Mantid
Loading...
Searching...
No Matches
DetectorEfficiencyCor.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"
11#include "MantidAPI/Run.h"
23#include "MantidTypes/SpectrumDefinition.h"
24
25#include <algorithm>
26#include <cmath>
27#include <functional>
28
29namespace Mantid::Algorithms {
30// Register the class into the algorithm factory
31DECLARE_ALGORITHM(DetectorEfficiencyCor)
32
33using namespace Kernel;
34using namespace API;
35using namespace Geometry;
36using namespace DataObjects;
37
38namespace {
39
40// E = KSquaredToE*K^2 KSquaredToE = (hbar^2)/(2*NeutronMass)
41const double KSquaredToE = 2.07212466; // units of meV Angstrom^-2
42
43const short NUMCOEFS = 25;
44// series expansion coefficients copied from a fortran source code file
45const double c_eff_f[] = {
46 0.7648360390553052, -0.3700950778935237, 0.1582704090813516, -6.0170218669705407E-02,
47 2.0465515957968953E-02, -6.2690181465706840E-03, 1.7408667184745830E-03, -4.4101378999425122E-04,
48 1.0252117967127217E-04, -2.1988904738111659E-05, 4.3729347905629990E-06, -8.0998753944849788E-07,
49 1.4031240949230472E-07, -2.2815971698619819E-08, 3.4943984983382137E-09, -5.0562696807254781E-10,
50 6.9315483353094009E-11, -9.0261598195695569E-12, 1.1192324844699897E-12, -1.3204992654891612E-13,
51 1.4100387524251801E-14, -8.6430862467068437E-16, -1.1129985821867194E-16, -4.5505266221823604E-16,
52 3.8885561437496108E-16};
53
54const double c_eff_g[] = {
55 2.033429926215546, -2.3123407369310212E-02, 7.0671915734894875E-03, -7.5970017538257162E-04,
56 7.4848652541832373E-05, 4.5642679186460588E-05, -2.3097291253000307E-05, 1.9697221715275770E-06,
57 2.4115259271262346E-06, -7.1302220919333692E-07, -2.5124427621592282E-07, 1.3246884875139919E-07,
58 3.4364196805913849E-08, -2.2891359549026546E-08, -6.7281240212491156E-09, 3.8292458615085678E-09,
59 1.6451021034313840E-09, -5.5868962123284405E-10, -4.2052310689211225E-10, 4.3217612266666094E-11,
60 9.9547699528024225E-11, 1.2882834243832519E-11, -1.9103066351000564E-11, -7.6805495297094239E-12,
61 1.8568853399347773E-12};
62
63// constants from the fortran code multiplied together sigref=143.23d0,
64// wref=3.49416d0, atmref=10.0d0 const = 2.0*sigref*wref/atmref
65const double g_helium_prefactor = 2.0 * 143.23 * 3.49416 / 10.0;
66
67// this should be a big number but not so big that there are rounding errors
68const double DIST_TO_UNIVERSE_EDGE = 1e3;
69
70// Name of pressure parameter
71const std::string PRESSURE_PARAM = "TubePressure";
72// Name of wall thickness parameter
73const std::string THICKNESS_PARAM = "TubeThickness";
74} // namespace
75
76// this default constructor calls default constructors and sets other member
77// data to impossible (flag) values
79 : Algorithm(), m_inputWS(), m_outputWS(), m_paraMap(nullptr), m_Ei(-1.0), m_ki(-1.0), m_shapeCache(), m_samplePos(),
80 m_spectraSkipped() {
81 m_shapeCache.clear();
82}
83
88 auto val = std::make_shared<CompositeValidator>();
89 val->add<WorkspaceUnitValidator>("DeltaE");
90 val->add<HistogramValidator>();
91 val->add<InstrumentValidator>();
92 declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input, val),
93 "The workspace to correct for detector efficiency");
94 declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
95 "The name of the workspace in which to store the result. Each histogram "
96 "from the input workspace maps to a histogram in this workspace that has "
97 "just one value which indicates if there was a bad detector.");
98 auto checkEi = std::make_shared<BoundedValidator<double>>();
99 checkEi->setLower(0.0);
100 declareProperty("IncidentEnergy", EMPTY_DBL(), checkEi,
101 "The energy of neutrons leaving the source as can be "
102 "calculated by :ref:`algm-GetEi`. If this value is provided, "
103 "uses property value, if it is not present, needs Ei log "
104 "value set on the workspace.");
105}
106
114 // gets and checks the values passed to the algorithm
116
117 // wave number that the neutrons originally had
118 m_ki = std::sqrt(m_Ei / KSquaredToE);
119
120 // Store some information about the instrument setup that will not change
121 m_samplePos = m_inputWS->getInstrument()->getSample()->getPos();
122
123 int64_t numHists = m_inputWS->getNumberHistograms();
124 auto numHists_d = static_cast<double>(numHists);
125 const auto progStep = static_cast<int64_t>(ceil(numHists_d / 100.0));
126 auto &spectrumInfo = m_inputWS->spectrumInfo();
127
129 for (int64_t i = 0; i < numHists; ++i) {
131
132 m_outputWS->setSharedX(i, m_inputWS->sharedX(i));
133 try {
134 correctForEfficiency(i, spectrumInfo);
135 } catch (Exception::NotFoundError &) {
136 // zero the Y data that can't be corrected
137 m_outputWS->mutableY(i) *= 0.0;
138 PARALLEL_CRITICAL(deteff_invalid) {
139 m_spectraSkipped.insert(m_spectraSkipped.end(), m_inputWS->getAxis(1)->spectraNo(i));
140 }
141 }
142 // make regular progress reports and check for canceling the algorithm
143 if (i % progStep == 0) {
144 progress(static_cast<double>(i) / numHists_d);
146 }
147
149 }
151
152 logErrors(numHists);
153 setProperty("OutputWorkspace", m_outputWS);
154}
161 // these first three properties are fully checked by validators
162 m_inputWS = getProperty("InputWorkspace");
163 m_paraMap = &(m_inputWS->constInstrumentParameters());
164
165 m_Ei = getProperty("IncidentEnergy");
166 // If we're not given an Ei, see if one has been set.
167 if (m_Ei == EMPTY_DBL()) {
168 if (m_inputWS->run().hasProperty("Ei")) {
169 m_Ei = m_inputWS->run().getPropertyValueAsType<double>("Ei");
170 g_log.debug() << "Using stored Ei value " << m_Ei << "\n";
171 } else {
172 throw std::invalid_argument("No Ei value has been set or stored within the run information.");
173 }
174 }
175
176 m_outputWS = getProperty("OutputWorkspace");
177 // If input and output workspaces are not the same, create a new workspace for
178 // the output
179 if (m_outputWS != m_inputWS) {
180 m_outputWS = create<MatrixWorkspace>(*m_inputWS);
181 }
182}
183
195void DetectorEfficiencyCor::correctForEfficiency(int64_t spectraIn, const SpectrumInfo &spectrumInfo) {
196 if (!spectrumInfo.hasDetectors(spectraIn))
197 throw Exception::NotFoundError("No detectors found", spectraIn);
198
199 if (spectrumInfo.isMonitor(spectraIn) || spectrumInfo.isMasked(spectraIn)) {
200 return;
201 }
202
203 auto &yout = m_outputWS->mutableY(spectraIn);
204 auto &eout = m_outputWS->mutableE(spectraIn);
205 // Need the original values so this is not a reference
206 const auto yValues = m_inputWS->y(spectraIn);
207 const auto eValues = m_inputWS->e(spectraIn);
208
209 // Storage for the reciprocal wave vectors that are calculated as the
210 // correction proceeds
211 std::vector<double> oneOverWaveVectors(yValues.size());
212 const auto &detectorInfo = m_inputWS->detectorInfo();
213 const auto &spectrumDefinition = spectrumInfo.spectrumDefinition(spectraIn);
214
215 for (const auto &index : spectrumDefinition) {
216 const auto detIndex = index.first;
217 const auto &det_member = detectorInfo.detector(detIndex);
218 Parameter_sptr par = m_paraMap->getRecursive(det_member.getComponentID(), PRESSURE_PARAM);
219 if (!par) {
220 throw Exception::NotFoundError(PRESSURE_PARAM, spectraIn);
221 }
222 const double atms = par->value<double>();
223 par = m_paraMap->getRecursive(det_member.getComponentID(), THICKNESS_PARAM);
224 if (!par) {
225 throw Exception::NotFoundError(THICKNESS_PARAM, spectraIn);
226 }
227 const double wallThickness = par->value<double>();
228 double detRadius(0.0);
229 V3D detAxis;
230 getDetectorGeometry(det_member, detRadius, detAxis);
231
232 // now get the sin of the angle, it's the magnitude of the cross product of
233 // unit vector along the detector tube axis and a unit vector directed from
234 // the sample to the detector centre
235 const V3D vectorFromSample = normalize(det_member.getPos() - m_samplePos);
236 Quat rot = det_member.getRotation();
237 // rotate the original cylinder object axis to get the detector axis in the
238 // actual instrument
239 rot.rotate(detAxis);
240 detAxis.normalize();
241 // Scalar product is quicker than cross product
242 double cosTheta = detAxis.scalar_prod(vectorFromSample);
243 double sinTheta = std::sqrt(1.0 - cosTheta * cosTheta);
244 // Detector constant
245 const double det_const = g_helium_prefactor * (detRadius - wallThickness) * atms / sinTheta;
246
247 auto yinItr = yValues.cbegin();
248 auto einItr = eValues.cbegin();
249 auto youtItr = yout.begin();
250 auto eoutItr = eout.begin();
251 auto xItr = m_inputWS->x(spectraIn).cbegin();
252 auto wavItr = oneOverWaveVectors.begin();
253
254 for (; youtItr != yout.end(); ++youtItr, ++eoutItr) {
255 if (index == spectrumDefinition[0]) {
256 *youtItr = 0.0;
257 *eoutItr = 0.0;
258 *wavItr = calculateOneOverK(*xItr, *(xItr + 1));
259 }
260 const double oneOverWave = *wavItr;
261 const auto nDets(static_cast<double>(spectrumDefinition.size()));
262 const double factor = 1.0 / nDets / detectorEfficiency(det_const * oneOverWave);
263 *youtItr += (*yinItr) * factor;
264 *eoutItr += (*einItr) * factor;
265 ++yinItr;
266 ++einItr;
267 ++xItr;
268 ++wavItr;
269 }
270 }
271}
272
282double DetectorEfficiencyCor::calculateOneOverK(double loBinBound, double uppBinBound) const {
283 double energy = m_Ei - 0.5 * (uppBinBound + loBinBound);
284 double oneOverKSquared = KSquaredToE / energy;
285 return std::sqrt(oneOverKSquared);
286}
287
293void DetectorEfficiencyCor::getDetectorGeometry(const Geometry::IDetector &det, double &detRadius, V3D &detAxis) {
294 std::shared_ptr<const IObject> shape_sptr = det.shape();
295 if (!shape_sptr->hasValidShape()) {
296 throw Exception::NotFoundError("Shape", "Detector has no shape");
297 }
298
299 std::map<const Geometry::IObject *, std::pair<double, Kernel::V3D>>::const_iterator it =
300 m_shapeCache.find(shape_sptr.get());
301 if (it == m_shapeCache.end()) {
302 double xDist = distToSurface(V3D(DIST_TO_UNIVERSE_EDGE, 0, 0), shape_sptr.get());
303 double zDist = distToSurface(V3D(0, 0, DIST_TO_UNIVERSE_EDGE), shape_sptr.get());
304 if (std::abs(zDist - xDist) < 1e-8) {
305 detRadius = zDist / 2.0;
306 detAxis = V3D(0, 1, 0);
307 // assume radi in z and x and the axis is in the y
308 PARALLEL_CRITICAL(deteff_shapecachea) {
309 m_shapeCache.emplace(shape_sptr.get(), std::make_pair(detRadius, detAxis));
310 }
311 return;
312 }
313 double yDist = distToSurface(V3D(0, DIST_TO_UNIVERSE_EDGE, 0), shape_sptr.get());
314 if (std::abs(yDist - zDist) < 1e-8) {
315 detRadius = yDist / 2.0;
316 detAxis = V3D(1, 0, 0);
317 // assume that y and z are radi of the cylinder's circular cross-section
318 // and the axis is perpendicular, in the x direction
319 PARALLEL_CRITICAL(deteff_shapecacheb) {
320 m_shapeCache.emplace(shape_sptr.get(), std::make_pair(detRadius, detAxis));
321 }
322 return;
323 }
324
325 if (std::abs(xDist - yDist) < 1e-8) {
326 detRadius = xDist / 2.0;
327 detAxis = V3D(0, 0, 1);
328 PARALLEL_CRITICAL(deteff_shapecachec) {
329 m_shapeCache.emplace(shape_sptr.get(), std::make_pair(detRadius, detAxis));
330 }
331 return;
332 }
333 } else {
334 std::pair<double, V3D> geometry = it->second;
335 detRadius = geometry.first;
336 detAxis = geometry.second;
337 }
338}
339
351double DetectorEfficiencyCor::distToSurface(const V3D &start, const IObject *shape) const {
352 // get a vector from the point that was passed to the origin
353 const V3D direction = normalize(-start);
354 // put the point and the vector (direction) together to get a line, here
355 // called a track
356 Track track(start, direction);
357 // split the track (line) up into the part that is inside the shape and the
358 // part that is outside
359 shape->interceptSurface(track);
360
361 if (track.count() != 1) { // the track missed the shape, probably the shape is
362 // not centered on the origin
363 throw std::invalid_argument("Fatal error interpreting the shape of a detector");
364 }
365 // the first part of the track will be the part inside the shape, return its
366 // length
367 return track.cbegin()->distInsideObject;
368}
369
376double DetectorEfficiencyCor::detectorEfficiency(const double alpha) const {
377 if (alpha < 9.0) {
378 return 0.25 * M_PI * alpha * chebevApprox(0.0, 10.0, c_eff_f, alpha);
379 }
380 if (alpha > 10.0) {
381 double y = 1.0 - 18.0 / alpha;
382 return 1.0 - chebevApprox(-1.0, 1.0, c_eff_g, y) / (alpha * alpha);
383 }
384 double eff_f = 0.25 * M_PI * alpha * chebevApprox(0.0, 10.0, c_eff_f, alpha);
385 double y = 1.0 - 18.0 / alpha;
386 double eff_g = 1.0 - chebevApprox(-1.0, 1.0, c_eff_g, y) / (alpha * alpha);
387 return (10.0 - alpha) * eff_f + (alpha - 9.0) * eff_g;
388}
389
401double DetectorEfficiencyCor::chebevApprox(double a, double b, const double exspansionCoefs[], double x) const {
402 double d = 0.0;
403 double dd = 0.0;
404 double y = (2.0 * x - a - b) / (b - a);
405 double y2 = 2.0 * y;
406 for (int j = NUMCOEFS - 1; j > 0; j -= 1) {
407 double sv = d;
408 d = y2 * d - dd + exspansionCoefs[j];
409 dd = sv;
410 }
411 return y * d - dd + 0.5 * exspansionCoefs[0];
412}
413
419void DetectorEfficiencyCor::logErrors(size_t totalNDetectors) const {
420 std::vector<int>::size_type nspecs = m_spectraSkipped.size();
421 if (!m_spectraSkipped.empty()) {
422 g_log.warning() << "There were " << nspecs
423 << " spectra that could not be corrected out of total: " << totalNDetectors << '\n';
424 g_log.warning() << "Their spectra were nullified\n";
425 g_log.debug() << " Nullified spectra numbers: ";
426 auto itend = m_spectraSkipped.end();
427 for (auto it = m_spectraSkipped.begin(); it != itend; ++it) {
428 g_log.debug() << *it << " ";
429 }
430 g_log.debug() << "\n";
431 }
432}
433
434} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
double energy
Definition: GetAllEi.cpp:157
const double DIST_TO_UNIVERSE_EDGE
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_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
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 progress(double p, const std::string &msg="", double estimatedTime=0.0, int progressPrecision=0)
Sends ProgressNotification.
Definition: Algorithm.cpp:231
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 hasDetectors(const size_t index) const
Returns true if the spectrum is associated with detectors in the instrument.
bool isMasked(const size_t index) const
Returns true if the detector(s) associated with the spectrum are masked.
const SpectrumDefinition & spectrumDefinition(const size_t index) const
Returns a const reference to the SpectrumDefinition of the spectrum.
A property class for workspaces.
A validator which checks that the unit of the workspace referred to by a WorkspaceProperty is the exp...
void getDetectorGeometry(const Geometry::IDetector &det, double &detRadius, Kernel::V3D &detAxis)
Sets the detector geometry cache if necessary.
double calculateOneOverK(double loBinBound, double uppBinBound) const
Calculate one over the wave vector for 2 bin bounds.
double distToSurface(const Kernel::V3D &start, const Geometry::IObject *shape) const
Computes the distance to the given shape from a starting point.
double m_ki
stores the wave number of incidient neutrons, calculated from the energy
void exec() override
Executes the algorithm.
std::list< int64_t > m_spectraSkipped
The spectra numbers that were skipped.
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...
API::MatrixWorkspace_const_sptr m_inputWS
the user selected workspace
API::MatrixWorkspace_sptr m_outputWS
output workspace, maybe the same as the input one
double m_Ei
stores the user selected value for incidient energy of the neutrons
void init() override
Declare algorithm properties.
void retrieveProperties()
Retrieve algorithm properties.
double detectorEfficiency(const double alpha) const
Computes the detector efficiency for a given paramater.
double chebevApprox(double a, double b, const double exspansionCoefs[], double x) const
Computes an approximate expansion of a Chebysev polynomial.
const Geometry::ParameterMap * m_paraMap
points the map that stores additional properties for detectors in that map
void correctForEfficiency(int64_t spectraIn, const API::SpectrumInfo &spectrumInfo)
Correct the given spectra index for efficiency.
void logErrors(size_t totalNDetectors) const
Log any errors with spectra that occurred.
Interface class for detector objects.
Definition: IDetector.h:43
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
std::shared_ptr< Parameter > getRecursive(const IComponent *comp, const std::string &name, const std::string &type="") const
Use get() recursively to see if can find param in all parents of comp and given type (std::string ver...
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
Exception for when an item is not found in a collection.
Definition: Exception.h:145
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
Class for quaternions.
Definition: Quat.h:39
void rotate(V3D &) const
Rotate a vector.
Definition: Quat.cpp:397
std::vector< double > getRotation(bool check_normalisation=false, bool throw_on_errors=false) const
returns the rotation matrix defined by this quaternion as an 9-point
Definition: Quat.cpp:453
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< Parameter > Parameter_sptr
Typedef for the shared pointer.
Definition: Parameter.h:195
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
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
Definition: EmptyValues.h:43
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54