Merge pull request #13143 from mganeva/toftof_moncoreff

Toftof moncoreff

Code/Mantid/Framework/Algorithms/CMakeLists.txt

@@ -155,6 +155,7 @@ set ( SRC_FILES
 	src/Minus.cpp
 	src/ModeratorTzero.cpp
 	src/ModeratorTzeroLinear.cpp
+        src/MonitorEfficiencyCorUser.cpp
 	src/MonteCarloAbsorption.cpp
 	src/MultipleScatteringCylinderAbsorption.cpp
 	src/Multiply.cpp
@@ -426,6 +427,7 @@ set ( INC_FILES
 	inc/MantidAlgorithms/Minus.h
 	inc/MantidAlgorithms/ModeratorTzero.h
 	inc/MantidAlgorithms/ModeratorTzeroLinear.h
+        inc/MantidAlgorithms/MonitorEfficiencyCorUser.h
 	inc/MantidAlgorithms/MonteCarloAbsorption.h
 	inc/MantidAlgorithms/MultipleScatteringCylinderAbsorption.h
 	inc/MantidAlgorithms/Multiply.h
@@ -697,6 +699,7 @@ set ( TEST_FILES
 	MinusTest.h
 	ModeratorTzeroLinearTest.h
 	ModeratorTzeroTest.h
+        MonitorEfficiencyCorUserTest.h
 	MonteCarloAbsorptionTest.h
 	MultipleScatteringCylinderAbsorptionTest.h
 	MultiplyRangeTest.h

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

+#ifndef MANTID_ALGORITHMS_MONITOREFFICIENCYCORUSER_H_
+#define MANTID_ALGORITHMS_MONITOREFFICIENCYCORUSER_H_
+
+#include "MantidKernel/System.h"
+#include "MantidAPI/Algorithm.h"
+
+namespace Mantid {
+namespace Algorithms {
+
+class DLLExport MonitorEfficiencyCorUser : public API::Algorithm {
+public:
+  /// (Empty) Constructor
+  MonitorEfficiencyCorUser(); // : Mantid::API::Algorithm() {}
+  /// Virtual destructor
+  virtual ~MonitorEfficiencyCorUser();
+  /// Algorithm's name
+  virtual const std::string name() const { return "MonitorEfficiencyCorUser"; }
+  /// Summary of algorithms purpose
+  virtual const std::string summary() const {
+    return "This algorithm normalises the counts by the monitor counts with "
+           "additional efficiency correction according to the formula set in "
+           "the instrument parameters file.";
+  }
+
+  /// Algorithm's version
+  virtual int version() const { return 1; }
+  /// Algorithm's category for identification
+  virtual const std::string category() const {
+    return "CorrectionFunctions\\NormalisationCorrections";
+  }
+
+private:
+  /// Initialisation code
+  void init();
+  /// Execution code
+  void exec();
+
+  double calculateFormulaValue(const std::string &, double);
+  std::string getValFromInstrumentDef(const std::string &);
+  void applyMonEfficiency(const size_t numberOfChannels, const MantidVec &yIn,
+                          const MantidVec &eIn, const double coeff,
+                          MantidVec &yOut, MantidVec &eOut);
+
+  /// The user selected (input) workspace
+  API::MatrixWorkspace_const_sptr m_inputWS;
+  /// The output workspace, maybe the same as the input one
+  API::MatrixWorkspace_sptr m_outputWS;
+  /// stores the incident energy of the neutrons
+  double m_Ei;
+  /// stores the total count of neutrons from the monitor
+  int m_monitorCounts;
+};
+
+} // namespace Algorithms
+} // namespace Mantid
+
+#endif /*MONITOREFFICIENCYCORUSER_H_*/

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

+#include "MantidAlgorithms/MonitorEfficiencyCorUser.h"
+#include "MantidAPI/WorkspaceValidators.h"
+#include "MantidKernel/BoundedValidator.h"
+#include "MantidKernel/CompositeValidator.h"
+#include "MantidGeometry/muParser_Silent.h"
+#include "MantidKernel/MultiThreaded.h"
+#include <ctime>
+
+namespace Mantid {
+namespace Algorithms {
+
+using namespace Kernel;
+using namespace API;
+using namespace Geometry;
+
+// Register the algorithm into the AlgorithmFactory
+DECLARE_ALGORITHM(MonitorEfficiencyCorUser)
+
+//----------------------------------------------------------------------------------------------
+/** Constructor
+ */
+MonitorEfficiencyCorUser::MonitorEfficiencyCorUser() {}
+
+//----------------------------------------------------------------------------------------------
+/** Destructor
+ */
+MonitorEfficiencyCorUser::~MonitorEfficiencyCorUser() {}
+
+//----------------------------------------------------------------------------------------------
+/** Initialize the algorithm's properties.
+ */
+
+void MonitorEfficiencyCorUser::init() {
+  declareProperty(
+      new WorkspaceProperty<>("InputWorkspace", "", Direction::Input,
+                              boost::make_shared<InstrumentValidator>()),
+      "The workspace to correct for monitor efficiency");
+  declareProperty(
+      new WorkspaceProperty<>("OutputWorkspace", "", Direction::Output),
+      "The name of the workspace in which to store the result.");
+}
+
+void MonitorEfficiencyCorUser::exec() {
+  m_inputWS = this->getProperty("InputWorkspace");
+
+  m_outputWS = this->getProperty("OutputWorkspace");
+
+  if (m_inputWS->getInstrument()->getName() != "TOFTOF") {
+    std::string message("The input workspace does not come from TOFTOF");
+    g_log.error(message);
+    throw std::invalid_argument(message);
+  }
+
+  // If input and output workspaces are not the same, create a new workspace for
+  // the output
+  if (m_outputWS != this->m_inputWS) {
+    m_outputWS = API::WorkspaceFactory::Instance().create(m_inputWS);
+  }
+  double val;
+  Strings::convert(m_inputWS->run().getProperty("Ei")->value(), val);
+  m_Ei = val;
+  Strings::convert(m_inputWS->run().getProperty("monitor_counts")->value(),
+                   m_monitorCounts);
+
+  // get Efficiency formula from the IDF - Parameters file
+  const std::string effFormula = getValFromInstrumentDef("formula_mon_eff");
+
+  // Calculate Efficiency for E = Ei
+  const double eff0 = m_monitorCounts * calculateFormulaValue(effFormula, m_Ei);
+
+  const size_t numberOfChannels = this->m_inputWS->blocksize();
+  // Calculate the number of spectra in this workspace
+  const int numberOfSpectra =
+      static_cast<int>(this->m_inputWS->getNumberHistograms());
+  API::Progress prog(this, 0.0, 1.0, numberOfSpectra);
+  int64_t numberOfSpectra_i =
+      static_cast<int64_t>(numberOfSpectra); // cast to make openmp happy
+
+  // Loop over the histograms (detector spectra)
+  PARALLEL_FOR2(m_outputWS, m_inputWS)
+  for (int64_t i = 0; i < numberOfSpectra_i; ++i) {
+    PARALLEL_START_INTERUPT_REGION
+    // MantidVec& xOut = m_outputWS->dataX(i);
+    MantidVec &yOut = m_outputWS->dataY(i);
+    MantidVec &eOut = m_outputWS->dataE(i);
+    const MantidVec &yIn = m_inputWS->readY(i);
+    const MantidVec &eIn = m_inputWS->readE(i);
+    m_outputWS->setX(i, m_inputWS->refX(i));
+
+    applyMonEfficiency(numberOfChannels, yIn, eIn, eff0, yOut, eOut);
+
+    prog.report("Detector Efficiency correction...");
+    PARALLEL_END_INTERUPT_REGION
+  } // end for i
+  PARALLEL_CHECK_INTERUPT_REGION
+
+  setProperty("OutputWorkspace", m_outputWS);
+}
+
+/**
+ * Apply the monitor efficiency to a single spectrum
+ * @param numberOfChannels Number of channels in a spectra (nbins - 1)
+ * @param yIn spectrum counts
+ * @param eIn spectrum errors
+ * @param effVec efficiency values (to be divided by the counts)
+ * @param yOut corrected spectrum counts
+ * @param eOut corrected spectrum errors
+ */
+void MonitorEfficiencyCorUser::applyMonEfficiency(
+    const size_t numberOfChannels, const MantidVec &yIn, const MantidVec &eIn,
+    const double effVec, MantidVec &yOut, MantidVec &eOut) {
+
+  for (unsigned int j = 0; j < numberOfChannels; ++j) {
+    yOut[j] = yIn[j] / effVec;
+    eOut[j] = eIn[j] / effVec;
+  }
+}
+
+/**
+ * Calculate the value of a formula
+ * @param formula :: Formula to correct
+ * @param energy :: value of the energy to use in the formula
+ * @return calculated corrected value of the monitor count
+ */
+double
+MonitorEfficiencyCorUser::calculateFormulaValue(const std::string &formula,
+                                                double energy) {
+  try {
+    mu::Parser p;
+    p.DefineVar("e", &energy);
+    p.SetExpr(formula);
+    double eff = p.Eval();
+    g_log.debug() << "Formula: " << formula << " with: " << energy
+                  << "evaluated to: " << eff << std::endl;
+    return eff;
+
+  } catch (mu::Parser::exception_type &e) {
+    throw Kernel::Exception::InstrumentDefinitionError(
+        "Error calculating formula from string. Muparser error message is: " +
+        e.GetMsg());
+  }
+}
+
+/**
+ * Returns the value associated to a parameter name in the IDF
+ * @param parameterName :: parameter name in the IDF
+ * @return the value associated to the parameter name
+ */
+std::string MonitorEfficiencyCorUser::getValFromInstrumentDef(
+    const std::string &parameterName) {
+
+  const ParameterMap &pmap = m_inputWS->constInstrumentParameters();
+  Instrument_const_sptr instrument = m_inputWS->getInstrument();
+  Parameter_sptr par =
+      pmap.getRecursive(instrument->getChild(0).get(), parameterName);
+  if (par) {
+    std::string ret = par->asString();
+    g_log.debug() << "Parsed parameter " << parameterName << ": " << ret
+                  << "\n";
+    return ret;
+  } else {
+    throw Kernel::Exception::InstrumentDefinitionError(
+        "There is no <" + parameterName + "> in the instrument definition!");
+  }
+}
+
+} // namespace Algorithms
+} // namespace Mantid

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

+#ifndef MANTID_ALGORITHMS_MONITOREFFICIENCYCORUSERTEST_H_
+#define MANTID_ALGORITHMS_MONITOREFFICIENCYCORUSERTEST_H_
+
+#include <cxxtest/TestSuite.h>
+
+#include "MantidAlgorithms/MonitorEfficiencyCorUser.h"
+#include "MantidDataObjects/Workspace2D.h"
+#include "MantidTestHelpers/WorkspaceCreationHelper.h"
+#include "MantidTestHelpers/ComponentCreationHelper.h"
+
+using Mantid::Algorithms::MonitorEfficiencyCorUser;
+
+using namespace Mantid;
+using namespace Mantid::API;
+using namespace Kernel;
+
+class MonitorEfficiencyCorUserTest : public CxxTest::TestSuite {
+public:
+  static MonitorEfficiencyCorUserTest *createSuite() {
+    return new MonitorEfficiencyCorUserTest();
+  }
+  static void destroySuite(MonitorEfficiencyCorUserTest *suite) {
+    delete suite;
+  }
+
+  // constructor
+  MonitorEfficiencyCorUserTest()
+      : m_inWSName("input_workspace"), m_outWSName("output_workspace") {
+    m_Ei = 3.27;
+    m_monitor_counts = 1000;
+    createInputWorkSpace();
+  }
+
+  // test init
+  void test_Init() {
+    TS_ASSERT_THROWS_NOTHING(alg.initialize())
+    TS_ASSERT(alg.isInitialized())
+  }
+
+  void testName() { TS_ASSERT_EQUALS(alg.name(), "MonitorEfficiencyCorUser"); }
+
+  void testVersion() { TS_ASSERT_EQUALS(alg.version(), 1); }
+
+  void test_exec() {
+
+    if (!alg.isInitialized())
+      alg.initialize();
+    TS_ASSERT_THROWS_NOTHING(
+        alg.setPropertyValue("InputWorkspace", m_inWSName));
+    TS_ASSERT_THROWS_NOTHING(
+        alg.setPropertyValue("OutputWorkspace", m_outWSName));
+    TS_ASSERT_THROWS_NOTHING(alg.execute(););
+    TS_ASSERT(alg.isExecuted());
+
+    // Retrieve the output workspace from data service.
+    MatrixWorkspace_sptr outWS;
+    TS_ASSERT_THROWS_NOTHING(
+        outWS = AnalysisDataService::Instance().retrieveWS<MatrixWorkspace>(
+            m_outWSName));
+    TS_ASSERT(outWS);
+    if (!outWS)
+      return;
+
+    // Retrieve the input workspace from data service.
+    MatrixWorkspace_sptr inWS;
+    TS_ASSERT_THROWS_NOTHING(
+        inWS = AnalysisDataService::Instance().retrieveWS<MatrixWorkspace>(
+            m_inWSName));
+    TS_ASSERT(inWS);
+    if (!inWS)
+      return;
+
+    // Test equality or proportionality between input and output workspaces
+    const size_t xsize = outWS->blocksize();
+    double proportionality_coeff = m_monitor_counts * sqrt(m_Ei / 25.3);
+    for (size_t i = 0; i < outWS->getNumberHistograms(); ++i) {
+      for (size_t j = 0; j < xsize; ++j) { // Same x-values
+        TS_ASSERT_DELTA(outWS->readX(i)[j], inWS->readX(i)[j], 1e-12);
+        // Output Y-values proportional to input
+        TS_ASSERT_DELTA(proportionality_coeff * outWS->readY(i)[j],
+                        inWS->readY(i)[j], 1e-12);
+
+        // Output Err-values proportional to input
+        TS_ASSERT_DELTA(proportionality_coeff * outWS->readE(i)[j],
+                        inWS->readE(i)[j], 1e-12);
+      }
+    }
+
+    // Remove workspace from the data service.
+    AnalysisDataService::Instance().remove(m_outWSName);
+    AnalysisDataService::Instance().remove(m_inWSName);
+  }
+
+private:
+  double m_Ei;
+  int m_monitor_counts;
+  const std::string m_inWSName, m_outWSName;
+  MonitorEfficiencyCorUser alg;
+
+  void createInputWorkSpace() {
+    int numHist = 1;
+    int numBins = 20;
+
+    DataObjects::Workspace2D_sptr dataws =
+        WorkspaceCreationHelper::create2DWorkspaceWithFullInstrument(
+            numHist, numBins, false, false, true, "TOFTOF");
+
+    static const double arr[] = {
+        -10.,  -9.25, -8.5,  -7.75, -7.,  -6.25, -5.5, -4.75, -4.,  -3.25, -2.5,
+        -1.75, -1.,   -0.25, 0.5,   1.25, 2.,    2.75, 3.5,   4.25, 5.};
+
+    std::vector<double> vec(arr, arr + sizeof(arr) / sizeof(arr[0]));
+
+    for (size_t wi = 0; wi < dataws->getNumberHistograms(); wi++) {
+      dataws->setX(wi, vec);
+    }
+
+    dataws->getAxis(0)->setUnit("TOF");
+    dataws->mutableRun().addProperty("Ei",
+                                     boost::lexical_cast<std::string>(m_Ei));
+    dataws->mutableRun().addProperty(
+        "monitor_counts", boost::lexical_cast<std::string>(m_monitor_counts));
+
+    dataws->instrumentParameters().addString(
+        dataws->getInstrument()->getChild(0).get(), "formula_mon_eff",
+        "sqrt(e/25.3)"); // TOFTOF
+
+    API::AnalysisDataService::Instance().addOrReplace(m_inWSName, dataws);
+  }
+};
+
+#endif /* MANTID_ALGORITHMS_MONITOREFFICIENCYCORUSERTEST_H_ */

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

+.. algorithm::
+
+.. summary::
+
+.. alias::
+
+.. properties::
+
+Description
+-----------
+
+This algorithm normalises the neutron counts by monitor counts with an additional efficiency correction.
+
+To date this algorithm only supports the TOFTOF instrument.
+
+The monitor counts is the total count and it is stored in the SampleLogs of the input workspace.
+
+This count is corrected taking into account the monitor efficiency. The formula used for the correction is stored in the Parameters file and requires the incident energy (Ei), which is stored in the SampleLogs of the input workspace.
+
+The corrected value of the monitor counts is used to normalise the input workspace.
+
+
+Restrictions
+###################################
+
+A formula named "formula\_eff" must be defined in the instrument
+parameters file. It is defined as "monitor counts * sqrt(Ei/25.3)"
+
+The incident energy Ei and the monitor counts are read in the SampleLogs of the input workspace.
+
+Usage
+-----
+
+**Example**
+
+.. testcode:: ExMonitorEfficiencyCorUser
+
+    import numpy as np
+    wsSample = LoadMLZ(Filename='TOFTOFTestdata.nxs')
+    wsNorm = MonitorEfficiencyCorUser(wsSample )
+    # Input and output workspaces have the same structure
+    print 'Number of histograms of the input and output workspaces:'
+    print wsSample.getNumberHistograms(), wsNorm.getNumberHistograms()
+    print 'Number of time channels of the input an output workspaces:'
+    print wsSample.blocksize(), wsNorm.blocksize()
+    # Check calculation of normalisation coefficient between input and output workspaces
+    wsCheck = Divide(wsSample,wsNorm)
+    print "Coefficient of proportionality between Input and Output of MonitorEfficiencyCorUser algorithm: %5.3f" % wsCheck.readY(102)[1]
+    # Read the values of the incident energy and of the monitor counts from the SampleLogs of wsSample
+    monitor_counts = float(mtd['wsSample'].getRun().getLogData('monitor_counts').value)
+    Ei = float(mtd['wsSample'].getRun().getLogData('Ei').value)
+    print "Coefficient from theoretical formula = monitor_counts * sqrt(Ei/25.3): %5.3f" % (monitor_counts*np.sqrt(Ei/25.3))
+ 
+ 
+Output:
+
+.. testoutput:: ExMonitorEfficiencyCorUser
+
+    Number of histograms of the input and output workspaces:
+    1006 1006
+    Number of time channels of the input an output workspaces:
+    1024 1024
+    Coefficient of proportionality between Input and Output of MonitorEfficiencyCorUser algorithm: 41038.432
+    Coefficient from theoretical formula = monitor_counts * sqrt(Ei/25.3): 41038.432 
+
+.. categories::
+
+.. sourcelink::