Mantid
Loading...
Searching...
No Matches
MDNormSCD.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
11#include "MantidAPI/Run.h"
24
25namespace Mantid::MDAlgorithms {
26
29using namespace Mantid::DataObjects;
30using namespace Mantid::API;
31using namespace Mantid::Kernel;
32
33namespace {
34// function to compare two intersections (h,k,l,Momentum) by Momentum
35bool compareMomentum(const std::array<double, 4> &v1, const std::array<double, 4> &v2) { return (v1[3] < v2[3]); }
36} // namespace
37
38// Register the algorithm into the AlgorithmFactory
39DECLARE_ALGORITHM(MDNormSCD)
40
41
45 : m_normWS(), m_inputWS(), m_hmin(0.0f), m_hmax(0.0f), m_kmin(0.0f), m_kmax(0.0f), m_lmin(0.0f), m_lmax(0.0f),
46 m_hIntegrated(true), m_kIntegrated(true), m_lIntegrated(true), m_rubw(3, 3), m_kiMin(0.0), m_kiMax(EMPTY_DBL()),
47 m_hIdx(-1), m_kIdx(-1), m_lIdx(-1), m_hX(), m_kX(), m_lX(), m_samplePos(), m_beamDir() {}
48
50int MDNormSCD::version() const { return 1; }
51
53const std::string MDNormSCD::category() const { return "MDAlgorithms\\Normalisation"; }
54
56const std::string MDNormSCD::summary() const {
57 return "Calculate normalization for an MDEvent workspace for single crystal "
58 "diffraction.";
59}
60
62const std::string MDNormSCD::name() const { return "MDNormSCD"; }
63
68 declareProperty(std::make_unique<WorkspaceProperty<IMDEventWorkspace>>("InputWorkspace", "", Direction::Input),
69 "An input MDWorkspace.");
70
71 std::string dimChars = getDimensionChars();
72 // --------------- Axis-aligned properties
73 // ---------------------------------------
74 for (size_t i = 0; i < dimChars.size(); i++) {
75 std::string dim(" ");
76 dim[0] = dimChars[i];
77 std::string propName = "AlignedDim" + dim;
79 "Binning parameters for the " + Strings::toString(i) +
80 "th dimension.\n"
81 "Enter it as a comma-separated list of values with the format: "
82 "'name,minimum,maximum,number_of_bins'. Leave blank for NONE.");
83 }
84
85 auto fluxValidator = std::make_shared<CompositeValidator>();
86 fluxValidator->add<WorkspaceUnitValidator>("Momentum");
87 fluxValidator->add<InstrumentValidator>();
88 fluxValidator->add<CommonBinsValidator>();
89 auto solidAngleValidator = fluxValidator->clone();
90
91 declareProperty(std::make_unique<WorkspaceProperty<>>("FluxWorkspace", "", Direction::Input, fluxValidator),
92 "An input workspace containing momentum dependent flux.");
94 std::make_unique<WorkspaceProperty<>>("SolidAngleWorkspace", "", Direction::Input, solidAngleValidator),
95 "An input workspace containing momentum integrated vanadium "
96 "(a measure of the solid angle).");
97
98 declareProperty(std::make_unique<PropertyWithValue<bool>>("SkipSafetyCheck", false, Direction::Input),
99 "If set to true, the algorithm does "
100 "not check history if the workspace was modified since the"
101 "ConvertToMD algorithm was run, and assume that the elastic "
102 "mode is used.");
103
104 declareProperty(std::make_unique<WorkspaceProperty<IMDHistoWorkspace>>("TemporaryNormalizationWorkspace", "",
106 "An input MDHistoWorkspace used to accumulate normalization "
107 "from multiple MDEventWorkspaces. "
108 "If unspecified a blank MDHistoWorkspace will be created.");
109
110 declareProperty(std::make_unique<WorkspaceProperty<IMDHistoWorkspace>>("TemporaryDataWorkspace", "", Direction::Input,
112 "An input MDHistoWorkspace used to accumulate data from "
113 "multiple MDEventWorkspaces. If "
114 "unspecified a blank MDHistoWorkspace will be created.");
115
116 declareProperty(std::make_unique<WorkspaceProperty<Workspace>>("OutputWorkspace", "", Direction::Output),
117 "A name for the output data MDHistoWorkspace.");
118 declareProperty(std::make_unique<WorkspaceProperty<Workspace>>("OutputNormalizationWorkspace", "", Direction::Output),
119 "A name for the output normalization MDHistoWorkspace.");
120}
121
126 cacheInputs();
127 auto outputWS = binInputWS();
128 convention = Kernel::ConfigService::Instance().getString("Q.convention");
129 outputWS->setDisplayNormalization(Mantid::API::NoNormalization);
130 setProperty<Workspace_sptr>("OutputWorkspace", outputWS);
131 createNormalizationWS(*outputWS);
132 m_normWS->setDisplayNormalization(Mantid::API::NoNormalization);
133 setProperty("OutputNormalizationWorkspace", m_normWS);
134
135 m_numExptInfos = outputWS->getNumExperimentInfo();
136 // loop over all experiment infos
137 for (uint16_t expInfoIndex = 0; expInfoIndex < m_numExptInfos; expInfoIndex++) {
138 // Check for other dimensions if we could measure anything in the original
139 // data
140 bool skipNormalization = false;
141 const std::vector<coord_t> otherValues = getValuesFromOtherDimensions(skipNormalization, expInfoIndex);
142 const auto affineTrans = findIntergratedDimensions(otherValues, skipNormalization);
144
145 if (!skipNormalization) {
146 calculateNormalization(otherValues, affineTrans, expInfoIndex);
147 } else {
148 g_log.warning("Binning limits are outside the limits of the MDWorkspace. "
149 "Not applying normalization.");
150 }
151 m_accumulate = true;
152 }
153}
154
159 m_inputWS = getProperty("InputWorkspace");
160 bool skipCheck = getProperty("SkipSafetyCheck");
161 if (!skipCheck && inputEnergyMode() != "Elastic") {
162 throw std::invalid_argument("Invalid energy transfer mode. Algorithm "
163 "currently only supports elastic data.");
164 }
165 // Min/max dimension values
166 const auto hdim(m_inputWS->getDimension(0)), kdim(m_inputWS->getDimension(1)), ldim(m_inputWS->getDimension(2));
167 m_hmin = hdim->getMinimum();
168 m_kmin = kdim->getMinimum();
169 m_lmin = ldim->getMinimum();
170 m_hmax = hdim->getMaximum();
171 m_kmax = kdim->getMaximum();
172 m_lmax = ldim->getMaximum();
173
174 const auto &exptInfoZero = *(m_inputWS->getExperimentInfo(0));
175 auto source = exptInfoZero.getInstrument()->getSource();
176 auto sample = exptInfoZero.getInstrument()->getSample();
177 if (source == nullptr || sample == nullptr) {
179 "Instrument not sufficiently defined: failed to get source and/or "
180 "sample");
181 }
182 m_samplePos = sample->getPos();
183 m_beamDir = normalize(m_samplePos - source->getPos());
184}
185
190std::string MDNormSCD::inputEnergyMode() const {
191 const auto &history = m_inputWS->getHistory();
192 const size_t nalgs = history.size();
193 const auto &lastAlgorithm = history.lastAlgorithm();
194
195 std::string emode;
196 if (lastAlgorithm->name() == "ConvertToMD") {
197 emode = lastAlgorithm->getPropertyValue("dEAnalysisMode");
198 } else if ((lastAlgorithm->name() == "Load" || history.lastAlgorithm()->name() == "LoadMD") &&
199 history.getAlgorithmHistory(nalgs - 2)->name() == "ConvertToMD") {
200 // get dEAnalysisMode
201 PropertyHistories histvec = history.getAlgorithmHistory(nalgs - 2)->getProperties();
202 for (auto &hist : histvec) {
203 if (hist->name() == "dEAnalysisMode") {
204 emode = hist->value();
205 break;
206 }
207 }
208 } else {
209 throw std::invalid_argument("The last algorithm in the history of the "
210 "input workspace is not ConvertToMD");
211 }
212 return emode;
213}
214
221 const auto &props = getProperties();
222 auto binMD = createChildAlgorithm("BinMD", 0.0, 0.3);
223 binMD->setPropertyValue("AxisAligned", "1");
224 for (auto prop : props) {
225 const auto &propName = prop->name();
226 if (propName != "FluxWorkspace" && propName != "SolidAngleWorkspace" &&
227 propName != "TemporaryNormalizationWorkspace" && propName != "OutputNormalizationWorkspace" &&
228 propName != "SkipSafetyCheck") {
229 binMD->setPropertyValue(propName, prop->value());
230 }
231 }
232 binMD->executeAsChildAlg();
233 Workspace_sptr outputWS = binMD->getProperty("OutputWorkspace");
234 return std::dynamic_pointer_cast<MDHistoWorkspace>(outputWS);
235}
236
242 // Copy the MDHisto workspace, and change signals and errors to 0.
243 std::shared_ptr<IMDHistoWorkspace> tmp = this->getProperty("TemporaryNormalizationWorkspace");
244 m_normWS = std::dynamic_pointer_cast<MDHistoWorkspace>(tmp);
245 if (!m_normWS) {
246 m_normWS = dataWS.clone();
247 m_normWS->setTo(0., 0., 0.);
248 } else {
249 m_accumulate = true;
250 }
251}
252
261std::vector<coord_t> MDNormSCD::getValuesFromOtherDimensions(bool &skipNormalization, uint16_t expInfoIndex) const {
262 const auto &currentRun = m_inputWS->getExperimentInfo(expInfoIndex)->run();
263
264 std::vector<coord_t> otherDimValues;
265 for (size_t i = 3; i < m_inputWS->getNumDims(); i++) {
266 const auto dimension = m_inputWS->getDimension(i);
267 auto dimMin = static_cast<float>(dimension->getMinimum());
268 auto dimMax = static_cast<float>(dimension->getMaximum());
269 auto *dimProp = dynamic_cast<Kernel::TimeSeriesProperty<double> *>(currentRun.getProperty(dimension->getName()));
270 if (dimProp) {
271 auto value = static_cast<coord_t>(dimProp->firstValue());
272 otherDimValues.emplace_back(value);
273 // in the original MD data no time was spent measuring between dimMin and
274 // dimMax
275 if (value < dimMin || value > dimMax) {
276 skipNormalization = true;
277 }
278 }
279 }
280 return otherDimValues;
281}
282
292Kernel::Matrix<coord_t> MDNormSCD::findIntergratedDimensions(const std::vector<coord_t> &otherDimValues,
293 bool &skipNormalization) {
294 // Get indices of the original dimensions in the output workspace,
295 // and if not found, the corresponding dimension is integrated
296 Kernel::Matrix<coord_t> affineMat = m_normWS->getTransformFromOriginal(0)->makeAffineMatrix();
297
298 const size_t nrm1 = affineMat.numRows() - 1;
299 const size_t ncm1 = affineMat.numCols() - 1;
300 for (size_t row = 0; row < nrm1; row++) // affine matrix, ignore last row
301 {
302 const auto dimen = m_normWS->getDimension(row);
303 const auto dimMin(dimen->getMinimum()), dimMax(dimen->getMaximum());
304 if (affineMat[row][0] == 1.) {
305 m_hIntegrated = false;
306 m_hIdx = row;
307 m_hmin = std::max(m_hmin, dimMin);
308 m_hmax = std::min(m_hmax, dimMax);
309 if (m_hmin > dimMax || m_hmax < dimMin) {
310 skipNormalization = true;
311 }
312 }
313 if (affineMat[row][1] == 1.) {
314 m_kIntegrated = false;
315 m_kIdx = row;
316 m_kmin = std::max(m_kmin, dimMin);
317 m_kmax = std::min(m_kmax, dimMax);
318 if (m_kmin > dimMax || m_kmax < dimMin) {
319 skipNormalization = true;
320 }
321 }
322 if (affineMat[row][2] == 1.) {
323 m_lIntegrated = false;
324 m_lIdx = row;
325 m_lmin = std::max(m_lmin, dimMin);
326 m_lmax = std::min(m_lmax, dimMax);
327 if (m_lmin > dimMax || m_lmax < dimMin) {
328 skipNormalization = true;
329 }
330 }
331
332 for (size_t col = 3; col < ncm1; col++) // affine matrix, ignore last column
333 {
334 if (affineMat[row][col] == 1.) {
335 double val = otherDimValues.at(col - 3);
336 if (val > dimMax || val < dimMin) {
337 skipNormalization = true;
338 }
339 }
340 }
341 }
342
343 return affineMat;
344}
345
350 if (!m_hIntegrated) {
351 auto &hDim = *m_normWS->getDimension(m_hIdx);
352 m_hX.resize(hDim.getNBoundaries());
353 for (size_t i = 0; i < m_hX.size(); ++i) {
354 m_hX[i] = hDim.getX(i);
355 }
356 }
357 if (!m_kIntegrated) {
358 auto &kDim = *m_normWS->getDimension(m_kIdx);
359 m_kX.resize(kDim.getNBoundaries());
360 for (size_t i = 0; i < m_kX.size(); ++i) {
361 m_kX[i] = kDim.getX(i);
362 }
363 }
364 if (!m_lIntegrated) {
365 auto &lDim = *m_normWS->getDimension(m_lIdx);
366 m_lX.resize(lDim.getNBoundaries());
367 for (size_t i = 0; i < m_lX.size(); ++i) {
368 m_lX[i] = lDim.getX(i);
369 }
370 }
371}
372
380void MDNormSCD::calculateNormalization(const std::vector<coord_t> &otherValues,
381 const Kernel::Matrix<coord_t> &affineTrans, uint16_t expInfoIndex) {
382 API::MatrixWorkspace_const_sptr integrFlux = getProperty("FluxWorkspace");
383 integrFlux->getXMinMax(m_kiMin, m_kiMax);
384 API::MatrixWorkspace_const_sptr solidAngleWS = getProperty("SolidAngleWorkspace");
385
386 const auto &currentExptInfo = *(m_inputWS->getExperimentInfo(expInfoIndex));
388 auto *rubwLog = dynamic_cast<VectorDoubleProperty *>(currentExptInfo.getLog("RUBW_MATRIX"));
389 if (!rubwLog) {
390 throw std::runtime_error("Wokspace does not contain a log entry for the RUBW matrix."
391 "Cannot continue.");
392 } else {
393 Kernel::DblMatrix rubwValue((*rubwLog)()); // includes the 2*pi factor but not goniometer for now :)
394 m_rubw = currentExptInfo.run().getGoniometerMatrix() * rubwValue;
395 m_rubw.Invert();
396 }
397 const double protonCharge = currentExptInfo.run().getProtonCharge();
398
399 const auto &spectrumInfo = currentExptInfo.spectrumInfo();
400
401 // Mappings
402 const auto ndets = static_cast<int64_t>(spectrumInfo.size());
403 const detid2index_map fluxDetToIdx = integrFlux->getDetectorIDToWorkspaceIndexMap();
404 const detid2index_map solidAngDetToIdx = solidAngleWS->getDetectorIDToWorkspaceIndexMap();
405
406 const size_t vmdDims = 4;
407 std::vector<std::atomic<signal_t>> signalArray(m_normWS->getNPoints());
408 std::vector<std::array<double, 4>> intersections;
409 std::vector<double> xValues, yValues;
410 std::vector<coord_t> pos, posNew;
411 double progStep = 0.7 / m_numExptInfos;
412 auto prog =
413 std::make_unique<API::Progress>(this, 0.3 + progStep * expInfoIndex, 0.3 + progStep * (expInfoIndex + 1.), ndets);
414
415PRAGMA_OMP(parallel for private(intersections, xValues, yValues, pos, posNew) if (Kernel::threadSafe(*integrFlux)))
416for (int64_t i = 0; i < ndets; i++) {
418
419 if (!spectrumInfo.hasDetectors(i) || spectrumInfo.isMonitor(i) || spectrumInfo.isMasked(i)) {
420 continue;
421 }
422 const auto &detector = spectrumInfo.detector(i);
423 double theta = detector.getTwoTheta(m_samplePos, m_beamDir);
424 double phi = detector.getPhi();
425 // If the dtefctor is a group, this should be the ID of the first detector
426 const auto detID = detector.getID();
427
428 // Intersections
429 this->calculateIntersections(intersections, theta, phi);
430 if (intersections.empty())
431 continue;
432
433 // get the flux spetrum number
434 size_t wsIdx = fluxDetToIdx.find(detID)->second;
435 // Get solid angle for this contribution
436 double solid = solidAngleWS->y(solidAngDetToIdx.find(detID)->second)[0] * protonCharge;
437
438 // -- calculate integrals for the intersection --
439 // momentum values at intersections
440 auto intersectionsBegin = intersections.begin();
441 // copy momenta to xValues
442 xValues.resize(intersections.size());
443 yValues.resize(intersections.size());
444 auto x = xValues.begin();
445 for (auto it = intersectionsBegin; it != intersections.end(); ++it, ++x) {
446 *x = (*it)[3];
447 }
448 // calculate integrals at momenta from xValues by interpolating between
449 // points in spectrum sp
450 // of workspace integrFlux. The result is stored in yValues
451 calcIntegralsForIntersections(xValues, *integrFlux, wsIdx, yValues);
452
453 // Compute final position in HKL
454 // pre-allocate for efficiency and copy non-hkl dim values into place
455 pos.resize(vmdDims + otherValues.size());
456 std::copy(otherValues.begin(), otherValues.end(), pos.begin() + vmdDims - 1);
457 pos.emplace_back(1.f);
458
459 for (auto it = intersectionsBegin + 1; it != intersections.end(); ++it) {
460 const auto &curIntSec = *it;
461 const auto &prevIntSec = *(it - 1);
462 // the full vector isn't used so compute only what is necessary
463 double delta = curIntSec[3] - prevIntSec[3];
464 if (delta < 1e-07)
465 continue; // Assume zero contribution if difference is small
466
467 // Average between two intersections for final position
468 std::transform(curIntSec.begin(), curIntSec.begin() + vmdDims - 1, prevIntSec.begin(), pos.begin(),
469 [](const double lhs, const double rhs) { return static_cast<coord_t>(0.5 * (lhs + rhs)); });
470 affineTrans.multiplyPoint(pos, posNew);
471 size_t linIndex = m_normWS->getLinearIndexAtCoord(posNew.data());
472 if (linIndex == size_t(-1))
473 continue;
474
475 // index of the current intersection
476 auto k = static_cast<size_t>(std::distance(intersectionsBegin, it));
477 // signal = integral between two consecutive intersections
478 signal_t signal = (yValues[k] - yValues[k - 1]) * solid;
479 Mantid::Kernel::AtomicOp(signalArray[linIndex], signal, std::plus<signal_t>());
480 }
481 prog->report();
482
484}
486if (m_accumulate) {
487 std::transform(signalArray.cbegin(), signalArray.cend(), m_normWS->getSignalArray(), m_normWS->mutableSignalArray(),
488 [](const std::atomic<signal_t> &a, const signal_t &b) { return a + b; });
489} else {
490 std::copy(signalArray.cbegin(), signalArray.cend(), m_normWS->mutableSignalArray());
491}
492}
493
502void MDNormSCD::calcIntegralsForIntersections(const std::vector<double> &xValues,
503 const API::MatrixWorkspace &integrFlux, size_t sp,
504 std::vector<double> &yValues) const {
505 assert(xValues.size() == yValues.size());
506
507 // the x-data from the workspace
508 const auto &xData = integrFlux.x(sp);
509 const double xStart = xData.front();
510 const double xEnd = xData.back();
511
512 // the values in integrFlux are expected to be integrals of a non-negative
513 // function
514 // ie they must make a non-decreasing function
515 const auto &yData = integrFlux.y(sp);
516 size_t spSize = yData.size();
517
518 const double yMin = 0.0;
519 const double yMax = yData.back();
520
521 size_t nData = xValues.size();
522 // all integrals below xStart must be 0
523 if (xValues[nData - 1] < xStart) {
524 std::fill(yValues.begin(), yValues.end(), yMin);
525 return;
526 }
527
528 // all integrals above xEnd must be equal tp yMax
529 if (xValues[0] > xEnd) {
530 std::fill(yValues.begin(), yValues.end(), yMax);
531 return;
532 }
533
534 size_t i = 0;
535 // integrals below xStart must be 0
536 while (i < nData - 1 && xValues[i] < xStart) {
537 yValues[i] = yMin;
538 i++;
539 }
540 size_t j = 0;
541 for (; i < nData; i++) {
542 // integrals above xEnd must be equal tp yMax
543 if (j >= spSize - 1) {
544 yValues[i] = yMax;
545 } else {
546 double xi = xValues[i];
547 while (j < spSize - 1 && xi > xData[j])
548 j++;
549 // if x falls onto an interpolation point return the corresponding y
550 if (xi == xData[j]) {
551 yValues[i] = yData[j];
552 } else if (j == spSize - 1) {
553 // if we get above xEnd it's yMax
554 yValues[i] = yMax;
555 } else if (j > 0) {
556 // interpolate between the consecutive points
557 double x0 = xData[j - 1];
558 double x1 = xData[j];
559 double y0 = yData[j - 1];
560 double y1 = yData[j];
561 yValues[i] = y0 + (y1 - y0) * (xi - x0) / (x1 - x0);
562 } else // j == 0
563 {
564 yValues[i] = yMin;
565 }
566 }
567 }
568}
569
578void MDNormSCD::calculateIntersections(std::vector<std::array<double, 4>> &intersections, const double theta,
579 const double phi) {
580 V3D q(-sin(theta) * cos(phi), -sin(theta) * sin(phi), 1. - cos(theta));
581 q = m_rubw * q;
582 if (convention == "Crystallography") {
583 q *= -1;
584 }
585
586 double hStart = q.X() * m_kiMin, hEnd = q.X() * m_kiMax;
587 double kStart = q.Y() * m_kiMin, kEnd = q.Y() * m_kiMax;
588 double lStart = q.Z() * m_kiMin, lEnd = q.Z() * m_kiMax;
589
590 double eps = 1e-7;
591
592 auto hNBins = m_hX.size();
593 auto kNBins = m_kX.size();
594 auto lNBins = m_lX.size();
595 intersections.clear();
596 intersections.reserve(hNBins + kNBins + lNBins + 8);
597
598 // calculate intersections with planes perpendicular to h
599 if (fabs(hStart - hEnd) > eps) {
600 double fmom = (m_kiMax - m_kiMin) / (hEnd - hStart);
601 double fk = (kEnd - kStart) / (hEnd - hStart);
602 double fl = (lEnd - lStart) / (hEnd - hStart);
603 if (!m_hIntegrated) {
604 for (size_t i = 0; i < hNBins; i++) {
605 double hi = m_hX[i];
606 if ((hi >= m_hmin) && (hi <= m_hmax) && ((hStart - hi) * (hEnd - hi) < 0)) {
607 // if hi is between hStart and hEnd, then ki and li will be between
608 // kStart, kEnd and lStart, lEnd and momi will be between m_kiMin and
609 // KnincidemtmMax
610 double ki = fk * (hi - hStart) + kStart;
611 double li = fl * (hi - hStart) + lStart;
612 if ((ki >= m_kmin) && (ki <= m_kmax) && (li >= m_lmin) && (li <= m_lmax)) {
613 double momi = fmom * (hi - hStart) + m_kiMin;
614 intersections.push_back({{hi, ki, li, momi}});
615 }
616 }
617 }
618 }
619
620 double momhMin = fmom * (m_hmin - hStart) + m_kiMin;
621 if ((momhMin > m_kiMin) && (momhMin < m_kiMax)) {
622 // khmin and lhmin
623 double khmin = fk * (m_hmin - hStart) + kStart;
624 double lhmin = fl * (m_hmin - hStart) + lStart;
625 if ((khmin >= m_kmin) && (khmin <= m_kmax) && (lhmin >= m_lmin) && (lhmin <= m_lmax)) {
626 intersections.push_back({{m_hmin, khmin, lhmin, momhMin}});
627 }
628 }
629 double momhMax = fmom * (m_hmax - hStart) + m_kiMin;
630 if ((momhMax > m_kiMin) && (momhMax < m_kiMax)) {
631 // khmax and lhmax
632 double khmax = fk * (m_hmax - hStart) + kStart;
633 double lhmax = fl * (m_hmax - hStart) + lStart;
634 if ((khmax >= m_kmin) && (khmax <= m_kmax) && (lhmax >= m_lmin) && (lhmax <= m_lmax)) {
635 intersections.push_back({{m_hmax, khmax, lhmax, momhMax}});
636 }
637 }
638 }
639
640 // calculate intersections with planes perpendicular to k
641 if (fabs(kStart - kEnd) > eps) {
642 double fmom = (m_kiMax - m_kiMin) / (kEnd - kStart);
643 double fh = (hEnd - hStart) / (kEnd - kStart);
644 double fl = (lEnd - lStart) / (kEnd - kStart);
645 if (!m_kIntegrated) {
646 for (size_t i = 0; i < kNBins; i++) {
647 double ki = m_kX[i];
648 if ((ki >= m_kmin) && (ki <= m_kmax) && ((kStart - ki) * (kEnd - ki) < 0)) {
649 // if ki is between kStart and kEnd, then hi and li will be between
650 // hStart, hEnd and lStart, lEnd
651 double hi = fh * (ki - kStart) + hStart;
652 double li = fl * (ki - kStart) + lStart;
653 if ((hi >= m_hmin) && (hi <= m_hmax) && (li >= m_lmin) && (li <= m_lmax)) {
654 double momi = fmom * (ki - kStart) + m_kiMin;
655 intersections.push_back({{hi, ki, li, momi}});
656 }
657 }
658 }
659 }
660
661 double momkMin = fmom * (m_kmin - kStart) + m_kiMin;
662 if ((momkMin > m_kiMin) && (momkMin < m_kiMax)) {
663 // hkmin and lkmin
664 double hkmin = fh * (m_kmin - kStart) + hStart;
665 double lkmin = fl * (m_kmin - kStart) + lStart;
666 if ((hkmin >= m_hmin) && (hkmin <= m_hmax) && (lkmin >= m_lmin) && (lkmin <= m_lmax)) {
667 intersections.push_back({{hkmin, m_kmin, lkmin, momkMin}});
668 }
669 }
670 double momkMax = fmom * (m_kmax - kStart) + m_kiMin;
671 if ((momkMax > m_kiMin) && (momkMax < m_kiMax)) {
672 // hkmax and lkmax
673 double hkmax = fh * (m_kmax - kStart) + hStart;
674 double lkmax = fl * (m_kmax - kStart) + lStart;
675 if ((hkmax >= m_hmin) && (hkmax <= m_hmax) && (lkmax >= m_lmin) && (lkmax <= m_lmax)) {
676 intersections.push_back({{hkmax, m_kmax, lkmax, momkMax}});
677 }
678 }
679 }
680
681 // calculate intersections with planes perpendicular to l
682 if (fabs(lStart - lEnd) > eps) {
683 double fmom = (m_kiMax - m_kiMin) / (lEnd - lStart);
684 double fh = (hEnd - hStart) / (lEnd - lStart);
685 double fk = (kEnd - kStart) / (lEnd - lStart);
686 if (!m_lIntegrated) {
687 for (size_t i = 0; i < lNBins; i++) {
688 double li = m_lX[i];
689 if ((li >= m_lmin) && (li <= m_lmax) && ((lStart - li) * (lEnd - li) < 0)) {
690 // if li is between lStart and lEnd, then hi and ki will be between
691 // hStart, hEnd and kStart, kEnd
692 double hi = fh * (li - lStart) + hStart;
693 double ki = fk * (li - lStart) + kStart;
694 if ((hi >= m_hmin) && (hi <= m_hmax) && (ki >= m_kmin) && (ki <= m_kmax)) {
695 double momi = fmom * (li - lStart) + m_kiMin;
696 intersections.push_back({{hi, ki, li, momi}});
697 }
698 }
699 }
700 }
701
702 double momlMin = fmom * (m_lmin - lStart) + m_kiMin;
703 if ((momlMin > m_kiMin) && (momlMin < m_kiMax)) {
704 // hlmin and klmin
705 double hlmin = fh * (m_lmin - lStart) + hStart;
706 double klmin = fk * (m_lmin - lStart) + kStart;
707 if ((hlmin >= m_hmin) && (hlmin <= m_hmax) && (klmin >= m_kmin) && (klmin <= m_kmax)) {
708 intersections.push_back({{hlmin, klmin, m_lmin, momlMin}});
709 }
710 }
711 double momlMax = fmom * (m_lmax - lStart) + m_kiMin;
712 if ((momlMax > m_kiMin) && (momlMax < m_kiMax)) {
713 // khmax and lhmax
714 double hlmax = fh * (m_lmax - lStart) + hStart;
715 double klmax = fk * (m_lmax - lStart) + kStart;
716 if ((hlmax >= m_hmin) && (hlmax <= m_hmax) && (klmax >= m_kmin) && (klmax <= m_kmax)) {
717 intersections.push_back({{hlmax, klmax, m_lmax, momlMax}});
718 }
719 }
720 }
721
722 // add endpoints
723 if ((hStart >= m_hmin) && (hStart <= m_hmax) && (kStart >= m_kmin) && (kStart <= m_kmax) && (lStart >= m_lmin) &&
724 (lStart <= m_lmax)) {
725 intersections.push_back({{hStart, kStart, lStart, m_kiMin}});
726 }
727 if ((hEnd >= m_hmin) && (hEnd <= m_hmax) && (kEnd >= m_kmin) && (kEnd <= m_kmax) && (lEnd >= m_lmin) &&
728 (lEnd <= m_lmax)) {
729 intersections.push_back({{hEnd, kEnd, lEnd, m_kiMax}});
730 }
731
732 // sort intersections by momentum
733 std::stable_sort(intersections.begin(), intersections.end(), compareMomentum);
734}
735
736} // namespace Mantid::MDAlgorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
gsl_vector * tmp
const std::vector< double > & rhs
double value
The value of the point.
Definition: FitMW.cpp:51
#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 PRAGMA_OMP(expression)
#define PARALLEL_CHECK_INTERRUPT_REGION
Adds a check after a Parallel region to see if it was interupted.
std::vector< history_type > history
history information
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
const std::vector< Kernel::Property * > & getProperties() const override
Get the list of managed properties.
Definition: Algorithm.cpp:2050
A validator which provides a TENTATIVE check that a workspace contains common bins in each spectrum.
Kernel::IValidator_sptr clone() const override
Clone the current state.
A validator which checks that a workspace has a valid instrument.
Base MatrixWorkspace Abstract Class.
const HistogramData::HistogramX & x(const size_t index) const
const HistogramData::HistogramY & y(const size_t index) const
A property class for workspaces.
A validator which checks that the unit of the workspace referred to by a WorkspaceProperty is the exp...
std::unique_ptr< MDHistoWorkspace > clone() const
Returns a clone of the workspace.
Exception for errors associated with the instrument definition.
Definition: Exception.h:220
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
Numerical Matrix class.
Definition: Matrix.h:42
T Invert()
LU inversion routine.
Definition: Matrix.cpp:924
void multiplyPoint(const std::vector< T > &in, std::vector< T > &out) const
Multiply M*Vec.
Definition: Matrix.cpp:375
size_t numRows() const
Return the number of rows in the matrix.
Definition: Matrix.h:144
size_t numCols() const
Return the number of columns in the matrix.
Definition: Matrix.h:147
The concrete, templated class for properties.
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
A specialised Property class for holding a series of time-value pairs.
Class for 3D vectors.
Definition: V3D.h:34
constexpr double X() const noexcept
Get x.
Definition: V3D.h:232
constexpr double Y() const noexcept
Get y.
Definition: V3D.h:233
constexpr double Z() const noexcept
Get z.
Definition: V3D.h:234
MDNormSCD : Generate MD normalization for single crystal diffraction.
Definition: MDNormSCD.h:20
const std::string category() const override
Algorithm's category for identification.
Definition: MDNormSCD.cpp:53
bool m_accumulate
internal flag to accumulate to an existing workspace
Definition: MDNormSCD.h:74
DataObjects::MDHistoWorkspace_sptr m_normWS
Normalization workspace.
Definition: MDNormSCD.h:52
std::vector< double > m_kX
Definition: MDNormSCD.h:66
bool m_hIntegrated
flag for integrated h,k,l dimensions
Definition: MDNormSCD.h:58
void cacheDimensionXValues()
Stores the X values from each H,K,L dimension as member variables.
Definition: MDNormSCD.cpp:349
coord_t m_hmin
limits for h,k,l dimensions
Definition: MDNormSCD.h:56
void exec() override
Execute the algorithm.
Definition: MDNormSCD.cpp:125
std::string inputEnergyMode() const
Currently looks for the ConvertToMD algorithm in the history.
Definition: MDNormSCD.cpp:190
Mantid::Kernel::DblMatrix m_rubw
(2*PiRUBW)^-1
Definition: MDNormSCD.h:60
DataObjects::MDHistoWorkspace_sptr binInputWS()
Runs the BinMD algorithm on the input to provide the output workspace All slicing algorithm propertie...
Definition: MDNormSCD.cpp:220
std::vector< double > m_hX
cached X values along dimensions h,k,l
Definition: MDNormSCD.h:66
const std::string name() const override
Algorithm's name for use in the GUI and help.
Definition: MDNormSCD.cpp:62
int version() const override
Algorithm's version for identification.
Definition: MDNormSCD.cpp:50
void calculateNormalization(const std::vector< coord_t > &otherValues, const Kernel::Matrix< coord_t > &affineTrans, uint16_t expInfoIndex)
Computed the normalization for the input workspace.
Definition: MDNormSCD.cpp:380
double m_kiMin
limits for momentum
Definition: MDNormSCD.h:62
std::vector< double > m_lX
Definition: MDNormSCD.h:66
uint16_t m_numExptInfos
number of experiment infos
Definition: MDNormSCD.h:76
size_t m_hIdx
index of h,k,l dimensions in the output workspaces
Definition: MDNormSCD.h:64
Kernel::V3D m_samplePos
Sample position.
Definition: MDNormSCD.h:68
void init() override
Initialize the algorithm's properties.
Definition: MDNormSCD.cpp:67
Kernel::Matrix< coord_t > findIntergratedDimensions(const std::vector< coord_t > &otherDimValues, bool &skipNormalization)
Checks the normalization workspace against the indices of the original dimensions.
Definition: MDNormSCD.cpp:292
API::IMDEventWorkspace_sptr m_inputWS
Input workspace.
Definition: MDNormSCD.h:54
void calculateIntersections(std::vector< std::array< double, 4 > > &intersections, const double theta, const double phi)
Calculate the points of intersection for the given detector with cuboid surrounding the detector posi...
Definition: MDNormSCD.cpp:578
const std::string summary() const override
Algorithm's summary for use in the GUI and help.
Definition: MDNormSCD.cpp:56
std::vector< coord_t > getValuesFromOtherDimensions(bool &skipNormalization, uint16_t expInfoIndex=0) const
Retrieve logged values from non-HKL dimensions.
Definition: MDNormSCD.cpp:261
Kernel::V3D m_beamDir
Beam direction.
Definition: MDNormSCD.h:70
void createNormalizationWS(const DataObjects::MDHistoWorkspace &dataWS)
Create & cached the normalization workspace.
Definition: MDNormSCD.cpp:241
void cacheInputs()
Set up starting values for cached variables.
Definition: MDNormSCD.cpp:158
std::string convention
ki-kf for Inelastic convention; kf-ki for Crystallography convention
Definition: MDNormSCD.h:72
void calcIntegralsForIntersections(const std::vector< double > &xValues, const API::MatrixWorkspace &integrFlux, size_t sp, std::vector< double > &yValues) const
Linearly interpolate between the points in integrFlux at xValues and save the results in yValues.
Definition: MDNormSCD.cpp:502
std::shared_ptr< Workspace > Workspace_sptr
shared pointer to Mantid::API::Workspace
Definition: Workspace_fwd.h:20
std::shared_ptr< const MatrixWorkspace > MatrixWorkspace_const_sptr
shared pointer to the matrix workspace base class (const version)
@ NoNormalization
Don't normalize = return raw counts.
Definition: IMDIterator.h:27
std::shared_ptr< MDHistoWorkspace > MDHistoWorkspace_sptr
A shared pointer to a MDHistoWorkspace.
std::string toString(const T &value)
Convert a number to a string.
Definition: Strings.cpp:703
std::vector< PropertyHistory_sptr > PropertyHistories
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
void AtomicOp(std::atomic< T > &f, T d, BinaryOp op)
Uses std::compare_exchange_weak to update the atomic value f = op(f, d) Used to improve parallel scal...
Definition: MultiThreaded.h:67
float coord_t
Typedef for the data type to use for coordinate axes in MD objects such as MDBox, MDEventWorkspace,...
Definition: MDTypes.h:27
std::unordered_map< detid_t, size_t > detid2index_map
Map with key = detector ID, value = workspace index.
double signal_t
Typedef for the signal recorded in a MDBox, etc.
Definition: MDTypes.h:36
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
Definition: EmptyValues.h:43
Describes the direction (within an algorithm) of a Property.
Definition: Property.h:50
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54