From da74966683e34160a95c9d4b2cd5c421d1bee4dd Mon Sep 17 00:00:00 2001 From: Andrei Savici Date: Tue, 9 Dec 2014 17:31:57 -0500 Subject: [PATCH] Start working on MDNormDirectSC. Refs #10698 --- .../Framework/MDAlgorithms/CMakeLists.txt | 3 + .../MantidMDAlgorithms/CalculateCoverageDGS.h | 2 +- .../inc/MantidMDAlgorithms/MDNormDirectSC.h | 97 ++ .../MDAlgorithms/src/MDNormDirectSC.cpp | 840 ++++++++++++++++++ .../MDAlgorithms/test/MDNormDirectSCTest.h | 119 +++ 5 files changed, 1060 insertions(+), 1 deletion(-) create mode 100644 Code/Mantid/Framework/MDAlgorithms/inc/MantidMDAlgorithms/MDNormDirectSC.h create mode 100644 Code/Mantid/Framework/MDAlgorithms/src/MDNormDirectSC.cpp create mode 100644 Code/Mantid/Framework/MDAlgorithms/test/MDNormDirectSCTest.h diff --git a/Code/Mantid/Framework/MDAlgorithms/CMakeLists.txt b/Code/Mantid/Framework/MDAlgorithms/CMakeLists.txt index 6443cf9d483c..9ef717f1dc21 100644 --- a/Code/Mantid/Framework/MDAlgorithms/CMakeLists.txt +++ b/Code/Mantid/Framework/MDAlgorithms/CMakeLists.txt @@ -40,6 +40,7 @@ set ( SRC_FILES src/LoadMD.cpp src/LoadSQW.cpp src/LogarithmMD.cpp + src/MDNormDirectSC.cpp src/MDNormSCD.cpp src/MaskMD.cpp src/MergeMD.cpp @@ -123,6 +124,7 @@ set ( INC_FILES inc/MantidMDAlgorithms/LoadMD.h inc/MantidMDAlgorithms/LoadSQW.h inc/MantidMDAlgorithms/LogarithmMD.h + inc/MantidMDAlgorithms/MDNormDirectSC.h inc/MantidMDAlgorithms/MDNormSCD.h inc/MantidMDAlgorithms/MaskMD.h inc/MantidMDAlgorithms/MergeMD.h @@ -205,6 +207,7 @@ set ( TEST_FILES LoadMDTest.h LoadSQWTest.h LogarithmMDTest.h + MDNormDirectSCTest.h MDNormSCDTest.h MDResolutionConvolutionFactoryTest.h MaskMDTest.h diff --git a/Code/Mantid/Framework/MDAlgorithms/inc/MantidMDAlgorithms/CalculateCoverageDGS.h b/Code/Mantid/Framework/MDAlgorithms/inc/MantidMDAlgorithms/CalculateCoverageDGS.h index 40e73e9a1f1c..20971c37f7a3 100644 --- a/Code/Mantid/Framework/MDAlgorithms/inc/MantidMDAlgorithms/CalculateCoverageDGS.h +++ b/Code/Mantid/Framework/MDAlgorithms/inc/MantidMDAlgorithms/CalculateCoverageDGS.h @@ -49,7 +49,7 @@ namespace MDAlgorithms /// limits for h,k,l,dE dimensions coord_t m_hmin, m_hmax, m_kmin, m_kmax, m_lmin, m_lmax, m_dEmin, m_dEmax; - /// cached values for incident energy and momentum + /// cached values for incident energy and momentum, final momentum min/max double m_Ei,m_ki, m_kfmin,m_kfmax; /// flag for integrated h,k,l,dE dimensions bool m_hIntegrated, m_kIntegrated, m_lIntegrated, m_dEIntegrated; diff --git a/Code/Mantid/Framework/MDAlgorithms/inc/MantidMDAlgorithms/MDNormDirectSC.h b/Code/Mantid/Framework/MDAlgorithms/inc/MantidMDAlgorithms/MDNormDirectSC.h new file mode 100644 index 000000000000..6b4bde6054fb --- /dev/null +++ b/Code/Mantid/Framework/MDAlgorithms/inc/MantidMDAlgorithms/MDNormDirectSC.h @@ -0,0 +1,97 @@ +#ifndef MANTID_MDALGORITHMS_MDNORMDIRECTSC_H_ +#define MANTID_MDALGORITHMS_MDNORMDIRECTSC_H_ + +#include "MantidAPI/Algorithm.h" +#include "MantidMDAlgorithms/SlicingAlgorithm.h" + +namespace Mantid +{ +namespace DataObjects +{ + class EventWorkspace; +} + namespace MDAlgorithms + { + + /** MDNormSCD : Generate MD normalization for single crystal diffraction + + Copyright © 2014 ISIS Rutherford Appleton Laboratory, NScD Oak Ridge National Laboratory & European Spallation Source + + This file is part of Mantid. + + Mantid is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 3 of the License, or + (at your option) any later version. + + Mantid is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . + + File change history is stored at: + Code Documentation is available at: + */ + class DLLExport MDNormDirectSC :public SlicingAlgorithm + { + public: + MDNormDirectSC(); + + virtual const std::string name() const; + virtual int version() const; + virtual const std::string category() const; + virtual const std::string summary() const; + + private: + void init(); + void exec(); + + void cacheInputs(); + std::string inputEnergyMode() const; + + MDEvents::MDHistoWorkspace_sptr binInputWS(); + void createNormalizationWS(const MDEvents::MDHistoWorkspace & dataWS); + std::vector getValuesFromOtherDimensions(bool & skipNormalization) const; + Kernel::Matrix findIntergratedDimensions(const std::vector & otherDimValues, + bool & skipNormalization); + void cacheDimensionXValues(); + void calculateNormalization(const std::vector &otherValues, + const Kernel::Matrix &affineTrans); + void calcIntegralsForIntersections( const std::vector &xValues, const API::MatrixWorkspace &integrFlux, + size_t sp, std::vector &yValues ) const; + std::vector removeGroupedIDs(const API::ExperimentInfo & exptInfo, + const std::vector &detIDs); + Geometry::IDetector_const_sptr getThetaPhi(const detid_t detID, + const API::ExperimentInfo & exptInfo, + double &theta, double &phi); + std::vector calculateIntersections(const double theta, const double phi); + + /// Normalization workspace + MDEvents::MDHistoWorkspace_sptr m_normWS; + /// Input workspace + API::IMDEventWorkspace_sptr m_inputWS; + /// limits for h,k,l, dE dimensions + coord_t m_hmin, m_hmax, m_kmin, m_kmax, m_lmin, m_lmax, m_dEmin, m_dEmax; + /// cached values for incident energy and momentum, final momentum min/max + double m_Ei,m_ki, m_kfmin,m_kfmax; + /// flag for integrated h,k,l, dE dimensions + bool m_hIntegrated, m_kIntegrated, m_lIntegrated, m_dEIntegrated; + /// (2*PiRUBW)^-1 + Mantid::Kernel::DblMatrix m_rubw; + ///index of h,k,l, dE dimensions in the output workspaces + size_t m_hIdx, m_kIdx, m_lIdx, m_eIdx; + /// cached X values along dimensions h,k,l. dE + std::vector m_hX, m_kX, m_lX, m_eX; + /// Sample position + Kernel::V3D m_samplePos; + /// Beam direction + Kernel::V3D m_beamDir; + }; + + } // namespace MDAlgorithms +} // namespace Mantid + +#endif /* MANTID_MDALGORITHMS_MDNORMDIRECTSC_H_ */ diff --git a/Code/Mantid/Framework/MDAlgorithms/src/MDNormDirectSC.cpp b/Code/Mantid/Framework/MDAlgorithms/src/MDNormDirectSC.cpp new file mode 100644 index 000000000000..8ba802941784 --- /dev/null +++ b/Code/Mantid/Framework/MDAlgorithms/src/MDNormDirectSC.cpp @@ -0,0 +1,840 @@ +#include "MantidMDAlgorithms/MDNormDirectSC.h" + +#include "MantidAPI/WorkspaceValidators.h" +#include "MantidDataObjects/EventWorkspace.h" +#include "MantidMDEvents/MDEventWorkspace.h" +#include "MantidMDEvents/MDHistoWorkspace.h" +#include "MantidKernel/TimeSeriesProperty.h" +#include "MantidKernel/VectorHelper.h" + +namespace Mantid +{ + namespace MDAlgorithms + { + + using Mantid::Kernel::Direction; + using Mantid::API::WorkspaceProperty; + using namespace Mantid::MDEvents; + using namespace Mantid::API; + using namespace Mantid::Kernel; + + namespace + { + //function to compare two intersections (h,k,l,Momentum) by Momentum + bool compareMomentum(const Mantid::Kernel::VMD &v1, const Mantid::Kernel::VMD &v2) + { + return (v1[3] < v2[3]); + } + } + + // Register the algorithm into the AlgorithmFactory + DECLARE_ALGORITHM(MDNormDirectSC) + + //---------------------------------------------------------------------------------------------- + /** + * Constructor + */ + MDNormDirectSC::MDNormDirectSC() : + m_normWS(), m_inputWS(), m_hmin(0.0f), m_hmax(0.0f), + m_kmin(0.0f), m_kmax(0.0f), m_lmin(0.0f), m_lmax(0.0f), m_dEmin(0.f), m_dEmax(0.f), + m_Ei(0.),m_ki(0.), m_kfmin(0.), m_kfmax(0.), + m_hIntegrated(true),m_kIntegrated(true), m_lIntegrated(true), m_dEIntegrated(false), m_rubw(3,3), + m_hIdx(-1), m_kIdx(-1), m_lIdx(-1),m_eIdx(-1), + m_hX(), m_kX(), m_lX(), m_eX(), m_samplePos(), m_beamDir() + { + } + + /// Algorithm's version for identification. @see Algorithm::version + int MDNormDirectSC::version() const { return 1; } + + /// Algorithm's category for identification. @see Algorithm::category + const std::string MDNormDirectSC::category() const { return "MDAlgorithms"; } + + /// Algorithm's summary for use in the GUI and help. @see Algorithm::summary + const std::string MDNormDirectSC::summary() const + { + return "Calculate normalization for an MDEvent workspace for single crystal diffraction."; + } + + /// Algorithm's name for use in the GUI and help. @see Algorithm::name + const std::string MDNormDirectSC::name() const { return "MDNormDirectSC"; } + + //---------------------------------------------------------------------------------------------- + /** + * Initialize the algorithm's properties. + */ + void MDNormDirectSC::init() + { + declareProperty(new WorkspaceProperty("InputWorkspace","",Direction::Input), + "An input MDWorkspace."); + + std::string dimChars = getDimensionChars(); + // --------------- Axis-aligned properties --------------------------------------- + for (size_t i = 0; i(propName,"",Direction::Input), + "Binning parameters for the " + Strings::toString(i) + "th dimension.\n" + "Enter it as a comma-separated list of values with the format: 'name,minimum,maximum,number_of_bins'. Leave blank for NONE."); + } + + auto solidAngleValidator = boost::make_shared(); + solidAngleValidator->add(); + solidAngleValidator->add(); + + + declareProperty(new WorkspaceProperty<>("SolidAngleWorkspace","",Direction::Input, solidAngleValidator), + "An input workspace containing integrated vanadium (a measure of the solid angle)."); + + declareProperty(new WorkspaceProperty("OutputWorkspace","",Direction::Output), + "A name for the output data MDHistoWorkspace."); + declareProperty(new WorkspaceProperty("OutputNormalizationWorkspace","",Direction::Output), + "A name for the output normalization MDHistoWorkspace."); + } + + //---------------------------------------------------------------------------------------------- + /** + * Execute the algorithm. + */ + void MDNormDirectSC::exec() + { + cacheInputs(); + auto outputWS = binInputWS(); + setProperty("OutputWorkspace", outputWS); + createNormalizationWS(*outputWS); + setProperty("OutputNormalizationWorkspace", m_normWS); + + // Check for other dimensions if we could measure anything in the original data + bool skipNormalization = false; + const std::vector otherValues = getValuesFromOtherDimensions(skipNormalization); + const auto affineTrans = findIntergratedDimensions(otherValues, skipNormalization); + cacheDimensionXValues(); + + if(!skipNormalization) + { + calculateNormalization(otherValues, affineTrans); + } + else + { + g_log.warning("Binning limits are outside the limits of the MDWorkspace. Not applying normalization."); + } + } + + /** + * Set up starting values for cached variables + */ + void MDNormDirectSC::cacheInputs() + { + m_inputWS = getProperty("InputWorkspace"); + if( inputEnergyMode() != "Direct" ) + { + throw std::invalid_argument("Invalid energy transfer mode. Algorithm only supports direct geometry spectrometers."); + } + // Min/max dimension values + const auto hdim(m_inputWS->getDimension(0)), kdim(m_inputWS->getDimension(1)), + ldim(m_inputWS->getDimension(2)),edim(m_inputWS->getDimension(3)); + m_hmin = hdim->getMinimum(); + m_kmin = kdim->getMinimum(); + m_lmin = ldim->getMinimum(); + m_dEmin = edim->getMinimum(); + m_hmax = hdim->getMaximum(); + m_kmax = kdim->getMaximum(); + m_lmax = ldim->getMaximum(); + m_dEmax = edim->getMaximum(); + + const auto & exptInfoZero = *(m_inputWS->getExperimentInfo(0)); + auto source = exptInfoZero.getInstrument()->getSource(); + auto sample = exptInfoZero.getInstrument()->getSample(); + if(source == NULL || sample == NULL) + { + throw Kernel::Exception::InstrumentDefinitionError("Instrument not sufficiently defined: failed to get source and/or sample"); + } + m_samplePos = sample->getPos(); + m_beamDir = m_samplePos - source->getPos(); + m_beamDir.normalize(); + } + + /** + * Currently looks for the ConvertToMD algorithm in the history + * @return A string donating the energy transfer mode of the input workspace + */ + std::string MDNormDirectSC::inputEnergyMode() const + { + const auto & hist = m_inputWS->getHistory(); + const size_t nalgs = hist.size(); + const auto & lastAlgorithm = hist.lastAlgorithm(); + + std::string emode(""); + if(lastAlgorithm->name() == "ConvertToMD") + { + emode = lastAlgorithm->getPropertyValue("dEAnalysisMode"); + } + else if ( (lastAlgorithm->name() == "Load" || hist.lastAlgorithm()->name() == "LoadMD") && + hist.getAlgorithmHistory(nalgs - 2)->name() == "ConvertToMD" ) + { + //get dEAnalysisMode + PropertyHistories histvec = hist.getAlgorithmHistory(nalgs - 2)->getProperties(); + for(auto it = histvec.begin(); it != histvec.end(); ++it) + { + if((*it)->name() == "dEAnalysisMode") + { + emode = (*it)->value(); + break; + } + } + } + else + { + throw std::invalid_argument("The last algorithm in the history of the input workspace is not ConvertToMD"); + } + return emode; + } + + /** + * Runs the BinMD algorithm on the input to provide the output workspace + * All slicing algorithm properties are passed along + * @return MDHistoWorkspace as a result of the binning + */ + MDHistoWorkspace_sptr MDNormDirectSC::binInputWS() + { + const auto & props = getProperties(); + IAlgorithm_sptr binMD = createChildAlgorithm("BinMD", 0.0, 0.3); + binMD->setPropertyValue("AxisAligned","1"); + for(auto it = props.begin(); it != props.end(); ++it) + { + const auto & propName = (*it)->name(); + if(propName != "FluxWorkspace" && propName != "SolidAngleWorkspace" && + propName != "OutputNormalizationWorkspace") + { + binMD->setPropertyValue(propName,(*it)->value()); + } + } + binMD->executeAsChildAlg(); + Workspace_sptr outputWS = binMD->getProperty("OutputWorkspace"); + return boost::dynamic_pointer_cast(outputWS); + } + + /** + * Create & cached the normalization workspace + * @param dataWS The binned workspace that will be used for the data + */ + void MDNormDirectSC::createNormalizationWS(const MDHistoWorkspace &dataWS) + { + // Copy the MDHisto workspace, and change signals and errors to 0. + m_normWS = boost::make_shared(dataWS); + m_normWS->setTo(0.,0.,0.); + } + + /** + * Retrieve logged values from non-HKL dimensions + * @param skipNormalization [InOut] Updated to false if any values are outside range measured by input workspace + * @return A vector of values from other dimensions to be include in normalized MD position calculation + */ + std::vector MDNormDirectSC::getValuesFromOtherDimensions(bool &skipNormalization) const + { + const auto & runZero = m_inputWS->getExperimentInfo(0)->run(); + + std::vector otherDimValues; + for(size_t i = 3; i < m_inputWS->getNumDims(); i++) + { + const auto dimension = m_inputWS->getDimension(i); + float dimMin = static_cast(dimension->getMinimum()); + float dimMax = static_cast(dimension->getMaximum()); + auto *dimProp = \ + dynamic_cast *>(runZero.getProperty(dimension->getName())); + if (dimProp) + { + coord_t value = static_cast(dimProp->firstValue()); + otherDimValues.push_back(value); + //in the original MD data no time was spent measuring between dimMin and dimMax + if (value < dimMin || value > dimMax) + { + skipNormalization = true; + } + } + } + return otherDimValues; + } + + /** + * Checks the normalization workspace against the indices of the original dimensions. + * If not found, the corresponding dimension is integrated + * @param otherDimValues Values from non-HKL dimensions + * @param skipNormalization [InOut] Sets the flag true if normalization values are outside of original inputs + * @return Affine trasform matrix + */ + Kernel::Matrix MDNormDirectSC::findIntergratedDimensions(const std::vector &otherDimValues, + bool &skipNormalization) + { + // Get indices of the original dimensions in the output workspace, + // and if not found, the corresponding dimension is integrated + Kernel::Matrix affineMat = m_normWS->getTransformFromOriginal(0)->makeAffineMatrix(); + + const size_t nrm1 = affineMat.numRows() - 1; + const size_t ncm1 = affineMat.numCols() - 1; + for( size_t row = 0; row < nrm1; row++ ) //affine matrix, ignore last row + { + const auto dimen = m_normWS->getDimension(row); + const auto dimMin(dimen->getMinimum()), dimMax(dimen->getMaximum()); + if(affineMat[row][0] == 1.) + { + m_hIntegrated = false; + m_hIdx = row; + if(m_hmin < dimMin) m_hmin = dimMin; + if(m_hmax > dimMax) m_hmax = dimMax; + if(m_hmin > dimMax || m_hmax < dimMin) + { + skipNormalization = true; + } + } + if(affineMat[row][1] == 1.) + { + m_kIntegrated = false; + m_kIdx = row; + if(m_kmin < dimMin) m_kmin = dimMin; + if(m_kmax > dimMax) m_kmax = dimMax; + if(m_kmin > dimMax || m_kmax < dimMin) + { + skipNormalization = true; + } + } + if(affineMat[row][2] == 1.) + { + m_lIntegrated = false; + m_lIdx = row; + if(m_lmin < dimMin) m_lmin = dimMin; + if(m_lmax > dimMax) m_lmax = dimMax; + if(m_lmin > dimMax || m_lmax < dimMin) + { + skipNormalization = true; + } + + for(size_t col = 3; col < ncm1; col++) //affine matrix, ignore last column + { + if(affineMat[row][col] == 1.) + { + double val = otherDimValues.at(col - 3); + if(val > dimMax || val < dimMin) + { + skipNormalization = true; + } + } + } + } + } + + return affineMat; + } + + /** + * Stores the X values from each H,K,L dimension as member variables + */ + void MDNormDirectSC::cacheDimensionXValues() + { + auto &hDim = *m_normWS->getDimension(m_hIdx); + m_hX.resize(hDim.getNBins()); + for(size_t i = 0; i < m_hX.size(); ++i) + { + m_hX[i] = hDim.getX(i); + } + auto &kDim = *m_normWS->getDimension(m_kIdx); + m_kX.resize( kDim.getNBins() ); + for(size_t i = 0; i < m_kX.size(); ++i) + { + m_kX[i] = kDim.getX(i); + } + auto &lDim = *m_normWS->getDimension(m_lIdx); + m_lX.resize( lDim.getNBins() ); + for(size_t i = 0; i < m_lX.size(); ++i) + { + m_lX[i] = lDim.getX(i); + } + } + + /** + * Computed the normalization for the input workspace. Results are stored in m_normWS + * @param otherValues + * @param affineTrans + */ + void MDNormDirectSC::calculateNormalization(const std::vector &otherValues, + const Kernel::Matrix &affineTrans) + { + /*API::MatrixWorkspace_const_sptr integrFlux = getProperty("FluxWorkspace"); + integrFlux->getXMinMax(m_kiMin, m_kiMax); + API::MatrixWorkspace_const_sptr solidAngleWS = getProperty("SolidAngleWorkspace"); + + const auto & exptInfoZero = *(m_inputWS->getExperimentInfo(0)); + typedef Kernel::PropertyWithValue > VectorDoubleProperty; + auto *rubwLog = dynamic_cast(exptInfoZero.getLog("RUBW_MATRIX")); + if(!rubwLog) + { + throw std::runtime_error("Wokspace does not contain a log entry for the RUBW matrix." + "Cannot continue."); + } + else + { + Kernel::DblMatrix rubwValue((*rubwLog)()); //includes the 2*pi factor but not goniometer for now :) + m_rubw = exptInfoZero.run().getGoniometerMatrix()*rubwValue; + m_rubw.Invert(); + } + const double protonCharge = exptInfoZero.run().getProtonCharge(); + + auto instrument = exptInfoZero.getInstrument(); + std::vector detIDs = instrument->getDetectorIDs(true); + // Prune out those that are part of a group and simply leave the head of the group + detIDs = removeGroupedIDs(exptInfoZero, detIDs); + + // Mappings + const int64_t ndets = static_cast(detIDs.size()); + const detid2index_map fluxDetToIdx = integrFlux->getDetectorIDToWorkspaceIndexMap(); + const detid2index_map solidAngDetToIdx = solidAngleWS->getDetectorIDToWorkspaceIndexMap(); + + auto *prog = new API::Progress(this, 0.3, 1.0, ndets); + PARALLEL_FOR1(integrFlux) + for(int64_t i = 0; i < ndets; i++) + { + PARALLEL_START_INTERUPT_REGION + + const auto detID = detIDs[i]; + double theta(0.0), phi(0.0); + bool skip(false); + try + { + auto spectrum = getThetaPhi(detID, exptInfoZero, theta, phi); + if(spectrum->isMonitor() || spectrum->isMasked()) continue; + } + catch(std::exception&) // detector might not exist or has no been included in grouping + { + skip = true;// Intel compiler has a problem with continue inside a catch inside openmp... + } + if(skip) continue; + + // Intersections + auto intersections = calculateIntersections(theta, phi); + if(intersections.empty()) continue; + + // get the flux spetrum number + size_t wsIdx = fluxDetToIdx.find(detID)->second; + // Get solid angle for this contribution + double solid = solidAngleWS->readY(solidAngDetToIdx.find(detID)->second)[0]*protonCharge; + + // -- calculate integrals for the intersection -- + // momentum values at intersections + auto intersectionsBegin = intersections.begin(); + std::vector xValues(intersections.size()), yValues(intersections.size()); + { + // copy momenta to xValues + auto x = xValues.begin(); + for (auto it = intersectionsBegin; it != intersections.end(); ++it, ++x) + { + *x = (*it)[3]; + } + } + // calculate integrals at momenta from xValues by interpolating between points in spectrum sp + // of workspace integrFlux. The result is stored in yValues + calcIntegralsForIntersections( xValues, *integrFlux, wsIdx, yValues ); + + // Compute final position in HKL + const size_t vmdDims = intersections.front().size(); + // pre-allocate for efficiency and copy non-hkl dim values into place + std::vector pos(vmdDims + otherValues.size()); + std::copy(otherValues.begin(), otherValues.end(), pos.begin() + vmdDims); + + for (auto it = intersectionsBegin + 1; it != intersections.end(); ++it) + { + const auto & curIntSec = *it; + const auto & prevIntSec = *(it-1); + // the full vector isn't used so compute only what is necessary + double delta = curIntSec[3] - prevIntSec[3]; + if(delta < 1e-07) continue; // Assume zero contribution if difference is small + + // Average between two intersections for final position + std::transform(curIntSec.getBareArray(), curIntSec.getBareArray() + vmdDims, + prevIntSec.getBareArray(), pos.begin(), + VectorHelper::SimpleAverage()); + std::vector posNew = affineTrans*pos; + size_t linIndex = m_normWS->getLinearIndexAtCoord(posNew.data()); + if(linIndex == size_t(-1)) continue; + + // index of the current intersection + size_t k = static_cast(std::distance(intersectionsBegin, it)); + // signal = integral between two consecutive intersections + double signal = (yValues[k] - yValues[k - 1])*solid; + + PARALLEL_CRITICAL(updateMD) + { + signal += m_normWS->getSignalAt(linIndex); + m_normWS->setSignalAt(linIndex, signal); + } + } + prog->report(); + + PARALLEL_END_INTERUPT_REGION + } + PARALLEL_CHECK_INTERUPT_REGION + + delete prog;*/ + } + + /** + * Linearly interpolate between the points in integrFlux at xValues and save the results in yValues. + * @param xValues :: X-values at which to interpolate + * @param integrFlux :: A workspace with the spectra to interpolate + * @param sp :: A workspace index for a spectrum in integrFlux to interpolate. + * @param yValues :: A vector to save the results. + */ + void MDNormDirectSC::calcIntegralsForIntersections( const std::vector &xValues, const API::MatrixWorkspace &integrFlux, size_t sp, std::vector &yValues ) const + { + assert( xValues.size() == yValues.size() ); + + // the x-data from the workspace + const auto &xData = integrFlux.readX(sp); + const double xStart = xData.front(); + const double xEnd = xData.back(); + + // the values in integrFlux are expected to be integrals of a non-negative function + // ie they must make a non-decreasing function + const auto &yData = integrFlux.readY(sp); + size_t spSize = yData.size(); + + const double yMin = 0.0; + const double yMax = yData.back(); + + size_t nData = xValues.size(); + // all integrals below xStart must be 0 + if (xValues[nData-1] < xStart) + { + std::fill( yValues.begin(), yValues.end(), yMin ); + return; + } + + // all integrals above xEnd must be equal tp yMax + if ( xValues[0] > xEnd ) + { + std::fill( yValues.begin(), yValues.end(), yMax ); + return; + } + + size_t i = 0; + // integrals below xStart must be 0 + while(i < nData - 1 && xValues[i] < xStart) + { + yValues[i] = yMin; + i++; + } + size_t j = 0; + for(; i < nData; i++) + { + // integrals above xEnd must be equal tp yMax + if (j >= spSize - 1) + { + yValues[i] = yMax; + } + else + { + double xi = xValues[i]; + while(j < spSize - 1 && xi > xData[j]) j++; + // if x falls onto an interpolation point return the corresponding y + if (xi == xData[j]) + { + yValues[i] = yData[j]; + } + else if (j == spSize - 1) + { + // if we get above xEnd it's yMax + yValues[i] = yMax; + } + else if (j > 0) + { + // interpolate between the consecutive points + double x0 = xData[j-1]; + double x1 = xData[j]; + double y0 = yData[j-1]; + double y1 = yData[j]; + yValues[i] = y0 + (y1 - y0)*(xi - x0)/(x1 - x0); + } + else // j == 0 + { + yValues[i] = yMin; + } + } + } + + } + + /** + * Checks for IDs that are actually part of the same group and just keeps one from the group. + * For a 1:1 map, none will be removed. + * @param exptInfo An ExperimentInfo object that defines the grouping + * @param detIDs A list of existing IDs + * @return A new list of IDs + */ + std::vector MDNormDirectSC::removeGroupedIDs(const ExperimentInfo &exptInfo, const std::vector &detIDs) + { + const size_t ntotal = detIDs.size(); + std::vector singleIDs; + singleIDs.reserve(ntotal/2); // reserve half. In the case of 1:1 it will double to the correct size once + std::set groupedIDs; + + for(auto iter = detIDs.begin(); iter != detIDs.end(); ++iter) + { + detid_t curID = *iter; + if(groupedIDs.count(curID) == 1) continue; // Already been processed + + try + { + const auto & members = exptInfo.getGroupMembers(curID); + singleIDs.push_back(members.front()); + std::copy(members.begin() + 1, members.end(), + std::inserter(groupedIDs, groupedIDs.begin())); + } + catch (std::runtime_error &) + { + singleIDs.push_back(curID); + } + } + + g_log.debug() << "Found " << singleIDs.size() << " spectra from " << detIDs.size() << " IDs\n"; + return singleIDs; + } + + /** + * Get the theta and phi angles for the given ID. If the detector was part of a group, + * as defined in the ExperimentInfo object, then the theta/phi are for the whole set. + * @param detID A reference to a single ID + * @param exptInfo A reference to the ExperimentInfo that defines that spectrum->detector mapping + * @param theta [Output] Set to the theta angle for the detector (set) + * @param phi [Output] Set to the phi angle for the detector (set) + * @return A poiner to the Detector object for this spectrum as a whole + * (may be a single pixel or group) + */ + Geometry::IDetector_const_sptr + MDNormDirectSC::getThetaPhi(const detid_t detID, + const ExperimentInfo &exptInfo, + double &theta, double &phi) + { + const auto spectrum = exptInfo.getDetectorByID(detID); + theta = spectrum->getTwoTheta(m_samplePos, m_beamDir); + phi = spectrum->getPhi(); + return spectrum; + } + + /** + * Calculate the points of intersection for the given detector with cuboid surrounding the + * detector position in HKL + * @param theta Polar angle withd detector + * @param phi Azimuthal angle with detector + * @return A list of intersections in HKL space + */ + std::vector MDNormDirectSC::calculateIntersections(const double theta, const double phi) + { + V3D qout(sin(theta)*cos(phi), sin(theta)*sin(phi), cos(theta)),qin(0.,0.,m_ki); + qout = m_rubw*qout; + qin = m_rubw*qin; + double hStart = qin.X()-qout.X()*m_kfmin, hEnd = qin.X()-qout.X()*m_kfmax; + double kStart = qin.Y()-qout.Y()*m_kfmin, kEnd = qin.Y()-qout.Y()*m_kfmax; + double lStart = qin.Z()-qout.Z()*m_kfmin, lEnd = qin.Z()-qout.Z()*m_kfmax; + double eps = 1e-10; + auto hNBins = m_hX.size(); + auto kNBins = m_kX.size(); + auto lNBins = m_lX.size(); + auto eNBins = m_eX.size(); + std::vector intersections; + intersections.reserve(hNBins + kNBins + lNBins + eNBins+8);//8 is 3*(min,max for each Q component)+kfmin+kfmax + + //calculate intersections with planes perpendicular to h + if (fabs(hStart-hEnd) > eps) + { + double fmom=(m_kfmax-m_kfmin)/(hEnd-hStart); + double fk=(kEnd-kStart)/(hEnd-hStart); + double fl=(lEnd-lStart)/(hEnd-hStart); + if(!m_hIntegrated) + { + for(size_t i = 0;i=m_hmin)&&(hi<=m_hmax) && ((hStart-hi)*(hEnd-hi)<0)) + { + // if hi is between hStart and hEnd, then ki and li will be between kStart, kEnd and lStart, lEnd and momi will be between m_kfmin and m_kfmax + double ki = fk*(hi-hStart)+kStart; + double li = fl*(hi-hStart)+lStart; + if ((ki>=m_kmin)&&(ki<=m_kmax)&&(li>=m_lmin)&&(li<=m_lmax)) + { + double momi = fmom*(hi-hStart)+m_kfmin; + Mantid::Kernel::VMD v(hi,ki,li,momi); + intersections.push_back(v); + } + } + } + } + double momhMin = fmom*(m_hmin-hStart)+m_kfmin; + if ((momhMin-m_kfmin)*(momhMin-m_kfmax)<0)//m_kfmin>m_kfmax + { + //khmin and lhmin + double khmin = fk*(m_hmin-hStart)+kStart; + double lhmin = fl*(m_hmin-hStart)+lStart; + if((khmin>=m_kmin)&&(khmin<=m_kmax)&&(lhmin>=m_lmin)&&(lhmin<=m_lmax)) + { + Mantid::Kernel::VMD v(m_hmin,khmin,lhmin,momhMin); + intersections.push_back(v); + } + } + double momhMax = fmom*(m_hmax-hStart)+m_kfmin; + if ((momhMax-m_kfmin)*(momhMax-m_kfmax)<=0) + { + //khmax and lhmax + double khmax = fk*(m_hmax-hStart)+kStart; + double lhmax = fl*(m_hmax-hStart)+lStart; + if((khmax>=m_kmin)&&(khmax<=m_kmax)&&(lhmax>=m_lmin)&&(lhmax<=m_lmax)) + { + Mantid::Kernel::VMD v(m_hmax,khmax,lhmax,momhMax); + intersections.push_back(v); + } + } + } + + //calculate intersections with planes perpendicular to k + if (fabs(kStart-kEnd) > eps) + { + double fmom=(m_kfmax-m_kfmin)/(kEnd-kStart); + double fh=(hEnd-hStart)/(kEnd-kStart); + double fl=(lEnd-lStart)/(kEnd-kStart); + if(!m_kIntegrated) + { + for(size_t i = 0;i=m_kmin)&&(ki<=m_kmax) && ((kStart-ki)*(kEnd-ki)<0)) + { + // if ki is between kStart and kEnd, then hi and li will be between hStart, hEnd and lStart, lEnd and momi will be between m_kfmin and m_kfmax + double hi = fh*(ki-kStart)+hStart; + double li = fl*(ki-kStart)+lStart; + if ((hi>=m_hmin)&&(hi<=m_hmax)&&(li>=m_lmin)&&(li<=m_lmax)) + { + double momi = fmom*(ki-kStart)+m_kfmin; + Mantid::Kernel::VMD v(hi,ki,li,momi); + intersections.push_back(v); + } + } + } + } + double momkMin = fmom*(m_kmin-kStart)+m_kfmin; + if ((momkMin-m_kfmin)*(momkMin-m_kfmax)<0) + { + //hkmin and lkmin + double hkmin = fh*(m_kmin-kStart)+hStart; + double lkmin = fl*(m_kmin-kStart)+lStart; + if((hkmin>=m_hmin)&&(hkmin<=m_hmax)&&(lkmin>=m_lmin)&&(lkmin<=m_lmax)) + { + Mantid::Kernel::VMD v(hkmin,m_kmin,lkmin,momkMin); + intersections.push_back(v); + } + } + double momkMax = fmom*(m_kmax-kStart)+m_kfmin; + if ((momkMax-m_kfmin)*(momkMax-m_kfmax)<=0) + { + //hkmax and lkmax + double hkmax = fh*(m_kmax-kStart)+hStart; + double lkmax = fl*(m_kmax-kStart)+lStart; + if((hkmax>=m_hmin)&&(hkmax<=m_hmax)&&(lkmax>=m_lmin)&&(lkmax<=m_lmax)) + { + Mantid::Kernel::VMD v(hkmax,m_kmax,lkmax,momkMax); + intersections.push_back(v); + } + } + } + + + //calculate intersections with planes perpendicular to l + if (fabs(lStart-lEnd) > eps) + { + double fmom=(m_kfmax-m_kfmin)/(lEnd-lStart); + double fh=(hEnd-hStart)/(lEnd-lStart); + double fk=(kEnd-kStart)/(lEnd-lStart); + if(!m_lIntegrated) + { + for(size_t i = 0;i=m_lmin)&&(li<=m_lmax) && ((lStart-li)*(lEnd-li)<0)) + { + double hi = fh*(li-lStart)+hStart; + double ki = fk*(li-lStart)+kStart; + if ((hi>=m_hmin)&&(hi<=m_hmax)&&(ki>=m_kmin)&&(ki<=m_kmax)) + { + double momi = fmom*(li-lStart)+m_kfmin; + Mantid::Kernel::VMD v(hi,ki,li,momi); + intersections.push_back(v); + } + } + } + } + double momlMin = fmom*(m_lmin-lStart)+m_kfmin; + if ((momlMin-m_kfmin)*(momlMin-m_kfmax)<=0) + { + //hlmin and klmin + double hlmin = fh*(m_lmin-lStart)+hStart; + double klmin = fk*(m_lmin-lStart)+kStart; + if((hlmin>=m_hmin)&&(hlmin<=m_hmax)&&(klmin>=m_kmin)&&(klmin<=m_kmax)) + { + Mantid::Kernel::VMD v(hlmin,klmin,m_lmin,momlMin); + intersections.push_back(v); + } + } + double momlMax = fmom*(m_lmax-lStart)+m_kfmin; + if ((momlMax-m_kfmin)*(momlMax-m_kfmax)<0) + { + //hlmax and klmax + double hlmax = fh*(m_lmax-lStart)+hStart; + double klmax = fk*(m_lmax-lStart)+kStart; + if((hlmax>=m_hmin)&&(hlmax<=m_hmax)&&(klmax>=m_kmin)&&(klmax<=m_kmax)) + { + Mantid::Kernel::VMD v(hlmax,klmax,m_lmax,momlMax); + intersections.push_back(v); + } + } + } + + //intersections with dE + if(!m_dEIntegrated) + { + for(size_t i = 0;i=m_hmin)&&(h<=m_hmax)&&(k>=m_kmin)&&(k<=m_kmax)&&(l>=m_lmin)&&(l<=m_lmax)) + { + Mantid::Kernel::VMD v(h,k,l,kfi); + intersections.push_back(v); + } + } + } + } + + //endpoints + if ((hStart>=m_hmin)&&(hStart<=m_hmax)&&(kStart>=m_kmin)&&(kStart<=m_kmax)&&(lStart>=m_lmin)&&(lStart<=m_lmax)) + { + Mantid::Kernel::VMD v(hStart,kStart,lStart,m_kfmin); + intersections.push_back(v); + } + if ((hEnd>=m_hmin)&&(hEnd<=m_hmax)&&(kEnd>=m_kmin)&&(kEnd<=m_kmax)&&(lEnd>=m_lmin)&&(lEnd<=m_lmax)) + { + Mantid::Kernel::VMD v(hEnd,kEnd,lEnd,m_kfmax); + intersections.push_back(v); + } + + //sort intersections by final momentum + typedef std::vector::iterator IterType; + std::stable_sort(intersections.begin(),intersections.end(),compareMomentum); + + + return intersections; + } + + } // namespace MDAlgorithms +} // namespace Mantid diff --git a/Code/Mantid/Framework/MDAlgorithms/test/MDNormDirectSCTest.h b/Code/Mantid/Framework/MDAlgorithms/test/MDNormDirectSCTest.h new file mode 100644 index 000000000000..6e6dae2b3e8c --- /dev/null +++ b/Code/Mantid/Framework/MDAlgorithms/test/MDNormDirectSCTest.h @@ -0,0 +1,119 @@ +#ifndef MANTID_MDALGORITHMS_MDNORMDIRECTSCTEST_H_ +#define MANTID_MDALGORITHMS_MDNORMDIRECTSCTEST_H_ + +#include + +#include "MantidMDAlgorithms/MDNormDirectSC.h" +#include "MantidMDAlgorithms/CreateMDWorkspace.h" +#include "MantidAPI/IMDHistoWorkspace.h" +#include "MantidAPI/WorkspaceFactory.h" +#include "MantidTestHelpers/WorkspaceCreationHelper.h" + +using Mantid::MDAlgorithms::MDNormDirectSC; +using namespace Mantid::API; + +class MDNormDirectSCTest : public CxxTest::TestSuite +{ +public: + // This pair of boilerplate methods prevent the suite being created statically + // This means the constructor isn't called when running other tests + static MDNormDirectSCTest *createSuite() { return new MDNormDirectSCTest(); } + static void destroySuite( MDNormDirectSCTest *suite ) { delete suite; } + + + void test_Init() + { + MDNormDirectSC alg; + TS_ASSERT_THROWS_NOTHING( alg.initialize() ) + TS_ASSERT( alg.isInitialized() ) + } + + void test_properties() + { + std::string mdWsName = "__temp_InputMDWorkspaceName"; + createMDWorkspace(mdWsName); + std::string fluxGoodWsName = "__temp_InputGoodFluxWorkspaceName"; + createGoodFluxWorkspace(fluxGoodWsName); + std::string fluxBadWsName = "__temp_InputBadFluxWorkspaceName"; + createBadFluxWorkspace(fluxBadWsName); + std::string saWsName = "__temp_InputSAWorkspaceName"; + createBadFluxWorkspace(saWsName); + + MDNormDirectSC alg; + TS_ASSERT_THROWS_NOTHING( alg.initialize() ) + TS_ASSERT( alg.isInitialized() ) + TS_ASSERT_THROWS_NOTHING( alg.setPropertyValue("InputWorkspace", mdWsName) ); + TS_ASSERT_THROWS_NOTHING( alg.setProperty("FluxWorkspace", fluxGoodWsName) ); + TS_ASSERT_THROWS_NOTHING( alg.setProperty("FluxWorkspace", fluxBadWsName) ); // it isn't bad any more + TS_ASSERT_THROWS_NOTHING( alg.setProperty("SolidAngleWorkspace", saWsName) ); + TS_ASSERT_THROWS_NOTHING( alg.setPropertyValue("OutputWorkspace", "OutWSName") ); + TS_ASSERT_THROWS_NOTHING( alg.setPropertyValue("OutputNormalizationWorkspace", "OutNormWSName") ); + + AnalysisDataService::Instance().clear(); + } + +private: + + void createMDWorkspace(const std::string& wsName) + { + const int ndims = 2; + std::string bins = "2,2"; + std::string extents = "0,1,0,1"; + std::vector names(ndims); + names[0] = "A"; names[1] = "B"; + std::vector units(ndims); + units[0] = "a"; units[1] = "b"; + + Mantid::MDAlgorithms::CreateMDWorkspace alg; + alg.initialize(); + TS_ASSERT_THROWS_NOTHING( alg.setProperty("Dimensions", ndims) ); + TS_ASSERT_THROWS_NOTHING( alg.setPropertyValue("Extents", extents) ); + TS_ASSERT_THROWS_NOTHING( alg.setProperty("Names", names) ); + TS_ASSERT_THROWS_NOTHING( alg.setProperty("Units", units) ); + TS_ASSERT_THROWS_NOTHING( alg.setPropertyValue("OutputWorkspace", wsName) ); + alg.execute(); + + } + + void createGoodFluxWorkspace(const std::string& wsName) + { + auto flux = WorkspaceCreationHelper::create2DWorkspaceWithFullInstrument( 2, 10 ); + auto &x = flux->dataX(0); + auto &y1 = flux->dataY(1); + + for(size_t i = 0; i < y1.size(); ++i) + { + y1[i] = 2 * x[i]; + } + flux->setX(1,x); + flux->getAxis(0)->setUnit("Momentum"); + + AnalysisDataService::Instance().addOrReplace( wsName, flux ); + } + + void createBadFluxWorkspace(const std::string& wsName) + { + auto flux = WorkspaceCreationHelper::create2DWorkspaceWithFullInstrument( 2, 10 ); + auto &x = flux->dataX(0); + auto &y1 = flux->dataY(1); + + for(size_t i = 0; i < y1.size(); ++i) + { + y1[i] = -2 * x[i]; + } + flux->setX(1,x); + flux->getAxis(0)->setUnit("Momentum"); + + AnalysisDataService::Instance().addOrReplace( wsName, flux ); + } + + void createSolidAngleWorkspace(const std::string& wsName) + { + auto sa = WorkspaceCreationHelper::create2DWorkspaceWithFullInstrument( 2, 10 ); + AnalysisDataService::Instance().addOrReplace( wsName, sa ); + } + +}; + + +#endif /* MANTID_MDALGORITHMS_MDNORMDIRECTSCTEST_H_ */