31#include "MantidHistogramData/LinearGenerator.h"
43#include "boost/math/distributions.hpp"
47#include <gsl/gsl_integration.h>
65 "An input MDEventWorkspace.");
67 auto radiiValidator = std::make_shared<ArrayBoundedValidator<double>>();
68 radiiValidator->setLower(0.0);
69 radiiValidator->setLowerExclusive(
true);
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).");
76 radiiValidator->setLowerExclusive(
false);
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 = "
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 "
91 "If smaller than PeakRadius, no background measurement is done.");
94 "A PeaksWorkspace containing the peaks to integrate.");
97 "The output PeaksWorkspace will be a copy of the input PeaksWorkspace "
98 "with the peaks' integrated intensities.");
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)");
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.");
110 "Default is false. If true, "
111 "BackgroundOuterRadius + AdaptiveQMultiplier * **|Q|** and "
112 "BackgroundInnerRadius + AdaptiveQMultiplier * **|Q|**");
116 declareProperty(
"FixQAxis",
false,
"Fix one axis of ellipsoid to be along direction of Q.");
118 declareProperty(
"Cylinder",
false,
"Default is sphere. Use next five parameters for cylinder.");
121 "Length of cylinder in which to integrate (in the same units as the "
125 "Percent of CylinderLength that is background (20 is 20%)");
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.");
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.");
142 std::vector<std::string>(1,
"profiles")),
143 "Save (Optionally) as Isaw peaks file with profiles included");
146 "PeakRadius + AdaptiveQMultiplier * **|Q|** "
147 "so each peak has a "
148 "different integration radius. Q includes the 2*pi factor.");
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).");
157 "If this options is enabled, then the the top 1% of the "
158 "background will be removed"
159 "before the background subtraction.");
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 "
172 "Perform integration on estimated centroid not peak position "
173 "(ignored if all three peak radii are specified).");
175 auto maxIterValidator = std::make_shared<BoundedValidator<int>>();
176 maxIterValidator->setLower(1);
178 "Number of iterations in covariance estimation (ignored if all "
179 "peak radii are specified). 2-3 should be sufficient.");
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.");
189 std::string general_grp =
"General Inputs";
190 std::string cylin_grp =
"Cylindrical Integration";
191 std::string ellip_grp =
"Ellipsoid Integration";
225 "Cylinder", [](std::string ellipsoid,
const std::string &cylinder) {
227 if (ellipsoid ==
"1" && cylinder ==
"1") {
228 return std::string{
"0"};
234 "Ellipsoid", [](std::string cylinder,
const std::string &ellipsoid) {
236 if (cylinder ==
"1" && ellipsoid ==
"1") {
237 return std::string{
"0"};
258 std::make_unique<Kernel::EnabledWhenProperty>(
"IntegrateIfOnEdge",
IS_EQUAL_TO,
"1"));
262 std::map<std::string, std::string> result;
264 std::vector<double> PeakRadius =
getProperty(
"PeakRadius");
265 std::vector<double> BackgroundInnerRadius =
getProperty(
"BackgroundInnerRadius");
266 std::vector<double> BackgroundOuterRadius =
getProperty(
"BackgroundOuterRadius");
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();
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();
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();
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();
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();
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();
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();
323 throw std::invalid_argument(
"For now, we expect the input MDEventWorkspace "
324 "to have 3 dimensions only.");
331 if (peakWS != inPeakWS)
332 peakWS = inPeakWS->clone();
344 g_log.
error(
"Can't execute MaskBTP algorithm for this instrument to set "
345 "edge for IntegrateIfOnEdge option");
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");
360 bool manualEllip =
false;
361 if (PeakRadius.size() > 1) {
364 if (BackgroundInnerRadius.size() == 1)
365 BackgroundInnerRadius.resize(3, BackgroundInnerRadius[0]);
366 if (BackgroundOuterRadius.size() == 1)
367 BackgroundOuterRadius.resize(3, BackgroundOuterRadius[0]);
370 double minInnerRadius = PeakRadius[0];
371 for (
size_t r = 0; r < BackgroundInnerRadius.size(); r++) {
373 minInnerRadius = PeakRadius[r];
375 if (BackgroundInnerRadius[r] < minInnerRadius)
376 BackgroundInnerRadius[r] = minInnerRadius;
381 bool majorAxisLengthFixed =
getProperty(
"FixMajorAxisLength");
385 double cylinderLength =
getProperty(
"CylinderLength");
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]);
399 size_t histogramNumber = peakWS->getNumberPeaks();
401 wsProfile2D = std::dynamic_pointer_cast<Workspace2D>(wsProfile);
404 wsFit2D = std::dynamic_pointer_cast<Workspace2D>(wsFit);
407 wsDiff2D = std::dynamic_pointer_cast<Workspace2D>(wsDiff);
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) {
420 IPeak &p = peakWS->getPeak(i);
421 std::ostringstream label;
424 newAxis1Raw->setLabel(i, label.str());
425 newAxis2Raw->setLabel(i, label.str());
426 newAxis3Raw->setLabel(i, label.str());
429 double percentBackground =
getProperty(
"PercentBackground");
431 size_t peakMax = numSteps;
434 peakMin =
static_cast<size_t>(
static_cast<double>(numSteps) * percentBackground / 100.);
435 peakMax = numSteps - peakMin - 1;
436 size_t numPeakCh = peakMax - peakMin + 1;
437 size_t numBkgCh = numSteps - numPeakCh;
438 ratio =
static_cast<double>(numPeakCh) /
static_cast<double>(numBkgCh);
441 bool replaceIntensity =
getProperty(
"ReplaceIntensity");
442 bool integrateEdge =
getProperty(
"IntegrateIfOnEdge");
445 std::string profileFunction =
getProperty(
"ProfileFunction");
446 std::string integrationOption =
getProperty(
"IntegrationOption");
448 if (cylinderBool && profileFunction !=
"NoFit") {
449 std::string outFile =
getProperty(
"InputWorkspace");
450 outFile.append(profileFunction);
451 outFile.append(
".dat");
453 outFile = save_path + outFile;
454 out.open(outFile.c_str(), std::ofstream::out);
457 double volumeBkg = 4.0 / 3.0 * M_PI * (std::pow(BackgroundOuterRadius[0], 3) - std::pow(BackgroundOuterRadius[0], 3));
459 double volumeRadius = 4.0 / 3.0 * M_PI * std::pow(PeakRadius[0], 3);
462 int nPeaks = peakWS->getNumberPeaks();
466 for (
int i = 0; i < nPeaks; ++i) {
471 IPeak &p = peakWS->getPeak(i);
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) {
499 bool dimensionsUsed[nd];
501 for (
size_t d = 0;
d < nd; ++
d) {
502 dimensionsUsed[
d] =
true;
503 center[
d] =
static_cast<coord_t>(pos[
d]);
509 double background_total = 0.0;
513 if (adaptiveQMultiplier != 0.0) {
515 for (
size_t d = 0;
d < nd; ++
d) {
516 lenQpeak += center[
d] * center[
d];
518 lenQpeak = std::sqrt(lenQpeak);
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
527 PeakRadiusVector[i] = 0.0;
528 BackgroundInnerRadiusVector[i] = 0.0;
529 BackgroundOuterRadiusVector[i] = 0.0;
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());
541 new PeakShapeSpherical(PeakRadiusVector[i], BackgroundInnerRadiusVector[i], BackgroundOuterRadiusVector[i],
544 const double scaleFactor = pow(PeakRadiusVector[i], 3) /
545 (pow(BackgroundOuterRadiusVector[i], 3) - pow(BackgroundInnerRadiusVector[i], 3));
547 if (BackgroundOuterRadius[0] > PeakRadius[0]) {
550 getRadiusSq,
static_cast<coord_t>(pow(BackgroundOuterRadiusVector[i], 2)), bgSignal, bgErrorSquared,
551 static_cast<coord_t>(pow(BackgroundInnerRadiusVector[i], 2)), useOnePercentBackgroundCorrection);
553 bgSignal *= scaleFactor;
554 bgErrorSquared *= scaleFactor * scaleFactor;
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);
564 if (PeakRadius.size() == 1) {
565 V3D mean(0.0, 0.0, 0.0);
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) {
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;
577 translation = mean - pos;
579 for (
size_t d = 0;
d < 3; ++
d) {
580 center[
d] =
static_cast<coord_t>(mean[
d]);
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));
594 eigenvects, eigenvals);
596 if (PeakRadius.size() == 1) {
597 if (BackgroundOuterRadius[0] > PeakRadius[0]) {
602 getRadiusSq,
static_cast<coord_t>(pow(BackgroundOuterRadiusVector[i], 2)), bgSignal, bgErrorSquared,
603 static_cast<coord_t>(pow(BackgroundInnerRadiusVector[i], 2)), useOnePercentBackgroundCorrection);
606 bgSignal *= scaleFactor;
607 bgErrorSquared *= scaleFactor * scaleFactor;
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;
622 new PeakShapeEllipsoid(eigenvects, peakRadii, backgroundInnerRadii, backgroundOuterRadii,
623 CoordinatesToUse, this->
name(), this->
version(), translation);
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); });
635 if (BackgroundOuterRadiusVector[0] > PeakRadiusVector[0]) {
638 eigenvects, eigenvals_background_inner);
640 eigenvects, eigenvals_background_outer);
649 static_cast<coord_t>(pow(BackgroundInnerRadiusVector[i], 2)), bgSignalInner,
650 bgErrorSquaredInner, 0.0, useOnePercentBackgroundCorrection);
652 static_cast<coord_t>(pow(BackgroundOuterRadiusVector[i], 2)), bgSignalOuter,
653 bgErrorSquaredOuter, 0.0, useOnePercentBackgroundCorrection);
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]);
663 bgErrorSquared *= pow(scale, 2);
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;
683 new PeakShapeEllipsoid(eigenvects, peakRadii, backgroundInnerRadii, backgroundOuterRadii,
689 signal, errorSquared, 0.0 ,
690 useOnePercentBackgroundCorrection);
696 Counts signal_fit(numSteps);
700 static_cast<coord_t>(cylinderLength), signal, errorSquared,
701 signal_fit.mutableRawData());
704 if (BackgroundOuterRadius[0] > PeakRadius[0]) {
709 static_cast<coord_t>(cylinderLength), bgSignal, bgErrorSquared,
710 signal_fit.mutableRawData());
712 Points points(signal_fit.size(), LinearGenerator(0, 1));
713 wsProfile2D->setHistogram(i, points, signal_fit);
720 if (BackgroundInnerRadius[0] != PeakRadius[0]) {
722 static_cast<coord_t>(cylinderLength), interiorSignal, interiorErrorSquared,
723 signal_fit.mutableRawData());
727 interiorSignal = signal;
728 interiorErrorSquared = errorSquared;
732 bgSignal -= interiorSignal;
736 bgErrorSquared -= interiorErrorSquared;
738 const double radiusRatio = (PeakRadius[0] / BackgroundOuterRadius[0]);
739 const double peakVolume = radiusRatio * radiusRatio * cylinderLength;
743 const double interiorRatio = (BackgroundInnerRadius[0] / BackgroundOuterRadius[0]);
746 const double bgVolume = 1.0 - interiorRatio * interiorRatio * cylinderLength;
750 const double scaleFactor = peakVolume / bgVolume;
751 bgSignal *= scaleFactor;
752 bgErrorSquared *= scaleFactor * scaleFactor;
754 Points points(signal_fit.size(), LinearGenerator(0, 1));
755 wsProfile2D->setHistogram(i, points, signal_fit);
758 if (profileFunction ==
"NoFit") {
760 for (
size_t j = 0; j < numSteps; j++) {
761 if (j < peakMin || j > peakMax)
762 background_total = background_total + wsProfile2D->y(i)[j];
764 signal = signal + wsProfile2D->y(i)[j];
766 errorSquared = std::fabs(signal);
772 std::string myFunc = std::string(
"name=LinearBackground;name=") + profileFunction;
773 auto maxPeak = std::max_element(signal_fit.begin(), signal_fit.end());
775 std::ostringstream strs;
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");
785 fitAlgorithm->setProperty(
"CalcErrors",
true);
786 fitAlgorithm->setProperty(
"Function", myFunc);
787 fitAlgorithm->setProperty(
"InputWorkspace", wsProfile2D);
788 fitAlgorithm->setProperty(
"WorkspaceIndex",
static_cast<int>(i));
790 fitAlgorithm->executeAsChildAlg();
798 out << std::setw(20) <<
"spectrum"
800 for (
size_t j = 0; j < ifun->nParams(); ++j)
801 out << std::setw(20) << ifun->parameterName(j) <<
" ";
802 out << std::setw(20) <<
"chi2"
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) <<
" ";
810 double chi2 = fitAlgorithm->getProperty(
"OutputChi2overDoF");
811 out << std::setw(20) << std::fixed << std::setprecision(10) << chi2 <<
"\n";
813 std::shared_ptr<const CompositeFunction> fun = std::dynamic_pointer_cast<const CompositeFunction>(ifun);
815 const auto &
x = wsProfile2D->x(i);
816 wsFit2D->setSharedX(i, wsProfile2D->sharedX(i));
817 wsDiff2D->setSharedX(i, wsProfile2D->sharedX(i));
821 fun->function(domain, yy);
824 wsFit2D->mutableY(i) = funcValues;
825 wsDiff2D->setSharedY(i, wsProfile2D->sharedY(i));
826 wsDiff2D->mutableY(i) -= wsFit2D->y(i);
830 if (integrationOption ==
"Sum") {
831 for (
size_t j = peakMin; j <= peakMax; j++)
832 if (std::isfinite(yy[j]))
835 gsl_integration_workspace *w = gsl_integration_workspace_alloc(1000);
843 gsl_integration_qags(&F,
x[peakMin],
x[peakMax], 0, 1e-7, 1000, w, &signal, &
error);
845 gsl_integration_workspace_free(w);
847 errorSquared = std::fabs(signal);
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;
856 checkOverlap(i, peakWS, CoordinatesToUse, 2.0 * std::max(PeakRadiusVector[i], BackgroundOuterRadiusVector[i]));
858 if (signal != 0. || replaceIntensity) {
859 double edgeMultiplier = 1.0;
860 double peakMultiplier = 1.0;
862 if (edgeDist < BackgroundOuterRadius[0]) {
863 double e1 = BackgroundOuterRadius[0] - edgeDist;
865 double f1 = M_PI * std::pow(e1, 2) / 3 * (3 * BackgroundOuterRadius[0] - e1);
866 edgeMultiplier = volumeBkg / (volumeBkg - f1);
868 if (edgeDist < PeakRadius[0]) {
869 double sigma = PeakRadius[0] / 3.0;
871 double e1 = std::exp(-std::pow(edgeDist, 2) / (2 *
sigma *
sigma)) * PeakRadius[0];
873 double f1 = M_PI * std::pow(e1, 2) / 3 * (3 * PeakRadius[0] - e1);
874 peakMultiplier = volumeRadius / (volumeRadius - f1);
878 p.
setIntensity(peakMultiplier * signal - edgeMultiplier * (ratio * background_total + bgSignal));
880 edgeMultiplier * (ratio * ratio * std::fabs(background_total) + bgErrorSquared)));
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";
891 peakWS->mutableRun().addProperty(
"PeaksIntegrated", 1,
true);
893 peakWS->mutableRun().addProperty(
"PeakRadius", PeakRadiusVector,
true);
894 peakWS->mutableRun().addProperty(
"BackgroundInnerRadius", BackgroundInnerRadiusVector,
true);
895 peakWS->mutableRun().addProperty(
"BackgroundOuterRadius", BackgroundOuterRadiusVector,
true);
898 const std::string outfile =
getProperty(
"ProfilesFile");
899 if (outfile.length() > 0) {
904 g_log.
error(
"Can't locate SaveIsawPeaks algorithm");
907 alg->setProperty(
"InputWorkspace", peakWS);
908 alg->setProperty(
"ProfileWorkspace", wsProfile2D);
909 alg->setPropertyValue(
"Filename", outfile);
933template <
typename MDE,
size_t nd>
936 const bool &qAxisIsFixed,
const bool &useCentroid,
const double &bgDensity,
937 std::vector<V3D> &eigenvects, std::vector<double> &eigenvals,
V3D &mean,
941 auto function = std::make_unique<Geometry::MDAlgorithms::MDBoxMaskFunction>(pos, radiusSquared);
946 std::vector<std::pair<V3D, double>> peak_events;
950 if (box && !box->getIsMasked()) {
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]);
963 if (displacement.
norm() <
static_cast<double>(sqrt(rboxSq)) +
static_cast<double>(sqrt(radiusSquared))) {
965 const std::vector<MDE> &events = box->getConstEvents();
966 auto bg = bgDensity / (
static_cast<double>(events.size()) * (box->getInverseVolume()));
968 for (
const auto &evnt : events) {
971 for (
size_t d = 0;
d < nd; ++
d) {
972 center_array[
d] = evnt.getCenter(
d);
975 auto *cen_ptr = center_array;
976 getRadiusSq.
apply(cen_ptr, out);
978 if (evnt.getSignal() > bg && out[0] < radiusSquared) {
981 for (
size_t d = 0;
d < nd; ++
d) {
982 center[
d] =
static_cast<double>(center_array[
d]);
984 peak_events.emplace_back(center, evnt.getSignal() - bg);
989 box->releaseEvents();
990 }
while (MDiter.
next());
991 calcCovar(peak_events, pos, radiusSquared, qAxisIsFixed, useCentroid, eigenvects, eigenvals, mean, maxIter);
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) {
1017 boost::math::chi_squared chisq(
static_cast<double>(nd));
1018 auto mdsq_max = boost::math::quantile(chisq, 0.997);
1020 double prev_cov_det = DBL_MAX;
1032 for (
int nIter = 0; nIter < maxIter; nIter++) {
1038 auto prev_pos = mean;
1040 for (
size_t ievent = 0; ievent < peak_events.size(); ievent++) {
1042 const auto event = peak_events[ievent];
1043 auto center =
event.first;
1046 center = Pinv * center;
1049 bool useEvent =
true;
1053 const auto displ = center - prev_pos;
1054 auto mdsq = displ.scalar_prod(invCov * displ);
1055 if (mdsq > mdsq_max) {
1064 const auto signal =
event.second;
1069 mean += (center - mean) * (signal / w_sum);
1073 auto wi = signal * (w_sum - signal) / w_sum;
1077 cov_mat[0][0] += wi * pow((center[0] - mean[0]), 2);
1080 for (
size_t row = istart; row < cov_mat.
numRows(); ++row) {
1081 for (
size_t col = istart; col < cov_mat.
numRows(); ++col) {
1084 auto cov = wi * (center[row] - mean[row]) * (center[col] - mean[col]);
1086 cov_mat[row][col] += cov;
1088 cov_mat[row][col] += cov;
1089 cov_mat[col][row] += cov;
1100 bool anyMasked = (nIter > 0) ? (nmasked > 0) :
true;
1103 bool isEllipVolGreater = cov_det > pow(
static_cast<double>(radiusSquared / 9), 3);
1105 bool isConverged = (cov_det > 0.95 * prev_cov_det);
1107 if (!anyMasked || isEllipVolGreater || isConverged) {
1110 prev_cov_det = cov_det;
1122 cov_mat = P * cov_mat * Pinv;
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]);
1132 if (min_eval >
static_cast<double>(radiusSquared / 9)) {
1135 evals = evals * (
static_cast<double>(radiusSquared) / 9);
1138 g_log.
warning() <<
" is not well constrained, it has been set to spherical" << std::endl;
1146 eigenvals.begin(), eigenvals.end(), [&](
auto x) { return x < 1e-6; }, 1e-6);
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]);
1176 }
while (abs(dotprod) > 1.0 - 1e-6);
1197 for (
size_t i = 0; i < detectorInfo.
size(); ++i) {
1202 const auto &det = detectorInfo.
detector(i);
1204 double ph1 = det.getPhi();
1205 V3D E1 =
V3D(-std::sin(tt1) * std::cos(ph1), -std::sin(tt1) * std::sin(ph1),
1206 1. - std::cos(tt1));
1207 E1 =
E1 * (1. /
E1.norm());
1222 double edgeDist = DBL_MAX;
1227 edgeDist = std::min(edgeDist, distv.
norm());
1233 const std::string &property,
const std::string &values) {
1235 if (property ==
"Tube" && peakWS->getInstrument()->getName() ==
"CORELLI") {
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");
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");
1251 alg->setProperty(property, values);
1252 if (!alg->execute())
1253 throw std::runtime_error(
"MaskDetectors Child Algorithm has not executed successfully");
1260 const IPeak &p1 = peakWS->getPeak(i);
1268 for (
int j = i + 1; j < peakWS->getNumberPeaks(); ++j) {
1270 const IPeak &p2 = peakWS->getPeak(j);
1279 g_log.
warning() <<
" Warning: Peak integration spheres for peaks " << i <<
" and " << j
1280 <<
" overlap. Distance between peaks is " << pos1.
distance(pos2) <<
'\n';
1295 std::shared_ptr<const API::CompositeFunction> fun =
1296 *
reinterpret_cast<std::shared_ptr<const API::CompositeFunction> *
>(params);
1299 fun->function(domain, yval);
#define DECLARE_ALGORITHM(classname)
#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...
#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.
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
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.
void progress(double p, const std::string &msg="", double estimatedTime=0.0, int progressPrecision=0)
Sends ProgressNotification.
@ OptionalSave
to specify a file to write to but an empty string is
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...
Helper class for reporting progress from algorithms.
A property class for workspaces.
Templated super-class of a multi-dimensional event "box".
void getCenter(coord_t *const center) const override
Get the center of the box.
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.
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.
Templated class for a multi-dimensional event "box".
std::shared_ptr< MDEventWorkspace< MDE, nd > > sptr
Typedef for a shared pointer of this kind of event workspace.
MDBoxBase< MDE, nd > * getBox()
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....
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.
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.
Support for a property that holds an array of values.
Exception for when an item is not found in a collection.
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.
void error(const std::string &msg)
Logs at error level.
void warning(const std::string &msg)
Logs at warning level.
void information(const std::string &msg)
Logs at information level.
T determinant() const
Calculate the determinant.
T Invert()
LU inversion routine.
void zeroMatrix()
Set the matrix to zero.
void identityMatrix()
Makes the matrix an idenity matrix.
int Diagonalise(Matrix< T > &, Matrix< T > &) const
(only Symmetric matrix)
void setRow(const size_t nRow, const std::vector< T > &newRow)
std::vector< T > Diagonal() const
Returns a vector of the diagonal.
size_t numRows() const
Return the number of rows in the matrix.
Matrix< T > & Transpose()
Transpose the matrix.
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...
double distance(const V3D &v) const noexcept
Calculates the distance between two vectors.
constexpr double scalar_prod(const V3D &v) const noexcept
Calculates the cross product.
void printSelf(std::ostream &) const
Prints a text representation of itself in format "[x,y,z]".
constexpr V3D cross_prod(const V3D &v) const noexcept
Cross product (this * argument)
double norm() const noexcept
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
std::shared_ptr< IFunction > IFunction_sptr
shared pointer to the function base class
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++ (!)
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.
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,...
double signal_t
Typedef for the signal recorded in a MDBox, etc.
@ Input
An input workspace.
@ Output
An output workspace.