diff --git a/Code/Mantid/Framework/Algorithms/CMakeLists.txt b/Code/Mantid/Framework/Algorithms/CMakeLists.txt index 1ca9d718127b410b64c47906db3b5c771355f93d..0ff5aa144fa66e474ddd7ed38a72d2382f8aa93b 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 0000000000000000000000000000000000000000..fdfca91cc8a2293a0b4790303c3a79b3b7bafd76 --- /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 0000000000000000000000000000000000000000..e58e67ebe457760eca11af4b2d0df23c4a188926 --- /dev/null +++ b/Code/Mantid/Framework/Algorithms/src/LorentzCorrection.cpp @@ -0,0 +1,150 @@ +#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"); + + 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& inX = inWS->readX(i); + + MantidVec& outY = outWS->dataY(i); + MantidVec& outE = outWS->dataE(i); + + IDetector_const_sptr detector; + try + { + detector = inWS->getDetector(i); + } catch (Exception::NotFoundError&) + { + // Catch if no detector. Next line tests whether this happened - test placed + // outside here because Mac Intel compiler doesn't like 'continue' in a catch + // in an openmp block. + } + // If no detector found, skip onto the next spectrum + if (!detector) + continue; + + const double twoTheta = inWS->detectorTwoTheta(detector); + const double sinTheta = std::sin(twoTheta / 2); + double sinThetaSq = sinTheta * sinTheta; + + for (size_t j = 0; j < inY.size(); ++j) + { + const double wL = isHist ? (0.5 * (inX[j] + inX[j + 1])) : inX[j]; + if(wL == 0) + { + std::stringstream buffer; + buffer << "Cannot have zero values Wavelength. At workspace index: " << i; + throw std::runtime_error(buffer.str()); + } + + double weight = sinThetaSq / (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 diff --git a/Code/Mantid/Framework/Algorithms/test/LorentzCorrectionTest.h b/Code/Mantid/Framework/Algorithms/test/LorentzCorrectionTest.h new file mode 100644 index 0000000000000000000000000000000000000000..bd553d7cea50cfdc9c2aba77cace6cdf0176527e --- /dev/null +++ b/Code/Mantid/Framework/Algorithms/test/LorentzCorrectionTest.h @@ -0,0 +1,164 @@ +#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_throws_if_wavelength_zero() + { + auto ws_lam = this->create_workspace(2 /*nBins*/); + ws_lam->dataX(0)[0] = 0; // Make wavelength zero + ws_lam->dataX(0)[1] = 0; // Make wavelength zero + LorentzCorrection alg; + alg.setChild(true); + alg.setRethrows(true); + alg.initialize(); + alg.setProperty("InputWorkspace", ws_lam); + alg.setPropertyValue("OutputWorkspace", "temp"); + TSM_ASSERT_THROWS("Should throw with zero wavelength values.", alg.execute(), std::runtime_error&); + } + + 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/AnvredCorrection-v1.rst b/Code/Mantid/docs/source/algorithms/AnvredCorrection-v1.rst index e489e5afe7d1c2c49d42e034308155d56a5dd410..a3fb3e4024d227aa9e08f4ba8782468907f49ef3 100644 --- a/Code/Mantid/docs/source/algorithms/AnvredCorrection-v1.rst +++ b/Code/Mantid/docs/source/algorithms/AnvredCorrection-v1.rst @@ -42,6 +42,8 @@ The pixel efficiency and incident spectrum correction are NOT CURRENTLY USED. The absorption correction, trans, depends on both lamda and the pixel, Which is a fairly expensive calulation when done for each event. +Also see :ref:`algm-LorentzCorrection` + Usage ----- 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 0000000000000000000000000000000000000000..8022a4b00d231a7efaba3e94ce7a5a1f35340c10 --- /dev/null +++ b/Code/Mantid/docs/source/algorithms/LorentzCorrection-v1.rst @@ -0,0 +1,48 @@ + +.. algorithm:: + +.. summary:: + +.. alias:: + +.. properties:: + +Description +----------- + +Calculates and applies the lorentz correction weights to a workspace. The Lorentz correction *L* is calculated according to: + +.. math:: + L = \sin(\theta)^{2}/\lambda^{4} + +Where :math:`\theta` is the scattering angle. + +The calculations performed in this Algorithm are a subset of those performed by the :ref:`algm-AnvredCorrection` + +Usage +----- + +**Example - LorentzCorrection** + +.. testcode:: LorentzCorrectionExample + + tof = Load(Filename='HRP39180.RAW') + lam = ConvertUnits(InputWorkspace=tof, Target='Wavelength') + corrected = LorentzCorrection(InputWorkspace=lam) + + y = corrected.readY(2) + e = corrected.readE(2) + # print first corrected yvalues + print y[1:5] + # print first corrected evalues + print e[1:5] + +Output: + +.. testoutput:: LorentzCorrectionExample + + [ 0.84604876 0.4213364 1.67862035 0. ] + [ 0.59824681 0.4213364 0.83931018 0. ] + +.. categories:: +