Mantid
Loading...
Searching...
No Matches
VesuvioCalculateMS.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 +
8// Use helpers for storing detector/resolution parameters
11
12#include "MantidAPI/Axis.h"
13#include "MantidAPI/Sample.h"
19
24
31
32#include <memory>
33
35using namespace API;
36using namespace Kernel;
37using namespace CurveFitting;
38using namespace CurveFitting::Functions;
39using Geometry::Track;
40
41namespace {
43const size_t MAX_SCATTER_PT_TRIES = 500;
45const double MASS_TO_MEV = 0.5 * PhysicalConstants::NeutronMass / PhysicalConstants::meV;
46} // end anonymous namespace
47
48// Register the algorithm into the AlgorithmFactory
49DECLARE_ALGORITHM(VesuvioCalculateMS)
50
51
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() {}
57
61void VesuvioCalculateMS::init() {
62 // Inputs
63 auto inputWSValidator = std::make_shared<CompositeValidator>();
64 inputWSValidator->add<WorkspaceUnitValidator>("TOF");
65 inputWSValidator->add<SampleShapeValidator>();
66 declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input, inputWSValidator),
67 "Input workspace to be corrected, in units of TOF.");
68
69 // -- Sample --
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");
73
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");
78
79 auto nonEmptyArray = std::make_shared<ArrayLengthValidator<double>>();
80 nonEmptyArray->setLengthMin(3);
81 declareProperty(std::make_unique<ArrayProperty<double>>("AtomicProperties", nonEmptyArray),
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.");
85 setPropertyGroup("NoOfMasses", "Sample");
86 setPropertyGroup("SampleDensity", "Sample");
87 setPropertyGroup("AtomicProperties", "Sample");
88
89 // -- Beam --
90 declareProperty("BeamRadius", -1.0, positiveNonZero, "Radius, in cm, of beam");
91
92 // -- Algorithm --
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");
97 setPropertyGroup("Seed", "Algorithm");
98 setPropertyGroup("NumScatters", "Algorithm");
99 setPropertyGroup("NumRuns", "Algorithm");
100 setPropertyGroup("NumEventsPerRun", "Algorithm");
101
102 // Outputs
103 declareProperty(std::make_unique<WorkspaceProperty<>>("TotalScatteringWS", "", Direction::Output),
104 "Workspace to store the calculated total scattering counts");
105 declareProperty(std::make_unique<WorkspaceProperty<>>("MultipleScatteringWS", "", Direction::Output),
106 "Workspace to store the calculated multiple scattering counts summed for "
107 "all orders");
108}
109
113void VesuvioCalculateMS::exec() {
114 m_inputWS = getProperty("InputWorkspace");
115 cacheInputs();
116
117 // Create new workspaces
120
121 // Setup progress
122 const size_t nhist = m_inputWS->getNumberHistograms();
123 m_progress = std::make_unique<Progress>(this, 0.0, 1.0, nhist * m_nruns * 2);
124 const auto &spectrumInfo = m_inputWS->spectrumInfo();
125
127 for (int64_t i = 0; i < static_cast<int64_t>(nhist); ++i) {
129
130 // set common X-values
131 totalsc->setSharedX(i, m_inputWS->sharedX(i));
132
133 multsc->setSharedX(i, m_inputWS->sharedX(i));
134
135 // Final detector position
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.";
139 g_log.information(os.str());
140 continue;
141 }
142
143 // Initialize random number generator
144 const int seed = getProperty("Seed");
146
147 // the output spectrum objects have references to where the data will be
148 // stored
149 calculateMS(rng, i, totalsc->getSpectrum(i), multsc->getSpectrum(i));
151 }
153
154 setProperty("TotalScatteringWS", totalsc);
155 setProperty("MultipleScatteringWS", multsc);
156}
157
161void VesuvioCalculateMS::cacheInputs() {
162 // Algorithm
163 int nscatters = getProperty("NumScatters");
164 m_nscatters = static_cast<size_t>(nscatters);
165 int nruns = getProperty("NumRuns");
166 m_nruns = static_cast<size_t>(nruns);
167 int nevents = getProperty("NumEventsPerRun");
168 m_nevents = static_cast<size_t>(nevents);
169
170 // -- Geometry --
171 const auto instrument = m_inputWS->getInstrument();
172 m_beamDir = instrument->getSample()->getPos() - instrument->getSource()->getPos();
174
175 const auto rframe = instrument->getReferenceFrame();
176 m_acrossIdx = rframe->pointingHorizontal();
177 m_upIdx = rframe->pointingUp();
178 m_beamIdx = rframe->pointingAlongBeam();
179
180 m_srcR2 = getProperty("BeamRadius");
181 // Convert to metres
182 m_srcR2 /= 100.0;
183
184 // Sample shape
185 m_sampleShape = &(m_inputWS->sample().getShape());
186 // We know the shape is valid from the property validator
187 // Use the bounding box as an approximation to determine the extents
188 // as this will match in both height and width for a cuboid & cylinder
189 // sample shape
191 V3D boxWidth = bounds.width();
192 // Use half-width/height for easier calculation later
193 m_halfSampleWidth = 0.5 * boxWidth[m_acrossIdx];
194 m_halfSampleHeight = 0.5 * boxWidth[m_upIdx];
195 m_halfSampleThick = 0.5 * boxWidth[m_beamIdx];
196
197 // -- Workspace --
198 const auto &inX = m_inputWS->x(0);
199 m_tmin = inX.front() * 1e-06;
200 m_tmax = inX.back() * 1e-06;
201 m_delt = (inX[1] - inX.front());
202
203 // -- Sample --
204 int nmasses = getProperty("NoOfMasses");
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());
213 }
214 const int natoms = nInputAtomProps / 3;
215 m_sampleProps = std::make_unique<SampleComptonProperties>(natoms);
216 m_sampleProps->density = getProperty("SampleDensity");
217
218 double totalMass(0.0); // total mass in grams
219 m_sampleProps->totalxsec = 0.0;
220 for (int i = 0; i < natoms; ++i) {
221 auto &comptonAtom = m_sampleProps->atoms[i];
222 comptonAtom.mass = sampleInfo[nExptdAtomProp * i];
223 totalMass += comptonAtom.mass * PhysicalConstants::AtomicMassUnit * 1000;
224
225 const double xsec = sampleInfo[nExptdAtomProp * i + 1];
226 comptonAtom.sclength = sqrt(xsec / (4.0 * M_PI));
227 const double factor = 1.0 + (PhysicalConstants::NeutronMassAMU / comptonAtom.mass);
228 m_sampleProps->totalxsec += (xsec / (factor * factor));
229
230 comptonAtom.profile = sampleInfo[nExptdAtomProp * i + 2];
231 }
232 const double numberDensity = m_sampleProps->density * 1e6 / totalMass; // formula units/m^3
233 m_sampleProps->mu = numberDensity * m_sampleProps->totalxsec * 1e-28;
234
235 // -- Detector geometry -- choose first detector that is not a monitor
236 const auto &spectrumInfo = m_inputWS->spectrumInfo();
237 int64_t index = -1;
238 for (size_t i = 0; i < m_inputWS->getNumberHistograms(); ++i) {
239 if (!spectrumInfo.hasDetectors(i))
240 continue;
241 if (!spectrumInfo.isMonitor(i)) {
242 index = i;
243 break;
244 }
245 }
246 // Bounding box in detector frame
247 if (index < 0) {
248 throw std::runtime_error("Failed to get detector");
249 }
250 // If is a detector group then take shape of first pixel
251 // All detectors in same bansk should be same shape anyway
252 // If the detector is a DetectorGroup, getID gives ID of first 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();
256
257 if (!pixelShape || !pixelShape->hasValidShape()) {
258 throw std::invalid_argument("Detector pixel has no defined shape!");
259 }
260 Geometry::BoundingBox detBounds = pixelShape->getBoundingBox();
261 V3D detBoxWidth = detBounds.width();
262 m_detWidth = detBoxWidth[m_acrossIdx];
263 m_detHeight = detBoxWidth[m_upIdx];
264 m_detThick = detBoxWidth[m_beamIdx];
265
266 // Foil resolution
267 auto foil = instrument->getComponentByName("foil-pos0");
268 if (!foil) {
269 throw std::runtime_error("Workspace has no gold foil component defined.");
270 }
271 auto param = m_inputWS->constInstrumentParameters().get(foil.get(), "hwhm_lorentz");
272 if (!param) {
273 throw std::runtime_error("Foil component has no hwhm_lorentz parameter defined.");
274 }
275 m_foilRes = param->value<double>();
276}
277
288void VesuvioCalculateMS::calculateMS(CurveFitting::MSVesuvioHelper::RandomVariateGenerator &rng, const size_t wsIndex,
289 API::ISpectrum &totalsc, API::ISpectrum &multsc) const {
290 // Detector information
292 detpar.t0 *= 1e6; // t0 in microseconds here
294
295 // Final counts averaged over all simulations
297 for (size_t i = 0; i < m_nruns; ++i) {
298 m_progress->report("MS calculation: idx=" + std::to_string(wsIndex) + ", run=" + std::to_string(i));
299
300 simulate(rng, detpar, respar, accumulator.newSimulation(m_nscatters, m_inputWS->blocksize()));
301
302 m_progress->report("MS calculation: idx=" + std::to_string(wsIndex) + ", run=" + std::to_string(i));
303 }
304
305 // Average over all runs and assign to output workspaces
307 avgCounts.normalise();
308
309 assignToOutput(avgCounts, totalsc, multsc);
310}
311
322void VesuvioCalculateMS::simulate(CurveFitting::MSVesuvioHelper::RandomVariateGenerator &rng,
323 const DetectorParams &detpar, const ResolutionParams &respar,
325 for (size_t i = 0; i < m_nevents; ++i) {
326 calculateCounts(rng, detpar, respar, simulCounts);
327 }
328}
329
338void VesuvioCalculateMS::assignToOutput(const CurveFitting::MSVesuvioHelper::SimulationWithErrors &avgCounts,
339 API::ISpectrum &totalsc, API::ISpectrum &multsc) const {
340 // Sum up all multiple scatter events
341 auto &msscatY = multsc.mutableY();
342 auto &msscatE = multsc.mutableE();
343 for (size_t i = 1; i < m_nscatters; ++i) //(i >= 1 for multiple scatters)
344 {
345 const auto &counts = avgCounts.sim.counts[i];
346
347 msscatY += counts;
348
349 const auto &scerrors = avgCounts.errors[i];
350 // sum errors in quadrature
351 std::transform(scerrors.begin(), scerrors.end(), msscatE.begin(), msscatE.begin(),
353 }
354 // for total scattering add on single-scatter events
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(),
362}
363
371double VesuvioCalculateMS::calculateCounts(CurveFitting::MSVesuvioHelper::RandomVariateGenerator &rng,
372 const DetectorParams &detpar, const ResolutionParams &respar,
374 double weightSum(0.0);
375
376 // moderator coord in lab frame
377 V3D srcPos = generateSrcPos(rng, detpar.l1);
378 if (fabs(srcPos[m_acrossIdx]) > m_halfSampleWidth || fabs(srcPos[m_upIdx]) > m_halfSampleHeight) {
379 return 0.0; // misses sample
380 }
381
382 // track various variables during calculation
383 std::vector<double> weights(m_nscatters, 1.0), // start at 1.0
384 tofs(m_nscatters,
385 0.0), // tof accumulates for each piece of the calculation
386 en1(m_nscatters, 0.0);
387
388 const double vel2 = sqrt(detpar.efixed / MASS_TO_MEV);
389 const double t2 = detpar.l2 / vel2;
390 en1[0] = generateE0(rng, detpar.l1, t2, weights[0]);
391 tofs[0] = generateTOF(rng, en1[0], respar.dtof,
392 respar.dl1); // correction for resolution in l1
393
394 // Neutron path
395 std::vector<V3D> scatterPts(m_nscatters), // track origin of each scatter
396 neutronDirs(m_nscatters); // neutron directions
397 V3D startPos(srcPos);
398 neutronDirs[0] = m_beamDir;
399
400 generateScatter(rng, startPos, neutronDirs[0], weights[0], scatterPts[0]);
401 double distFromStart = startPos.distance(scatterPts[0]);
402 // Compute TOF for first scatter event
403 const double vel0 = sqrt(en1[0] / MASS_TO_MEV);
404 tofs[0] += (distFromStart * 1e6 / vel0);
405
406 // multiple scatter events within sample, i.e not including zeroth
407 for (size_t i = 1; i < m_nscatters; ++i) {
408 weights[i] = weights[i - 1];
409 tofs[i] = tofs[i - 1];
410
411 // Generate a new direction of travel
412 const V3D &prevSc = scatterPts[i - 1];
413 V3D &curSc = scatterPts[i];
414 const V3D &oldDir = neutronDirs[i - 1];
415 V3D &newDir = neutronDirs[i];
416 size_t ntries(0);
417 do {
418 const double randth = acos(2.0 * rng.flat() - 1.0);
419 const double randphi = 2.0 * M_PI * rng.flat();
420 newDir.azimuth_polar_SNS(1.0, randphi, randth);
421
422 // Update weight
423 const double wgt = weights[i];
424 if (generateScatter(rng, prevSc, newDir, weights[i], curSc))
425 break;
426 else {
427 weights[i] = wgt; // put it back to what it was
428 ++ntries;
429 }
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?");
435 }
436
437 const double scang = newDir.angle(oldDir);
438 auto e1range = calculateE1Range(scang, en1[i - 1]);
439
440 en1[i] = e1range.first + rng.flat() * (e1range.second - e1range.first);
441
442 const double d2sig = partialDiffXSec(en1[i - 1], en1[i], scang);
443 double weight = d2sig * 4.0 * M_PI * (e1range.second - e1range.first) / m_sampleProps->totalxsec;
444 // accumulate total weight
445 weightSum += weight;
446 weights[i] *= weight; // account for this scatter on top of previous
447
448 // Increment time of flight...
449 const double veli = sqrt(en1[i] / MASS_TO_MEV);
450 tofs[i] += (curSc.distance(prevSc) * 1e6 / veli);
451 }
452
453 // force all orders in to current detector
454 const auto &inX = m_inputWS->x(0);
455 for (size_t i = 0; i < m_nscatters; ++i) {
456 double scang(0.0), distToExit(0.0);
457
458 V3D detPos = generateDetectorPos(rng, detpar.pos, en1[i], scatterPts[i], neutronDirs[i], scang, distToExit);
459
460 // Weight by probability neutron leaves sample
461 double &curWgt = weights[i];
462 curWgt *= exp(-m_sampleProps->mu * distToExit);
463 // Weight by cross-section for the final energy
464 const double efinal = generateE1(rng, detpar.theta, detpar.efixed, m_foilRes);
465 curWgt *= partialDiffXSec(en1[i], efinal, scang) / m_sampleProps->totalxsec;
466 // final TOF
467 const double veli = sqrt(efinal / MASS_TO_MEV);
468 tofs[i] += detpar.t0 + (scatterPts[i].distance(detPos) * 1e6) / veli;
469 // "Bin" weight into appropriate place
470 std::vector<double> &counts = simulation.counts[i];
471 const double finalTOF = tofs[i];
472
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;
476 break;
477 }
478 }
479 }
480
481 return weightSum;
482}
483
491V3D VesuvioCalculateMS::generateSrcPos(CurveFitting::MSVesuvioHelper::RandomVariateGenerator &rng,
492 const double l1) const {
493 double radius(-1.0), widthPos(0.0), heightPos(0.0);
494 do {
495 widthPos = -m_srcR2 + 2.0 * m_srcR2 * rng.flat();
496 heightPos = -m_srcR2 + 2.0 * m_srcR2 * rng.flat();
497 using std::sqrt;
498 radius = sqrt(widthPos * widthPos + heightPos * heightPos);
499 } while (radius > m_srcR2);
500 // assign to output
501 V3D srcPos;
502 srcPos[m_acrossIdx] = widthPos;
503 srcPos[m_upIdx] = heightPos;
504 srcPos[m_beamIdx] = -l1;
505 return srcPos;
506}
507
517double VesuvioCalculateMS::generateE0(CurveFitting::MSVesuvioHelper::RandomVariateGenerator &rng, const double l1,
518 const double t2, double &weight) const {
519 const double tof = m_tmin + (m_tmax - m_tmin) * rng.flat();
520 const double t1 = (tof - t2);
521 const double vel0 = l1 / t1;
522 const double en0 = MASS_TO_MEV * vel0 * vel0;
523
524 weight = 2.0 * en0 / t1 / pow(en0, 0.9);
525 weight *= 1e-4; // Reduce weight to ~1
526
527 return en0;
528}
529
540double VesuvioCalculateMS::generateTOF(CurveFitting::MSVesuvioHelper::RandomVariateGenerator &rng, const double en0,
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);
546 // Generate a random y position in th distribution
547 const double yv = rng.flat();
548
549 double xt(xmin);
550 double tof = rng.gaussian(0.0, dtof);
551 while (true) {
552 xt += dx;
553 // Y=1-(0.5*X**2/T0**2+X/T0+1)*EXP(-X/T0)
554 double y = 1.0 - (0.5 * xt * xt / (dt1 * dt1) + xt / dt1 + 1) * exp(-xt / dt1);
555 if (fabs(y - yv) < 1e-4) {
556 tof += xt - 3 * dt1;
557 break;
558 }
559 if (y > yv) {
560 dx = -fabs(0.5 * dx);
561 } else {
562 dx = fabs(0.5 * dx);
563 }
564 }
565 return tof;
566}
567
579bool VesuvioCalculateMS::generateScatter(CurveFitting::MSVesuvioHelper::RandomVariateGenerator &rng,
580 const Kernel::V3D &startPos, const Kernel::V3D &direc, double &weight,
581 V3D &scatterPt) const {
582 Track scatterTrack(startPos, direc);
583 if (m_sampleShape->interceptSurface(scatterTrack) != 1) {
584 return false;
585 }
586 // Find distance inside object and compute probability of scattering
587 const auto &link = scatterTrack.cbegin();
588 double totalObjectDist = link->distInsideObject;
589 const double scatterProb = 1.0 - exp(-m_sampleProps->mu * totalObjectDist);
590 // Select a random point on the track that is the actual scatter point
591 // from the scattering probability distribution
592 const double dist = -log(1.0 - rng.flat() * scatterProb) / m_sampleProps->mu;
593 const double fraction = dist / totalObjectDist;
594 // Scatter point is then entry point + fraction of width in each direction
595 scatterPt = link->entryPoint;
596 V3D edgeDistances = (link->exitPoint - link->entryPoint);
597 scatterPt += edgeDistances * fraction;
598 // Update weight
599 weight *= scatterProb;
600 return true;
601}
602
608std::pair<double, double> VesuvioCalculateMS::calculateE1Range(const double theta, const double en0) const {
609 const double k0 = sqrt(en0 / PhysicalConstants::E_mev_toNeutronWavenumberSq);
610 const double sth(sin(theta)), cth(cos(theta));
611
612 double e1min(1e10), e1max(-1e10); // large so that anything else is smaller
613 const auto &atoms = m_sampleProps->atoms;
614 for (const auto &atom : atoms) {
615 const double mass = atom.mass;
616
617 const double fraction = (cth + sqrt(mass * mass - sth * sth)) / (1.0 + mass);
618 const double k1 = fraction * k0;
619 const double en1 = PhysicalConstants::E_mev_toNeutronWavenumberSq * k1 * k1;
620 const double qr = sqrt(k0 * k0 + k1 * k1 - 2.0 * k0 * k1 * cth);
621 const double wr = en0 - en1;
622 const double width = PhysicalConstants::E_mev_toNeutronWavenumberSq * atom.profile * qr / mass;
623 const double e1a = en0 - wr - 10.0 * width;
624 const double e1b = en0 - wr + 10.0 * width;
625 if (e1a < e1min)
626 e1min = e1a;
627 if (e1b > e1max)
628 e1max = e1b;
629 }
630 if (e1min < 0.0)
631 e1min = 0.0;
632 return std::make_pair(e1min, e1max);
633}
634
642double VesuvioCalculateMS::partialDiffXSec(const double en0, const double en1, const double theta) const {
643 const double rt2pi = sqrt(2.0 * M_PI);
644
645 const double k0 = sqrt(en0 / PhysicalConstants::E_mev_toNeutronWavenumberSq);
646 const double k1 = sqrt(en1 / PhysicalConstants::E_mev_toNeutronWavenumberSq);
647 const double q = sqrt(k0 * k0 + k1 * k1 - 2.0 * k0 * k1 * cos(theta));
648 const double w = en0 - en1;
649
650 double pdcs(0.0);
651 const auto &atoms = m_sampleProps->atoms;
652 if (q > 0.0) // avoid continuous checking in loop
653 {
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);
660
661 const double sclength = atom.sclength;
662 pdcs += sclength * sclength * (k1 / k0) * sqw;
663 }
664 } else {
665 for (const auto &atom : atoms) {
666 const double sclength = atom.sclength;
667 pdcs += sclength * sclength;
668 }
669 }
670
671 return pdcs;
672}
673
689V3D VesuvioCalculateMS::generateDetectorPos(CurveFitting::MSVesuvioHelper::RandomVariateGenerator &rng,
690 const V3D &nominalPos, const double energy, const V3D &scatterPt,
691 const V3D &direcBeforeSc, double &scang, double &distToExit) const {
692 // Inverse attenuation length (m-1) for vesuvio det.
693 const double mu = 7430.0 / sqrt(energy);
694 // Probability of detection in path thickness.
695 const double ps = 1.0 - exp(-mu * m_detThick);
696 V3D detPos;
697 scang = 0.0;
698 distToExit = 0.0;
699 size_t ntries(0);
700 do {
701 // Beam direction by moving to front of "box"define by detector dimensions
702 // and then
703 // computing expected distance travelled based on probability
704 detPos[m_beamIdx] = (nominalPos[m_beamIdx] - 0.5 * m_detThick) - (log(1.0 - rng.flat() * ps) / mu);
705 // perturb away from nominal position
706 detPos[m_acrossIdx] = nominalPos[m_acrossIdx] + (rng.flat() - 0.5) * m_detWidth;
707 detPos[m_upIdx] = nominalPos[m_upIdx] + (rng.flat() - 0.5) * m_detHeight;
708
709 // Distance to exit the sample for this order
710 const V3D scToDet = normalize(detPos - scatterPt);
711 Geometry::Track scatterToDet(scatterPt, scToDet);
712 if (m_sampleShape->interceptSurface(scatterToDet) > 0) {
713 scang = direcBeforeSc.angle(scToDet);
714 const auto &link = scatterToDet.cbegin();
715 distToExit = link->distInsideObject;
716 break;
717 }
718 // if point is very close surface then there may be no valid intercept so
719 // try again
720 ++ntries;
721 } while (ntries < MAX_SCATTER_PT_TRIES);
722 if (ntries == MAX_SCATTER_PT_TRIES) {
723 // Assume it is very close to the surface so that the distance travelled
724 // would
725 // be a neglible contribution
726 distToExit = 0.0;
727 }
728 return detPos;
729}
730
739double VesuvioCalculateMS::generateE1(CurveFitting::MSVesuvioHelper::RandomVariateGenerator &rng, const double angle,
740 const double e1nom, const double e1res) const {
741 if (e1res == 0.0)
742 return e1nom;
743
744 const double randv = rng.flat();
745 if (e1nom < 5000.0) {
746 if (angle > 90.0)
748 else
750 } else {
752 }
753}
754
755} // namespace Mantid::CurveFitting::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
double energy
Definition: GetAllEi.cpp:157
std::map< DeltaEMode::Type, std::string > index
Definition: DeltaEMode.cpp:19
#define fabs(x)
Definition: Matrix.cpp:22
#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
Kernel::Logger & g_log
Definition: Algorithm.h:451
A "spectrum" is an object that holds the data for a particular spectrum, in particular:
Definition: ISpectrum.h:39
HistogramData::HistogramY & mutableY() &
Definition: ISpectrum.h:178
HistogramData::HistogramE & mutableE() &
Definition: ISpectrum.h:182
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.
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
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.
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 ...
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.
Definition: BoundingBox.h:34
Kernel::V3D width() const
Returns the width of the box.
Definition: BoundingBox.h:98
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.
Definition: Track.h:165
LType::const_iterator cbegin() const
Returns an interator to the start of the set of links (const version)
Definition: Track.h:206
Support for a property that holds an array of values.
Definition: ArrayProperty.h:28
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.
Definition: Logger.cpp:105
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
double distance(const V3D &v) const noexcept
Calculates the distance between two vectors.
Definition: V3D.h:287
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...
Definition: V3D.cpp:95
double normalize()
Make a normalized vector (return norm value)
Definition: V3D.cpp:130
double angle(const V3D &) const
Angle between this and another vector.
Definition: V3D.cpp:165
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.
Definition: MultiThreaded.h:22
MANTID_KERNEL_DLL V3D normalize(V3D v)
Normalizes a V3D.
Definition: V3D.h:341
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 theta
scattering angle in radians
double l1
source-sample distance in metres
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)
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.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54
Functor used for computing the sum of the square values of a vector, using the accumulate algorithm.
Definition: VectorHelper.h:150