34const double TOL = 1.0e-8;
38using namespace DataObjects;
39using namespace HistogramData;
45 :
Algorithm(), m_inputWS(), m_outputWS(), m_paraMap(
nullptr), m_shapeCache(), m_samplePos(), m_spectraSkipped(),
56 auto wsValidator = std::make_shared<CompositeValidator>();
62 "Name of the input workspace");
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);
69 "Constant factor with which to scale the calculated"
70 "detector efficiency. Same factor applies to all efficiencies.");
72 auto mustBePosArr = std::make_shared<Kernel::ArrayBoundedValidator<double>>();
73 mustBePosArr->setLower(0.0);
75 "Provide overriding the default tube pressure. The pressure must "
76 "be specified in atm.");
78 "Provide overriding the default tube thickness. The thickness must "
79 "be specified in metres.");
81 "Provide overriding the default tube temperature. The temperature must "
82 "be specified in Kelvin.");
105 std::dynamic_pointer_cast<const DataObjects::EventWorkspace>(
m_inputWS);
106 if (eventW !=
nullptr) {
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();
116 for (
int i = 0; i < static_cast<int>(numHists); ++i) {
122 }
catch (std::out_of_range &) {
156 if (spectrumInfo.
isMonitor(spectraIndex) || spectrumInfo.
isMasked(spectraIndex)) {
160 const auto &det = spectrumInfo.
detector(spectraIndex);
162 const double scale = this->
getProperty(
"ScaleFactor");
164 const auto &yValues =
m_inputWS->y(spectraIndex);
165 auto &yOut =
m_outputWS->mutableY(spectraIndex);
167 const auto wavelength =
m_inputWS->points(spectraIndex);
170 std::vector<double> effCorrection(wavelength.size());
174 std::transform(yValues.cbegin(), yValues.cend(), effCorrection.cbegin(), yOut.begin(),
175 [&](
const double y,
const double effcorr) { return y * effcorr; });
177 const auto &eValues =
m_inputWS->e(spectraIndex);
178 auto &eOut =
m_outputWS->mutableE(spectraIndex);
181 std::transform(eValues.cbegin(), eValues.cend(), effCorrection.cbegin(), eOut.begin(),
182 [&](
const double e,
const double effcorr) { return e * effcorr; });
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);
199 double detRadius(0.0);
202 double detDiameter = 2.0 * detRadius;
203 double twiceTubeThickness = 2.0 * tubethickness;
215 double cosTheta = detAxis.
scalar_prod(vectorFromSample);
216 double sinTheta = std::sqrt(1.0 - cosTheta * cosTheta);
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");
224 const double pathlength = straight_path / sinTheta;
235 std::shared_ptr<const Geometry::IObject> shape_sptr = det.
shape();
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!");
248 if (std::abs(zDist - xDist) < 1e-8) {
249 detRadius = zDist / 2.0;
256 if (std::abs(yDist - zDist) < 1e-8) {
257 detRadius = yDist / 2.0;
265 if (std::abs(xDist - yDist) < 1e-8) {
266 detRadius = xDist / 2.0;
272 std::pair<double, Kernel::V3D> geometry = it->second;
273 detRadius = geometry.first;
274 detAxis = geometry.second;
299 if (track.
count() != 1) {
302 throw std::invalid_argument(
"Fatal error interpreting the shape of a detector");
306 return track.
cbegin()->distInsideObject;
317 return (scale_factor / (1.0 - std::exp(-alpha)));
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) {
346 std::vector<double> wsProp = this->
getProperty(wsPropName);
348 if (wsProp.empty()) {
351 if (wsProp.size() == 1) {
354 return wsProp.at(currentIndex);
369 if (matrixOutputWS != matrixInputWS) {
370 matrixOutputWS = matrixInputWS->clone();
373 auto m_outputWS = std::dynamic_pointer_cast<DataObjects::EventWorkspace>(matrixOutputWS);
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);
380 for (
int i = 0; i < static_cast<int>(numHistograms); ++i) {
383 const auto &det = spectrumInfo.detector(i);
384 if (spectrumInfo.isMonitor(i) || spectrumInfo.isMasked(i)) {
388 double exp_constant = 0.0;
391 }
catch (std::out_of_range &) {
396 spectrumInfo.setMasked(i,
true);
402 switch (evlist.getEventType()) {
408 eventHelper(evlist.getWeightedEvents(), exp_constant);
411 eventHelper(evlist.getWeightedEventsNoTime(), exp_constant);
439 const double expConstant,
const double scale)
const {
441 std::transform(xes.cbegin(), xes.cend(), effCorrection.begin(),
442 [&](
double wavelen) { return detectorEfficiency(expConstant * wavelen, scale); });
451 const double scale = this->
getProperty(
"ScaleFactor");
453 for (
auto &event : events) {
455 event.m_weight *= de;
456 event.m_errorSquared *= de * de;
#define DECLARE_ALGORITHM(classname)
const double DIST_TO_UNIVERSE_EDGE
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...
#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.
void clear() override
Clears all properties under management.
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
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 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.
void execEvent()
Execute for events.
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.
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.
virtual int interceptSurface(Geometry::Track &) const =0
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)
Support for a property that holds an array of values.
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 information(const std::string &msg)
Logs at information level.
The concrete, templated class for properties.
void rotate(V3D &) const
Rotate a vector.
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< 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.
MANTID_KERNEL_DLL V3D normalize(V3D v)
Normalizes a V3D.
std::string to_string(const wide_integer< Bits, Signed > &n)
@ Input
An input workspace.
@ Output
An output workspace.