16#include "MantidTypes/SpectrumDefinition.h"
18#define BOOST_ALLOW_DEPRECATED_HEADERS
19#include <boost/graph/adjacency_list.hpp>
20#undef BOOST_ALLOW_DEPRECATED_HEADERS
22#include <boost/graph/connected_components.hpp>
27const double TWO_PI = 2 * M_PI;
29bool isDifferenceLargerThanTolerance(
const double angle1,
const double angle2,
const double tolerance) {
30 auto difference = std::abs(angle1 - angle2);
34 if (difference > TWO_PI) {
35 difference = std::fmod(difference, TWO_PI);
39 if (difference > M_PI) {
40 difference = TWO_PI - difference;
68SXPeak::SXPeak(
double t,
double phi,
double intensity,
const std::vector<int> &spectral,
const size_t wsIndex,
70 : m_tof(t),
m_phi(phi), m_intensity(intensity), m_spectra(spectral), m_wsIndex(wsIndex) {
73 throw std::invalid_argument(
"SXPeak: Cannot have an intensity < 0");
75 if (spectral.empty()) {
76 throw std::invalid_argument(
"SXPeak: Cannot have zero sized spectral list");
79 throw std::invalid_argument(
"SXPeak: Spectrum at ws index " +
std::to_string(wsIndex) +
" doesn't have detectors");
82 const auto l1 = spectrumInfo.
l1();
88 throw std::invalid_argument(
"SXPeak: Cannot have detector distance < 0");
97 unit->initialize(l1, 0,
pmap);
100 }
catch (std::exception &) {
108 const auto beamDir =
normalize(samplePos - sourcePos);
110 const auto detDir =
normalize(detPos - samplePos);
136 if (std::abs(x_1 - x_2) > xTolerance) {
140 if (isDifferenceLargerThanTolerance(
m_phi,
rhs.m_phi, phiTolerance)) {
143 if (isDifferenceLargerThanTolerance(
m_twoTheta,
rhs.m_twoTheta, thetaTolerance)) {
164 double wvi = 1.0 / wi;
205 : m_y(
y), m_startIndex(0), m_stopIndex(m_y.size() - 1), m_maxIndex(0) {}
249 const HistogramData::HistogramY & )
const {
254 : m_backgroundMultiplier(backgroundMultiplier) {}
257 const HistogramData::HistogramY &
y)
const {
258 auto background = 0.5 * (1.0 +
y.front() +
y.back());
260 return intensity < background;
271 : m_backgroundStrategy(backgroundStrategy), m_minValue(
minValue), m_maxValue(
maxValue),
272 m_spectrumInfo(spectrumInfo), m_units(units) {}
275 const int workspaceIndex)
const {
280 auto lowit = boundsIterator.first;
281 auto highit = boundsIterator.second;
284 if (lowit ==
x.end() || highit ==
x.begin()) {
300 using std::placeholders::_1;
303 : std::find_if(lowit,
x.end(), std::bind(std::greater<double>(), _1,
m_maxValue));
305 return std::make_pair(lowit, highit);
320 const auto numberOfDetectors = spectrumDefinition.size();
322 if (numberOfDetectors == 1) {
327 if (!detectorGroup) {
328 throw std::runtime_error(
"Could not cast to detector group");
330 phi = detectorGroup->getPhi();
340 auto leftBinPosition =
x.begin() + peakLocation;
341 const double leftBinEdge = *leftBinPosition;
342 const double rightBinEdge = *std::next(leftBinPosition);
343 return 0.5 * (leftBinEdge + rightBinEdge);
357 return unit->singleToTOF(xValue);
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);
372 auto maxY = (
y.size() > 1) ? std::max_element(
y.begin() + distmin,
y.begin() + distmax) :
y.begin();
375 double intensity = (*maxY);
381 const auto distance = std::distance(
y.begin(), maxY);
386 std::vector<int> specs(1, workspaceIndex);
387 std::vector<SXPeak> peaks;
388 peaks.emplace_back(tof, phi, *maxY, specs, workspaceIndex,
m_spectrumInfo);
397 throw std::invalid_argument(
"The AllPeaksStrategy has to be initialized "
398 "with the AbsoluteBackgroundStrategy.");
403 Bound low,
Bound high,
const int workspaceIndex)
const {
413std::vector<std::unique_ptr<PeakContainer>>
422 bool isRecording =
false;
424 std::unique_ptr<PeakContainer> currentPeak =
nullptr;
425 std::vector<std::unique_ptr<PeakContainer>> peaks;
429 if (high !=
x.end()) {
432 auto distanceMin = std::distance(
x.begin(), low);
433 auto distanceMax = std::distance(
x.begin(), high);
435 const auto lowY =
y.begin() + distanceMin;
436 auto highY = distanceMax < static_cast<int>(
y.size()) ?
y.begin() + distanceMax :
y.end();
438 for (
auto it = lowY; it != highY; ++it) {
439 const auto signal = *it;
447 if (!isRecording && (!std::isfinite(signal) || !isAboveThreshold)) {
449 }
else if (!isRecording && isAboveThreshold && std::isfinite(signal)) {
451 currentPeak = std::make_unique<PeakContainer>(
y);
452 currentPeak->startRecord(it);
454 }
else if (isRecording && !isAboveThreshold) {
455 currentPeak->stopRecord(it);
456 peaks.emplace_back(std::move(currentPeak));
457 currentPeak =
nullptr;
461 currentPeak->record(it);
467 if (highY ==
y.end()) {
470 currentPeak->stopRecord(highY);
471 peaks.emplace_back(std::move(currentPeak));
472 currentPeak =
nullptr;
479 const std::vector<std::unique_ptr<PeakContainer>> &foundPeaks,
480 const int workspaceIndex)
const {
483 if (foundPeaks.empty()) {
488 peaks = std::vector<FindSXPeaksHelper::SXPeak>();
490 for (
const auto &peak : foundPeaks) {
492 auto maxY = peak->getMaxIterator();
493 const auto distance = std::distance(
y.begin(), maxY);
498 std::vector<int> specs(1, workspaceIndex);
499 (*peaks).emplace_back(tof, phi, *maxY, specs, workspaceIndex,
m_spectrumInfo);
511 : m_compareStrategy(compareStrategy) {}
523 std::vector<SXPeak> finalPeaks;
524 for (
const auto ¤tPeak : peaks) {
525 auto pos = std::find_if(finalPeaks.begin(), finalPeaks.end(), [¤tPeak,
this](
SXPeak &peak) {
526 auto result = this->m_compareStrategy->compare(currentPeak, peak);
533 if (pos == finalPeaks.end()) {
534 finalPeaks.emplace_back(currentPeak);
559using PeakGraph = adjacency_list<vecS, vecS, undirectedS, SXPeak *>;
560using Vertex = boost::graph_traits<PeakGraph>::vertex_descriptor;
561using Edge = boost::graph_traits<PeakGraph>::edge_descriptor;
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); });
576 PeakGraph::vertex_iterator vertexIt, vertexEnd;
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.");
590 std::string message = std::string(
"There are ") +
std::to_string(numberOfPeaksFound) +
591 std::string(
" peaks. Investigating peak number ");
594 for (
auto peak : peaks) {
598 auto vertex = add_vertex(peak, graph);
602 std::tie(vertexIt, vertexEnd) = vertices(graph);
605 if (peakCounter > 50) {
609 for (; vertexIt != vertexEnd; ++vertexIt) {
612 if (vertex == *vertexIt) {
617 if (boost::edge(vertex, *vertexIt, graph).second) {
622 const auto toCheck = graph[*vertexIt];
625 add_edge(vertex, *vertexIt, graph);
631 std::vector<int> components(boost::num_vertices(graph));
632 const int numberOfPeaks = connected_components(graph, &components[0]);
634 std::vector<std::vector<SXPeak *>> peakGroups(numberOfPeaks);
636 for (
auto i = 0u; i < components.size(); ++i) {
637 auto index = components[i];
638 peakGroups[
index].emplace_back(graph[i]);
645 std::vector<SXPeak> peaks;
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();
662 peaks.emplace_back(*maxPeak);
679 const double twoThetaResolution,
const XAxisUnit units)
680 : m_xUnitResolution(xUnitResolution), m_phiResolution(phiResolution), m_twoThetaResolution(twoThetaResolution),
683 constexpr double rad2deg = M_PI / 180.;
const std::vector< double > & rhs
API::SpectrumInfo is an intermediate step towards a SpectrumInfo that is part of Instrument-2....
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)
double m_twoThetaResolution
bool compare(const SXPeak &lhs, const SXPeak &rhs) const override
const double m_xUnitResolution
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)
void startRecord(yIt item)
yIt getMaxIterator() const
const HistogramData::HistogramY & m_y
void stopRecord(yIt item)
size_t getNumberOfPointsInPeak() const
double calculatePhi(size_t workspaceIndex) const
Calculates the average phi value if the workspace contains multiple detectors per spectrum,...
const BackgroundStrategy * m_backgroundStrategy
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
const API::SpectrumInfo & m_spectrumInfo
const CompareStrategy * m_compareStrategy
ReducePeakListStrategy(const CompareStrategy *compareStrategy)
bool compare(const SXPeak &lhs, const SXPeak &rhs) const override
const double m_resolution
RelativeCompareStrategy(const double resolution)
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.
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.
void warning(const std::string &msg)
Logs at warning level.
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.
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
std::unordered_map< UnitParams, double > UnitParametersMap
MANTID_KERNEL_DLL V3D normalize(V3D v)
Normalizes a V3D.
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.
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
Generate a tableworkspace to store the calibration results.
std::string to_string(const wide_integer< Bits, Signed > &n)
AbsoluteBackgroundStrategy(const double background)
bool isBelowBackground(const double intensity, const HistogramData::HistogramY &) const override
virtual bool isBelowBackground(const double intensity, const HistogramData::HistogramY &y) const =0
PerSpectrumBackgroundStrategy(const double backgroundMultiplier)
double m_backgroundMultiplier
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.
detid_t m_detId
Detector ID.
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.
std::string m_qConvention
Q Convention.
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 ...