Mantid
Loading...
Searching...
No Matches
ConvertEmptyToTof.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
9#include "MantidAPI/Axis.h"
13#include "MantidAPI/Run.h"
21
22#include <cmath>
23#include <map>
24#include <numeric> // std::accumulate
25#include <utility> // std::pair
26
27namespace Mantid::Algorithms {
28
29using namespace Kernel;
30using namespace API;
31
32// Register the algorithm into the AlgorithmFactory
33DECLARE_ALGORITHM(ConvertEmptyToTof)
34
35
36const std::string ConvertEmptyToTof::name() const { return "ConvertEmptyToTof"; }
37
39int ConvertEmptyToTof::version() const { return 1; }
40
42const std::string ConvertEmptyToTof::category() const { return "Transforms\\Units"; }
43
47
48 auto wsValidator = std::make_shared<WorkspaceUnitValidator>("Empty");
50 wsValidator),
51 "Name of the input workspace");
52 declareProperty(std::make_unique<WorkspaceProperty<API::MatrixWorkspace>>("OutputWorkspace", "", Direction::Output),
53 "Name of the output workspace, can be the same as the input");
54 declareProperty(std::make_unique<Kernel::ArrayProperty<int>>("ListOfSpectraIndices"),
55 "A list of spectra workspace indices as a string with ranges; e.g. "
56 "5-10,15,20-23. \n"
57 "Optional: if not specified, then the Start/EndIndex fields "
58 "are used alone. "
59 "If specified, the range and the list are combined (without "
60 "duplicating indices). For example, a range of 10 to 20 and "
61 "a list '12,15,26,28' gives '10-20,26,28'.");
62
63 declareProperty(std::make_unique<Kernel::ArrayProperty<int>>("ListOfChannelIndices"),
64 "A list of spectra indices as a string with ranges; e.g. "
65 "5-10,15,20-23. \n"
66 "Optional: if not specified, then the Start/EndIndex fields "
67 "are used alone. "
68 "If specified, the range and the list are combined (without "
69 "duplicating indices). For example, a range of 10 to 20 and "
70 "a list '12,15,26,28' gives '10-20,26,28'.");
71
72 // OR Specify EPP
73 auto mustBePositive = std::make_shared<BoundedValidator<int>>();
74 mustBePositive->setLower(0);
75 declareProperty("ElasticPeakPosition", EMPTY_INT(), mustBePositive,
76 "Value of elastic peak position if none of the above are filled in.");
77 declareProperty("ElasticPeakPositionSpectrum", EMPTY_INT(), mustBePositive,
78 "Workspace Index used for elastic peak position above.");
79}
80
84
85 m_inputWS = this->getProperty("InputWorkspace");
86 m_outputWS = this->getProperty("OutputWorkspace");
87 m_spectrumInfo = &m_inputWS->spectrumInfo();
88 std::vector<int> spectraIndices = getProperty("ListOfSpectraIndices");
89 std::vector<int> channelIndices = getProperty("ListOfChannelIndices");
90 const int epp = getProperty("ElasticPeakPosition");
91 const int eppSpectrum = getProperty("ElasticPeakPositionSpectrum");
92
93 std::vector<double> tofAxis;
94 double channelWidth{0.0};
95 if (m_inputWS->run().hasProperty("channel_width")) {
96 channelWidth = m_inputWS->run().getPropertyValueAsType<double>("channel_width");
97 } else {
98 const std::string mesg = "No property channel_width found in the input workspace....";
99 throw std::runtime_error(mesg);
100 }
101
102 // If the ElasticPeakPosition and the ElasticPeakPositionSpectrum were
103 // specified
104 if (epp != EMPTY_INT() && eppSpectrum != EMPTY_INT()) {
105 g_log.information("Using the specified ElasticPeakPosition and "
106 "ElasticPeakPositionSpectrum");
107
108 double wavelength{0.0};
109 if (m_inputWS->run().hasProperty("wavelength")) {
110 wavelength = m_inputWS->run().getPropertyValueAsType<double>("wavelength");
111 } else {
112 std::string mesg = "No property wavelength found in the input workspace....";
113 throw std::runtime_error(mesg);
114 }
115
116 const double l1 = m_spectrumInfo->l1();
117 const double l2 = m_spectrumInfo->l2(eppSpectrum);
118 const double epTof = (calculateTOF(l1, wavelength) + calculateTOF(l2, wavelength)) * 1e6; // microsecs
119
120 tofAxis = makeTofAxis(epp, epTof, m_inputWS->blocksize() + 1, channelWidth);
121
122 }
123 // If the spectraIndices and channelIndices were specified
124 else {
125
126 // validations
127 validateWorkspaceIndices(spectraIndices);
128 validateChannelIndices(channelIndices);
129
130 // Map of spectra index, epp
131 const std::map<int, int> eppMap = findElasticPeakPositions(spectraIndices, channelIndices);
132
133 if (g_log.is(Logger::Priority::PRIO_DEBUG)) {
134 for (auto &e : eppMap) {
135 g_log.debug() << "Spectra idx =" << e.first << ", epp=" << e.second << '\n';
136 }
137 }
138
139 const std::pair<int, double> eppAndEpTof = findAverageEppAndEpTof(eppMap);
140
141 tofAxis = makeTofAxis(eppAndEpTof.first, eppAndEpTof.second, m_inputWS->blocksize() + 1, channelWidth);
142 }
143
144 // If input and output workspaces are not the same, create a new workspace for
145 // the output
146 if (m_outputWS != m_inputWS) {
147 m_outputWS = m_inputWS->clone();
148 }
149
150 setTofInWS(tofAxis, m_outputWS);
151
152 setProperty("OutputWorkspace", m_outputWS);
153}
154
161 auto nHist = m_inputWS->getNumberHistograms();
162 if (v.empty()) {
163 g_log.information("No spectrum number given. Using all spectra to calculate "
164 "the elastic peak.");
165 // use all spectra indices
166 v.reserve(nHist);
167
168 for (unsigned int i = 0; i < nHist; ++i)
169 v[i] = i;
170 } else {
171 const auto found = std::find_if(
172 v.cbegin(), v.cend(), [nHist](const auto index) { return index < 0 || static_cast<size_t>(index) >= nHist; });
173 if (found != v.cend()) {
174 throw std::runtime_error("Spectra index out of limits: " + std::to_string(*found));
175 }
176 }
177}
178
185 auto blockSize = m_inputWS->blocksize() + 1;
186 if (v.empty()) {
187 g_log.information("No channel index given. Using all channels (full "
188 "spectrum!) to calculate the elastic peak.");
189 // use all channel indices
190 v.reserve(blockSize);
191 for (unsigned int i = 0; i < blockSize; ++i)
192 v.emplace_back(i);
193 } else {
194 const auto found = std::find_if(v.cbegin(), v.cend(), [blockSize](const auto index) {
195 return index < 0 || static_cast<size_t>(index) >= blockSize;
196 });
197 if (found != v.cend()) {
198 throw std::runtime_error("Channel index out of limits: " + std::to_string(*found));
199 }
200 }
201}
202
208std::map<int, int> ConvertEmptyToTof::findElasticPeakPositions(const std::vector<int> &spectraIndices,
209 const std::vector<int> &channelIndices) {
210
211 std::map<int, int> eppMap;
212
213 // make sure we not looking for channel indices outside the bounds
214 assert(static_cast<size_t>(*(channelIndices.end() - 1)) < m_inputWS->blocksize() + 1);
215
216 g_log.information() << "Peak detection, search for peak \n";
217
218 for (auto spectrumIndex : spectraIndices) {
219
220 auto &thisSpecY = m_inputWS->y(spectrumIndex);
221
222 int minChannelIndex = *(channelIndices.begin());
223 int maxChannelIndex = *(channelIndices.end() - 1);
224
225 double center, sigma, height, minX, maxX;
226 minX = static_cast<double>(minChannelIndex);
227 maxX = static_cast<double>(maxChannelIndex);
228 estimateFWHM(thisSpecY, center, sigma, height, minX, maxX);
229
230 g_log.debug() << "Peak estimate :: center=" << center << "\t sigma=" << sigma << "\t height=" << height
231 << "\t minX=" << minX << "\t maxX=" << maxX << '\n';
232
233 bool doFit = doFitGaussianPeak(spectrumIndex, center, sigma, height, minX, maxX);
234 if (!doFit) {
235 g_log.error() << "doFitGaussianPeak failed...\n";
236 throw std::runtime_error("Gaussin Peak Fit failed....");
237 }
238
239 g_log.debug() << "Peak Fitting :: center=" << center << "\t sigma=" << sigma << "\t height=" << height
240 << "\t minX=" << minX << "\t maxX=" << maxX << '\n';
241
242 // round up the center to the closest int
243 eppMap[spectrumIndex] = roundUp(center);
244 }
245 return eppMap;
246}
247
252void ConvertEmptyToTof::estimateFWHM(const Mantid::HistogramData::HistogramY &spec, double &center, double &sigma,
253 double &height, double &minX, double &maxX) {
254
255 auto maxValueIt = std::max_element(spec.begin() + static_cast<size_t>(minX),
256 spec.begin() + static_cast<size_t>(maxX)); // max value
257 double maxValue = *maxValueIt;
258 size_t maxIndex = std::distance(spec.begin(), maxValueIt); // index of max value
259
260 auto minFwhmIt =
261 std::find_if(MantidVec::const_reverse_iterator(maxValueIt), MantidVec::const_reverse_iterator(spec.cbegin()),
262 [maxValue](double value) { return value < 0.5 * maxValue; });
263 auto maxFwhmIt = std::find_if(maxValueIt, spec.end(), [maxValue](double value) { return value < 0.5 * maxValue; });
264
265 // double fwhm = thisSpecX[maxFwhmIndex] - thisSpecX[minFwhmIndex + 1];
266 double fwhm = static_cast<double>(std::distance(minFwhmIt.base(), maxFwhmIt) + 1);
267
268 // parameters for the gaussian peak fit
269 center = static_cast<double>(maxIndex);
270 sigma = fwhm;
272
273 g_log.debug() << "Peak estimate : center=" << center << "\t sigma=" << sigma << "\t h=" << height << '\n';
274
275 // determination of the range used for the peak definition
276 size_t ipeak_min = std::max(
277 std::size_t{0}, maxIndex - static_cast<size_t>(2.5 * static_cast<double>(std::distance(maxValueIt, maxFwhmIt))));
278 size_t ipeak_max = std::min(
279 spec.size(), maxIndex + static_cast<size_t>(2.5 * static_cast<double>(std::distance(maxFwhmIt, maxValueIt))));
280 size_t i_delta_peak = ipeak_max - ipeak_min;
281
282 g_log.debug() << "Peak estimate xmin/max: " << ipeak_min - 1 << "\t" << ipeak_max + 1 << '\n';
283
284 minX = static_cast<double>(ipeak_min - 2 * i_delta_peak);
285 maxX = static_cast<double>(ipeak_max + 2 * i_delta_peak);
286}
287
301bool ConvertEmptyToTof::doFitGaussianPeak(int workspaceindex, double &center, double &sigma, double &height,
302 double startX, double endX) {
303
304 g_log.debug("Calling doFitGaussianPeak...");
305
306 // 1. Estimate
307 sigma = sigma * 0.5;
308
309 // 2. Use factory to generate Gaussian
310 auto temppeak = API::FunctionFactory::Instance().createFunction("Gaussian");
311 auto gaussianpeak = std::dynamic_pointer_cast<API::IPeakFunction>(temppeak);
312 gaussianpeak->setHeight(height);
313 gaussianpeak->setCentre(center);
314 gaussianpeak->setFwhm(sigma);
315
316 // 3. Constraint
317 double centerleftend = center - sigma * 0.5;
318 double centerrightend = center + sigma * 0.5;
319 std::ostringstream os;
320 os << centerleftend << " < PeakCentre < " << centerrightend;
321 auto centerbound = std::unique_ptr<API::IConstraint>(
322 API::ConstraintFactory::Instance().createInitialized(gaussianpeak.get(), os.str(), false));
323 gaussianpeak->addConstraint(std::move(centerbound));
324
325 g_log.debug("Calling createChildAlgorithm : Fit...");
326 // 4. Fit
327 auto fitalg = createChildAlgorithm("Fit", -1, -1, true);
328 fitalg->initialize();
329
330 fitalg->setProperty("Function", std::dynamic_pointer_cast<API::IFunction>(gaussianpeak));
331 fitalg->setProperty("InputWorkspace", m_inputWS);
332 fitalg->setProperty("WorkspaceIndex", workspaceindex);
333 fitalg->setProperty("Minimizer", "Levenberg-MarquardtMD");
334 fitalg->setProperty("CostFunction", "Least squares");
335 fitalg->setProperty("MaxIterations", 1000);
336 fitalg->setProperty("Output", "FitGaussianPeak");
337 fitalg->setProperty("StartX", startX);
338 fitalg->setProperty("EndX", endX);
339
340 // 5. Result
341 bool successfulfit = fitalg->execute();
342 if (!fitalg->isExecuted() || !successfulfit) {
343 // Early return due to bad fit
344 g_log.warning() << "Fitting Gaussian peak for peak around " << gaussianpeak->centre() << '\n';
345 return false;
346 }
347
348 // 6. Get result
349 center = gaussianpeak->centre();
350 height = gaussianpeak->height();
351 double fwhm = gaussianpeak->fwhm();
352
353 return fwhm > 0.0;
354}
355
362std::pair<int, double> ConvertEmptyToTof::findAverageEppAndEpTof(const std::map<int, int> &eppMap) {
363
364 double l1 = m_spectrumInfo->l1();
365 double wavelength{0.0};
366 if (m_inputWS->run().hasProperty("wavelength")) {
367 wavelength = m_inputWS->run().getPropertyValueAsType<double>("wavelength");
368 } else {
369 std::string mesg = "No property wavelength found in the input workspace....";
370 throw std::runtime_error(mesg);
371 }
372
373 std::vector<double> epTofList;
374 std::vector<int> eppList;
375
376 double firstL2 = m_spectrumInfo->l2(eppMap.begin()->first);
377 for (const auto &epp : eppMap) {
378
379 double l2 = m_spectrumInfo->l2(epp.first);
380 if (!areEqual(l2, firstL2, 0.0001)) {
381 g_log.error() << "firstL2=" << firstL2 << " , "
382 << "l2=" << l2 << '\n';
383 throw std::runtime_error("All the pixels for selected spectra must have "
384 "the same distance from the sample!");
385 } else {
386 firstL2 = l2;
387 }
388
389 epTofList.emplace_back((calculateTOF(l1, wavelength) + calculateTOF(l2, wavelength)) * 1e6); // microsecs
390 eppList.emplace_back(epp.first);
391
392 g_log.debug() << "WS index = " << epp.first << ", l1 = " << l1 << ", l2 = " << l2
393 << ", TOF(l1+l2) = " << *(epTofList.end() - 1) << '\n';
394 }
395
396 double averageEpTof =
397 std::accumulate(epTofList.begin(), epTofList.end(), 0.0) / static_cast<double>(epTofList.size());
398 int averageEpp = roundUp(static_cast<double>(std::accumulate(eppList.begin(), eppList.end(), 0)) /
399 static_cast<double>(eppList.size()));
400
401 g_log.debug() << "Average epp=" << averageEpp << " , Average epTof=" << averageEpTof << '\n';
402 return std::make_pair(averageEpp, averageEpTof);
403}
404
405double ConvertEmptyToTof::calculateTOF(double distance, double wavelength) {
406 if (wavelength <= 0) {
407 throw std::runtime_error("Wavelenght is <= 0");
408 }
409
410 double velocity = PhysicalConstants::h / (PhysicalConstants::NeutronMass * wavelength * 1e-10); // m/s
411
412 return distance / velocity;
413}
414
418bool ConvertEmptyToTof::areEqual(double a, double b, double epsilon) { return fabs(a - b) < epsilon; }
419
420int ConvertEmptyToTof::roundUp(double value) { return static_cast<int>(std::floor(value + 0.5)); }
421
425std::vector<double> ConvertEmptyToTof::makeTofAxis(int epp, double epTof, size_t size, double channelWidth) {
426 std::vector<double> axis(size);
427
428 g_log.debug() << "Building the TOF X Axis: epp=" << epp << ", epTof=" << epTof << ", Channel Width=" << channelWidth
429 << '\n';
430 for (size_t i = 0; i < size; ++i) {
431 axis[i] = epTof + channelWidth * static_cast<double>(static_cast<int>(i) - epp) -
432 channelWidth / 2; // to make sure the bin is in the middle of the elastic peak
433 }
434 g_log.debug() << "TOF X Axis: [start,end] = [" << *axis.begin() << "," << *(axis.end() - 1) << "]\n";
435 return axis;
436}
437
438void ConvertEmptyToTof::setTofInWS(const std::vector<double> &tofAxis, const API::MatrixWorkspace_sptr &outputWS) {
439
440 const size_t numberOfSpectra = m_inputWS->getNumberHistograms();
441
442 g_log.debug() << "Setting the TOF X Axis for numberOfSpectra=" << numberOfSpectra << '\n';
443
444 HistogramData::BinEdges edges(tofAxis);
445 Progress prog(this, 0.0, 0.2, numberOfSpectra);
446
447 for (size_t i = 0; i < numberOfSpectra; ++i) {
448 // Replace bin edges with tof axis
449 outputWS->setBinEdges(i, edges);
450
451 prog.report();
452 } // end for i
453
454 outputWS->getAxis(0)->unit() = UnitFactory::Instance().create("TOF");
455}
456
457} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
double maxValue
double value
The value of the point.
Definition: FitMW.cpp:51
double sigma
Definition: GetAllEi.cpp:156
double height
Definition: GetAllEi.cpp:155
std::map< DeltaEMode::Type, std::string > index
Definition: DeltaEMode.cpp:19
#define fabs(x)
Definition: Matrix.cpp:22
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
virtual std::shared_ptr< Algorithm > createChildAlgorithm(const std::string &name, const double startProgress=-1., const double endProgress=-1., const bool enableLogging=true, const int &version=-1)
Create a Child Algorithm.
Definition: Algorithm.cpp:842
Kernel::Logger & g_log
Definition: Algorithm.h:451
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
double l2(const size_t index) const
Returns L2 (distance from sample to spectrum).
double l1() const
Returns L1 (distance from source to sample).
A property class for workspaces.
bool doFitGaussianPeak(int, double &, double &, double &, double, double)
Fit peak without background i.e, with background removed inspired from FitPowderDiffPeaks....
API::MatrixWorkspace_sptr m_outputWS
void validateChannelIndices(std::vector< int > &v)
Check if the channel indices are in the limits of the number of the block size in the input workspace...
void validateWorkspaceIndices(std::vector< int > &v)
Check if spectra indices are in the limits of the number of histograms in the input workspace.
void exec() override
Execute the algorithm.
const std::string category() const override
Algorithm's category for identification.
void estimateFWHM(const Mantid::HistogramData::HistogramY &, double &, double &, double &, double &, double &)
Estimated the FWHM for Gaussian peak fitting.
std::map< int, int > findElasticPeakPositions(const std::vector< int > &, const std::vector< int > &)
Looks for the EPP positions in the spectraIndices.
bool areEqual(double, double, double)
Compare two double with a precision epsilon.
int version() const override
Algorithm's version for identification.
std::vector< double > makeTofAxis(int, double, size_t, double)
Builds the X time axis.
void init() override
Initialize the algorithm's properties.
std::pair< int, double > findAverageEppAndEpTof(const std::map< int, int > &)
Finds the TOF for a given epp.
DataObjects::Workspace2D_sptr m_inputWS
void setTofInWS(const std::vector< double > &, const API::MatrixWorkspace_sptr &)
const API::SpectrumInfo * m_spectrumInfo
Support for a property that holds an array of values.
Definition: ArrayProperty.h:28
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void debug(const std::string &msg)
Logs at debug level.
Definition: Logger.cpp:114
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
bool is(int level) const
Returns true if at least the given log level is set.
Definition: Logger.cpp:146
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
Definition: ProgressBase.h:51
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
static constexpr double NeutronMass
Mass of the neutron in kg.
static constexpr double h
Planck constant in J*s.
constexpr int EMPTY_INT() noexcept
Returns what we consider an "empty" integer within a property.
Definition: EmptyValues.h:25
STL namespace.
std::string to_string(const wide_integer< Bits, Signed > &n)
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54