36using namespace Kernel;
37using namespace CurveFitting;
38using namespace CurveFitting::Functions;
43const size_t MAX_SCATTER_PT_TRIES = 500;
53 :
Algorithm(), m_acrossIdx(0), m_upIdx(1), m_beamIdx(3), m_beamDir(), m_srcR2(0.0), m_halfSampleHeight(0.0),
54 m_halfSampleWidth(0.0), m_halfSampleThick(0.0), m_sampleShape(
nullptr), m_sampleProps(
nullptr), m_detHeight(-1.0),
55 m_detWidth(-1.0), m_detThick(-1.0), m_tmin(-1.0), m_tmax(-1.0), m_delt(-1.0), m_foilRes(-1.0), m_nscatters(0),
56 m_nruns(0), m_nevents(0), m_progress(
nullptr), m_inputWS() {}
61void VesuvioCalculateMS::init() {
63 auto inputWSValidator = std::make_shared<CompositeValidator>();
67 "Input workspace to be corrected, in units of TOF.");
70 auto positiveInt = std::make_shared<Kernel::BoundedValidator<int>>();
71 positiveInt->setLower(1);
72 declareProperty(
"NoOfMasses", -1, positiveInt,
"The number of masses contained within the sample");
74 auto positiveNonZero = std::make_shared<BoundedValidator<double>>();
75 positiveNonZero->setLower(0.0);
76 positiveNonZero->setLowerExclusive(
true);
77 declareProperty(
"SampleDensity", -1.0, positiveNonZero,
"The density of the sample in gm/cm^3");
79 auto nonEmptyArray = std::make_shared<ArrayLengthValidator<double>>();
80 nonEmptyArray->setLengthMin(3);
82 "Atomic properties of masses within the sample. "
83 "The expected format is 3 consecutive values per mass: "
84 "mass(amu), cross-section (barns) & s.d of Compton profile.");
90 declareProperty(
"BeamRadius", -1.0, positiveNonZero,
"Radius, in cm, of beam");
93 declareProperty(
"Seed", 123456789, positiveInt,
"Seed the random number generator with this value");
94 declareProperty(
"NumScatters", 3, positiveInt,
"Number of scattering orders to calculate");
95 declareProperty(
"NumRuns", 10, positiveInt,
"Number of simulated runs per spectrum");
96 declareProperty(
"NumEventsPerRun", 50000, positiveInt,
"Number of events per run");
104 "Workspace to store the calculated total scattering counts");
106 "Workspace to store the calculated multiple scattering counts summed for "
113void VesuvioCalculateMS::exec() {
122 const size_t nhist =
m_inputWS->getNumberHistograms();
124 const auto &spectrumInfo =
m_inputWS->spectrumInfo();
127 for (int64_t i = 0; i < static_cast<int64_t>(nhist); ++i) {
131 totalsc->setSharedX(i,
m_inputWS->sharedX(i));
133 multsc->setSharedX(i,
m_inputWS->sharedX(i));
136 if (!spectrumInfo.hasDetectors(i)) {
137 std::ostringstream os;
138 os <<
"No valid detector object found for spectrum at workspace index '" << i <<
"'. No correction calculated.";
149 calculateMS(rng, i, totalsc->getSpectrum(i), multsc->getSpectrum(i));
161void VesuvioCalculateMS::cacheInputs() {
166 m_nruns =
static_cast<size_t>(nruns);
168 m_nevents =
static_cast<size_t>(nevents);
171 const auto instrument =
m_inputWS->getInstrument();
172 m_beamDir = instrument->getSample()->getPos() - instrument->getSource()->getPos();
175 const auto rframe = instrument->getReferenceFrame();
177 m_upIdx = rframe->pointingUp();
199 m_tmin = inX.front() * 1e-06;
200 m_tmax = inX.back() * 1e-06;
201 m_delt = (inX[1] - inX.front());
205 std::vector<double> sampleInfo =
getProperty(
"AtomicProperties");
206 const auto nInputAtomProps =
static_cast<int>(sampleInfo.size());
207 const int nExptdAtomProp(3);
208 if (nInputAtomProps != nExptdAtomProp * nmasses) {
209 std::ostringstream os;
210 os <<
"Inconsistent AtomicProperties list defined. Expected " << nExptdAtomProp * nmasses
211 <<
" values, however, only " << sampleInfo.size() <<
" have been given.";
212 throw std::invalid_argument(os.str());
214 const int natoms = nInputAtomProps / 3;
215 m_sampleProps = std::make_unique<SampleComptonProperties>(natoms);
218 double totalMass(0.0);
220 for (
int i = 0; i < natoms; ++i) {
222 comptonAtom.mass = sampleInfo[nExptdAtomProp * i];
225 const double xsec = sampleInfo[nExptdAtomProp * i + 1];
226 comptonAtom.sclength = sqrt(xsec / (4.0 * M_PI));
230 comptonAtom.profile = sampleInfo[nExptdAtomProp * i + 2];
232 const double numberDensity =
m_sampleProps->density * 1e6 / totalMass;
236 const auto &spectrumInfo =
m_inputWS->spectrumInfo();
238 for (
size_t i = 0; i <
m_inputWS->getNumberHistograms(); ++i) {
239 if (!spectrumInfo.hasDetectors(i))
241 if (!spectrumInfo.isMonitor(i)) {
248 throw std::runtime_error(
"Failed to get detector");
253 const auto &detectorInfo =
m_inputWS->detectorInfo();
254 const size_t detIndex = detectorInfo.indexOf(spectrumInfo.detector(
index).getID());
255 const auto pixelShape = detectorInfo.detector(detIndex).shape();
257 if (!pixelShape || !pixelShape->hasValidShape()) {
258 throw std::invalid_argument(
"Detector pixel has no defined shape!");
261 V3D detBoxWidth = detBounds.
width();
267 auto foil = instrument->getComponentByName(
"foil-pos0");
269 throw std::runtime_error(
"Workspace has no gold foil component defined.");
271 auto param =
m_inputWS->constInstrumentParameters().get(foil.get(),
"hwhm_lorentz");
273 throw std::runtime_error(
"Foil component has no hwhm_lorentz parameter defined.");
297 for (
size_t i = 0; i <
m_nruns; ++i) {
345 const auto &counts = avgCounts.
sim.
counts[i];
349 const auto &scerrors = avgCounts.
errors[i];
351 std::transform(scerrors.begin(), scerrors.end(), msscatE.begin(), msscatE.begin(),
355 auto &totalscY = totalsc.
mutableY();
356 auto &totalscE = totalsc.
mutableE();
357 const auto &counts0 = avgCounts.
sim.
counts.front();
358 std::transform(counts0.begin(), counts0.end(), msscatY.begin(), totalscY.begin(), std::plus<double>());
359 const auto &errors0 = avgCounts.
errors.front();
360 std::transform(errors0.begin(), errors0.end(), msscatE.begin(), totalscE.begin(),
374 double weightSum(0.0);
388 const double vel2 = sqrt(detpar.
efixed / MASS_TO_MEV);
389 const double t2 = detpar.
l2 / vel2;
397 V3D startPos(srcPos);
400 generateScatter(rng, startPos, neutronDirs[0], weights[0], scatterPts[0]);
401 double distFromStart = startPos.
distance(scatterPts[0]);
403 const double vel0 = sqrt(en1[0] / MASS_TO_MEV);
404 tofs[0] += (distFromStart * 1e6 / vel0);
408 weights[i] = weights[i - 1];
409 tofs[i] = tofs[i - 1];
412 const V3D &prevSc = scatterPts[i - 1];
413 V3D &curSc = scatterPts[i];
414 const V3D &oldDir = neutronDirs[i - 1];
415 V3D &newDir = neutronDirs[i];
418 const double randth = acos(2.0 * rng.
flat() - 1.0);
419 const double randphi = 2.0 * M_PI * rng.
flat();
423 const double wgt = weights[i];
430 }
while (ntries < MAX_SCATTER_PT_TRIES);
431 if (ntries == MAX_SCATTER_PT_TRIES) {
432 throw std::runtime_error(
"Cannot generate valid trajectory from within "
433 "the sample that intersects the sample. Does it "
434 "have a valid shape?");
437 const double scang = newDir.
angle(oldDir);
440 en1[i] = e1range.first + rng.
flat() * (e1range.second - e1range.first);
443 double weight = d2sig * 4.0 * M_PI * (e1range.second - e1range.first) /
m_sampleProps->totalxsec;
446 weights[i] *= weight;
449 const double veli = sqrt(en1[i] / MASS_TO_MEV);
450 tofs[i] += (curSc.
distance(prevSc) * 1e6 / veli);
456 double scang(0.0), distToExit(0.0);
461 double &curWgt = weights[i];
467 const double veli = sqrt(efinal / MASS_TO_MEV);
468 tofs[i] += detpar.
t0 + (scatterPts[i].distance(detPos) * 1e6) / veli;
470 std::vector<double> &counts = simulation.
counts[i];
471 const double finalTOF = tofs[i];
473 for (
size_t it = 0; it < inX.size(); ++it) {
474 if (inX[it] - 0.5 *
m_delt < finalTOF && finalTOF < inX[it] + 0.5 *
m_delt) {
475 counts[it] += curWgt;
492 const double l1)
const {
493 double radius(-1.0), widthPos(0.0), heightPos(0.0);
498 radius = sqrt(widthPos * widthPos + heightPos * heightPos);
518 const double t2,
double &weight)
const {
520 const double t1 = (tof - t2);
521 const double vel0 = l1 / t1;
522 const double en0 = MASS_TO_MEV * vel0 * vel0;
524 weight = 2.0 * en0 / t1 / pow(en0, 0.9);
541 const double dtof,
const double dl1)
const {
542 const double vel1 = sqrt(en0 / MASS_TO_MEV);
543 const double dt1 = (dl1 / vel1) * 1e6;
544 const double xmin(0.0), xmax(15.0 * dt1);
545 double dx = 0.5 * (xmax - xmin);
547 const double yv = rng.
flat();
550 double tof = rng.
gaussian(0.0, dtof);
554 double y = 1.0 - (0.5 * xt * xt / (dt1 * dt1) + xt / dt1 + 1) * exp(-xt / dt1);
555 if (
fabs(
y - yv) < 1e-4) {
560 dx = -
fabs(0.5 * dx);
581 V3D &scatterPt)
const {
582 Track scatterTrack(startPos, direc);
587 const auto &link = scatterTrack.
cbegin();
588 double totalObjectDist = link->distInsideObject;
589 const double scatterProb = 1.0 - exp(-
m_sampleProps->mu * totalObjectDist);
593 const double fraction = dist / totalObjectDist;
595 scatterPt = link->entryPoint;
596 V3D edgeDistances = (link->exitPoint - link->entryPoint);
597 scatterPt += edgeDistances * fraction;
599 weight *= scatterProb;
608std::pair<double, double> VesuvioCalculateMS::calculateE1Range(
const double theta,
const double en0)
const {
610 const double sth(sin(theta)), cth(cos(theta));
612 double e1min(1e10), e1max(-1e10);
614 for (
const auto &atom : atoms) {
615 const double mass = atom.mass;
617 const double fraction = (cth + sqrt(mass * mass - sth * sth)) / (1.0 + mass);
618 const double k1 = fraction * k0;
620 const double qr = sqrt(k0 * k0 + k1 * k1 - 2.0 * k0 * k1 * cth);
621 const double wr = en0 - en1;
623 const double e1a = en0 - wr - 10.0 * width;
624 const double e1b = en0 - wr + 10.0 * width;
632 return std::make_pair(e1min, e1max);
642double VesuvioCalculateMS::partialDiffXSec(
const double en0,
const double en1,
const double theta)
const {
643 const double rt2pi = sqrt(2.0 * M_PI);
647 const double q = sqrt(k0 * k0 + k1 * k1 - 2.0 * k0 * k1 * cos(theta));
648 const double w = en0 - en1;
654 for (
const auto &atom : atoms) {
655 const double jstddev = atom.profile;
656 const double mass = atom.mass;
657 const double y = mass * w / (4.18036 * q) - 0.5 * q;
658 const double jy = exp(-0.5 *
y *
y / (jstddev * jstddev)) / (jstddev * rt2pi);
659 const double sqw = mass * jy / (4.18036 * q);
661 const double sclength = atom.sclength;
662 pdcs += sclength * sclength * (k1 / k0) * sqw;
665 for (
const auto &atom : atoms) {
666 const double sclength = atom.sclength;
667 pdcs += sclength * sclength;
690 const V3D &nominalPos,
const double energy,
const V3D &scatterPt,
691 const V3D &direcBeforeSc,
double &scang,
double &distToExit)
const {
693 const double mu = 7430.0 / sqrt(
energy);
713 scang = direcBeforeSc.
angle(scToDet);
714 const auto &link = scatterToDet.
cbegin();
715 distToExit = link->distInsideObject;
721 }
while (ntries < MAX_SCATTER_PT_TRIES);
722 if (ntries == MAX_SCATTER_PT_TRIES) {
740 const double e1nom,
const double e1res)
const {
744 const double randv = rng.
flat();
745 if (e1nom < 5000.0) {
#define DECLARE_ALGORITHM(classname)
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_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.
A "spectrum" is an object that holds the data for a particular spectrum, in particular:
HistogramData::HistogramY & mutableY() &
HistogramData::HistogramE & mutableE() &
Verify that a workspace has valid sample shape.
A property class for workspaces.
A validator which checks that the unit of the workspace referred to by a WorkspaceProperty is the exp...
static DetectorParams getDetectorParameters(const API::MatrixWorkspace_const_sptr &ws, const size_t index)
Creates a POD struct containing the required detector parameters for this spectrum.
Calculates the multiple scattering & total scattering contributions for a flat-plate or cylindrical s...
void cacheInputs()
Caches inputs insuitable form for speed in later calculations.
std::pair< double, double > calculateE1Range(const double theta, const double en0) const
Kernel::V3D generateDetectorPos(CurveFitting::MSVesuvioHelper::RandomVariateGenerator &rng, const Kernel::V3D &nominalPos, const double energy, const Kernel::V3D &scatterPt, const Kernel::V3D &direcBeforeSc, double &scang, double &distToExit) const
Generate a random position within the final detector in the lab frame.
std::unique_ptr< SampleComptonProperties > m_sampleProps
double generateE0(CurveFitting::MSVesuvioHelper::RandomVariateGenerator &rng, const double l1, const double t2, double &weight) const
Generate an incident energy based on a randomly-selected TOF value It is assigned a weight = (2....
void calculateMS(CurveFitting::MSVesuvioHelper::RandomVariateGenerator &rng, const size_t wsIndex, API::ISpectrum &totalsc, API::ISpectrum &multsc) const
Calculate the total scattering and contributions from higher-order scattering for given spectrum.
std::unique_ptr< API::Progress > m_progress
double generateE1(CurveFitting::MSVesuvioHelper::RandomVariateGenerator &rng, const double angle, const double e1nom, const double e1res) const
Generate the final energy of the analyser.
double calculateCounts(CurveFitting::MSVesuvioHelper::RandomVariateGenerator &rng, const DetectorParams &detpar, const Functions::ResolutionParams &respar, MSVesuvioHelper::Simulation &simulation) const
API::MatrixWorkspace_sptr m_inputWS
Kernel::V3D generateSrcPos(CurveFitting::MSVesuvioHelper::RandomVariateGenerator &rng, const double l1) const
Sample from the moderator assuming it can be seen as a cylindrical ring with inner and outer radius.
Geometry::IObject const * m_sampleShape
bool generateScatter(CurveFitting::MSVesuvioHelper::RandomVariateGenerator &rng, const Kernel::V3D &startPos, const Kernel::V3D &direc, double &weight, Kernel::V3D &scatterPt) const
Generate a scatter event and update the weight according to the amount the beam would be attenuted by...
double partialDiffXSec(const double en0, const double en1, const double theta) const
Compute the partial differential cross section for this energy and theta.
void assignToOutput(const MSVesuvioHelper::SimulationWithErrors &avgCounts, API::ISpectrum &totalsc, API::ISpectrum &multsc) const
Assign the averaged counts to the correct workspaces.
double generateTOF(CurveFitting::MSVesuvioHelper::RandomVariateGenerator &rng, const double en0, const double dtof, const double dl1) const
Generate an initial tof from this distribution: 1-(0.5*X**2/T0**2+X/T0+1)*EXP(-X/T0),...
void simulate(CurveFitting::MSVesuvioHelper::RandomVariateGenerator &rng, const DetectorParams &detpar, const Functions::ResolutionParams &respar, MSVesuvioHelper::Simulation &simulCounts) const
Perform a single simulation of a given number of events for up to a maximum number of scatterings on ...
double m_halfSampleHeight
static ResolutionParams getResolutionParameters(const API::MatrixWorkspace_const_sptr &ws, const size_t index)
Creates a POD struct containing the required resolution parameters for this spectrum.
double flat()
Returns a flat random number between 0.0 & 1.0.
double gaussian(const double mean, const double sigma)
Returns a random number distributed by a normal distribution.
A simple structure that defines an axis-aligned cuboid shaped bounding box for a geometrical object.
Kernel::V3D width() const
Returns the width of the box.
virtual int interceptSurface(Geometry::Track &) const =0
virtual const BoundingBox & getBoundingBox() const =0
Return cached value of axis-aligned bounding box.
Defines a track as a start point and a direction.
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 setPropertyGroup(const std::string &name, const std::string &group)
Set the group for a given property.
void information(const std::string &msg)
Logs at information level.
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
double distance(const V3D &v) const noexcept
Calculates the distance between two vectors.
void azimuth_polar_SNS(const double R, const double azimuth, const double polar) noexcept
Sets the vector position based on azimuth and polar angle, in RADIANS, in the SNS instrument coordina...
double normalize()
Make a normalized vector (return norm value)
double angle(const V3D &) const
Angle between this and another vector.
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
double finalEnergyUranium(const double randv)
Generate the final energy of a neutron for uranium foil analyser at 293K with number density of 1....
double finalEnergyAuYap(const double randv)
Generate the final energy of a neutron for gold foil analyser at 293K with number density of 7....
double finalEnergyAuDD(const double randv)
Generate the final energy of a neutron for gold foil analyser at 293K in double-difference mode:
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.
static constexpr double E_mev_toNeutronWavenumberSq
Transformation coefficient to transform neutron energy into neutron wavevector: K-neutron[m^-10] = sq...
static constexpr double NeutronMassAMU
Mass of the neutron in AMU.
static constexpr double NeutronMass
Mass of the neutron in kg.
static constexpr double AtomicMassUnit
AMU in kg.
static constexpr double meV
1 meV in Joules.
std::string to_string(const wide_integer< Bits, Signed > &n)
Simple data structure to store nominal detector values It avoids some functions taking a huge number ...
double efixed
final energy
double theta
scattering angle in radians
double l1
source-sample distance in metres
double t0
time delay in seconds
Kernel::V3D pos
Full 3D position.
double l2
sample-detector distance in metres
Simple data structure to store resolution parameter values It avoids some functions taking a huge num...
double dtof
spread in tof measurement (us)
double dl1
spread in source-sample distance (m)
Simulation & newSimulation(const size_t order, const size_t ntimes)
SimulationWithErrors average() const
std::vector< std::vector< double > > errors
void normalise()
Normalise the counts so that the integral over the single-scatter events is 1.
std::vector< std::vector< double > > counts
@ Input
An input workspace.
@ Output
An output workspace.
Functor used for computing the sum of the square values of a vector, using the accumulate algorithm.