refs #10208. Add new algorithm

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

Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/LorentzCorrection.h

+#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

Code/Mantid/Framework/Algorithms/src/LorentzCorrection.cpp

+#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

Code/Mantid/Framework/Algorithms/test/LorentzCorrectionTest.h

+#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_ */

Code/Mantid/docs/source/algorithms/LorentzCorrection-v1.rst

+
+.. 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::
+