Bin2DPowderDiffraction.cpp

#include <fstream>
#include <algorithm>
#include "MantidAlgorithms/Bin2DPowderDiffraction.h"
#include "MantidAPI/BinEdgeAxis.h"
#include "MantidAPI/FileProperty.h"
#include "MantidAPI/InstrumentValidator.h"
#include "MantidAPI/SpectrumDetectorMapping.h"
#include "MantidAPI/SpectrumInfo.h"
#include "MantidAPI/SpectraAxisValidator.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidAPI/WorkspaceUnitValidator.h"
#include "MantidKernel/CompositeValidator.h"
#include "MantidAPI/WorkspaceProperty.h"
#include "MantidKernel/UnitFactory.h"
#include "MantidDataObjects/EventList.h"
#include "MantidDataObjects/FractionalRebinning.h"
#include "MantidDataObjects/RebinnedOutput.h"
#include "MantidDataObjects/WorkspaceCreation.h"
#include "MantidGeometry/Math/ConvexPolygon.h"
#include "MantidGeometry/Math/PolygonIntersection.h"
#include "MantidGeometry/Math/Quadrilateral.h"
#include "MantidKernel/ArrayProperty.h"
#include "MantidKernel/PropertyWithValue.h"
#include "MantidKernel/RebinParamsValidator.h"
#include "MantidKernel/VectorHelper.h"

