-
Notifications
You must be signed in to change notification settings - Fork 122
/
MaxMin.cpp
152 lines (130 loc) · 5.72 KB
/
MaxMin.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
// NScD Oak Ridge National Laboratory, European Spallation Source,
// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidAlgorithms/MaxMin.h"
#include "MantidAPI/HistogramValidator.h"
#include "MantidDataObjects/TableWorkspace.h"
#include "MantidDataObjects/Workspace2D.h"
#include "MantidDataObjects/WorkspaceCreation.h"
#include "MantidHistogramData/Histogram.h"
#include "MantidKernel/BoundedValidator.h"
namespace Mantid {
namespace Algorithms {
// Register the class into the algorithm factory
DECLARE_ALGORITHM(MaxMin)
using namespace Kernel;
using namespace API;
using namespace Mantid::DataObjects;
using namespace Mantid::HistogramData;
/** Initialisation method.
*
*/
void MaxMin::init() {
declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input,
std::make_shared<HistogramValidator>()),
"The name of the Workspace2D to take as input");
declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
"The name of the workspace in which to store the result");
declareProperty("ShowMin", false, "Flag to show minimum instead of maximum (default=false)");
declareProperty("RangeLower", EMPTY_DBL(), "The X value to search from (default min)");
declareProperty("RangeUpper", EMPTY_DBL(), "The X value to search to (default max)");
auto mustBePositive = std::make_shared<BoundedValidator<int>>();
mustBePositive->setLower(0);
declareProperty("StartWorkspaceIndex", 0, mustBePositive, "Start spectrum number (default 0)");
declareProperty("EndWorkspaceIndex", EMPTY_INT(), mustBePositive, "End spectrum number (default max)");
}
/** Executes the algorithm
*
* @throw runtime_error Thrown if algorithm cannot execute
*/
void MaxMin::exec() {
// Try and retrieve the optional properties
/// The value in X to start the search from
double MinRange = getProperty("RangeLower");
/// The value in X to finish the search at
double MaxRange = getProperty("RangeUpper");
/// The spectrum to start the integration from
int MinSpec = getProperty("StartWorkspaceIndex");
/// The spectrum to finish the integration at
int MaxSpec = getProperty("EndWorkspaceIndex");
/// The flag to show minimum
bool showMin = getProperty("ShowMin");
// Get the input workspace
MatrixWorkspace_const_sptr localworkspace = getProperty("InputWorkspace");
const auto numberOfSpectra = static_cast<int>(localworkspace->getNumberHistograms());
// Check 'StartSpectrum' is in range 0-numberOfSpectra
if (MinSpec > numberOfSpectra) {
g_log.warning("StartSpectrum out of range! Set to 0.");
MinSpec = 0;
}
if (isEmpty(MaxSpec))
MaxSpec = numberOfSpectra - 1;
if (MaxSpec > numberOfSpectra - 1 || MaxSpec < MinSpec) {
g_log.warning("EndSpectrum out of range! Set to max detector number");
MaxSpec = numberOfSpectra;
}
if (MinRange > MaxRange) {
g_log.warning("Range_upper is less than Range_lower. Will integrate up to "
"frame maximum.");
MaxRange = 0.0;
}
// Create the 1D workspace for the output
MatrixWorkspace_sptr outputWorkspace;
outputWorkspace = create<HistoWorkspace>(*localworkspace, MaxSpec - MinSpec + 1, BinEdges(2));
Progress progress(this, 0.0, 1.0, (MaxSpec - MinSpec + 1));
PARALLEL_FOR_IF(Kernel::threadSafe(*localworkspace, *outputWorkspace))
// Loop over spectra
for (int i = MinSpec; i <= MaxSpec; ++i) {
PARALLEL_START_INTERUPT_REGION
int newindex = i - MinSpec;
// Copy over spectrum and detector number info
outputWorkspace->getSpectrum(newindex).copyInfoFrom(localworkspace->getSpectrum(i));
// Retrieve the spectrum into a vector
auto &X = localworkspace->x(i);
auto &Y = localworkspace->y(i);
// Find the range [min,max]
MantidVec::const_iterator lowit, highit;
if (MinRange == EMPTY_DBL())
lowit = X.begin();
else
lowit = std::lower_bound(X.begin(), X.end(), MinRange);
if (MaxRange == EMPTY_DBL())
highit = X.end();
else {
using std::placeholders::_1;
highit = std::find_if(lowit, X.end(), std::bind(std::greater<double>(), _1, MaxRange));
}
// If range specified doesn't overlap with this spectrum then bail out
if (lowit == X.end() || highit == X.begin() || lowit == highit)
continue;
--highit; // Upper limit is the bin before, i.e. the last value smaller than
// MaxRange
MantidVec::difference_type distmin = std::distance(X.begin(), lowit);
MantidVec::difference_type distmax = std::distance(X.begin(), highit);
MantidVec::const_iterator maxY;
// Find the max/min element
if (showMin) {
maxY = std::min_element(Y.begin() + distmin, Y.begin() + distmax);
} else {
maxY = std::max_element(Y.begin() + distmin, Y.begin() + distmax);
}
MantidVec::difference_type d = std::distance(Y.begin(), maxY);
// X boundaries for the max/min element
outputWorkspace->mutableX(newindex)[0] = *(X.begin() + d);
outputWorkspace->mutableX(newindex)[1] = *(X.begin() + d + 1); // This is safe since X is of dimension Y+1
outputWorkspace->mutableY(newindex)[0] = *maxY;
progress.report();
PARALLEL_END_INTERUPT_REGION
}
PARALLEL_CHECK_INTERUPT_REGION
// Assign it to the output workspace property
setProperty("OutputWorkspace", outputWorkspace);
}
} // namespace Algorithms
} // namespace Mantid