Mantid
Loading...
Searching...
No Matches
IntegrateFlux.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 +
12#include "MantidHistogramData/LinearGenerator.h"
14
15#include <memory>
16#include <numeric>
17
18namespace Mantid::MDAlgorithms {
19
22using namespace Mantid::HistogramData;
23using Mantid::Types::Event::TofEvent;
24
25// Register the algorithm into the AlgorithmFactory
26DECLARE_ALGORITHM(IntegrateFlux)
27
28namespace {
29
31class NoEventWorkspaceDeleting {
32public:
34 void operator()(const API::MatrixWorkspace * /*unused*/) {}
35};
36} // namespace
37
38//----------------------------------------------------------------------------------------------
39
41const std::string IntegrateFlux::name() const { return "IntegrateFlux"; }
42
44int IntegrateFlux::version() const { return 1; }
45
47const std::string IntegrateFlux::category() const { return "MDAlgorithms\\Normalisation"; }
48
50const std::string IntegrateFlux::summary() const {
51 return "Integrates spectra in a matrix workspace at a set of points.";
52}
53
54//----------------------------------------------------------------------------------------------
60 "InputWorkspace", "", Direction::Input, std::make_shared<API::WorkspaceUnitValidator>("Momentum")),
61 "An input workspace. Must have units of Momentum");
62 auto validator = std::make_shared<Kernel::BoundedValidator<int>>();
63 validator->setLower(2);
64 declareProperty("NPoints", 1000, validator, "Number of points per output spectrum.");
65 declareProperty(std::make_unique<WorkspaceProperty<API::Workspace>>("OutputWorkspace", "", Direction::Output),
66 "An output workspace.");
67}
68
69//----------------------------------------------------------------------------------------------
73 API::MatrixWorkspace_sptr inputWS = getProperty("InputWorkspace");
74 size_t nX = static_cast<size_t>(static_cast<int>(getProperty("NPoints")));
75
76 auto outputWS = createOutputWorkspace(*inputWS, nX);
77
78 integrateSpectra(*inputWS, *outputWS);
79
80 setProperty("OutputWorkspace", outputWS);
81}
82
90std::shared_ptr<API::MatrixWorkspace> IntegrateFlux::createOutputWorkspace(const API::MatrixWorkspace &inputWS,
91 size_t nX) const {
92 size_t nSpec = inputWS.getNumberHistograms();
93
94 if (nSpec == 0) {
95 throw std::runtime_error("Input workspace has no data.");
96 }
97
98 // make sure the output spectrum size isn't too large
99 auto maxPoints = getMaxNumberOfPoints(inputWS);
100 if (nX > maxPoints) {
101 nX = maxPoints;
102 }
103
104 // and not 0 or 1 as they are to be used for interpolation
105 if (nX < 2) {
106 throw std::runtime_error("Failed to create output."
107 "Output spectra should have at least two points.");
108 }
109
110 // crate empty output workspace
112 std::shared_ptr<const API::MatrixWorkspace>(&inputWS, NoEventWorkspaceDeleting()), nSpec, nX, nX);
113
114 // claculate the integration points and save them in the x-vactors of
115 // integrFlux
116 double xMin = inputWS.getXMin();
117 double xMax = inputWS.getXMax();
118 double dx = (xMax - xMin) / static_cast<double>(nX - 1);
119 const auto &x = ws->x(0);
120
121 ws->setPoints(0, Points(x.size(), LinearGenerator(xMin, dx)));
122 for (size_t sp = 1; sp < nSpec; ++sp)
123 ws->setSharedX(sp, ws->sharedX(0));
124
125 return ws;
126}
127
136 auto eventWS = dynamic_cast<const DataObjects::EventWorkspace *>(&inputWS);
137
138 if (eventWS) {
139 auto eventType = eventWS->getEventType();
140 switch (eventType) {
142 integrateSpectraEvents<DataObjects::WeightedEventNoTime>(*eventWS, integrWS);
143 return;
144 case (API::WEIGHTED):
145 integrateSpectraEvents<DataObjects::WeightedEvent>(*eventWS, integrWS);
146 return;
147 case (API::TOF):
148 integrateSpectraEvents<TofEvent>(*eventWS, integrWS);
149 return;
150 }
151 } else {
152 integrateSpectraMatrix(inputWS, integrWS);
153 }
154}
155
162template <class EventType>
164 API::MatrixWorkspace &integrWS) const {
165 inputWS.sortAll(DataObjects::TOF_SORT, nullptr);
166 size_t nSpec = inputWS.getNumberHistograms();
167 assert(nSpec == integrWS.getNumberHistograms());
168
169 auto &X = integrWS.x(0);
170 // loop overr the spectra and integrate
171 for (size_t sp = 0; sp < nSpec; ++sp) {
172 const std::vector<EventType> *el;
174 auto &outY = integrWS.mutableY(sp);
175
176 double sum = 0;
177 auto x = X.begin() + 1;
178 size_t i = 1;
179 // the integral is a running sum of the event weights in the spectrum
180 for (auto evnt = el->begin(); evnt != el->end(); ++evnt) {
181 double tof = evnt->tof();
182 while (x != X.end() && *x < tof) {
183 outY[i] = sum;
184 ++x;
185 ++i;
186 }
187 if (x == X.end())
188 break;
189 sum += evnt->weight();
190 outY[i] = sum;
191 }
192
193 while (x != X.end()) {
194 outY[i] = sum;
195 ++x;
196 ++i;
197 }
198 }
199}
200
208 bool isHistogram = inputWS.isHistogramData();
209
210 if (isHistogram) {
211 integrateSpectraHistograms(inputWS, integrWS);
212 } else {
213 integrateSpectraPointData(inputWS, integrWS);
214 }
215}
216
224 API::MatrixWorkspace &integrWS) const {
225 size_t nSpec = inputWS.getNumberHistograms();
226 assert(nSpec == integrWS.getNumberHistograms());
227
228 auto &X = integrWS.x(0);
229
230 // loop over the spectra and integrate
231 for (size_t sp = 0; sp < nSpec; ++sp) {
232 auto &inX = inputWS.x(sp);
233 auto inY = inputWS.counts(sp); // make a copy
234
235 // integral at the first point is always 0
236 auto outY = integrWS.mutableY(sp).begin();
237 *outY = 0.0;
238 ++outY;
239 // initialize summation
240 double sum = 0;
241 // cache some iterators
242 auto inXbegin = inX.begin();
243 auto inXend = inX.end();
244 auto x0 = inXbegin; // iterator over x in input workspace
245 // loop over the iteration points starting from the second one
246 for (auto outX = X.begin() + 1; outX != X.end(); ++outX, ++outY) {
247 // there are no data to integrate
248 if (x0 == inXend) {
249 *outY = sum;
250 continue;
251 }
252
253 // in each iteration we find the integral of the input spectrum
254 // between bounds [lowerBound,upperBound]
255 const double &lowerBound = *(outX - 1);
256 double upperBound = *outX;
257
258 // interval [*x0, *x1] is the smalest interval in inX that contains
259 // the integration interval [lowerBound,upperBound]
260 auto x1 = std::lower_bound(x0, inXend, upperBound);
261
262 // reached end of input data
263 if (x1 == inXend) {
264 --x1;
265 if (x1 == x0) {
266 *outY = sum;
267 x0 = inXend;
268 continue;
269 }
270 upperBound = *x1;
271 }
272
273 // if starting point in input x is smaller (not equal) than the lower
274 // integration bound
275 // then there is a partial bin at the beginning of the interval
276 if (*x0 < lowerBound) {
277 // first find the part of bin [*x0,*(x0+1)] which hasn't been integrated
278 // yet
279 // the left boundary == lowerBound
280 // the right boundary == min( upperBound, *(x0+1) )
281 const double leftX = lowerBound;
282 const double rightX = std::min(upperBound, *(x0 + 1));
283
284 auto i = static_cast<size_t>(std::distance(inXbegin, x0));
285 // add bin's fraction between leftX and rightX
286 sum += inY[i] * (rightX - leftX) / (*(x0 + 1) - *x0);
287
288 // if rightX == upperBound there is nothing left to integrate, move to
289 // the next integration point
290 if (rightX == upperBound) {
291 *outY = sum;
292 continue;
293 }
294
295 ++x0;
296 }
297
298 // accumulate values in bins that fit entirely into the integration
299 // interval [lowerBound,upperBound]
300 auto i0 = static_cast<size_t>(std::distance(inXbegin, x0));
301 auto i1 = static_cast<size_t>(std::distance(inXbegin, x1));
302 if (*x1 > upperBound)
303 --i1;
304 for (auto i = i0; i < i1; ++i) {
305 sum += inY[i];
306 }
307
308 // if x1 is greater than upperBound there is a partial "bin" that has to
309 // be added
310 if (*x1 > upperBound) {
311 // find the part of "bin" [*(x1-1),*x1] which needs to be integrated
312 // the left boundary == *(x1-1)
313 // the right boundary == upperBound
314 const double leftX = *(x1 - 1);
315 const double rightX = upperBound;
316
317 auto i = static_cast<size_t>(std::distance(inXbegin, x1));
318 // add the area under the line between leftX and rightX
319 sum += inY[i - 1] * (rightX - leftX) / (*x1 - *(x1 - 1));
320
321 // advance in the input workspace
322 x0 = x1 - 1;
323 } else {
324 // advance in the input workspace
325 x0 = x1;
326 }
327
328 // store the current sum
329 *outY = sum;
330 }
331 }
332}
333
341 API::MatrixWorkspace &integrWS) const {
342 size_t nSpec = inputWS.getNumberHistograms();
343 assert(nSpec == integrWS.getNumberHistograms());
344
345 auto &X = integrWS.x(0);
346
347 // loop overr the spectra and integrate
348 for (size_t sp = 0; sp < nSpec; ++sp) {
349 auto &inX = inputWS.x(sp);
350 auto &inY = inputWS.y(sp);
351
352 // integral at the first point is always 0
353 auto outY = integrWS.mutableY(sp).begin();
354 *outY = 0.0;
355 ++outY;
356 // initialize summation
357 double sum = 0;
358 // cache some iterators
359 auto inXbegin = inX.begin();
360 auto inXend = inX.end();
361 auto x0 = inXbegin; // iterator over x in input workspace
362
363 // loop over the iteration points starting from the second one
364 for (auto outX = X.begin() + 1; outX != X.end(); ++outX, ++outY) {
365 // there are no data to integrate
366 if (x0 == inXend) {
367 *outY = sum;
368 continue;
369 }
370
371 // in each iteration we find the integral of the input spectrum
372 // between bounds [lowerBound,upperBound]
373 const double &lowerBound = *(outX - 1);
374 double upperBound = *outX;
375
376 // interval [*x0, *x1] is the smalest interval in inX that contains
377 // the integration interval [lowerBound,upperBound]
378 auto x1 = std::lower_bound(x0, inXend, upperBound);
379
380 // reached end of input data
381 if (x1 == inXend) {
382 --x1;
383 if (x1 == x0) {
384 *outY = sum;
385 x0 = inXend;
386 continue;
387 }
388 upperBound = *x1;
389 }
390
391 // if starting point in input x is smaller (not equal) than the lower
392 // integration bound
393 // then there is a partial "bin" at the beginning of the interval
394 if (*x0 < lowerBound) {
395 // first find the part of "bin" [*x0,*(x0+1)] which hasn't been
396 // integrated yet
397 // the left boundary == lowerBound
398 // the right boundary == min( upperBound, *(x0+1) )
399 const double leftX = lowerBound;
400 const double rightX = std::min(upperBound, *(x0 + 1));
401
402 auto i = static_cast<size_t>(std::distance(inXbegin, x0));
403 // gradient of "bin" [*x0,*(x0+1)]
404 double dy_dx = (inY[i + 1] - inY[i]) / (*(x0 + 1) - *x0);
405
406 // add the area under the line between leftX and rightX
407 sum += (inY[i] + 0.5 * dy_dx * (leftX + rightX - 2 * (*(x0)))) * (rightX - leftX);
408
409 // if rightX == upperBound there is nothing left to integrate, move to
410 // the next integration point
411 if (rightX == upperBound) {
412 *outY = sum;
413 continue;
414 }
415
416 ++x0;
417 }
418
419 // accumulate values in bins that fit entirely into the integration
420 // interval [lowerBound,upperBound]
421 auto i0 = static_cast<size_t>(std::distance(inXbegin, x0));
422 auto i1 = static_cast<size_t>(std::distance(inXbegin, x1));
423 if (*x1 > upperBound)
424 --i1;
425
426 for (auto i = i0; i < i1; ++i) {
427 sum += (inY[i] + inY[i + 1]) / 2 * (inX[i + 1] - inX[i]);
428 }
429
430 // if x1 is greater than upperBound there is a partial "bin" that has to
431 // be added
432 if (*x1 > upperBound) {
433 // find the part of "bin" [*(x1-1),*x1] which needs to be integrated
434 // the left boundary == *(x1-1)
435 // the right boundary == upperBound
436 const double leftX = *(x1 - 1);
437 const double rightX = upperBound;
438
439 auto i = static_cast<size_t>(std::distance(inXbegin, x1));
440 // gradient of "bin" [*(x1-1),*x1]
441 double dy_dx = (inY[i] - inY[i - 1]) / (*x1 - *(x1 - 1));
442
443 // add the area under the line between leftX and rightX
444 sum += (inY[i - 1] + 0.5 * dy_dx * (rightX - *(x1 - 1))) * (rightX - leftX);
445
446 // advance in the input workspace
447 x0 = x1 - 1;
448 } else {
449 // advance in the input workspace
450 x0 = x1;
451 }
452
453 // store the current sum
454 *outY = sum;
455 }
456 }
457}
458
464 // if it's events we shouldn't care about binning
465 auto eventWS = dynamic_cast<const DataObjects::EventWorkspace *>(&inputWS);
466 if (eventWS) {
467 return eventWS->getSpectrum(0).getNumberEvents();
468 }
469
470 return inputWS.blocksize();
471}
472
473} // namespace Mantid::MDAlgorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
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
Base MatrixWorkspace Abstract Class.
virtual std::size_t blocksize() const =0
Returns the size of each block of data returned by the dataY accessors.
virtual double getXMin() const
virtual std::size_t getNumberHistograms() const =0
Returns the number of histograms in the workspace.
const HistogramData::HistogramX & x(const size_t index) const
HistogramData::Counts counts(const size_t index) const
virtual double getXMax() const
HistogramData::HistogramY & mutableY(const size_t index) &
const HistogramData::HistogramY & y(const size_t index) const
virtual bool isHistogramData() const
Returns true if the workspace contains data in histogram form (as opposed to point-like)
A property class for workspaces.
std::size_t getNumberEvents() const override
Return the number of events in the list.
Definition: EventList.cpp:1143
This class is intended to fulfill the design specified in <https://github.com/mantidproject/documents...
EventList & getSpectrum(const size_t index) override
Return the underlying ISpectrum ptr at the given workspace index.
std::size_t getNumberHistograms() const override
Get the number of histograms, usually the same as the number of pixels or detectors.
void sortAll(EventSortType sortType, Mantid::API::Progress *prog) const
Mantid::API::EventType getEventType() const override
Get the EventType of the most-specialized EventList in the workspace.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
const std::string name() const override
Algorithms name for identification.
void integrateSpectraHistograms(const API::MatrixWorkspace &inputWS, API::MatrixWorkspace &integrWS) const
Integrate spectra in inputWS at x-values in integrWS and save the results in y-vectors of integrWS.
int version() const override
Algorithm's version for identification.
void integrateSpectraEvents(const DataObjects::EventWorkspace &inputWS, API::MatrixWorkspace &integrWS) const
Integrate spectra in inputWS at x-values in integrWS and save the results in y-vectors of integrWS.
void init() override
Initialize the algorithm's properties.
void integrateSpectra(const API::MatrixWorkspace &inputWS, API::MatrixWorkspace &integrWS) const
Integrate spectra in inputWS at x-values in integrWS and save the results in y-vectors of integrWS.
std::shared_ptr< API::MatrixWorkspace > createOutputWorkspace(const API::MatrixWorkspace &inputWS, size_t nX) const
Create an empty output workspace with required dimensions and defined x-values.
size_t getMaxNumberOfPoints(const API::MatrixWorkspace &inputWS) const
Calculate the maximun number of points in the integration grid.
const std::string category() const override
Algorithm's category for identification.
void integrateSpectraMatrix(const API::MatrixWorkspace &inputWS, API::MatrixWorkspace &integrWS) const
Integrate spectra in inputWS at x-values in integrWS and save the results in y-vectors of integrWS.
void integrateSpectraPointData(const API::MatrixWorkspace &inputWS, API::MatrixWorkspace &integrWS) const
Integrate spectra in inputWS at x-values in integrWS and save the results in y-vectors of integrWS.
void exec() override
Execute the algorithm.
const std::string summary() const override
Algorithm's summary for use in the GUI and help.
@ WEIGHTED_NOTIME
Definition: IEventList.h:18
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
DLLExport void getEventsFrom(EventList &el, std::vector< Types::Event::TofEvent > *&events)
Describes the direction (within an algorithm) of a Property.
Definition: Property.h:50
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54