Mantid
Loading...
Searching...
No Matches
FindSXPeaksHelper.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 +
11#include "MantidKernel/Logger.h"
13#include "MantidKernel/Unit.h"
15
16#include "MantidTypes/SpectrumDefinition.h"
17
18#define BOOST_ALLOW_DEPRECATED_HEADERS
19#include <boost/graph/adjacency_list.hpp>
20#undef BOOST_ALLOW_DEPRECATED_HEADERS
21
22#include <boost/graph/connected_components.hpp>
23#include <cmath>
24
25namespace {
26
27const double TWO_PI = 2 * M_PI;
28
29bool isDifferenceLargerThanTolerance(const double angle1, const double angle2, const double tolerance) {
30 auto difference = std::abs(angle1 - angle2);
31
32 // If we have more than 360 degree angle difference then we need to wrap it
33 // back to 360
34 if (difference > TWO_PI) {
35 difference = std::fmod(difference, TWO_PI);
36 }
37
38 // If we have more than 180 degrees then we must take the smaller angle
39 if (difference > M_PI) {
40 difference = TWO_PI - difference;
41 }
42
43 return difference > tolerance;
44}
45} // namespace
46
47using namespace boost;
49
51
52Mantid::Kernel::Logger g_log("FindSXPeaksHelper");
53
54/* ------------------------------------------------------------------------------------------
55 * Single Crystal peak representation
56 * ------------------------------------------------------------------------------------------
57 */
58
68SXPeak::SXPeak(double t, double phi, double intensity, const std::vector<int> &spectral, const size_t wsIndex,
69 const API::SpectrumInfo &spectrumInfo)
70 : m_tof(t), m_phi(phi), m_intensity(intensity), m_spectra(spectral), m_wsIndex(wsIndex) {
71 // Sanity checks
72 if (intensity < 0) {
73 throw std::invalid_argument("SXPeak: Cannot have an intensity < 0");
74 }
75 if (spectral.empty()) {
76 throw std::invalid_argument("SXPeak: Cannot have zero sized spectral list");
77 }
78 if (!spectrumInfo.hasDetectors(m_wsIndex)) {
79 throw std::invalid_argument("SXPeak: Spectrum at ws index " + std::to_string(wsIndex) + " doesn't have detectors");
80 }
81
82 const auto l1 = spectrumInfo.l1();
83 const auto l2 = spectrumInfo.l2(m_wsIndex);
84
85 m_twoTheta = spectrumInfo.twoTheta(m_wsIndex);
86 m_LTotal = l1 + l2;
87 if (m_LTotal < 0) {
88 throw std::invalid_argument("SXPeak: Cannot have detector distance < 0");
89 }
90 m_detId = spectrumInfo.detector(m_wsIndex).getID();
91 m_nPixels = 1;
92
94 const auto unit = Mantid::Kernel::UnitFactory::Instance().create("dSpacing");
96 spectrumInfo.getDetectorValues(tof, *unit, Kernel::DeltaEMode::Elastic, false, m_wsIndex, pmap);
97 unit->initialize(l1, 0, pmap);
98 try {
99 m_dSpacing = unit->singleFromTOF(m_tof);
100 } catch (std::exception &) {
101 m_dSpacing = 0;
102 }
103
104 const auto samplePos = spectrumInfo.samplePosition();
105 const auto sourcePos = spectrumInfo.sourcePosition();
106 const auto detPos = spectrumInfo.position(m_wsIndex);
107 // Normalized beam direction
108 const auto beamDir = normalize(samplePos - sourcePos);
109 // Normalized detector direction
110 const auto detDir = normalize(detPos - samplePos);
111 m_unitWaveVector = beamDir - detDir;
112 m_qConvention = Kernel::ConfigService::Instance().getString("Q.convention");
113}
114
120bool SXPeak::compare(const SXPeak &rhs, double tolerance) const {
121 if (std::abs(m_tof / m_nPixels - rhs.m_tof / rhs.m_nPixels) > tolerance * m_tof / m_nPixels)
122 return false;
123 if (std::abs(m_phi / m_nPixels - rhs.m_phi / rhs.m_nPixels) > tolerance * m_phi / m_nPixels)
124 return false;
125 if (std::abs(m_twoTheta / m_nPixels - rhs.m_twoTheta / rhs.m_nPixels) > tolerance * m_twoTheta / m_nPixels)
126 return false;
127 return true;
128}
129
130bool SXPeak::compare(const SXPeak &rhs, const double xTolerance, const double phiTolerance, const double thetaTolerance,
131 const XAxisUnit units) const {
132
133 const auto x_1 = (units == XAxisUnit::TOF) ? m_tof : m_dSpacing;
134 const auto x_2 = (units == XAxisUnit::TOF) ? rhs.m_tof : rhs.m_dSpacing;
135
136 if (std::abs(x_1 - x_2) > xTolerance) {
137 return false;
138 }
139
140 if (isDifferenceLargerThanTolerance(m_phi, rhs.m_phi, phiTolerance)) {
141 return false;
142 }
143 if (isDifferenceLargerThanTolerance(m_twoTheta, rhs.m_twoTheta, thetaTolerance)) {
144 return false;
145 }
146 return true;
147}
148
154 double qSign = 1.0;
155 if (m_qConvention == "Crystallography") {
156 qSign = -1.0;
157 }
158 double vi = m_LTotal / (m_tof * 1e-6);
159 // wavenumber = h_bar / mv
161 // in angstroms
162 wi *= 1e10;
163 // wavevector=1/wavenumber = 2pi/wavelength
164 double wvi = 1.0 / wi;
165 // Now calculate the wavevector of the scattered neutron
166 return m_unitWaveVector * wvi * qSign;
167}
168
174 m_tof += rhs.m_tof;
175 m_phi += rhs.m_phi;
176 m_twoTheta += rhs.m_twoTheta;
177 m_intensity += rhs.m_intensity;
178 m_LTotal += rhs.m_LTotal;
179 m_nPixels += 1;
180 m_spectra.insert(m_spectra.end(), rhs.m_spectra.cbegin(), rhs.m_spectra.cend());
181 return *this;
182}
183
186 m_tof /= m_nPixels;
187 m_phi /= m_nPixels;
191 m_nPixels = 1;
192}
193
197const double &SXPeak::getIntensity() const { return m_intensity; }
198
203
204PeakContainer::PeakContainer(const HistogramData::HistogramY &y)
205 : m_y(y), m_startIndex(0), m_stopIndex(m_y.size() - 1), m_maxIndex(0) {}
206
208 m_startIndex = std::distance(m_y.begin(), item);
210 m_maxSignal = *item;
211}
212
214 if (*item > m_maxSignal) {
215 m_maxIndex = std::distance(m_y.begin(), item);
216 m_maxSignal = *item;
217 }
218}
219
221 if (*item > m_maxSignal) {
222 m_maxIndex = std::distance(m_y.begin(), item);
223 m_maxSignal = *item;
224 }
225 m_stopIndex = std::distance(m_y.begin(), item);
226
227 // Peak end is one back though
228 --m_stopIndex;
229}
230
232 // If we didn't record anything then the start iterator is at the end
233 if (m_startIndex >= m_y.size()) {
234 return 0;
235 }
236
237 return m_stopIndex - m_startIndex;
238}
239
240yIt PeakContainer::getMaxIterator() const { return m_y.begin() + m_maxIndex; }
241
242/* ------------------------------------------------------------------------------------------
243 * Background
244 * ------------------------------------------------------------------------------------------
245 */
246AbsoluteBackgroundStrategy::AbsoluteBackgroundStrategy(const double background) : m_background(background) {}
247
249 const HistogramData::HistogramY & /*y*/) const {
250 return intensity < m_background;
251}
252
254 : m_backgroundMultiplier(backgroundMultiplier) {}
255
257 const HistogramData::HistogramY &y) const {
258 auto background = 0.5 * (1.0 + y.front() + y.back());
259 background *= m_backgroundMultiplier;
260 return intensity < background;
261}
262
263/* ------------------------------------------------------------------------------------------
264 * Peak Finding Strategy
265 * ------------------------------------------------------------------------------------------
266 */
267
269 const API::SpectrumInfo &spectrumInfo, const double minValue,
270 const double maxValue, const XAxisUnit units)
271 : m_backgroundStrategy(backgroundStrategy), m_minValue(minValue), m_maxValue(maxValue),
272 m_spectrumInfo(spectrumInfo), m_units(units) {}
273
274PeakList PeakFindingStrategy::findSXPeaks(const HistogramData::HistogramX &x, const HistogramData::HistogramY &y,
275 const int workspaceIndex) const {
276 // ---------------------------------------
277 // Get the lower and upper bound iterators
278 // ---------------------------------------
279 auto boundsIterator = getBounds(x);
280 auto lowit = boundsIterator.first;
281 auto highit = boundsIterator.second;
282
283 // If range specified doesn't overlap with this spectrum then bail out
284 if (lowit == x.end() || highit == x.begin()) {
285 return PeakList();
286 }
287
288 // Upper limit is the bin before, i.e. the last value smaller than MaxRange
289 --highit;
290
291 // ---------------------------------------
292 // Perform the search of the peaks
293 // ---------------------------------------
294 return dofindSXPeaks(x, y, lowit, highit, workspaceIndex);
295}
296
297BoundsIterator PeakFindingStrategy::getBounds(const HistogramData::HistogramX &x) const {
298 // Find the range [min,max]
299 auto lowit = (m_minValue == EMPTY_DBL()) ? x.begin() : std::lower_bound(x.begin(), x.end(), m_minValue);
300 using std::placeholders::_1;
301 auto highit = (m_maxValue == EMPTY_DBL())
302 ? x.end()
303 : std::find_if(lowit, x.end(), std::bind(std::greater<double>(), _1, m_maxValue));
304
305 return std::make_pair(lowit, highit);
306}
307
315double PeakFindingStrategy::calculatePhi(size_t workspaceIndex) const {
316 double phi;
317
318 // Get the detectors for the workspace index
319 const auto &spectrumDefinition = m_spectrumInfo.spectrumDefinition(workspaceIndex);
320 const auto numberOfDetectors = spectrumDefinition.size();
321 const auto &det = m_spectrumInfo.detector(workspaceIndex);
322 if (numberOfDetectors == 1) {
323 phi = det.getPhi();
324 } else {
325 // Have to average the value for phi
326 auto detectorGroup = dynamic_cast<const Mantid::Geometry::DetectorGroup *>(&det);
327 if (!detectorGroup) {
328 throw std::runtime_error("Could not cast to detector group");
329 }
330 phi = detectorGroup->getPhi();
331 }
332
333 if (phi < 0) {
334 phi += 2.0 * M_PI;
335 }
336 return phi;
337}
338
339double PeakFindingStrategy::getXValue(const HistogramData::HistogramX &x, const size_t peakLocation) const {
340 auto leftBinPosition = x.begin() + peakLocation;
341 const double leftBinEdge = *leftBinPosition;
342 const double rightBinEdge = *std::next(leftBinPosition);
343 return 0.5 * (leftBinEdge + rightBinEdge);
344}
345
346double PeakFindingStrategy::convertToTOF(const double xValue, const size_t workspaceIndex) const {
347 if (m_units == XAxisUnit::TOF) {
348 // we're already using TOF units
349 return xValue;
350 } else {
351 const auto unit = UnitFactory::Instance().create("dSpacing");
354 m_spectrumInfo.getDetectorValues(*unit, tof, Kernel::DeltaEMode::Elastic, false, workspaceIndex, pmap);
355 // we're using d-spacing, convert the point to TOF
356 unit->initialize(m_spectrumInfo.l1(), 0, pmap);
357 return unit->singleToTOF(xValue);
358 }
359}
360
362 const API::SpectrumInfo &spectrumInfo, const double minValue,
363 const double maxValue, const XAxisUnit units)
364 : PeakFindingStrategy(backgroundStrategy, spectrumInfo, minValue, maxValue, units) {}
365
366PeakList StrongestPeaksStrategy::dofindSXPeaks(const HistogramData::HistogramX &x, const HistogramData::HistogramY &y,
367 Bound low, Bound high, const int workspaceIndex) const {
368 auto distmin = std::distance(x.begin(), low);
369 auto distmax = std::distance(x.begin(), high);
370
371 // Find the max element
372 auto maxY = (y.size() > 1) ? std::max_element(y.begin() + distmin, y.begin() + distmax) : y.begin();
373
374 // Perform a check against the background
375 double intensity = (*maxY);
376 if (m_backgroundStrategy->isBelowBackground(intensity, y)) {
377 return PeakList();
378 }
379
380 // Create the SXPeak information
381 const auto distance = std::distance(y.begin(), maxY);
382 const auto xValue = getXValue(x, distance);
383 const auto tof = convertToTOF(xValue, workspaceIndex);
384 const double phi = calculatePhi(workspaceIndex);
385
386 std::vector<int> specs(1, workspaceIndex);
387 std::vector<SXPeak> peaks;
388 peaks.emplace_back(tof, phi, *maxY, specs, workspaceIndex, m_spectrumInfo);
389 return peaks;
390}
391
392AllPeaksStrategy::AllPeaksStrategy(const BackgroundStrategy *backgroundStrategy, const API::SpectrumInfo &spectrumInfo,
393 const double minValue, const double maxValue, const XAxisUnit units)
394 : PeakFindingStrategy(backgroundStrategy, spectrumInfo, minValue, maxValue, units) {
395 // We only allow the AbsoluteBackgroundStrategy for now
396 if (!dynamic_cast<const AbsoluteBackgroundStrategy *>(m_backgroundStrategy)) {
397 throw std::invalid_argument("The AllPeaksStrategy has to be initialized "
398 "with the AbsoluteBackgroundStrategy.");
399 }
400}
401
402PeakList AllPeaksStrategy::dofindSXPeaks(const HistogramData::HistogramX &x, const HistogramData::HistogramY &y,
403 Bound low, Bound high, const int workspaceIndex) const {
404 // Get all peaks from the container
405 auto foundPeaks = getAllPeaks(x, y, low, high, m_backgroundStrategy);
406
407 // Convert the found peaks to SXPeaks
408 auto peaks = convertToSXPeaks(x, y, foundPeaks, workspaceIndex);
409
410 return peaks;
411}
412
413std::vector<std::unique_ptr<PeakContainer>>
414AllPeaksStrategy::getAllPeaks(const HistogramData::HistogramX &x, const HistogramData::HistogramY &y, Bound low,
415 Bound high,
416 const Mantid::Crystal::FindSXPeaksHelper::BackgroundStrategy *backgroundStrategy) const {
417 // We iterate over the data and only consider data which is above the
418 // threshold.
419 // Once data starts to be above the threshold we start to record it and add it
420 // to a peak. Once it falls below, it concludes recording of that particular
421 // peak
422 bool isRecording = false;
423
424 std::unique_ptr<PeakContainer> currentPeak = nullptr;
425 std::vector<std::unique_ptr<PeakContainer>> peaks;
426
427 // We want to the upper boundary to be inclusive hence we need to increment it
428 // by one
429 if (high != x.end()) {
430 ++high;
431 }
432 auto distanceMin = std::distance(x.begin(), low);
433 auto distanceMax = std::distance(x.begin(), high);
434
435 const auto lowY = y.begin() + distanceMin;
436 auto highY = distanceMax < static_cast<int>(y.size()) ? y.begin() + distanceMax : y.end();
437
438 for (auto it = lowY; it != highY; ++it) {
439 const auto signal = *it;
440 const auto isAboveThreshold = !backgroundStrategy->isBelowBackground(signal, y);
441
442 // There are four scenarios:
443 // 1. Not recording + below threshold => continue
444 // 2. Not recording + above treshold => start recording
445 // 3. Recording + below threshold => stop recording
446 // 4. Recording + above threshold => continue recording
447 if (!isRecording && (!std::isfinite(signal) || !isAboveThreshold)) {
448 continue;
449 } else if (!isRecording && isAboveThreshold && std::isfinite(signal)) {
450 // only start recording if is finite as NaN values will be found to be above threshold
451 currentPeak = std::make_unique<PeakContainer>(y);
452 currentPeak->startRecord(it);
453 isRecording = true;
454 } else if (isRecording && !isAboveThreshold) {
455 currentPeak->stopRecord(it);
456 peaks.emplace_back(std::move(currentPeak));
457 currentPeak = nullptr;
458 isRecording = false;
459 } else {
460 // this will continue to record NaN if previous point was above the background
461 currentPeak->record(it);
462 }
463 }
464
465 // Handle a peak on the edge if it exists
466 if (isRecording) {
467 if (highY == y.end()) {
468 --highY;
469 }
470 currentPeak->stopRecord(highY);
471 peaks.emplace_back(std::move(currentPeak));
472 currentPeak = nullptr;
473 }
474
475 return peaks;
476}
477
478PeakList AllPeaksStrategy::convertToSXPeaks(const HistogramData::HistogramX &x, const HistogramData::HistogramY &y,
479 const std::vector<std::unique_ptr<PeakContainer>> &foundPeaks,
480 const int workspaceIndex) const {
481 PeakList peaks;
482
483 if (foundPeaks.empty()) {
484 return peaks;
485 }
486
487 // Add a vector to the boost optional
488 peaks = std::vector<FindSXPeaksHelper::SXPeak>();
489
490 for (const auto &peak : foundPeaks) {
491 // Get the index of the bin
492 auto maxY = peak->getMaxIterator();
493 const auto distance = std::distance(y.begin(), maxY);
494 const auto xValue = getXValue(x, distance);
495 const auto tof = convertToTOF(xValue, workspaceIndex);
496 const double phi = calculatePhi(workspaceIndex);
497
498 std::vector<int> specs(1, workspaceIndex);
499 (*peaks).emplace_back(tof, phi, *maxY, specs, workspaceIndex, m_spectrumInfo);
500 }
501
502 return peaks;
503}
504
505/* ------------------------------------------------------------------------------------------
506 * PeakList Reduction Strategy
507 * ------------------------------------------------------------------------------------------
508 */
509
511 : m_compareStrategy(compareStrategy) {}
512
514 : ReducePeakListStrategy(compareStrategy) {}
515
516std::vector<SXPeak> SimpleReduceStrategy::reduce(const std::vector<SXPeak> &peaks,
517 Mantid::Kernel::ProgressBase & /*progress*/) const {
518 // If the peaks are empty then do nothing
519 if (peaks.empty()) {
520 return peaks;
521 }
522
523 std::vector<SXPeak> finalPeaks;
524 for (const auto &currentPeak : peaks) {
525 auto pos = std::find_if(finalPeaks.begin(), finalPeaks.end(), [&currentPeak, this](SXPeak &peak) {
526 auto result = this->m_compareStrategy->compare(currentPeak, peak);
527 // bool result = currentPeak.compare(peak,
528 // resolution);
529 if (result)
530 peak += currentPeak;
531 return result;
532 });
533 if (pos == finalPeaks.end()) {
534 finalPeaks.emplace_back(currentPeak);
535 }
536 }
537
538 return finalPeaks;
539}
540
542 : ReducePeakListStrategy(compareStrategy) {}
543
544std::vector<SXPeak> FindMaxReduceStrategy::reduce(const std::vector<SXPeak> &peaks,
545 Mantid::Kernel::ProgressBase &progress) const {
546 // If the peaks are empty then do nothing
547 if (peaks.empty()) {
548 return peaks;
549 }
550
551 // Groups the peaks into elements which are considered alike
552 auto peakGroups = getPeakGroups(peaks, progress);
553
554 // Now reduce the peaks groups
555 return getFinalPeaks(peakGroups);
556}
557
558// Define some graph elements
559using PeakGraph = adjacency_list<vecS, vecS, undirectedS, SXPeak *>;
560using Vertex = boost::graph_traits<PeakGraph>::vertex_descriptor;
561using Edge = boost::graph_traits<PeakGraph>::edge_descriptor;
562
563std::vector<std::vector<SXPeak *>> FindMaxReduceStrategy::getPeakGroups(const std::vector<SXPeak> &peakList,
564 Mantid::Kernel::ProgressBase &progress) const {
565
566 // Create a vector of addresses. Note that the peaks live on the stack. This
567 // here only works, because the peaks are always in a stack frame below.
568 std::vector<SXPeak *> peaks;
569 peaks.reserve(peakList.size());
570 std::transform(peakList.cbegin(), peakList.cend(), std::back_inserter(peaks),
571 [](const auto &peak) { return &const_cast<SXPeak &>(peak); });
572
573 // Add the peaks to a graph
574 Edge edge;
575 PeakGraph graph;
576 PeakGraph::vertex_iterator vertexIt, vertexEnd;
577
578 // Provide a warning if there are more than 500 peaks found.
579 const size_t numberOfPeaksFound = peaks.size();
580 if (numberOfPeaksFound > 500) {
581 std::string warningMessage = std::string("There are ") + std::to_string(numberOfPeaksFound) +
582 std::string(" peaks being processed. This might take a long time. "
583 "Please check that the cutoff of the background that "
584 "you have selected is high enough, else the algorithm will"
585 " mistake background noise for peaks. The instrument view "
586 "allows you to easily inspect the typical background level.");
587 g_log.warning(warningMessage);
588 }
589
590 std::string message = std::string("There are ") + std::to_string(numberOfPeaksFound) +
591 std::string(" peaks. Investigating peak number ");
592 int peakCounter = 0;
593
594 for (auto peak : peaks) {
595 ++peakCounter;
596
597 // 1. Add the vertex
598 auto vertex = add_vertex(peak, graph);
599
600 // 2. Iterate over all elements already in the graph and check if they need
601 // to edstablish an edge between them.
602 std::tie(vertexIt, vertexEnd) = vertices(graph);
603
604 // Provide a progress report such that users can escape the graph generation
605 if (peakCounter > 50) {
606 progress.doReport(message + std::to_string(peakCounter));
607 }
608
609 for (; vertexIt != vertexEnd; ++vertexIt) {
610 // 2.1 Check if we are looking at the new vertex itself. We don't want
611 // self-loops
612 if (vertex == *vertexIt) {
613 continue;
614 }
615
616 // 2.2 Check if the edge exists already
617 if (boost::edge(vertex, *vertexIt, graph).second) {
618 continue;
619 }
620
621 // 2.3 Check if the two vertices should have an edge
622 const auto toCheck = graph[*vertexIt];
623 if (m_compareStrategy->compare(*peak, *toCheck)) {
624 // We need to create an edge
625 add_edge(vertex, *vertexIt, graph);
626 }
627 }
628 }
629
630 // Create disjoined graphs from graph above
631 std::vector<int> components(boost::num_vertices(graph));
632 const int numberOfPeaks = connected_components(graph, &components[0]);
633
634 std::vector<std::vector<SXPeak *>> peakGroups(numberOfPeaks);
635
636 for (auto i = 0u; i < components.size(); ++i) {
637 auto index = components[i];
638 peakGroups[index].emplace_back(graph[i]);
639 }
640
641 return peakGroups;
642}
643
644std::vector<SXPeak> FindMaxReduceStrategy::getFinalPeaks(const std::vector<std::vector<SXPeak *>> &peakGroups) const {
645 std::vector<SXPeak> peaks;
646 // For each peak groupf find one peak
647 // Currently we select the peak with the largest signal (this strategy could
648 // be changed to something like a weighted mean or similar)
649 for (auto &group : peakGroups) {
650 SXPeak *maxPeak = nullptr;
651 double maxIntensity = std::numeric_limits<double>::min();
652 for (auto *element : group) {
653 if (element->getIntensity() > maxIntensity) {
654 maxIntensity = element->getIntensity();
655 maxPeak = element;
656 }
657 }
658
659 // Add the max peak if valid
660 if (maxPeak) {
661 // check not null as is case when intensity is NaN (though this shouldn't occur now)
662 peaks.emplace_back(*maxPeak);
663 }
664 }
665 return peaks;
666}
667
668/* ------------------------------------------------------------------------------------------
669 * Comparison Strategy
670 * ------------------------------------------------------------------------------------------
671 */
672RelativeCompareStrategy::RelativeCompareStrategy(const double resolution) : m_resolution(resolution) {}
673
674bool RelativeCompareStrategy::compare(const SXPeak &lhs, const SXPeak &rhs) const {
675 return lhs.compare(rhs, m_resolution);
676}
677
678AbsoluteCompareStrategy::AbsoluteCompareStrategy(const double xUnitResolution, const double phiResolution,
679 const double twoThetaResolution, const XAxisUnit units)
680 : m_xUnitResolution(xUnitResolution), m_phiResolution(phiResolution), m_twoThetaResolution(twoThetaResolution),
681 m_units(units) {
682 // Convert the input from degree to radians
683 constexpr double rad2deg = M_PI / 180.;
686}
687
688bool AbsoluteCompareStrategy::compare(const SXPeak &lhs, const SXPeak &rhs) const {
690}
691
692} // namespace Mantid::Crystal::FindSXPeaksHelper
const std::vector< double > & rhs
double maxValue
double minValue
const double m_phi
Definition: LoadEMU.cpp:227
double tolerance
API::SpectrumInfo is an intermediate step towards a SpectrumInfo that is part of Instrument-2....
Definition: SpectrumInfo.h:53
Kernel::V3D sourcePosition() const
Returns the source position.
Kernel::V3D samplePosition() const
Returns the sample position.
bool hasDetectors(const size_t index) const
Returns true if the spectrum is associated with detectors in the instrument.
double twoTheta(const size_t index) const
Returns the scattering angle 2 theta in radians (angle w.r.t.
Kernel::V3D position(const size_t index) const
Returns the position of the spectrum with given index.
double l2(const size_t index) const
Returns L2 (distance from sample to spectrum).
const SpectrumDefinition & spectrumDefinition(const size_t index) const
Returns a const reference to the SpectrumDefinition of the spectrum.
const Geometry::IDetector & detector(const size_t index) const
Return a const reference to the detector or detector group of the spectrum with given index.
double l1() const
Returns L1 (distance from source to sample).
void getDetectorValues(const Kernel::Unit &inputUnit, const Kernel::Unit &outputUnit, const Kernel::DeltaEMode::Type emode, const bool signedTheta, int64_t wsIndex, Kernel::UnitParametersMap &pmap) const
Get the detector values relevant to unit conversion for a workspace index.
AbsoluteCompareStrategy(const double tofResolution, const double phiResolution, const double twoThetaResolution, const XAxisUnit units=XAxisUnit::TOF)
bool compare(const SXPeak &lhs, const SXPeak &rhs) const override
std::vector< std::unique_ptr< PeakContainer > > getAllPeaks(const HistogramData::HistogramX &x, const HistogramData::HistogramY &y, Bound low, Bound high, const Mantid::Crystal::FindSXPeaksHelper::BackgroundStrategy *backgroundStrategy) const
PeakList dofindSXPeaks(const HistogramData::HistogramX &x, const HistogramData::HistogramY &y, Bound low, Bound high, const int workspaceIndex) const override
AllPeaksStrategy(const BackgroundStrategy *backgroundStrategy, const API::SpectrumInfo &spectrumInfo, const double minValue=EMPTY_DBL(), const double maxValue=EMPTY_DBL(), const XAxisUnit units=XAxisUnit::TOF)
PeakList convertToSXPeaks(const HistogramData::HistogramX &x, const HistogramData::HistogramY &y, const std::vector< std::unique_ptr< PeakContainer > > &peaks, const int workspaceIndex) const
virtual bool compare(const SXPeak &lhs, const SXPeak &rhs) const =0
FindMaxReduceStrategy(const CompareStrategy *compareStrategy)
std::vector< SXPeak > reduce(const std::vector< SXPeak > &peaks, Mantid::Kernel::ProgressBase &progress) const override
std::vector< std::vector< SXPeak * > > getPeakGroups(const std::vector< SXPeak > &peakList, Mantid::Kernel::ProgressBase &progress) const
std::vector< SXPeak > getFinalPeaks(const std::vector< std::vector< SXPeak * > > &peakGroups) const
PeakContainer(const HistogramData::HistogramY &y)
double calculatePhi(size_t workspaceIndex) const
Calculates the average phi value if the workspace contains multiple detectors per spectrum,...
virtual PeakList dofindSXPeaks(const HistogramData::HistogramX &x, const HistogramData::HistogramY &y, Bound low, Bound high, const int workspaceIndex) const =0
double getXValue(const HistogramData::HistogramX &x, const size_t peakLocation) const
PeakFindingStrategy(const BackgroundStrategy *backgroundStrategy, const API::SpectrumInfo &spectrumInfo, const double minValue=EMPTY_DBL(), const double maxValue=EMPTY_DBL(), const XAxisUnit units=XAxisUnit::TOF)
BoundsIterator getBounds(const HistogramData::HistogramX &x) const
PeakList findSXPeaks(const HistogramData::HistogramX &x, const HistogramData::HistogramY &y, const int workspaceIndex) const
double convertToTOF(const double xValue, const size_t workspaceIndex) const
ReducePeakListStrategy(const CompareStrategy *compareStrategy)
bool compare(const SXPeak &lhs, const SXPeak &rhs) const override
std::vector< SXPeak > reduce(const std::vector< SXPeak > &peaks, Mantid::Kernel::ProgressBase &progress) const override
SimpleReduceStrategy(const CompareStrategy *compareStrategy)
PeakList dofindSXPeaks(const HistogramData::HistogramX &x, const HistogramData::HistogramY &y, Bound low, Bound high, const int workspaceIndex) const override
StrongestPeaksStrategy(const BackgroundStrategy *backgroundStrategy, const API::SpectrumInfo &spectrumInfo, const double minValue=EMPTY_DBL(), const double maxValue=EMPTY_DBL(), const XAxisUnit units=XAxisUnit::TOF)
Holds a collection of detectors.
Definition: DetectorGroup.h:28
virtual double getPhi() const =0
Gives the phi of this detector object in radians.
virtual detid_t getID() const =0
Get the detector ID.
The Logger class is in charge of the publishing messages from the framework through various channels.
Definition: Logger.h:52
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
virtual void doReport(const std::string &msg="")=0
Pure virtual method that does the progress reporting, to be overridden.
Manage the lifetime of a class intended to be a singleton.
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
Time of flight in microseconds.
Definition: Unit.h:283
Class for 3D vectors.
Definition: V3D.h:34
Mantid::HistogramData::HistogramY::const_iterator yIt
boost::optional< std::vector< SXPeak > > PeakList
boost::graph_traits< PeakGraph >::vertex_descriptor Vertex
adjacency_list< vecS, vecS, undirectedS, SXPeak * > PeakGraph
Mantid::Kernel::Logger g_log("FindSXPeaksHelper")
boost::graph_traits< PeakGraph >::edge_descriptor Edge
HistogramData::HistogramX::const_iterator Bound
XAxisUnit
enum to determine the units of the workspaces X axis we are searching in
std::pair< Bound, Bound > BoundsIterator
constexpr double rad2deg
std::unordered_map< UnitParams, double > UnitParametersMap
Definition: Unit.h:30
MANTID_KERNEL_DLL V3D normalize(V3D v)
Normalizes a V3D.
Definition: V3D.h:341
static constexpr double NeutronMass
Mass of the neutron in kg.
static constexpr double h_bar
Planck constant in J*s, divided by 2 PI.
int32_t detid_t
Typedef for a detector ID.
Definition: SpectrumInfo.h:21
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
Definition: EmptyValues.h:43
Definition: NDArray.h:49
Generate a tableworkspace to store the calibration results.
std::string to_string(const wide_integer< Bits, Signed > &n)
bool isBelowBackground(const double intensity, const HistogramData::HistogramY &) const override
virtual bool isBelowBackground(const double intensity, const HistogramData::HistogramY &y) const =0
bool isBelowBackground(const double intensity, const HistogramData::HistogramY &y) const override
int m_nPixels
Number of contributing pixels.
detid_t getDetectorId() const
Getter for the detector Id.
double m_phi
Phi angle for the centre detector of the peak.
double m_tof
TOF for the peak centre.
std::vector< int > m_spectra
Contributing spectra to this peak.
Kernel::V3D m_unitWaveVector
Unit vector in the direction of the wavevector.
double m_LTotal
Detector-sample distance.
SXPeak(double t, double phi, double intensity, const std::vector< int > &spectral, const size_t wsIndex, const Mantid::API::SpectrumInfo &spectrumInfo)
Constructor.
bool compare(const SXPeak &rhs, double tolerance) const
Object comparison.
double m_intensity
Measured intensity of centre of the peak.
double m_dSpacing
d-spacing at the peak centre
void reduce()
Normalise by number of pixels.
SXPeak & operator+=(const SXPeak &rhs)
Operator addition overload.
Mantid::Kernel::V3D getQ() const
Getter for LabQ.
const double & getIntensity() const
Getter for the intensity.
double m_twoTheta
2 * theta angle for then centre detector of the peak
size_t m_wsIndex
Detector workspace index.
Peak indexing algorithm, which works by assigning multiple possible HKL values to each peak and then ...
Definition: IndexSXPeaks.h:29