Mantid
Loading...
Searching...
No Matches
MultipleScatteringCorrection.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2021 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 +
7
12#include "MantidAPI/Sample.h"
22#include "MantidHistogramData/Interpolate.h"
26
27namespace Mantid::Algorithms {
28
29using namespace API;
30using namespace Geometry;
31using namespace Kernel;
32using namespace Mantid::DataObjects;
33
34using HistogramData::interpolateLinearInplace;
35
36DECLARE_ALGORITHM(MultipleScatteringCorrection)
37
38namespace {
39// the maximum number of elements to combine at once in the pairwise summation
40constexpr int64_t MAX_INTEGRATION_LENGTH{1000};
41
42static constexpr double RAD2DEG = 180.0 / M_PI; // save some flops??
43
44inline int64_t findMiddle(const int64_t start, const int64_t stop) {
45 auto half = static_cast<int64_t>(floor(.5 * (static_cast<double>(stop - start))));
46 return start + half;
47}
48
49inline size_t calcLinearIdxFromUpperTriangular(const size_t N, const size_t row_idx, const size_t col_idx) {
50 // calculate the linear index from the upper triangular matrix from a (N x N) matrix
51 // row_idx < col_idx due to upper triangular matrix
52 assert(row_idx < col_idx); // only relevant during Debug build
53 return N * (N - 1) / 2 - (N - row_idx) * (N - row_idx - 1) / 2 + col_idx - row_idx - 1;
54}
55
56// being added to make it clearer that we are creating a unit vector
57inline const V3D getDirection(const V3D &posInitial, const V3D &posFinal) { return normalize(posFinal - posInitial); }
58
59// make code slightly clearer
60inline double getDistanceInsideObject(const IObject &shape, Track &track) {
61 if (shape.interceptSurface(track) > 0) {
62 return track.totalDistInsideObject();
63 } else {
64 return 0.0;
65 }
66}
67
68inline double checkzero(const double x) { return std::abs(x) < std::numeric_limits<float>::min() ? 0.0 : x; }
69} // namespace
70
76 // 1- The input workspace must have an instrument and units of wavelength
77 auto wsValidator = std::make_shared<CompositeValidator>();
78 wsValidator->add<WorkspaceUnitValidator>("Wavelength");
79 wsValidator->add<InstrumentValidator>();
80 // 2- The input workspace must have a sample defined (shape and material)
81 wsValidator->add<SampleValidator, unsigned int>((SampleValidator::Shape | SampleValidator::Material));
82
83 declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input, wsValidator),
84 "The X values for the input workspace must be in units of wavelength");
85
86 auto positiveInt = std::make_shared<BoundedValidator<int64_t>>();
87 positiveInt->setLower(1);
88 declareProperty("NumberOfWavelengthPoints", static_cast<int64_t>(EMPTY_INT()), positiveInt,
89 "The number of wavelength points for which the numerical integral is\n"
90 "calculated (default: all points)");
91
92 auto moreThanZero = std::make_shared<BoundedValidator<double>>();
93 moreThanZero->setLower(0.001);
94 declareProperty("ElementSize", 1.0, moreThanZero, "The size of one side of an integration element cube in mm");
95
96 declareProperty("ContainerElementSize", EMPTY_DBL(),
97 "The size of one side of an integration element cube in mm for container."
98 "Default to be the same as ElementSize.");
99
100 std::vector<std::string> methodOptions{"SampleOnly", "SampleAndContainer"};
101 declareProperty("Method", "SampleOnly", std::make_shared<StringListValidator>(methodOptions),
102 "Correction method, use either SampleOnly or SampleAndContainer.");
103
104 declareProperty(std::make_unique<WorkspaceProperty<WorkspaceGroup>>("OutputWorkspace", "", Direction::Output),
105 "Output workspace name. "
106 "A Workspace2D containing the correction matrix that can be directly applied to the corresponding "
107 "Event workspace for multiple scattering correction.");
108}
109
115std::map<std::string, std::string> MultipleScatteringCorrection::validateInputs() {
116 std::map<std::string, std::string> result;
117
118 // check 0: input workspace must have a valida sample
119 // NOTE: technically the workspace validator should be able to catch the undefined sample
120 // error. Keeping this check here in case the validator is changed in the future.
121 m_inputWS = getProperty("InputWorkspace");
122 const auto &sample = m_inputWS->sample();
123 if (!sample.getShape().hasValidShape()) {
124 result["InputWorkspace"] = "The input workspace must have a valid sample shape";
125 }
126
127 // check 1: input workspace must have a valid sample environment
128 std::string method = getProperty("Method");
129 if (method == "SampleAndContainer") {
130 const auto &containerShape = m_inputWS->sample().getEnvironment().getContainer();
131 if (!containerShape.hasValidShape()) {
132 result["Method"] = "SampleAndContainer requires a valid container shape.";
133 }
134 }
135
136 return result;
137}
138
144 // Parse input properties and assign corresponding values to the member
145 // variables
146 parseInputs();
147
148 std::string method = getProperty("Method");
149 if (method == "SampleOnly") {
150 //-- Setup output workspace
151 // set the OutputWorkspace to a workspace group with one workspace:
152 // ${OutputWorkspace}_sampleOnly
153 API::MatrixWorkspace_sptr ws_sampleOnly = create<HistoWorkspace>(*m_inputWS);
154 ws_sampleOnly->setYUnit(""); // Need to explicitly set YUnit to nothing
155 ws_sampleOnly->setDistribution(true); // The output of this is a distribution
156 ws_sampleOnly->setYUnitLabel("Multiple Scattering Correction factor");
157 //-- Fill the workspace with sample only correction factors
158 const auto &sampleShape = m_inputWS->sample().getShape();
159 calculateSingleComponent(ws_sampleOnly, sampleShape, m_sampleElementSize);
160 //-- Package output to workspace group
161 const std::string outWSName = getProperty("OutputWorkspace");
162 std::vector<std::string> names;
163 names.emplace_back(outWSName + "_sampleOnly");
164 API::AnalysisDataService::Instance().addOrReplace(names.back(), ws_sampleOnly);
165 // group
166 auto group = createChildAlgorithm("GroupWorkspaces");
167 group->initialize();
168 group->setProperty("InputWorkspaces", names);
169 group->setProperty("OutputWorkspace", outWSName);
170 group->execute();
171 API::WorkspaceGroup_sptr outWS = group->getProperty("OutputWorkspace");
172 // NOTE:
173 // The output here is a workspace group of one, and it is an intended design as
174 // the MantidTotalScattering would like to have consistent output type regardless
175 // of the correction method.
176 setProperty("OutputWorkspace", outWS);
177 } else if (method == "SampleAndContainer") {
178 //-- Setup output workspace
179 // set the OutputWorkspace to a workspace group with two workspaces:
180 // ${OutputWorkspace}_containerOnly
181 // ${OutputWorkspace}_sampleAndContainer
182 // 1. container only
183 API::MatrixWorkspace_sptr ws_containerOnly = create<HistoWorkspace>(*m_inputWS);
184 ws_containerOnly->setYUnit(""); // Need to explicitly set YUnit to nothing
185 ws_containerOnly->setDistribution(true); // The output of this is a distribution
186 ws_containerOnly->setYUnitLabel("Multiple Scattering Correction factor");
187 const auto &containerShape = m_inputWS->sample().getEnvironment().getContainer();
188 calculateSingleComponent(ws_containerOnly, containerShape, m_containerElementSize);
189 // 2. sample and container
190 API::MatrixWorkspace_sptr ws_sampleAndContainer = create<HistoWorkspace>(*m_inputWS);
191 ws_sampleAndContainer->setYUnit(""); // Need to explicitly set YUnit to nothing
192 ws_sampleAndContainer->setDistribution(true); // The output of this is a distribution
193 ws_sampleAndContainer->setYUnitLabel("Multiple Scattering Correction factor");
194 calculateSampleAndContainer(ws_sampleAndContainer);
195 //-- Package output to workspace group
196 const std::string outWSName = getProperty("OutputWorkspace");
197 std::vector<std::string> names;
198 names.emplace_back(outWSName + "_containerOnly");
199 API::AnalysisDataService::Instance().addOrReplace(names.back(), ws_containerOnly);
200 names.emplace_back(outWSName + "_sampleAndContainer");
201 API::AnalysisDataService::Instance().addOrReplace(names.back(), ws_sampleAndContainer);
202 // group
203 auto group = createChildAlgorithm("GroupWorkspaces");
204 group->initialize();
205 group->setProperty("InputWorkspaces", names);
206 group->setProperty("OutputWorkspace", outWSName);
207 group->execute();
208 API::WorkspaceGroup_sptr outWS = group->getProperty("OutputWorkspace");
209 //
210 setProperty("OutputWorkspace", outWS);
211 } else {
212 // With validator guarding the gate, this should never happen. However, just incase it
213 // does, we should throw an exception.
214 throw std::invalid_argument("Invalid method: " + method);
215 }
216}
217
223 // Get input workspace
224 m_inputWS = getProperty("InputWorkspace");
225
226 // Get the beam direction
227 m_beamDirection = m_inputWS->getInstrument()->getBeamDirection();
228
229 // Get the total number of wavelength points, default to use all if not specified
230 m_num_lambda = isDefault("NumberOfWavelengthPoints") ? static_cast<int64_t>(m_inputWS->blocksize())
231 : getProperty("NumberOfWavelengthPoints");
232 // -- while we're here, compute the step in bin number between two adjacent points
233 const auto specSize = static_cast<int64_t>(m_inputWS->blocksize());
234 m_xStep = std::max(int64_t(1), specSize / m_num_lambda); // Bin step between points to calculate
235 // -- notify the user of the bin step
236 std::ostringstream msg;
237 msg << "Numerical integration performed every " << m_xStep << " wavelength points";
238 g_log.information(msg.str());
239 g_log.information() << msg.str();
240
241 // Get the element size
242 m_sampleElementSize = getProperty("ElementSize"); // in mm
243 m_sampleElementSize = m_sampleElementSize * 1e-3; // convert to m
244 m_containerElementSize = getProperty("ContainerElementSize");
246}
247
256 const Geometry::IObject &shape, const double elementSize) {
257 const auto material = shape.material();
258 // Cache distances
259 // NOTE: cannot use IObject_sprt for sample shape as the getShape() method dereferenced
260 // the shared pointer upon returning.
261 MultipleScatteringCorrectionDistGraber distGraber(shape, elementSize);
262 distGraber.cacheLS1(m_beamDirection);
263
264 const int64_t numVolumeElements = distGraber.m_numVolumeElements;
265
266 // calculate distance within material from source to scattering point
267 std::vector<double> LS1s(numVolumeElements, 0.0);
268 calculateLS1s(distGraber, LS1s, shape);
269
270 // L12 is independent from the detector, therefore can be cached outside
271 // - L12 is a upper off-diagonal matrix
272 // NOTE: if the sample size/volume is too large, we might need to use openMP
273 // to parallelize the calculation
274 const int64_t len_l12 = numVolumeElements * (numVolumeElements - 1) / 2;
275 std::vector<double> L12s(len_l12, 0.0);
276 calculateL12s(distGraber, L12s, shape);
277
278 // L2D needs to be calculated w.r.t the detector
279
280 // compute the prefactor for multiple scattering correction factor Delta
281 // Delta = totScatterCoeff * A2/A1
282 // NOTE: Unit is important
283 // rho in 1/A^3, and sigma_s in 1/barns (1e-8 A^(-2))
284 // so rho * sigma_s = 1e-8 A^(-1) = 100 meters
285 // A2/A1 gives length in meters
286 const double rho = material.numberDensityEffective();
287 const double sigma_s = material.totalScatterXSection();
288 const double unit_scaling = 1e2;
289 const double totScatterCoeff = rho * sigma_s * unit_scaling;
290
291 // Calculate one detector at a time
292 const auto &spectrumInfo = m_inputWS->spectrumInfo();
293 const auto numHists = static_cast<int64_t>(m_inputWS->getNumberHistograms());
294 const auto specSize = static_cast<int64_t>(m_inputWS->blocksize());
295 Progress prog(this, 0.0, 1.0, numHists);
296 // -- loop over the spectra/detectors
298 for (int64_t workspaceIndex = 0; workspaceIndex < numHists; ++workspaceIndex) {
300 // locate the spectrum and its detector
301 if (!spectrumInfo.hasDetectors(workspaceIndex)) {
302 g_log.information() << "Spectrum " << workspaceIndex << " does not have a detector defined for it\n";
303 continue;
304 }
305 const auto &det = spectrumInfo.detector(workspaceIndex);
306
307 // compute L2D
308 std::vector<double> L2Ds(numVolumeElements, 0.0);
309 calculateL2Ds(distGraber, det, L2Ds, shape);
310
311 const auto wavelengths = m_inputWS->points(workspaceIndex);
312 // these need to have the minus sign applied still
313
314 const auto LinearCoefAbs = material.linearAbsorpCoef(wavelengths.cbegin(), wavelengths.cend());
315
316 auto &output = outws->mutableY(workspaceIndex);
317 // -- loop over the wavelength points every m_xStep
318 for (int64_t wvBinsIndex = 0; wvBinsIndex < specSize; wvBinsIndex += m_xStep) {
319 double A1 = 0.0;
320 double A2 = 0.0;
321
322 pairWiseSum(A1, A2, -LinearCoefAbs[wvBinsIndex], distGraber, LS1s, L12s, L2Ds, 0, numVolumeElements);
323
324 // compute the correction factor
325 // NOTE: prefactor, totScatterCoeff, is pre-calculated outside the loop (see above)
326 output[wvBinsIndex] = totScatterCoeff / (4 * M_PI) * (A2 / A1);
327
328 // debug output
329#ifndef NDEBUG
330 std::ostringstream msg_debug;
331 msg_debug << "Det_" << workspaceIndex << "@spectrum_" << wvBinsIndex << '\n'
332 << "\trho = " << rho << ", sigma_s = " << sigma_s << '\n'
333 << "\tA1 = " << A1 << '\n'
334 << "\tA2 = " << A2 << '\n'
335 << "\tms_factor = " << output[wvBinsIndex] << '\n';
336 g_log.notice(msg_debug.str());
337#endif
338
339 // Make certain that last point is calculated
340 if (m_xStep > 1 && wvBinsIndex + m_xStep >= specSize && wvBinsIndex + 1 != specSize) {
341 wvBinsIndex = specSize - m_xStep - 1;
342 }
343 }
344
345 // Interpolate linearly between points separated by m_xStep,
346 // last point required
347 if (m_xStep > 1) {
348 auto histNew = outws->histogram(workspaceIndex);
349 interpolateLinearInplace(histNew, m_xStep);
350 outws->setHistogram(workspaceIndex, histNew);
351 }
352
353 prog.report();
354
356 }
358
359 g_log.notice() << "finished integration.\n";
360}
361
369 // retrieve container related properties as they are relevant now
370 const auto &sample = m_inputWS->sample();
371 const auto &sampleMaterial = sample.getShape().material();
372 const auto &containerMaterial = sample.getEnvironment().getContainer().material();
373
374 // get the sample and container shapes
375 const auto &sampleShape = sample.getShape();
376 const auto &containerShape = sample.getEnvironment().getContainer();
377
378 MultipleScatteringCorrectionDistGraber distGraberSample(sampleShape, m_sampleElementSize);
379 distGraberSample.cacheLS1(m_beamDirection);
380 MultipleScatteringCorrectionDistGraber distGraberContainer(containerShape, m_containerElementSize);
381 distGraberContainer.cacheLS1(m_beamDirection);
382
383 // useful info to have
384 const int64_t numVolumeElementsSample = distGraberSample.m_numVolumeElements;
385 const int64_t numVolumeElementsContainer = distGraberContainer.m_numVolumeElements;
386 const int64_t numVolumeElements = numVolumeElementsSample + numVolumeElementsContainer;
387 g_log.information() << "numVolumeElementsSample=" << numVolumeElementsSample
388 << ", numVolumeElementsContainer=" << numVolumeElementsContainer << "\n";
389
390 // Total elements = [container_elements] + [sample_elements]
391 // Schematic for scattering element i (*)
392 // | \ / |
393 // | container \ sample / container |
394 // | \ / |
395 // | ---LS1_container[i] --- \ LS1_sample[i] * L2D_sample[i] / ---L2D_container[i] --- |
396 // | \ / |
397
398 // LS1 can be cached here, but L2D must be calculated within the loop of spectra
399 std::vector<double> LS1_container(numVolumeElements, 0.0);
400 std::vector<double> LS1_sample(numVolumeElements, 0.0);
401 calculateLS1s(distGraberContainer, distGraberSample, LS1_container, LS1_sample, containerShape, sampleShape);
402
403 // cache L12 for both sample and container
404 // L12 is a upper off-diagonal matrix from the hybrid of container and sample.
405 // e.g. container: i, j, k
406 // sample: l, m, n
407 // hybrid: i, j, k, l, m, n
408 // L12:
409 // i j k l m n
410 // -------------------------------------------------
411 // i | x L12[0] L12[1] L12[2] L12[3] L12[4]
412 // j | x x L12[5] L12[6] L12[7] L12[8]
413 // k | x x x L12[9] L12[10] L12[11]
414 // l | x x x x L12[12] L12[13]
415 // m | x x x x x L12[14]
416 // n | x x x x x x
417 const int64_t len_l12 = numVolumeElements * (numVolumeElements - 1) / 2;
418 std::vector<double> L12_container(len_l12, 0.0);
419 std::vector<double> L12_sample(len_l12, 0.0);
420 calculateL12s(distGraberContainer, distGraberSample, L12_container, L12_sample, containerShape, sampleShape);
421#ifndef NDEBUG
422 for (size_t i = 0; i < size_t(numVolumeElements); ++i) {
423 for (size_t j = i + 1; j < size_t(numVolumeElements); ++j) {
424 const auto idx = calcLinearIdxFromUpperTriangular(numVolumeElements, i, j);
425 const auto l12 = L12_container[idx] + L12_sample[idx];
426 if (l12 < 1e-9) {
427 g_log.notice() << "L12_container(" << i << "," << j << ")=" << L12_container[idx] << '\n'
428 << "L12_sample(" << i << "," << j << ")=" << L12_sample[idx] << '\n';
429 }
430 }
431 }
432#endif
433
434 // cache the elementsVolumes
435 std::vector<double> elementVolumes(distGraberContainer.m_elementVolumes.begin(),
436 distGraberContainer.m_elementVolumes.end());
437 elementVolumes.insert(elementVolumes.end(), distGraberSample.m_elementVolumes.begin(),
438 distGraberSample.m_elementVolumes.end());
439#ifndef NDEBUG
440 for (size_t i = 0; i < elementVolumes.size(); ++i) {
441 if (elementVolumes[i] < 1e-16) {
442 g_log.notice() << "Element_" << i << " has near zero volume: " << elementVolumes[i] << '\n';
443 }
444 }
445 g_log.notice() << "V_container = "
446 << std::accumulate(distGraberContainer.m_elementVolumes.begin(),
447 distGraberContainer.m_elementVolumes.end(), 0.0)
448 << '\n'
449 << "V_sample = "
450 << std::accumulate(distGraberSample.m_elementVolumes.begin(), distGraberSample.m_elementVolumes.end(),
451 0.0)
452 << '\n';
453#endif
454
455 // NOTE: Unit is important
456 // rho in 1/A^3, and sigma_s in 1/barns (1e-8 A^(-2))
457 // so rho * sigma_s = 1e-8 A^(-1) = 100 meters
458 // A2/A1 gives length in meters
459 const double unit_scaling = 1e2;
460 const double rho_sample = sampleMaterial.numberDensityEffective();
461 const double sigma_s_sample = sampleMaterial.totalScatterXSection();
462 const double rho_container = containerMaterial.numberDensityEffective();
463 const double sigma_s_container = containerMaterial.totalScatterXSection();
464 const double totScatterCoef_container = rho_container * sigma_s_container * unit_scaling;
465 const double totScatterCoef_sample = rho_sample * sigma_s_sample * unit_scaling;
466
467 // Compute the multiple scattering factor: one detector at a time
468 const auto &spectrumInfo = m_inputWS->spectrumInfo();
469 const auto numHists = static_cast<int64_t>(m_inputWS->getNumberHistograms());
470 const auto specSize = static_cast<int64_t>(m_inputWS->blocksize());
471 Progress prog(this, 0.0, 1.0, numHists);
472 // -- loop over the spectra/detectors
474 for (int64_t workspaceIndex = 0; workspaceIndex < numHists; ++workspaceIndex) {
476 // locate the spectrum and its detector
477 if (!spectrumInfo.hasDetectors(workspaceIndex)) {
478 g_log.information() << "Spectrum " << workspaceIndex << " does not have a detector defined for it\n";
479 continue;
480 }
481 const auto &det = spectrumInfo.detector(workspaceIndex);
482
483 // calculate L2D (2_element -> detector)
484 // 1. container
485 std::vector<double> L2D_container(numVolumeElements, 0.0);
486 std::vector<double> L2D_sample(numVolumeElements, 0.0);
487 calculateL2Ds(distGraberContainer, distGraberSample, det, L2D_container, L2D_sample, containerShape, sampleShape);
488
489 // prepare material wise linear coefficient
490 const auto wavelengths = m_inputWS->points(workspaceIndex);
491 const auto sampleLinearCoefAbs = sampleMaterial.linearAbsorpCoef(wavelengths.cbegin(), wavelengths.cend());
492 const auto containerLinearCoefAbs = containerMaterial.linearAbsorpCoef(wavelengths.cbegin(), wavelengths.cend());
493
494 auto &output = outws->mutableY(workspaceIndex);
495 for (int64_t wvBinsIndex = 0; wvBinsIndex < specSize; wvBinsIndex += m_xStep) {
496 double A1 = 0.0;
497 double A2 = 0.0;
498
499 // compute the multiple scattering correction factor Delta
500 pairWiseSum(A1, A2, // output values
501 -containerLinearCoefAbs[wvBinsIndex], -sampleLinearCoefAbs[wvBinsIndex], // absorption coefficient
502 numVolumeElementsContainer, numVolumeElements, // number of elements for checking type
503 totScatterCoef_container, totScatterCoef_sample, // volumes
504 elementVolumes, // source -> 1st scattering element
505 LS1_container, LS1_sample, // 1st -> 2nd scattering element
506 L12_container, L12_sample, // 2nd scattering element -> detector
507 L2D_container, L2D_sample, // starting element idx, ending element idx
508 0, numVolumeElements);
509
510 output[wvBinsIndex] = (A2 / A1) / (4.0 * M_PI);
511
512 // debug output
513#ifndef NDEBUG
514 std::ostringstream msg_debug;
515 msg_debug << "Det_" << workspaceIndex << "@spectrum_" << wvBinsIndex << '\n'
516 << "-containerLinearCoefAbs[wvBinsIndex] = " << -containerLinearCoefAbs[wvBinsIndex] << '\n'
517 << "-sampleLinearCoefAbs[wvBinsIndex] = " << -sampleLinearCoefAbs[wvBinsIndex] << '\n'
518 << "numVolumeElementsContainer = " << numVolumeElementsContainer << '\n'
519 << "numVolumeElements = " << numVolumeElements << '\n'
520 << "totScatterCoef_container = " << totScatterCoef_container << '\n'
521 << "totScatterCoef_sample = " << totScatterCoef_sample << '\n'
522 << "\tA1 = " << A1 << '\n'
523 << "\tA2 = " << A2 << '\n'
524 << "\tms_factor = " << output[wvBinsIndex] << '\n';
525 g_log.notice(msg_debug.str());
526#endif
527
528 // Make certain that last point is calculated
529 if (m_xStep > 1 && wvBinsIndex + m_xStep >= specSize && wvBinsIndex + 1 != specSize) {
530 wvBinsIndex = specSize - m_xStep - 1;
531 }
532 }
533 // Interpolate linearly between points separated by m_xStep,
534 // last point required
535 if (m_xStep > 1) {
536 auto histNew = outws->histogram(workspaceIndex);
537 interpolateLinearInplace(histNew, m_xStep);
538 outws->setHistogram(workspaceIndex, histNew);
539 }
540 prog.report();
542 }
544 g_log.notice() << "finished integration.\n";
545}
546
555 std::vector<double> &LS1s, //
556 const Geometry::IObject &shape) const {
557 const auto &sourcePos = m_inputWS->getInstrument()->getSource()->getPos();
558 const int64_t numVolumeElements = distGraber.m_numVolumeElements;
559 Track trackerLS1(V3D{0, 0, 1}, V3D{0, 0, 1});
560
561 for (int64_t idx = 0; idx < numVolumeElements; ++idx) {
562 const auto pos = distGraber.m_elementPositions[idx];
563 const auto vec = getDirection(pos, sourcePos);
564 //
565 trackerLS1.reset(pos, vec);
566 trackerLS1.clearIntersectionResults();
567 LS1s[idx] = getDistanceInsideObject(shape, trackerLS1);
568 }
569}
570
582 const MultipleScatteringCorrectionDistGraber &distGraberSample, //
583 std::vector<double> &LS1sContainer, //
584 std::vector<double> &LS1sSample, //
585 const Geometry::IObject &shapeContainer, //
586 const Geometry::IObject &shapeSample) const {
587 // Total elements = [container_elements] + [sample_elements]
588 // Schematic for scattering element i (*)
589 // | \ / |
590 // | container \ sample / container |
591 // | \ / |
592 // | ---LS1_container[i] --- \ LS1_sample[i] * L2D_sample[i] / ---L2D_container[i] --- |
593 // | \ / |
594 const int64_t numVolumeElementsSample = distGraberSample.m_numVolumeElements;
595 const int64_t numVolumeElementsContainer = distGraberContainer.m_numVolumeElements;
596 const int64_t numVolumeElements = numVolumeElementsSample + numVolumeElementsContainer;
597 const auto sourcePos = m_inputWS->getInstrument()->getSource()->getPos();
598 Track trackerLS1(V3D{0, 0, 1}, V3D{0, 0, 1}); // reusable tracker for calculating LS1
599 for (int64_t idx = 0; idx < numVolumeElements; ++idx) {
600 const auto pos = idx < numVolumeElementsContainer
601 ? distGraberContainer.m_elementPositions[idx]
602 : distGraberSample.m_elementPositions[idx - numVolumeElementsContainer];
603 const auto vec = getDirection(pos, sourcePos);
604 //
605 trackerLS1.reset(pos, vec);
606 trackerLS1.clearIntersectionResults();
607 LS1sContainer[idx] = getDistanceInsideObject(shapeContainer, trackerLS1);
608 trackerLS1.reset(pos, vec);
609 trackerLS1.clearIntersectionResults();
610 LS1sSample[idx] = getDistanceInsideObject(shapeSample, trackerLS1);
611#ifndef NDEBUG
612 // debug
613 std::ostringstream msg_debug;
614 msg_debug << "idx=" << idx << ", pos=" << pos << ", vec=" << vec << '\n';
615 if (idx < numVolumeElementsContainer) {
616 msg_debug << "Container element " << idx << '\n';
617 } else {
618 msg_debug << "Sample element " << idx - numVolumeElementsContainer << '\n';
619 }
620 msg_debug << "LS1_container=" << LS1sContainer[idx] << ", LS1_sample=" << LS1sSample[idx] << '\n';
621 g_log.notice(msg_debug.str());
622#endif
623 }
624}
625
634 std::vector<double> &L12s, //
635 const Geometry::IObject &shape) {
636 const int64_t numVolumeElements = distGraber.m_numVolumeElements;
637 // L12 is independent from the detector, therefore can be cached outside
638 // - L12 is a upper off-diagonal matrix
639 // NOTE: if the sample size/volume is too large, we might need to use openMP
640 // to parallelize the calculation
641
643 for (int64_t indexTo = 0; indexTo < numVolumeElements; ++indexTo) {
645
646 const auto posTo = distGraber.m_elementPositions[indexTo];
647 Track track(posTo, V3D{0, 0, 1}); // take object creation out of the loop
648 for (int64_t indexFrom = indexTo + 1; indexFrom < numVolumeElements; ++indexFrom) {
649 // where in the final result to update, e.g.
650 // x 1 2 3
651 // x x 4 5
652 // x x x 6
653 // x x x x
654 const int64_t idx = calcLinearIdxFromUpperTriangular(numVolumeElements, indexTo, indexFrom);
655
656 const auto posFrom = distGraber.m_elementPositions[indexFrom];
657 const V3D unitVector = getDirection(posFrom, posTo);
658
659 // reset information in the Track and calculate distance
660 track.reset(posFrom, unitVector);
662 const auto rayLengthOne1 = getDistanceInsideObject(shape, track);
663
664 track.reset(posTo, unitVector);
666 const auto rayLengthOne2 = getDistanceInsideObject(shape, track);
667 // getDistanceInsideObject returns the line segment inside the shape from the given ray (defined by track)
668 // therefore, the distance between the two point (element) can be found by the difference of the two line
669 // segments.
670 L12s[idx] = checkzero(rayLengthOne1 - rayLengthOne2);
671 }
672
674 }
676}
677
679 const MultipleScatteringCorrectionDistGraber &distGraberSample, //
680 std::vector<double> &L12sContainer, //
681 std::vector<double> &L12sSample, //
682 const Geometry::IObject &shapeContainer, //
683 const Geometry::IObject &shapeSample) {
684 const int64_t numVolumeElementsSample = distGraberSample.m_numVolumeElements;
685 const int64_t numVolumeElementsContainer = distGraberContainer.m_numVolumeElements;
686 const int64_t numVolumeElements = numVolumeElementsSample + numVolumeElementsContainer;
687 // L12 is a upper off-diagonal matrix from the hybrid of container and sample.
688 // e.g. container: a, b, c
689 // sample: alpha, beta, gamma
690 // hybrid: a, b, c, alpha, beta, gamma
691 // L12:
692 // a b c alpha beta gamma
693 // -------------------------------------------------
694 // a | x L12[0] L12[1] L12[2] L12[3] L12[4]
695 // b | x x L12[5] L12[6] L12[7] L12[8]
696 // c | x x x L12[9] L12[10] L12[11]
697 // alpha | x x x x L12[12] L12[13]
698 // beta | x x x x x L12[14]
699 // gamma | x x x x x x
700
702 for (int64_t indexTo = 0; indexTo < numVolumeElements; ++indexTo) {
704 // find the position of first scattering element (container or sample)
705 const auto posTo = indexTo < numVolumeElementsContainer
706 ? distGraberContainer.m_elementPositions[indexTo]
707 : distGraberSample.m_elementPositions[indexTo - numVolumeElementsContainer];
708 // NOTE:
709 // We are using Track to ensure that the distance calculated are within the material
710 Track track(posTo, V3D{0, 0, 1}); // reusable track, eco-friendly
711 for (int64_t indexFrom = indexTo + 1; indexFrom < numVolumeElements; ++indexFrom) {
712 const int64_t idx = calcLinearIdxFromUpperTriangular(numVolumeElements, indexTo, indexFrom);
713 const auto posFrom = indexFrom < numVolumeElementsContainer
714 ? distGraberContainer.m_elementPositions[indexFrom]
715 : distGraberSample.m_elementPositions[indexFrom - numVolumeElementsContainer];
716 const V3D unitVector = getDirection(posFrom, posTo);
717
718 // combined
720 track.reset(posFrom, unitVector);
721 const auto rayLen1_container = getDistanceInsideObject(shapeContainer, track);
723 track.reset(posFrom, unitVector);
724 const auto rayLen1_sample = getDistanceInsideObject(shapeSample, track);
726 track.reset(posTo, unitVector);
727 const auto rayLen2_container = getDistanceInsideObject(shapeContainer, track);
729 track.reset(posTo, unitVector);
730 const auto rayLen2_sample = getDistanceInsideObject(shapeSample, track);
731 //
732 L12sContainer[idx] = checkzero(rayLen1_container - rayLen2_container);
733 L12sSample[idx] = checkzero(rayLen1_sample - rayLen2_sample);
734 }
736 }
738}
739
749 const IDetector &detector, std::vector<double> &L2Ds,
750 const Geometry::IObject &shape) const {
751 V3D detectorPos(detector.getPos());
752 if (detector.nDets() > 1) {
753 // We need to make sure this is right for grouped detectors - should use
754 // average theta & phi
755 detectorPos.spherical(detectorPos.norm(), detector.getTwoTheta(V3D(), V3D(0, 0, 1)) * RAD2DEG,
756 detector.getPhi() * RAD2DEG);
757 }
758
759 // calculate the distance between the detector and the shape (sample or container)
760 Track TwoToDetector(distGraber.m_elementPositions[0], V3D{1, 0, 0}); // take object creation out of the loop
761 for (size_t i = 0; i < distGraber.m_elementPositions.size(); ++i) {
762 const auto &elementPos = distGraber.m_elementPositions[i];
763 TwoToDetector.reset(elementPos, getDirection(elementPos, detectorPos));
764 TwoToDetector.clearIntersectionResults();
765 // find distance in sample
766 L2Ds[i] = getDistanceInsideObject(shape, TwoToDetector);
767 }
768}
769
771 const MultipleScatteringCorrectionDistGraber &distGraberSample, //
772 const IDetector &detector, //
773 std::vector<double> &container_L2Ds, //
774 std::vector<double> &sample_L2Ds, //
775 const Geometry::IObject &shapeContainer, //
776 const Geometry::IObject &shapeSample) const {
777 V3D detectorPos(detector.getPos());
778 if (detector.nDets() > 1) {
779 // We need to make sure this is right for grouped detectors - should use
780 // average theta & phi
781 detectorPos.spherical(detectorPos.norm(), detector.getTwoTheta(V3D(), V3D(0, 0, 1)) * RAD2DEG,
782 detector.getPhi() * RAD2DEG);
783 }
784
785 const int64_t numVolumeElementsSample = distGraberSample.m_numVolumeElements;
786 const int64_t numVolumeElementsContainer = distGraberContainer.m_numVolumeElements;
787 const int64_t numVolumeElements = numVolumeElementsSample + numVolumeElementsContainer;
788 // Total elements = [container_elements] + [sample_elements]
789 // Schematic for scattering element i (*)
790 // | \ / |
791 // | container \ sample / container |
792 // | \ / |
793 // | ---LS1_container[i] --- \ LS1_sample[i] * L2D_sample[i] / ---L2D_container[i] --- |
794 // | \ / |
795 // L2D must be calculated within the loop of spectra, so we cannot use OpenMP here
796 Track trackerL2D(V3D{0, 0, 1}, V3D{0, 0, 1}); // reusable tracker for calculating L2D
797 for (int64_t idx = 0; idx < numVolumeElements; ++idx) {
798 const auto pos = idx < numVolumeElementsContainer
799 ? distGraberContainer.m_elementPositions[idx]
800 : distGraberSample.m_elementPositions[idx - numVolumeElementsContainer];
801 const auto vec = getDirection(pos, detectorPos);
802 //
803 trackerL2D.reset(pos, vec);
804 trackerL2D.clearIntersectionResults();
805 container_L2Ds[idx] = getDistanceInsideObject(shapeContainer, trackerL2D);
806 trackerL2D.reset(pos, vec);
807 trackerL2D.clearIntersectionResults();
808 sample_L2Ds[idx] = getDistanceInsideObject(shapeSample, trackerL2D);
809 }
810}
811
812void MultipleScatteringCorrection::pairWiseSum(double &A1, double &A2, //
813 const double linearCoefAbs, //
814 const MultipleScatteringCorrectionDistGraber &distGraber, //
815 const std::vector<double> &LS1s, //
816 const std::vector<double> &L12s, //
817 const std::vector<double> &L2Ds, //
818 const int64_t startIndex, const int64_t endIndex) const {
819 if (endIndex - startIndex > MAX_INTEGRATION_LENGTH) {
820 int64_t middle = findMiddle(startIndex, endIndex);
821
822 // recursive to process upper and lower part
823 pairWiseSum(A1, A2, linearCoefAbs, distGraber, LS1s, L12s, L2Ds, startIndex, middle);
824 pairWiseSum(A1, A2, linearCoefAbs, distGraber, LS1s, L12s, L2Ds, middle, endIndex);
825 } else {
826 // perform the integration
827 const auto &elementVolumes = distGraber.m_elementVolumes;
828 const auto nElements = distGraber.m_numVolumeElements;
829 for (int64_t i = startIndex; i < endIndex; ++i) {
830 // compute A1
831 double exponent = (LS1s[i] + L2Ds[i]) * linearCoefAbs;
832 A1 += exp(exponent) * elementVolumes[i];
833 // compute A2
834 double a2 = 0.0;
835 for (int64_t j = 0; j < int64_t(nElements); ++j) {
836 if (i == j) {
837 // skip self (second order scattering must happen in a different element)
838 continue;
839 }
840 // L12 is a pre-computed vector, therefore we can use the index directly
841 size_t idx_l12 = i < j ? calcLinearIdxFromUpperTriangular(nElements, i, j)
842 : calcLinearIdxFromUpperTriangular(nElements, j, i);
843 // compute a2 component
844 const auto l12 = L12s[idx_l12];
845 if (l12 > 0.0) {
846 exponent = (LS1s[i] + L12s[idx_l12] + L2Ds[j]) * linearCoefAbs;
847 a2 += exp(exponent) * elementVolumes[j] / (L12s[idx_l12] * L12s[idx_l12]);
848 }
849 }
850 A2 += a2 * elementVolumes[i];
851 }
852 }
853}
854
856 double &A1, double &A2, //
857 const double linearCoefAbsContainer, const double linearCoefAbsSample, //
858 const int64_t numVolumeElementsContainer, const int64_t numVolumeElementsTotal, //
859 const double totScatterCoefContainer, // rho * sigma_s * unit_scaling
860 const double totScatterCoefSample, // unit_scaling = 100
861 const std::vector<double> &elementVolumes, //
862 const std::vector<double> &LS1sContainer, const std::vector<double> &LS1sSample, // source -> 1_element
863 const std::vector<double> &L12sContainer, const std::vector<double> &L12sSample, // 1_element -> 2_element
864 const std::vector<double> &L2DsContainer, const std::vector<double> &L2DsSample, // 2_element -> detector
865 const int64_t startIndex, const int64_t endIndex) const {
866 if (endIndex - startIndex > MAX_INTEGRATION_LENGTH) {
867 int64_t middle = findMiddle(startIndex, endIndex);
868 // recursive to process upper and lower part
869 pairWiseSum(A1, A2, linearCoefAbsContainer, linearCoefAbsSample, numVolumeElementsContainer, numVolumeElementsTotal,
870 totScatterCoefContainer, totScatterCoefSample, elementVolumes, LS1sContainer, LS1sSample, L12sContainer,
871 L12sSample, L2DsContainer, L2DsSample, startIndex, middle);
872 pairWiseSum(A1, A2, linearCoefAbsContainer, linearCoefAbsSample, numVolumeElementsContainer, numVolumeElementsTotal,
873 totScatterCoefContainer, totScatterCoefSample, elementVolumes, LS1sContainer, LS1sSample, L12sContainer,
874 L12sSample, L2DsContainer, L2DsSample, middle, endIndex);
875 } else {
876 // perform the integration
877 for (int64_t i = startIndex; i < endIndex; ++i) {
878 const double factor_i = i > numVolumeElementsContainer ? totScatterCoefSample : totScatterCoefContainer;
879 // compute A1
880 double exponent = (LS1sContainer[i] + L2DsContainer[i]) * linearCoefAbsContainer +
881 (LS1sSample[i] + L2DsSample[i]) * linearCoefAbsSample;
882 A1 += exp(exponent) * factor_i * elementVolumes[i];
883 // compute A2
884 double a2 = 0.0;
885 for (int64_t j = 0; j < numVolumeElementsTotal; ++j) {
886 if (i == j) {
887 // skip self (second order scattering must happen in a different element)
888 continue;
889 }
890 const double factor_j = j > numVolumeElementsContainer ? totScatterCoefSample : totScatterCoefContainer;
891 // L12 is a pre-computed vector, therefore we can use the index directly
892 size_t idx_l12 = i < j ? calcLinearIdxFromUpperTriangular(numVolumeElementsTotal, i, j)
893 : calcLinearIdxFromUpperTriangular(numVolumeElementsTotal, j, i);
894 const double l12 = L12sContainer[idx_l12] + L12sSample[idx_l12];
895 if (l12 > 0.0) {
896 exponent = (LS1sContainer[i] + L12sContainer[idx_l12] + L2DsContainer[j]) * linearCoefAbsContainer + //
897 (LS1sSample[i] + L12sSample[idx_l12] + L2DsSample[j]) * linearCoefAbsSample;
898 a2 += exp(exponent) * factor_j * elementVolumes[j] / (l12 * l12);
899 }
900 }
901 A2 += a2 * factor_i * elementVolumes[i];
902 }
903 }
904}
905
906} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
#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.
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
virtual std::shared_ptr< Algorithm > createChildAlgorithm(const std::string &name, const double startProgress=-1., const double endProgress=-1., const bool enableLogging=true, const int &version=-1)
Create a Child Algorithm.
Definition: Algorithm.cpp:842
Kernel::Logger & g_log
Definition: Algorithm.h:451
bool isDefault(const std::string &name) const
Definition: Algorithm.cpp:2084
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...
MultipleScatteringCorrectionDistGraber : This is a helper class to calculate the distance from source...
std::vector< Kernel::V3D > m_elementPositions
Cached element positions.
void cacheLS1(const Mantid::Kernel::V3D &beamDirection)
pre-calculate the distance from source to L1 for all the voxels in the sample
void calculateLS1s(const MultipleScatteringCorrectionDistGraber &distGraber, std::vector< double > &LS1s, const Geometry::IObject &shape) const
compute LS1s within given shape for single component case
void calculateSingleComponent(const API::MatrixWorkspace_sptr &outws, const Geometry::IObject &shape, const double elementSize)
calculate the correction factor per detector for sample only case
void parseInputs()
parse and assign corresponding values from input properties
void pairWiseSum(double &A1, double &A2, const double linearCoefAbs, const MultipleScatteringCorrectionDistGraber &distGraber, const std::vector< double > &LS1s, const std::vector< double > &L12s, const std::vector< double > &L2Ds, const int64_t startIndex, const int64_t endIndex) const
double m_sampleElementSize
The size of the integration element for sample in meters.
void calculateSampleAndContainer(const API::MatrixWorkspace_sptr &outws)
calculate the multiple scattering factor (0, 1) for sample and container case
API::MatrixWorkspace_sptr m_inputWS
A pointer to the input workspace.
void calculateL2Ds(const MultipleScatteringCorrectionDistGraber &distGraber, const IDetector &detector, std::vector< double > &L2Ds, const Geometry::IObject &shape) const
Calculate distance between exiting element to the detector for single component case.
int64_t m_num_lambda
The number of points in wavelength, the rest is interpolated linearly.
int64_t m_xStep
The step in bin number between adjacent points for linear interpolation.
void init() override
interface initialisation method
std::map< std::string, std::string > validateInputs() override
validate the inputs
void calculateL12s(const MultipleScatteringCorrectionDistGraber &distGraber, std::vector< double > &L12s, const Geometry::IObject &shape)
calculate L12 for single component case
double m_containerElementSize
the size of the integration element for container in meters
virtual Kernel::V3D getPos() const =0
Get the position of the IComponent. Tree structure is traverse through the.
Interface class for detector objects.
Definition: IDetector.h:43
virtual double getPhi() const =0
Gives the phi of this detector object in radians.
virtual std::size_t nDets() const =0
Get the number of physical detectors this object represents.
virtual double getTwoTheta(const Kernel::V3D &observer, const Kernel::V3D &axis) const =0
Gives the angle of this detector object with respect to an axis.
IObject : Interface for geometry objects.
Definition: IObject.h:41
virtual int interceptSurface(Geometry::Track &) const =0
virtual const Kernel::Material & material() const =0
Defines a track as a start point and a direction.
Definition: Track.h:165
void clearIntersectionResults()
Clear the current set of intersection results.
Definition: Track.cpp:55
double totalDistInsideObject() const
Returns the sum of all the links distInsideObject in the track.
Definition: Track.cpp:254
void reset(const Kernel::V3D &startPoint, const Kernel::V3D &direction)
Set a starting point and direction.
Definition: Track.cpp:45
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void notice(const std::string &msg)
Logs at notice level.
Definition: Logger.cpp:95
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
Definition: ProgressBase.h:51
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
void spherical(const double R, const double theta, const double phi) noexcept
Sets the vector position based on spherical coordinates.
Definition: V3D.cpp:57
double norm() const noexcept
Definition: V3D.h:263
std::shared_ptr< WorkspaceGroup > WorkspaceGroup_sptr
shared pointer to Mantid::API::WorkspaceGroup
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
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
constexpr int EMPTY_INT() noexcept
Returns what we consider an "empty" integer within a property.
Definition: EmptyValues.h:25
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
Definition: EmptyValues.h:43
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54