From 3ce66bb52273ba4311046c2ee31137147badf75e Mon Sep 17 00:00:00 2001 From: Dan Nixon Date: Tue, 26 May 2015 14:59:42 +0100 Subject: [PATCH] Initial version of resolution algorithm from Fortran Refs #11833 --- .../Framework/Algorithms/CMakeLists.txt | 3 + .../VesuvioL1ThetaResolution.h | 67 ++++ .../src/VesuvioL1ThetaResolution.cpp | 306 ++++++++++++++++++ .../test/VesuvioL1ThetaResolutionTest.h | 60 ++++ .../VesuvioL1ThetaResolution-v1.rst | 44 +++ 5 files changed, 480 insertions(+) create mode 100644 Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/VesuvioL1ThetaResolution.h create mode 100644 Code/Mantid/Framework/Algorithms/src/VesuvioL1ThetaResolution.cpp create mode 100644 Code/Mantid/Framework/Algorithms/test/VesuvioL1ThetaResolutionTest.h create mode 100644 Code/Mantid/docs/source/algorithms/VesuvioL1ThetaResolution-v1.rst diff --git a/Code/Mantid/Framework/Algorithms/CMakeLists.txt b/Code/Mantid/Framework/Algorithms/CMakeLists.txt index 012f793358d5..99e02a77ca29 100644 --- a/Code/Mantid/Framework/Algorithms/CMakeLists.txt +++ b/Code/Mantid/Framework/Algorithms/CMakeLists.txt @@ -250,6 +250,7 @@ set ( SRC_FILES src/UnwrapMonitor.cpp src/UnwrapSNS.cpp src/UpdateScriptRepository.cpp + src/VesuvioL1ThetaResolution.cpp src/WeightedMean.cpp src/WeightedMeanOfWorkspace.cpp src/WeightingStrategy.cpp @@ -511,6 +512,7 @@ set ( INC_FILES inc/MantidAlgorithms/UnwrapMonitor.h inc/MantidAlgorithms/UnwrapSNS.h inc/MantidAlgorithms/UpdateScriptRepository.h + inc/MantidAlgorithms/VesuvioL1ThetaResolution.h inc/MantidAlgorithms/WeightedMean.h inc/MantidAlgorithms/WeightedMeanOfWorkspace.h inc/MantidAlgorithms/WeightingStrategy.h @@ -757,6 +759,7 @@ set ( TEST_FILES UnGroupWorkspaceTest.h UnaryOperationTest.h UnwrapSNSTest.h + VesuvioL1ThetaResolutionTest.h WeightedMeanOfWorkspaceTest.h WeightedMeanTest.h WeightingStrategyTest.h diff --git a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/VesuvioL1ThetaResolution.h b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/VesuvioL1ThetaResolution.h new file mode 100644 index 000000000000..e88d2da10e9e --- /dev/null +++ b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/VesuvioL1ThetaResolution.h @@ -0,0 +1,67 @@ +#ifndef MANTID_ALGORITHMS_VESUVIOL1THETARESOLUTION_H_ +#define MANTID_ALGORITHMS_VESUVIOL1THETARESOLUTION_H_ + +#include "MantidKernel/System.h" +#include "MantidAPI/Algorithm.h" + +#include + +namespace Mantid { +namespace Algorithms { + +/** VesuvioL1ThetaResolution + + Calculates the resolution function for L1 and Theta. + + Copyright © 2015 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 VesuvioL1ThetaResolution : public API::Algorithm { +public: + VesuvioL1ThetaResolution(); + virtual ~VesuvioL1ThetaResolution(); + + 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 loadInstrument(); + + void calculateDetector(Mantid::Geometry::IDetector_const_sptr detector, std::vector& l1Values, std::vector& thetaValues); + Mantid::API::MatrixWorkspace_sptr processDistribution(Mantid::API::MatrixWorkspace_sptr ws, const double binWidth); + double random(); + + Mantid::API::MatrixWorkspace_sptr m_instWorkspace; + Mantid::Geometry::IComponent_const_sptr m_sample; + Mantid::API::MatrixWorkspace_sptr m_l1DistributionWs; + Mantid::API::MatrixWorkspace_sptr m_thetaDistributionWs; + + boost::mt19937 m_generator; +}; + +} // namespace Algorithms +} // namespace Mantid + +#endif /* MANTID_ALGORITHMS_VESUVIOL1THETARESOLUTION_H_ */ diff --git a/Code/Mantid/Framework/Algorithms/src/VesuvioL1ThetaResolution.cpp b/Code/Mantid/Framework/Algorithms/src/VesuvioL1ThetaResolution.cpp new file mode 100644 index 000000000000..9553c33119d7 --- /dev/null +++ b/Code/Mantid/Framework/Algorithms/src/VesuvioL1ThetaResolution.cpp @@ -0,0 +1,306 @@ +#include "MantidAlgorithms/VesuvioL1ThetaResolution.h" + +#include "MantidAPI/AlgorithmManager.h" +#include "MantidKernel/BoundedValidator.h" +#include "MantidKernel/Statistics.h" + +#include +#include +#include + +namespace Mantid { +namespace Algorithms { + +using namespace Mantid::Kernel; +using namespace Mantid::API; +using namespace Mantid::Geometry; + +namespace { + Mantid::Kernel::Logger g_log("VesuvioL1ThetaResolution"); +} + +// Register the algorithm into the AlgorithmFactory +DECLARE_ALGORITHM(VesuvioL1ThetaResolution) + +//---------------------------------------------------------------------------------------------- +/** Constructor + */ +VesuvioL1ThetaResolution::VesuvioL1ThetaResolution() {} + +//---------------------------------------------------------------------------------------------- +/** Destructor + */ +VesuvioL1ThetaResolution::~VesuvioL1ThetaResolution() {} + +//---------------------------------------------------------------------------------------------- + +/// Algorithms name for identification. @see Algorithm::name +const std::string VesuvioL1ThetaResolution::name() const { return "VesuvioL1ThetaResolution"; } + +/// Algorithm's version for identification. @see Algorithm::version +int VesuvioL1ThetaResolution::version() const { return 1; } + +/// Algorithm's category for identification. @see Algorithm::category +const std::string VesuvioL1ThetaResolution::category() const { + return "CorrectionFunctions"; +} + +/// Algorithm's summary for use in the GUI and help. @see Algorithm::summary +const std::string VesuvioL1ThetaResolution::summary() const { + return "Calculates resolution of l1 and theta"; +} + +//---------------------------------------------------------------------------------------------- +/** Initialize the algorithm's properties. + */ +void VesuvioL1ThetaResolution::init() { + auto positiveInt = boost::make_shared>(); + positiveInt->setLower(1); + + declareProperty("SpectrumMin", 3, + "Index of minimum spectrum"); + declareProperty("SpectrumMax", 198, + "Index of maximum spectrum"); + + declareProperty("NumEvents", 10000, + "Number of scattering events"); + declareProperty("Seed", 123456789, + "Seed for random number generator"); + + declareProperty("L1BinWidth", 0.01, + "Bin width for L1 distribution."); + declareProperty("ThetaBinWidth", 0.0001, + "Bin width for theta distribution."); + + declareProperty( + new WorkspaceProperty<>("L1Distribution", "", Direction::Output, PropertyMode::Optional), + "Distribution of lengths of the final flight path."); + + declareProperty( + new WorkspaceProperty<>("ThetaDistribution", "", Direction::Output, PropertyMode::Optional), + "Distribution of scattering angles."); +} + +//---------------------------------------------------------------------------------------------- +/** Execute the algorithm. + */ +void VesuvioL1ThetaResolution::exec() { + // Set up random number generator + m_generator.seed(static_cast(static_cast(getProperty("Seed")))); + + // Load the instrument workspace + loadInstrument(); + + const std::string l1DistributionWsName = getPropertyValue("L1Distribution"); + const std::string thetaDistributionWsName = getPropertyValue("ThetaDistribution"); + const size_t numHist = m_instWorkspace->getNumberHistograms(); + const int numEvents = getProperty("NumEvents"); + + // Create output workspaces for distributions if required + if(!l1DistributionWsName.empty()) + m_l1DistributionWs = WorkspaceFactory::Instance().create("Workspace2D", numHist, numEvents, numEvents); + if(!thetaDistributionWsName.empty()) + m_thetaDistributionWs = WorkspaceFactory::Instance().create("Workspace2D", numHist, numEvents, numEvents); + + // Set up progress reporting + Progress prog(this, 0.0, 1.0, numHist); + + // Loop for all detectors + for(size_t i = 0; i < numHist; i++) { + std::vector l1; + std::vector theta; + IDetector_const_sptr det = m_instWorkspace->getDetector(i); + + // Report progress + std::stringstream report; + report << "Detector " << det->getID(); + prog.report(report.str()); + + // Do simulation + calculateDetector(det, l1, theta); + + // Calculate statistics for L1 and theta + Statistics l1Stats = getStatistics(l1); + Statistics thetaStats = getStatistics(theta); + + g_log.information() << "Detector ID: " << det->getID() << std::endl + << "l0: mean=" << l1Stats.mean << ", std.dev.=" + << l1Stats.standard_deviation << std::endl + << "theta: mean=" << thetaStats.mean << ", std.dev.=" + << thetaStats.standard_deviation << std::endl; + + // Process data for L1 distribution + if(m_l1DistributionWs) { + std::vector& x = m_l1DistributionWs->dataX(i); + std::vector y(numEvents, 1.0); + + std::sort(l1.begin(), l1.end()); + std::copy(l1.begin(), l1.end(), x.begin()); + + m_l1DistributionWs->dataY(i) = y; + } + + // Process data for theta distribution + if(m_thetaDistributionWs) { + std::vector& x = m_thetaDistributionWs->dataX(i); + std::vector y(numEvents, 1.0); + + std::sort(theta.begin(), theta.end()); + std::copy(theta.begin(), theta.end(), x.begin()); + + m_thetaDistributionWs->dataY(i) = y; + } + } + + // Process the L1 distribution workspace + if(m_l1DistributionWs) { + const double binWidth = getProperty("L1BinWidth"); + setProperty("L1Distribution", processDistribution(m_l1DistributionWs, binWidth)); + } + + // Process the theta distribution workspace + if(m_thetaDistributionWs) { + const double binWidth = getProperty("ThetaBinWidth"); + setProperty("ThetaDistribution", processDistribution(m_thetaDistributionWs, binWidth)); + } + +} + +//---------------------------------------------------------------------------------------------- +/** Loads the instrument into a workspace. + */ +void VesuvioL1ThetaResolution::loadInstrument() { + MatrixWorkspace_sptr tempWS = WorkspaceFactory::Instance().create("Workspace2D", 1, 1, 1); + const std::string vesuvioIPF = tempWS->getInstrumentFilename("VESUVIO"); + + IAlgorithm_sptr loadInst = AlgorithmManager::Instance().create("LoadEmptyInstrument"); + loadInst->initialize(); + loadInst->setChild(true); + loadInst->setLogging(false); + loadInst->setProperty("OutputWorkspace", "__evs"); + loadInst->setProperty("Filename", vesuvioIPF); + loadInst->execute(); + m_instWorkspace = loadInst->getProperty("OutputWorkspace"); + + //TODO: load par file + + const int specIdxMin = static_cast(m_instWorkspace->getIndexFromSpectrumNumber(getProperty("SpectrumMin"))); + const int specIdxMax = static_cast(m_instWorkspace->getIndexFromSpectrumNumber(getProperty("SpectrumMax"))); + + IAlgorithm_sptr crop = AlgorithmManager::Instance().create("CropWorkspace"); + crop->initialize(); + crop->setChild(true); + crop->setLogging(false); + crop->setProperty("InputWorkspace", m_instWorkspace); + crop->setProperty("OutputWorkspace", "__evs"); + crop->setProperty("StartWorkspaceIndex", specIdxMin); + crop->setProperty("EndWorkspaceIndex", specIdxMax); + crop->execute(); + m_instWorkspace = crop->getProperty("OutputWorkspace"); + + m_sample = m_instWorkspace->getInstrument()->getSample(); +} + +//---------------------------------------------------------------------------------------------- +/** Loads the instrument into a workspace. + */ +void VesuvioL1ThetaResolution::calculateDetector(IDetector_const_sptr detector, std::vector& l1Values, std::vector& thetaValues) { + const int numEvents = getProperty("NumEvents"); + l1Values.reserve(numEvents); + thetaValues.reserve(numEvents); + + //TODO + const double detHeight = 25.0; + const double detWidth = 2.5; + const double beamWidth = 3.0; + + // Scattering angle in rad + const double theta = m_instWorkspace->detectorSignedTwoTheta(detector); + if(theta == 0.0) + return; + + // Final flight path in cm + const double l1av = detector->getDistance(*m_sample) * 100.0; + + const double x0 = l1av * sin(theta); + const double y0 = l1av * cos(theta); + + // Get as many events as defined by NumEvents + // This loop is not iteration limited but highly unlikely to ever become infinate + while(l1Values.size() < static_cast(numEvents)) { + const double xs = -beamWidth/2 + beamWidth*random(); + const double ys = 0.0; + const double zs = -beamWidth/2 + beamWidth*random(); + const double rs = sqrt(pow(xs, 2) + pow(xs, 2)); + + if(rs <= beamWidth/2) { + const double a = -detWidth/2 + detWidth*random(); + const double xd = x0 - a*cos(theta); + const double yd = y0 + a*sin(theta); + const double zd = -detHeight/2 + detHeight*random(); + + const double l1 = sqrt(pow(xd-xs, 2) + pow(yd-ys, 2) + pow(zd-zs, 2)); + double angle = acos(yd / l1); + + if(xd < 0.0) + angle *= -1; + + //TODO: convert angle to degrees + + l1Values.push_back(l1); + thetaValues.push_back(angle); + } + + interruption_point(); + } +} + +//---------------------------------------------------------------------------------------------- +/** Rebins the distributions and sets error values. + */ +MatrixWorkspace_sptr VesuvioL1ThetaResolution::processDistribution(MatrixWorkspace_sptr ws, const double binWidth) { + const size_t numHist = ws->getNumberHistograms(); + + double xMin(DBL_MAX); + double xMax(DBL_MIN); + for(size_t i = 0; i < numHist; i++) { + const std::vector x = ws->readX(i); + if(x[0] < xMin) + xMin = x[0]; + if(x[x.size()-1] > xMax) + xMax = x[x.size()-1]; + } + + std::stringstream binParams; + binParams << xMin << "," << binWidth << "," << xMax; + + IAlgorithm_sptr rebin = AlgorithmManager::Instance().create("Rebin"); + rebin->initialize(); + rebin->setChild(true); + rebin->setLogging(false); + rebin->setProperty("InputWorkspace", ws); + rebin->setProperty("OutputWorkspace", "__rebin"); + rebin->setProperty("Params", binParams.str()); + rebin->execute(); + ws = rebin->getProperty("OutputWorkspace"); + + for(size_t i = 0; i < numHist; i++) { + const std::vector y = ws->readY(i); + std::vector& e = ws->dataE(i); + + std::transform(y.begin(), y.end(), e.begin(), sqrt); + } + + return ws; +} + +//---------------------------------------------------------------------------------------------- +/** Generates a random number. + */ +double VesuvioL1ThetaResolution::random() { + typedef boost::uniform_real uniform_double; + return boost::variate_generator(m_generator, uniform_double(0.0, 1.0))(); +} + +} // namespace Algorithms +} // namespace Mantid diff --git a/Code/Mantid/Framework/Algorithms/test/VesuvioL1ThetaResolutionTest.h b/Code/Mantid/Framework/Algorithms/test/VesuvioL1ThetaResolutionTest.h new file mode 100644 index 000000000000..055657534192 --- /dev/null +++ b/Code/Mantid/Framework/Algorithms/test/VesuvioL1ThetaResolutionTest.h @@ -0,0 +1,60 @@ +#ifndef MANTID_ALGORITHMS_VESUVIOL1THETARESOLUTIONTEST_H_ +#define MANTID_ALGORITHMS_VESUVIOL1THETARESOLUTIONTEST_H_ + +#include + +#include "MantidAlgorithms/VesuvioL1ThetaResolution.h" + +using Mantid::Algorithms::VesuvioL1ThetaResolution; +using namespace Mantid::API; + +class VesuvioL1ThetaResolutionTest : 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 VesuvioL1ThetaResolutionTest *createSuite() { return new VesuvioL1ThetaResolutionTest(); } + static void destroySuite( VesuvioL1ThetaResolutionTest *suite ) { delete suite; } + + + void test_Init() + { + VesuvioL1ThetaResolution alg; + TS_ASSERT_THROWS_NOTHING( alg.initialize() ) + TS_ASSERT( alg.isInitialized() ) + } + + void test_exec() + { + // Name of the output workspace. + std::string outWSName("VesuvioL1ThetaResolutionTest_OutputWS"); + + VesuvioL1ThetaResolution alg; + TS_ASSERT_THROWS_NOTHING( alg.initialize() ) + TS_ASSERT( alg.isInitialized() ) + TS_ASSERT_THROWS_NOTHING( alg.setPropertyValue("REPLACE_PROPERTY_NAME_HERE!!!!", "value") ); + TS_ASSERT_THROWS_NOTHING( alg.setPropertyValue("OutputWorkspace", outWSName) ); + TS_ASSERT_THROWS_NOTHING( alg.execute(); ); + TS_ASSERT( alg.isExecuted() ); + + // Retrieve the workspace from data service. TODO: Change to your desired type + Workspace_sptr ws; + TS_ASSERT_THROWS_NOTHING( ws = AnalysisDataService::Instance().retrieveWS(outWSName) ); + TS_ASSERT(ws); + if (!ws) return; + + // TODO: Check the results + + // Remove workspace from the data service. + AnalysisDataService::Instance().remove(outWSName); + } + + void test_Something() + { + TSM_ASSERT( "You forgot to write a test!", 0); + } + + +}; + + +#endif /* MANTID_ALGORITHMS_VESUVIOL1THETARESOLUTIONTEST_H_ */ \ No newline at end of file diff --git a/Code/Mantid/docs/source/algorithms/VesuvioL1ThetaResolution-v1.rst b/Code/Mantid/docs/source/algorithms/VesuvioL1ThetaResolution-v1.rst new file mode 100644 index 000000000000..e35a37172506 --- /dev/null +++ b/Code/Mantid/docs/source/algorithms/VesuvioL1ThetaResolution-v1.rst @@ -0,0 +1,44 @@ + +.. algorithm:: + +.. summary:: + +.. alias:: + +.. properties:: + +Description +----------- + +TODO: Enter a full rst-markup description of your algorithm here. + + +Usage +----- +.. Try not to use files in your examples, + but if you cannot avoid it then the (small) files must be added to + autotestdata\UsageData and the following tag unindented + .. include:: ../usagedata-note.txt + +**Example - VesuvioL1ThetaResolution** + +.. testcode:: VesuvioL1ThetaResolutionExample + + # Create a host workspace + ws = CreateWorkspace(DataX=range(0,3), DataY=(0,2)) + or + ws = CreateSampleWorkspace() + + wsOut = VesuvioL1ThetaResolution() + + # Print the result + print "The output workspace has %i spectra" % wsOut.getNumberHistograms() + +Output: + +.. testoutput:: VesuvioL1ThetaResolutionExample + + The output workspace has ?? spectra + +.. categories:: +