Mantid
Loading...
Searching...
No Matches
CalculateCoverageDGS.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 +
9#include "MantidAPI/Run.h"
10#include "MantidAPI/Sample.h"
24
25#include <boost/lexical_cast.hpp>
26
27namespace Mantid::MDAlgorithms {
28using namespace Mantid::Kernel;
30using namespace Mantid::API;
31using namespace Mantid::DataObjects;
32using namespace Mantid::Geometry;
33
34namespace {
35// function to compare two intersections (h,k,l,Momentum) by scattered momentum
36bool compareMomentum(const Mantid::Kernel::VMD &v1, const Mantid::Kernel::VMD &v2) { return (v1[3] < v2[3]); }
37} // namespace
38// Register the algorithm into the AlgorithmFactory
39DECLARE_ALGORITHM(CalculateCoverageDGS)
40
41
44 : m_hmin(0.f), m_hmax(0.f), m_kmin(0.f), m_kmax(0.f), m_lmin(0.f), m_lmax(0.f), m_dEmin(0.f), m_dEmax(0.f),
45 m_Ei(0.), m_ki(0.), m_kfmin(0.), m_kfmax(0.), m_hIntegrated(false), m_kIntegrated(false), m_lIntegrated(false),
46 m_dEIntegrated(false), m_hX(), m_kX(), m_lX(), m_eX(), m_hIdx(-1), m_kIdx(-1), m_lIdx(-1), m_eIdx(-1),
47 m_rubw(3, 3), m_normWS() {}
48
50const std::string CalculateCoverageDGS::name() const { return "CalculateCoverageDGS"; }
51
53int CalculateCoverageDGS::version() const { return 1; }
54
56const std::string CalculateCoverageDGS::category() const { return "Inelastic\\Planning;MDAlgorithms\\Planning"; }
57
59const std::string CalculateCoverageDGS::summary() const {
60 return "Calculate the reciprocal space coverage for direct geometry "
61 "spectrometers";
62}
63
68 const double energyToK = 8.0 * M_PI * M_PI * PhysicalConstants::NeutronMass * PhysicalConstants::meV * 1e-20 /
70
71 auto &hDim = *m_normWS->getDimension(m_hIdx);
72 m_hX.resize(hDim.getNBins());
73 for (size_t i = 0; i < m_hX.size(); ++i) {
74 m_hX[i] = hDim.getX(i);
75 }
76 auto &kDim = *m_normWS->getDimension(m_kIdx);
77 m_kX.resize(kDim.getNBins());
78 for (size_t i = 0; i < m_kX.size(); ++i) {
79 m_kX[i] = kDim.getX(i);
80 }
81 auto &lDim = *m_normWS->getDimension(m_lIdx);
82 m_lX.resize(lDim.getNBins());
83 for (size_t i = 0; i < m_lX.size(); ++i) {
84 m_lX[i] = lDim.getX(i);
85 }
86 // NOTE: store k final instead
87 auto &eDim = *m_normWS->getDimension(m_eIdx);
88 m_eX.resize(eDim.getNBins());
89 for (size_t i = 0; i < m_eX.size(); ++i) {
90 m_eX[i] = std::sqrt(energyToK * (m_Ei - eDim.getX(i)));
91 }
92}
93
97 declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Mantid::Kernel::Direction::Input,
98 std::make_shared<InstrumentValidator>()),
99 "An input workspace.");
100
101 // clang-format off
102 auto mustBe3D = std::make_shared<ArrayLengthValidator<double> >(3);
103 // clang-format on
104 auto mustBePositive = std::make_shared<BoundedValidator<double>>();
105 mustBePositive->setLower(0.0);
106
107 std::vector<double> Q1(3, 0.), Q2(3, 0), Q3(3, 0);
108 Q1[0] = 1.;
109 Q2[1] = 1.;
110 Q3[2] = 1.;
111 declareProperty(std::make_unique<ArrayProperty<double>>("Q1Basis", std::move(Q1), mustBe3D->clone()),
112 "Q1 projection direction in the x,y,z format. Q1, Q2, Q3 "
113 "must not be coplanar");
114 declareProperty(std::make_unique<ArrayProperty<double>>("Q2Basis", std::move(Q2), mustBe3D->clone()),
115 "Q2 projection direction in the x,y,z format. Q1, Q2, Q3 "
116 "must not be coplanar");
117 declareProperty(std::make_unique<ArrayProperty<double>>("Q3Basis", std::move(Q3), std::move(mustBe3D)),
118 "Q3 projection direction in the x,y,z format. Q1, Q2, Q3 "
119 "must not be coplanar");
120 declareProperty(std::make_unique<PropertyWithValue<double>>("IncidentEnergy", EMPTY_DBL(), std::move(mustBePositive),
122 "Incident energy. If set, will override Ei in the input workspace");
123
124 std::vector<std::string> options{"Q1", "Q2", "Q3", "DeltaE"};
125
126 for (int i = 1; i <= 4; i++) {
127 std::string dim("Dimension");
128 dim += boost::lexical_cast<std::string>(i);
129 declareProperty(dim, options[i - 1], std::make_shared<StringListValidator>(options),
130 "Dimension to bin or integrate");
131 declareProperty(dim + "Min", EMPTY_DBL(), dim + " minimum value. If empty will take minimum possible value.");
132 declareProperty(dim + "Max", EMPTY_DBL(), dim + " maximum value. If empty will take maximum possible value.");
133 declareProperty(dim + "Step", EMPTY_DBL(),
134 dim + " step size. If empty the dimension will be "
135 "integrated between minimum and maximum values");
136 }
138 std::make_unique<WorkspaceProperty<Workspace>>("OutputWorkspace", "", Mantid::Kernel::Direction::Output),
139 "A name for the output data MDHistoWorkspace.");
140}
141
145 const double energyToK = 8.0 * M_PI * M_PI * PhysicalConstants::NeutronMass * PhysicalConstants::meV * 1e-20 /
147 // get the limits
148 Mantid::API::MatrixWorkspace_const_sptr inputWS = getProperty("InputWorkspace");
149 convention = Kernel::ConfigService::Instance().getString("Q.convention");
150 // cache two theta and phi
151 const auto &detectorInfo = inputWS->detectorInfo();
152 std::vector<double> tt, phi;
153 for (size_t i = 0; i < detectorInfo.size(); ++i) {
154 if (!detectorInfo.isMasked(i) && !detectorInfo.isMonitor(i)) {
155 const auto &detector = detectorInfo.detector(i);
156 tt.emplace_back(detector.getTwoTheta(V3D(0, 0, 0), V3D(0, 0, 1)));
157 phi.emplace_back(detector.getPhi());
158 }
159 }
160
161 double ttmax = *(std::max_element(tt.begin(), tt.end()));
162 m_Ei = getProperty("IncidentEnergy");
163 if (m_Ei == EMPTY_DBL()) {
164 if (inputWS->run().hasProperty("Ei")) {
165 Kernel::Property *eiprop = inputWS->run().getProperty("Ei");
166 m_Ei = boost::lexical_cast<double>(eiprop->value());
167 if (m_Ei <= 0) {
168 throw std::invalid_argument("Ei stored in the workspace is not positive");
169 }
170 } else {
171 throw std::invalid_argument("Could not find Ei in the workspace. Please enter a positive value");
172 }
173 }
174
175 // Check if there is duplicate elements in dimension selection, and calculate
176 // the affine matrix
177 Kernel::Matrix<coord_t> affineMat(4, 4);
178 double q1min(0.), q1max(0.), q2min(0.), q2max(0.), q3min(0.), q3max(0.);
179 double q1step(0.), q2step(0.), q3step(0.), dEstep(0.);
180 size_t q1NumBins = 1, q2NumBins = 1, q3NumBins = 1, dENumBins = 1;
181 for (int i = 1; i <= 4; i++) {
182 std::string dim("Dimension");
183 dim += boost::lexical_cast<std::string>(i);
184 std::string dimensioni = getProperty(dim);
185 if (dimensioni == "Q1") {
186 affineMat[i - 1][0] = 1.;
187 q1min = getProperty(dim + "Min");
188 q1max = getProperty(dim + "Max");
189 q1step = getProperty(dim + "Step");
190 m_hIdx = i - 1;
191 }
192 if (dimensioni == "Q2") {
193 affineMat[i - 1][1] = 1.;
194 q2min = getProperty(dim + "Min");
195 q2max = getProperty(dim + "Max");
196 q2step = getProperty(dim + "Step");
197 m_kIdx = i - 1;
198 }
199 if (dimensioni == "Q3") {
200 affineMat[i - 1][2] = 1.;
201 q3min = getProperty(dim + "Min");
202 q3max = getProperty(dim + "Max");
203 q3step = getProperty(dim + "Step");
204 m_lIdx = i - 1;
205 }
206 if (dimensioni == "DeltaE") {
207 affineMat[i - 1][3] = 1.;
208 m_eIdx = i - 1;
209 double dmin = getProperty(dim + "Min");
210 if (dmin == EMPTY_DBL()) {
211 m_dEmin = -static_cast<coord_t>(m_Ei);
212 } else {
213 m_dEmin = static_cast<coord_t>(dmin);
214 }
215 double dmax = getProperty(dim + "Max");
216 if (dmax == EMPTY_DBL()) {
217 m_dEmax = static_cast<coord_t>(m_Ei);
218 } else {
219 m_dEmax = static_cast<coord_t>(dmax);
220 }
221 dEstep = getProperty(dim + "Step");
222 if (dEstep == EMPTY_DBL()) {
223 dENumBins = 1;
224 } else {
225 dENumBins = static_cast<size_t>((m_dEmax - m_dEmin) / dEstep);
226 if (dEstep * static_cast<double>(dENumBins) + m_dEmin < m_dEmax) {
227 dENumBins += 1;
228 m_dEmax = static_cast<coord_t>(dEstep * static_cast<double>(dENumBins)) + m_dEmin;
229 }
230 }
231 }
232 }
233
234 if (affineMat.determinant() == 0.) {
235 g_log.debug() << affineMat;
236 throw std::invalid_argument("Please make sure each dimension is selected only once.");
237 }
238
239 // Qmax is at kf=kfmin or kf=kfmax
240 m_ki = std::sqrt(energyToK * m_Ei);
241 if (m_Ei > m_dEmin) {
242 m_kfmin = std::sqrt(energyToK * (m_Ei - m_dEmin));
243 } else {
244 m_kfmin = 0.;
245 }
246 if (m_Ei > m_dEmax) {
247 m_kfmax = std::sqrt(energyToK * (m_Ei - m_dEmax));
248 } else {
249 m_kfmax = 0;
250 }
251
252 double QmaxTemp = sqrt(m_ki * m_ki + m_kfmin * m_kfmin - 2 * m_ki * m_kfmin * cos(ttmax));
253 double Qmax = QmaxTemp;
254 QmaxTemp = sqrt(m_ki * m_ki + m_kfmax * m_kfmax - 2 * m_ki * m_kfmax * cos(ttmax));
255 if (QmaxTemp > Qmax)
256 Qmax = QmaxTemp;
257
258 // get goniometer
259 DblMatrix gon = DblMatrix(3, 3, true);
260 if (inputWS->run().getGoniometer().isDefined()) {
261 gon = inputWS->run().getGoniometerMatrix();
262 }
263
264 // get the UB
265 DblMatrix UB = DblMatrix(3, 3, true) * (0.5 / M_PI);
266 if (inputWS->sample().hasOrientedLattice()) {
267 UB = inputWS->sample().getOrientedLattice().getUB();
268 }
269
270 // get the W matrix
271 DblMatrix W = DblMatrix(3, 3);
272 std::vector<double> Q1Basis = getProperty("Q1Basis");
273 std::vector<double> Q2Basis = getProperty("Q2Basis");
274 std::vector<double> Q3Basis = getProperty("Q3Basis");
275 W.setColumn(0, Q1Basis);
276 W.setColumn(1, Q2Basis);
277 W.setColumn(2, Q3Basis);
278
279 m_rubw = gon * UB * W * (2.0 * M_PI);
280
281 // calculate maximum original limits
283 ol.setUB(m_rubw);
284 m_hmin = static_cast<coord_t>(-Qmax * ol.a());
285 m_hmax = static_cast<coord_t>(Qmax * ol.a());
286 m_kmin = static_cast<coord_t>(-Qmax * ol.b());
287 m_kmax = static_cast<coord_t>(Qmax * ol.b());
288 m_lmin = static_cast<coord_t>(-Qmax * ol.c());
289 m_lmax = static_cast<coord_t>(Qmax * ol.c());
290 m_rubw.Invert();
291 // adjust Q steps/dimensions
292 if (q1min == EMPTY_DBL()) {
293 q1min = m_hmin;
294 }
295 if (q1max == EMPTY_DBL()) {
296 q1max = m_hmax;
297 }
298 if (q1min >= q1max) {
299 throw std::invalid_argument("Q1max has to be greater than Q1min");
300 }
301 if (q1step == EMPTY_DBL()) {
302 q1NumBins = 1;
303 } else {
304 q1NumBins = static_cast<size_t>((q1max - q1min) / q1step);
305 if (q1step * static_cast<double>(q1NumBins) + q1min < q1max) {
306 q1NumBins += 1;
307 q1max = q1step * static_cast<double>(q1NumBins) + q1min;
308 }
309 }
310
311 if (q2min == EMPTY_DBL()) {
312 q2min = m_kmin;
313 }
314 if (q2max == EMPTY_DBL()) {
315 q2max = m_kmax;
316 }
317 if (q2min >= q2max) {
318 throw std::invalid_argument("Q2max has to be greater than Q2min");
319 }
320 if (q2step == EMPTY_DBL()) {
321 q2NumBins = 1;
322 } else {
323 q2NumBins = static_cast<size_t>((q2max - q2min) / q2step);
324 if (q2step * static_cast<double>(q2NumBins) + q2min < q2max) {
325 q2NumBins += 1;
326 q2max = q2step * static_cast<double>(q2NumBins) + q2min;
327 }
328 }
329
330 if (q3min == EMPTY_DBL()) {
331 q3min = m_lmin;
332 }
333 if (q3max == EMPTY_DBL()) {
334 q3max = m_lmax;
335 }
336 if (q3min >= q3max) {
337 throw std::invalid_argument("Q3max has to be greater than Q3min");
338 }
339 if (q3step == EMPTY_DBL()) {
340 q3NumBins = 1;
341 } else {
342 q3NumBins = static_cast<size_t>((q3max - q3min) / q3step);
343 if (q3step * static_cast<double>(q3NumBins) + q3min < q3max) {
344 q3NumBins += 1;
345 q3max = q3step * static_cast<double>(q3NumBins) + q3min;
346 }
347 }
348
349 // create the output workspace
350 std::vector<Mantid::Geometry::MDHistoDimension_sptr> binDimensions;
351
352 // Define frames
353 Mantid::Geometry::GeneralFrame frame1("Q1", "");
354 Mantid::Geometry::GeneralFrame frame2("Q2", "");
355 Mantid::Geometry::GeneralFrame frame3("Q3", "");
356 Mantid::Geometry::GeneralFrame frame4("meV", "");
357 auto out1 = std::make_shared<MDHistoDimension>("Q1", "Q1", frame1, static_cast<coord_t>(q1min),
358 static_cast<coord_t>(q1max), q1NumBins);
359 auto out2 = std::make_shared<MDHistoDimension>("Q2", "Q2", frame2, static_cast<coord_t>(q2min),
360 static_cast<coord_t>(q2max), q2NumBins);
361 auto out3 = std::make_shared<MDHistoDimension>("Q3", "Q3", frame3, static_cast<coord_t>(q3min),
362 static_cast<coord_t>(q3max), q3NumBins);
363 auto out4 = std::make_shared<MDHistoDimension>("DeltaE", "DeltaE", frame4, static_cast<coord_t>(m_dEmin),
364 static_cast<coord_t>(m_dEmax), dENumBins);
365
366 for (size_t row = 0; row <= 3; row++) {
367 if (affineMat[row][0] == 1.) {
368 binDimensions.emplace_back(out1);
369 }
370 if (affineMat[row][1] == 1.) {
371 binDimensions.emplace_back(out2);
372 }
373 if (affineMat[row][2] == 1.) {
374 binDimensions.emplace_back(out3);
375 }
376 if (affineMat[row][3] == 1.) {
377 binDimensions.emplace_back(out4);
378 }
379 }
380
381 m_normWS = std::make_shared<MDHistoWorkspace>(binDimensions);
382 m_normWS->setTo(0., 0., 0.);
383 setProperty("OutputWorkspace", std::dynamic_pointer_cast<Workspace>(m_normWS));
384
386
387 const auto ndets = static_cast<int64_t>(tt.size());
388
390 for (int64_t i = 0; i < ndets; i++) {
392 auto intersections = calculateIntersections(tt[i], phi[i]);
393 if (intersections.empty())
394 continue;
395 auto intersectionsBegin = intersections.begin();
396 for (auto it = intersectionsBegin + 1; it != intersections.end(); ++it) {
397 const auto &curIntSec = *it;
398 const auto &prevIntSec = *(it - 1);
399 // the full vector isn't used so compute only what is necessary
400 double delta = curIntSec[3] - prevIntSec[3];
401 if (delta < 1e-10)
402 continue; // Assume zero contribution if difference is small
403 // Average between two intersections for final position
404 std::vector<coord_t> pos(4);
405 std::transform(curIntSec.getBareArray(), curIntSec.getBareArray() + 4, prevIntSec.getBareArray(), pos.begin(),
406 [](const double lhs, const double rhs) { return static_cast<coord_t>(0.5 * (lhs + rhs)); });
407 // transform kf to energy transfer
408 pos[3] = static_cast<coord_t>(m_Ei - pos[3] * pos[3] / energyToK);
409
410 std::vector<coord_t> posNew = affineMat * pos;
411 size_t linIndex = m_normWS->getLinearIndexAtCoord(posNew.data());
412 if (linIndex == size_t(-1))
413 continue;
414 PARALLEL_CRITICAL(updateMD) { m_normWS->setSignalAt(linIndex, 1.); }
415 }
417 }
419}
420
429std::vector<Kernel::VMD> CalculateCoverageDGS::calculateIntersections(const double theta, const double phi) {
430 V3D qout(sin(theta) * cos(phi), sin(theta) * sin(phi), cos(theta)), qin(0., 0., m_ki);
431 qout = m_rubw * qout;
432 qin = m_rubw * qin;
433 if (convention == "Crystallography") {
434 qout *= -1;
435 qin *= -1;
436 }
437 double hStart = qin.X() - qout.X() * m_kfmin, hEnd = qin.X() - qout.X() * m_kfmax;
438 double kStart = qin.Y() - qout.Y() * m_kfmin, kEnd = qin.Y() - qout.Y() * m_kfmax;
439 double lStart = qin.Z() - qout.Z() * m_kfmin, lEnd = qin.Z() - qout.Z() * m_kfmax;
440 double eps = 1e-10;
441 auto hNBins = m_hX.size();
442 auto kNBins = m_kX.size();
443 auto lNBins = m_lX.size();
444 auto eNBins = m_eX.size();
445 std::vector<Kernel::VMD> intersections;
446 intersections.reserve(hNBins + kNBins + lNBins + eNBins + 8); // 8 is 3*(min,max for each Q component)+kfmin+kfmax
447
448 // calculate intersections with planes perpendicular to h
449 if (fabs(hStart - hEnd) > eps) {
450 double fmom = (m_kfmax - m_kfmin) / (hEnd - hStart);
451 double fk = (kEnd - kStart) / (hEnd - hStart);
452 double fl = (lEnd - lStart) / (hEnd - hStart);
453 if (!m_hIntegrated) {
454 for (size_t i = 0; i < hNBins; i++) {
455 double hi = m_hX[i];
456 if ((hi >= m_hmin) && (hi <= m_hmax) && ((hStart - hi) * (hEnd - hi) < 0)) {
457 // if hi is between hStart and hEnd, then ki and li will be between
458 // kStart, kEnd and lStart, lEnd and momi will be between m_kfmin and
459 // m_kfmax
460 double ki = fk * (hi - hStart) + kStart;
461 double li = fl * (hi - hStart) + lStart;
462 if ((ki >= m_kmin) && (ki <= m_kmax) && (li >= m_lmin) && (li <= m_lmax)) {
463 double momi = fmom * (hi - hStart) + m_kfmin;
464 intersections.emplace_back(hi, ki, li, momi);
465 }
466 }
467 }
468 }
469 double momhMin = fmom * (m_hmin - hStart) + m_kfmin;
470 if ((momhMin - m_kfmin) * (momhMin - m_kfmax) < 0) // m_kfmin>m_kfmax
471 {
472 // khmin and lhmin
473 double khmin = fk * (m_hmin - hStart) + kStart;
474 double lhmin = fl * (m_hmin - hStart) + lStart;
475 if ((khmin >= m_kmin) && (khmin <= m_kmax) && (lhmin >= m_lmin) && (lhmin <= m_lmax)) {
476 intersections.emplace_back(m_hmin, khmin, lhmin, momhMin);
477 }
478 }
479 double momhMax = fmom * (m_hmax - hStart) + m_kfmin;
480 if ((momhMax - m_kfmin) * (momhMax - m_kfmax) <= 0) {
481 // khmax and lhmax
482 double khmax = fk * (m_hmax - hStart) + kStart;
483 double lhmax = fl * (m_hmax - hStart) + lStart;
484 if ((khmax >= m_kmin) && (khmax <= m_kmax) && (lhmax >= m_lmin) && (lhmax <= m_lmax)) {
485 intersections.emplace_back(m_hmax, khmax, lhmax, momhMax);
486 }
487 }
488 }
489
490 // calculate intersections with planes perpendicular to k
491 if (fabs(kStart - kEnd) > eps) {
492 double fmom = (m_kfmax - m_kfmin) / (kEnd - kStart);
493 double fh = (hEnd - hStart) / (kEnd - kStart);
494 double fl = (lEnd - lStart) / (kEnd - kStart);
495 if (!m_kIntegrated) {
496 for (size_t i = 0; i < kNBins; i++) {
497 double ki = m_kX[i];
498 if ((ki >= m_kmin) && (ki <= m_kmax) && ((kStart - ki) * (kEnd - ki) < 0)) {
499 // if ki is between kStart and kEnd, then hi and li will be between
500 // hStart, hEnd and lStart, lEnd and momi will be between m_kfmin and
501 // m_kfmax
502 double hi = fh * (ki - kStart) + hStart;
503 double li = fl * (ki - kStart) + lStart;
504 if ((hi >= m_hmin) && (hi <= m_hmax) && (li >= m_lmin) && (li <= m_lmax)) {
505 double momi = fmom * (ki - kStart) + m_kfmin;
506 intersections.emplace_back(hi, ki, li, momi);
507 }
508 }
509 }
510 }
511 double momkMin = fmom * (m_kmin - kStart) + m_kfmin;
512 if ((momkMin - m_kfmin) * (momkMin - m_kfmax) < 0) {
513 // hkmin and lkmin
514 double hkmin = fh * (m_kmin - kStart) + hStart;
515 double lkmin = fl * (m_kmin - kStart) + lStart;
516 if ((hkmin >= m_hmin) && (hkmin <= m_hmax) && (lkmin >= m_lmin) && (lkmin <= m_lmax)) {
517 intersections.emplace_back(hkmin, m_kmin, lkmin, momkMin);
518 }
519 }
520 double momkMax = fmom * (m_kmax - kStart) + m_kfmin;
521 if ((momkMax - m_kfmin) * (momkMax - m_kfmax) <= 0) {
522 // hkmax and lkmax
523 double hkmax = fh * (m_kmax - kStart) + hStart;
524 double lkmax = fl * (m_kmax - kStart) + lStart;
525 if ((hkmax >= m_hmin) && (hkmax <= m_hmax) && (lkmax >= m_lmin) && (lkmax <= m_lmax)) {
526 intersections.emplace_back(hkmax, m_kmax, lkmax, momkMax);
527 }
528 }
529 }
530
531 // calculate intersections with planes perpendicular to l
532 if (fabs(lStart - lEnd) > eps) {
533 double fmom = (m_kfmax - m_kfmin) / (lEnd - lStart);
534 double fh = (hEnd - hStart) / (lEnd - lStart);
535 double fk = (kEnd - kStart) / (lEnd - lStart);
536 if (!m_lIntegrated) {
537 for (size_t i = 0; i < lNBins; i++) {
538 double li = m_lX[i];
539 if ((li >= m_lmin) && (li <= m_lmax) && ((lStart - li) * (lEnd - li) < 0)) {
540 double hi = fh * (li - lStart) + hStart;
541 double ki = fk * (li - lStart) + kStart;
542 if ((hi >= m_hmin) && (hi <= m_hmax) && (ki >= m_kmin) && (ki <= m_kmax)) {
543 double momi = fmom * (li - lStart) + m_kfmin;
544 intersections.emplace_back(hi, ki, li, momi);
545 }
546 }
547 }
548 }
549 double momlMin = fmom * (m_lmin - lStart) + m_kfmin;
550 if ((momlMin - m_kfmin) * (momlMin - m_kfmax) <= 0) {
551 // hlmin and klmin
552 double hlmin = fh * (m_lmin - lStart) + hStart;
553 double klmin = fk * (m_lmin - lStart) + kStart;
554 if ((hlmin >= m_hmin) && (hlmin <= m_hmax) && (klmin >= m_kmin) && (klmin <= m_kmax)) {
555 intersections.emplace_back(hlmin, klmin, m_lmin, momlMin);
556 }
557 }
558 double momlMax = fmom * (m_lmax - lStart) + m_kfmin;
559 if ((momlMax - m_kfmin) * (momlMax - m_kfmax) < 0) {
560 // hlmax and klmax
561 double hlmax = fh * (m_lmax - lStart) + hStart;
562 double klmax = fk * (m_lmax - lStart) + kStart;
563 if ((hlmax >= m_hmin) && (hlmax <= m_hmax) && (klmax >= m_kmin) && (klmax <= m_kmax)) {
564 intersections.emplace_back(hlmax, klmax, m_lmax, momlMax);
565 }
566 }
567 }
568
569 // intersections with dE
570 if (!m_dEIntegrated) {
571 for (size_t i = 0; i < eNBins; i++) {
572 double kfi = m_eX[i];
573 if ((kfi - m_kfmin) * (kfi - m_kfmax) <= 0) {
574 double h = qin.X() - qout.X() * kfi;
575 double k = qin.Y() - qout.Y() * kfi;
576 double l = qin.Z() - qout.Z() * kfi;
577 if ((h >= m_hmin) && (h <= m_hmax) && (k >= m_kmin) && (k <= m_kmax) && (l >= m_lmin) && (l <= m_lmax)) {
578 intersections.emplace_back(h, k, l, kfi);
579 }
580 }
581 }
582 }
583
584 // endpoints
585 if ((hStart >= m_hmin) && (hStart <= m_hmax) && (kStart >= m_kmin) && (kStart <= m_kmax) && (lStart >= m_lmin) &&
586 (lStart <= m_lmax)) {
587 intersections.emplace_back(hStart, kStart, lStart, m_kfmin);
588 }
589 if ((hEnd >= m_hmin) && (hEnd <= m_hmax) && (kEnd >= m_kmin) && (kEnd <= m_kmax) && (lEnd >= m_lmin) &&
590 (lEnd <= m_lmax)) {
591 intersections.emplace_back(hEnd, kEnd, lEnd, m_kfmax);
592 }
593
594 // sort intersections by final momentum
595 using IterType = std::vector<Mantid::Kernel::VMD>::iterator;
596 std::stable_sort<IterType, bool (*)(const Mantid::Kernel::VMD &, const Mantid::Kernel::VMD &)>(
597 intersections.begin(), intersections.end(), compareMomentum);
598
599 return intersections;
600}
601
602} // namespace Mantid::MDAlgorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
const std::vector< double > & rhs
#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_CRITICAL(name)
#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
Kernel::Logger & g_log
Definition: Algorithm.h:451
A property class for workspaces.
GeneralFrame : Any MDFrame that isn't related to momemtum transfer.
Definition: GeneralFrame.h:21
Class to implement UB matrix.
void setUB(const Kernel::DblMatrix &newUB)
Sets the UB matrix and recalculates lattice parameters.
Support for a property that holds an array of values.
Definition: ArrayProperty.h:28
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void debug(const std::string &msg)
Logs at debug level.
Definition: Logger.cpp:114
Numerical Matrix class.
Definition: Matrix.h:42
T determinant() const
Calculate the determinant.
Definition: Matrix.cpp:1048
T Invert()
LU inversion routine.
Definition: Matrix.cpp:924
void setColumn(const size_t nCol, const std::vector< T > &newCol)
Definition: Matrix.cpp:675
The concrete, templated class for properties.
Base class for properties.
Definition: Property.h:94
virtual std::string value() const =0
Returns the value of the property as a string.
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
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
CalculateCoverageDGS : Calculate coverage for single crystal direct geometry scattering.
bool m_hIntegrated
flag for integrated h,k,l,dE dimensions
const std::string category() const override
Algorithm's category for identification.
Mantid::DataObjects::MDHistoWorkspace_sptr m_normWS
Normalization workspace (this is the coverage workspace)
const std::string name() const override
Algorithms name for identification.
coord_t m_hmin
limits for h,k,l,dE dimensions
void cacheDimensionXValues()
Stores the X values from each H,K,L dimension as member variables.
double m_Ei
cached values for incident energy and momentum, final momentum min/max
Mantid::Kernel::DblMatrix m_rubw
(2*PiRUBW)^-1
int version() const override
Algorithm's version for identification.
const std::string summary() const override
Algorithm's summary for use in the GUI and help.
std::vector< double > m_hX
cached X values along dimensions h,k,l, dE
size_t m_hIdx
index of h,k,l,dE dimensions in the output workspaces
void exec() override
Execute the algorithm.
void init() override
Initialize the algorithm's properties.
std::vector< Kernel::VMD > calculateIntersections(const double theta, const double phi)
Calculate the points of intersection for the given detector with cuboid surrounding the detector posi...
std::shared_ptr< const MatrixWorkspace > MatrixWorkspace_const_sptr
shared pointer to the matrix workspace base class (const version)
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::Matrix< double > DblMatrix
Definition: Matrix.h:206
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
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