Mantid
Loading...
Searching...
No Matches
IntegratePeaksMD.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 +
9#include "MantidAPI/Column.h"
16#include "MantidAPI/Run.h"
17#include "MantidAPI/TableRow.h"
18#include "MantidAPI/TextAxis.h"
25#include "MantidHistogramData/LinearGenerator.h"
27#include "MantidKernel/System.h"
28#include "MantidKernel/Utils.h"
30
31#include <algorithm>
32#include <cmath>
33#include <fstream>
34#include <gsl/gsl_integration.h>
35
36namespace Mantid::MDAlgorithms {
37
38// Register the algorithm into the AlgorithmFactory
39DECLARE_ALGORITHM(IntegratePeaksMD)
40
41using namespace Mantid::Kernel;
42using namespace Mantid::API;
43using namespace Mantid::DataObjects;
44using namespace Mantid::DataObjects;
45using namespace Mantid::Geometry;
46using namespace Mantid::HistogramData;
47
51 declareProperty(std::make_unique<WorkspaceProperty<IMDEventWorkspace>>("InputWorkspace", "", Direction::Input),
52 "An input MDEventWorkspace.");
53
54 std::vector<std::string> propOptions{"Q (lab frame)", "Q (sample frame)", "HKL"};
55 declareProperty("CoordinatesToUse", "Q (lab frame)", std::make_shared<StringListValidator>(propOptions),
56 "Ignored: algorithm uses the InputWorkspace's coordinates.");
57
58 declareProperty(std::make_unique<PropertyWithValue<double>>("PeakRadius", 1.0, Direction::Input),
59 "Fixed radius around each peak position in which to integrate (in the "
60 "same units as the workspace).");
61
62 declareProperty(std::make_unique<PropertyWithValue<double>>("BackgroundInnerRadius", 0.0, Direction::Input),
63 "Inner radius to use to evaluate the background of the peak.\n"
64 "If smaller than PeakRadius, then we assume BackgroundInnerRadius = "
65 "PeakRadius.");
66
67 declareProperty(std::make_unique<PropertyWithValue<double>>("BackgroundOuterRadius", 0.0, Direction::Input),
68 "Outer radius to use to evaluate the background of the peak.\n"
69 "The signal density around the peak (BackgroundInnerRadius < r < "
70 "BackgroundOuterRadius) is used to estimate the background under the "
71 "peak.\n"
72 "If smaller than PeakRadius, no background measurement is done.");
73
74 declareProperty(std::make_unique<WorkspaceProperty<PeaksWorkspace>>("PeaksWorkspace", "", Direction::Input),
75 "A PeaksWorkspace containing the peaks to integrate.");
76
77 declareProperty(std::make_unique<WorkspaceProperty<PeaksWorkspace>>("OutputWorkspace", "", Direction::Output),
78 "The output PeaksWorkspace will be a copy of the input PeaksWorkspace "
79 "with the peaks' integrated intensities.");
80
81 declareProperty("ReplaceIntensity", true,
82 "Always replace intensity in PeaksWorkspacem (default).\n"
83 "If false, then do not replace intensity if calculated value "
84 "is 0 (used for SNSSingleCrystalReduction)");
85
86 declareProperty("IntegrateIfOnEdge", true,
87 "Only warning if all of peak outer radius is not on detector (default).\n"
88 "If false, do not integrate if the outer radius is not on a detector.");
89
90 declareProperty("AdaptiveQRadius", false,
91 "Default is false. If true, all input radii are multiplied "
92 "by the magnitude of Q at the peak center so each peak has a "
93 "different integration radius.");
94
95 declareProperty("Cylinder", false, "Default is sphere. Use next five parameters for cylinder.");
96
97 declareProperty(std::make_unique<PropertyWithValue<double>>("CylinderLength", 0.0, Direction::Input),
98 "Length of cylinder in which to integrate (in the same units as the "
99 "workspace).");
100
101 declareProperty(std::make_unique<PropertyWithValue<double>>("PercentBackground", 0.0, Direction::Input),
102 "Percent of CylinderLength that is background (20 is 20%)");
103
104 std::vector<std::string> peakNames = FunctionFactory::Instance().getFunctionNames<IPeakFunction>();
105 peakNames.emplace_back("NoFit");
106 declareProperty("ProfileFunction", "Gaussian", std::make_shared<StringListValidator>(peakNames),
107 "Fitting function for profile that is used only with "
108 "Cylinder integration.");
109
110 std::vector<std::string> integrationOptions(2);
111 integrationOptions[0] = "Sum";
112 integrationOptions[1] = "GaussianQuadrature";
113 auto integrationvalidator = std::make_shared<StringListValidator>(integrationOptions);
114 declareProperty("IntegrationOption", "GaussianQuadrature", integrationvalidator,
115 "Integration method for calculating intensity "
116 "used only with Cylinder integration.");
117
118 declareProperty(std::make_unique<FileProperty>("ProfilesFile", "", FileProperty::OptionalSave,
119 std::vector<std::string>(1, "profiles")),
120 "Save (Optionally) as Isaw peaks file with profiles included");
121}
122
123//----------------------------------------------------------------------------------------------
128template <typename MDE, size_t nd> void IntegratePeaksMD::integrate(typename MDEventWorkspace<MDE, nd>::sptr ws) {
129 if (nd != 3)
130 throw std::invalid_argument("For now, we expect the input MDEventWorkspace "
131 "to have 3 dimensions only.");
132
134 Mantid::DataObjects::PeaksWorkspace_sptr inPeakWS = getProperty("PeaksWorkspace");
135
137 Mantid::DataObjects::PeaksWorkspace_sptr peakWS = getProperty("OutputWorkspace");
138 if (peakWS != inPeakWS)
139 peakWS = inPeakWS->clone();
140
142 std::string CoordinatesToUseStr = getPropertyValue("CoordinatesToUse");
144 g_log.warning() << " Warning" << CoordinatesToUse << '\n';
145 if (CoordinatesToUse == Kernel::QLab && CoordinatesToUseStr != "Q (lab frame)")
146 g_log.warning() << "Warning: used Q (lab frame) coordinates for MD "
147 "workspace, not CoordinatesToUse from input \n";
148 else if (CoordinatesToUse == Kernel::QSample && CoordinatesToUseStr != "Q (sample frame)")
149 g_log.warning() << "Warning: used Q (sample frame) coordinates for MD "
150 "workspace, not CoordinatesToUse from input \n";
151 else if (CoordinatesToUse == Kernel::HKL && CoordinatesToUseStr != "HKL")
152 g_log.warning() << "Warning: used HKL coordinates for MD workspace, not "
153 "CoordinatesToUse from input \n";
154
156 double PeakRadius = getProperty("PeakRadius");
158 double BackgroundOuterRadius = getProperty("BackgroundOuterRadius");
160 double BackgroundInnerRadius = getProperty("BackgroundInnerRadius");
161 if (BackgroundInnerRadius < PeakRadius)
162 BackgroundInnerRadius = PeakRadius;
164 double cylinderLength = getProperty("CylinderLength");
165 Workspace2D_sptr wsProfile2D, wsFit2D, wsDiff2D;
166 size_t numSteps = 0;
167 bool cylinderBool = getProperty("Cylinder");
168 bool adaptiveQRadius = getProperty("AdaptiveQRadius");
169 std::vector<double> PeakRadiusVector(peakWS->getNumberPeaks(), PeakRadius);
170 std::vector<double> BackgroundInnerRadiusVector(peakWS->getNumberPeaks(), BackgroundInnerRadius);
171 std::vector<double> BackgroundOuterRadiusVector(peakWS->getNumberPeaks(), BackgroundOuterRadius);
172 if (cylinderBool) {
173 numSteps = 100;
174 size_t histogramNumber = peakWS->getNumberPeaks();
175 Workspace_sptr wsProfile = WorkspaceFactory::Instance().create("Workspace2D", histogramNumber, numSteps, numSteps);
176 wsProfile2D = std::dynamic_pointer_cast<Workspace2D>(wsProfile);
177 AnalysisDataService::Instance().addOrReplace("ProfilesData", wsProfile2D);
178 Workspace_sptr wsFit = WorkspaceFactory::Instance().create("Workspace2D", histogramNumber, numSteps, numSteps);
179 wsFit2D = std::dynamic_pointer_cast<Workspace2D>(wsFit);
180 AnalysisDataService::Instance().addOrReplace("ProfilesFit", wsFit2D);
181 Workspace_sptr wsDiff = WorkspaceFactory::Instance().create("Workspace2D", histogramNumber, numSteps, numSteps);
182 wsDiff2D = std::dynamic_pointer_cast<Workspace2D>(wsDiff);
183 AnalysisDataService::Instance().addOrReplace("ProfilesFitDiff", wsDiff2D);
184 auto newAxis1 = std::make_unique<TextAxis>(peakWS->getNumberPeaks());
185 auto newAxis2 = std::make_unique<TextAxis>(peakWS->getNumberPeaks());
186 auto newAxis3 = std::make_unique<TextAxis>(peakWS->getNumberPeaks());
187 auto newAxis1Raw = newAxis1.get();
188 auto newAxis2Raw = newAxis2.get();
189 auto newAxis3Raw = newAxis3.get();
190 wsProfile2D->replaceAxis(1, std::move(newAxis1));
191 wsFit2D->replaceAxis(1, std::move(newAxis2));
192 wsDiff2D->replaceAxis(1, std::move(newAxis3));
193 for (int i = 0; i < peakWS->getNumberPeaks(); ++i) {
194 // Get a direct ref to that peak.
195 IPeak &p = peakWS->getPeak(i);
196 std::ostringstream label;
197 label << Utils::round(p.getH()) << "_" << Utils::round(p.getK()) << "_" << Utils::round(p.getL()) << "_"
198 << p.getRunNumber();
199 newAxis1Raw->setLabel(i, label.str());
200 newAxis2Raw->setLabel(i, label.str());
201 newAxis3Raw->setLabel(i, label.str());
202 }
203 }
204 double percentBackground = getProperty("PercentBackground");
205 size_t peakMin = 0;
206 size_t peakMax = numSteps;
207 double ratio = 0.0;
208 if (cylinderBool) {
209 peakMin = static_cast<size_t>(static_cast<double>(numSteps) * percentBackground / 100.);
210 peakMax = numSteps - peakMin - 1;
211 size_t numPeakCh = peakMax - peakMin + 1; // number of peak channels
212 size_t numBkgCh = numSteps - numPeakCh; // number of background channels
213 ratio = static_cast<double>(numPeakCh) / static_cast<double>(numBkgCh);
214 }
216 bool replaceIntensity = getProperty("ReplaceIntensity");
217 bool integrateEdge = getProperty("IntegrateIfOnEdge");
218
219 std::string profileFunction = getProperty("ProfileFunction");
220 std::string integrationOption = getProperty("IntegrationOption");
221 std::ofstream out;
222 if (cylinderBool && profileFunction != "NoFit") {
223 std::string outFile = getProperty("InputWorkspace");
224 outFile.append(profileFunction);
225 outFile.append(".dat");
226 std::string save_path = ConfigService::Instance().getString("defaultsave.directory");
227 outFile = save_path + outFile;
228 out.open(outFile.c_str(), std::ofstream::out);
229 }
230 //
231 // If the following OMP pragma is included, this algorithm seg faults
232 // sporadically when processing multiple TOPAZ runs in a script, on
233 // Scientific Linux 6.2. Typically, it seg faults after 2 to 6 runs are
234 // processed, though occasionally it will process all 8 requested in the
235 // script without crashing. Since the lower level codes already use OpenMP,
236 // parallelizing at this level is only marginally useful, giving about a
237 // 5-10% speedup. Perhaps is should just be removed permanantly, but for
238 // now it is commented out to avoid the seg faults. Refs #5533
239 // PRAGMA_OMP(parallel for schedule(dynamic, 10) )
240 for (int i = 0; i < peakWS->getNumberPeaks(); ++i) {
241 // Get a direct ref to that peak.
242 IPeak &p = peakWS->getPeak(i);
243
244 // Get the peak center as a position in the dimensions of the workspace
245 V3D pos;
246 if (CoordinatesToUse == Kernel::QLab) //"Q (lab frame)"
247 pos = p.getQLabFrame();
248 else if (CoordinatesToUse == Kernel::QSample) //"Q (sample frame)"
249 pos = p.getQSampleFrame();
250 else if (CoordinatesToUse == Kernel::HKL) //"HKL"
251 pos = p.getHKL();
252
253 // Get the instrument and its detectors
254 inst = peakWS->getInstrument();
255 // Do not integrate if sphere is off edge of detector
256 if (BackgroundOuterRadius > PeakRadius) {
257 if (!detectorQ(p.getQLabFrame(), BackgroundOuterRadius)) {
258 g_log.warning() << "Warning: sphere/cylinder for integration is off "
259 "edge of detector for peak "
260 << i << '\n';
261 if (!integrateEdge)
262 continue;
263 }
264 } else {
265 if (!detectorQ(p.getQLabFrame(), PeakRadius)) {
266 g_log.warning() << "Warning: sphere/cylinder for integration is off "
267 "edge of detector for peak "
268 << i << '\n';
269 if (!integrateEdge)
270 continue;
271 }
272 }
273
274 // Build the sphere transformation
275 bool dimensionsUsed[nd];
276 coord_t center[nd];
277 for (size_t d = 0; d < nd; ++d) {
278 dimensionsUsed[d] = true; // Use all dimensions
279 center[d] = static_cast<coord_t>(pos[d]);
280 }
281 signal_t signal = 0;
282 signal_t errorSquared = 0;
283 signal_t bgSignal = 0;
284 signal_t bgErrorSquared = 0;
285 double background_total = 0.0;
286 if (!cylinderBool) {
287 // modulus of Q
288 coord_t lenQpeak = 1.0;
289 if (adaptiveQRadius) {
290 lenQpeak = 0.0;
291 for (size_t d = 0; d < nd; d++) {
292 lenQpeak += center[d] * center[d];
293 }
294 lenQpeak = std::sqrt(lenQpeak);
295 }
296 PeakRadiusVector[i] = lenQpeak * PeakRadius;
297 BackgroundInnerRadiusVector[i] = lenQpeak * BackgroundInnerRadius;
298 BackgroundOuterRadiusVector[i] = lenQpeak * BackgroundOuterRadius;
299 CoordTransformDistance sphere(nd, center, dimensionsUsed);
300
301 if (auto *shapeablePeak = dynamic_cast<Peak *>(&p)) {
302
303 PeakShape *sphereShape =
304 new PeakShapeSpherical(PeakRadiusVector[i], BackgroundInnerRadiusVector[i], BackgroundOuterRadiusVector[i],
305 CoordinatesToUse, this->name(), this->version());
306 shapeablePeak->setPeakShape(sphereShape);
307 }
308
309 // Perform the integration into whatever box is contained within.
310 ws->getBox()->integrateSphere(sphere, static_cast<coord_t>(lenQpeak * PeakRadius * lenQpeak * PeakRadius), signal,
311 errorSquared);
312
313 // Integrate around the background radius
314
315 if (BackgroundOuterRadius > PeakRadius) {
316 // Get the total signal inside "BackgroundOuterRadius"
317 ws->getBox()->integrateSphere(
318 sphere, static_cast<coord_t>(lenQpeak * BackgroundOuterRadius * lenQpeak * BackgroundOuterRadius), bgSignal,
319 bgErrorSquared);
320
321 // Evaluate the signal inside "BackgroundInnerRadius"
322 signal_t interiorSignal = 0;
323 signal_t interiorErrorSquared = 0;
324
325 // Integrate this 3rd radius, if needed
326 if (BackgroundInnerRadius != PeakRadius) {
327 ws->getBox()->integrateSphere(
328 sphere, static_cast<coord_t>(lenQpeak * BackgroundInnerRadius * lenQpeak * BackgroundInnerRadius),
329 interiorSignal, interiorErrorSquared);
330 } else {
331 // PeakRadius == BackgroundInnerRadius, so use the previous value
332 interiorSignal = signal;
333 interiorErrorSquared = errorSquared;
334 }
335 // Subtract the peak part to get the intensity in the shell
336 // (BackgroundInnerRadius < r < BackgroundOuterRadius)
337 bgSignal -= interiorSignal;
338 // We can subtract the error (instead of adding) because the two values
339 // are 100% dependent; this is the same as integrating a shell.
340 bgErrorSquared -= interiorErrorSquared;
341
342 // Relative volume of peak vs the BackgroundOuterRadius sphere
343 const double radiusRatio = (PeakRadius / BackgroundOuterRadius);
344 const double peakVolume = radiusRatio * radiusRatio * radiusRatio;
345
346 // Relative volume of the interior of the shell vs overall background
347 const double interiorRatio = (BackgroundInnerRadius / BackgroundOuterRadius);
348 // Volume of the bg shell, relative to the volume of the
349 // BackgroundOuterRadius sphere
350 const double bgVolume = 1.0 - interiorRatio * interiorRatio * interiorRatio;
351
352 // Finally, you will multiply the bg intensity by this to get the
353 // estimated background under the peak volume
354 const double scaleFactor = peakVolume / bgVolume;
355 bgSignal *= scaleFactor;
356 bgErrorSquared *= scaleFactor * scaleFactor;
357 }
358 } else {
359 CoordTransformDistance cylinder(nd, center, dimensionsUsed, 2);
360
361 // Perform the integration into whatever box is contained within.
362 Counts signal_fit(numSteps);
363 signal_fit = 0;
364
365 ws->getBox()->integrateCylinder(cylinder, static_cast<coord_t>(PeakRadius), static_cast<coord_t>(cylinderLength),
366 signal, errorSquared, signal_fit.mutableRawData());
367
368 Points points(signal_fit.size(), LinearGenerator(0, 1));
369 wsProfile2D->setHistogram(i, points, signal_fit);
370
371 // Integrate around the background radius
372 if (BackgroundOuterRadius > PeakRadius) {
373 // Get the total signal inside "BackgroundOuterRadius"
374 signal_fit = 0;
375
376 ws->getBox()->integrateCylinder(cylinder, static_cast<coord_t>(BackgroundOuterRadius),
377 static_cast<coord_t>(cylinderLength), bgSignal, bgErrorSquared,
378 signal_fit.mutableRawData());
379
380 wsProfile2D->setHistogram(i, points, signal_fit);
381
382 // Evaluate the signal inside "BackgroundInnerRadius"
383 signal_t interiorSignal = 0;
384 signal_t interiorErrorSquared = 0;
385
386 // Integrate this 3rd radius, if needed
387 if (BackgroundInnerRadius != PeakRadius) {
388 ws->getBox()->integrateCylinder(cylinder, static_cast<coord_t>(BackgroundInnerRadius),
389 static_cast<coord_t>(cylinderLength), interiorSignal, interiorErrorSquared,
390 signal_fit.mutableRawData());
391 } else {
392 // PeakRadius == BackgroundInnerRadius, so use the previous value
393 interiorSignal = signal;
394 interiorErrorSquared = errorSquared;
395 }
396 // Subtract the peak part to get the intensity in the shell
397 // (BackgroundInnerRadius < r < BackgroundOuterRadius)
398 bgSignal -= interiorSignal;
399 // We can subtract the error (instead of adding) because the two values
400 // are 100% dependent; this is the same as integrating a shell.
401 bgErrorSquared -= interiorErrorSquared;
402 // Relative volume of peak vs the BackgroundOuterRadius cylinder
403 const double radiusRatio = (PeakRadius / BackgroundOuterRadius);
404 const double peakVolume = radiusRatio * radiusRatio * cylinderLength;
405
406 // Relative volume of the interior of the shell vs overall background
407 const double interiorRatio = (BackgroundInnerRadius / BackgroundOuterRadius);
408 // Volume of the bg shell, relative to the volume of the
409 // BackgroundOuterRadius cylinder
410 const double bgVolume = 1.0 - interiorRatio * interiorRatio * cylinderLength;
411
412 // Finally, you will multiply the bg intensity by this to get the
413 // estimated background under the peak volume
414 const double scaleFactor = peakVolume / bgVolume;
415 bgSignal *= scaleFactor;
416 bgErrorSquared *= scaleFactor * scaleFactor;
417 } else {
418 wsProfile2D->setHistogram(i, points, signal_fit);
419 }
420
421 if (profileFunction == "NoFit") {
422 auto &y = wsProfile2D->y(i);
423 // sum signal between range
424 signal = y.sum(peakMin, peakMax);
425 // sum background outside of range
426 background_total += y.sum(0, peakMin);
427 background_total += y.sum(peakMax);
428 errorSquared = std::abs(signal);
429 } else {
430 auto findpeaks = createChildAlgorithm("FindPeaks", -1, -1, false);
431 findpeaks->setProperty("InputWorkspace", wsProfile2D);
432 findpeaks->setProperty<int>("FWHM", 7);
433 findpeaks->setProperty<int>("Tolerance", 4);
434 // FindPeaks will do the checking on the validity of WorkspaceIndex
435 findpeaks->setProperty("WorkspaceIndex", static_cast<int>(i));
436
437 // Get the specified peak positions, which is optional
438 findpeaks->setProperty<std::string>("PeakFunction", profileFunction);
439 // FindPeaks will use linear or flat if they are better
440 findpeaks->setProperty<std::string>("BackgroundType", "Quadratic");
441 findpeaks->setProperty<bool>("HighBackground", true);
442 findpeaks->setProperty<bool>("RawPeakParameters", true);
443 std::vector<double> peakPosToFit;
444 peakPosToFit.emplace_back(static_cast<double>(numSteps) / 2.0);
445 findpeaks->setProperty("PeakPositions", peakPosToFit);
446 findpeaks->setProperty<int>("MinGuessedPeakWidth", 4);
447 findpeaks->setProperty<int>("MaxGuessedPeakWidth", 4);
448 try {
449 findpeaks->executeAsChildAlg();
450 } catch (...) {
451 g_log.error("Can't execute FindPeaks algorithm");
452 continue;
453 }
454
455 API::ITableWorkspace_sptr paramws = findpeaks->getProperty("PeaksList");
456 if (paramws->rowCount() < 1)
457 continue;
458 std::ostringstream fun_str;
459 fun_str << "name=" << profileFunction;
460
461 size_t numcols = paramws->columnCount();
462 std::vector<std::string> paramsName = paramws->getColumnNames();
463 std::vector<double> paramsValue;
464 API::TableRow row = paramws->getRow(0);
465 int spectrum;
466 row >> spectrum;
467 for (size_t j = 1; j < numcols; ++j) {
468 double parvalue;
469 row >> parvalue;
470 if (j == numcols - 4)
471 fun_str << ";name=Quadratic";
472 // erase f0. or f1.
473 // if (j > 0 && j < numcols-1) fun_str << "," <<
474 // paramsName[j].erase(0,3) <<"="<<parvalue;
475 if (j > 0 && j < numcols - 1)
476 fun_str << "," << paramsName[j] << "=" << parvalue;
477 paramsValue.emplace_back(parvalue);
478 }
479 if (i == 0) {
480 for (size_t j = 0; j < numcols; ++j)
481 out << std::setw(20) << paramsName[j] << " ";
482 out << "\n";
483 }
484 out << std::setw(20) << i;
485 for (size_t j = 0; j < numcols - 1; ++j)
486 out << std::setw(20) << std::fixed << std::setprecision(10) << paramsValue[j] << " ";
487 out << "\n";
488
489 // Evaluate fit at points
490
491 IFunction_sptr ifun = FunctionFactory::Instance().createInitialized(fun_str.str());
492 std::shared_ptr<const CompositeFunction> fun = std::dynamic_pointer_cast<const CompositeFunction>(ifun);
493 const auto &x = wsProfile2D->x(i);
494 wsFit2D->setSharedX(i, wsProfile2D->sharedX(i));
495 wsDiff2D->setSharedX(i, wsProfile2D->sharedX(i));
496
497 FunctionDomain1DVector domain(x.rawData());
498 FunctionValues yy(domain);
499 fun->function(domain, yy);
500 auto funcValues = yy.toVector();
501
502 wsFit2D->mutableY(i) = funcValues;
503 wsDiff2D->setSharedY(i, wsProfile2D->sharedY(i));
504 wsDiff2D->mutableY(i) -= wsFit2D->y(i);
505
506 // Calculate intensity
507 signal = 0.0;
508 if (integrationOption == "Sum") {
509
510 for (size_t j = peakMin; j <= peakMax; j++)
511 if (std::isfinite(yy[j]))
512 signal += yy[j];
513 } else {
514 gsl_integration_workspace *w = gsl_integration_workspace_alloc(1000);
515
516 double error;
517
518 gsl_function F;
519 F.function = &Mantid::MDAlgorithms::f_eval;
520 F.params = &fun;
521
522 gsl_integration_qags(&F, x[peakMin], x[peakMax], 0, 1e-7, 1000, w, &signal, &error);
523
524 gsl_integration_workspace_free(w);
525 }
526 errorSquared = std::fabs(signal);
527 // Get background counts
528 for (size_t j = 0; j < numSteps; j++) {
529 double background =
530 paramsValue[numcols - 3] * x[j] * x[j] + paramsValue[numcols - 4] * x[j] + paramsValue[numcols - 5];
531 if (j < peakMin || j > peakMax)
532 background_total = background_total + background;
533 }
534 }
535 }
536 // Save it back in the peak object.
537 if (signal != 0. || replaceIntensity) {
538 p.setIntensity(signal - ratio * background_total - bgSignal);
539 p.setSigmaIntensity(sqrt(errorSquared + ratio * ratio * std::fabs(background_total) + bgErrorSquared));
540 }
541
542 g_log.information() << "Peak " << i << " at " << pos << ": signal " << signal << " (sig^2 " << errorSquared
543 << "), with background " << bgSignal + ratio * background_total << " (sig^2 "
544 << bgErrorSquared + ratio * ratio * std::fabs(background_total) << ") subtracted.\n";
545 }
546 // This flag is used by the PeaksWorkspace to evaluate whether it has been
547 // integrated.
548 peakWS->mutableRun().addProperty("PeaksIntegrated", 1, true);
549 // These flags are specific to the algorithm.
550 peakWS->mutableRun().addProperty("PeakRadius", PeakRadiusVector, true);
551 peakWS->mutableRun().addProperty("BackgroundInnerRadius", BackgroundInnerRadiusVector, true);
552 peakWS->mutableRun().addProperty("BackgroundOuterRadius", BackgroundOuterRadiusVector, true);
553
554 // save profiles in peaks file
555 const std::string outfile = getProperty("ProfilesFile");
556 if (outfile.length() > 0) {
557 IAlgorithm_sptr alg;
558 try {
559 alg = createChildAlgorithm("SaveIsawPeaks", -1, -1, false);
560 } catch (Exception::NotFoundError &) {
561 g_log.error("Can't locate SaveIsawPeaks algorithm");
562 throw;
563 }
564 alg->setProperty("InputWorkspace", peakWS);
565 alg->setProperty("ProfileWorkspace", wsProfile2D);
566 alg->setPropertyValue("Filename", outfile);
567 alg->execute();
568 }
569 // Save the output
570 setProperty("OutputWorkspace", peakWS);
571}
572
578bool IntegratePeaksMD::detectorQ(const Mantid::Kernel::V3D &QLabFrame, double r) {
579 bool in = true;
580 const int nAngles = 8;
581 double dAngles = static_cast<coord_t>(nAngles);
582 // check 64 points in theta and phi at outer radius
583 for (int i = 0; i < nAngles; ++i) {
584 double theta = (2 * M_PI) / dAngles * i;
585 for (int j = 0; j < nAngles; ++j) {
586 double phi = (2 * M_PI) / dAngles * j;
587 // Calculate an edge position at this point on the sphere surface.
588 // Spherical coordinates to cartesian.
589 V3D edge = V3D(QLabFrame.X() + r * std::cos(theta) * std::sin(phi),
590 QLabFrame.Y() + r * std::sin(theta) * std::sin(phi), QLabFrame.Z() + r * std::cos(phi));
591 // Create the peak using the Q in the lab frame with all its info:
592 try {
593 Peak p(inst, edge);
594 in = (in && p.findDetector());
595 if (!in) {
596 return false;
597 }
598 } catch (...) {
599 return false;
600 }
601 }
602 }
603 return in;
604}
605
606//----------------------------------------------------------------------------------------------
610 inWS = getProperty("InputWorkspace");
611
613}
614
615double f_eval(double x, void *params) {
616 std::shared_ptr<const API::CompositeFunction> fun =
617 *reinterpret_cast<std::shared_ptr<const API::CompositeFunction> *>(params);
619 FunctionValues yval(domain);
620 fun->function(domain, yval);
621 return yval[0];
622}
623
624} // namespace Mantid::MDAlgorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
double error
Definition: IndexPeaks.cpp:133
#define CALL_MDEVENT_FUNCTION(funcname, workspace)
Macro that makes it possible to call a templated method for a MDEventWorkspace using a IMDEventWorksp...
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
std::string getPropertyValue(const std::string &name) const override
Get the value of a property as a string.
Definition: Algorithm.cpp:2026
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
@ OptionalSave
to specify a file to write to but an empty string is
Definition: FileProperty.h:50
Implements FunctionDomain1D with its own storage in form of a std::vector.
A class to store values calculated by a function.
std::vector< double > toVector() const
Return the calculated values as a vector.
An interface to a peak function, which extend the interface of IFunctionWithLocation by adding method...
Definition: IPeakFunction.h:51
TableRow represents a row in a TableWorkspace.
Definition: TableRow.h:39
A property class for workspaces.
Unique CoordCenterVectorParam type declaration for ndimensional coordinate centers.
void integrateSphere(Mantid::API::CoordTransform &radiusTransform, const coord_t radiusSquared, signal_t &signal, signal_t &errorSquared, const coord_t innerRadiusSquared=0.0, const bool useOnePercentBackgroundCorrection=true) const override=0
Sphere (peak) integration.
void integrateCylinder(Mantid::API::CoordTransform &radiusTransform, const coord_t radius, const coord_t length, signal_t &signal, signal_t &errorSquared, std::vector< signal_t > &signal_fit) const override=0
Cylinder (peak) integration.
std::shared_ptr< MDEventWorkspace< MDE, nd > > sptr
Typedef for a shared pointer of this kind of event workspace.
Kernel::SpecialCoordinateSystem getSpecialCoordinateSystem() const override
Get the coordinate system.
PeakShapeSpherical : PeakShape for a spherical peak.
Structure describing a single-crystal peak.
Definition: Peak.h:34
bool findDetector()
After creating a peak using the Q in the lab frame, the detPos is set to the direction of the detecto...
Definition: Peak.cpp:613
Structure describing a single-crystal peak.
Definition: IPeak.h:26
virtual double getH() const =0
virtual void setIntensity(double m_Intensity)=0
virtual double getK() const =0
virtual int getRunNumber() const =0
virtual Mantid::Kernel::V3D getQSampleFrame() const =0
virtual double getL() const =0
virtual void setSigmaIntensity(double m_SigmaIntensity)=0
virtual Mantid::Kernel::V3D getQLabFrame() const =0
virtual Mantid::Kernel::V3D getHKL() const =0
PeakShape : Abstract type to describes the shape of a peak.
Definition: PeakShape.h:20
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 warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
The concrete, templated class for properties.
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
Class for 3D vectors.
Definition: V3D.h:34
constexpr double X() const noexcept
Get x.
Definition: V3D.h:232
constexpr double Y() const noexcept
Get y.
Definition: V3D.h:233
constexpr double Z() const noexcept
Get z.
Definition: V3D.h:234
void integrate(typename DataObjects::MDEventWorkspace< MDE, nd >::sptr ws)
Integrate the peaks of the workspace using parameters saved in the algorithm class.
int version() const override
Algorithm's version for identification.
bool detectorQ(const Mantid::Kernel::V3D &QLabFrame, double r)
Calculate if this Q is on a detector.
void exec() override
Run the algorithm.
Mantid::API::IMDEventWorkspace_sptr inWS
Input MDEventWorkspace.
void init() override
Initialise the properties.
Geometry::Instrument_const_sptr inst
Instrument reference.
const std::string name() const override
Algorithm's name for identification.
std::shared_ptr< IAlgorithm > IAlgorithm_sptr
shared pointer to Mantid::API::IAlgorithm
std::shared_ptr< ITableWorkspace > ITableWorkspace_sptr
shared pointer to Mantid::API::ITableWorkspace
std::shared_ptr< Workspace > Workspace_sptr
shared pointer to Mantid::API::Workspace
Definition: Workspace_fwd.h:20
std::shared_ptr< IFunction > IFunction_sptr
shared pointer to the function base class
Definition: IFunction.h:732
std::shared_ptr< Workspace2D > Workspace2D_sptr
shared pointer to Mantid::DataObjects::Workspace2D
std::shared_ptr< PeaksWorkspace > PeaksWorkspace_sptr
Typedef for a shared pointer to a peaks workspace.
long round(double x)
Custom rounding method for a double->long because none is portable in C++ (!)
Definition: Utils.h:37
SpecialCoordinateSystem
Special coordinate systems for Q3D.
double f_eval(double x, void *params)
float coord_t
Typedef for the data type to use for coordinate axes in MD objects such as MDBox, MDEventWorkspace,...
Definition: MDTypes.h:27
double signal_t
Typedef for the signal recorded in a MDBox, etc.
Definition: MDTypes.h:36
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54