23#include "MantidTypes/SpectrumDefinition.h"
33using namespace Kernel;
35using namespace Geometry;
36using namespace DataObjects;
41const double KSquaredToE = 2.07212466;
43const short NUMCOEFS = 25;
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};
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};
65const double g_helium_prefactor = 2.0 * 143.23 * 3.49416 / 10.0;
71const std::string PRESSURE_PARAM =
"TubePressure";
73const std::string THICKNESS_PARAM =
"TubeThickness";
79 :
Algorithm(), m_inputWS(), m_outputWS(), m_paraMap(nullptr), m_Ei(-1.0), m_ki(-1.0), m_shapeCache(), m_samplePos(),
88 auto val = std::make_shared<CompositeValidator>();
93 "The workspace to correct for detector efficiency");
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);
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.");
118 m_ki = std::sqrt(
m_Ei / KSquaredToE);
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();
129 for (int64_t i = 0; i < numHists; ++i) {
143 if (i % progStep == 0) {
144 progress(
static_cast<double>(i) / numHists_d);
168 if (
m_inputWS->run().hasProperty(
"Ei")) {
169 m_Ei =
m_inputWS->run().getPropertyValueAsType<
double>(
"Ei");
172 throw std::invalid_argument(
"No Ei value has been set or stored within the run information.");
206 const auto yValues =
m_inputWS->y(spectraIn);
207 const auto eValues =
m_inputWS->e(spectraIn);
211 std::vector<double> oneOverWaveVectors(yValues.size());
212 const auto &detectorInfo =
m_inputWS->detectorInfo();
215 for (
const auto &
index : spectrumDefinition) {
216 const auto detIndex =
index.first;
217 const auto &det_member = detectorInfo.detector(detIndex);
222 const double atms = par->value<
double>();
227 const double wallThickness = par->value<
double>();
228 double detRadius(0.0);
242 double cosTheta = detAxis.
scalar_prod(vectorFromSample);
243 double sinTheta = std::sqrt(1.0 - cosTheta * cosTheta);
245 const double det_const = g_helium_prefactor * (detRadius - wallThickness) * atms / sinTheta;
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();
254 for (; youtItr != yout.end(); ++youtItr, ++eoutItr) {
255 if (
index == spectrumDefinition[0]) {
260 const double oneOverWave = *wavItr;
261 const auto nDets(
static_cast<double>(spectrumDefinition.size()));
263 *youtItr += (*yinItr) * factor;
264 *eoutItr += (*einItr) * factor;
283 double energy =
m_Ei - 0.5 * (uppBinBound + loBinBound);
284 double oneOverKSquared = KSquaredToE /
energy;
285 return std::sqrt(oneOverKSquared);
294 std::shared_ptr<const IObject> shape_sptr = det.
shape();
295 if (!shape_sptr->hasValidShape()) {
299 std::map<const Geometry::IObject *, std::pair<double, Kernel::V3D>>::const_iterator it =
304 if (std::abs(zDist - xDist) < 1e-8) {
305 detRadius = zDist / 2.0;
306 detAxis =
V3D(0, 1, 0);
309 m_shapeCache.emplace(shape_sptr.get(), std::make_pair(detRadius, detAxis));
314 if (std::abs(yDist - zDist) < 1e-8) {
315 detRadius = yDist / 2.0;
316 detAxis =
V3D(1, 0, 0);
320 m_shapeCache.emplace(shape_sptr.get(), std::make_pair(detRadius, detAxis));
325 if (std::abs(xDist - yDist) < 1e-8) {
326 detRadius = xDist / 2.0;
327 detAxis =
V3D(0, 0, 1);
329 m_shapeCache.emplace(shape_sptr.get(), std::make_pair(detRadius, detAxis));
334 std::pair<double, V3D> geometry = it->second;
335 detRadius = geometry.first;
336 detAxis = geometry.second;
356 Track track(start, direction);
361 if (track.
count() != 1) {
363 throw std::invalid_argument(
"Fatal error interpreting the shape of a detector");
367 return track.
cbegin()->distInsideObject;
378 return 0.25 * M_PI * alpha *
chebevApprox(0.0, 10.0, c_eff_f, alpha);
381 double y = 1.0 - 18.0 / alpha;
382 return 1.0 -
chebevApprox(-1.0, 1.0, c_eff_g,
y) / (alpha * alpha);
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;
404 double y = (2.0 *
x - a - b) / (b - a);
406 for (
int j = NUMCOEFS - 1; j > 0; j -= 1) {
408 d = y2 *
d - dd + exspansionCoefs[j];
411 return y *
d - dd + 0.5 * exspansionCoefs[0];
423 <<
" spectra that could not be corrected out of total: " << totalNDetectors <<
'\n';
425 g_log.
debug() <<
" Nullified spectra numbers: ";
#define DECLARE_ALGORITHM(classname)
const double DIST_TO_UNIVERSE_EDGE
std::map< DeltaEMode::Type, std::string > index
#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...
#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.
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
void progress(double p, const std::string &msg="", double estimatedTime=0.0, int progressPrecision=0)
Sends ProgressNotification.
void interruption_point()
This is called during long-running operations, and check if the algorithm has requested that it be ca...
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....
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...
Kernel::V3D m_samplePos
Sample position.
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.
virtual const std::shared_ptr< const IObject > shape() const =0
Returns the shape of the Object.
IObject : Interface for geometry objects.
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.
int count() const
Returns the number of links.
LType::const_iterator cbegin() const
Returns an interator to the start of the set of links (const version)
Exception for when an item is not found in a collection.
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.
void warning(const std::string &msg)
Logs at warning level.
void rotate(V3D &) const
Rotate a vector.
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
constexpr double scalar_prod(const V3D &v) const noexcept
Calculates the cross product.
double normalize()
Make a normalized vector (return norm value)
std::shared_ptr< Parameter > Parameter_sptr
Typedef for the shared pointer.
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.
MANTID_KERNEL_DLL V3D normalize(V3D v)
Normalizes a V3D.
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
@ Input
An input workspace.
@ Output
An output workspace.