Mantid
Loading...
Searching...
No Matches
ModeratorTzero.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 +
8#include "MantidAPI/Axis.h"
9#include "MantidAPI/Run.h"
18#include "MantidHistogramData/Histogram.h"
21
22namespace Mantid::Algorithms {
23
24// Register the algorithm into the AlgorithmFactory
25DECLARE_ALGORITHM(ModeratorTzero)
26
27using namespace Mantid::Kernel;
28using namespace Mantid::API;
29using namespace Mantid::Geometry;
30using namespace Mantid::DataObjects;
31using namespace Mantid::HistogramData;
32
34 : Mantid::API::Algorithm(),
35 m_convfactor(0.5e+12 * Mantid::PhysicalConstants::NeutronMass / Mantid::PhysicalConstants::meV), m_niter(1),
36 m_tolTOF(0.), m_t1min(200.0) {}
37
39void ModeratorTzero::setFormula(const std::string &formula) { m_formula = formula; }
40
42
43 auto wsValidator = std::make_shared<WorkspaceUnitValidator>("TOF");
45 std::make_unique<WorkspaceProperty<MatrixWorkspace>>("InputWorkspace", "", Direction::Input, wsValidator),
46 "The name of the input workspace, containing events and/or "
47 "histogram data, in units of time-of-flight");
48 // declare the output workspace
49 declareProperty(std::make_unique<WorkspaceProperty<MatrixWorkspace>>("OutputWorkspace", "", Direction::Output),
50 "The name of the output workspace");
51
52 // declare the instrument scattering mode
53 std::vector<std::string> EModeOptions{"Indirect", "Direct", "Elastic"};
54 this->declareProperty("EMode", "Indirect", std::make_shared<StringListValidator>(EModeOptions),
55 "The energy mode (default: Indirect)");
56
58 "Tolerance in the calculation of the emission time, in "
59 "microseconds (default:1)");
61 "Number of iterations (default:1)");
62} // end of void ModeratorTzero::init()
63
65 m_tolTOF = getProperty("tolTOF"); // Tolerance in the calculation of the
66 // emission time, in microseconds
67 m_niter = getProperty("Niter"); // number of iterations
68 const MatrixWorkspace_sptr inputWS = getProperty("InputWorkspace");
69 const std::string emode = getProperty("EMode");
70
71 // Check if Ei stored in workspace
72 if (emode == "Direct") {
73 if (!inputWS->run().hasProperty("Ei")) {
74 g_log.error("No incident energy value (Ei) has been set or stored.");
75 return;
76 }
77 }
78
79 // extract formula from instrument parameters
80 auto t0_formula = inputWS->getInstrument()->getStringParameter("t0_formula");
81 if (t0_formula.empty()) {
82 g_log.error("Unable to retrieve t0_formula among instrument parameters.");
83 return;
84 }
85 m_formula = t0_formula[0];
86
87 // Run execEvent if eventWorkSpace
88 EventWorkspace_const_sptr eventWS = std::dynamic_pointer_cast<const EventWorkspace>(inputWS);
89 if (eventWS != nullptr) {
90 execEvent(emode);
91 return;
92 }
93
94 // Run exec for matrix workspace
95 MatrixWorkspace_sptr outputWS = getProperty("OutputWorkspace");
96 // Check whether input == output to see whether a new workspace is required.
97 if (outputWS != inputWS) {
98 // Create new workspace for output from old
99 outputWS = create<MatrixWorkspace>(*inputWS);
100 }
101
102 // calculate tof shift once for all neutrons if emode==Direct
103 double t0_direct(-1);
104 if (emode == "Direct") {
105 auto Ei = inputWS->run().getPropertyValueAsType<double>("Ei");
106 mu::Parser parser;
107 parser.DefineVar("incidentEnergy", &Ei); // associate E1 to this parser
108 parser.SetExpr(m_formula);
109 t0_direct = parser.Eval();
110 }
111
112 const auto &spectrumInfo = inputWS->spectrumInfo();
113 const double Lss = spectrumInfo.l1();
114
115 const auto numHists = static_cast<size_t>(inputWS->getNumberHistograms());
116 Progress prog(this, 0.0, 1.0, numHists); // report progress of algorithm
117 PARALLEL_FOR_IF(Kernel::threadSafe(*inputWS, *outputWS))
118 // iterate over the spectra
119 for (int i = 0; i < static_cast<int>(numHists); ++i) {
121 outputWS->setHistogram(i, inputWS->histogram(i));
122
123 // One parser for each parallel processor needed (except Edirect mode)
124 double E1;
125 mu::Parser parser;
126 parser.DefineVar("incidentEnergy", &E1); // associate E1 to this parser
127 parser.SetExpr(m_formula);
128
129 double L1(Lss); // distance from source to sample
130 double L2(-1); // distance from sample to detector
131 if (spectrumInfo.hasDetectors(i)) {
132 if (spectrumInfo.isMonitor(i)) {
133 // redefine the sample as the monitor
134 L1 = Lss + spectrumInfo.l2(i); // L2 in SpectrumInfo defined negative
135 L2 = 0;
136 } else {
137 L2 = spectrumInfo.l2(i);
138 }
139 } else {
140 g_log.error() << "Unable to calculate distances to/from detector" << i << '\n';
141 }
142
143 if (L2 >= 0) {
144 // fast neutrons are shifted by min_t0_next, irrespective of tof
145 double v1_max = L1 / m_t1min;
146 E1 = m_convfactor * v1_max * v1_max;
147 double min_t0_next = parser.Eval();
148
149 if (emode == "Indirect") {
150 double t2(-1.0); // time from sample to detector. (-1) signals error
151 if (spectrumInfo.isMonitor(i)) {
152 t2 = 0.0;
153 } else {
154 static const double convFact = 1.0e-6 * sqrt(2 * PhysicalConstants::meV / PhysicalConstants::NeutronMass);
155 std::vector<double> wsProp = spectrumInfo.detector(i).getNumberParameter("Efixed");
156 if (!wsProp.empty()) {
157 double E2 = wsProp.at(0); //[E2]=meV
158 double v2 = convFact * sqrt(E2); //[v2]=meter/microsec
159 t2 = L2 / v2;
160 } else {
161 // t2 is kept to -1 if no Efixed is found
162 g_log.debug() << "Efixed not found for detector " << i << '\n';
163 }
164 }
165 // shift the time of flights by the emission time from the moderator
166 if (t2 >= 0) // t2 < 0 when no detector info is available
167 {
168 auto &outbins = outputWS->mutableX(i);
169 for (auto &tof : outbins) {
170 if (tof < m_t1min + t2)
171 tof -= min_t0_next;
172 else
173 tof -= CalculateT0indirect(tof, L1, t2, E1, parser);
174 }
175 }
176 } // end of if(emode=="Indirect")
177 else if (emode == "Elastic") {
178 auto &outbins = outputWS->mutableX(i);
179 for (auto &tof : outbins) {
180 if (tof < m_t1min * (L1 + L2) / L1)
181 tof -= min_t0_next;
182 else
183 tof -= CalculateT0elastic(tof, L1 + L2, E1, parser);
184 }
185 } // end of else if(emode=="Elastic")
186 else if (emode == "Direct") {
187 outputWS->mutableX(i) += -t0_direct;
188 } // end of else if(emode="Direct")
189 } // end of if(L2 >= 0)
190
191 prog.report();
193 } // end of for (int i = 0; i < static_cast<int>(numHists); ++i)
195
196 // Copy units
197 if (inputWS->getAxis(0)->unit().get()) {
198 outputWS->getAxis(0)->unit() = inputWS->getAxis(0)->unit();
199 }
200 try {
201 if (inputWS->getAxis(1)->unit().get()) {
202 outputWS->getAxis(1)->unit() = inputWS->getAxis(1)->unit();
203 }
204 } catch (Exception::IndexError &) {
205 // OK, so this isn't a Workspace2D
206 }
207
208 // Assign it to the output workspace property
209} // end of void ModeratorTzero::exec
210
211void ModeratorTzero::execEvent(const std::string &emode) {
212 g_log.information("Processing event workspace");
213
214 const MatrixWorkspace_const_sptr matrixInputWS = getProperty("InputWorkspace");
215
216 // generate the output workspace pointer
217 API::MatrixWorkspace_sptr matrixOutputWS = getProperty("OutputWorkspace");
218 if (matrixOutputWS != matrixInputWS) {
219 matrixOutputWS = matrixInputWS->clone();
220 setProperty("OutputWorkspace", matrixOutputWS);
221 }
222 auto outputWS = std::dynamic_pointer_cast<EventWorkspace>(matrixOutputWS);
223
224 // calculate tof shift once for all neutrons if emode==Direct
225 double t0_direct(-1);
226 if (emode == "Direct") {
227 auto Ei = outputWS->run().getPropertyValueAsType<double>("Ei");
228 mu::Parser parser;
229 parser.DefineVar("incidentEnergy", &Ei); // associate E1 to this parser
230 parser.SetExpr(m_formula);
231 t0_direct = parser.Eval();
232 }
233
234 const auto &spectrumInfo = outputWS->spectrumInfo();
235 const double Lss = spectrumInfo.l1();
236
237 // Loop over the spectra
238 const auto numHists = static_cast<size_t>(outputWS->getNumberHistograms());
239 Progress prog(this, 0.0, 1.0, numHists); // report progress of algorithm
241 for (int i = 0; i < static_cast<int>(numHists); ++i) {
243 auto wsIndex = static_cast<size_t>(i);
244 EventList &evlist = outputWS->getSpectrum(wsIndex);
245 if (evlist.getNumberEvents() > 0) // don't bother with empty lists
246 {
247 double L1(Lss); // distance from source to sample
248 double L2(-1); // distance from sample to detector
249
250 if (spectrumInfo.hasDetectors(i)) {
251 if (spectrumInfo.isMonitor(i)) {
252 // redefine the sample as the monitor
253 L1 = Lss + spectrumInfo.l2(i); // L2 in SpectrumInfo defined negative
254 L2 = 0;
255 } else {
256 L2 = spectrumInfo.l2(i);
257 }
258 } else {
259 g_log.error() << "Unable to calculate distances to/from detector" << i << '\n';
260 }
261
262 if (L2 >= 0) {
263 // One parser for each parallel processor needed (except Edirect mode)
264 double E1;
265 mu::Parser parser;
266 parser.DefineVar("incidentEnergy", &E1); // associate E1 to this parser
267 parser.SetExpr(m_formula);
268
269 // fast neutrons are shifted by min_t0_next, irrespective of tof
270 double v1_max = L1 / m_t1min;
271 E1 = m_convfactor * v1_max * v1_max;
272 double min_t0_next = parser.Eval();
273
274 if (emode == "Indirect") {
275 double t2(-1.0); // time from sample to detector. (-1) signals error
276 if (spectrumInfo.isMonitor(i)) {
277 t2 = 0.0;
278 } else {
279 static const double convFact = 1.0e-6 * sqrt(2 * PhysicalConstants::meV / PhysicalConstants::NeutronMass);
280 std::vector<double> wsProp = spectrumInfo.detector(i).getNumberParameter("Efixed");
281 if (!wsProp.empty()) {
282 double E2 = wsProp.at(0); //[E2]=meV
283 double v2 = convFact * sqrt(E2); //[v2]=meter/microsec
284 t2 = L2 / v2;
285 } else {
286 // t2 is kept to -1 if no Efixed is found
287 g_log.debug() << "Efixed not found for detector " << i << '\n';
288 }
289 }
290 if (t2 >= 0) // t2 < 0 when no detector info is available
291 {
292 // fix the histogram bins
293 auto &x = evlist.mutableX();
294 for (double &tof : x) {
295 if (tof < m_t1min + t2)
296 tof -= min_t0_next;
297 else
298 tof -= CalculateT0indirect(tof, L1, t2, E1, parser);
299 }
300
301 MantidVec tofs = evlist.getTofs();
302 for (double &tof : tofs) {
303 if (tof < m_t1min + t2)
304 tof -= min_t0_next;
305 else
306 tof -= CalculateT0indirect(tof, L1, t2, E1, parser);
307 }
308 evlist.setTofs(tofs);
310 } // end of if( t2>= 0)
311 } // end of if(emode=="Indirect")
312 else if (emode == "Elastic") {
313 // Apply t0 correction to histogram bins
314 auto &x = evlist.mutableX();
315 for (double &tof : x) {
316 if (tof < m_t1min * (L1 + L2) / L1)
317 tof -= min_t0_next;
318 else
319 tof -= CalculateT0elastic(tof, L1 + L2, E1, parser);
320 }
321
322 MantidVec tofs = evlist.getTofs();
323 for (double &tof : tofs) {
324 // add a [-0.1,0.1] microsecond noise to avoid artifacts
325 // resulting from original tof data
326 if (tof < m_t1min * (L1 + L2) / L1)
327 tof -= min_t0_next;
328 else
329 tof -= CalculateT0elastic(tof, L1 + L2, E1, parser);
330 }
331 evlist.setTofs(tofs);
333 } // end of else if(emode=="Elastic")
334 else if (emode == "Direct") {
335 // fix the histogram bins
336 evlist.mutableX() -= t0_direct;
337
338 MantidVec tofs = evlist.getTofs();
339
340 std::transform(tofs.begin(), tofs.end(), tofs.begin(), [&t0_direct](double tof) { return tof - t0_direct; });
341 evlist.setTofs(tofs);
343 } // end of else if(emode=="Direct")
344 } // end of if(L2 >= 0)
345 } // end of if (evlist.getNumberEvents() > 0)
346 prog.report();
348 } // end of for (int i = 0; i < static_cast<int>(numHists); ++i)
350 outputWS->clearMRU(); // Clears the Most Recent Used lists */
351} // end of void ModeratorTzero::execEvent()
352
355double ModeratorTzero::CalculateT0indirect(const double &tof, const double &L1, const double &t2, double &E1,
356 mu::Parser &parser) {
357
358 double t0_curr, t0_next;
359 t0_curr = m_tolTOF; // current iteration emission time
360 t0_next = 0.0; // next iteration emission time, initialized to zero
361 size_t iiter(0); // current iteration number
362 // iterate until convergence in t0 reached
363 while (std::fabs(t0_curr - t0_next) >= m_tolTOF && iiter < m_niter) {
364 t0_curr = t0_next;
365 double t1 = tof - t0_curr - t2;
366 double v1 = L1 / t1;
367 E1 = m_convfactor * v1 * v1; // Energy in meV if v1 in meter/microsecond
368 t0_next = parser.Eval();
369 iiter++;
370 }
371 return t0_next;
372}
373
376double ModeratorTzero::CalculateT0elastic(const double &tof, const double &L12, double &E1, mu::Parser &parser) {
377 double t0_curr, t0_next;
378 t0_curr = m_tolTOF; // current iteration emission time
379 t0_next = 0.0; // next iteration emission time, initialized to zero
380 size_t iiter(0); // current iteration number
381 // iterate until convergence in t0 reached
382 while (std::fabs(t0_curr - t0_next) >= m_tolTOF && iiter < m_niter) {
383 t0_curr = t0_next;
384 double v1 = L12 / (tof - t0_curr); // v1 = v2 = v since emode is elastic
385 E1 = m_convfactor * v1 * v1; // Energy in meV if v1 in meter/microsecond
386 t0_next = parser.Eval();
387 iiter++;
388 }
389 return t0_next;
390}
392
393} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
#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.
Base class from which all concrete algorithm classes should be derived.
Definition: Algorithm.h:85
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
HistogramData::HistogramX & mutableX() &
Definition: ISpectrum.h:176
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
A property class for workspaces.
const double m_t1min
tof limit for fast neutrons
void execEvent(const std::string &emode)
Execution code for event workspace.
double CalculateT0elastic(const double &tof, const double &L12, double &E1, mu::Parser &parser)
Calculate emission time from the moderator for a given detector (L1, t2) and TOF when Emode==Elastic.
std::string m_formula
string containing the heuristic regression for the moderator emission time versus neutron energy
void init() override
Virtual method - must be overridden by concrete algorithm.
ModeratorTzero()
Default constructor.
double m_tolTOF
tolerance for calculating E1, in micro-seconds
void exec() override
Execution code for histogram workspace.
size_t m_niter
Maximum number of iterations when calculating the emission time from the moderator.
void setFormula(const std::string &formula)
set attribute m_formula
double CalculateT0indirect(const double &tof, const double &L1, const double &t2, double &E1, mu::Parser &parser)
Calculate emission time from the moderator for a given detector (L1, t2) and TOF when Emode==Inelasti...
A class for holding :
Definition: EventList.h:56
void setTofs(const MantidVec &tofs) override
Set a list of TOFs to the current event list.
Definition: EventList.cpp:3226
void setSortOrder(const EventSortType order) const
Manually set the event list sort order value.
Definition: EventList.cpp:943
std::size_t getNumberEvents() const override
Return the number of events in the list.
Definition: EventList.cpp:1143
void getTofs(std::vector< double > &tofs) const override
Fill a vector with the list of TOFs.
Definition: EventList.cpp:2734
Exception for index errors.
Definition: Exception.h:284
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 error(const std::string &msg)
Logs at error level.
Definition: Logger.cpp:77
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
Definition: ProgressBase.h:51
The concrete, templated class for properties.
std::complex< double > MANTID_API_DLL E1(std::complex< double > z)
Integral for Gamma.
std::shared_ptr< const MatrixWorkspace > MatrixWorkspace_const_sptr
shared pointer to the matrix workspace base class (const version)
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
std::shared_ptr< const EventWorkspace > EventWorkspace_const_sptr
shared pointer to a const Workspace2D
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
A namespace containing physical constants that are required by algorithms and unit routines.
Definition: Atom.h:14
static constexpr double NeutronMass
Mass of the neutron in kg.
static constexpr double meV
1 meV in Joules.
Helper class which provides the Collimation Length for SANS instruments.
std::vector< double > MantidVec
typedef for the data storage used in Mantid matrix workspaces
Definition: cow_ptr.h:172
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54