Mantid
Loading...
Searching...
No Matches
MayersSampleCorrection.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 +
9#include "MantidAPI/Sample.h"
21
22#include <cmath>
23
24namespace Mantid::Algorithms {
25
28using Kernel::Direction;
29using Kernel::V3D;
30
31// Register the algorithm into the AlgorithmFactory
32DECLARE_ALGORITHM(MayersSampleCorrection)
33
34
38
40const std::string MayersSampleCorrection::name() const { return "MayersSampleCorrection"; }
41
43int MayersSampleCorrection::version() const { return 1; }
44
46const std::string MayersSampleCorrection::category() const { return "CorrectionFunctions\\AbsorptionCorrections"; }
47
49const std::string MayersSampleCorrection::summary() const {
50 return "Corrects the input data for the effects of attenuation & multiple "
51 "scattering";
52}
53
58 // Inputs
60 std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input, createInputWSValidator()),
61 "Input workspace with X units in TOF. The workspace must "
62 "also have a sample with a cylindrical shape and an "
63 "instrument with a defined source and sample position.");
64 declareProperty("MultipleScattering", false,
65 "If True then also correct for the effects of multiple scattering."
66 "Please note that the MS correction assumes the scattering is elastic.",
68 declareProperty("MSEvents", 10000,
69 "Controls the number of second-scatter "
70 "events generated. Only applicable where "
71 "MultipleScattering=True.",
73 declareProperty("MSRuns", 10,
74 "Controls the number of simulations, each containing "
75 "MSEvents, performed. The final MS correction is "
76 "computed as the average over the runs. Only applicable"
77 "where MultipleScattering=True.",
79 // Outputs
80 declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
81 "An output workspace.");
82}
83
85 using API::Progress;
87
88 MatrixWorkspace_sptr inputWS = getProperty("InputWorkspace");
89 const bool mscatOn = getProperty("MultipleScattering");
90 const int msEvents = getProperty("MSEvents");
91 const int msRuns = getProperty("MSRuns");
92
93 MatrixWorkspace_sptr outputWS = WorkspaceFactory::Instance().create(inputWS);
94 // Instrument constants
95 auto instrument = inputWS->getInstrument();
96 const auto frame = instrument->getReferenceFrame();
97
98 // Sample
99 const auto &sampleShape = inputWS->sample().getShape();
100 // Current Object code computes quite an inaccurate bounding box so we do
101 // something better for the time being
102 const double big(100.); // seems to be a sweet spot...
103 double minX(-big), maxX(big), minY(-big), maxY(big), minZ(-big), maxZ(big);
104 sampleShape.getBoundingBox(maxX, maxY, maxZ, minX, minY, minZ);
105 V3D boxWidth(maxX - minX, maxY - minY, maxZ - minZ);
106 const double radius(0.5 * boxWidth[frame->pointingHorizontal()]), height(boxWidth[frame->pointingUp()]);
107 const auto &sampleMaterial = sampleShape.material();
108
109 const size_t nhist(inputWS->getNumberHistograms());
110 Progress prog(this, 0.0, 1.0, nhist);
111 prog.setNotifyStep(0.01);
112
113 const auto &spectrumInfo = inputWS->spectrumInfo();
114
115 PARALLEL_FOR_IF(Kernel::threadSafe(*inputWS, *outputWS))
116 for (int64_t i = 0; i < static_cast<int64_t>(nhist); ++i) {
118
119 if (!spectrumInfo.hasDetectors(i) || spectrumInfo.isMonitor(i) || spectrumInfo.isMasked(i) ||
120 spectrumInfo.l2(i) == 0) {
121 continue;
122 }
123
125 params.mscat = mscatOn;
126 params.l1 = spectrumInfo.l1();
127 params.l2 = spectrumInfo.l2(i);
128 params.twoTheta = spectrumInfo.twoTheta(i);
129 // Code requires angle above/below scattering plane not our definition of
130 // phi
131 const auto pos = spectrumInfo.position(i);
132 // Theta here is the angle between beam and neutron path not necessarily
133 // twoTheta if sample not at origin
134 double _(0.0), theta(0.0), phi(0.0);
135 pos.getSpherical(_, theta, phi);
136 params.azimuth = asin(sin(theta) * sin(phi));
137 params.rho = sampleMaterial.numberDensityEffective();
138 params.sigmaAbs = sampleMaterial.absorbXSection();
139 params.sigmaSc = sampleMaterial.totalScatterXSection();
140 params.cylRadius = radius;
141 params.cylHeight = height;
142 params.msNEvents = static_cast<size_t>(msEvents);
143 params.msNRuns = static_cast<size_t>(msRuns);
144
145 MayersSampleCorrectionStrategy correction(params, inputWS->histogram(i));
146 outputWS->setHistogram(i, correction.getCorrectedHisto());
147 prog.report();
148
150 }
152
153 setProperty("OutputWorkspace", outputWS);
154}
155
156//------------------------------------------------------------------------------
157// Private members
158//------------------------------------------------------------------------------
167 auto validator = std::make_shared<CompositeValidator>();
168
170 validator->add<InstrumentValidator, unsigned int>(requires);
171
173 validator->add<SampleValidator, unsigned int>(requires);
174
175 // NOTE: Mayers correction requires the input to be of TOF
176 validator->add<WorkspaceUnitValidator>("TOF");
177
178 return validator;
179}
180
181} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
double height
Definition: GetAllEi.cpp:155
#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.
double radius
Definition: Rasterize.cpp:31
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
A validator which checks that a workspace has a valid instrument.
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
A validator which checks that sample has the required properties.
A property class for workspaces.
A validator which checks that the unit of the workspace referred to by a WorkspaceProperty is the exp...
Applies the procedure found in section 4 of https://inis.iaea.org/search/search.aspx?...
Mantid::HistogramData::Histogram getCorrectedHisto()
Return the correction factors.
Corrects for the effects of absorption and multiple scattering using the algorithm of Jerry Mayers.
void exec() override
Virtual method - must be overridden by concrete algorithm.
void init() override
Initialize the algorithm's properties.
int version() const override
Algorithm's version for identification.
const std::string name() const override
Algorithms name for identification.
const std::string category() const override
Algorithm's category for identification.
Kernel::IValidator_sptr createInputWSValidator() const
const std::string summary() const override
Algorithm's summary for use in the GUI and help.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
Definition: ProgressBase.h:51
void setNotifyStep(double notifyStepPct)
Override the frequency at which notifications are sent out.
Manage the lifetime of a class intended to be a singleton.
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
Class for 3D vectors.
Definition: V3D.h:34
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
std::shared_ptr< const Mantid::Geometry::IDetector > IDetector_const_sptr
Shared pointer to IDetector (const version)
Definition: IDetector.h:102
std::shared_ptr< IValidator > IValidator_sptr
A shared_ptr to an IValidator.
Definition: IValidator.h:26
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
Stores parameters for a single calculation for a given angle and sample details.
size_t msNRuns
Number of runs to average ms correction over.
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54