25#include "MantidHistogramData/LinearGenerator.h"
34#include <gsl/gsl_integration.h>
52 "An input MDEventWorkspace.");
54 std::vector<std::string> propOptions{
"Q (lab frame)",
"Q (sample frame)",
"HKL"};
55 declareProperty(
"CoordinatesToUse",
"Q (lab frame)", std::make_shared<StringListValidator>(propOptions),
56 "Ignored: algorithm uses the InputWorkspace's coordinates.");
59 "Fixed radius around each peak position in which to integrate (in the "
60 "same units as the workspace).");
63 "Inner radius to use to evaluate the background of the peak.\n"
64 "If smaller than PeakRadius, then we assume BackgroundInnerRadius = "
68 "Outer radius to use to evaluate the background of the peak.\n"
69 "The signal density around the peak (BackgroundInnerRadius < r < "
70 "BackgroundOuterRadius) is used to estimate the background under the "
72 "If smaller than PeakRadius, no background measurement is done.");
75 "A PeaksWorkspace containing the peaks to integrate.");
78 "The output PeaksWorkspace will be a copy of the input PeaksWorkspace "
79 "with the peaks' integrated intensities.");
82 "Always replace intensity in PeaksWorkspacem (default).\n"
83 "If false, then do not replace intensity if calculated value "
84 "is 0 (used for SNSSingleCrystalReduction)");
87 "Only warning if all of peak outer radius is not on detector (default).\n"
88 "If false, do not integrate if the outer radius is not on a detector.");
91 "Default is false. If true, all input radii are multiplied "
92 "by the magnitude of Q at the peak center so each peak has a "
93 "different integration radius.");
95 declareProperty(
"Cylinder",
false,
"Default is sphere. Use next five parameters for cylinder.");
98 "Length of cylinder in which to integrate (in the same units as the "
102 "Percent of CylinderLength that is background (20 is 20%)");
105 peakNames.emplace_back(
"NoFit");
106 declareProperty(
"ProfileFunction",
"Gaussian", std::make_shared<StringListValidator>(peakNames),
107 "Fitting function for profile that is used only with "
108 "Cylinder integration.");
110 std::vector<std::string> integrationOptions(2);
111 integrationOptions[0] =
"Sum";
112 integrationOptions[1] =
"GaussianQuadrature";
113 auto integrationvalidator = std::make_shared<StringListValidator>(integrationOptions);
114 declareProperty(
"IntegrationOption",
"GaussianQuadrature", integrationvalidator,
115 "Integration method for calculating intensity "
116 "used only with Cylinder integration.");
119 std::vector<std::string>(1,
"profiles")),
120 "Save (Optionally) as Isaw peaks file with profiles included");
130 throw std::invalid_argument(
"For now, we expect the input MDEventWorkspace "
131 "to have 3 dimensions only.");
138 if (peakWS != inPeakWS)
139 peakWS = inPeakWS->clone();
144 g_log.
warning() <<
" Warning" << CoordinatesToUse <<
'\n';
145 if (CoordinatesToUse ==
Kernel::QLab && CoordinatesToUseStr !=
"Q (lab frame)")
146 g_log.
warning() <<
"Warning: used Q (lab frame) coordinates for MD "
147 "workspace, not CoordinatesToUse from input \n";
148 else if (CoordinatesToUse ==
Kernel::QSample && CoordinatesToUseStr !=
"Q (sample frame)")
149 g_log.
warning() <<
"Warning: used Q (sample frame) coordinates for MD "
150 "workspace, not CoordinatesToUse from input \n";
151 else if (CoordinatesToUse ==
Kernel::HKL && CoordinatesToUseStr !=
"HKL")
152 g_log.
warning() <<
"Warning: used HKL coordinates for MD workspace, not "
153 "CoordinatesToUse from input \n";
158 double BackgroundOuterRadius =
getProperty(
"BackgroundOuterRadius");
160 double BackgroundInnerRadius =
getProperty(
"BackgroundInnerRadius");
161 if (BackgroundInnerRadius < PeakRadius)
162 BackgroundInnerRadius = PeakRadius;
164 double cylinderLength =
getProperty(
"CylinderLength");
168 bool adaptiveQRadius =
getProperty(
"AdaptiveQRadius");
169 std::vector<double> PeakRadiusVector(peakWS->getNumberPeaks(), PeakRadius);
170 std::vector<double> BackgroundInnerRadiusVector(peakWS->getNumberPeaks(), BackgroundInnerRadius);
171 std::vector<double> BackgroundOuterRadiusVector(peakWS->getNumberPeaks(), BackgroundOuterRadius);
174 size_t histogramNumber = peakWS->getNumberPeaks();
176 wsProfile2D = std::dynamic_pointer_cast<Workspace2D>(wsProfile);
179 wsFit2D = std::dynamic_pointer_cast<Workspace2D>(wsFit);
182 wsDiff2D = std::dynamic_pointer_cast<Workspace2D>(wsDiff);
184 auto newAxis1 = std::make_unique<TextAxis>(peakWS->getNumberPeaks());
185 auto newAxis2 = std::make_unique<TextAxis>(peakWS->getNumberPeaks());
186 auto newAxis3 = std::make_unique<TextAxis>(peakWS->getNumberPeaks());
187 auto newAxis1Raw = newAxis1.get();
188 auto newAxis2Raw = newAxis2.get();
189 auto newAxis3Raw = newAxis3.get();
190 wsProfile2D->replaceAxis(1, std::move(newAxis1));
191 wsFit2D->replaceAxis(1, std::move(newAxis2));
192 wsDiff2D->replaceAxis(1, std::move(newAxis3));
193 for (
int i = 0; i < peakWS->getNumberPeaks(); ++i) {
195 IPeak &p = peakWS->getPeak(i);
196 std::ostringstream label;
199 newAxis1Raw->setLabel(i, label.str());
200 newAxis2Raw->setLabel(i, label.str());
201 newAxis3Raw->setLabel(i, label.str());
204 double percentBackground =
getProperty(
"PercentBackground");
206 size_t peakMax = numSteps;
209 peakMin =
static_cast<size_t>(
static_cast<double>(numSteps) * percentBackground / 100.);
210 peakMax = numSteps - peakMin - 1;
211 size_t numPeakCh = peakMax - peakMin + 1;
212 size_t numBkgCh = numSteps - numPeakCh;
213 ratio =
static_cast<double>(numPeakCh) /
static_cast<double>(numBkgCh);
216 bool replaceIntensity =
getProperty(
"ReplaceIntensity");
217 bool integrateEdge =
getProperty(
"IntegrateIfOnEdge");
219 std::string profileFunction =
getProperty(
"ProfileFunction");
220 std::string integrationOption =
getProperty(
"IntegrationOption");
222 if (cylinderBool && profileFunction !=
"NoFit") {
223 std::string outFile =
getProperty(
"InputWorkspace");
224 outFile.append(profileFunction);
225 outFile.append(
".dat");
227 outFile = save_path + outFile;
228 out.open(outFile.c_str(), std::ofstream::out);
240 for (
int i = 0; i < peakWS->getNumberPeaks(); ++i) {
242 IPeak &p = peakWS->getPeak(i);
254 inst = peakWS->getInstrument();
256 if (BackgroundOuterRadius > PeakRadius) {
258 g_log.
warning() <<
"Warning: sphere/cylinder for integration is off "
259 "edge of detector for peak "
266 g_log.
warning() <<
"Warning: sphere/cylinder for integration is off "
267 "edge of detector for peak "
275 bool dimensionsUsed[nd];
277 for (
size_t d = 0;
d < nd; ++
d) {
278 dimensionsUsed[
d] =
true;
279 center[
d] =
static_cast<coord_t>(pos[
d]);
285 double background_total = 0.0;
289 if (adaptiveQRadius) {
291 for (
size_t d = 0;
d < nd;
d++) {
292 lenQpeak += center[
d] * center[
d];
294 lenQpeak = std::sqrt(lenQpeak);
296 PeakRadiusVector[i] = lenQpeak * PeakRadius;
297 BackgroundInnerRadiusVector[i] = lenQpeak * BackgroundInnerRadius;
298 BackgroundOuterRadiusVector[i] = lenQpeak * BackgroundOuterRadius;
301 if (
auto *shapeablePeak =
dynamic_cast<Peak *
>(&p)) {
304 new PeakShapeSpherical(PeakRadiusVector[i], BackgroundInnerRadiusVector[i], BackgroundOuterRadiusVector[i],
306 shapeablePeak->setPeakShape(sphereShape);
315 if (BackgroundOuterRadius > PeakRadius) {
318 sphere,
static_cast<coord_t>(lenQpeak * BackgroundOuterRadius * lenQpeak * BackgroundOuterRadius), bgSignal,
326 if (BackgroundInnerRadius != PeakRadius) {
328 sphere,
static_cast<coord_t>(lenQpeak * BackgroundInnerRadius * lenQpeak * BackgroundInnerRadius),
329 interiorSignal, interiorErrorSquared);
332 interiorSignal = signal;
333 interiorErrorSquared = errorSquared;
337 bgSignal -= interiorSignal;
340 bgErrorSquared -= interiorErrorSquared;
343 const double radiusRatio = (PeakRadius / BackgroundOuterRadius);
344 const double peakVolume = radiusRatio * radiusRatio * radiusRatio;
347 const double interiorRatio = (BackgroundInnerRadius / BackgroundOuterRadius);
350 const double bgVolume = 1.0 - interiorRatio * interiorRatio * interiorRatio;
354 const double scaleFactor = peakVolume / bgVolume;
355 bgSignal *= scaleFactor;
356 bgErrorSquared *= scaleFactor * scaleFactor;
362 Counts signal_fit(numSteps);
366 signal, errorSquared, signal_fit.mutableRawData());
368 Points points(signal_fit.size(), LinearGenerator(0, 1));
369 wsProfile2D->setHistogram(i, points, signal_fit);
372 if (BackgroundOuterRadius > PeakRadius) {
377 static_cast<coord_t>(cylinderLength), bgSignal, bgErrorSquared,
378 signal_fit.mutableRawData());
380 wsProfile2D->setHistogram(i, points, signal_fit);
387 if (BackgroundInnerRadius != PeakRadius) {
389 static_cast<coord_t>(cylinderLength), interiorSignal, interiorErrorSquared,
390 signal_fit.mutableRawData());
393 interiorSignal = signal;
394 interiorErrorSquared = errorSquared;
398 bgSignal -= interiorSignal;
401 bgErrorSquared -= interiorErrorSquared;
403 const double radiusRatio = (PeakRadius / BackgroundOuterRadius);
404 const double peakVolume = radiusRatio * radiusRatio * cylinderLength;
407 const double interiorRatio = (BackgroundInnerRadius / BackgroundOuterRadius);
410 const double bgVolume = 1.0 - interiorRatio * interiorRatio * cylinderLength;
414 const double scaleFactor = peakVolume / bgVolume;
415 bgSignal *= scaleFactor;
416 bgErrorSquared *= scaleFactor * scaleFactor;
418 wsProfile2D->setHistogram(i, points, signal_fit);
421 if (profileFunction ==
"NoFit") {
422 auto &
y = wsProfile2D->y(i);
424 signal =
y.sum(peakMin, peakMax);
426 background_total +=
y.sum(0, peakMin);
427 background_total +=
y.sum(peakMax);
428 errorSquared = std::abs(signal);
431 findpeaks->setProperty(
"InputWorkspace", wsProfile2D);
432 findpeaks->setProperty<
int>(
"FWHM", 7);
433 findpeaks->setProperty<
int>(
"Tolerance", 4);
435 findpeaks->setProperty(
"WorkspaceIndex",
static_cast<int>(i));
438 findpeaks->setProperty<std::string>(
"PeakFunction", profileFunction);
440 findpeaks->setProperty<std::string>(
"BackgroundType",
"Quadratic");
441 findpeaks->setProperty<
bool>(
"HighBackground",
true);
442 findpeaks->setProperty<
bool>(
"RawPeakParameters",
true);
443 std::vector<double> peakPosToFit;
444 peakPosToFit.emplace_back(
static_cast<double>(numSteps) / 2.0);
445 findpeaks->setProperty(
"PeakPositions", peakPosToFit);
446 findpeaks->setProperty<
int>(
"MinGuessedPeakWidth", 4);
447 findpeaks->setProperty<
int>(
"MaxGuessedPeakWidth", 4);
449 findpeaks->executeAsChildAlg();
451 g_log.
error(
"Can't execute FindPeaks algorithm");
456 if (paramws->rowCount() < 1)
458 std::ostringstream fun_str;
459 fun_str <<
"name=" << profileFunction;
461 size_t numcols = paramws->columnCount();
462 std::vector<std::string> paramsName = paramws->getColumnNames();
463 std::vector<double> paramsValue;
467 for (
size_t j = 1; j < numcols; ++j) {
470 if (j == numcols - 4)
471 fun_str <<
";name=Quadratic";
475 if (j > 0 && j < numcols - 1)
476 fun_str <<
"," << paramsName[j] <<
"=" << parvalue;
477 paramsValue.emplace_back(parvalue);
480 for (
size_t j = 0; j < numcols; ++j)
481 out << std::setw(20) << paramsName[j] <<
" ";
484 out << std::setw(20) << i;
485 for (
size_t j = 0; j < numcols - 1; ++j)
486 out << std::setw(20) << std::fixed << std::setprecision(10) << paramsValue[j] <<
" ";
492 std::shared_ptr<const CompositeFunction> fun = std::dynamic_pointer_cast<const CompositeFunction>(ifun);
493 const auto &
x = wsProfile2D->x(i);
494 wsFit2D->setSharedX(i, wsProfile2D->sharedX(i));
495 wsDiff2D->setSharedX(i, wsProfile2D->sharedX(i));
499 fun->function(domain, yy);
502 wsFit2D->mutableY(i) = funcValues;
503 wsDiff2D->setSharedY(i, wsProfile2D->sharedY(i));
504 wsDiff2D->mutableY(i) -= wsFit2D->y(i);
508 if (integrationOption ==
"Sum") {
510 for (
size_t j = peakMin; j <= peakMax; j++)
511 if (std::isfinite(yy[j]))
514 gsl_integration_workspace *w = gsl_integration_workspace_alloc(1000);
522 gsl_integration_qags(&F,
x[peakMin],
x[peakMax], 0, 1e-7, 1000, w, &signal, &
error);
524 gsl_integration_workspace_free(w);
526 errorSquared = std::fabs(signal);
528 for (
size_t j = 0; j < numSteps; j++) {
530 paramsValue[numcols - 3] *
x[j] *
x[j] + paramsValue[numcols - 4] *
x[j] + paramsValue[numcols - 5];
531 if (j < peakMin || j > peakMax)
532 background_total = background_total + background;
537 if (signal != 0. || replaceIntensity) {
538 p.
setIntensity(signal - ratio * background_total - bgSignal);
539 p.
setSigmaIntensity(sqrt(errorSquared + ratio * ratio * std::fabs(background_total) + bgErrorSquared));
542 g_log.
information() <<
"Peak " << i <<
" at " << pos <<
": signal " << signal <<
" (sig^2 " << errorSquared
543 <<
"), with background " << bgSignal + ratio * background_total <<
" (sig^2 "
544 << bgErrorSquared + ratio * ratio * std::fabs(background_total) <<
") subtracted.\n";
548 peakWS->mutableRun().addProperty(
"PeaksIntegrated", 1,
true);
550 peakWS->mutableRun().addProperty(
"PeakRadius", PeakRadiusVector,
true);
551 peakWS->mutableRun().addProperty(
"BackgroundInnerRadius", BackgroundInnerRadiusVector,
true);
552 peakWS->mutableRun().addProperty(
"BackgroundOuterRadius", BackgroundOuterRadiusVector,
true);
555 const std::string outfile =
getProperty(
"ProfilesFile");
556 if (outfile.length() > 0) {
561 g_log.
error(
"Can't locate SaveIsawPeaks algorithm");
564 alg->setProperty(
"InputWorkspace", peakWS);
565 alg->setProperty(
"ProfileWorkspace", wsProfile2D);
566 alg->setPropertyValue(
"Filename", outfile);
580 const int nAngles = 8;
581 double dAngles =
static_cast<coord_t>(nAngles);
583 for (
int i = 0; i < nAngles; ++i) {
584 double theta = (2 * M_PI) / dAngles * i;
585 for (
int j = 0; j < nAngles; ++j) {
586 double phi = (2 * M_PI) / dAngles * j;
589 V3D edge =
V3D(QLabFrame.
X() + r * std::cos(theta) * std::sin(phi),
590 QLabFrame.
Y() + r * std::sin(theta) * std::sin(phi), QLabFrame.
Z() + r * std::cos(phi));
616 std::shared_ptr<const API::CompositeFunction> fun =
617 *
reinterpret_cast<std::shared_ptr<const API::CompositeFunction> *
>(params);
620 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...
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
std::string getPropertyValue(const std::string &name) const override
Get the value of a property as a string.
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.
@ 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...
TableRow represents a row in a TableWorkspace.
A property class for workspaces.
void integrateSphere(Mantid::API::CoordTransform &radiusTransform, const coord_t radiusSquared, signal_t &signal, signal_t &errorSquared, const coord_t innerRadiusSquared=0.0, const bool useOnePercentBackgroundCorrection=true) const override=0
Sphere (peak) integration.
void integrateCylinder(Mantid::API::CoordTransform &radiusTransform, const coord_t radius, const coord_t length, signal_t &signal, signal_t &errorSquared, std::vector< signal_t > &signal_fit) const override=0
Cylinder (peak) integration.
std::shared_ptr< MDEventWorkspace< MDE, nd > > sptr
Typedef for a shared pointer of this kind of event workspace.
MDBoxBase< MDE, nd > * getBox()
Kernel::SpecialCoordinateSystem getSpecialCoordinateSystem() const override
Get the coordinate system.
PeakShapeSpherical : PeakShape for a spherical peak.
Structure describing a single-crystal peak.
bool findDetector()
After creating a peak using the Q in the lab frame, the detPos is set to the direction of the detecto...
Structure describing a single-crystal peak.
virtual double getH() const =0
virtual void setIntensity(double m_Intensity)=0
virtual double getK() const =0
virtual int getRunNumber() const =0
virtual Mantid::Kernel::V3D getQSampleFrame() const =0
virtual double getL() const =0
virtual void setSigmaIntensity(double m_SigmaIntensity)=0
virtual Mantid::Kernel::V3D getQLabFrame() const =0
virtual Mantid::Kernel::V3D getHKL() const =0
PeakShape : Abstract type to describes the shape of a peak.
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 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.
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...
constexpr double X() const noexcept
Get x.
constexpr double Y() const noexcept
Get y.
constexpr double Z() const noexcept
Get z.
void integrate(typename DataObjects::MDEventWorkspace< MDE, nd >::sptr ws)
Integrate the peaks of the workspace using parameters saved in the algorithm class.
int version() const override
Algorithm's version for identification.
bool detectorQ(const Mantid::Kernel::V3D &QLabFrame, double r)
Calculate if this Q is on a detector.
void exec() override
Run the algorithm.
Mantid::API::IMDEventWorkspace_sptr inWS
Input MDEventWorkspace.
void init() override
Initialise the properties.
Geometry::Instrument_const_sptr inst
Instrument reference.
const std::string name() const override
Algorithm's name for identification.
std::shared_ptr< IAlgorithm > IAlgorithm_sptr
shared pointer to Mantid::API::IAlgorithm
std::shared_ptr< ITableWorkspace > ITableWorkspace_sptr
shared pointer to Mantid::API::ITableWorkspace
std::shared_ptr< Workspace > Workspace_sptr
shared pointer to Mantid::API::Workspace
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.
double f_eval(double x, void *params)
float coord_t
Typedef for the data type to use for coordinate axes in MD objects such as MDBox, MDEventWorkspace,...
double signal_t
Typedef for the signal recorded in a MDBox, etc.
@ Input
An input workspace.
@ Output
An output workspace.