Mantid
Loading...
Searching...
No Matches
FindSXPeaks.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#include "MantidAPI/Axis.h"
10#include "MantidAPI/Run.h"
13#include "MantidIndexing/IndexInfo.h"
18#include "MantidKernel/Unit.h"
20
21#include <vector>
22
23using namespace Mantid::DataObjects;
24using namespace Mantid::API;
26
27namespace Mantid::Crystal {
28
29const std::string FindSXPeaks::strongestPeakStrategy = "StrongestPeakOnly";
30const std::string FindSXPeaks::allPeaksStrategy = "AllPeaks";
31
32const std::string FindSXPeaks::relativeResolutionStrategy = "RelativeResolution";
33const std::string FindSXPeaks::absoluteResolutionPeaksStrategy = "AbsoluteResolution";
34
35// Register the class into the algorithm factory
37
38using namespace Kernel;
39using namespace API;
41
43 : API::Algorithm(), m_MinRange(DBL_MAX), m_MaxRange(-DBL_MAX), m_MinWsIndex(0), m_MaxWsIndex(0) {}
44
49 auto wsValidation = std::make_shared<CompositeValidator>();
50 wsValidation->add<HistogramValidator>();
51
52 auto unitValidation = std::make_shared<CompositeValidator>(CompositeRelation::OR);
53 unitValidation->add<WorkspaceUnitValidator>("TOF");
54 unitValidation->add<WorkspaceUnitValidator>("dSpacing");
55
56 wsValidation->add(unitValidation);
57
58 declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input, wsValidation),
59 "The name of the Workspace2D to take as input");
60 declareProperty("RangeLower", EMPTY_DBL(), "The X value to search from (default 0)");
61 declareProperty("RangeUpper", EMPTY_DBL(), "The X value to search to (default total number of bins)");
62 auto mustBePositive = std::make_shared<BoundedValidator<int>>();
63 mustBePositive->setLower(0);
64 declareProperty("StartWorkspaceIndex", 0, mustBePositive, "Start workspace index (default 0)");
65 declareProperty("EndWorkspaceIndex", EMPTY_INT(), mustBePositive,
66 "End workspace index (default to total number of histograms)");
67
68 // ---------------------------------------------------------------
69 // Peak strategies + Threshold
70 // ---------------------------------------------------------------
71 auto mustBePositiveDouble = std::make_shared<BoundedValidator<double>>();
72 mustBePositiveDouble->setLower(0.0);
73
74 std::vector<std::string> peakFindingStrategy = {strongestPeakStrategy, allPeaksStrategy};
75 declareProperty("PeakFindingStrategy", strongestPeakStrategy,
76 std::make_shared<StringListValidator>(peakFindingStrategy),
77 "Different options for peak finding."
78 "1. StrongestPeakOnly: Looks only for the the strongest peak in each "
79 "spectrum (provided there is "
80 "one). This options is more performant than the AllPeaks option.\n"
81 "2. AllPeaks: This strategy will find all peaks in each "
82 "spectrum. This is slower than StrongestPeakOnly. Note that the "
83 "recommended ResolutionStrategy in this mode is AbsoluteResolution.\n");
84
85 // Declare
86 declareProperty("SignalBackground", 10.0, mustBePositiveDouble,
87 "Multiplication factor for the signal background. Peaks which are"
88 " below the estimated background are discarded. The background is "
89 "estimated"
90 " to be an average of the first and the last signal and multiplied"
91 " by the SignalBackground property.\n");
92
93 declareProperty("AbsoluteBackground", 30.0, mustBePositiveDouble,
94 "Peaks which are below the specified absolute background are discarded."
95 " The background is gloabally specified for all spectra. Inspect your "
96 "data in the InstrumentView to get a good feeling for the background "
97 "threshold.\n"
98 "Background thresholds which are too low will mistake noise for peaks.");
99
100 // Enable
101 setPropertySettings("SignalBackground", std::make_unique<EnabledWhenProperty>(
104
105 setPropertySettings("AbsoluteBackground",
106 std::make_unique<EnabledWhenProperty>(
108
109 // Group
110 const std::string peakGroup = "Peak Finding Settings";
111 setPropertyGroup("PeakFindingStrategy", peakGroup);
112 setPropertyGroup("SignalBackground", peakGroup);
113 setPropertyGroup("AbsoluteBackground", peakGroup);
114
115 // ---------------------------------------------------------------
116 // Resolution
117 // ---------------------------------------------------------------
118 // Declare
119 std::vector<std::string> resolutionStrategy = {relativeResolutionStrategy, absoluteResolutionPeaksStrategy};
120 declareProperty("ResolutionStrategy", relativeResolutionStrategy,
121 std::make_shared<StringListValidator>(resolutionStrategy),
122 "Different options for the resolution."
123 "1. RelativeResolution: This defines a relative tolerance "
124 "needed to avoid peak duplication in number of pixels. "
125 "This selection will enable the Resolution property and "
126 "disable the XResolution, PhiResolution, ThetaResolution.\n"
127 "1. AbsoluteResolution: This defines an absolute tolerance "
128 "needed to avoid peak duplication in number of pixels. "
129 "This selection will disable the Resolution property and "
130 "enable the XResolution, PhiResolution, "
131 "ThetaResolution.\n");
132
133 declareProperty("Resolution", 0.01, mustBePositiveDouble,
134 "Tolerance needed to avoid peak duplication in number of pixels");
135
136 declareProperty("XResolution", 0., mustBePositiveDouble,
137 "Absolute tolerance in time-of-flight or d-spacing needed to avoid peak "
138 "duplication in number of pixels. The values are specified "
139 "in either microseconds or angstroms.");
140
141 declareProperty("PhiResolution", 1., mustBePositiveDouble,
142 "Absolute tolerance in the phi "
143 "coordinate needed to avoid peak "
144 "duplication in number of pixels. The "
145 "values are specified in degrees.");
146
147 declareProperty("TwoThetaResolution", 1., mustBePositiveDouble,
148 "Absolute tolerance of two theta value needed to avoid peak "
149 "duplication in number of pixels. The values are specified "
150 "in degrees.");
151
152 // Enable
153 setPropertySettings("Resolution", std::make_unique<EnabledWhenProperty>(
156
157 setPropertySettings("XResolution", std::make_unique<EnabledWhenProperty>(
160
161 setPropertySettings("PhiResolution", std::make_unique<EnabledWhenProperty>(
164
165 setPropertySettings("TwoThetaResolution", std::make_unique<EnabledWhenProperty>(
168
169 declareProperty(std::make_unique<WorkspaceProperty<PeaksWorkspace>>("OutputWorkspace", "", Direction::Output),
170 "The name of the PeaksWorkspace in which to store the list "
171 "of peaks found");
172
173 // Group
174 const std::string resolutionGroup = "Resolution Settings";
175 setPropertyGroup("ResolutionStrategy", resolutionGroup);
176 setPropertyGroup("Resolution", resolutionGroup);
177 setPropertyGroup("XResolution", resolutionGroup);
178 setPropertyGroup("PhiResolution", resolutionGroup);
179 setPropertyGroup("TwoThetaResolution", resolutionGroup);
180
181 // Create the output peaks workspace
182 m_peaks.reset(new PeaksWorkspace);
183}
184
185/*
186 * Validate the input parameters
187 * @returns map with keys corresponding to properties with errors and values
188 * containing the error messages.
189 */
190std::map<std::string, std::string> FindSXPeaks::validateInputs() {
191 // create the map
192 std::map<std::string, std::string> validationOutput;
193 const std::string resolutionStrategy = getProperty("ResolutionStrategy");
194 const auto xResolutionProperty = getPointerToProperty("XResolution");
195
196 // Check that the user has set a valid value for the x resolution when
197 // in absolute resolution mode.
198 if (resolutionStrategy == FindSXPeaks::absoluteResolutionPeaksStrategy && xResolutionProperty->isDefault()) {
199 validationOutput["XResolution"] = "XResolution must be set to a value greater than 0";
200 }
201
202 return validationOutput;
203}
204
210 // Try and retrieve the optional properties
211 m_MinRange = getProperty("RangeLower");
212 m_MaxRange = getProperty("RangeUpper");
213
214 // the assignment below is intended and if removed will break the unit tests
215 m_MinWsIndex = static_cast<int>(getProperty("StartWorkspaceIndex"));
216 m_MaxWsIndex = static_cast<int>(getProperty("EndWorkspaceIndex"));
217
218 // Get the input workspace
219 MatrixWorkspace_const_sptr localworkspace = getProperty("InputWorkspace");
220
221 // copy the instrument across. Cannot generate peaks without doing this
222 // first.
223 m_peaks->setInstrument(localworkspace->getInstrument());
224
225 size_t numberOfSpectra = localworkspace->getNumberHistograms();
226
227 // Check 'StartSpectrum' is in range 0-numberOfSpectra
228 if (m_MinWsIndex > numberOfSpectra) {
229 g_log.warning("StartSpectrum out of range! Set to 0.");
230 m_MinWsIndex = 0;
231 }
233 throw std::invalid_argument("Cannot have StartWorkspaceIndex > EndWorkspaceIndex");
234 }
236 m_MaxWsIndex = numberOfSpectra - 1;
237 if (m_MaxWsIndex > numberOfSpectra - 1 || m_MaxWsIndex < m_MinWsIndex) {
238 g_log.warning("EndSpectrum out of range! Set to max detector number");
239 m_MaxWsIndex = numberOfSpectra;
240 }
241 if (m_MinRange > m_MaxRange) {
242 g_log.warning("Range_upper is less than Range_lower. Will integrate up to "
243 "frame maximum.");
244 m_MaxRange = 0.0;
245 }
246
247 Progress progress(this, 0.0, 1.0, m_MaxWsIndex - m_MinWsIndex + 2);
248
249 // Calculate the primary flight path.
250 const auto &spectrumInfo = localworkspace->spectrumInfo();
251
252 // Get the background strategy
253 auto backgroundStrategy = getBackgroundStrategy();
254
255 // Get the peak finding strategy
256 const auto xUnit = getWorkspaceXAxisUnit(localworkspace);
257 auto peakFindingStrategy =
258 getPeakFindingStrategy(backgroundStrategy.get(), spectrumInfo, m_MinRange, m_MaxRange, xUnit);
259
260 peakvector entries;
261 entries.reserve(m_MaxWsIndex - m_MinWsIndex);
262 // Count the peaks so that we can resize the peak vector at the end.
263 PARALLEL_FOR_IF(Kernel::threadSafe(*localworkspace))
264 for (auto wsIndex = static_cast<int>(m_MinWsIndex); wsIndex <= static_cast<int>(m_MaxWsIndex); ++wsIndex) {
266
267 // If no detector found / monitor, skip onto the next spectrum
268 const auto wsIndexSize_t = static_cast<size_t>(wsIndex);
269 if (!spectrumInfo.hasDetectors(wsIndexSize_t) || spectrumInfo.isMonitor(wsIndexSize_t)) {
270 continue;
271 }
272
273 // Retrieve the spectrum into a vector
274 const auto &x = localworkspace->x(wsIndex);
275 const auto &y = localworkspace->y(wsIndex);
276
277 // Run the peak finding strategy
278 auto foundPeaks = peakFindingStrategy->findSXPeaks(x, y, wsIndex);
279 if (!foundPeaks) {
280 continue;
281 }
282
283 PARALLEL_CRITICAL(entries) { std::copy(foundPeaks->cbegin(), foundPeaks->cend(), std::back_inserter(entries)); }
284 progress.report();
286 }
288
289 // Now reduce the list with duplicate entries
290 reducePeakList(entries, progress);
291
292 setProperty("OutputWorkspace", m_peaks);
293 progress.report();
294}
295
303 MatrixWorkspace_const_sptr localworkspace = getProperty("InputWorkspace");
304 auto &goniometerMatrix = localworkspace->run().getGoniometer().getR();
305 auto compareStrategy = getCompareStrategy();
306 auto reductionStrategy = getReducePeakListStrategy(compareStrategy.get());
307 auto finalv = reductionStrategy->reduce(pcv, progress);
308
309 for (auto &finalPeak : finalv) {
310 finalPeak.reduce();
311 try {
312 IPeak_uptr ipeak = m_peaks->createPeak(finalPeak.getQ());
313 Peak_uptr peak(static_cast<Peak *>(ipeak.release()));
314 if (peak) {
315 peak->setIntensity(finalPeak.getIntensity());
316 peak->setDetectorID(finalPeak.getDetectorId());
317 peak->setGoniometerMatrix(goniometerMatrix);
318 peak->setRunNumber(localworkspace->getRunNumber());
319 m_peaks->addPeak(*peak);
320 }
321 } catch (std::exception &e) {
322 g_log.error() << e.what() << '\n';
323 }
324 }
325}
326
336 const auto xAxis = workspace->getAxis(0);
337 const auto unitID = xAxis->unit()->unitID();
338
339 if (unitID == "TOF") {
340 return XAxisUnit::TOF;
341 } else {
342 return XAxisUnit::DSPACING;
343 }
344}
345
346std::unique_ptr<BackgroundStrategy> FindSXPeaks::getBackgroundStrategy() const {
347 const std::string peakFindingStrategy = getProperty("PeakFindingStrategy");
348 if (peakFindingStrategy == strongestPeakStrategy) {
349 const double signalBackground = getProperty("SignalBackground");
350 return std::make_unique<PerSpectrumBackgroundStrategy>(signalBackground);
351 } else if (peakFindingStrategy == allPeaksStrategy) {
352 const double background = getProperty("AbsoluteBackground");
353 return std::make_unique<AbsoluteBackgroundStrategy>(background);
354 } else {
355 throw std::invalid_argument("The selected background strategy has not been implemented yet.");
356 }
357}
358
359std::unique_ptr<FindSXPeaksHelper::PeakFindingStrategy>
360FindSXPeaks::getPeakFindingStrategy(const BackgroundStrategy *backgroundStrategy, const API::SpectrumInfo &spectrumInfo,
361 const double minValue, const double maxValue, const XAxisUnit tofUnits) const {
362 // Get the peak finding stratgy
363 std::string peakFindingStrategy = getProperty("PeakFindingStrategy");
364 if (peakFindingStrategy == strongestPeakStrategy) {
365 return std::make_unique<StrongestPeaksStrategy>(backgroundStrategy, spectrumInfo, minValue, maxValue, tofUnits);
366 } else if (peakFindingStrategy == allPeaksStrategy) {
367 return std::make_unique<AllPeaksStrategy>(backgroundStrategy, spectrumInfo, minValue, maxValue, tofUnits);
368 } else {
369 throw std::invalid_argument("The selected peak finding strategy has not been implemented yet.");
370 }
371}
372
373std::unique_ptr<FindSXPeaksHelper::ReducePeakListStrategy>
375 const std::string peakFindingStrategy = getProperty("PeakFindingStrategy");
376 auto useSimpleReduceStrategy = peakFindingStrategy == strongestPeakStrategy;
377 if (useSimpleReduceStrategy) {
378 return std::make_unique<FindSXPeaksHelper::SimpleReduceStrategy>(compareStrategy);
379 } else {
380 return std::make_unique<FindSXPeaksHelper::FindMaxReduceStrategy>(compareStrategy);
381 }
382}
383
384std::unique_ptr<FindSXPeaksHelper::CompareStrategy> FindSXPeaks::getCompareStrategy() const {
385 const std::string resolutionStrategy = getProperty("ResolutionStrategy");
386 auto useRelativeResolutionStrategy = resolutionStrategy == relativeResolutionStrategy;
387 if (useRelativeResolutionStrategy) {
388 double resolution = getProperty("Resolution");
389 return std::make_unique<FindSXPeaksHelper::RelativeCompareStrategy>(resolution);
390 } else {
391 double xUnitResolution = getProperty("XResolution");
392 double phiResolution = getProperty("PhiResolution");
393 double twoThetaResolution = getProperty("TwoThetaResolution");
394 const auto tofUnits = getWorkspaceXAxisUnit(getProperty("InputWorkspace"));
395 return std::make_unique<FindSXPeaksHelper::AbsoluteCompareStrategy>(xUnitResolution, phiResolution,
396 twoThetaResolution, tofUnits);
397 }
398}
399
400} // namespace Mantid::Crystal
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
double maxValue
double minValue
IPeaksWorkspace_sptr workspace
Definition: IndexPeaks.cpp:114
#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.
Base class from which all concrete algorithm classes should be derived.
Definition: Algorithm.h:85
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
Kernel::Property * getPointerToProperty(const std::string &name) const override
Get a property by name.
Definition: Algorithm.cpp:2033
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
void progress(double p, const std::string &msg="", double estimatedTime=0.0, int progressPrecision=0)
Sends ProgressNotification.
Definition: Algorithm.cpp:231
static bool isEmpty(const NumT toCheck)
checks that the value was not set by users, uses the value in empty double/int.
A validator which checks that a workspace contains histogram data (the default) or point data as requ...
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
API::SpectrumInfo is an intermediate step towards a SpectrumInfo that is part of Instrument-2....
Definition: SpectrumInfo.h:53
A property class for workspaces.
A validator which checks that the unit of the workspace referred to by a WorkspaceProperty is the exp...
BackgroundStrategy : Abstract class used for identifying elements of a IMDWorkspace that are not cons...
Search detector space for single crystal peaks.
Definition: FindSXPeaks.h:38
FindSXPeaksHelper::XAxisUnit getWorkspaceXAxisUnit(const Mantid::API::MatrixWorkspace_const_sptr &workspace) const
Check what x units this workspace has.
static const std::string absoluteResolutionPeaksStrategy
Definition: FindSXPeaks.h:60
Mantid::DataObjects::PeaksWorkspace_sptr m_peaks
Definition: FindSXPeaks.h:99
size_t m_MaxWsIndex
The spectrum to finish the integration at.
Definition: FindSXPeaks.h:97
void exec() override
Executes the algorithm.
std::unique_ptr< FindSXPeaksHelper::BackgroundStrategy > getBackgroundStrategy() const
Selects a background strategy.
FindSXPeaks()
Default constructor.
Definition: FindSXPeaks.cpp:42
std::unique_ptr< FindSXPeaksHelper::CompareStrategy > getCompareStrategy() const
Selects a comparison strategy.
std::map< std::string, std::string > validateInputs() override
Perform validation of ALL the input properties of the algorithm.
double m_MinRange
The value in X to start the search from.
Definition: FindSXPeaks.h:91
double m_MaxRange
The value in X to finish the search at.
Definition: FindSXPeaks.h:93
static const std::string strongestPeakStrategy
Definition: FindSXPeaks.h:57
static const std::string relativeResolutionStrategy
Definition: FindSXPeaks.h:59
static const std::string allPeaksStrategy
Definition: FindSXPeaks.h:58
void init() override
Initialisation method.
Definition: FindSXPeaks.cpp:48
void reducePeakList(const peakvector &, Mantid::API::Progress &progress)
Reduce the peak list by removing duplicates then convert SXPeaks objects to PeakObjects and add them ...
size_t m_MinWsIndex
The spectrum to start the integration from.
Definition: FindSXPeaks.h:95
std::unique_ptr< FindSXPeaksHelper::PeakFindingStrategy > getPeakFindingStrategy(const FindSXPeaksHelper::BackgroundStrategy *backgroundStrategy, const API::SpectrumInfo &spectrumInfo, const double minValue, const double maxValue, const FindSXPeaksHelper::XAxisUnit tofUnits=FindSXPeaksHelper::XAxisUnit::TOF) const
Selects a peak finding strategy.
std::unique_ptr< FindSXPeaksHelper::ReducePeakListStrategy > getReducePeakListStrategy(const FindSXPeaksHelper::CompareStrategy *compareStrategy) const
Selects a peak finding strategy.
Structure describing a single-crystal peak.
Definition: Peak.h:34
The class PeaksWorkspace stores information about a set of SCD peaks.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void setPropertySettings(const std::string &name, std::unique_ptr< IPropertySettings > settings)
void setPropertyGroup(const std::string &name, const std::string &group)
Set the group for a given property.
void error(const std::string &msg)
Logs at error level.
Definition: Logger.cpp:77
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
std::shared_ptr< const MatrixWorkspace > MatrixWorkspace_const_sptr
shared pointer to the matrix workspace base class (const version)
XAxisUnit
enum to determine the units of the workspaces X axis we are searching in
std::vector< FindSXPeaksHelper::SXPeak > peakvector
Definition: FindSXPeaks.h:23
std::unique_ptr< Peak > Peak_uptr
Definition: Peak.h:170
std::unique_ptr< IPeak > IPeak_uptr
Definition: IPeak.h:103
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
constexpr int EMPTY_INT() noexcept
Returns what we consider an "empty" integer within a property.
Definition: EmptyValues.h:25
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