Mantid
Loading...
Searching...
No Matches
MonteCarloAbsorption.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"
29
30using namespace Mantid::API;
31using namespace Mantid::Geometry;
32using namespace Mantid::Kernel;
35
37namespace {
38
39constexpr int DEFAULT_NEVENTS = 1000;
40constexpr int DEFAULT_SEED = 123456789;
41constexpr int DEFAULT_LATITUDINAL_DETS = 5;
42constexpr int DEFAULT_LONGITUDINAL_DETS = 10;
43
45inline double toWavelength(double energy) {
46 static const double factor =
48 return factor / sqrt(energy);
49}
50
51struct EFixedProvider {
52 explicit EFixedProvider(const ExperimentInfo &expt) : m_expt(expt), m_emode(expt.getEMode()), m_value(0.0) {
53 if (m_emode == DeltaEMode::Direct) {
54 m_value = m_expt.getEFixed();
55 }
56 }
57 inline DeltaEMode::Type emode() const { return m_emode; }
58 inline double value(const Mantid::detid_t detID) const {
59 if (m_emode != DeltaEMode::Indirect)
60 return m_value;
61 else
62 return m_expt.getEFixed(detID);
63 }
64
65private:
66 const ExperimentInfo &m_expt;
67 const DeltaEMode::Type m_emode;
68 double m_value;
69};
70} // namespace
72
73namespace Mantid::Algorithms {
74
75DECLARE_ALGORITHM(MonteCarloAbsorption)
76
77//------------------------------------------------------------------------------
78// Private methods
79//------------------------------------------------------------------------------
80
81
85 // The input workspace must have an instrument and units of wavelength
86 auto wsValidator = std::make_shared<CompositeValidator>();
87 wsValidator->add<WorkspaceUnitValidator>("Wavelength");
88 wsValidator->add<InstrumentValidator>();
89
90 declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input, wsValidator),
91 "The name of the input workspace. The input workspace must "
92 "have X units of wavelength.");
93 declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
94 "The name to use for the output workspace.");
95
96 auto positiveInt = std::make_shared<Kernel::BoundedValidator<int>>();
97 positiveInt->setLower(1);
98 declareProperty("NumberOfWavelengthPoints", EMPTY_INT(), positiveInt,
99 "The number of wavelength points for which a simulation is "
100 "attempted if ResimulateTracksForDifferentWavelengths=true");
101 declareProperty("EventsPerPoint", DEFAULT_NEVENTS, positiveInt,
102 "The number of \"neutron\" events to generate per simulated point");
103 declareProperty("SeedValue", DEFAULT_SEED, positiveInt, "Seed the random number generator with this value");
104
105 auto interpolateOpt = createInterpolateOption();
106 declareProperty(interpolateOpt->property(), interpolateOpt->propertyDoc());
107 declareProperty("SparseInstrument", false,
108 "Enable simulation on special "
109 "instrument with a sparse grid of "
110 "detectors interpolating the "
111 "results to the real instrument.");
112 auto threeOrMore = std::make_shared<Kernel::BoundedValidator<int>>();
113 threeOrMore->setLower(3);
114 declareProperty("NumberOfDetectorRows", DEFAULT_LATITUDINAL_DETS, threeOrMore,
115 "Number of detector rows in the detector grid of the sparse instrument.");
116 setPropertySettings("NumberOfDetectorRows",
117 std::make_unique<EnabledWhenProperty>("SparseInstrument", ePropertyCriterion::IS_NOT_DEFAULT));
118 auto twoOrMore = std::make_shared<Kernel::BoundedValidator<int>>();
119 twoOrMore->setLower(2);
120 declareProperty("NumberOfDetectorColumns", DEFAULT_LONGITUDINAL_DETS, twoOrMore,
121 "Number of detector columns in the detector grid "
122 "of the sparse instrument.");
123 setPropertySettings("NumberOfDetectorColumns",
124 std::make_unique<EnabledWhenProperty>("SparseInstrument", ePropertyCriterion::IS_NOT_DEFAULT));
125
126 // Control the number of attempts made to generate a random point in the
127 // object
128 declareProperty("MaxScatterPtAttempts", 5000, positiveInt,
129 "Maximum number of tries made to generate a scattering point "
130 "within the sample (+ optional container etc). Objects with "
131 "holes in them, e.g. a thin annulus can cause problems "
132 "if this number is too low.\n"
133 "If a scattering point cannot be generated by increasing "
134 "this value then there is most likely a problem with "
135 "the sample geometry.");
136 declareProperty("ResimulateTracksForDifferentWavelengths", false, "Resimulate tracks for each wavelength point.");
137 setPropertySettings("NumberOfWavelengthPoints",
138 std::make_unique<EnabledWhenProperty>("ResimulateTracksForDifferentWavelengths",
139 ePropertyCriterion::IS_NOT_DEFAULT));
140 auto scatteringOptionValidator = std::make_shared<StringListValidator>();
141 scatteringOptionValidator->addAllowedValue("SampleAndEnvironment");
142 scatteringOptionValidator->addAllowedValue("SampleOnly");
143 scatteringOptionValidator->addAllowedValue("EnvironmentOnly");
144 declareProperty("SimulateScatteringPointIn", "SampleAndEnvironment",
145 "Simulate the scattering point in the vicinity of the sample or its "
146 "environment or both (default).",
147 scatteringOptionValidator);
148}
149
154 const MatrixWorkspace_sptr inputWS = getProperty("InputWorkspace");
155 const int nevents = getProperty("EventsPerPoint");
156 const bool resimulateTracks = getProperty("ResimulateTracksForDifferentWavelengths");
157 const int seed = getProperty("SeedValue");
158 InterpolationOption interpolateOpt;
159 interpolateOpt.set(getPropertyValue("Interpolation"), true, resimulateTracks);
160 const bool useSparseInstrument = getProperty("SparseInstrument");
161 const int maxScatterPtAttempts = getProperty("MaxScatterPtAttempts");
163 const auto pointsInProperty = getPropertyValue("SimulateScatteringPointIn");
164 if (pointsInProperty == "SampleOnly") {
166 } else if (pointsInProperty == "EnvironmentOnly") {
168 }
169 auto outputWS = doSimulation(*inputWS, static_cast<size_t>(nevents), resimulateTracks, seed, interpolateOpt,
170 useSparseInstrument, static_cast<size_t>(maxScatterPtAttempts), simulatePointsIn);
171 setProperty("OutputWorkspace", std::move(outputWS));
172}
173
178std::map<std::string, std::string> MonteCarloAbsorption::validateInputs() {
179 std::map<std::string, std::string> issues;
180 const bool resimulateTracksForDiffWavelengths = getProperty("ResimulateTracksForDifferentWavelengths");
181 // Only interpolate between wavelength points if resimulating tracks
182 if (resimulateTracksForDiffWavelengths) {
183 const int nlambda = getProperty("NumberOfWavelengthPoints");
184 InterpolationOption interpOpt;
185 const std::string interpValue = getPropertyValue("Interpolation");
186 interpOpt.set(interpValue, true, resimulateTracksForDiffWavelengths);
187 const auto nlambdaIssue = interpOpt.validateInputSize(nlambda);
188 if (!nlambdaIssue.empty()) {
189 issues["NumberOfWavelengthPoints"] = nlambdaIssue;
190 }
191 }
192 return issues;
193}
194
204std::shared_ptr<IMCInteractionVolume>
205MonteCarloAbsorption::createInteractionVolume(const API::Sample &sample, const size_t maxScatterPtAttempts,
207 auto interactionVol = std::make_shared<MCInteractionVolume>(sample, maxScatterPtAttempts, pointsIn);
208 return interactionVol;
209}
210
225std::shared_ptr<IMCAbsorptionStrategy>
227 Kernel::DeltaEMode::Type EMode, const size_t nevents,
228 const size_t maxScatterPtAttempts, const bool regenerateTracksForEachLambda) {
229 auto MCAbs = std::make_shared<MCAbsorptionStrategy>(interactionVol, beamProfile, EMode, nevents, maxScatterPtAttempts,
230 regenerateTracksForEachLambda);
231 return MCAbs;
232}
233
243std::shared_ptr<SparseWorkspace> MonteCarloAbsorption::createSparseWorkspace(const API::MatrixWorkspace &modelWS,
244 const size_t wavelengthPoints,
245 const size_t rows, const size_t columns) {
246 auto sparseWS = std::make_shared<SparseWorkspace>(modelWS, wavelengthPoints, rows, columns);
247 return sparseWS;
248}
249
255std::unique_ptr<InterpolationOption> MonteCarloAbsorption::createInterpolateOption() {
256 auto interpolationOpt = std::make_unique<InterpolationOption>();
257 return interpolationOpt;
258}
259
275 const bool resimulateTracksForDiffWavelengths, const int seed,
276 const InterpolationOption &interpolateOpt,
277 const bool useSparseInstrument,
278 const size_t maxScatterPtAttempts,
280 auto outputWS = createOutputWorkspace(inputWS);
281 const auto inputNbins = static_cast<int>(inputWS.blocksize());
282
283 int nlambda;
284 if (resimulateTracksForDiffWavelengths) {
285 nlambda = getProperty("NumberOfWavelengthPoints");
286 if (isEmpty(nlambda) || nlambda > inputNbins) {
287 if (!isEmpty(nlambda)) {
288 g_log.warning() << "The requested number of wavelength points is larger "
289 "than the spectra size. "
290 "Defaulting to spectra size.\n";
291 }
292 nlambda = inputNbins;
293 }
294 } else {
295 nlambda = inputNbins;
296 }
297 SparseWorkspace_sptr sparseWS;
298 if (useSparseInstrument) {
299 const int latitudinalDets = getProperty("NumberOfDetectorRows");
300 const int longitudinalDets = getProperty("NumberOfDetectorColumns");
301 sparseWS = createSparseWorkspace(inputWS, nlambda, latitudinalDets, longitudinalDets);
302 }
303 MatrixWorkspace &simulationWS = useSparseInstrument ? *sparseWS : *outputWS;
304 const MatrixWorkspace &instrumentWS = useSparseInstrument ? simulationWS : inputWS;
305 // Cache information about the workspace that will be used repeatedly
306 auto instrument = instrumentWS.getInstrument();
307 const auto nhists = static_cast<int64_t>(instrumentWS.getNumberHistograms());
308
309 EFixedProvider efixed(instrumentWS);
310 auto beamProfile = BeamProfileFactory::createBeamProfile(*instrument, inputWS.sample());
311
312 // Configure progress
313 Progress prog(this, 0.0, 1.0, nhists);
314 prog.setNotifyStep(0.01);
315 const std::string reportMsg = "Computing corrections";
316
317 // Configure strategy
318 auto interactionVolume = createInteractionVolume(inputWS.sample(), maxScatterPtAttempts, pointsIn);
319 auto strategy = createStrategy(*interactionVolume, *beamProfile, efixed.emode(), nevents, maxScatterPtAttempts,
320 resimulateTracksForDiffWavelengths);
321
322 const auto &spectrumInfo = simulationWS.spectrumInfo();
323
325 for (int64_t i = 0; i < nhists; ++i) {
327
328 auto &outE = simulationWS.mutableE(i);
329 // The input was cloned so clear the errors out
330 outE = 0.0;
331
332 if (!spectrumInfo.hasDetectors(i) || spectrumInfo.isMasked(i)) {
333 continue;
334 }
335 // Per spectrum values
336 const auto &detPos = spectrumInfo.position(i);
337 const double lambdaFixed = toWavelength(efixed.value(spectrumInfo.detector(i).getID()));
338 MersenneTwister rng(seed + int(i));
339
340 const auto lambdas = simulationWS.points(i).rawData();
341
342 const auto nbins = lambdas.size();
343 const size_t lambdaStepSize = nbins / nlambda;
344
345 std::vector<double> packedLambdas;
346 std::vector<double> packedAttFactors;
347 std::vector<double> packedAttFactorErrors;
348
349 for (size_t j = 0; j < nbins; j += lambdaStepSize) {
350 packedLambdas.push_back(lambdas[j]);
351 packedAttFactors.push_back(0);
352 packedAttFactorErrors.push_back(0);
353 // Ensure we have the last point for the interpolation
354 if (lambdaStepSize > 1 && j + lambdaStepSize >= nbins && j + 1 != nbins) {
355 j = nbins - lambdaStepSize - 1;
356 }
357 }
358 MCInteractionStatistics detStatistics(spectrumInfo.detector(i).getID(), inputWS.sample());
359
360 strategy->calculate(rng, detPos, packedLambdas, lambdaFixed, packedAttFactors, packedAttFactorErrors,
361 detStatistics);
362
363 if (g_log.is(Kernel::Logger::Priority::PRIO_DEBUG)) {
364 g_log.debug(detStatistics.generateScatterPointStats());
365 }
366
367 for (size_t j = 0; j < packedLambdas.size(); j++) {
368 auto idx = simulationWS.yIndexOfX(packedLambdas[j], i);
369 simulationWS.getSpectrum(i).dataY()[idx] = packedAttFactors[j];
370 simulationWS.getSpectrum(i).dataE()[idx] = packedAttFactorErrors[j];
371 }
372
373 // Interpolate through points not simulated. Simulation WS only has
374 // reduced X values if using sparse instrument so no interpolation required
375
376 if (!useSparseInstrument && lambdaStepSize > 1) {
377 auto histnew = simulationWS.histogram(i);
378
379 if (lambdaStepSize < nbins) {
380 interpolateOpt.applyInplace(histnew, lambdaStepSize);
381 } else {
382 std::fill(histnew.mutableY().begin() + 1, histnew.mutableY().end(), histnew.y()[0]);
383 }
384 outputWS->setHistogram(i, histnew);
385 }
386
387 prog.report(reportMsg);
388
390 }
392
393 if (useSparseInstrument) {
394 interpolateFromSparse(*outputWS, *sparseWS, interpolateOpt);
395 }
396
397 return outputWS;
398}
399
401 MatrixWorkspace_uptr outputWS = DataObjects::create<Workspace2D>(inputWS);
402 // The algorithm computes the signal values at bin centres so they should
403 // be treated as a distribution
404 outputWS->setDistribution(true);
405 outputWS->setYUnit("");
406 outputWS->setYUnitLabel("Attenuation factor");
407 return outputWS;
408}
409
412 const auto &spectrumInfo = targetWS.spectrumInfo();
413 const auto refFrame = targetWS.getInstrument()->getReferenceFrame();
414 PARALLEL_FOR_IF(Kernel::threadSafe(targetWS, sparseWS))
415 for (int64_t i = 0; i < static_cast<decltype(i)>(spectrumInfo.size()); ++i) {
417 if (spectrumInfo.hasDetectors(i)) {
418 double lat, lon;
419 std::tie(lat, lon) = spectrumInfo.geographicalAngles(i);
420 const auto spatiallyInterpHisto = sparseWS.bilinearInterpolateFromDetectorGrid(lat, lon);
421 if (spatiallyInterpHisto.size() > 1) {
422 auto targetHisto = targetWS.histogram(i);
423 interpOpt.applyInPlace(spatiallyInterpHisto, targetHisto);
424 targetWS.setHistogram(i, targetHisto);
425 } else {
426 targetWS.mutableY(i) = spatiallyInterpHisto.y().front();
427 }
428 }
430 }
432}
433} // namespace Mantid::Algorithms
const std::string & m_value
Definition: Algorithm.cpp:71
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
double value
The value of the point.
Definition: FitMW.cpp:51
double energy
Definition: GetAllEi.cpp:157
#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.
std::string getPropertyValue(const std::string &name) const override
Get the value of a property as a string.
Definition: Algorithm.cpp:2026
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
static bool isEmpty(const NumT toCheck)
checks that the value was not set by users, uses the value in empty double/int.
This class is shared by a few Workspace types and holds information related to a particular experimen...
const SpectrumInfo & spectrumInfo() const
Return a reference to the SpectrumInfo object.
Geometry::Instrument_const_sptr getInstrument() const
Returns the parameterized instrument.
const Sample & sample() const
Sample accessors.
virtual MantidVec & dataY()=0
virtual MantidVec & dataE()=0
A validator which checks that a workspace has a valid instrument.
Base MatrixWorkspace Abstract Class.
virtual ISpectrum & getSpectrum(const size_t index)=0
Return the underlying ISpectrum ptr at the given workspace index.
HistogramData::Points points(const size_t index) const
virtual std::size_t blocksize() const =0
Returns the size of each block of data returned by the dataY accessors.
virtual std::size_t getNumberHistograms() const =0
Returns the number of histograms in the workspace.
void setHistogram(const size_t index, T &&...data) &
std::size_t yIndexOfX(const double xValue, const std::size_t &index=0, const double tolerance=0.0) const
Returns the y index which corresponds to the X Value provided.
HistogramData::Histogram histogram(const size_t index) const
Returns the Histogram at the given workspace index.
HistogramData::HistogramE & mutableE(const size_t index) &
HistogramData::HistogramY & mutableY(const size_t index) &
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
This class stores information about the sample used in particular run.
Definition: Sample.h:33
A property class for workspaces.
A validator which checks that the unit of the workspace referred to by a WorkspaceProperty is the exp...
static std::unique_ptr< IBeamProfile > createBeamProfile(const Geometry::Instrument &instrument, const API::Sample &sample)
Base class for all classes defining a beam profile.
Definition: IBeamProfile.h:26
Defines a base class for objects describing a volume where interactions of Tracks and Objects can tak...
Class to provide a consistent interface to an interpolation option on algorithms.
std::string validateInputSize(const size_t size) const
Validate the size of input histogram.
void applyInPlace(const HistogramData::Histogram &in, HistogramData::Histogram &out) const
Apply the interpolation method to the output histogram.
void set(const Value &kind, const bool calculateErrors, const bool independentErrors)
Set the interpolation option.
void applyInplace(HistogramData::Histogram &inOut, size_t stepSize) const
Apply the interpolation method to the given histogram.
Stores statistics relating to the tracks generated in MCInteractionVolume for a specific detector.
std::string generateScatterPointStats()
Log a debug string summarising which parts of the environment the simulated scatter points occurred i...
Calculates attenuation due to absorption and scattering in a sample + its environment using a Monte C...
virtual std::shared_ptr< IMCAbsorptionStrategy > createStrategy(IMCInteractionVolume &interactionVol, const IBeamProfile &beamProfile, Kernel::DeltaEMode::Type EMode, const size_t nevents, const size_t maxScatterPtAttempts, const bool regenerateTracksForEachLambda)
Factory method to return an instance of the required absorption strategy class.
void interpolateFromSparse(API::MatrixWorkspace &targetWS, const SparseWorkspace &sparseWS, const Mantid::Algorithms::InterpolationOption &interpOpt)
virtual std::shared_ptr< SparseWorkspace > createSparseWorkspace(const API::MatrixWorkspace &modelWS, const size_t wavelengthPoints, const size_t rows, const size_t columns)
Factory method to return an instance of the required SparseInstrument class.
API::MatrixWorkspace_uptr createOutputWorkspace(const API::MatrixWorkspace &inputWS) const
std::map< std::string, std::string > validateInputs() override
Validate the input properties.
virtual std::unique_ptr< InterpolationOption > createInterpolateOption()
Factory method to return an instance of the required InterpolationOption class.
virtual std::shared_ptr< IMCInteractionVolume > createInteractionVolume(const API::Sample &sample, const size_t maxScatterPtAttempts, const MCInteractionVolume::ScatteringPointVicinity pointsIn)
Factory method to return an instance of the required interaction volume class.
API::MatrixWorkspace_uptr doSimulation(const API::MatrixWorkspace &inputWS, const size_t nevents, const bool simulateTracksForEachWavelength, const int seed, const InterpolationOption &interpolateOpt, const bool useSparseInstrument, const size_t maxScatterPtAttempts, const MCInteractionVolume::ScatteringPointVicinity pointsIn)
Run the simulation over the whole input workspace.
Defines functions and utilities to create and deal with sparse instruments.
virtual HistogramData::Histogram bilinearInterpolateFromDetectorGrid(const double lat, const double lon) const
Spatially interpolate a single histogram from nearby detectors using bilinear interpolation method.
Concrete workspace implementation.
Definition: Workspace2D.h:29
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
bool is(int level) const
Returns true if at least the given log level is set.
Definition: Logger.cpp:146
This implements the the Mersenne Twister 19937 pseudo-random number generator algorithm as a specialz...
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.
std::unique_ptr< MatrixWorkspace > MatrixWorkspace_uptr
unique pointer to Mantid::API::MatrixWorkspace
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
std::shared_ptr< SparseWorkspace > SparseWorkspace_sptr
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
A namespace containing physical constants that are required by algorithms and unit routines.
Definition: Atom.h:14
static constexpr double NeutronMass
Mass of the neutron in kg.
static constexpr double h
Planck constant in J*s.
static constexpr double meV
1 meV in Joules.
constexpr int EMPTY_INT() noexcept
Returns what we consider an "empty" integer within a property.
Definition: EmptyValues.h:25
int32_t detid_t
Typedef for a detector ID.
Definition: SpectrumInfo.h:21
Type
Define the available energy transfer modes It is important to assign enums proper numbers,...
Definition: DeltaEMode.h:29
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54