Mantid
Loading...
Searching...
No Matches
MDNormDirectSC.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(MDNormDirectSC)
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_dEmin(0.f), m_dEmax(0.f), m_Ei(0.), m_ki(0.), m_kfmin(0.), m_kfmax(0.), m_hIntegrated(true),
47 m_kIntegrated(true), m_lIntegrated(true), m_dEIntegrated(true), m_rubw(3, 3), m_hIdx(-1), m_kIdx(-1), m_lIdx(-1),
48 m_eIdx(-1), m_hX(), m_kX(), m_lX(), m_eX(), m_samplePos(), m_beamDir() {}
49
51int MDNormDirectSC::version() const { return 1; }
52
54const std::string MDNormDirectSC::category() const { return "MDAlgorithms\\Normalisation"; }
55
57const std::string MDNormDirectSC::summary() const {
58 return "Calculate normalization for an MDEvent workspace for single crystal "
59 "direct geometry inelastic measurement.";
60}
61
63const std::string MDNormDirectSC::name() const { return "MDNormDirectSC"; }
64
69 declareProperty(std::make_unique<WorkspaceProperty<IMDEventWorkspace>>("InputWorkspace", "", Direction::Input),
70 "An input MDWorkspace.");
71
72 std::string dimChars = getDimensionChars();
73 // --------------- Axis-aligned properties
74 // ---------------------------------------
75 for (size_t i = 0; i < dimChars.size(); i++) {
76 std::string dim(" ");
77 dim[0] = dimChars[i];
78 std::string propName = "AlignedDim" + dim;
80 "Binning parameters for the " + Strings::toString(i) +
81 "th dimension.\n"
82 "Enter it as a comma-separated list of values with the format: "
83 "'name,minimum,maximum,number_of_bins'. Leave blank for NONE.");
84 }
85
86 auto solidAngleValidator = std::make_shared<CompositeValidator>();
87 solidAngleValidator->add<InstrumentValidator>();
88 solidAngleValidator->add<CommonBinsValidator>();
89
90 declareProperty(std::make_unique<WorkspaceProperty<>>("SolidAngleWorkspace", "", Direction::Input,
91 PropertyMode::Optional, solidAngleValidator),
92 "An input workspace containing integrated vanadium (a measure of the "
93 "solid angle).");
94
95 declareProperty(std::make_unique<PropertyWithValue<bool>>("SkipSafetyCheck", false, Direction::Input),
96 "If set to true, the algorithm does "
97 "not check history if the workspace was modified since the"
98 "ConvertToMD algorithm was run, and assume that the direct "
99 "geometry inelastic mode is used.");
100
101 declareProperty(std::make_unique<WorkspaceProperty<IMDHistoWorkspace>>("TemporaryNormalizationWorkspace", "",
103 "An input MDHistoWorkspace used to accumulate normalization "
104 "from multiple MDEventWorkspaces. If unspecified a blank "
105 "MDHistoWorkspace will be created.");
106
107 declareProperty(std::make_unique<WorkspaceProperty<IMDHistoWorkspace>>("TemporaryDataWorkspace", "", Direction::Input,
109 "An input MDHistoWorkspace used to accumulate data from "
110 "multiple MDEventWorkspaces. If unspecified a blank "
111 "MDHistoWorkspace will be created.");
112
113 declareProperty(std::make_unique<WorkspaceProperty<Workspace>>("OutputWorkspace", "", Direction::Output),
114 "A name for the output data MDHistoWorkspace.");
115 declareProperty(std::make_unique<WorkspaceProperty<Workspace>>("OutputNormalizationWorkspace", "", Direction::Output),
116 "A name for the output normalization MDHistoWorkspace.");
117}
118
119//----------------------------------------------------------------------------------------------
124 cacheInputs();
125 auto outputWS = binInputWS();
126 convention = Kernel::ConfigService::Instance().getString("Q.convention");
127 outputWS->setDisplayNormalization(Mantid::API::NoNormalization);
128 setProperty<Workspace_sptr>("OutputWorkspace", outputWS);
129 createNormalizationWS(*outputWS);
130 m_normWS->setDisplayNormalization(Mantid::API::NoNormalization);
131 setProperty("OutputNormalizationWorkspace", m_normWS);
132
133 m_numExptInfos = outputWS->getNumExperimentInfo();
134 // loop over all experiment infos
135 for (uint16_t expInfoIndex = 0; expInfoIndex < m_numExptInfos; expInfoIndex++) {
136 // Check for other dimensions if we could measure anything in the original
137 // data
138 bool skipNormalization = false;
139 const std::vector<coord_t> otherValues = getValuesFromOtherDimensions(skipNormalization, expInfoIndex);
140 const auto affineTrans = findIntergratedDimensions(otherValues, skipNormalization);
142
143 if (!skipNormalization) {
144 calculateNormalization(otherValues, affineTrans, expInfoIndex);
145 } else {
146 g_log.warning("Binning limits are outside the limits of the MDWorkspace. "
147 "Not applying normalization.");
148 }
149 // if more than one experiment info, keep accumulating
150 m_accumulate = true;
151 }
152
153 // Set the display normalization based on the input workspace
154 outputWS->setDisplayNormalization(m_inputWS->displayNormalizationHisto());
155}
156
161 m_inputWS = getProperty("InputWorkspace");
162 bool skipCheck = getProperty("SkipSafetyCheck");
163 if (!skipCheck && (inputEnergyMode() != "Direct")) {
164 throw std::invalid_argument("Invalid energy transfer mode. Algorithm only "
165 "supports direct geometry spectrometers.");
166 }
167 // Min/max dimension values
168 const auto hdim(m_inputWS->getDimension(0)), kdim(m_inputWS->getDimension(1)), ldim(m_inputWS->getDimension(2)),
169 edim(m_inputWS->getDimension(3));
170 m_hmin = hdim->getMinimum();
171 m_kmin = kdim->getMinimum();
172 m_lmin = ldim->getMinimum();
173 m_dEmin = edim->getMinimum();
174 m_hmax = hdim->getMaximum();
175 m_kmax = kdim->getMaximum();
176 m_lmax = ldim->getMaximum();
177 m_dEmax = edim->getMaximum();
178
179 const auto &exptInfoZero = *(m_inputWS->getExperimentInfo(0));
180 auto source = exptInfoZero.getInstrument()->getSource();
181 auto sample = exptInfoZero.getInstrument()->getSample();
182 if (source == nullptr || sample == nullptr) {
184 "Instrument not sufficiently defined: failed to get source and/or "
185 "sample");
186 }
187 m_samplePos = sample->getPos();
188 m_beamDir = normalize(m_samplePos - source->getPos());
189
190 double originaldEmin = exptInfoZero.run().getBinBoundaries().front();
191 double originaldEmax = exptInfoZero.run().getBinBoundaries().back();
192 if (exptInfoZero.run().hasProperty("Ei")) {
193 m_Ei = exptInfoZero.run().getPropertyValueAsType<double>("Ei");
194 if (m_Ei <= 0) {
195 throw std::invalid_argument("Ei stored in the workspace is not positive");
196 }
197 } else {
198 throw std::invalid_argument("Could not find Ei value in the workspace.");
199 }
200 double eps = 1e-7;
201 if (m_Ei - originaldEmin < eps) {
202 originaldEmin = m_Ei - eps;
203 }
204 if (m_Ei - originaldEmax < eps) {
205 originaldEmax = m_Ei - 1e-7;
206 }
207 if (originaldEmin == originaldEmax) {
208 throw std::runtime_error("The limits of the original workspace used in "
209 "ConvertToMD are incorrect");
210 }
211 const double energyToK = 8.0 * M_PI * M_PI * PhysicalConstants::NeutronMass * PhysicalConstants::meV * 1e-20 /
213 m_ki = std::sqrt(energyToK * m_Ei);
214 m_kfmin = std::sqrt(energyToK * (m_Ei - originaldEmin));
215 m_kfmax = std::sqrt(energyToK * (m_Ei - originaldEmax));
216}
217
223 const auto &hist = m_inputWS->getHistory();
224 const size_t nalgs = hist.size();
225 const auto &lastAlgHist = hist.getAlgorithmHistory(nalgs - 1);
226 const auto &penultimateAlgHist = hist.getAlgorithmHistory(nalgs - 2);
227
228 std::string emode;
229 if (lastAlgHist->name() == "ConvertToMD") {
230 emode = lastAlgHist->getPropertyValue("dEAnalysisMode");
231 } else if ((lastAlgHist->name() == "Load" || lastAlgHist->name() == "LoadMD") &&
232 penultimateAlgHist->name() == "ConvertToMD") {
233 // get dEAnalysisMode
234 emode = penultimateAlgHist->getPropertyValue("dEAnalysisMode");
235 } else {
236 throw std::invalid_argument("The last algorithm in the history of the "
237 "input workspace is not ConvertToMD");
238 }
239 return emode;
240}
241
248 const auto &props = getProperties();
249 auto binMD = createChildAlgorithm("BinMD", 0.0, 0.3);
250 binMD->setPropertyValue("AxisAligned", "1");
251 for (auto prop : props) {
252 const auto &propName = prop->name();
253 if (propName != "SolidAngleWorkspace" && propName != "TemporaryNormalizationWorkspace" &&
254 propName != "OutputNormalizationWorkspace" && propName != "SkipSafetyCheck") {
255 binMD->setPropertyValue(propName, prop->value());
256 }
257 }
258 binMD->executeAsChildAlg();
259 Workspace_sptr outputWS = binMD->getProperty("OutputWorkspace");
260 return std::dynamic_pointer_cast<MDHistoWorkspace>(outputWS);
261}
262
268 // Copy the MDHisto workspace, and change signals and errors to 0.
269 std::shared_ptr<IMDHistoWorkspace> tmp = this->getProperty("TemporaryNormalizationWorkspace");
270 m_normWS = std::dynamic_pointer_cast<MDHistoWorkspace>(tmp);
271 if (!m_normWS) {
272 m_normWS = dataWS.clone();
273 m_normWS->setTo(0., 0., 0.);
274 } else {
275 m_accumulate = true;
276 }
277}
278
287std::vector<coord_t> MDNormDirectSC::getValuesFromOtherDimensions(bool &skipNormalization,
288 uint16_t expInfoIndex) const {
289 const auto &currentRun = m_inputWS->getExperimentInfo(expInfoIndex)->run();
290
291 std::vector<coord_t> otherDimValues;
292 for (size_t i = 4; i < m_inputWS->getNumDims(); i++) {
293 const auto dimension = m_inputWS->getDimension(i);
294 auto dimMin = static_cast<float>(dimension->getMinimum());
295 auto dimMax = static_cast<float>(dimension->getMaximum());
296 auto *dimProp = dynamic_cast<Kernel::TimeSeriesProperty<double> *>(currentRun.getProperty(dimension->getName()));
297 if (dimProp) {
298 auto value = static_cast<coord_t>(dimProp->firstValue());
299 otherDimValues.emplace_back(value);
300 // in the original MD data no time was spent measuring between dimMin and
301 // dimMax
302 if (value < dimMin || value > dimMax) {
303 skipNormalization = true;
304 }
305 }
306 }
307 return otherDimValues;
308}
309
320 bool &skipNormalization) {
321 // Get indices of the original dimensions in the output workspace,
322 // and if not found, the corresponding dimension is integrated
323 Kernel::Matrix<coord_t> affineMat = m_normWS->getTransformFromOriginal(0)->makeAffineMatrix();
324
325 const size_t nrm1 = affineMat.numRows() - 1;
326 const size_t ncm1 = affineMat.numCols() - 1;
327 for (size_t row = 0; row < nrm1; row++) // affine matrix, ignore last row
328 {
329 const auto dimen = m_normWS->getDimension(row);
330 const auto dimMin(dimen->getMinimum()), dimMax(dimen->getMaximum());
331 if (affineMat[row][0] == 1.) {
332 m_hIntegrated = false;
333 m_hIdx = row;
334 m_hmin = std::max(m_hmin, dimMin);
335 m_hmax = std::min(m_hmax, dimMax);
336 if (m_hmin > dimMax || m_hmax < dimMin) {
337 skipNormalization = true;
338 }
339 }
340 if (affineMat[row][1] == 1.) {
341 m_kIntegrated = false;
342 m_kIdx = row;
343 m_kmin = std::max(m_kmin, dimMin);
344 m_kmax = std::min(m_kmax, dimMax);
345 if (m_kmin > dimMax || m_kmax < dimMin) {
346 skipNormalization = true;
347 }
348 }
349 if (affineMat[row][2] == 1.) {
350 m_lIntegrated = false;
351 m_lIdx = row;
352 m_lmin = std::max(m_lmin, dimMin);
353 m_lmax = std::min(m_lmax, dimMax);
354 if (m_lmin > dimMax || m_lmax < dimMin) {
355 skipNormalization = true;
356 }
357 }
358
359 if (affineMat[row][3] == 1.) {
360 m_dEIntegrated = false;
361 m_eIdx = row;
362 m_dEmin = std::max(m_dEmin, dimMin);
363 m_dEmax = std::min(m_dEmax, dimMax);
364 if (m_dEmin > dimMax || m_dEmax < dimMin) {
365 skipNormalization = true;
366 }
367 }
368 for (size_t col = 4; col < ncm1; col++) // affine matrix, ignore last column
369 {
370 if (affineMat[row][col] == 1.) {
371 double val = otherDimValues.at(col - 3);
372 if (val > dimMax || val < dimMin) {
373 skipNormalization = true;
374 }
375 }
376 }
377 }
378
379 return affineMat;
380}
381
387 constexpr double energyToK = 8.0 * M_PI * M_PI * PhysicalConstants::NeutronMass * PhysicalConstants::meV * 1e-20 /
389 if (!m_hIntegrated) {
390 auto &hDim = *m_normWS->getDimension(m_hIdx);
391 m_hX.resize(hDim.getNBoundaries());
392 for (size_t i = 0; i < m_hX.size(); ++i) {
393 m_hX[i] = hDim.getX(i);
394 }
395 }
396 if (!m_kIntegrated) {
397 auto &kDim = *m_normWS->getDimension(m_kIdx);
398 m_kX.resize(kDim.getNBoundaries());
399 for (size_t i = 0; i < m_kX.size(); ++i) {
400 m_kX[i] = kDim.getX(i);
401 }
402 }
403 if (!m_lIntegrated) {
404 auto &lDim = *m_normWS->getDimension(m_lIdx);
405 m_lX.resize(lDim.getNBoundaries());
406 for (size_t i = 0; i < m_lX.size(); ++i) {
407 m_lX[i] = lDim.getX(i);
408 }
409 }
410 if (!m_dEIntegrated) {
411 // NOTE: store k final instead
412 auto &eDim = *m_normWS->getDimension(m_eIdx);
413 m_eX.resize(eDim.getNBoundaries());
414 for (size_t i = 0; i < m_eX.size(); ++i) {
415 double temp = m_Ei - eDim.getX(i);
416 temp = std::max(temp, 0.);
417 m_eX[i] = std::sqrt(energyToK * temp);
418 }
419 }
420}
421
429void MDNormDirectSC::calculateNormalization(const std::vector<coord_t> &otherValues,
430 const Kernel::Matrix<coord_t> &affineTrans, uint16_t expInfoIndex) {
431 constexpr double energyToK = 8.0 * M_PI * M_PI * PhysicalConstants::NeutronMass * PhysicalConstants::meV * 1e-20 /
433 const auto &currentExptInfo = *(m_inputWS->getExperimentInfo(expInfoIndex));
435 auto *rubwLog = dynamic_cast<VectorDoubleProperty *>(currentExptInfo.getLog("RUBW_MATRIX"));
436 if (!rubwLog) {
437 throw std::runtime_error("Wokspace does not contain a log entry for the RUBW matrix."
438 "Cannot continue.");
439 } else {
440 Kernel::DblMatrix rubwValue((*rubwLog)()); // includes the 2*pi factor but not goniometer for now :)
441 m_rubw = currentExptInfo.run().getGoniometerMatrix() * rubwValue;
442 m_rubw.Invert();
443 }
444 const double protonCharge = currentExptInfo.run().getProtonCharge();
445
446 const auto &spectrumInfo = currentExptInfo.spectrumInfo();
447
448 // Mapping
449 const auto ndets = static_cast<int64_t>(spectrumInfo.size());
450 bool haveSA = false;
451 API::MatrixWorkspace_const_sptr solidAngleWS = getProperty("SolidAngleWorkspace");
452 detid2index_map solidAngDetToIdx;
453 if (solidAngleWS != nullptr) {
454 haveSA = true;
455 solidAngDetToIdx = solidAngleWS->getDetectorIDToWorkspaceIndexMap();
456 }
457
458 const size_t vmdDims = 4;
459 std::vector<std::atomic<signal_t>> signalArray(m_normWS->getNPoints());
460 std::vector<std::array<double, 4>> intersections;
461 std::vector<coord_t> pos, posNew;
462 double progStep = 0.7 / m_numExptInfos;
463 auto prog =
464 std::make_unique<API::Progress>(this, 0.3 + progStep * expInfoIndex, 0.3 + progStep * (expInfoIndex + 1.), ndets);
465
466PRAGMA_OMP(parallel for private(intersections, pos, posNew))
467for (int64_t i = 0; i < ndets; i++) {
469
470 if (!spectrumInfo.hasDetectors(i) || spectrumInfo.isMonitor(i) || spectrumInfo.isMasked(i)) {
471 continue;
472 }
473 const auto &detector = spectrumInfo.detector(i);
474 double theta = detector.getTwoTheta(m_samplePos, m_beamDir);
475 double phi = detector.getPhi();
476 // If the detector is a group, this should be the ID of the first detector
477 const auto detID = detector.getID();
478
479 // Intersections
480 this->calculateIntersections(intersections, theta, phi);
481 if (intersections.empty())
482 continue;
483
484 // Get solid angle for this contribution
485 double solid = protonCharge;
486 if (haveSA) {
487 solid = solidAngleWS->y(solidAngDetToIdx.find(detID)->second)[0] * protonCharge;
488 }
489 // Compute final position in HKL
490 // pre-allocate for efficiency and copy non-hkl dim values into place
491 pos.resize(vmdDims + otherValues.size() + 1);
492 std::copy(otherValues.begin(), otherValues.end(), pos.begin() + vmdDims);
493 pos.emplace_back(1.f);
494 auto intersectionsBegin = intersections.begin();
495 for (auto it = intersectionsBegin + 1; it != intersections.end(); ++it) {
496 const auto &curIntSec = *it;
497 const auto &prevIntSec = *(it - 1);
498 // the full vector isn't used so compute only what is necessary
499 double delta = (curIntSec[3] * curIntSec[3] - prevIntSec[3] * prevIntSec[3]) / energyToK;
500 if (delta < 1e-10)
501 continue; // Assume zero contribution if difference is small
502
503 // Average between two intersections for final position
504 std::transform(curIntSec.data(), curIntSec.data() + vmdDims, prevIntSec.data(), pos.begin(),
505 [](const double rhs, const double lhs) { return static_cast<coord_t>(0.5 * (rhs + lhs)); });
506
507 // transform kf to energy transfer
508 pos[3] = static_cast<coord_t>(m_Ei - pos[3] * pos[3] / energyToK);
509 affineTrans.multiplyPoint(pos, posNew);
510 size_t linIndex = m_normWS->getLinearIndexAtCoord(posNew.data());
511 if (linIndex == size_t(-1))
512 continue;
513
514 // signal = integral between two consecutive intersections *solid angle
515 // *PC
516 double signal = solid * delta;
517 Mantid::Kernel::AtomicOp(signalArray[linIndex], signal, std::plus<signal_t>());
518 }
519 prog->report();
520
522}
524if (m_accumulate) {
525 std::transform(signalArray.cbegin(), signalArray.cend(), m_normWS->getSignalArray(), m_normWS->mutableSignalArray(),
526 [](const std::atomic<signal_t> &a, const signal_t &b) { return a + b; });
527} else {
528 std::copy(signalArray.cbegin(), signalArray.cend(), m_normWS->mutableSignalArray());
529}
530}
531
540void MDNormDirectSC::calculateIntersections(std::vector<std::array<double, 4>> &intersections, const double theta,
541 const double phi) {
542 V3D qout(sin(theta) * cos(phi), sin(theta) * sin(phi), cos(theta)), qin(0., 0., m_ki);
543
544 qout = m_rubw * qout;
545 qin = m_rubw * qin;
546 if (convention == "Crystallography") {
547 qout *= -1;
548 qin *= -1;
549 }
550 double hStart = qin.X() - qout.X() * m_kfmin, hEnd = qin.X() - qout.X() * m_kfmax;
551 double kStart = qin.Y() - qout.Y() * m_kfmin, kEnd = qin.Y() - qout.Y() * m_kfmax;
552 double lStart = qin.Z() - qout.Z() * m_kfmin, lEnd = qin.Z() - qout.Z() * m_kfmax;
553 double eps = 1e-10;
554 auto hNBins = m_hX.size();
555 auto kNBins = m_kX.size();
556 auto lNBins = m_lX.size();
557 auto eNBins = m_eX.size();
558 intersections.clear();
559 intersections.reserve(hNBins + kNBins + lNBins + eNBins + 8); // 8 is 3*(min,max for each Q component)+kfmin+kfmax
560
561 // calculate intersections with planes perpendicular to h
562 if (fabs(hStart - hEnd) > eps) {
563 double fmom = (m_kfmax - m_kfmin) / (hEnd - hStart);
564 double fk = (kEnd - kStart) / (hEnd - hStart);
565 double fl = (lEnd - lStart) / (hEnd - hStart);
566 if (!m_hIntegrated) {
567 for (size_t i = 0; i < hNBins; i++) {
568 double hi = m_hX[i];
569 if ((hi >= m_hmin) && (hi <= m_hmax) && ((hStart - hi) * (hEnd - hi) < 0)) {
570 // if hi is between hStart and hEnd, then ki and li will be between
571 // kStart, kEnd and lStart, lEnd and momi will be between m_kfmin and
572 // m_kfmax
573 double ki = fk * (hi - hStart) + kStart;
574 double li = fl * (hi - hStart) + lStart;
575 if ((ki >= m_kmin) && (ki <= m_kmax) && (li >= m_lmin) && (li <= m_lmax)) {
576 double momi = fmom * (hi - hStart) + m_kfmin;
577 intersections.push_back({{hi, ki, li, momi}});
578 }
579 }
580 }
581 }
582 double momhMin = fmom * (m_hmin - hStart) + m_kfmin;
583 if ((momhMin - m_kfmin) * (momhMin - m_kfmax) < 0) // m_kfmin>m_kfmax
584 {
585 // khmin and lhmin
586 double khmin = fk * (m_hmin - hStart) + kStart;
587 double lhmin = fl * (m_hmin - hStart) + lStart;
588 if ((khmin >= m_kmin) && (khmin <= m_kmax) && (lhmin >= m_lmin) && (lhmin <= m_lmax)) {
589 intersections.push_back({{m_hmin, khmin, lhmin, momhMin}});
590 }
591 }
592 double momhMax = fmom * (m_hmax - hStart) + m_kfmin;
593 if ((momhMax - m_kfmin) * (momhMax - m_kfmax) <= 0) {
594 // khmax and lhmax
595 double khmax = fk * (m_hmax - hStart) + kStart;
596 double lhmax = fl * (m_hmax - hStart) + lStart;
597 if ((khmax >= m_kmin) && (khmax <= m_kmax) && (lhmax >= m_lmin) && (lhmax <= m_lmax)) {
598 intersections.push_back({{m_hmax, khmax, lhmax, momhMax}});
599 }
600 }
601 }
602
603 // calculate intersections with planes perpendicular to k
604 if (fabs(kStart - kEnd) > eps) {
605 double fmom = (m_kfmax - m_kfmin) / (kEnd - kStart);
606 double fh = (hEnd - hStart) / (kEnd - kStart);
607 double fl = (lEnd - lStart) / (kEnd - kStart);
608 if (!m_kIntegrated) {
609 for (size_t i = 0; i < kNBins; i++) {
610 double ki = m_kX[i];
611 if ((ki >= m_kmin) && (ki <= m_kmax) && ((kStart - ki) * (kEnd - ki) < 0)) {
612 // if ki is between kStart and kEnd, then hi and li will be between
613 // hStart, hEnd and lStart, lEnd and momi will be between m_kfmin and
614 // m_kfmax
615 double hi = fh * (ki - kStart) + hStart;
616 double li = fl * (ki - kStart) + lStart;
617 if ((hi >= m_hmin) && (hi <= m_hmax) && (li >= m_lmin) && (li <= m_lmax)) {
618 double momi = fmom * (ki - kStart) + m_kfmin;
619 intersections.push_back({{hi, ki, li, momi}});
620 }
621 }
622 }
623 }
624 double momkMin = fmom * (m_kmin - kStart) + m_kfmin;
625 if ((momkMin - m_kfmin) * (momkMin - m_kfmax) < 0) {
626 // hkmin and lkmin
627 double hkmin = fh * (m_kmin - kStart) + hStart;
628 double lkmin = fl * (m_kmin - kStart) + lStart;
629 if ((hkmin >= m_hmin) && (hkmin <= m_hmax) && (lkmin >= m_lmin) && (lkmin <= m_lmax)) {
630 intersections.push_back({{hkmin, m_kmin, lkmin, momkMin}});
631 }
632 }
633 double momkMax = fmom * (m_kmax - kStart) + m_kfmin;
634 if ((momkMax - m_kfmin) * (momkMax - m_kfmax) <= 0) {
635 // hkmax and lkmax
636 double hkmax = fh * (m_kmax - kStart) + hStart;
637 double lkmax = fl * (m_kmax - kStart) + lStart;
638 if ((hkmax >= m_hmin) && (hkmax <= m_hmax) && (lkmax >= m_lmin) && (lkmax <= m_lmax)) {
639 intersections.push_back({{hkmax, m_kmax, lkmax, momkMax}});
640 }
641 }
642 }
643
644 // calculate intersections with planes perpendicular to l
645 if (fabs(lStart - lEnd) > eps) {
646 double fmom = (m_kfmax - m_kfmin) / (lEnd - lStart);
647 double fh = (hEnd - hStart) / (lEnd - lStart);
648 double fk = (kEnd - kStart) / (lEnd - lStart);
649 if (!m_lIntegrated) {
650 for (size_t i = 0; i < lNBins; i++) {
651 double li = m_lX[i];
652 if ((li >= m_lmin) && (li <= m_lmax) && ((lStart - li) * (lEnd - li) < 0)) {
653 double hi = fh * (li - lStart) + hStart;
654 double ki = fk * (li - lStart) + kStart;
655 if ((hi >= m_hmin) && (hi <= m_hmax) && (ki >= m_kmin) && (ki <= m_kmax)) {
656 double momi = fmom * (li - lStart) + m_kfmin;
657 intersections.push_back({{hi, ki, li, momi}});
658 }
659 }
660 }
661 }
662 double momlMin = fmom * (m_lmin - lStart) + m_kfmin;
663 if ((momlMin - m_kfmin) * (momlMin - m_kfmax) <= 0) {
664 // hlmin and klmin
665 double hlmin = fh * (m_lmin - lStart) + hStart;
666 double klmin = fk * (m_lmin - lStart) + kStart;
667 if ((hlmin >= m_hmin) && (hlmin <= m_hmax) && (klmin >= m_kmin) && (klmin <= m_kmax)) {
668 intersections.push_back({{hlmin, klmin, m_lmin, momlMin}});
669 }
670 }
671 double momlMax = fmom * (m_lmax - lStart) + m_kfmin;
672 if ((momlMax - m_kfmin) * (momlMax - m_kfmax) < 0) {
673 // hlmax and klmax
674 double hlmax = fh * (m_lmax - lStart) + hStart;
675 double klmax = fk * (m_lmax - lStart) + kStart;
676 if ((hlmax >= m_hmin) && (hlmax <= m_hmax) && (klmax >= m_kmin) && (klmax <= m_kmax)) {
677 intersections.push_back({{hlmax, klmax, m_lmax, momlMax}});
678 }
679 }
680 }
681
682 // intersections with dE
683 if (!m_dEIntegrated) {
684 for (size_t i = 0; i < eNBins; i++) {
685 double kfi = m_eX[i];
686 if ((kfi - m_kfmin) * (kfi - m_kfmax) <= 0) {
687 double h = qin.X() - qout.X() * kfi;
688 double k = qin.Y() - qout.Y() * kfi;
689 double l = qin.Z() - qout.Z() * kfi;
690 if ((h >= m_hmin) && (h <= m_hmax) && (k >= m_kmin) && (k <= m_kmax) && (l >= m_lmin) && (l <= m_lmax)) {
691 intersections.push_back({{h, k, l, kfi}});
692 }
693 }
694 }
695 }
696
697 // endpoints
698 if ((hStart >= m_hmin) && (hStart <= m_hmax) && (kStart >= m_kmin) && (kStart <= m_kmax) && (lStart >= m_lmin) &&
699 (lStart <= m_lmax)) {
700 intersections.push_back({{hStart, kStart, lStart, m_kfmin}});
701 }
702 if ((hEnd >= m_hmin) && (hEnd <= m_hmax) && (kEnd >= m_kmin) && (kEnd <= m_kmax) && (lEnd >= m_lmin) &&
703 (lEnd <= m_lmax)) {
704 intersections.push_back({{hEnd, kEnd, lEnd, m_kfmax}});
705 }
706
707 // sort intersections by final momentum
708 std::stable_sort(intersections.begin(), intersections.end(), compareMomentum);
709}
710
711} // 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.
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.
A validator which checks that a workspace has a valid instrument.
A property class for workspaces.
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.
DataObjects::MDHistoWorkspace_sptr binInputWS()
Runs the BinMD algorithm on the input to provide the output workspace All slicing algorithm propertie...
Mantid::Kernel::DblMatrix m_rubw
(2*PiRUBW)^-1
API::IMDEventWorkspace_sptr m_inputWS
Input workspace.
double m_Ei
cached values for incident energy and momentum, final momentum min/max
int version() const override
Algorithm's version for identification.
Kernel::V3D m_samplePos
Sample position.
Kernel::V3D m_beamDir
Beam direction.
void calculateNormalization(const std::vector< coord_t > &otherValues, const Kernel::Matrix< coord_t > &affineTrans, uint16_t expInfoIndex)
Computed the normalization for the input workspace.
void createNormalizationWS(const DataObjects::MDHistoWorkspace &dataWS)
Create & cached the normalization workspace.
DataObjects::MDHistoWorkspace_sptr m_normWS
Normalization workspace.
bool m_hIntegrated
flag for integrated h,k,l, dE dimensions
Kernel::Matrix< coord_t > findIntergratedDimensions(const std::vector< coord_t > &otherDimValues, bool &skipNormalization)
Checks the normalization workspace against the indices of the original dimensions.
const std::string name() const override
Algorithm's name for use in the GUI and help.
coord_t m_hmin
limits for h,k,l, dE dimensions
std::vector< coord_t > getValuesFromOtherDimensions(bool &skipNormalization, uint16_t expInfoIndex=0) const
Retrieve logged values from non-HKL dimensions.
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...
void cacheInputs()
Set up starting values for cached variables.
void init() override
Initialize the algorithm's properties.
const std::string summary() const override
Algorithm's summary for use in the GUI and help.
uint16_t m_numExptInfos
number of experiment infos
std::vector< double > m_hX
cached X values along dimensions h,k,l. dE
const std::string category() const override
Algorithm's category for identification.
size_t m_hIdx
index of h,k,l, dE dimensions in the output workspaces
std::string convention
ki-kf for Inelastic convention; kf-ki for Crystallography convention
void cacheDimensionXValues()
Stores the X values from each H,K,L,E dimension as member variables Energy dimension is transformed t...
bool m_accumulate
internal flag to accumulate to an existing workspace
std::string inputEnergyMode() const
Currently looks for the ConvertToMD algorithm in the history.
void exec() override
Execute the algorithm.
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
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
static constexpr double NeutronMass
Mass of the neutron in kg.
static constexpr double h
Planck constant in J*s.
static constexpr double meV
1 meV in Joules.
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
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