Mantid
Loading...
Searching...
No Matches
IntegratePeaksMD2.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/Progress.h"
17#include "MantidAPI/Run.h"
18#include "MantidAPI/TableRow.h"
19#include "MantidAPI/TextAxis.h"
31#include "MantidHistogramData/LinearGenerator.h"
37#include "MantidKernel/System.h"
38#include "MantidKernel/Utils.h"
42
43#include "boost/math/distributions.hpp"
44
45#include <cmath>
46#include <fstream>
47#include <gsl/gsl_integration.h>
48
49namespace Mantid::MDAlgorithms {
50
51// Register the algorithm into the AlgorithmFactory
52DECLARE_ALGORITHM(IntegratePeaksMD2)
53
54using namespace Mantid::HistogramData;
55using namespace Mantid::Kernel;
56using namespace Mantid::API;
57using namespace Mantid::DataObjects;
58using namespace Mantid::DataObjects;
59using namespace Mantid::Geometry;
60
64 declareProperty(std::make_unique<WorkspaceProperty<IMDEventWorkspace>>("InputWorkspace", "", Direction::Input),
65 "An input MDEventWorkspace.");
66
67 auto radiiValidator = std::make_shared<ArrayBoundedValidator<double>>();
68 radiiValidator->setLower(0.0);
69 radiiValidator->setLowerExclusive(true);
70 declareProperty(std::make_unique<ArrayProperty<double>>("PeakRadius", std::vector<double>({1.0}), radiiValidator,
72 "Fixed radius around each peak position in which to integrate, or the "
73 "semi-axis lengths (a,b,c) describing an ellipsoid shape used for "
74 "integration (in the same units as the workspace).");
75
76 radiiValidator->setLowerExclusive(false);
77 declareProperty(std::make_unique<ArrayProperty<double>>("BackgroundInnerRadius", std::vector<double>({0.0}),
78 radiiValidator, Direction::Input),
79 "Inner radius, or three values for semi-axis lengths (a,b,c) of the "
80 "ellipsoid shape, used to evaluate the background of the peak.\n"
81 "If smaller than PeakRadius, then we assume BackgroundInnerRadius = "
82 "PeakRadius.");
83
84 declareProperty(std::make_unique<ArrayProperty<double>>("BackgroundOuterRadius", std::vector<double>({0.0}),
85 radiiValidator, Direction::Input),
86 "Outer radius, or three values for semi-axis lengths (a,b,c) of the "
87 "ellipsoid shape, to use to evaluate the background of the peak.\n"
88 "The signal density around the peak (BackgroundInnerRadius < r < "
89 "BackgroundOuterRadius) is used to estimate the background under the "
90 "peak.\n"
91 "If smaller than PeakRadius, no background measurement is done.");
92
93 declareProperty(std::make_unique<WorkspaceProperty<IPeaksWorkspace>>("PeaksWorkspace", "", Direction::Input),
94 "A PeaksWorkspace containing the peaks to integrate.");
95
96 declareProperty(std::make_unique<WorkspaceProperty<IPeaksWorkspace>>("OutputWorkspace", "", Direction::Output),
97 "The output PeaksWorkspace will be a copy of the input PeaksWorkspace "
98 "with the peaks' integrated intensities.");
99
100 declareProperty("ReplaceIntensity", true,
101 "Always replace intensity in PeaksWorkspacem (default).\n"
102 "If false, then do not replace intensity if calculated value "
103 "is 0 (used for SNSSingleCrystalReduction)");
104
105 declareProperty("IntegrateIfOnEdge", true,
106 "Only warning if all of peak outer radius is not on detector (default).\n"
107 "If false, do not integrate if the outer radius is not on a detector.");
108
109 declareProperty("AdaptiveQBackground", false,
110 "Default is false. If true, "
111 "BackgroundOuterRadius + AdaptiveQMultiplier * **|Q|** and "
112 "BackgroundInnerRadius + AdaptiveQMultiplier * **|Q|**");
113
114 declareProperty("Ellipsoid", false, "Default is sphere.");
115
116 declareProperty("FixQAxis", false, "Fix one axis of ellipsoid to be along direction of Q.");
117
118 declareProperty("Cylinder", false, "Default is sphere. Use next five parameters for cylinder.");
119
120 declareProperty(std::make_unique<PropertyWithValue<double>>("CylinderLength", 0.0, Direction::Input),
121 "Length of cylinder in which to integrate (in the same units as the "
122 "workspace).");
123
124 declareProperty(std::make_unique<PropertyWithValue<double>>("PercentBackground", 0.0, Direction::Input),
125 "Percent of CylinderLength that is background (20 is 20%)");
126
127 std::vector<std::string> peakNames = FunctionFactory::Instance().getFunctionNames<IPeakFunction>();
128 peakNames.emplace_back("NoFit");
129 declareProperty("ProfileFunction", "Gaussian", std::make_shared<StringListValidator>(peakNames),
130 "Fitting function for profile that is used only with "
131 "Cylinder integration.");
132
133 std::vector<std::string> integrationOptions(2);
134 integrationOptions[0] = "Sum";
135 integrationOptions[1] = "GaussianQuadrature";
136 auto integrationvalidator = std::make_shared<StringListValidator>(integrationOptions);
137 declareProperty("IntegrationOption", "GaussianQuadrature", integrationvalidator,
138 "Integration method for calculating intensity "
139 "used only with Cylinder integration.");
140
141 declareProperty(std::make_unique<FileProperty>("ProfilesFile", "", FileProperty::OptionalSave,
142 std::vector<std::string>(1, "profiles")),
143 "Save (Optionally) as Isaw peaks file with profiles included");
144
145 declareProperty("AdaptiveQMultiplier", 0.0,
146 "PeakRadius + AdaptiveQMultiplier * **|Q|** "
147 "so each peak has a "
148 "different integration radius. Q includes the 2*pi factor.");
149
150 declareProperty("CorrectIfOnEdge", false,
151 "Only warning if all of peak outer radius is not on detector (default).\n"
152 "If false, correct for volume off edge for both background and "
153 "intensity (the peak is assumed uniform Gaussian so this only applies "
154 "to spherical integration).");
155
156 declareProperty("UseOnePercentBackgroundCorrection", true,
157 "If this options is enabled, then the the top 1% of the "
158 "background will be removed"
159 "before the background subtraction.");
160
161 // continued ellipsoid args
162 declareProperty("FixMajorAxisLength", true,
163 "This option is ignored if all peak radii are specified. "
164 "Otherwise, if True the ellipsoid radidi (proportional to "
165 "the sqrt of the eigenvalues of the covariance matrix) are "
166 "scaled such that the major axis radius is equal to the "
167 "PeakRadius. If False then the ellipsoid radii are set to "
168 "3 times the sqrt of the eigenvalues of the covariance "
169 "matrix");
170
171 declareProperty("UseCentroid", false,
172 "Perform integration on estimated centroid not peak position "
173 "(ignored if all three peak radii are specified).");
174
175 auto maxIterValidator = std::make_shared<BoundedValidator<int>>();
176 maxIterValidator->setLower(1);
177 declareProperty("MaxIterations", 1, maxIterValidator,
178 "Number of iterations in covariance estimation (ignored if all "
179 "peak radii are specified). 2-3 should be sufficient.");
180
182 "MaskEdgeTubes", true,
183 "Mask tubes on the edge of all banks in the PeaksWorkspace instrument (note the edge pixels at top/bottom of all "
184 "tubes will always be masked even if this property is False). Note the algorithm will treat "
185 "any masked pixels as edges (including pixels already masked prior to the execution of this algorithm) - this "
186 "means a custom mask can be applied to the PeaksWorkspace before integration.");
187
188 // Group Properties
189 std::string general_grp = "General Inputs";
190 std::string cylin_grp = "Cylindrical Integration";
191 std::string ellip_grp = "Ellipsoid Integration";
192
193 setPropertyGroup("InputWorkspace", general_grp);
194 setPropertyGroup("PeakRadius", general_grp);
195 setPropertyGroup("BackgroundInnerRadius", general_grp);
196 setPropertyGroup("BackgroundOuterRadius", general_grp);
197 setPropertyGroup("PeaksWorkspace", general_grp);
198 setPropertyGroup("OutputWorkspace", general_grp);
199 setPropertyGroup("ReplaceIntensity", general_grp);
200 setPropertyGroup("IntegrateIfOnEdge", general_grp);
201 setPropertyGroup("AdaptiveQBackground", general_grp);
202
203 setPropertyGroup("Ellipsoid", ellip_grp);
204 setPropertyGroup("FixQAxis", ellip_grp);
205
206 setPropertyGroup("Cylinder", cylin_grp);
207 setPropertyGroup("CylinderLength", cylin_grp);
208 setPropertyGroup("PercentBackground", cylin_grp);
209 setPropertyGroup("ProfileFunction", cylin_grp);
210 setPropertyGroup("IntegrationOption", cylin_grp);
211 setPropertyGroup("ProfilesFile", cylin_grp);
212
213 setPropertyGroup("AdaptiveQMultiplier", general_grp);
214 setPropertyGroup("CorrectIfOnEdge", general_grp);
215 setPropertyGroup("UseOnePercentBackgroundCorrection", general_grp);
216
217 setPropertyGroup("FixMajorAxisLength", ellip_grp);
218 setPropertyGroup("UseCentroid", ellip_grp);
219 setPropertyGroup("MaxIterations", ellip_grp);
220
221 setPropertyGroup("MaskEdgeTubes", general_grp);
222
223 // SetValue when another property value changes
224 setPropertySettings("Ellipsoid", std::make_unique<SetValueWhenProperty>(
225 "Cylinder", [](std::string ellipsoid, const std::string &cylinder) {
226 // Set Ellipsoid to 0, if user has set Cylinder to 1
227 if (ellipsoid == "1" && cylinder == "1") {
228 return std::string{"0"};
229 } else {
230 return ellipsoid;
231 };
232 }));
233 setPropertySettings("Cylinder", std::make_unique<SetValueWhenProperty>(
234 "Ellipsoid", [](std::string cylinder, const std::string &ellipsoid) {
235 // Set Cylinder to 0, if user has set Ellipsoid to 1
236 if (cylinder == "1" && ellipsoid == "1") {
237 return std::string{"0"};
238 } else {
239 return cylinder;
240 };
241 }));
242
243 // Set these Properties as visible only when Cylinder = 1
244 setPropertySettings("CylinderLength", std::make_unique<VisibleWhenProperty>("Cylinder", IS_EQUAL_TO, "1"));
245 setPropertySettings("PercentBackground", std::make_unique<VisibleWhenProperty>("Cylinder", IS_EQUAL_TO, "1"));
246 setPropertySettings("ProfileFunction", std::make_unique<VisibleWhenProperty>("Cylinder", IS_EQUAL_TO, "1"));
247 setPropertySettings("IntegrationOption", std::make_unique<VisibleWhenProperty>("Cylinder", IS_EQUAL_TO, "1"));
248 setPropertySettings("ProfilesFile", std::make_unique<VisibleWhenProperty>("Cylinder", IS_EQUAL_TO, "1"));
249
250 // Set these Properties as visible only when Ellipsoid = 1
251 setPropertySettings("FixQAxis", std::make_unique<VisibleWhenProperty>("Ellipsoid", IS_EQUAL_TO, "1"));
252 setPropertySettings("FixMajorAxisLength", std::make_unique<VisibleWhenProperty>("Ellipsoid", IS_EQUAL_TO, "1"));
253 setPropertySettings("UseCentroid", std::make_unique<VisibleWhenProperty>("Ellipsoid", IS_EQUAL_TO, "1"));
254 setPropertySettings("MaxIterations", std::make_unique<VisibleWhenProperty>("Ellipsoid", IS_EQUAL_TO, "1"));
255
256 // Disable / greyed out these Properties based on the value of another
257 setPropertySettings("CorrectIfOnEdge",
258 std::make_unique<Kernel::EnabledWhenProperty>("IntegrateIfOnEdge", IS_EQUAL_TO, "1"));
259}
260
261std::map<std::string, std::string> IntegratePeaksMD2::validateInputs() {
262 std::map<std::string, std::string> result;
263
264 std::vector<double> PeakRadius = getProperty("PeakRadius");
265 std::vector<double> BackgroundInnerRadius = getProperty("BackgroundInnerRadius");
266 std::vector<double> BackgroundOuterRadius = getProperty("BackgroundOuterRadius");
267 bool ellipsoid = getProperty("Ellipsoid");
268 bool cylinder = getProperty("Cylinder");
269
270 if (PeakRadius.size() != 1 && PeakRadius.size() != 3) {
271 std::stringstream errmsg;
272 errmsg << "Only one or three values should be specified";
273 result["PeakRadius"] = errmsg.str();
274 }
275
276 if (!ellipsoid && PeakRadius.size() != 1) {
277 std::stringstream errmsg;
278 errmsg << "One value must be specified when Ellipsoid is false";
279 result["PeakRadius"] = errmsg.str();
280 }
281
282 if (BackgroundInnerRadius.size() != 1 && BackgroundInnerRadius.size() != 3) {
283 std::stringstream errmsg;
284 errmsg << "Only one or three values should be specified";
285 result["BackgroundInnerRadius"] = errmsg.str();
286 }
287
288 if (!ellipsoid && BackgroundInnerRadius.size() != 1) {
289 std::stringstream errmsg;
290 errmsg << "One value must be specified when Ellipsoid is false";
291 result["BackgroundInnerRadius"] = errmsg.str();
292 }
293
294 if (BackgroundOuterRadius.size() != 1 && BackgroundOuterRadius.size() != 3) {
295 std::stringstream errmsg;
296 errmsg << "Only one or three values should be specified";
297 result["BackgroundOuterRadius"] = errmsg.str();
298 }
299
300 if (!ellipsoid && BackgroundOuterRadius.size() != 1) {
301 std::stringstream errmsg;
302 errmsg << "One value must be specified when Ellipsoid is false";
303 result["BackgroundOuterRadius"] = errmsg.str();
304 }
305
306 if (ellipsoid && cylinder) {
307 std::stringstream errmsg;
308 errmsg << "Ellipsoid and Cylinder cannot both be true";
309 result["Ellipsoid"] = errmsg.str();
310 result["Cylinder"] = errmsg.str();
311 }
312
313 return result;
314}
315
316//----------------------------------------------------------------------------------------------
321template <typename MDE, size_t nd> void IntegratePeaksMD2::integrate(typename MDEventWorkspace<MDE, nd>::sptr ws) {
322 if (nd != 3)
323 throw std::invalid_argument("For now, we expect the input MDEventWorkspace "
324 "to have 3 dimensions only.");
325
327 IPeaksWorkspace_sptr inPeakWS = getProperty("PeaksWorkspace");
328
330 IPeaksWorkspace_sptr peakWS = getProperty("OutputWorkspace");
331 if (peakWS != inPeakWS)
332 peakWS = inPeakWS->clone();
333 // This only fails in the unit tests which say that MaskBTP is not registered
334 bool maskTubes = getProperty("MaskEdgeTubes");
335 try {
336 PeaksWorkspace_sptr p = std::dynamic_pointer_cast<PeaksWorkspace>(inPeakWS);
337 if (p) {
338 if (maskTubes) {
339 runMaskDetectors(p, "Tube", "edges");
340 }
341 runMaskDetectors(p, "Pixel", "edges");
342 }
343 } catch (...) {
344 g_log.error("Can't execute MaskBTP algorithm for this instrument to set "
345 "edge for IntegrateIfOnEdge option");
346 }
347
348 calculateE1(inPeakWS->detectorInfo()); // fill E1Vec for use in detectorQ
350
352 std::vector<double> PeakRadius = getProperty("PeakRadius");
354 std::vector<double> BackgroundOuterRadius = getProperty("BackgroundOuterRadius");
356 std::vector<double> BackgroundInnerRadius = getProperty("BackgroundInnerRadius");
358 bool useOnePercentBackgroundCorrection = getProperty("UseOnePercentBackgroundCorrection");
359
360 bool manualEllip = false;
361 if (PeakRadius.size() > 1) {
362 manualEllip = true;
363 // make sure the background radii are 3 values (they default to 1)
364 if (BackgroundInnerRadius.size() == 1)
365 BackgroundInnerRadius.resize(3, BackgroundInnerRadius[0]);
366 if (BackgroundOuterRadius.size() == 1)
367 BackgroundOuterRadius.resize(3, BackgroundOuterRadius[0]);
368 }
369
370 double minInnerRadius = PeakRadius[0];
371 for (size_t r = 0; r < BackgroundInnerRadius.size(); r++) {
372 if (manualEllip) {
373 minInnerRadius = PeakRadius[r];
374 }
375 if (BackgroundInnerRadius[r] < minInnerRadius)
376 BackgroundInnerRadius[r] = minInnerRadius;
377 }
378 // Ellipsoid
379 bool isEllipse = getProperty("Ellipsoid");
380 bool qAxisIsFixed = getProperty("FixQAxis");
381 bool majorAxisLengthFixed = getProperty("FixMajorAxisLength");
382 bool useCentroid = getProperty("UseCentroid");
383 int maxCovarIter = getProperty("MaxIterations");
385 double cylinderLength = getProperty("CylinderLength");
386 Workspace2D_sptr wsProfile2D, wsFit2D, wsDiff2D;
387 size_t numSteps = 0;
388 bool cylinderBool = getProperty("Cylinder");
389 bool adaptiveQBackground = getProperty("AdaptiveQBackground");
390 double adaptiveQMultiplier = getProperty("AdaptiveQMultiplier");
391 double adaptiveQBackgroundMultiplier = 0.0;
392 if (adaptiveQBackground)
393 adaptiveQBackgroundMultiplier = adaptiveQMultiplier;
394 std::vector<double> PeakRadiusVector(peakWS->getNumberPeaks(), PeakRadius[0]);
395 std::vector<double> BackgroundInnerRadiusVector(peakWS->getNumberPeaks(), BackgroundInnerRadius[0]);
396 std::vector<double> BackgroundOuterRadiusVector(peakWS->getNumberPeaks(), BackgroundOuterRadius[0]);
397 if (cylinderBool) {
398 numSteps = 100;
399 size_t histogramNumber = peakWS->getNumberPeaks();
400 Workspace_sptr wsProfile = WorkspaceFactory::Instance().create("Workspace2D", histogramNumber, numSteps, numSteps);
401 wsProfile2D = std::dynamic_pointer_cast<Workspace2D>(wsProfile);
402 AnalysisDataService::Instance().addOrReplace("ProfilesData", wsProfile2D);
403 Workspace_sptr wsFit = WorkspaceFactory::Instance().create("Workspace2D", histogramNumber, numSteps, numSteps);
404 wsFit2D = std::dynamic_pointer_cast<Workspace2D>(wsFit);
405 AnalysisDataService::Instance().addOrReplace("ProfilesFit", wsFit2D);
406 Workspace_sptr wsDiff = WorkspaceFactory::Instance().create("Workspace2D", histogramNumber, numSteps, numSteps);
407 wsDiff2D = std::dynamic_pointer_cast<Workspace2D>(wsDiff);
408 AnalysisDataService::Instance().addOrReplace("ProfilesFitDiff", wsDiff2D);
409 auto newAxis1 = std::make_unique<TextAxis>(peakWS->getNumberPeaks());
410 auto newAxis2 = std::make_unique<TextAxis>(peakWS->getNumberPeaks());
411 auto newAxis3 = std::make_unique<TextAxis>(peakWS->getNumberPeaks());
412 auto newAxis1Raw = newAxis1.get();
413 auto newAxis2Raw = newAxis2.get();
414 auto newAxis3Raw = newAxis3.get();
415 wsProfile2D->replaceAxis(1, std::move(newAxis1));
416 wsFit2D->replaceAxis(1, std::move(newAxis2));
417 wsDiff2D->replaceAxis(1, std::move(newAxis3));
418 for (int i = 0; i < peakWS->getNumberPeaks(); ++i) {
419 // Get a direct ref to that peak.
420 IPeak &p = peakWS->getPeak(i);
421 std::ostringstream label;
422 label << Utils::round(p.getH()) << "_" << Utils::round(p.getK()) << "_" << Utils::round(p.getL()) << "_"
423 << p.getRunNumber();
424 newAxis1Raw->setLabel(i, label.str());
425 newAxis2Raw->setLabel(i, label.str());
426 newAxis3Raw->setLabel(i, label.str());
427 }
428 }
429 double percentBackground = getProperty("PercentBackground");
430 size_t peakMin = 0;
431 size_t peakMax = numSteps;
432 double ratio = 0.0;
433 if (cylinderBool) {
434 peakMin = static_cast<size_t>(static_cast<double>(numSteps) * percentBackground / 100.);
435 peakMax = numSteps - peakMin - 1;
436 size_t numPeakCh = peakMax - peakMin + 1; // number of peak channels
437 size_t numBkgCh = numSteps - numPeakCh; // number of background channels
438 ratio = static_cast<double>(numPeakCh) / static_cast<double>(numBkgCh);
439 }
441 bool replaceIntensity = getProperty("ReplaceIntensity");
442 bool integrateEdge = getProperty("IntegrateIfOnEdge");
443 bool correctEdge = getProperty("CorrectIfOnEdge");
444
445 std::string profileFunction = getProperty("ProfileFunction");
446 std::string integrationOption = getProperty("IntegrationOption");
447 std::ofstream out;
448 if (cylinderBool && profileFunction != "NoFit") {
449 std::string outFile = getProperty("InputWorkspace");
450 outFile.append(profileFunction);
451 outFile.append(".dat");
452 std::string save_path = ConfigService::Instance().getString("defaultsave.directory");
453 outFile = save_path + outFile;
454 out.open(outFile.c_str(), std::ofstream::out);
455 }
456 // volume of Background sphere with inner volume subtracted
457 double volumeBkg = 4.0 / 3.0 * M_PI * (std::pow(BackgroundOuterRadius[0], 3) - std::pow(BackgroundOuterRadius[0], 3));
458 // volume of PeakRadius sphere
459 double volumeRadius = 4.0 / 3.0 * M_PI * std::pow(PeakRadius[0], 3);
460
461 // Initialize progress reporting
462 int nPeaks = peakWS->getNumberPeaks();
463 Progress progress(this, 0., 1., nPeaks);
464 bool doParallel = cylinderBool ? false : Kernel::threadSafe(*ws, *peakWS);
465 PARALLEL_FOR_IF(doParallel)
466 for (int i = 0; i < nPeaks; ++i) {
468 progress.report();
469
470 // Get a direct ref to that peak.
471 IPeak &p = peakWS->getPeak(i);
472
473 // Get the peak center as a position in the dimensions of the workspace
474 V3D pos;
475 if (CoordinatesToUse == Mantid::Kernel::QLab) //"Q (lab frame)"
476 pos = p.getQLabFrame();
477 else if (CoordinatesToUse == Mantid::Kernel::QSample) //"Q (sample frame)"
478 pos = p.getQSampleFrame();
479 else if (CoordinatesToUse == Mantid::Kernel::HKL) //"HKL"
480 pos = p.getHKL();
481
482 // Do not integrate if sphere is off edge of detector
483
484 const double edgeDist = calculateDistanceToEdge(p.getQLabFrame());
485 if (edgeDist < std::max(BackgroundOuterRadius[0], PeakRadius[0])) {
486 g_log.warning() << "Warning: sphere/cylinder for integration is off edge "
487 "of detector for peak "
488 << i << "; radius of edge = " << edgeDist << '\n';
489 if (!integrateEdge) {
490 if (replaceIntensity) {
491 p.setIntensity(0.0);
492 p.setSigmaIntensity(0.0);
493 }
494 continue;
495 }
496 }
497
498 // Build the sphere transformation
499 bool dimensionsUsed[nd];
500 coord_t center[nd];
501 for (size_t d = 0; d < nd; ++d) {
502 dimensionsUsed[d] = true; // Use all dimensions
503 center[d] = static_cast<coord_t>(pos[d]);
504 }
505 signal_t signal = 0;
506 signal_t errorSquared = 0;
507 signal_t bgSignal = 0;
508 signal_t bgErrorSquared = 0;
509 double background_total = 0.0;
510 if (!cylinderBool) {
511 // modulus of Q
512 coord_t lenQpeak = 0.0;
513 if (adaptiveQMultiplier != 0.0) {
514 lenQpeak = 0.0;
515 for (size_t d = 0; d < nd; ++d) {
516 lenQpeak += center[d] * center[d];
517 }
518 lenQpeak = std::sqrt(lenQpeak);
519 }
520 double adaptiveRadius = adaptiveQMultiplier * lenQpeak + *std::max_element(PeakRadius.begin(), PeakRadius.end());
521 if (adaptiveRadius <= 0.0) {
522 g_log.error() << "Error: Radius for integration sphere of peak " << i << " is negative = " << adaptiveRadius
523 << '\n';
524 adaptiveRadius = 0.;
525 p.setIntensity(0.0);
526 p.setSigmaIntensity(0.0);
527 PeakRadiusVector[i] = 0.0;
528 BackgroundInnerRadiusVector[i] = 0.0;
529 BackgroundOuterRadiusVector[i] = 0.0;
530 continue;
531 }
532 PeakRadiusVector[i] = adaptiveRadius;
533 BackgroundInnerRadiusVector[i] = adaptiveQBackgroundMultiplier * lenQpeak +
534 *std::max_element(BackgroundInnerRadius.begin(), BackgroundInnerRadius.end());
535 BackgroundOuterRadiusVector[i] = adaptiveQBackgroundMultiplier * lenQpeak +
536 *std::max_element(BackgroundOuterRadius.begin(), BackgroundOuterRadius.end());
537 // define the radius squared for a sphere intially
538 CoordTransformDistance getRadiusSq(nd, center, dimensionsUsed);
539 // set spherical shape
540 PeakShape *sphereShape =
541 new PeakShapeSpherical(PeakRadiusVector[i], BackgroundInnerRadiusVector[i], BackgroundOuterRadiusVector[i],
542 CoordinatesToUse, this->name(), this->version());
543 p.setPeakShape(sphereShape);
544 const double scaleFactor = pow(PeakRadiusVector[i], 3) /
545 (pow(BackgroundOuterRadiusVector[i], 3) - pow(BackgroundInnerRadiusVector[i], 3));
546 // Integrate spherical background shell if specified
547 if (BackgroundOuterRadius[0] > PeakRadius[0]) {
548 // Get the total signal inside background shell
549 ws->getBox()->integrateSphere(
550 getRadiusSq, static_cast<coord_t>(pow(BackgroundOuterRadiusVector[i], 2)), bgSignal, bgErrorSquared,
551 static_cast<coord_t>(pow(BackgroundInnerRadiusVector[i], 2)), useOnePercentBackgroundCorrection);
552 // correct bg signal by Vpeak/Vshell (same for sphere and ellipse)
553 bgSignal *= scaleFactor;
554 bgErrorSquared *= scaleFactor * scaleFactor;
555 }
556 // if ellipsoid find covariance and centroid in spherical region
557 // using one-pass algorithm from https://doi.org/10.1145/359146.359153
558 if (isEllipse) {
559 // flat bg to subtract
560 const auto bgDensity = bgSignal / (4 * M_PI * pow(PeakRadiusVector[i], 3) / 3);
561 std::vector<V3D> eigenvects;
562 std::vector<double> eigenvals;
563 V3D translation(0.0, 0.0, 0.0); // translation from peak pos to centroid
564 if (PeakRadius.size() == 1) {
565 V3D mean(0.0, 0.0, 0.0); // vector to hold centroid
566 findEllipsoid<MDE, nd>(ws, getRadiusSq, pos, static_cast<coord_t>(pow(PeakRadiusVector[i], 2)), qAxisIsFixed,
567 useCentroid, bgDensity, eigenvects, eigenvals, mean, maxCovarIter);
568 if (!majorAxisLengthFixed) {
569 // replace radius for this peak with 3*stdev along major axis
570 auto max_stdev = sqrt(*std::max_element(eigenvals.begin(), eigenvals.end()));
571 BackgroundOuterRadiusVector[i] = 3 * max_stdev * (BackgroundOuterRadiusVector[i] / PeakRadiusVector[i]);
572 BackgroundInnerRadiusVector[i] = 3 * max_stdev * (BackgroundInnerRadiusVector[i] / PeakRadiusVector[i]);
573 PeakRadiusVector[i] = 3 * max_stdev;
574 }
575 if (useCentroid) {
576 // calculate translation to apply when drawing
577 translation = mean - pos;
578 // update integration center with mean
579 for (size_t d = 0; d < 3; ++d) {
580 center[d] = static_cast<coord_t>(mean[d]);
581 }
582 }
583 } else {
584 // Use the manually specified radii instead of finding them via
585 // findEllipsoid
586 std::transform(PeakRadius.begin(), PeakRadius.end(), std::back_inserter(eigenvals),
587 [](double r) { return std::pow(r, 2.0); });
588 eigenvects.push_back(V3D(1.0, 0.0, 0.0));
589 eigenvects.push_back(V3D(0.0, 1.0, 0.0));
590 eigenvects.push_back(V3D(0.0, 0.0, 1.0));
591 }
592 // transform ellispoid onto sphere of radius = R
593 getRadiusSq = CoordTransformDistance(nd, center, dimensionsUsed, 1, /* outD */
594 eigenvects, eigenvals);
595 // Integrate ellipsoid background shell if specified
596 if (PeakRadius.size() == 1) {
597 if (BackgroundOuterRadius[0] > PeakRadius[0]) {
598 // Get the total signal inside "BackgroundOuterRadius"
599 bgSignal = 0;
600 bgErrorSquared = 0;
601 ws->getBox()->integrateSphere(
602 getRadiusSq, static_cast<coord_t>(pow(BackgroundOuterRadiusVector[i], 2)), bgSignal, bgErrorSquared,
603 static_cast<coord_t>(pow(BackgroundInnerRadiusVector[i], 2)), useOnePercentBackgroundCorrection);
604 // correct bg signal by Vpeak/Vshell (same as previously
605 // calculated for sphere)
606 bgSignal *= scaleFactor;
607 bgErrorSquared *= scaleFactor * scaleFactor;
608 }
609 // set peak shape
610 // get radii in same proprtion as eigenvalues
611 auto max_stdev = pow(*std::max_element(eigenvals.begin(), eigenvals.end()), 0.5);
612 std::vector<double> peakRadii(3, 0.0);
613 std::vector<double> backgroundInnerRadii(3, 0.0);
614 std::vector<double> backgroundOuterRadii(3, 0.0);
615 for (size_t irad = 0; irad < peakRadii.size(); irad++) {
616 auto scale = pow(eigenvals[irad], 0.5) / max_stdev;
617 peakRadii[irad] = PeakRadiusVector[i] * scale;
618 backgroundInnerRadii[irad] = BackgroundInnerRadiusVector[i] * scale;
619 backgroundOuterRadii[irad] = BackgroundOuterRadiusVector[i] * scale;
620 }
621 PeakShape *ellipsoidShape =
622 new PeakShapeEllipsoid(eigenvects, peakRadii, backgroundInnerRadii, backgroundOuterRadii,
623 CoordinatesToUse, this->name(), this->version(), translation);
624 p.setPeakShape(ellipsoidShape);
625 } else {
626 // Use the manually specified radii instead of finding them via
627 // findEllipsoid
628 std::vector<double> eigenvals_background_inner;
629 std::vector<double> eigenvals_background_outer;
630 std::transform(BackgroundInnerRadius.begin(), BackgroundInnerRadius.end(),
631 std::back_inserter(eigenvals_background_inner), [](double r) { return std::pow(r, 2.0); });
632 std::transform(BackgroundOuterRadius.begin(), BackgroundOuterRadius.end(),
633 std::back_inserter(eigenvals_background_outer), [](double r) { return std::pow(r, 2.0); });
634
635 if (BackgroundOuterRadiusVector[0] > PeakRadiusVector[0]) {
636 // transform ellispoid onto sphere of radius = R
637 auto getRadiusSqInner = CoordTransformDistance(nd, center, dimensionsUsed, 1, /* outD */
638 eigenvects, eigenvals_background_inner);
639 auto getRadiusSqOuter = CoordTransformDistance(nd, center, dimensionsUsed, 1, /* outD */
640 eigenvects, eigenvals_background_outer);
641 // Get the total signal inside "BackgroundOuterRadius"
642 bgSignal = 0;
643 bgErrorSquared = 0;
644 signal_t bgSignalInner = 0;
645 signal_t bgSignalOuter = 0;
646 signal_t bgErrorSquaredInner = 0;
647 signal_t bgErrorSquaredOuter = 0;
648 ws->getBox()->integrateSphere(getRadiusSqInner,
649 static_cast<coord_t>(pow(BackgroundInnerRadiusVector[i], 2)), bgSignalInner,
650 bgErrorSquaredInner, 0.0, useOnePercentBackgroundCorrection);
651 ws->getBox()->integrateSphere(getRadiusSqOuter,
652 static_cast<coord_t>(pow(BackgroundOuterRadiusVector[i], 2)), bgSignalOuter,
653 bgErrorSquaredOuter, 0.0, useOnePercentBackgroundCorrection);
654 // correct bg signal by Vpeak/Vshell (same as previously
655 // calculated for sphere)
656 bgSignal = bgSignalOuter - bgSignalInner;
657 bgErrorSquared = bgErrorSquaredInner + bgErrorSquaredOuter;
658 g_log.debug() << "unscaled background signal from ellipsoid integration = " << bgSignal << '\n';
659 const double scale = (PeakRadius[0] * PeakRadius[1] * PeakRadius[2]) /
660 (BackgroundOuterRadius[0] * BackgroundOuterRadius[1] * BackgroundOuterRadius[2] -
661 BackgroundInnerRadius[0] * BackgroundInnerRadius[1] * BackgroundInnerRadius[2]);
662 bgSignal *= scale;
663 bgErrorSquared *= pow(scale, 2);
664 }
665 // set peak shape
666 // get radii in same proprtion as eigenvalues
667 auto max_stdev = pow(*std::max_element(eigenvals.begin(), eigenvals.end()), 0.5);
668 auto max_stdev_inner =
669 pow(*std::max_element(eigenvals_background_inner.begin(), eigenvals_background_inner.end()), 0.5);
670 auto max_stdev_outer =
671 pow(*std::max_element(eigenvals_background_outer.begin(), eigenvals_background_outer.end()), 0.5);
672 std::vector<double> peakRadii(3, 0.0);
673 std::vector<double> backgroundInnerRadii(3, 0.0);
674 std::vector<double> backgroundOuterRadii(3, 0.0);
675 for (size_t irad = 0; irad < peakRadii.size(); irad++) {
676 peakRadii[irad] = PeakRadiusVector[i] * pow(eigenvals[irad], 0.5) / max_stdev;
677 backgroundInnerRadii[irad] =
678 BackgroundInnerRadiusVector[i] * pow(eigenvals_background_inner[irad], 0.5) / max_stdev_inner;
679 backgroundOuterRadii[irad] =
680 BackgroundOuterRadiusVector[i] * pow(eigenvals_background_outer[irad], 0.5) / max_stdev_outer;
681 }
682 PeakShape *ellipsoidShape =
683 new PeakShapeEllipsoid(eigenvects, peakRadii, backgroundInnerRadii, backgroundOuterRadii,
684 CoordinatesToUse, this->name(), this->version());
685 p.setPeakShape(ellipsoidShape);
686 }
687 }
688 ws->getBox()->integrateSphere(getRadiusSq, static_cast<coord_t>(PeakRadiusVector[i] * PeakRadiusVector[i]),
689 signal, errorSquared, 0.0 /* innerRadiusSquared */,
690 useOnePercentBackgroundCorrection);
691 //
692 } else {
693 CoordTransformDistance cylinder(nd, center, dimensionsUsed, 2);
694
695 // Perform the integration into whatever box is contained within.
696 Counts signal_fit(numSteps);
697 signal_fit = 0;
698
699 ws->getBox()->integrateCylinder(cylinder, static_cast<coord_t>(PeakRadius[0]),
700 static_cast<coord_t>(cylinderLength), signal, errorSquared,
701 signal_fit.mutableRawData());
702
703 // Integrate around the background radius
704 if (BackgroundOuterRadius[0] > PeakRadius[0]) {
705 // Get the total signal inside "BackgroundOuterRadius"
706 signal_fit = 0;
707
708 ws->getBox()->integrateCylinder(cylinder, static_cast<coord_t>(BackgroundOuterRadius[0]),
709 static_cast<coord_t>(cylinderLength), bgSignal, bgErrorSquared,
710 signal_fit.mutableRawData());
711
712 Points points(signal_fit.size(), LinearGenerator(0, 1));
713 wsProfile2D->setHistogram(i, points, signal_fit);
714
715 // Evaluate the signal inside "BackgroundInnerRadius"
716 signal_t interiorSignal = 0;
717 signal_t interiorErrorSquared = 0;
718
719 // Integrate this 3rd radius, if needed
720 if (BackgroundInnerRadius[0] != PeakRadius[0]) {
721 ws->getBox()->integrateCylinder(cylinder, static_cast<coord_t>(BackgroundInnerRadius[0]),
722 static_cast<coord_t>(cylinderLength), interiorSignal, interiorErrorSquared,
723 signal_fit.mutableRawData());
724 } else {
725 // PeakRadius == BackgroundInnerRadius, so use the previous
726 // value
727 interiorSignal = signal;
728 interiorErrorSquared = errorSquared;
729 }
730 // Subtract the peak part to get the intensity in the shell
731 // (BackgroundInnerRadius < r < BackgroundOuterRadius)
732 bgSignal -= interiorSignal;
733 // We can subtract the error (instead of adding) because the two
734 // values are 100% dependent; this is the same as integrating a
735 // shell.
736 bgErrorSquared -= interiorErrorSquared;
737 // Relative volume of peak vs the BackgroundOuterRadius cylinder
738 const double radiusRatio = (PeakRadius[0] / BackgroundOuterRadius[0]);
739 const double peakVolume = radiusRatio * radiusRatio * cylinderLength;
740
741 // Relative volume of the interior of the shell vs overall
742 // background
743 const double interiorRatio = (BackgroundInnerRadius[0] / BackgroundOuterRadius[0]);
744 // Volume of the bg shell, relative to the volume of the
745 // BackgroundOuterRadius cylinder
746 const double bgVolume = 1.0 - interiorRatio * interiorRatio * cylinderLength;
747
748 // Finally, you will multiply the bg intensity by this to get the
749 // estimated background under the peak volume
750 const double scaleFactor = peakVolume / bgVolume;
751 bgSignal *= scaleFactor;
752 bgErrorSquared *= scaleFactor * scaleFactor;
753 } else {
754 Points points(signal_fit.size(), LinearGenerator(0, 1));
755 wsProfile2D->setHistogram(i, points, signal_fit);
756 }
757
758 if (profileFunction == "NoFit") {
759 signal = 0.;
760 for (size_t j = 0; j < numSteps; j++) {
761 if (j < peakMin || j > peakMax)
762 background_total = background_total + wsProfile2D->y(i)[j];
763 else
764 signal = signal + wsProfile2D->y(i)[j];
765 }
766 errorSquared = std::fabs(signal);
767 } else {
768
769 auto fitAlgorithm = createChildAlgorithm("Fit", -1, -1, false);
770 // fitAlgorithm->setProperty("CreateOutput", true);
771 // fitAlgorithm->setProperty("Output", "FitPeaks1D");
772 std::string myFunc = std::string("name=LinearBackground;name=") + profileFunction;
773 auto maxPeak = std::max_element(signal_fit.begin(), signal_fit.end());
774
775 std::ostringstream strs;
776 strs << maxPeak[0];
777 std::string strMax = strs.str();
778 if (profileFunction == "Gaussian") {
779 myFunc += ", PeakCentre=50, Height=" + strMax;
780 fitAlgorithm->setProperty("Constraints", "40<f1.PeakCentre<60");
781 } else if (profileFunction == "BackToBackExponential" || profileFunction == "IkedaCarpenterPV") {
782 myFunc += ", X0=50, I=" + strMax;
783 fitAlgorithm->setProperty("Constraints", "40<f1.X0<60");
784 }
785 fitAlgorithm->setProperty("CalcErrors", true);
786 fitAlgorithm->setProperty("Function", myFunc);
787 fitAlgorithm->setProperty("InputWorkspace", wsProfile2D);
788 fitAlgorithm->setProperty("WorkspaceIndex", static_cast<int>(i));
789 try {
790 fitAlgorithm->executeAsChildAlg();
791 } catch (...) {
792 g_log.error("Can't execute Fit algorithm");
793 continue;
794 }
795
796 IFunction_sptr ifun = fitAlgorithm->getProperty("Function");
797 if (i == 0) {
798 out << std::setw(20) << "spectrum"
799 << " ";
800 for (size_t j = 0; j < ifun->nParams(); ++j)
801 out << std::setw(20) << ifun->parameterName(j) << " ";
802 out << std::setw(20) << "chi2"
803 << " ";
804 out << "\n";
805 }
806 out << std::setw(20) << i << " ";
807 for (size_t j = 0; j < ifun->nParams(); ++j) {
808 out << std::setw(20) << std::fixed << std::setprecision(10) << ifun->getParameter(j) << " ";
809 }
810 double chi2 = fitAlgorithm->getProperty("OutputChi2overDoF");
811 out << std::setw(20) << std::fixed << std::setprecision(10) << chi2 << "\n";
812
813 std::shared_ptr<const CompositeFunction> fun = std::dynamic_pointer_cast<const CompositeFunction>(ifun);
814
815 const auto &x = wsProfile2D->x(i);
816 wsFit2D->setSharedX(i, wsProfile2D->sharedX(i));
817 wsDiff2D->setSharedX(i, wsProfile2D->sharedX(i));
818
819 FunctionDomain1DVector domain(x.rawData());
820 FunctionValues yy(domain);
821 fun->function(domain, yy);
822 auto funcValues = yy.toVector();
823
824 wsFit2D->mutableY(i) = funcValues;
825 wsDiff2D->setSharedY(i, wsProfile2D->sharedY(i));
826 wsDiff2D->mutableY(i) -= wsFit2D->y(i);
827
828 // Calculate intensity
829 signal = 0.0;
830 if (integrationOption == "Sum") {
831 for (size_t j = peakMin; j <= peakMax; j++)
832 if (std::isfinite(yy[j]))
833 signal += yy[j];
834 } else {
835 gsl_integration_workspace *w = gsl_integration_workspace_alloc(1000);
836
837 double error;
838
839 gsl_function F;
840 F.function = &Mantid::MDAlgorithms::f_eval2;
841 F.params = &fun;
842
843 gsl_integration_qags(&F, x[peakMin], x[peakMax], 0, 1e-7, 1000, w, &signal, &error);
844
845 gsl_integration_workspace_free(w);
846 }
847 errorSquared = std::fabs(signal);
848 // Get background counts
849 for (size_t j = 0; j < numSteps; j++) {
850 double background = ifun->getParameter(0) + ifun->getParameter(1) * x[j];
851 if (j < peakMin || j > peakMax)
852 background_total = background_total + background;
853 }
854 }
855 }
856 checkOverlap(i, peakWS, CoordinatesToUse, 2.0 * std::max(PeakRadiusVector[i], BackgroundOuterRadiusVector[i]));
857 // Save it back in the peak object.
858 if (signal != 0. || replaceIntensity) {
859 double edgeMultiplier = 1.0;
860 double peakMultiplier = 1.0;
861 if (correctEdge) {
862 if (edgeDist < BackgroundOuterRadius[0]) {
863 double e1 = BackgroundOuterRadius[0] - edgeDist;
864 // volume of cap of sphere with h = edge
865 double f1 = M_PI * std::pow(e1, 2) / 3 * (3 * BackgroundOuterRadius[0] - e1);
866 edgeMultiplier = volumeBkg / (volumeBkg - f1);
867 }
868 if (edgeDist < PeakRadius[0]) {
869 double sigma = PeakRadius[0] / 3.0;
870 // assume gaussian peak
871 double e1 = std::exp(-std::pow(edgeDist, 2) / (2 * sigma * sigma)) * PeakRadius[0];
872 // volume of cap of sphere with h = edge
873 double f1 = M_PI * std::pow(e1, 2) / 3 * (3 * PeakRadius[0] - e1);
874 peakMultiplier = volumeRadius / (volumeRadius - f1);
875 }
876 }
877
878 p.setIntensity(peakMultiplier * signal - edgeMultiplier * (ratio * background_total + bgSignal));
879 p.setSigmaIntensity(sqrt(peakMultiplier * errorSquared +
880 edgeMultiplier * (ratio * ratio * std::fabs(background_total) + bgErrorSquared)));
881 }
882
883 g_log.information() << "Peak " << i << " at " << pos << ": signal " << signal << " (sig^2 " << errorSquared
884 << "), with background " << bgSignal + ratio * background_total << " (sig^2 "
885 << bgErrorSquared + ratio * ratio * std::fabs(background_total) << ") subtracted.\n";
887 }
889 // This flag is used by the PeaksWorkspace to evaluate whether it has
890 // been integrated.
891 peakWS->mutableRun().addProperty("PeaksIntegrated", 1, true);
892 // These flags are specific to the algorithm.
893 peakWS->mutableRun().addProperty("PeakRadius", PeakRadiusVector, true);
894 peakWS->mutableRun().addProperty("BackgroundInnerRadius", BackgroundInnerRadiusVector, true);
895 peakWS->mutableRun().addProperty("BackgroundOuterRadius", BackgroundOuterRadiusVector, true);
896
897 // save profiles in peaks file
898 const std::string outfile = getProperty("ProfilesFile");
899 if (outfile.length() > 0) {
900 IAlgorithm_sptr alg;
901 try {
902 alg = createChildAlgorithm("SaveIsawPeaks", -1, -1, false);
903 } catch (Exception::NotFoundError &) {
904 g_log.error("Can't locate SaveIsawPeaks algorithm");
905 throw;
906 }
907 alg->setProperty("InputWorkspace", peakWS);
908 alg->setProperty("ProfileWorkspace", wsProfile2D);
909 alg->setPropertyValue("Filename", outfile);
910 alg->execute();
911 }
912 // Save the output
913 setProperty("OutputWorkspace", peakWS);
914}
915
933template <typename MDE, size_t nd>
935 const CoordTransform &getRadiusSq, const V3D &pos, const coord_t &radiusSquared,
936 const bool &qAxisIsFixed, const bool &useCentroid, const double &bgDensity,
937 std::vector<V3D> &eigenvects, std::vector<double> &eigenvals, V3D &mean,
938 const int maxIter) {
939
940 // get leaf-only iterators over all boxes in ws
941 auto function = std::make_unique<Geometry::MDAlgorithms::MDBoxMaskFunction>(pos, radiusSquared);
942 MDBoxBase<MDE, nd> *baseBox = ws->getBox();
943 MDBoxIterator<MDE, nd> MDiter(baseBox, 1000, true, function.get());
944
945 // get initial vector of events inside sphere
946 std::vector<std::pair<V3D, double>> peak_events;
947
948 do {
949 auto *box = dynamic_cast<MDBox<MDE, nd> *>(MDiter.getBox());
950 if (box && !box->getIsMasked()) {
951
952 // simple check whether box is defintely not contained
953 coord_t boxCenter[nd];
954 box->getCenter(boxCenter);
955 V3D displacement; // vector between peak pos and box center
956 coord_t rboxSq = 0; // dist from center to vertex sq
957 for (size_t d = 0; d < nd; ++d) {
958 auto dim = box->getExtents(d);
959 rboxSq += static_cast<coord_t>(0.25 * dim.getSize() * dim.getSize());
960 displacement[d] = pos[d] - static_cast<double>(boxCenter[d]);
961 }
962
963 if (displacement.norm() < static_cast<double>(sqrt(rboxSq)) + static_cast<double>(sqrt(radiusSquared))) {
964 // box MIGHT intersect peak spherical region so go through events
965 const std::vector<MDE> &events = box->getConstEvents();
966 auto bg = bgDensity / (static_cast<double>(events.size()) * (box->getInverseVolume()));
967 // For each event
968 for (const auto &evnt : events) {
969
970 coord_t center_array[nd];
971 for (size_t d = 0; d < nd; ++d) {
972 center_array[d] = evnt.getCenter(d);
973 }
974 coord_t out[1];
975 auto *cen_ptr = center_array; // pointer to first element
976 getRadiusSq.apply(cen_ptr, out);
977
978 if (evnt.getSignal() > bg && out[0] < radiusSquared) {
979 // need in V3D for matrix maths later
980 V3D center;
981 for (size_t d = 0; d < nd; ++d) {
982 center[d] = static_cast<double>(center_array[d]);
983 }
984 peak_events.emplace_back(center, evnt.getSignal() - bg);
985 }
986 }
987 }
988 }
989 box->releaseEvents();
990 } while (MDiter.next());
991 calcCovar(peak_events, pos, radiusSquared, qAxisIsFixed, useCentroid, eigenvects, eigenvals, mean, maxIter);
992}
993
1009void IntegratePeaksMD2::calcCovar(const std::vector<std::pair<V3D, double>> &peak_events, const V3D &pos,
1010 const coord_t &radiusSquared, const bool &qAxisIsFixed, const bool &useCentroid,
1011 std::vector<V3D> &eigenvects, std::vector<double> &eigenvals, V3D &mean,
1012 const int maxIter) {
1013
1014 size_t nd = 3;
1015
1016 // to calc threshold mdsq to exclude events over 3 stdevs away
1017 boost::math::chi_squared chisq(static_cast<double>(nd));
1018 auto mdsq_max = boost::math::quantile(chisq, 0.997);
1019 Matrix<double> invCov; // required to calc mdsq
1020 double prev_cov_det = DBL_MAX;
1021
1022 // initialise mean with pos
1023 mean = pos;
1024 Matrix<double> Pinv(nd, nd);
1025 if (qAxisIsFixed) {
1026 // transformation from Qlab to Qhat, vhat and uhat,
1027 getPinv(pos, Pinv);
1028 mean = Pinv * mean;
1029 }
1030 Matrix<double> cov_mat(nd, nd);
1031
1032 for (int nIter = 0; nIter < maxIter; nIter++) {
1033
1034 // reset on each loop
1035 cov_mat.zeroMatrix();
1036 double w_sum = 0; // sum of weights
1037 size_t nmasked = 0; // num masked events outside 3stdevs
1038 auto prev_pos = mean;
1039
1040 for (size_t ievent = 0; ievent < peak_events.size(); ievent++) {
1041
1042 const auto event = peak_events[ievent];
1043 auto center = event.first;
1044 if (qAxisIsFixed) {
1045 // transform coords to Q, uhat, vhat basis
1046 center = Pinv * center;
1047 }
1048
1049 bool useEvent = true;
1050 if (nIter > 0) {
1051 // check if point within 3 stdevs of mean (in MD frame)
1052 // prev_pos is the mean if useCentroid
1053 const auto displ = center - prev_pos;
1054 auto mdsq = displ.scalar_prod(invCov * displ);
1055 if (mdsq > mdsq_max) {
1056 // exclude points outside 3 stdevs
1057 useEvent = false;
1058 nmasked += 1;
1059 }
1060 }
1061
1062 // if prev cov_mat chec
1063 if (useEvent) {
1064 const auto signal = event.second;
1065 w_sum += signal;
1066
1067 if (useCentroid) {
1068 // update mean
1069 mean += (center - mean) * (signal / w_sum);
1070 }
1071
1072 // weight for variance
1073 auto wi = signal * (w_sum - signal) / w_sum;
1074 size_t istart = 0;
1075 if (qAxisIsFixed) {
1076 // variance along Q (skipped in next nested loops below)
1077 cov_mat[0][0] += wi * pow((center[0] - mean[0]), 2);
1078 istart = 1;
1079 }
1080 for (size_t row = istart; row < cov_mat.numRows(); ++row) {
1081 for (size_t col = istart; col < cov_mat.numRows(); ++col) {
1082 // symmeteric matrix
1083 if (row <= col) {
1084 auto cov = wi * (center[row] - mean[row]) * (center[col] - mean[col]);
1085 if (row == col) {
1086 cov_mat[row][col] += cov;
1087 } else {
1088 cov_mat[row][col] += cov;
1089 cov_mat[col][row] += cov;
1090 }
1091 }
1092 }
1093 }
1094 }
1095 }
1096 // normalise the covariance matrix
1097 cov_mat /= w_sum; // normalise by sum of weights
1098
1099 // check if another iteration is required
1100 bool anyMasked = (nIter > 0) ? (nmasked > 0) : true;
1101 // check if ellipsoid volume greater than sphere
1102 auto cov_det = cov_mat.determinant();
1103 bool isEllipVolGreater = cov_det > pow(static_cast<double>(radiusSquared / 9), 3);
1104 // check for convergence of variances
1105 bool isConverged = (cov_det > 0.95 * prev_cov_det);
1106
1107 if (!anyMasked || isEllipVolGreater || isConverged) {
1108 break;
1109 } else {
1110 prev_cov_det = cov_det;
1111 // required to eval Mahalanobis distance
1112 invCov = Matrix<double>(cov_mat);
1113 invCov.Invert();
1114 }
1115 }
1116
1117 if (qAxisIsFixed) {
1118 // transform back to MD basis
1119 Matrix<double> P(Pinv);
1120 P.Transpose();
1121 mean = P * mean;
1122 cov_mat = P * cov_mat * Pinv;
1123 }
1124 Matrix<double> evecs; // hold eigenvectors
1125 Matrix<double> evals; // hold eigenvals in diag
1126 cov_mat.Diagonalise(evecs, evals);
1127
1128 auto min_eval = evals[0][0];
1129 for (size_t d = 1; d < nd; ++d) {
1130 min_eval = std::min(min_eval, evals[d][d]);
1131 }
1132 if (min_eval > static_cast<double>(radiusSquared / 9)) {
1133 // haven't found good covar - set to spherical region
1134 evals.identityMatrix();
1135 evals = evals * (static_cast<double>(radiusSquared) / 9);
1136 g_log.warning() << "Covariance of peak at ";
1137 pos.printSelf(g_log.warning());
1138 g_log.warning() << " is not well constrained, it has been set to spherical" << std::endl;
1139 }
1140
1141 // convert to vectors for output
1142 eigenvals = evals.Diagonal();
1143 // set min eigenval to be small but non-zero (1e-6)
1144 // when no discernible peak above background
1145 std::replace_if(
1146 eigenvals.begin(), eigenvals.end(), [&](auto x) { return x < 1e-6; }, 1e-6);
1147
1148 // populate V3D vector of eigenvects (needed for ellipsoid shape)
1149 eigenvects = std::vector<V3D>(nd);
1150 for (size_t ivect = 0; ivect < nd; ++ivect) {
1151 eigenvects[ivect] = V3D(evecs[0][ivect], evecs[1][ivect], evecs[2][ivect]);
1152 }
1153}
1154
1165 // loop over 3 mutually-orthogonal vectors until get one with
1166 // a component perp to Q (within tolerance)
1167 double dotprod = 1;
1168 size_t ii = 0;
1169 V3D qhat = q / q.norm();
1170 V3D tmp;
1171 do {
1172 tmp = V3D(0, 0, 0); // reset u
1173 tmp[ii] = 1.0;
1174 dotprod = qhat.scalar_prod(tmp);
1175 ii++;
1176 } while (abs(dotprod) > 1.0 - 1e-6);
1177 // populate Pinv with basis vector rows
1178 Pinv.setRow(0, qhat);
1179 tmp = qhat.cross_prod(tmp);
1180 Pinv.setRow(1, tmp / tmp.norm());
1181 tmp = qhat.cross_prod(tmp);
1182 Pinv.setRow(2, tmp / tmp.norm());
1183}
1184
1185/*
1186 * Define edges for each instrument by masking. For CORELLI, tubes 1 and
1187 *16, and pixels 0 and 255. Get Q in the lab frame for every peak, call it
1188 *C For every point on the edge, the trajectory in reciprocal space is a
1189 *straight line, going through O=V3D(0,0,0). Calculate a point at a fixed
1190 *momentum, say k=1. Q in the lab frame
1191 *E=V3D(-k*sin(tt)*cos(ph),-k*sin(tt)*sin(ph),k-k*cos(ph)).
1192 * Normalize E to 1: E=E*(1./E.norm())
1193 *
1194 * @param inst: instrument
1195 */
1197 for (size_t i = 0; i < detectorInfo.size(); ++i) {
1198 if (detectorInfo.isMonitor(i))
1199 continue; // skip monitor
1200 if (!detectorInfo.isMasked(i))
1201 continue; // edge is masked so don't check if not masked
1202 const auto &det = detectorInfo.detector(i);
1203 double tt1 = det.getTwoTheta(V3D(0, 0, 0), V3D(0, 0, 1)); // two theta
1204 double ph1 = det.getPhi(); // phi
1205 V3D E1 = V3D(-std::sin(tt1) * std::cos(ph1), -std::sin(tt1) * std::sin(ph1),
1206 1. - std::cos(tt1)); // end of trajectory
1207 E1 = E1 * (1. / E1.norm()); // normalize
1208 E1Vec.emplace_back(E1);
1209 }
1210}
1211
1222 double edgeDist = DBL_MAX;
1223 for (auto &E1 : E1Vec) {
1224 V3D distv = QLabFrame - E1 * (QLabFrame.scalar_prod(E1)); // distance to the
1225 // trajectory as a
1226 // vector
1227 edgeDist = std::min(edgeDist, distv.norm()); // want smallest dist to peak
1228 }
1229 return edgeDist;
1230}
1231
1233 const std::string &property, const std::string &values) {
1234 // For CORELLI do not count as edge if next to another detector bank
1235 if (property == "Tube" && peakWS->getInstrument()->getName() == "CORELLI") {
1236 auto alg = createChildAlgorithm("MaskBTP");
1237 alg->setProperty<Workspace_sptr>("Workspace", peakWS);
1238 alg->setProperty("Bank", "1,7,12,17,22,27,30,59,63,69,74,79,84,89");
1239 alg->setProperty(property, "1");
1240 if (!alg->execute())
1241 throw std::runtime_error("MaskDetectors Child Algorithm has not executed successfully");
1242 auto alg2 = createChildAlgorithm("MaskBTP");
1243 alg2->setProperty<Workspace_sptr>("Workspace", peakWS);
1244 alg2->setProperty("Bank", "6,11,16,21,26,29,58,62,68,73,78,83,88,91");
1245 alg2->setProperty(property, "16");
1246 if (!alg2->execute())
1247 throw std::runtime_error("MaskDetectors Child Algorithm has not executed successfully");
1248 } else {
1249 auto alg = createChildAlgorithm("MaskBTP");
1250 alg->setProperty<Workspace_sptr>("Workspace", peakWS);
1251 alg->setProperty(property, values);
1252 if (!alg->execute())
1253 throw std::runtime_error("MaskDetectors Child Algorithm has not executed successfully");
1254 }
1255}
1256
1258 Mantid::Kernel::SpecialCoordinateSystem CoordinatesToUse, double radius) {
1259 // Get a direct ref to that peak.
1260 const IPeak &p1 = peakWS->getPeak(i);
1261 V3D pos1;
1262 if (CoordinatesToUse == Kernel::QLab) //"Q (lab frame)"
1263 pos1 = p1.getQLabFrame();
1264 else if (CoordinatesToUse == Kernel::QSample) //"Q (sample frame)"
1265 pos1 = p1.getQSampleFrame();
1266 else if (CoordinatesToUse == Kernel::HKL) //"HKL"
1267 pos1 = p1.getHKL();
1268 for (int j = i + 1; j < peakWS->getNumberPeaks(); ++j) {
1269 // Get a direct ref to rest of peaks peak.
1270 const IPeak &p2 = peakWS->getPeak(j);
1271 V3D pos2;
1272 if (CoordinatesToUse == Kernel::QLab) //"Q (lab frame)"
1273 pos2 = p2.getQLabFrame();
1274 else if (CoordinatesToUse == Kernel::QSample) //"Q (sample frame)"
1275 pos2 = p2.getQSampleFrame();
1276 else if (CoordinatesToUse == Kernel::HKL) //"HKL"
1277 pos2 = p2.getHKL();
1278 if (pos1.distance(pos2) < radius) {
1279 g_log.warning() << " Warning: Peak integration spheres for peaks " << i << " and " << j
1280 << " overlap. Distance between peaks is " << pos1.distance(pos2) << '\n';
1281 }
1282 }
1283}
1284
1285//----------------------------------------------------------------------------------------------
1289 inWS = getProperty("InputWorkspace");
1290
1292}
1293
1294double f_eval2(double x, void *params) {
1295 std::shared_ptr<const API::CompositeFunction> fun =
1296 *reinterpret_cast<std::shared_ptr<const API::CompositeFunction> *>(params);
1297 FunctionDomain1DVector domain(x);
1298 FunctionValues yval(domain);
1299 fun->function(domain, yval);
1300 return yval[0];
1301}
1302
1303} // namespace Mantid::MDAlgorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
gsl_vector * tmp
double sigma
Definition: GetAllEi.cpp:156
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...
#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.
double radius
Definition: Rasterize.cpp:31
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
void progress(double p, const std::string &msg="", double estimatedTime=0.0, int progressPrecision=0)
Sends ProgressNotification.
Definition: Algorithm.cpp:231
Unique SingleValueParameter Declaration for InputNDimensions.
virtual void apply(const coord_t *inputVector, coord_t *outVector) const =0
@ 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
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
A property class for workspaces.
Unique CoordCenterVectorParam type declaration for ndimensional coordinate centers.
Templated super-class of a multi-dimensional event "box".
Definition: MDBoxBase.h:50
void getCenter(coord_t *const center) const override
Get the center of the box.
Definition: MDBoxBase.h:210
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.
MDBoxIterator: iterate through MDBoxBase hierarchy down to a given maximum depth.
Definition: MDBoxIterator.h:28
bool next() override
Advance to the next cell.
MDBoxBase< MDE, nd > * getBox() const
Return a pointer to the current box pointed to by the iterator.
Definition: MDBoxIterator.h:39
Templated class for a multi-dimensional event "box".
Definition: MDBox.h:45
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.
PeakShapeEllipsoid : PeakShape representing a 3D ellipsoid.
PeakShapeSpherical : PeakShape for a spherical peak.
Geometry::DetectorInfo is an intermediate step towards a DetectorInfo that is part of Instrument-2....
Definition: DetectorInfo.h:49
bool isMasked(const size_t index) const
Returns true if the detector is masked.
const Geometry::IDetector & detector(const size_t index) const
Return a const reference to the detector with given index.
size_t size() const
Returns the size of the DetectorInfo, i.e., the number of detectors in the instrument.
bool isMonitor(const size_t index) const
Returns true if the detector is a monitor.
virtual double getTwoTheta(const Kernel::V3D &observer, const Kernel::V3D &axis) const =0
Gives the angle of this detector object with respect to an axis.
Structure describing a single-crystal peak.
Definition: IPeak.h:26
virtual double getH() const =0
virtual void setIntensity(double m_Intensity)=0
virtual void setPeakShape(Mantid::Geometry::PeakShape *shape)=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
Support for a property that holds an array of values.
Definition: ArrayProperty.h:28
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 setPropertySettings(const std::string &name, std::unique_ptr< IPropertySettings > settings)
void setPropertyGroup(const std::string &name, const std::string &group)
Set the group for a given property.
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
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
Numerical Matrix class.
Definition: Matrix.h:42
T determinant() const
Calculate the determinant.
Definition: Matrix.cpp:1048
T Invert()
LU inversion routine.
Definition: Matrix.cpp:924
void zeroMatrix()
Set the matrix to zero.
Definition: Matrix.cpp:646
void identityMatrix()
Makes the matrix an idenity matrix.
Definition: Matrix.cpp:661
int Diagonalise(Matrix< T > &, Matrix< T > &) const
(only Symmetric matrix)
Definition: Matrix.cpp:1319
void setRow(const size_t nRow, const std::vector< T > &newRow)
Definition: Matrix.cpp:686
std::vector< T > Diagonal() const
Returns a vector of the diagonal.
Definition: Matrix.cpp:1265
size_t numRows() const
Return the number of rows in the matrix.
Definition: Matrix.h:144
Matrix< T > & Transpose()
Transpose the matrix.
Definition: Matrix.cpp:793
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
double distance(const V3D &v) const noexcept
Calculates the distance between two vectors.
Definition: V3D.h:287
constexpr double scalar_prod(const V3D &v) const noexcept
Calculates the cross product.
Definition: V3D.h:274
void printSelf(std::ostream &) const
Prints a text representation of itself in format "[x,y,z]".
Definition: V3D.cpp:357
constexpr V3D cross_prod(const V3D &v) const noexcept
Cross product (this * argument)
Definition: V3D.h:278
double norm() const noexcept
Definition: V3D.h:263
std::vector< Kernel::V3D > E1Vec
save for all detector pixels
void runMaskDetectors(const Mantid::DataObjects::PeaksWorkspace_sptr &peakWS, const std::string &property, const std::string &values)
void calculateE1(const Geometry::DetectorInfo &detectorInfo)
Calculate if this Q is on a detector.
void getPinv(const Mantid::Kernel::V3D &q, Mantid::Kernel::Matrix< double > &Pinv)
Get the inverse of the matrix P.
void calcCovar(const std::vector< std::pair< Mantid::Kernel::V3D, double > > &peak_events, const Mantid::Kernel::V3D &pos, const coord_t &radiusSquared, const bool &qAxisIsFixed, const bool &useCentroid, std::vector< Mantid::Kernel::V3D > &eigenvects, std::vector< double > &eigenvals, Mantid::Kernel::V3D &mean, const int maxIter)
Calculate the covariance matrix of a spherical region and store the eigenvectors and eigenvalues that...
void exec() override
Run the algorithm.
double calculateDistanceToEdge(const Mantid::Kernel::V3D &QLabFrame)
Calculate if this Q is on a detector The distance from C to OE is given by dv=C-E*(C....
int version() const override
Algorithm's version for identification.
void checkOverlap(int i, const Mantid::API::IPeaksWorkspace_sptr &peakWS, Mantid::Kernel::SpecialCoordinateSystem CoordinatesToUse, double radius)
Check if peaks overlap.
void init() override
Initialise the properties.
void integrate(typename DataObjects::MDEventWorkspace< MDE, nd >::sptr ws)
Integrate the peaks of the workspace using parameters saved in the algorithm class.
void findEllipsoid(const typename DataObjects::MDEventWorkspace< MDE, nd >::sptr ws, const Mantid::API::CoordTransform &getRadiusSq, const Mantid::Kernel::V3D &pos, const coord_t &radiusSquared, const bool &qAxisBool, const bool &useCentroid, const double &bgDensity, std::vector< Mantid::Kernel::V3D > &eigenvects, std::vector< double > &eigenvals, Mantid::Kernel::V3D &mean, const int maxIter=1)
Calculate the covariance matrix of a spherical region and store the eigenvectors and eigenvalues that...
Mantid::API::IMDEventWorkspace_sptr inWS
Input MDEventWorkspace.
const std::string name() const override
Algorithm's name for identification.
std::map< std::string, std::string > validateInputs() override
Perform validation of ALL the input properties of the algorithm.
std::shared_ptr< IPeaksWorkspace > IPeaksWorkspace_sptr
shared pointer to Mantid::API::IPeaksWorkspace
std::shared_ptr< IAlgorithm > IAlgorithm_sptr
shared pointer to Mantid::API::IAlgorithm
std::complex< double > MANTID_API_DLL E1(std::complex< double > z)
Integral for Gamma.
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.
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
double f_eval2(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