Mantid
Loading...
Searching...
No Matches
He3TubeEfficiency.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"
23
24#include <cmath>
25#include <stdexcept>
26
27// this should be a big number but not so big that there are rounding errors
28const double DIST_TO_UNIVERSE_EDGE = 1e3;
29
30// Scalar constant for exponential in units of K / (m Angstroms atm)
31const double EXP_SCALAR_CONST = 2175.486863864;
32
33// Tolerance for diameter/thickness comparison
34const double TOL = 1.0e-8;
35
36namespace Mantid::Algorithms {
37
38using namespace DataObjects;
39using namespace HistogramData;
40// Register the class into the algorithm factory
41DECLARE_ALGORITHM(He3TubeEfficiency)
42
43
45 : Algorithm(), m_inputWS(), m_outputWS(), m_paraMap(nullptr), m_shapeCache(), m_samplePos(), m_spectraSkipped(),
46 m_progress(nullptr) {
47 m_shapeCache.clear();
48}
49
54 using namespace Mantid::Kernel;
55
56 auto wsValidator = std::make_shared<CompositeValidator>();
57 wsValidator->add<API::WorkspaceUnitValidator>("Wavelength");
58 wsValidator->add<API::HistogramValidator>();
59 wsValidator->add<API::InstrumentValidator>();
61 "InputWorkspace", "", Kernel::Direction::Input, wsValidator),
62 "Name of the input workspace");
63 this->declareProperty(
64 std::make_unique<API::WorkspaceProperty<API::MatrixWorkspace>>("OutputWorkspace", "", Kernel::Direction::Output),
65 "Name of the output workspace, can be the same as the input");
66 auto mustBePositive = std::make_shared<Kernel::BoundedValidator<double>>();
67 mustBePositive->setLower(0.0);
68 this->declareProperty(std::make_unique<Kernel::PropertyWithValue<double>>("ScaleFactor", 1.0, mustBePositive),
69 "Constant factor with which to scale the calculated"
70 "detector efficiency. Same factor applies to all efficiencies.");
71
72 auto mustBePosArr = std::make_shared<Kernel::ArrayBoundedValidator<double>>();
73 mustBePosArr->setLower(0.0);
74 this->declareProperty(std::make_unique<Kernel::ArrayProperty<double>>("TubePressure", mustBePosArr),
75 "Provide overriding the default tube pressure. The pressure must "
76 "be specified in atm.");
77 this->declareProperty(std::make_unique<Kernel::ArrayProperty<double>>("TubeThickness", mustBePosArr),
78 "Provide overriding the default tube thickness. The thickness must "
79 "be specified in metres.");
80 this->declareProperty(std::make_unique<Kernel::ArrayProperty<double>>("TubeTemperature", mustBePosArr),
81 "Provide overriding the default tube temperature. The temperature must "
82 "be specified in Kelvin.");
83}
84
89 // Get the workspaces
90 m_inputWS = this->getProperty("InputWorkspace");
91 m_outputWS = this->getProperty("OutputWorkspace");
92
93 if (m_outputWS != m_inputWS) {
94 m_outputWS = create<API::MatrixWorkspace>(*m_inputWS);
95 }
96
97 // Get the detector parameters
98 m_paraMap = &(m_inputWS->constInstrumentParameters());
99
100 // Store some information about the instrument setup that will not change
101 m_samplePos = m_inputWS->getInstrument()->getSample()->getPos();
102
103 // Check if it is an event workspace
105 std::dynamic_pointer_cast<const DataObjects::EventWorkspace>(m_inputWS);
106 if (eventW != nullptr) {
107 this->execEvent();
108 return;
109 }
110
111 std::size_t numHists = m_inputWS->getNumberHistograms();
112 m_progress = std::make_unique<API::Progress>(this, 0.0, 1.0, numHists);
113 const auto &spectrumInfo = m_inputWS->spectrumInfo();
114
116 for (int i = 0; i < static_cast<int>(numHists); ++i) {
118
119 m_outputWS->setSharedX(i, m_inputWS->sharedX(i));
120 try {
121 this->correctForEfficiency(i, spectrumInfo);
122 } catch (std::out_of_range &) {
123 // if we don't have all the data there will be spectra we can't correct,
124 // avoid leaving the workspace part corrected
125 m_outputWS->mutableY(i) = 0;
126 PARALLEL_CRITICAL(deteff_invalid) { m_spectraSkipped.emplace_back(m_inputWS->getAxis(1)->spectraNo(i)); }
127 }
128
129 // make regular progress reports
130 m_progress->report();
131 // check for canceling the algorithm
132 if (i % 1000 == 0) {
134 }
135
137 }
139
140 this->logErrors();
141 this->setProperty("OutputWorkspace", m_outputWS);
142}
143
155void He3TubeEfficiency::correctForEfficiency(std::size_t spectraIndex, const API::SpectrumInfo &spectrumInfo) {
156 if (spectrumInfo.isMonitor(spectraIndex) || spectrumInfo.isMasked(spectraIndex)) {
157 return;
158 }
159
160 const auto &det = spectrumInfo.detector(spectraIndex);
161 const double exp_constant = this->calculateExponential(spectraIndex, det);
162 const double scale = this->getProperty("ScaleFactor");
163
164 const auto &yValues = m_inputWS->y(spectraIndex);
165 auto &yOut = m_outputWS->mutableY(spectraIndex);
166
167 const auto wavelength = m_inputWS->points(spectraIndex);
168
169 // todo is it going to die? returning local var by reference
170 std::vector<double> effCorrection(wavelength.size());
171
172 computeEfficiencyCorrection(effCorrection, wavelength, exp_constant, scale);
173
174 std::transform(yValues.cbegin(), yValues.cend(), effCorrection.cbegin(), yOut.begin(),
175 [&](const double y, const double effcorr) { return y * effcorr; });
176
177 const auto &eValues = m_inputWS->e(spectraIndex);
178 auto &eOut = m_outputWS->mutableE(spectraIndex);
179
180 // reset wavelength iterator
181 std::transform(eValues.cbegin(), eValues.cend(), effCorrection.cbegin(), eOut.begin(),
182 [&](const double e, const double effcorr) { return e * effcorr; });
183}
184
193double He3TubeEfficiency::calculateExponential(std::size_t spectraIndex, const Geometry::IDetector &idet) {
194 // Get the parameters for the current associated tube
195 double pressure = this->getParameter("TubePressure", spectraIndex, "tube_pressure", idet);
196 double tubethickness = this->getParameter("TubeThickness", spectraIndex, "tube_thickness", idet);
197 double temperature = this->getParameter("TubeTemperature", spectraIndex, "tube_temperature", idet);
198
199 double detRadius(0.0);
200 Kernel::V3D detAxis;
201 this->getDetectorGeometry(idet, detRadius, detAxis);
202 double detDiameter = 2.0 * detRadius;
203 double twiceTubeThickness = 2.0 * tubethickness;
204
205 // now get the sin of the angle, it's the magnitude of the cross product of
206 // unit vector along the detector tube axis and a unit vector directed from
207 // the sample to the detector center
208 const Kernel::V3D vectorFromSample = normalize(idet.getPos() - m_samplePos);
209 Kernel::Quat rot = idet.getRotation();
210 // rotate the original cylinder object axis to get the detector axis in the
211 // actual instrument
212 rot.rotate(detAxis);
213 detAxis.normalize();
214 // Scalar product is quicker than cross product
215 double cosTheta = detAxis.scalar_prod(vectorFromSample);
216 double sinTheta = std::sqrt(1.0 - cosTheta * cosTheta);
217
218 const double straight_path = detDiameter - twiceTubeThickness;
219 if (std::fabs(straight_path - 0.0) < TOL) {
220 throw std::out_of_range("Twice tube thickness cannot be greater than "
221 "or equal to the tube diameter");
222 }
223
224 const double pathlength = straight_path / sinTheta;
225 return EXP_SCALAR_CONST * (pressure / temperature) * pathlength;
226}
227
234void He3TubeEfficiency::getDetectorGeometry(const Geometry::IDetector &det, double &detRadius, Kernel::V3D &detAxis) {
235 std::shared_ptr<const Geometry::IObject> shape_sptr = det.shape();
236 if (!shape_sptr) {
237 throw std::runtime_error("Detector geometry error: detector with id: " + std::to_string(det.getID()) +
238 " does not have shape. Is this a detectors group?\n"
239 "The algorithm works for instruments with one-to-one "
240 "spectra-to-detector maps only!");
241 }
242 // std::map<const Geometry::Object *, std::pair<double,
243 // Kernel::V3D>>::const_iterator
244 auto it = m_shapeCache.find(shape_sptr.get());
245 if (it == m_shapeCache.end()) {
246 double xDist = distToSurface(Kernel::V3D(DIST_TO_UNIVERSE_EDGE, 0, 0), shape_sptr.get());
247 double zDist = distToSurface(Kernel::V3D(0, 0, DIST_TO_UNIVERSE_EDGE), shape_sptr.get());
248 if (std::abs(zDist - xDist) < 1e-8) {
249 detRadius = zDist / 2.0;
250 detAxis = Kernel::V3D(0, 1, 0);
251 // assume radii in z and x and the axis is in the y
252 PARALLEL_CRITICAL(deteff_shapecachea) { m_shapeCache.insert({shape_sptr.get(), {detRadius, detAxis}}); }
253 return;
254 }
255 double yDist = distToSurface(Kernel::V3D(0, DIST_TO_UNIVERSE_EDGE, 0), shape_sptr.get());
256 if (std::abs(yDist - zDist) < 1e-8) {
257 detRadius = yDist / 2.0;
258 detAxis = Kernel::V3D(1, 0, 0);
259 // assume that y and z are radii of the cylinder's circular cross-section
260 // and the axis is perpendicular, in the x direction
261 PARALLEL_CRITICAL(deteff_shapecacheb) { m_shapeCache.insert({shape_sptr.get(), {detRadius, detAxis}}); }
262 return;
263 }
264
265 if (std::abs(xDist - yDist) < 1e-8) {
266 detRadius = xDist / 2.0;
267 detAxis = Kernel::V3D(0, 0, 1);
268 PARALLEL_CRITICAL(deteff_shapecachec) { m_shapeCache.insert({shape_sptr.get(), {detRadius, detAxis}}); }
269 return;
270 }
271 } else {
272 std::pair<double, Kernel::V3D> geometry = it->second;
273 detRadius = geometry.first;
274 detAxis = geometry.second;
275 }
276}
277
289double He3TubeEfficiency::distToSurface(const Kernel::V3D &start, const Geometry::IObject *shape) const {
290 // get a vector from the point that was passed to the origin
291 const Kernel::V3D direction = normalize(-start);
292 // put the point and the vector (direction) together to get a line,
293 // here called a track
294 Geometry::Track track(start, direction);
295 // split the track (line) up into the part that is inside the shape and the
296 // part that is outside
297 shape->interceptSurface(track);
298
299 if (track.count() != 1) {
300 // the track missed the shape, probably the shape is not centered on
301 // the origin
302 throw std::invalid_argument("Fatal error interpreting the shape of a detector");
303 }
304 // the first part of the track will be the part inside the shape,
305 // return its length
306 return track.cbegin()->distInsideObject;
307}
308
316double He3TubeEfficiency::detectorEfficiency(const double alpha, const double scale_factor) const {
317 return (scale_factor / (1.0 - std::exp(-alpha)));
318}
319
324 std::vector<int>::size_type nspecs = m_spectraSkipped.size();
325 if (nspecs > 0) {
326 this->g_log.warning() << "There were " << nspecs << " spectra that could not be corrected. ";
327 this->g_log.debug() << "Unaffected spectra numbers: ";
328 for (size_t i = 0; i < nspecs; ++i) {
329 this->g_log.debug() << m_spectraSkipped[i] << " ";
330 }
331 this->g_log.debug() << '\n';
332 }
333}
334
344double He3TubeEfficiency::getParameter(const std::string &wsPropName, std::size_t currentIndex,
345 const std::string &detPropName, const Geometry::IDetector &idet) {
346 std::vector<double> wsProp = this->getProperty(wsPropName);
347
348 if (wsProp.empty()) {
349 return idet.getNumberParameter(detPropName).at(0);
350 } else {
351 if (wsProp.size() == 1) {
352 return wsProp.at(0);
353 } else {
354 return wsProp.at(currentIndex);
355 }
356 }
357}
358
363 this->g_log.information("Processing event workspace");
364
365 const API::MatrixWorkspace_const_sptr matrixInputWS = getProperty("InputWorkspace");
366
367 // generate the output workspace pointer
368 API::MatrixWorkspace_sptr matrixOutputWS = getProperty("OutputWorkspace");
369 if (matrixOutputWS != matrixInputWS) {
370 matrixOutputWS = matrixInputWS->clone();
371 setProperty("OutputWorkspace", matrixOutputWS);
372 }
373 auto m_outputWS = std::dynamic_pointer_cast<DataObjects::EventWorkspace>(matrixOutputWS);
374
375 std::size_t numHistograms = m_outputWS->getNumberHistograms();
376 auto &spectrumInfo = m_outputWS->mutableSpectrumInfo();
377 m_progress = std::make_unique<API::Progress>(this, 0.0, 1.0, numHistograms);
378
380 for (int i = 0; i < static_cast<int>(numHistograms); ++i) {
382
383 const auto &det = spectrumInfo.detector(i);
384 if (spectrumInfo.isMonitor(i) || spectrumInfo.isMasked(i)) {
385 continue;
386 }
387
388 double exp_constant = 0.0;
389 try {
390 exp_constant = this->calculateExponential(i, det);
391 } catch (std::out_of_range &) {
392 // Parameters are bad so skip correction
393 PARALLEL_CRITICAL(deteff_invalid) {
394 m_spectraSkipped.emplace_back(m_outputWS->getAxis(1)->spectraNo(i));
395 m_outputWS->getSpectrum(i).clearData();
396 spectrumInfo.setMasked(i, true);
397 }
398 }
399
400 // Do the correction
401 auto &evlist = m_outputWS->getSpectrum(i);
402 switch (evlist.getEventType()) {
403 case API::TOF:
404 // Switch to weights if needed.
405 evlist.switchTo(API::WEIGHTED);
406 // Fall through
407 case API::WEIGHTED:
408 eventHelper(evlist.getWeightedEvents(), exp_constant);
409 break;
411 eventHelper(evlist.getWeightedEventsNoTime(), exp_constant);
412 break;
413 }
414
415 m_progress->report();
416
417 // check for canceling the algorithm
418 if (i % 1000 == 0) {
420 }
421
423 }
425
426 m_outputWS->clearMRU();
427
428 this->logErrors();
429}
430
438void He3TubeEfficiency::computeEfficiencyCorrection(std::vector<double> &effCorrection, const Points &xes,
439 const double expConstant, const double scale) const {
440
441 std::transform(xes.cbegin(), xes.cend(), effCorrection.begin(),
442 [&](double wavelen) { return detectorEfficiency(expConstant * wavelen, scale); });
443}
444
450template <class T> void He3TubeEfficiency::eventHelper(std::vector<T> &events, double expval) {
451 const double scale = this->getProperty("ScaleFactor");
452
453 for (auto &event : events) {
454 auto de = static_cast<float>(this->detectorEfficiency(expval * event.tof(), scale));
455 event.m_weight *= de;
456 event.m_errorSquared *= de * de;
457 }
458}
459
460} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
const double DIST_TO_UNIVERSE_EDGE
const double TOL
const double EXP_SCALAR_CONST
#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_CRITICAL(name)
#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
void clear() override
Clears all properties under management.
Definition: Algorithm.cpp:2125
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
void interruption_point()
This is called during long-running operations, and check if the algorithm has requested that it be ca...
Definition: Algorithm.cpp:1687
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.
API::SpectrumInfo is an intermediate step towards a SpectrumInfo that is part of Instrument-2....
Definition: SpectrumInfo.h:53
bool isMonitor(const size_t index) const
Returns true if the detector(s) associated with the spectrum are monitors.
bool isMasked(const size_t index) const
Returns true if the detector(s) associated with the spectrum are masked.
const Geometry::IDetector & detector(const size_t index) const
Return a const reference to the detector or detector group of the spectrum with given index.
A property class for workspaces.
A validator which checks that the unit of the workspace referred to by a WorkspaceProperty is the exp...
Corrects the input workspace for helium3 tube efficiency based on an exponential parameterization.
double detectorEfficiency(const double alpha, const double scale_factor=1.0) const
Calculate the detector efficiency.
std::map< const Geometry::IObject *, std::pair< double, Kernel::V3D > > m_shapeCache
A lookup of previously seen shape objects used to save calculation time as most detectors have the sa...
double getParameter(const std::string &wsPropName, std::size_t currentIndex, const std::string &detPropName, const Geometry::IDetector &idet)
Retrieve the detector parameters from workspace or detector properties.
Kernel::V3D m_samplePos
Sample position.
void init() override
Declare algorithm properties.
void eventHelper(std::vector< T > &events, double expval)
Helper for event handling.
const Geometry::ParameterMap * m_paraMap
Map that stores additional properties for detectors.
API::MatrixWorkspace_const_sptr m_inputWS
The user selected (input) workspace.
std::vector< specnum_t > m_spectraSkipped
The spectra numbers that were skipped.
void correctForEfficiency(std::size_t spectraIndex, const API::SpectrumInfo &spectrumInfo)
Correct the given spectra index for efficiency.
void getDetectorGeometry(const Geometry::IDetector &det, double &detRadius, Kernel::V3D &detAxis)
Sets the detector geometry cache if necessary.
double calculateExponential(std::size_t spectraIndex, const Geometry::IDetector &idet)
Function to calculate exponential contribution.
double distToSurface(const Kernel::V3D &start, const Geometry::IObject *shape) const
Computes the distance to the given shape from a starting point.
void logErrors() const
Log any errors with spectra that occurred.
void computeEfficiencyCorrection(std::vector< double > &effCorrection, const HistogramData::Points &wavelength, const double expConstant, const double scale) const
Calculates the efficiency correction from the points.
void exec() override
Executes the algorithm.
API::MatrixWorkspace_sptr m_outputWS
The output workspace, maybe the same as the input one.
std::unique_ptr< API::Progress > m_progress
Algorithm progress keeper.
virtual Kernel::V3D getPos() const =0
Get the position of the IComponent. Tree structure is traverse through the.
virtual Kernel::Quat getRotation() const =0
Get the absolute orientation of the IComponent.
virtual std::vector< double > getNumberParameter(const std::string &pname, bool recursive=true) const =0
Get a parameter defined as a double.
Interface class for detector objects.
Definition: IDetector.h:43
virtual detid_t getID() const =0
Get the detector ID.
virtual const std::shared_ptr< const IObject > shape() const =0
Returns the shape of the Object.
IObject : Interface for geometry objects.
Definition: IObject.h:41
virtual int interceptSurface(Geometry::Track &) const =0
Defines a track as a start point and a direction.
Definition: Track.h:165
int count() const
Returns the number of links.
Definition: Track.h:219
LType::const_iterator cbegin() const
Returns an interator to the start of the set of links (const version)
Definition: Track.h:206
Support for a property that holds an array of values.
Definition: ArrayProperty.h:28
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 warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
The concrete, templated class for properties.
Class for quaternions.
Definition: Quat.h:39
void rotate(V3D &) const
Rotate a vector.
Definition: Quat.cpp:397
Class for 3D vectors.
Definition: V3D.h:34
constexpr double scalar_prod(const V3D &v) const noexcept
Calculates the cross product.
Definition: V3D.h:274
double normalize()
Make a normalized vector (return norm value)
Definition: V3D.cpp:130
std::shared_ptr< const MatrixWorkspace > MatrixWorkspace_const_sptr
shared pointer to the matrix workspace base class (const version)
@ WEIGHTED_NOTIME
Definition: IEventList.h:18
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
MANTID_KERNEL_DLL V3D normalize(V3D v)
Normalizes a V3D.
Definition: V3D.h:341
std::string to_string(const wide_integer< Bits, Signed > &n)
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54