39constexpr int DEFAULT_NEVENTS = 1000;
40constexpr int DEFAULT_SEED = 123456789;
41constexpr int DEFAULT_LATITUDINAL_DETS = 5;
42constexpr int DEFAULT_LONGITUDINAL_DETS = 10;
45inline double toWavelength(
double energy) {
46 static const double factor =
48 return factor / sqrt(
energy);
51struct EFixedProvider {
52 explicit EFixedProvider(
const ExperimentInfo &expt) : m_expt(expt), m_emode(expt.getEMode()),
m_value(0.0) {
62 return m_expt.getEFixed(detID);
86 auto wsValidator = std::make_shared<CompositeValidator>();
91 "The name of the input workspace. The input workspace must "
92 "have X units of wavelength.");
94 "The name to use for the output workspace.");
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");
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));
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);
156 const bool resimulateTracks =
getProperty(
"ResimulateTracksForDifferentWavelengths");
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") {
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));
179 std::map<std::string, std::string> issues;
180 const bool resimulateTracksForDiffWavelengths =
getProperty(
"ResimulateTracksForDifferentWavelengths");
182 if (resimulateTracksForDiffWavelengths) {
183 const int nlambda =
getProperty(
"NumberOfWavelengthPoints");
186 interpOpt.
set(interpValue,
true, resimulateTracksForDiffWavelengths);
188 if (!nlambdaIssue.empty()) {
189 issues[
"NumberOfWavelengthPoints"] = nlambdaIssue;
204std::shared_ptr<IMCInteractionVolume>
207 auto interactionVol = std::make_shared<MCInteractionVolume>(sample, maxScatterPtAttempts, pointsIn);
208 return interactionVol;
225std::shared_ptr<IMCAbsorptionStrategy>
228 const size_t maxScatterPtAttempts,
const bool regenerateTracksForEachLambda) {
229 auto MCAbs = std::make_shared<MCAbsorptionStrategy>(interactionVol, beamProfile, EMode, nevents, maxScatterPtAttempts,
230 regenerateTracksForEachLambda);
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);
256 auto interpolationOpt = std::make_unique<InterpolationOption>();
257 return interpolationOpt;
275 const bool resimulateTracksForDiffWavelengths,
const int seed,
277 const bool useSparseInstrument,
278 const size_t maxScatterPtAttempts,
281 const auto inputNbins =
static_cast<int>(inputWS.
blocksize());
284 if (resimulateTracksForDiffWavelengths) {
286 if (
isEmpty(nlambda) || nlambda > inputNbins) {
288 g_log.
warning() <<
"The requested number of wavelength points is larger "
289 "than the spectra size. "
290 "Defaulting to spectra size.\n";
292 nlambda = inputNbins;
295 nlambda = inputNbins;
298 if (useSparseInstrument) {
299 const int latitudinalDets =
getProperty(
"NumberOfDetectorRows");
300 const int longitudinalDets =
getProperty(
"NumberOfDetectorColumns");
303 MatrixWorkspace &simulationWS = useSparseInstrument ? *sparseWS : *outputWS;
304 const MatrixWorkspace &instrumentWS = useSparseInstrument ? simulationWS : inputWS;
309 EFixedProvider
efixed(instrumentWS);
313 Progress prog(
this, 0.0, 1.0, nhists);
315 const std::string reportMsg =
"Computing corrections";
319 auto strategy =
createStrategy(*interactionVolume, *beamProfile,
efixed.emode(), nevents, maxScatterPtAttempts,
320 resimulateTracksForDiffWavelengths);
325 for (int64_t i = 0; i < nhists; ++i) {
328 auto &outE = simulationWS.
mutableE(i);
332 if (!spectrumInfo.hasDetectors(i) || spectrumInfo.isMasked(i)) {
336 const auto &detPos = spectrumInfo.position(i);
337 const double lambdaFixed = toWavelength(
efixed.value(spectrumInfo.detector(i).getID()));
340 const auto lambdas = simulationWS.
points(i).rawData();
342 const auto nbins = lambdas.size();
343 const size_t lambdaStepSize = nbins / nlambda;
345 std::vector<double> packedLambdas;
346 std::vector<double> packedAttFactors;
347 std::vector<double> packedAttFactorErrors;
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);
354 if (lambdaStepSize > 1 && j + lambdaStepSize >= nbins && j + 1 != nbins) {
355 j = nbins - lambdaStepSize - 1;
360 strategy->calculate(rng, detPos, packedLambdas, lambdaFixed, packedAttFactors, packedAttFactorErrors,
363 if (
g_log.
is(Kernel::Logger::Priority::PRIO_DEBUG)) {
367 for (
size_t j = 0; j < packedLambdas.size(); j++) {
368 auto idx = simulationWS.
yIndexOfX(packedLambdas[j], i);
376 if (!useSparseInstrument && lambdaStepSize > 1) {
377 auto histnew = simulationWS.
histogram(i);
379 if (lambdaStepSize < nbins) {
382 std::fill(histnew.mutableY().begin() + 1, histnew.mutableY().end(), histnew.y()[0]);
384 outputWS->setHistogram(i, histnew);
393 if (useSparseInstrument) {
404 outputWS->setDistribution(
true);
405 outputWS->setYUnit(
"");
406 outputWS->setYUnitLabel(
"Attenuation factor");
413 const auto refFrame = targetWS.
getInstrument()->getReferenceFrame();
415 for (int64_t i = 0; i < static_cast<decltype(i)>(spectrumInfo.size()); ++i) {
417 if (spectrumInfo.hasDetectors(i)) {
419 std::tie(lat, lon) = spectrumInfo.geographicalAngles(i);
421 if (spatiallyInterpHisto.size() > 1) {
422 auto targetHisto = targetWS.
histogram(i);
423 interpOpt.
applyInPlace(spatiallyInterpHisto, targetHisto);
426 targetWS.
mutableY(i) = spatiallyInterpHisto.y().front();
const std::string & m_value
#define DECLARE_ALGORITHM(classname)
double value
The value of the point.
#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_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.
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
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.
This class stores information about the sample used in particular run.
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.
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 exec() override
Execution code.
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.
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.
bool is(int level) const
Returns true if at least the given log level is set.
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.
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.
A namespace containing physical constants that are required by algorithms and unit routines.
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.
int32_t detid_t
Typedef for a detector ID.
Type
Define the available energy transfer modes It is important to assign enums proper numbers,...
@ Input
An input workspace.
@ Output
An output workspace.