namespace Mantid {
namespace Algorithms {

using namespace Kernel;
using namespace API;
using namespace Geometry;
using namespace DataObjects;
using namespace Mantid::HistogramData;


// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(Bin2DPowderDiffraction)

//----------------------------------------------------------------------------------------------

/// Algorithms name for identification. @see Algorithm::name
const std::string Bin2DPowderDiffraction::name() const { return "Bin2DPowderDiffraction"; }

/// Algorithm's version for identification. @see Algorithm::version
int Bin2DPowderDiffraction::version() const { return 1; }

/// Algorithm's category for identification. @see Algorithm::category
const std::string Bin2DPowderDiffraction::category() const {
    return "Diffraction\\Focussing";
}

/// Algorithm's summary for use in the GUI and help. @see Algorithm::summary
const std::string Bin2DPowderDiffraction::summary() const {
    return "Bins TOF powder diffraction event data in 2D space.";
}

//----------------------------------------------------------------------------------------------
/** Initialize the algorithm's properties.
 */
void Bin2DPowderDiffraction::init() {
    auto wsValidator = boost::make_shared<CompositeValidator>();
    wsValidator->add<WorkspaceUnitValidator>("Wavelength");
    wsValidator->add<SpectraAxisValidator>();
    wsValidator->add<InstrumentValidator>();

    declareProperty(make_unique<WorkspaceProperty<EventWorkspace>>(
                        "InputWorkspace", "", Direction::Input,
                        wsValidator),
                    "An input EventWorkspace. X-axis units must be wavelength.");

    declareProperty(
                make_unique<WorkspaceProperty<API::Workspace>>("OutputWorkspace", "",
                                                               Direction::Output),
                "An output workspace.");

    const std::string docString =
            "A comma separated list of first bin boundary, width, last bin boundary. "
            "Optionally "
            "this can be followed by a comma and more widths and last boundary "
            "pairs. "
            "Negative width values indicate logarithmic binning.";
    auto rebinValidator = boost::make_shared<RebinParamsValidator>(true);
    declareProperty(make_unique<ArrayProperty<double>>("Axis1Binning",
                                                       rebinValidator),
                    docString);
    declareProperty(make_unique<ArrayProperty<double>>("Axis2Binning",
                                                       rebinValidator),
                    docString);

    const std::vector<std::string> exts{".txt", ".dat"};
    declareProperty(make_unique<FileProperty>("BinEdgesFile", "",
                                              FileProperty::OptionalLoad, exts),
                    "Optional: The ascii file containing the list of bin edges. "
                    "Either this or Axis1- and Axis2Binning needs to be specified.");

    declareProperty(
          Kernel::make_unique<PropertyWithValue<bool>>("NormalizeByBinArea", true),
          "Normalize the binned workspace by the bin area.");
}

//----------------------------------------------------------------------------------------------
/** Execute the algorithm.
 */
void Bin2DPowderDiffraction::exec() {
    m_inputWS = this->getProperty("InputWorkspace");
    m_numberOfSpectra = static_cast<int>(m_inputWS->getNumberHistograms());
    g_log.debug() << "Number of spectra in input workspace: " << m_numberOfSpectra
                  << "\n";

    const auto &oldXEdges = m_inputWS->x(0);
    BinEdges newXBins(oldXEdges.size());
    BinEdges newYBins(oldXEdges.size());

    MatrixWorkspace_sptr outputWS = createOutputWorkspace(newXBins, newYBins);

    bool normalizeByBinArea = this->getProperty("NormalizeByBinArea");
    if (normalizeByBinArea)
        normalizeToBinArea(outputWS);

    setProperty("OutputWorkspace", outputWS);
}
//----------------------------------------------------------------------------------------------
/**
 * @brief Bin2DPowderDiffraction::validateInputs Validate inputs
 * @return
 */
std::map<std::string, std::string> Bin2DPowderDiffraction::validateInputs() {
    std::map<std::string, std::string> result;

    int numWays = 0;

    const std::string beFileName = getProperty("BinEdgesFile");
    if (!beFileName.empty())
        numWays += 1;

    const std::vector<double> axis1Binning = getProperty("Axis1Binning");
    const std::vector<double> axis2Binning = getProperty("Axis2Binning");

    if (!axis1Binning.empty() && !axis2Binning.empty())
        numWays += 1;


    std::string message;
    if (numWays == 0) {
        message = "You must specify either Axis1Binning and Axis2Binning, "
                  "or a BinEdgesFile.";
    }
    if (numWays > 1) {
        message = "You must specify either Axis1Binning and Axis2Binning, "
                  "or a BinEdgesFile, but not both.";
    }

    if (!message.empty()) {
        result["Axis1Binning"] = message;
        result["BinEdgesFile"] = message;
    }

    return result;
}
//----------------------------------------------------------------------------------------------
/**
 * @brief createOutputWorkspace create an output workspace and setup axis
 * @param parent
 * @param newXBins
 * @param newYBins
 * @return
 */

MatrixWorkspace_sptr Bin2DPowderDiffraction::createOutputWorkspace(
        HistogramData::BinEdges &newXBins, HistogramData::BinEdges &newYBins) {

    using VectorHelper::createAxisFromRebinParams;
    bool binsFromFile(false);
    int newYSize = 0;
    size_t newXSize = 0;
    MatrixWorkspace_sptr outputWS;
    const auto &spectrumInfo = m_inputWS->spectrumInfo();

    const std::string beFileName = getProperty("BinEdgesFile");
    if (!beFileName.empty())
        binsFromFile=true;

    auto &newY = newYBins.mutableRawData();
    std::vector<std::vector<double>> fileXbins;

    // First create the output Workspace filled with zeros
    if (binsFromFile) {
        newY.clear();
        ReadBinsFromFile(newY, fileXbins);
        newYSize = static_cast<int>(newY.size());
        // unify xbins
        newXSize = UnifyXBins(fileXbins);
        g_log.debug() << "Maximal size of Xbins = " << newXSize;
        outputWS =  WorkspaceFactory::Instance().create(m_inputWS, newYSize-1, newXSize, newXSize-1);
        g_log.debug() << "Outws has " << outputWS->getNumberHistograms() << " histograms and " << outputWS->blocksize() << " bins." << std::endl;

        size_t idx = 0;
        for (auto Xbins : fileXbins) {
            g_log.debug() << "Xbins size: " << Xbins.size() << std::endl;
            BinEdges binEdges (Xbins);
            outputWS->setBinEdges(idx, binEdges);
            idx++;
        }

    } else {
        static_cast<void>(createAxisFromRebinParams(getProperty("Axis1Binning"),
                                                    newXBins.mutableRawData()));
        HistogramData::BinEdges binEdges(newXBins);
        newYSize =
                createAxisFromRebinParams(getProperty("Axis2Binning"), newY);
        newXSize = binEdges.size();
        outputWS =  WorkspaceFactory::Instance().create(m_inputWS, newYSize - 1, newXSize, newXSize-1);
        for (auto idx=0; idx<newYSize-1; idx++)
            outputWS->setBinEdges(idx, binEdges);
        NumericAxis *const abscissa = new BinEdgeAxis(newXBins.mutableRawData());
        outputWS->replaceAxis(0, abscissa);
    }

    outputWS->getAxis(0)->unit() = UnitFactory::Instance().create("dSpacing");

    NumericAxis *const verticalAxis = new BinEdgeAxis(newY);
    // Meta data
    verticalAxis->unit() = UnitFactory::Instance().create("dSpacingPerpendicular");
    verticalAxis->title() = "d_p";
    outputWS->replaceAxis(1, verticalAxis);

    Progress prog(this, 0.0, 1.0, m_numberOfSpectra);
    int64_t numSpectra = static_cast<int64_t>(m_numberOfSpectra);
    std::vector<std::vector<double>> newYValues(newYSize-1, std::vector<double>(newXSize-1, 0.0));
    std::vector<std::vector<double>> newEValues(newYSize-1, std::vector<double>(newXSize-1, 0.0));

    // fill the workspace with data
    g_log.notice() << "newYSize = " << newYSize << std::endl;
    g_log.notice() << "newXSize = " << newXSize << std::endl;
    std::vector<double> dp_vec (verticalAxis->getValues());

    PARALLEL_FOR_IF(Kernel::threadSafe(*m_inputWS, *outputWS))
    for (int64_t snum=0; snum<numSpectra; ++snum) {
        PARALLEL_START_INTERUPT_REGION
        double theta = 0.5*spectrumInfo.twoTheta(snum);
        EventList &evList = m_inputWS->getSpectrum(snum);

        // Switch to weighted if needed.
        if (evList.getEventType() == TOF)
            evList.switchTo(WEIGHTED);

        std::vector<WeightedEvent> events = evList.getWeightedEvents();

        for(const auto &ev:events){
            double d, dp;
            convertToDSpacing(ev.tof(), theta, &d, &dp);
            std::vector<double>::iterator upy = std::lower_bound(dp_vec.begin(), dp_vec.end(), dp);
            // long int h_index = (upy-dp_vec.begin()) -1;
            long int h_index = std::distance(dp_vec.begin(), upy) -1;
            if ((h_index < newYSize -1) && h_index > -1) {
                if (h_index == newYSize-1)
                    g_log.error("h_index is equal to the size of the Y axis!");
                auto xs = binsFromFile ? fileXbins[h_index] : newXBins.rawData();
                std::vector<double>::iterator lowx = std::lower_bound(xs.begin(), xs.end(), d);
                long int index = std::distance(xs.begin(), lowx) -1;
                if ((index <  static_cast<int>(newXSize-1)) && (index > -1)) {
                    // writing to the same vectors is not threat-safe
                    PARALLEL_CRITICAL(newValues) {
                        newYValues[h_index][index] += ev.weight();
                        newEValues[h_index][index] += ev.errorSquared();
                    }
                }
                // g_log.notice() << "wl = " << ev.tof() << ", theta = " << theta << ", d = " << d << ", dp = " << dp << std::endl;
                // g_log.notice() << "h_index = " << static_cast<int>(h_index) << ", x_index = " << static_cast<int>(index) << std::endl;

            }

        }
        prog.report("Binning event data...");
        PARALLEL_END_INTERUPT_REGION

    }
    PARALLEL_CHECK_INTERUPT_REGION
    size_t idx=0;
    for (const auto &yVec: newYValues) {
        outputWS->setCounts(idx, yVec);
        idx++;
    }
    idx=0;
    for (auto &eVec: newEValues) {
        std::transform(eVec.begin(), eVec.end(), eVec.begin(),
                           static_cast<double (*)(double)>(sqrt));
        outputWS->setCountStandardDeviations(idx, eVec);
        idx++;
    }
    return outputWS;
}

//----------------------------------------------------------------------------------------------
/**
 * @brief Bin2DPowderDiffraction::ReadBinsFromFile
 * @param[out] Ybins vector of doubles to save the dOrth bin parameters
 * @param[out] Xbins vector of vectors of doubles to sabe the dSpacing bin edges
 */
void Bin2DPowderDiffraction::ReadBinsFromFile(std::vector<double> &Ybins, std::vector<std::vector<double> > &Xbins) const
{
    const std::string beFileName = getProperty("BinEdgesFile");
    std::ifstream file (beFileName);
    std::string line;
    std::string::size_type n;
    std::string::size_type sz;
    std::vector<double> tmp;
    int dpno = 0;
    while(getline(file, line)){
        n = line.find("dp =");
        if (n != std::string::npos) {
            if (!tmp.empty()){
                Xbins.push_back(tmp);
                tmp.clear();
            }
            double dp1 = std::stod (line.substr(4),&sz);        // 4 is needed to crop 'dp='
            double dp2 = std::stod (line.substr(sz + 4));
            if (dpno < 1){
                Ybins.push_back(dp1);
                Ybins.push_back(dp2);
            } else {
                Ybins.push_back(dp2);
            }
            dpno++;
        } else if(line.find("#")==std::string::npos) {
            std::stringstream ss(line);
            double d;
            while (ss >> d)
            {
                tmp.push_back(d);
            }

        }
    }
    Xbins.push_back(tmp);
    g_log.information() << "Number of Ybins: " <<  Ybins.size() << std::endl;
    g_log.information() << "Number of Xbins sets: " <<  Xbins.size() << std::endl;

}
//----------------------------------------------------------------------------------------------
/**
 * @brief Bin2DPowderDiffraction::UnifyXBins unifies size of the vectors in Xbins.
 * Just fills std::nans at the end of the shorter bins.
 * Required to avoid garbage values in the X values after ws->setHistogram.
 * returns the maximal size
 *
 * @param Xbins[in] --- bins to unify. Will be overwritten.
 */
size_t Bin2DPowderDiffraction::UnifyXBins(std::vector<std::vector<double> > &Xbins) const
{
    // get maximal vector size
    size_t max_size = 0;
    for (const auto& v : Xbins) {
            max_size = std::max(v.size(), max_size);
    }
    // resize all vectors to maximum size, fill last vector element at the end
    for (auto &v: Xbins) {
        if (v.size() < max_size)
            v.resize(max_size, v.back());
    }
    return max_size;

}

void Bin2DPowderDiffraction::normalizeToBinArea(MatrixWorkspace_sptr outWS)
{
    NumericAxis *verticalAxis = dynamic_cast<NumericAxis *>(outWS->getAxis(1));
    const std::vector<double> yValues = verticalAxis->getValues();
    auto nhist = outWS->getNumberHistograms();
    g_log.debug() << "Number of hists: " << nhist << " Length of YAxis: " << verticalAxis->length() << std::endl;

    for (size_t idx=0; idx<nhist; ++idx){
        double factor = 1.0/(yValues[idx+1] - yValues[idx]);
        // divide by the xBinWidth
        outWS->convertToFrequencies(idx);
        auto &freqs = outWS->mutableY(idx);
        std::transform(freqs.begin(), freqs.end(), freqs.begin(),
                       std::bind1st(std::multiplies<double>(), factor));
        auto &errors = outWS->mutableE(idx);
        std::transform(errors.begin(), errors.end(), errors.begin(),
                       std::bind1st(std::multiplies<double>(), factor));

    }
}

//TODO: take care of theta=0.
void convertToDSpacing(double wavelength, double theta, double *d, double *dp)
{
    *d = wavelength*0.5/sin(theta);
    *dp = sqrt(wavelength*wavelength - 2.0*log(cos(theta)));
}

} // namespace Algorithms
} // namespace Mantid