Mantid
Loading...
Searching...
No Matches
DetectorEfficiencyCorUser.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
4// NScD Oak Ridge National Laboratory, European Spallation Source,
5// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6// SPDX - License - Identifier: GPL - 3.0 +
10#include "MantidAPI/Run.h"
19
20namespace Mantid::Algorithms {
21
22using namespace Kernel;
23using namespace API;
24using namespace Geometry;
25
26// Register the algorithm into the AlgorithmFactory
27DECLARE_ALGORITHM(DetectorEfficiencyCorUser)
28
29
30const std::string DetectorEfficiencyCorUser::name() const { return "DetectorEfficiencyCorUser"; }
31
33int DetectorEfficiencyCorUser::version() const { return 1; }
34
36const std::string DetectorEfficiencyCorUser::category() const {
37 return "CorrectionFunctions\\EfficiencyCorrections;Inelastic\\Corrections";
38}
39
43 auto val = std::make_shared<CompositeValidator>();
44 // val->add<WorkspaceUnitValidator>("Energy");
45 val->add<WorkspaceUnitValidator>("DeltaE");
46 val->add<HistogramValidator>();
47 val->add<InstrumentValidator>();
48 declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input, val),
49 "The workspace to correct for detector efficiency");
50 declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
51 "The name of the workspace in which to store the result.");
52 auto checkEi = std::make_shared<BoundedValidator<double>>();
53 checkEi->setLower(0.0);
54 declareProperty("IncidentEnergy", EMPTY_DBL(), checkEi, "The energy of neutrons leaving the source.");
55}
56
60 // get input properties (WSs, Ei)
62
63 const size_t numberOfChannels = this->m_inputWS->blocksize();
64 // Calculate the number of spectra in this workspace
65 const auto numberOfSpectra = static_cast<int>(this->m_inputWS->size() / numberOfChannels);
66 API::Progress prog(this, 0.0, 1.0, numberOfSpectra);
67 auto numberOfSpectra_i = static_cast<int64_t>(numberOfSpectra); // cast to make openmp happy
68
69 // Loop over the histograms (detector spectra)
71 for (int64_t i = 0; i < numberOfSpectra_i; ++i) {
73 const auto effFormula = retrieveFormula(i);
74 // Calculate Efficiency for E = Ei
75 double e;
76 auto parser = generateParser(effFormula, &e);
77 e = m_Ei;
78 const double eff0 = evaluate(parser);
79 correctHistogram(i, eff0, e, parser);
80
81 prog.report("Detector Efficiency correction...");
82
84 } // end for i
86
87 this->setProperty("OutputWorkspace", this->m_outputWS);
88}
89
98void DetectorEfficiencyCorUser::correctHistogram(const size_t index, const double eff0, double &e,
99 const mu::Parser &parser) {
100 const auto &xIn = m_inputWS->points(index);
101 const auto &yIn = m_inputWS->y(index);
102 const auto &eIn = m_inputWS->e(index);
103 auto &yOut = m_outputWS->mutableY(index);
104 auto &eOut = m_outputWS->mutableE(index);
105 for (size_t i = 0; i < xIn.size(); ++i) {
106 e = m_Ei - xIn[i];
107 const double eff = evaluate(parser);
108 const double corr = eff / eff0;
109 yOut[i] = yIn[i] / corr;
110 eOut[i] = eIn[i] / corr;
111 }
112}
113
120double DetectorEfficiencyCorUser::evaluate(const mu::Parser &parser) const {
121 try {
122 return parser.Eval();
123 } catch (mu::Parser::exception_type &e) {
125 "Error calculating formula from string. Muparser error message is: " + e.GetMsg());
126 }
127}
128
129mu::Parser DetectorEfficiencyCorUser::generateParser(const std::string &formula, double *e) const {
130 mu::Parser p;
131 p.DefineVar("e", e);
132 p.SetExpr(formula);
133 return p;
134}
135
141std::string DetectorEfficiencyCorUser::retrieveFormula(const size_t workspaceIndex) {
142 const std::string formulaParamName("formula_eff");
143 const auto &paramMap = m_inputWS->constInstrumentParameters();
144 auto det = m_inputWS->getDetector(workspaceIndex);
145 auto param = paramMap.getRecursive(det.get(), formulaParamName, "string");
146 if (!param) {
147 throw Kernel::Exception::InstrumentDefinitionError("No <" + formulaParamName + "> parameter found for component '" +
148 det->getFullName() + "' in the instrument definition.");
149 }
150 const auto ret = param->asString();
151 g_log.debug() << "Found formula for workspace index " << workspaceIndex << ": " << ret << "\n";
152 return ret;
153}
154
161 // Get the workspaces
162 m_inputWS = this->getProperty("InputWorkspace");
163
164 m_outputWS = this->getProperty("OutputWorkspace");
165
166 // If input and output workspaces are not the same, create a new workspace for
167 // the output
168 if (m_outputWS != this->m_inputWS) {
169 m_outputWS.reset(Mantid::DataObjects::create<DataObjects::Workspace2D>(*m_inputWS).release());
170 }
171
172 // these first three properties are fully checked by validators
173 m_Ei = this->getProperty("IncidentEnergy");
174 // If we're not given an Ei, see if one has been set.
175 if (m_Ei == EMPTY_DBL()) {
176 Mantid::Kernel::Property *prop = m_inputWS->run().getProperty("Ei");
177 double val;
178 if (!prop || !Strings::convert(prop->value(), val)) {
179 throw std::invalid_argument("No Ei value has been set or stored within the run information.");
180 }
181 m_Ei = val;
182 g_log.debug() << "Using stored Ei value " << m_Ei << "\n";
183 } else {
184 g_log.debug() << "Using user input Ei value: " << m_Ei << "\n";
185 }
186}
187} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
std::map< DeltaEMode::Type, std::string > index
Definition: DeltaEMode.cpp:19
#define PARALLEL_START_INTERRUPT_REGION
Begins a block to skip processing is the algorithm has been interupted Note the end of the block if n...
Definition: MultiThreaded.h:94
#define PARALLEL_END_INTERRUPT_REGION
Ends a block to skip processing is the algorithm has been interupted Note the start of the block if n...
#define PARALLEL_FOR_IF(condition)
Empty definitions - to enable set your complier to enable openMP.
#define PARALLEL_CHECK_INTERRUPT_REGION
Adds a check after a Parallel region to see if it was interupted.
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
Definition: Algorithm.cpp:1913
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Definition: Algorithm.cpp:2076
Kernel::Logger & g_log
Definition: Algorithm.h:451
A validator which checks that a workspace contains histogram data (the default) or point data as requ...
A validator which checks that a workspace has a valid instrument.
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
A property class for workspaces.
A validator which checks that the unit of the workspace referred to by a WorkspaceProperty is the exp...
double m_Ei
stores the user selected value for incidient energy of the neutrons
mu::Parser generateParser(const std::string &formula, double *e) const
const std::string category() const override
Algorithm's category for identification.
API::MatrixWorkspace_sptr m_outputWS
The output workspace, maybe the same as the input one.
void retrieveProperties()
Loads and checks the values passed to the algorithm.
double evaluate(const mu::Parser &parser) const
Calculate the value of a formula parsed by muParser.
API::MatrixWorkspace_const_sptr m_inputWS
The user selected (input) workspace.
void init() override
Initialize the algorithm's properties.
int version() const override
Algorithm's version for identification.
void correctHistogram(const size_t index, const double eff0, double &e, const mu::Parser &parser)
Apply efficiency corrections to a histogram in the output workspace.
std::string retrieveFormula(const size_t workspaceIndex)
Returns the efficiency correction formula associated to a detector.
Exception for errors associated with the instrument definition.
Definition: Exception.h:220
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void debug(const std::string &msg)
Logs at debug level.
Definition: Logger.cpp:114
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
Definition: ProgressBase.h:51
Base class for properties.
Definition: Property.h:94
virtual std::string value() const =0
Returns the value of the property as a string.
int convert(const std::string &A, T &out)
Convert a string into a number.
Definition: Strings.cpp:665
std::enable_if< std::is_pointer< Arg >::value, bool >::type threadSafe(Arg workspace)
Thread-safety check Checks the workspace to ensure it is suitable for multithreaded access.
Definition: MultiThreaded.h:22
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
Definition: EmptyValues.h:43
STL namespace.
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54