From 69ee6ac724f0cb89dbe5fd9493e17934b8fc219f Mon Sep 17 00:00:00 2001 From: Owen Arnold <owen.arnold@stfc.ac.uk> Date: Sat, 13 Sep 2014 09:07:30 +0100 Subject: [PATCH] refs #10208. Add new algorithm --- .../Framework/Algorithms/CMakeLists.txt | 3 + .../inc/MantidAlgorithms/LorentzCorrection.h | 56 +++++++ .../Algorithms/src/LorentzCorrection.cpp | 124 +++++++++++++++ .../Algorithms/test/LorentzCorrectionTest.h | 150 ++++++++++++++++++ .../algorithms/LorentzCorrection-v1.rst | 44 +++++ 5 files changed, 377 insertions(+) create mode 100644 Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/LorentzCorrection.h create mode 100644 Code/Mantid/Framework/Algorithms/src/LorentzCorrection.cpp create mode 100644 Code/Mantid/Framework/Algorithms/test/LorentzCorrectionTest.h create mode 100644 Code/Mantid/docs/source/algorithms/LorentzCorrection-v1.rst diff --git a/Code/Mantid/Framework/Algorithms/CMakeLists.txt b/Code/Mantid/Framework/Algorithms/CMakeLists.txt index 1ca9d718127..0ff5aa144fa 100644 --- a/Code/Mantid/Framework/Algorithms/CMakeLists.txt +++ b/Code/Mantid/Framework/Algorithms/CMakeLists.txt @@ -134,6 +134,7 @@ set ( SRC_FILES src/InterpolatingRebin.cpp src/InvertMask.cpp src/Logarithm.cpp + src/LorentzCorrection.cpp src/MaskBins.cpp src/MaskBinsFromTable.cpp src/MaskDetectorsIf.cpp @@ -379,6 +380,7 @@ set ( INC_FILES inc/MantidAlgorithms/InterpolatingRebin.h inc/MantidAlgorithms/InvertMask.h inc/MantidAlgorithms/Logarithm.h + inc/MantidAlgorithms/LorentzCorrection.h inc/MantidAlgorithms/MaskBins.h inc/MantidAlgorithms/MaskBinsFromTable.h inc/MantidAlgorithms/MaskDetectorsIf.h @@ -625,6 +627,7 @@ set ( TEST_FILES InterpolatingRebinTest.h InvertMaskTest.h LogarithmTest.h + LorentzCorrectionTest.h MaskBinsFromTableTest.h MaskBinsTest.h MaxMinTest.h diff --git a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/LorentzCorrection.h b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/LorentzCorrection.h new file mode 100644 index 00000000000..fdfca91cc8a --- /dev/null +++ b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/LorentzCorrection.h @@ -0,0 +1,56 @@ +#ifndef MANTID_ALGORITHMS_LORENTZCORRECTION_H_ +#define MANTID_ALGORITHMS_LORENTZCORRECTION_H_ + +#include "MantidKernel/System.h" +#include "MantidAPI/Algorithm.h" + +namespace Mantid +{ +namespace Algorithms +{ + + /** LorentzCorrection : Algorithm Performs a lorentz correction (sin(theta)^2)/(wavelength^4) on a MatrixWorkspace in units of wavelength + + Copyright © 2014 ISIS Rutherford Appleton Laboratory & NScD Oak Ridge National Laboratory + + 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 <http://www.gnu.org/licenses/>. + + File change history is stored at: <https://github.com/mantidproject/mantid> + Code Documentation is available at: <http://doxygen.mantidproject.org> + */ + class DLLExport LorentzCorrection : public API::Algorithm + { + public: + LorentzCorrection(); + virtual ~LorentzCorrection(); + + 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(); + + + }; + + +} // namespace Algorithms +} // namespace Mantid + +#endif /* MANTID_ALGORITHMS_LORENTZCORRECTION_H_ */ \ No newline at end of file diff --git a/Code/Mantid/Framework/Algorithms/src/LorentzCorrection.cpp b/Code/Mantid/Framework/Algorithms/src/LorentzCorrection.cpp new file mode 100644 index 00000000000..0517fcc4e16 --- /dev/null +++ b/Code/Mantid/Framework/Algorithms/src/LorentzCorrection.cpp @@ -0,0 +1,124 @@ +#include "MantidAlgorithms/LorentzCorrection.h" +#include "MantidAPI/MatrixWorkspace.h" +#include "MantidAPI/Progress.h" +#include "MantidAPI/WorkspaceValidators.h" +#include "MantidKernel/MultiThreaded.h" +#include <boost/make_shared.hpp> +#include <boost/shared_ptr.hpp> +#include <cmath> + +namespace Mantid +{ +namespace Algorithms +{ + using namespace Mantid::Kernel; + using namespace Mantid::API; + using namespace Mantid::Geometry; + using namespace Mantid::Kernel; + using Mantid::API::WorkspaceProperty; + + // Register the algorithm into the AlgorithmFactory + DECLARE_ALGORITHM(LorentzCorrection) + + + + //---------------------------------------------------------------------------------------------- + /** Constructor + */ + LorentzCorrection::LorentzCorrection() + { + } + + //---------------------------------------------------------------------------------------------- + /** Destructor + */ + LorentzCorrection::~LorentzCorrection() + { + } + + + //---------------------------------------------------------------------------------------------- + + + /// Algorithm's version for identification. @see Algorithm::version + int LorentzCorrection::version() const { return 1;}; + + /// Algorithm's category for identification. @see Algorithm::category + const std::string LorentzCorrection::category() const { return "Crystal";} + + /// Algorithm's summary for use in the GUI and help. @see Algorithm::summary + const std::string LorentzCorrection::summary() const { return "Performs a white beam Lorentz Correction";}; + + const std::string LorentzCorrection::name() const {return "LorentzCorrection";} + + //---------------------------------------------------------------------------------------------- + /** Initialize the algorithm's properties. + */ + void LorentzCorrection::init() + { + + declareProperty(new WorkspaceProperty<MatrixWorkspace>("InputWorkspace","",Direction::Input, PropertyMode::Mandatory, boost::make_shared<WorkspaceUnitValidator>("Wavelength")), "Input workspace to correct in Wavelength."); + declareProperty(new WorkspaceProperty<MatrixWorkspace>("OutputWorkspace","",Direction::Output), "An output workspace."); + } + + //---------------------------------------------------------------------------------------------- + /** Execute the algorithm. + */ + void LorentzCorrection::exec() + { + MatrixWorkspace_sptr inWS = this->getProperty("InputWorkspace"); + const std::string inWSName = this->getPropertyValue("InputWorkspace"); + const std::string outWSName = this->getPropertyValue("OutputWorkspace"); + + auto cloneAlg = this->createChildAlgorithm("CloneWorkspace", 0, 0.1); + cloneAlg->initialize(); + cloneAlg->setProperty("InputWorkspace", inWS); + cloneAlg->execute(); + Workspace_sptr temp = cloneAlg->getProperty("OutputWorkspace"); + MatrixWorkspace_sptr outWS = boost::dynamic_pointer_cast<MatrixWorkspace>(temp); + + const auto numHistos = inWS->getNumberHistograms(); + Progress prog(this, 0, 1, numHistos); + const bool isHist = inWS->isHistogramData(); + + PARALLEL_FOR1(inWS) + for (int64_t i = 0; i < int64_t(numHistos); ++i) + { + PARALLEL_START_INTERUPT_REGION + + + const MantidVec& inY = inWS->readY(i); + const MantidVec& inE = inWS->readE(i); + const MantidVec& inX = inWS->readX(i); + + MantidVec& outY = outWS->dataY(i); + MantidVec& outE = outWS->dataE(i); + + const auto & detector = inWS->getDetector(i); + + const double twoTheta = inWS->detectorTwoTheta(detector); + + for (int64_t j = 0; j < inY.size(); ++j) + { + const double wL = isHist ? (0.5 * (inX[j] + inX[j + 1])) : inX[j]; + double sinTheta = std::sin( twoTheta/2 ); + double weight = sinTheta * sinTheta / (wL * wL * wL * wL); + outY[j] *= weight; + outE[j] *= weight; + } + + prog.report(); + + PARALLEL_END_INTERUPT_REGION + } + PARALLEL_CHECK_INTERUPT_REGION + + + + this->setProperty("OutputWorkspace", outWS); + } + + + +} // namespace Algorithms +} // namespace Mantid \ No newline at end of file diff --git a/Code/Mantid/Framework/Algorithms/test/LorentzCorrectionTest.h b/Code/Mantid/Framework/Algorithms/test/LorentzCorrectionTest.h new file mode 100644 index 00000000000..1a5fdafa23b --- /dev/null +++ b/Code/Mantid/Framework/Algorithms/test/LorentzCorrectionTest.h @@ -0,0 +1,150 @@ +#ifndef MANTID_ALGORITHMS_LORENTZCORRECTIONTEST_H_ +#define MANTID_ALGORITHMS_LORENTZCORRECTIONTEST_H_ + +#include <cxxtest/TestSuite.h> + +#include "MantidAlgorithms/LorentzCorrection.h" +#include "MantidTestHelpers/WorkspaceCreationHelper.h" +#include "MantidDataObjects/Workspace2D.h" +#include "MantidKernel/V3D.h" +#include "MantidGeometry/Instrument.h" +#include "MantidGeometry/Instrument/ReferenceFrame.h" +#include "MantidGeometry/Instrument/Detector.h" +#include "MantidGeometry/Instrument/ObjComponent.h" +#include <cmath> + +using Mantid::Algorithms::LorentzCorrection; +using namespace Mantid::Kernel; +using namespace Mantid::API; +using namespace Mantid::Geometry; +using namespace WorkspaceCreationHelper; + +class LorentzCorrectionTest: public CxxTest::TestSuite +{ + +private: + + /* + * Calculate what the weight should be. + */ + double calculate_weight_at(MatrixWorkspace_sptr& ws, const int bin_index) + { + const Mantid::MantidVec& xData = ws->readX(0); + + auto detector = ws->getDetector(0); + double twotheta = ws->detectorTwoTheta(detector); + double lam = (xData[bin_index] + xData[bin_index + 1]) / 2; + double weight = std::sin(twotheta / 2); + weight = weight * weight; + weight = weight / (lam * lam * lam * lam); + return weight; + } + + /** + Create a workspace in wavelength with a simple instrument defined with a single detector. + */ + MatrixWorkspace_sptr create_workspace(const int nBins) + { + Instrument_sptr instrument = boost::make_shared<Instrument>(); + instrument->setReferenceFrame(boost::make_shared<ReferenceFrame>(Y, X, Left, "0,0,0")); + + ObjComponent *source = new ObjComponent("source"); + source->setPos(V3D(0, 0, 0)); + instrument->add(source); + instrument->markAsSource(source); + + ObjComponent *sample = new ObjComponent("some-surface-holder"); + source->setPos(V3D(15, 0, 0)); + instrument->add(sample); + instrument->markAsSamplePos(sample); + + Detector* det = new Detector("my-detector", 1, NULL); + det->setPos(20, (20 - sample->getPos().X()), 0); + instrument->add(det); + instrument->markAsDetector(det); + + const int nSpectra = 1; + const double deltaX = 10; + const double startX = 0; + auto workspace = Create2DWorkspaceBinned(nSpectra, nBins, startX, deltaX); // Creates histogram data + + auto & ydata = workspace->dataY(0); + auto & edata = workspace->dataE(0); + for (size_t i = 0; i < ydata.size(); ++i) + { + ydata[i] = 1; + edata[i] = 1; + } + + workspace->getAxis(0)->setUnit("Wavelength"); + workspace->setYUnit("Counts"); + workspace->setInstrument(instrument); + workspace->getSpectrum(0)->addDetectorID(det->getID()); + return workspace; + } + +public: + // This pair of boilerplate methods prevent the suite being created statically + // This means the constructor isn't called when running other tests + static LorentzCorrectionTest *createSuite() + { + return new LorentzCorrectionTest(); + } + static void destroySuite(LorentzCorrectionTest *suite) + { + delete suite; + } + + void test_init() + { + LorentzCorrection alg; + TS_ASSERT_THROWS_NOTHING( alg.initialize()) + TS_ASSERT( alg.isInitialized()) + } + + void test_check_input_units() + { + const int nHisto = 1; + const int nBins = 1; + auto ws_tof = create2DWorkspaceWithFullInstrument(nHisto, nBins); + + LorentzCorrection alg; + alg.setChild(true); + alg.initialize(); + TSM_ASSERT_THROWS("Workspace must be in units of wavelength", + alg.setProperty("InputWorkspace", ws_tof), std::invalid_argument&); + } + + void test_execute() + { + auto ws_lam = this->create_workspace(2 /*nBins*/); + + LorentzCorrection alg; + alg.initialize(); + alg.setChild(true); + alg.setProperty("InputWorkspace", ws_lam); + alg.setPropertyValue("OutputWorkspace", "temp"); + alg.execute(); + MatrixWorkspace_sptr out_ws = alg.getProperty("OutputWorkspace"); + TS_ASSERT(out_ws != NULL); + + const std::string unitID = out_ws->getAxis(0)->unit()->unitID(); + TS_ASSERT_EQUALS(unitID, "Wavelength"); + + const Mantid::MantidVec& yData = out_ws->readY(0); + const Mantid::MantidVec& eData = out_ws->readE(0); + + int index = 0; + double weight = calculate_weight_at(out_ws, index /*y index*/); + TS_ASSERT_EQUALS(yData[index], weight); + TS_ASSERT_EQUALS(eData[index], weight); + + index++; // go to 1 + weight = calculate_weight_at(out_ws, index /*y index*/); + TS_ASSERT_EQUALS(yData[index], weight); + TS_ASSERT_EQUALS(eData[index], weight); + } + +}; + +#endif /* MANTID_ALGORITHMS_LORENTZCORRECTIONTEST_H_ */ diff --git a/Code/Mantid/docs/source/algorithms/LorentzCorrection-v1.rst b/Code/Mantid/docs/source/algorithms/LorentzCorrection-v1.rst new file mode 100644 index 00000000000..fff2112bf5e --- /dev/null +++ b/Code/Mantid/docs/source/algorithms/LorentzCorrection-v1.rst @@ -0,0 +1,44 @@ + +.. algorithm:: + +.. summary:: + +.. alias:: + +.. properties:: + +Description +----------- + +Performs a lorentz correction on a white beam workspace in units of wavelength. + + L = sin\theta^{2}/\lam^{4} + + +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 - LorentzCorrection** + +.. testcode:: LorentzCorrectionExample + + ws = CreateSampleWorkspace("Histogram",NumBanks=1,BankPixelWidth=1) + wavelength = ConvertUnits(ws, Target="Wavelength") + + corrected = LorentzCorrection(InputWorkspace=wavelength) + + # Print some of the weights. Constant theta angle. + print "The output workspace has %i spectra" % wsOut.getNumberHistograms() + +Output: + +.. testoutput:: LorentzCorrectionExample + + The output workspace has ?? spectra + +.. categories:: + -- GitLab