Mantid
Loading...
Searching...
No Matches
PeakIntensityVsRadius.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 +
8#include "MantidAPI/Axis.h"
11#include "MantidAPI/TextAxis.h"
14
16
18
19using namespace Mantid::Kernel;
20using namespace Mantid::API;
21using namespace Mantid::DataObjects;
22
23namespace Mantid::Crystal {
24
25// Register the algorithm into the AlgorithmFactory
26DECLARE_ALGORITHM(PeakIntensityVsRadius)
27
28//----------------------------------------------------------------------------------------------
30const std::string PeakIntensityVsRadius::name() const { return "PeakIntensityVsRadius"; }
31
33int PeakIntensityVsRadius::version() const { return 1; }
34
36const std::string PeakIntensityVsRadius::category() const { return "Crystal\\Integration"; }
37
38//----------------------------------------------------------------------------------------------
39
40//----------------------------------------------------------------------------------------------
44 declareProperty(std::make_unique<WorkspaceProperty<IMDEventWorkspace>>("InputWorkspace", "", Direction::Input),
45 "An input MDEventWorkspace containing the SCD data.");
46 declareProperty(std::make_unique<WorkspaceProperty<PeaksWorkspace>>("PeaksWorkspace", "", Direction::Input),
47 "The list of peaks to integrate, matching the InputWorkspace.");
48
49 declareProperty("RadiusStart", 0.0, "Radius at which to start integrating.");
50 declareProperty("RadiusEnd", 1.0, "Radius at which to stop integrating.");
51 auto mustBePositive = std::make_shared<BoundedValidator<int>>();
52 mustBePositive->setLower(0);
53 declareProperty("NumSteps", 10, mustBePositive, "Number of steps, between start and end, to calculate radius.");
54
55 declareProperty("BackgroundInnerFactor", 0.0,
56 "For background subtraction: the peak radius will be multiplied\n"
57 "by this factor and passed to the BackgroundInnerRadius parameter.\n"
58 "Default 0.0 (no background).");
59 declareProperty("BackgroundOuterFactor", 0.0,
60 "For background subtraction: the peak radius will be multiplied\n"
61 "by this factor and passed to the BackgroundOuterRadius parameter.\n"
62 "Default 0.0 (no background).");
63
64 setPropertyGroup("BackgroundInnerFactor", "Variable Background Shell");
65 setPropertyGroup("BackgroundOuterFactor", "Variable Background Shell");
66
67 declareProperty("BackgroundInnerRadius", 0.0,
68 "For background subtraction:\n"
69 "Specify a fixed BackgroundInnerRadius, which does not "
70 "change with PeakRadius.\n"
71 "Default 0.0 (no background).");
72
73 declareProperty("BackgroundOuterRadius", 0.0,
74 "For background subtraction:\n"
75 "Specify a fixed BackgroundOuterRadius, which does not "
76 "change with PeakRadius.\n"
77 "Default 0.0 (no background).");
78
79 setPropertyGroup("BackgroundInnerRadius", "Fixed Background Shell");
80 setPropertyGroup("BackgroundOuterRadius", "Fixed Background Shell");
81
82 declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
83 "An output workspace2D containing intensity vs radius.");
84 declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace2", "NumberPeaksIntegrated", Direction::Output),
85 "An output workspace2D containing number of peaks at levels "
86 "of I/sigI vs radius.");
87}
88
90std::map<std::string, std::string> PeakIntensityVsRadius::validateInputs() {
91 std::map<std::string, std::string> out;
92 double BackgroundInnerFactor = getProperty("BackgroundInnerFactor");
93 double BackgroundOuterFactor = getProperty("BackgroundOuterFactor");
94 double BackgroundInnerRadius = getProperty("BackgroundInnerRadius");
95 double BackgroundOuterRadius = getProperty("BackgroundOuterRadius");
96 if ((BackgroundInnerRadius > 0) && (BackgroundInnerFactor > 0)) {
97 out["BackgroundInnerRadius"] = "Do not specify both BackgroundInnerRadius and BackgroundInnerFactor.";
98 out["BackgroundInnerFactor"] = "Do not specify both BackgroundInnerRadius and BackgroundInnerFactor.";
99 }
100 if ((BackgroundOuterRadius > 0) && (BackgroundOuterFactor > 0)) {
101 out["BackgroundOuterRadius"] = "Do not specify both BackgroundOuterRadius and BackgroundOuterFactor.";
102 out["BackgroundOuterFactor"] = "Do not specify both BackgroundOuterRadius and BackgroundOuterFactor.";
103 }
104 return out;
105}
106
107//----------------------------------------------------------------------------------------------
111 IMDEventWorkspace_sptr inWS = getProperty("InputWorkspace");
112 PeaksWorkspace_sptr peaksWS = getProperty("PeaksWorkspace");
113 double RadiusStart = getProperty("RadiusStart");
114 double RadiusEnd = getProperty("RadiusEnd");
115 double BackgroundInnerFactor = getProperty("BackgroundInnerFactor");
116 double BackgroundOuterFactor = getProperty("BackgroundOuterFactor");
117 double BackgroundInnerRadius = getProperty("BackgroundInnerRadius");
118 double BackgroundOuterRadius = getProperty("BackgroundOuterRadius");
119 int NumSteps = getProperty("NumSteps");
120
121 // Create a workspace with one spectrum per peak, and one point per radius
122 // step
124 WorkspaceFactory::Instance().create("Workspace2D", peaksWS->getNumberPeaks(), NumSteps, NumSteps);
125
126 // Create a text axis for axis(1), with H K L of each peak
127 auto ax = std::make_unique<TextAxis>(outWS->getNumberHistograms());
128 for (int i = 0; i < peaksWS->getNumberPeaks(); i++) {
129 V3D hkl = peaksWS->getPeak(i).getHKL();
130 hkl.round(); // Round HKL to make the string prettier
131 ax->setLabel(size_t(i), hkl.toString());
132 }
133 outWS->replaceAxis(1, std::move(ax));
134
135 MatrixWorkspace_sptr outWS2 = WorkspaceFactory::Instance().create("Workspace2D", 4, NumSteps, NumSteps);
136 // Create a text axis for axis(1), with H K L of each peak
137 auto ax2 = std::make_unique<TextAxis>(outWS2->getNumberHistograms());
138 ax2->setLabel(0, "I/SigI=2");
139 ax2->setLabel(1, "I/SigI=3");
140 ax2->setLabel(2, "I/SigI=5");
141 ax2->setLabel(3, "I/SigI=10");
142 outWS2->replaceAxis(1, std::move(ax2));
143
144 Progress prog(this, 0.0, 1.0, NumSteps);
145 double progStep = 1.0 / double(NumSteps);
146 for (int step = 0; step < NumSteps; step++) {
147 // Step from RadiusStart to RadiusEnd, inclusively
148 double radius = RadiusStart + double(step) * (RadiusEnd - RadiusStart) / (double(NumSteps - 1));
149 g_log.debug() << "Integrating radius " << radius << '\n';
150 prog.report("Radius " + Kernel::Strings::toString(radius));
151
152 double OuterRadius = 0;
153 if (BackgroundOuterRadius > 0)
154 OuterRadius = BackgroundOuterRadius;
155 if (BackgroundOuterFactor > 0)
156 OuterRadius = BackgroundOuterFactor * radius;
157 double InnerRadius = 0;
158 if (BackgroundInnerRadius > 0)
159 InnerRadius = BackgroundInnerRadius;
160 if (BackgroundInnerFactor > 0)
161 InnerRadius = BackgroundInnerFactor * radius;
162
163 // Run the integrate algo with this background
164 auto alg = createChildAlgorithm("IntegratePeaksMD", progStep * double(step), progStep *double(step + 1), false);
165 alg->setProperty("InputWorkspace", inWS);
166 alg->setProperty("PeaksWorkspace", peaksWS);
167 alg->setProperty<std::vector<double>>("PeakRadius", {radius});
168 alg->setProperty<std::vector<double>>("BackgroundOuterRadius", {OuterRadius});
169 alg->setProperty<std::vector<double>>("BackgroundInnerRadius", {InnerRadius});
170 alg->setPropertyValue("OutputWorkspace", "__tmp__PeakIntensityVsRadius");
171 alg->execute();
172 if (alg->isExecuted()) {
173 size_t ISigI2 = 0;
174 size_t ISigI3 = 0;
175 size_t ISigI5 = 0;
176 size_t ISigI10 = 0;
177 for (int i = 0; i < 4; i++) {
178 outWS2->mutableX(i)[step] = radius;
179 }
180 // Retrieve the integrated workspace
181 IPeaksWorkspace_sptr outPeaks = alg->getProperty("OutputWorkspace");
182 for (int i = 0; i < outPeaks->getNumberPeaks(); i++) {
183 auto wi = size_t(i); // workspace index in output
184 Geometry::IPeak &p = outPeaks->getPeak(i);
185 outWS->mutableX(wi)[step] = radius;
186 outWS->mutableY(wi)[step] = p.getIntensity();
187 outWS->mutableE(wi)[step] = p.getSigmaIntensity();
188 double ISigI = p.getIntensity() / p.getSigmaIntensity();
189 if (ISigI > 10.0) {
190 ISigI2++;
191 ISigI3++;
192 ISigI5++;
193 ISigI10++;
194 } else if (ISigI > 5.0) {
195 ISigI2++;
196 ISigI3++;
197 ISigI5++;
198 } else if (ISigI > 3.0) {
199 ISigI2++;
200 ISigI3++;
201 } else if (ISigI > 2.0) {
202 ISigI2++;
203 }
204 }
205 outWS2->mutableY(0)[step] = static_cast<double>(ISigI2);
206 outWS2->mutableY(1)[step] = static_cast<double>(ISigI3);
207 outWS2->mutableY(2)[step] = static_cast<double>(ISigI5);
208 outWS2->mutableY(3)[step] = static_cast<double>(ISigI10);
209 } else {
210 // TODO: Clear the point
211 }
212 }
213
214 // Fix units and labels
215 outWS->setYUnit("Integrated Intensity");
216 outWS->getAxis(0)->title() = "Radius";
217 outWS2->setYUnit("Number Peaks");
218 outWS2->getAxis(0)->title() = "Radius";
219
220 setProperty("OutputWorkspace", outWS);
221 setProperty("OutputWorkspace2", outWS2);
222}
223
224} // namespace Mantid::Crystal
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
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
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
A property class for workspaces.
Calculate the integrated intensity of peaks vs integration radius.
int version() const override
Algorithm's version for identification.
std::map< std::string, std::string > validateInputs() override
Check for reasonable values.
void exec() override
Execute the algorithm.
const std::string category() const override
Algorithm's category for identification.
void init() override
Initialize the algorithm's properties.
Structure describing a single-crystal peak.
Definition: IPeak.h:26
virtual double getSigmaIntensity() const =0
virtual double getIntensity() const =0
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
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 report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
Definition: ProgressBase.h:51
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
Class for 3D vectors.
Definition: V3D.h:34
std::string toString() const
Definition: V3D.cpp:340
void round() noexcept
Round each component to the nearest integer.
Definition: V3D.cpp:140
std::shared_ptr< IMDEventWorkspace > IMDEventWorkspace_sptr
Shared pointer to Mantid::API::IMDEventWorkspace.
std::shared_ptr< IPeaksWorkspace > IPeaksWorkspace_sptr
shared pointer to Mantid::API::IPeaksWorkspace
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
std::shared_ptr< PeaksWorkspace > PeaksWorkspace_sptr
Typedef for a shared pointer to a peaks workspace.
std::string toString(const T &value)
Convert a number to a string.
Definition: Strings.cpp:703
STL namespace.
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54