DetectorEfficiencyCorUser.cpp

Go to the documentation of this file.

// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
//   NScD Oak Ridge National Laboratory, European Spallation Source,
//   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
#include "MantidAlgorithms/DetectorEfficiencyCorUser.h"
#include "MantidAPI/HistogramValidator.h"
#include "MantidAPI/InstrumentValidator.h"
#include "MantidAPI/Run.h"
#include "MantidAPI/SpectrumInfo.h"
#include "MantidAPI/WorkspaceUnitValidator.h"
#include "MantidDataObjects/Workspace2D.h"
#include "MantidDataObjects/WorkspaceCreation.h"
#include "MantidGeometry/Instrument.h"
#include "MantidGeometry/Instrument/ComponentInfo.h"
#include "MantidGeometry/muParser_Silent.h"
#include "MantidKernel/BoundedValidator.h"
#include "MantidKernel/CompositeValidator.h"
#include "MantidKernel/Strings.h"
 
namespace Mantid::Algorithms {
 
using namespace Kernel;
using namespace API;
using namespace Geometry;
 
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(DetectorEfficiencyCorUser)
 
 
const std::string DetectorEfficiencyCorUser::name() const { return "DetectorEfficiencyCorUser"; }
 
int DetectorEfficiencyCorUser::version() const { return 1; }
 
const std::string DetectorEfficiencyCorUser::category() const {
  return "CorrectionFunctions\\EfficiencyCorrections;Inelastic\\Corrections";
}
 
void DetectorEfficiencyCorUser::init() {
  auto val = std::make_shared<CompositeValidator>();
  // val->add<WorkspaceUnitValidator>("Energy");
  val->add<WorkspaceUnitValidator>("DeltaE");
  val->add<HistogramValidator>();
  val->add<InstrumentValidator>();
  declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input, val),
                  "The workspace to correct for detector efficiency");
  declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
                  "The name of the workspace in which to store the result.");
  auto checkEi = std::make_shared<BoundedValidator<double>>();
  checkEi->setLower(0.0);
  declareProperty("IncidentEnergy", EMPTY_DBL(), checkEi, "The energy of neutrons leaving the source.");
}
 
void DetectorEfficiencyCorUser::exec() {
  // get input properties (WSs, Ei)
  retrieveProperties();
 
  const size_t numberOfChannels = this->m_inputWS->blocksize();
  // Calculate the number of spectra in this workspace
  const auto numberOfSpectra = static_cast<int>(this->m_inputWS->size() / numberOfChannels);
  API::Progress prog(this, 0.0, 1.0, numberOfSpectra);
  auto numberOfSpectra_i = static_cast<int64_t>(numberOfSpectra); // cast to make openmp happy
 
  // Loop over the histograms (detector spectra)
  PARALLEL_FOR_IF(Kernel::threadSafe(*m_outputWS, *m_inputWS))
  for (int64_t i = 0; i < numberOfSpectra_i; ++i) {
    PARALLEL_START_INTERRUPT_REGION
    const auto effFormula = retrieveFormula(i);
    // Calculate Efficiency for E = Ei
    double e;
    auto parser = generateParser(effFormula, &e);
    e = m_Ei;
    const double eff0 = evaluate(parser);
    correctHistogram(i, eff0, e, parser);
 
    prog.report("Detector Efficiency correction...");
 
    PARALLEL_END_INTERRUPT_REGION
  } // end for i
  PARALLEL_CHECK_INTERRUPT_REGION
 
  this->setProperty("OutputWorkspace", this->m_outputWS);
}
 
void DetectorEfficiencyCorUser::correctHistogram(const size_t index, const double eff0, double &e,
                                                 const mu::Parser &parser) {
  const auto &xIn = m_inputWS->points(index);
  const auto &yIn = m_inputWS->y(index);
  const auto &eIn = m_inputWS->e(index);
  auto &yOut = m_outputWS->mutableY(index);
  auto &eOut = m_outputWS->mutableE(index);
  for (size_t i = 0; i < xIn.size(); ++i) {
    e = m_Ei - xIn[i];
    const double eff = evaluate(parser);
    const double corr = eff / eff0;
    yOut[i] = yIn[i] / corr;
    eOut[i] = eIn[i] / corr;
  }
}
 
double DetectorEfficiencyCorUser::evaluate(const mu::Parser &parser) const {
  try {
    return parser.Eval();
  } catch (mu::Parser::exception_type &e) {
    throw Kernel::Exception::InstrumentDefinitionError(
        "Error calculating formula from string. Muparser error message is: " + e.GetMsg());
  }
}
 
mu::Parser DetectorEfficiencyCorUser::generateParser(const std::string &formula, double *e) const {
  mu::Parser p;
  p.DefineVar("e", e);
  p.SetExpr(formula);
  return p;
}
 
std::string DetectorEfficiencyCorUser::retrieveFormula(const size_t workspaceIndex) {
  const std::string formulaParamName("formula_eff");
  const auto &spectrumInfo = m_inputWS->spectrumInfo();
  if (!spectrumInfo.hasDetectors(workspaceIndex)) {
    throw Kernel::Exception::NotFoundError("MatrixWorkspace::getDetector(): No detectors for this workspace index.",
                                           std::to_string(workspaceIndex));
  }
  const auto &componentInfo = m_inputWS->componentInfo();
  const auto detIndex = spectrumInfo.firstDetectorIndex(workspaceIndex);
  // A spectrum of several detectors has no parameters of its own
  const auto param = spectrumInfo.hasUniqueDetector(workspaceIndex)
                         ? componentInfo.parameter(detIndex, formulaParamName, "string")
                         : nullptr;
  if (!param) {
    throw Kernel::Exception::InstrumentDefinitionError("No <" + formulaParamName + "> parameter found for component '" +
                                                       componentInfo.fullName(detIndex) +
                                                       "' in the instrument definition.");
  }
  const auto ret = param->asString();
  g_log.debug() << "Found formula for workspace index " << workspaceIndex << ": " << ret << "\n";
  return ret;
}
 
void DetectorEfficiencyCorUser::retrieveProperties() {
  // Get the workspaces
  m_inputWS = this->getProperty("InputWorkspace");
 
  m_outputWS = this->getProperty("OutputWorkspace");
 
  // If input and output workspaces are not the same, create a new workspace for
  // the output
  if (m_outputWS != this->m_inputWS) {
    m_outputWS.reset(Mantid::DataObjects::create<DataObjects::Workspace2D>(*m_inputWS).release());
  }
 
  // these first three properties are fully checked by validators
  m_Ei = this->getProperty("IncidentEnergy");
  // If we're not given an Ei, see if one has been set.
  if (m_Ei == EMPTY_DBL()) {
    Mantid::Kernel::Property const *prop = m_inputWS->run().getProperty("Ei");
    double val;
    if (!prop || !Strings::convert(prop->value(), val)) {
      throw std::invalid_argument("No Ei value has been set or stored within the run information.");
    }
    m_Ei = val;
    g_log.debug() << "Using stored Ei value " << m_Ei << "\n";
  } else {
    g_log.debug() << "Using user input Ei value: " << m_Ei << "\n";
  }
}
} // namespace Mantid::Algorithms