Mantid
Loading...
Searching...
No Matches
PeakIntegration.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 "MantidAPI/Sample.h"
21
22#include <boost/math/special_functions/round.hpp>
23#include <cmath>
24
25namespace Mantid::Crystal {
26
27// Register the class into the algorithm factory
28DECLARE_ALGORITHM(PeakIntegration)
29
30using namespace Kernel;
31using namespace Geometry;
32using namespace API;
33using namespace DataObjects;
34
39
40 declareProperty(std::make_unique<WorkspaceProperty<PeaksWorkspace>>("InPeaksWorkspace", "", Direction::Input),
41 "Name of the peaks workspace.");
42 declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input,
43 std::make_shared<InstrumentValidator>()),
44 "A 2D workspace with X values of time of flight");
45 declareProperty(std::make_unique<WorkspaceProperty<PeaksWorkspace>>("OutPeaksWorkspace", "", Direction::Output),
46 "Name of the output peaks workspace with integrated intensities.");
47 declareProperty("IkedaCarpenterTOF", false,
48 "Integrate TOF using IkedaCarpenter fit.\n"
49 "Default is false which is best for corrected data.");
50 declareProperty("MatchingRunNo", true,
51 "Integrate only peaks where run "
52 "number of peak matches run number of "
53 "sample.\n"
54 "Default is true.");
55 declareProperty("NBadEdgePixels", 0, "Number of bad Edge Pixels");
56}
57
65
67 PeaksWorkspace_sptr inPeaksW = getProperty("InPeaksWorkspace");
68
70 PeaksWorkspace_sptr peaksW = getProperty("OutPeaksWorkspace");
71 if (peaksW != inPeaksW)
72 peaksW = inPeaksW->clone();
73
74 double qspan = 0.12;
75 m_IC = getProperty("IkedaCarpenterTOF");
76 bool matchRun = getProperty("MatchingRunNo");
77 if (peaksW->mutableSample().hasOrientedLattice()) {
78 OrientedLattice latt = peaksW->mutableSample().getOrientedLattice();
79 qspan = 1. / std::max(latt.a(), std::max(latt.b(), latt.c())); // 1/6*2Pi about 1
80
81 } else {
82 qspan = 0.12;
83 }
84
85 // To get the workspace index from the detector ID
86 const auto pixel_to_wi = inputW->getDetectorIDToWorkspaceIndexMap();
87
88 // Sort events if EventWorkspace so it will run in parallel
89 EventWorkspace_const_sptr inWS = std::dynamic_pointer_cast<const EventWorkspace>(inputW);
90 if (inWS) {
91 inWS->sortAll(TOF_SORT, nullptr);
92 }
93
94 // Get some stuff from the input workspace
95 const auto YLength = static_cast<int>(inputW->blocksize());
96 outputW = API::WorkspaceFactory::Instance().create(inputW, peaksW->getNumberPeaks(), YLength + 1, YLength);
97 // Copy geometry over.
98 API::WorkspaceFactory::Instance().initializeFromParent(*inputW, *outputW, true);
99 size_t Numberwi = inputW->getNumberHistograms();
100 int NumberPeaks = peaksW->getNumberPeaks();
101 int MinPeaks = 0;
102
103 std::vector<int> badPeaks;
104 for (int i = NumberPeaks - 1; i >= 0; i--) {
105 Peak &peak = peaksW->getPeaks()[i];
106 int pixelID = peak.getDetectorID();
107
108 // Find the workspace index for this detector ID
109 auto wiEntry = pixel_to_wi.find(pixelID);
110 if (wiEntry != pixel_to_wi.end()) {
111 size_t wi = wiEntry->second;
112 if ((matchRun && peak.getRunNumber() != inputW->getRunNumber()) || wi >= Numberwi)
113 badPeaks.emplace_back(i);
114 } else // This is for appending peak workspaces when running
115 // SNSSingleCrystalReduction one bank at at time
116 if (i + 1 > MinPeaks)
117 MinPeaks = i + 1;
118 }
119 peaksW->removePeaks(std::move(badPeaks));
120 NumberPeaks = peaksW->getNumberPeaks();
121 if (NumberPeaks <= 0) {
122 g_log.error("RunNumbers of InPeaksWorkspace and InputWorkspace do not match");
123 return;
124 }
125
126 Progress prog(this, MinPeaks, 1.0, NumberPeaks);
128 for (int i = MinPeaks; i < NumberPeaks; i++) {
129
131
132 // Direct ref to that peak
133 Peak &peak = peaksW->getPeaks()[i];
134
135 double col = peak.getCol();
136 double row = peak.getRow();
137
138 // Average integer postion
139 int XPeak = boost::math::iround(col);
140 int YPeak = boost::math::iround(row);
141
142 double TOFPeakd = peak.getTOF();
143 std::string bankName = peak.getBankName();
144
145 std::shared_ptr<const IComponent> parent = inputW->getInstrument()->getComponentByName(bankName);
146
147 if (!parent)
148 continue;
149
150 int TOFPeak = 0, TOFmin = 0, TOFmax = 0;
151 TOFmax = fitneighbours(i, bankName, XPeak, YPeak, i, qspan, peaksW, pixel_to_wi);
152
153 double I = 0., sigI = 0.;
154 // Find point of peak centre
155 // Get references to the current spectrum
156 const auto &X = outputW->x(i);
157 const auto &Y = outputW->y(i);
158 const auto &E = outputW->e(i);
159 auto numbins = static_cast<int>(Y.size());
160 if (TOFmin > numbins)
161 TOFmin = numbins;
162 if (TOFmax > numbins)
163 TOFmax = numbins;
164
165 TOFPeak = VectorHelper::getBinIndex(X.rawData(), TOFPeakd);
166 const double peakLoc = X[TOFPeak];
167 int iTOF;
168 for (iTOF = TOFmin; iTOF < TOFmax; iTOF++) {
169 if (Y[iTOF] > 0.0 && Y[iTOF + 1] > 0.0)
170 break;
171 }
172 TOFmin = iTOF;
173 for (iTOF = TOFmax; iTOF > TOFmin; iTOF--) {
174 if (Y[iTOF] > 0.0 && Y[iTOF - 1] > 0.0)
175 break;
176 }
177 TOFmax = iTOF;
178 if (TOFmax <= TOFmin)
179 continue;
180 const int n = TOFmax - TOFmin + 1;
181 // double pktime = 0.0;
182
183 // for (iTOF = TOFmin; iTOF < TOFmax; iTOF++) pktime+= X[iTOF];
184 if (n >= 8 && m_IC) // Number of fitting parameters large enough if
185 // Ikeda-Carpenter fit
186 {
187 for (iTOF = TOFmin; iTOF <= TOFmax; iTOF++) {
188 if (((Y[iTOF] - Y[TOFPeak] / 2.) * (Y[iTOF + 1] - Y[TOFPeak] / 2.)) < 0.)
189 break;
190 }
191 double Gamma = fabs(X[TOFPeak] - X[iTOF]);
192 double SigmaSquared = Gamma * Gamma;
193 const double peakHeight = Y[TOFPeak] * Gamma; // Intensity*HWHM
194
195 IAlgorithm_sptr fit_alg;
196 try {
197 fit_alg = createChildAlgorithm("Fit", -1, -1, false);
198 } catch (Exception::NotFoundError &) {
199 g_log.error("Can't locate Fit algorithm");
200 throw;
201 }
202 // Initialize Ikeda-Carpender function variables
203 double Alpha0 = 1.6;
204 double Alpha1 = 1.5;
205 double Beta0 = 31.9;
206 double Kappa = 46.0;
207 std::ostringstream fun_str;
208 fun_str << "name=IkedaCarpenterPV,I=" << peakHeight << ",Alpha0=" << Alpha0 << ",Alpha1=" << Alpha1
209 << ",Beta0=" << Beta0 << ",Kappa=" << Kappa << ",SigmaSquared=" << SigmaSquared << ",Gamma=" << Gamma
210 << ",X0=" << peakLoc;
211 fit_alg->setPropertyValue("Function", fun_str.str());
212 if (Alpha0 != 1.6 || Alpha1 != 1.5 || Beta0 != 31.9 || Kappa != 46.0) {
213 std::ostringstream tie_str;
214 tie_str << "Alpha0=" << Alpha0 << ",Alpha1=" << Alpha1 << ",Beta0=" << Beta0 << ",Kappa=" << Kappa;
215 fit_alg->setProperty("Ties", tie_str.str());
216 }
217 fit_alg->setProperty("InputWorkspace", outputW);
218 fit_alg->setProperty("WorkspaceIndex", i);
219 fit_alg->setProperty("StartX", X[TOFmin]);
220 fit_alg->setProperty("EndX", X[TOFmax]);
221 fit_alg->setProperty("MaxIterations", 5000);
222 fit_alg->setProperty("CreateOutput", true);
223 fit_alg->setProperty("Output", "fit");
224 fit_alg->executeAsChildAlg();
225 MatrixWorkspace_sptr fitWS = fit_alg->getProperty("OutputWorkspace");
226
227 /*double chisq = fit_alg->getProperty("OutputChi2overDoF");
228 if(chisq > 0 && chisq < 400 && !haveFit && PeakIntensity < 30.0) // use
229 fit of strong peaks for weak peaks
230 {
231 std::vector<double> params = fit_alg->getProperty("Parameters");
232 Alpha0 = params[1];
233 Alpha1 = params[2];
234 Beta0 = params[3];
235 Kappa = params[4];
236 haveFit = true;
237 }
238 std::string funct = fit_alg->getPropertyValue("Function");
239 */
240
241 // Evaluate fit at points
242 const auto &y = fitWS->y(1);
243
244 // Calculate intensity
245 for (iTOF = 0; iTOF < n; iTOF++)
246 if (std::isfinite(y[iTOF]))
247 I += y[iTOF];
248 } else
249 for (iTOF = TOFmin; iTOF <= TOFmax; iTOF++)
250 I += Y[iTOF];
251
252 if (!m_IC) {
253 sigI = peak.getSigmaIntensity();
254 } else {
255 // Calculate errors correctly for nonPoisson distributions
256 for (iTOF = TOFmin; iTOF <= TOFmax; iTOF++)
257 sigI += E[iTOF] * E[iTOF];
258 sigI = sqrt(sigI);
259 }
260
261 peak.setIntensity(I);
262 peak.setSigmaIntensity(sigI);
263
264 prog.report();
266 }
267
269
270 // Save the output
271 setProperty("OutPeaksWorkspace", peaksW);
272}
273
275 inputW = getProperty("InputWorkspace");
276 if (inputW->y(0).size() <= 1)
277 throw std::runtime_error("Must Rebin data with more than 1 bin");
278 // Check if detectors are RectangularDetectors
279 Instrument_const_sptr inst = inputW->getInstrument();
280 std::shared_ptr<RectangularDetector> det;
281 for (int i = 0; i < inst->nelements(); i++) {
282 det = std::dynamic_pointer_cast<RectangularDetector>((*inst)[i]);
283 if (det)
284 break;
285 }
286}
287
288int PeakIntegration::fitneighbours(int ipeak, const std::string &det_name, int x0, int y0, int idet, double qspan,
289 PeaksWorkspace_sptr &Peaks, const detid2index_map &pixel_to_wi) {
290 UNUSED_ARG(det_name);
291 UNUSED_ARG(x0);
292 UNUSED_ARG(y0);
293 Peak &peak = Peaks->getPeak(ipeak);
294 // Number of slices
295 int TOFmax = 0;
296
297 auto slice_alg = createChildAlgorithm("IntegratePeakTimeSlices");
298 slice_alg->setProperty<MatrixWorkspace_sptr>("InputWorkspace", inputW);
299 std::ostringstream tab_str;
300 tab_str << "LogTable" << ipeak;
301
302 slice_alg->setPropertyValue("OutputWorkspace", tab_str.str());
303 slice_alg->setProperty<PeaksWorkspace_sptr>("Peaks", Peaks);
304 slice_alg->setProperty("PeakIndex", ipeak);
305 slice_alg->setProperty("PeakQspan", qspan);
306
307 int nPixels = std::max<int>(0, getProperty("NBadEdgePixels"));
308
309 slice_alg->setProperty("NBadEdgePixels", nPixels);
310 slice_alg->executeAsChildAlg();
311
312 auto &Xout = outputW->mutableX(idet);
313 auto &Yout = outputW->mutableY(idet);
314 auto &Eout = outputW->mutableE(idet);
315 TableWorkspace_sptr logtable = slice_alg->getProperty("OutputWorkspace");
316
317 peak.setIntensity(slice_alg->getProperty("Intensity"));
318 peak.setSigmaIntensity(slice_alg->getProperty("SigmaIntensity"));
319
320 TOFmax = static_cast<int>(logtable->rowCount());
321 for (int iTOF = 0; iTOF < TOFmax; iTOF++) {
322 Xout[iTOF] = logtable->getRef<double>(std::string("Time"), iTOF);
323 if (m_IC) // Ikeda-Carpenter fit
324 {
325 Yout[iTOF] = logtable->getRef<double>(std::string("TotIntensity"), iTOF);
326 Eout[iTOF] = logtable->getRef<double>(std::string("TotIntensityError"), iTOF);
327 } else {
328 Yout[iTOF] = logtable->getRef<double>(std::string("ISAWIntensity"), iTOF);
329 Eout[iTOF] = logtable->getRef<double>(std::string("ISAWIntensityError"), iTOF);
330 }
331 }
332
333 outputW->getSpectrum(idet).clearDetectorIDs();
334 // Find the pixel ID at that XY position on the rectangular detector
335 int pixelID = peak.getDetectorID(); // det->getAtXY(x0,y0)->getID();
336
337 // Find the corresponding workspace index, if any
338 auto wiEntry = pixel_to_wi.find(pixelID);
339 if (wiEntry != pixel_to_wi.end()) {
340 size_t wi = wiEntry->second;
341 // Set detectorIDs
342 outputW->getSpectrum(idet).addDetectorIDs(inputW->getSpectrum(wi).getDetectorIDs());
343 }
344
345 return TOFmax - 1;
346}
347
348} // namespace Mantid::Crystal
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
#define fabs(x)
Definition: Matrix.cpp:22
#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_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.
#define UNUSED_ARG(x)
Function arguments are sometimes unused in certain implmentations but are required for documentation ...
Definition: System.h:64
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
A property class for workspaces.
bool m_IC
Ikeida Carpenter fit of TOF.
void retrieveProperties()
Read in all the input parameters.
int fitneighbours(int ipeak, const std::string &det_name, int x0, int y0, int idet, double qspan, DataObjects::PeaksWorkspace_sptr &Peaks, const detid2index_map &pixel_to_wi)
API::MatrixWorkspace_sptr outputW
A pointer to the output workspace.
void exec() override
Executes the algorithm.
void init() override
Initialisation method.
API::MatrixWorkspace_sptr inputW
A pointer to the input workspace.
void setSigmaIntensity(double m_sigmaIntensity) override
Set the error on the integrated peak intensity.
Definition: BasePeak.cpp:206
int getRunNumber() const override
Return the run number this peak was measured at.
Definition: BasePeak.cpp:78
void setIntensity(double m_intensity) override
Set the integrated peak intensity.
Definition: BasePeak.cpp:198
Structure describing a single-crystal peak.
Definition: Peak.h:34
int getCol() const override
For RectangularDetectors only, returns the column (x) of the pixel of the detector or -1 if not found...
Definition: Peak.cpp:341
int getDetectorID() const
Get the ID of the detector at the center of the peak
Definition: Peak.cpp:271
Class to implement UB matrix.
double a(int nd) const
Get lattice parameter a1-a3 as function of index (0-2)
Definition: UnitCell.cpp:94
double c() const
Get lattice parameter.
Definition: UnitCell.cpp:128
double b() const
Get lattice parameter.
Definition: UnitCell.cpp:123
Exception for when an item is not found in a collection.
Definition: Exception.h:145
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void error(const std::string &msg)
Logs at error level.
Definition: Logger.cpp:77
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< IAlgorithm > IAlgorithm_sptr
shared pointer to Mantid::API::IAlgorithm
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
std::shared_ptr< const EventWorkspace > EventWorkspace_const_sptr
shared pointer to a const Workspace2D
std::shared_ptr< TableWorkspace > TableWorkspace_sptr
shared pointer to Mantid::DataObjects::TableWorkspace
std::shared_ptr< PeaksWorkspace > PeaksWorkspace_sptr
Typedef for a shared pointer to a peaks workspace.
std::shared_ptr< const Instrument > Instrument_const_sptr
Shared pointer to an const instrument object.
MANTID_KERNEL_DLL int getBinIndex(const std::vector< double > &bins, const double value)
Return the index into a vector of bin boundaries for a particular X value.
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
std::unordered_map< detid_t, size_t > detid2index_map
Map with key = detector ID, value = workspace index.
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54