Mantid
Loading...
Searching...
No Matches
IntegratePeaksMDHKL.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 +
18#include <algorithm>
19#include <boost/math/special_functions/round.hpp>
20#include <limits>
21
22namespace Mantid::MDAlgorithms {
23
25// using Mantid::API::WorkspaceProperty;
26using namespace Mantid::DataObjects;
27using namespace Mantid::API;
28using namespace Mantid::Kernel;
29using namespace Mantid::Geometry;
30
31// Register the algorithm into the AlgorithmFactory
32DECLARE_ALGORITHM(IntegratePeaksMDHKL)
33
34//----------------------------------------------------------------------------------------------
38void IntegratePeaksMDHKL::init() {
39 declareProperty(std::make_unique<WorkspaceProperty<IMDWorkspace>>("InputWorkspace", "", Direction::Input),
40 "An input Sample MDHistoWorkspace or MDEventWorkspace in HKL.");
41 declareProperty("DeltaHKL", 0.5, "Distance from integer HKL to integrate peak.");
42 declareProperty("GridPoints", 201, "Number of grid points for each dimension of HKL box.");
43 declareProperty("NeighborPoints", 10,
44 "Number of points in 5^3 surrounding "
45 "points above intensity threshold for "
46 "point to be part of peak.");
47 auto fluxValidator = std::make_shared<CompositeValidator>();
48 fluxValidator->add<WorkspaceUnitValidator>("Momentum");
49 fluxValidator->add<InstrumentValidator>();
50 fluxValidator->add<CommonBinsValidator>();
51 auto solidAngleValidator = fluxValidator->clone();
52
53 declareProperty(std::make_unique<WorkspaceProperty<>>("FluxWorkspace", "", Direction::Input, PropertyMode::Optional,
54 fluxValidator),
55 "An optional input workspace containing momentum dependent flux for "
56 "normalization.");
57 declareProperty(std::make_unique<WorkspaceProperty<>>("SolidAngleWorkspace", "", Direction::Input,
58 PropertyMode::Optional, solidAngleValidator),
59 "An optional input workspace containing momentum integrated "
60 "vanadium for normalization "
61 "(a measure of the solid angle).");
62
63 declareProperty(std::make_unique<WorkspaceProperty<PeaksWorkspace>>("PeaksWorkspace", "", Direction::Input),
64 "A PeaksWorkspace containing the peaks to integrate.");
65
66 declareProperty(std::make_unique<WorkspaceProperty<PeaksWorkspace>>("OutputWorkspace", "", Direction::Output),
67 "The output PeaksWorkspace will be a copy of the input PeaksWorkspace "
68 "with the peaks' integrated intensities.");
69 declareProperty(std::make_unique<PropertyWithValue<double>>("BackgroundInnerRadius", EMPTY_DBL(), Direction::Input),
70 "Optional:Inner radius to use to evaluate the background of the peak.\n"
71 "If omitted background is region of HKL box - peak. ");
72
73 declareProperty(std::make_unique<PropertyWithValue<double>>("BackgroundOuterRadius", EMPTY_DBL(), Direction::Input),
74 "Optional:Outer radius to use to evaluate the background of the peak.\n"
75 "The signal density around the peak (BackgroundInnerRadius < r < "
76 "BackgroundOuterRadius) is used to estimate the background under the "
77 "peak.\n"
78 "If omitted background is region of HKL box - peak.");
79}
80
81//----------------------------------------------------------------------------------------------
86 IMDWorkspace_sptr m_inputWS = getProperty("InputWorkspace");
88 boost::optional<Mantid::Kernel::SpecialCoordinateSystem> coordinateSystem = converter(m_inputWS.get());
89 if (*coordinateSystem != Mantid::Kernel::SpecialCoordinateSystem::HKL) {
90 std::stringstream errmsg;
91 errmsg << "Input MDWorkspace's coordinate system is not HKL.";
92 throw std::invalid_argument(errmsg.str());
93 }
94
96 PeaksWorkspace_sptr inPeakWS = getProperty("PeaksWorkspace");
97 const double box = getProperty("DeltaHKL");
98 const int gridPts = getProperty("GridPoints");
99 const int neighborPts = getProperty("NeighborPoints");
101 PeaksWorkspace_sptr peakWS = getProperty("OutputWorkspace");
102 if (peakWS != inPeakWS)
103 peakWS = inPeakWS->clone();
104
105 MatrixWorkspace_sptr flux = getProperty("FluxWorkspace");
106 MatrixWorkspace_sptr sa = getProperty("SolidAngleWorkspace");
107
108 IMDEventWorkspace_sptr m_eventWS = std::dynamic_pointer_cast<IMDEventWorkspace>(m_inputWS);
109 IMDHistoWorkspace_sptr m_histoWS = std::dynamic_pointer_cast<IMDHistoWorkspace>(m_inputWS);
110 int npeaks = peakWS->getNumberPeaks();
111
112 auto prog = std::make_unique<Progress>(this, 0.3, 1.0, npeaks);
114 for (int i = 0; i < npeaks; i++) {
116
117 IPeak &p = peakWS->getPeak(i);
118 // round to integer
119 int h = static_cast<int>(boost::math::iround(p.getH()));
120 int k = static_cast<int>(boost::math::iround(p.getK()));
121 int l = static_cast<int>(boost::math::iround(p.getL()));
122 MDHistoWorkspace_sptr histoBox;
123 if (m_histoWS) {
124 histoBox = cropHisto(h, k, l, box, m_histoWS);
125 } else if (sa && flux) {
126 histoBox = normalize(h, k, l, box, gridPts, flux, sa, m_eventWS);
127 } else {
128 histoBox = binEvent(h, k, l, box, gridPts, m_eventWS);
129 }
130 double intensity = 0.0;
131 double errorSquared = 0.0;
132 integratePeak(neighborPts, histoBox, intensity, errorSquared);
133 p.setIntensity(intensity);
134 p.setSigmaIntensity(sqrt(errorSquared));
135 prog->report();
137 }
139 // Save the output
140 setProperty("OutputWorkspace", peakWS);
141}
142
143MDHistoWorkspace_sptr IntegratePeaksMDHKL::normalize(int h, int k, int l, double box, int gridPts,
144 const MatrixWorkspace_sptr &flux, const MatrixWorkspace_sptr &sa,
145 const IMDEventWorkspace_sptr &ws) {
146 auto normAlg = createChildAlgorithm("MDNormSCD");
147 normAlg->setProperty("InputWorkspace", ws);
148 normAlg->setProperty("AlignedDim0", "[H,0,0]," + boost::lexical_cast<std::string>(h - box) + "," +
149 boost::lexical_cast<std::string>(h + box) + "," + std::to_string(gridPts));
150 normAlg->setProperty("AlignedDim1", "[0,K,0]," + boost::lexical_cast<std::string>(k - box) + "," +
151 boost::lexical_cast<std::string>(k + box) + "," + std::to_string(gridPts));
152 normAlg->setProperty("AlignedDim2", "[0,0,L]," + boost::lexical_cast<std::string>(l - box) + "," +
153 boost::lexical_cast<std::string>(l + box) + "," + std::to_string(gridPts));
154 normAlg->setProperty("FluxWorkspace", flux);
155 normAlg->setProperty("SolidAngleWorkspace", sa);
156 normAlg->setProperty("OutputWorkspace", "mdout");
157 normAlg->setProperty("OutputNormalizationWorkspace", "mdnorm");
158 normAlg->executeAsChildAlg();
159 Workspace_sptr mdout = normAlg->getProperty("OutputWorkspace");
160 Workspace_sptr mdnorm = normAlg->getProperty("OutputNormalizationWorkspace");
161
162 auto alg = createChildAlgorithm("DivideMD");
163 alg->setProperty("LHSWorkspace", mdout);
164 alg->setProperty("RHSWorkspace", mdnorm);
165 alg->setPropertyValue("OutputWorkspace", "out");
166 alg->execute();
167 IMDWorkspace_sptr out = alg->getProperty("OutputWorkspace");
168 return std::dynamic_pointer_cast<MDHistoWorkspace>(out);
169}
170
171void IntegratePeaksMDHKL::integratePeak(const int neighborPts, const MDHistoWorkspace_sptr &out, double &intensity,
172 double &errorSquared) {
173 std::vector<int> gridPts;
175 double BackgroundOuterRadius2 = getProperty("BackgroundOuterRadius");
176 if (BackgroundOuterRadius2 != EMPTY_DBL())
177 BackgroundOuterRadius2 = pow(BackgroundOuterRadius2, 2.0);
179 double BackgroundInnerRadius2 = getProperty("BackgroundInnerRadius");
180 if (BackgroundInnerRadius2 != EMPTY_DBL())
181 BackgroundInnerRadius2 = pow(BackgroundInnerRadius2, 2.0);
182 const size_t dimensionality = out->getNumDims();
183 for (size_t i = 0; i < dimensionality; ++i) {
184 gridPts.emplace_back(static_cast<int>(out->getDimension(i)->getNBins()));
185 }
186
187 const auto F = out->getSignalArray();
188 double Fmax = 0.;
189 double Fmin = std::numeric_limits<double>::max();
190 for (int i = 0; i < gridPts[0] * gridPts[1] * gridPts[2]; i++) {
191 if (std::isnormal(F[i])) {
192 Fmin = std::min(Fmin, F[i]);
193 Fmax = std::max(Fmax, F[i]);
194 }
195 }
196 const auto SqError = out->getErrorSquaredArray();
197
198 double minIntensity = Fmin + 0.01 * (Fmax - Fmin);
199 int measuredPoints = 0;
200 int peakPoints = 0;
201 double peakSum = 0.0;
202 double measuredSum = 0.0;
203 double errSqSum = 0.0;
204 double measuredErrSqSum = 0.0;
205 int backgroundPoints = 0;
206 double backgroundSum = 0.0;
207 double backgroundErrSqSum = 0.0;
208 double Hcenter = gridPts[0] * 0.5;
209 double Kcenter = gridPts[1] * 0.5;
210 double Lcenter = gridPts[2] * 0.5;
211
212 for (int Hindex = 0; Hindex < gridPts[0]; Hindex++) {
213 for (int Kindex = 0; Kindex < gridPts[1]; Kindex++) {
214 for (int Lindex = 0; Lindex < gridPts[2]; Lindex++) {
215 int iHKL = Hindex + gridPts[0] * (Kindex + gridPts[1] * Lindex);
216 if (BackgroundOuterRadius2 != EMPTY_DBL()) {
217 double radius2 = pow((double(Hindex) - Hcenter) / gridPts[0], 2) +
218 pow((double(Kindex) - Kcenter) / gridPts[1], 2) +
219 pow((double(Lindex) - Lcenter) / gridPts[2], 2);
220 if (radius2 < BackgroundOuterRadius2 && BackgroundInnerRadius2 < radius2) {
221 backgroundPoints = backgroundPoints + 1;
222 backgroundSum = backgroundSum + F[iHKL];
223 backgroundErrSqSum = backgroundErrSqSum + SqError[iHKL];
224 }
225 }
226 if (std::isfinite(F[iHKL])) {
227 measuredPoints = measuredPoints + 1;
228 measuredSum = measuredSum + F[iHKL];
229 measuredErrSqSum = measuredErrSqSum + SqError[iHKL];
230 if (F[iHKL] > minIntensity) {
231 int neighborPoints = 0;
232 for (int Hj = -2; Hj < 3; Hj++) {
233 for (int Kj = -2; Kj < 3; Kj++) {
234 for (int Lj = -2; Lj < 3; Lj++) {
235 int jHKL = Hindex + Hj + gridPts[0] * (Kindex + Kj + gridPts[1] * (Lindex + Lj));
236 if (Lindex + Lj >= 0 && Lindex + Lj < gridPts[2] && Kindex + Kj >= 0 && Kindex + Kj < gridPts[1] &&
237 Hindex + Hj >= 0 && Hindex + Hj < gridPts[0] && F[jHKL] > minIntensity) {
238 neighborPoints = neighborPoints + 1;
239 }
240 }
241 }
242 }
243 if (neighborPoints >= neighborPts) {
244 peakPoints = peakPoints + 1;
245 peakSum = peakSum + F[iHKL];
246 errSqSum = errSqSum + SqError[iHKL];
247 }
248 }
249 } else {
250 double minR = sqrt(std::pow(float(Hindex) / float(gridPts[0]) - 0.5, 2) +
251 std::pow(float(Kindex) / float(gridPts[1]) - 0.5, 2) +
252 std::pow(float(Lindex) / float(gridPts[0]) - 0.5, 2));
253 if (minR < 0.05) {
254 intensity = 0.0;
255 errorSquared = 0.0;
256 return;
257 }
258 }
259 }
260 }
261 }
262 if (BackgroundOuterRadius2 != EMPTY_DBL()) {
263 double ratio = 0.0;
264 if (backgroundPoints > 0) {
265 ratio = float(peakPoints) / float(backgroundPoints);
266 }
267 intensity = peakSum - ratio * (backgroundSum);
268 errorSquared = errSqSum + ratio * ratio * (backgroundErrSqSum);
269 } else {
270 double ratio = float(peakPoints) / float(measuredPoints - peakPoints);
271 intensity = peakSum - ratio * (measuredSum - peakSum);
272 errorSquared = errSqSum + ratio * ratio * (measuredErrSqSum - errSqSum);
273 }
274 return;
275}
276
282MDHistoWorkspace_sptr IntegratePeaksMDHKL::binEvent(int h, int k, int l, double box, int gridPts,
283 const IMDWorkspace_sptr &ws) {
284 auto binMD = createChildAlgorithm("BinMD", 0.0, 0.3);
285 binMD->setProperty("InputWorkspace", ws);
286 binMD->setProperty("AlignedDim0", "[H,0,0]," + boost::lexical_cast<std::string>(h - box) + "," +
287 boost::lexical_cast<std::string>(h + box) + "," + std::to_string(gridPts));
288 binMD->setProperty("AlignedDim1", "[0,K,0]," + boost::lexical_cast<std::string>(k - box) + "," +
289 boost::lexical_cast<std::string>(k + box) + "," + std::to_string(gridPts));
290 binMD->setProperty("AlignedDim2", "[0,0,L]," + boost::lexical_cast<std::string>(l - box) + "," +
291 boost::lexical_cast<std::string>(l + box) + "," + std::to_string(gridPts));
292 binMD->setPropertyValue("AxisAligned", "1");
293 binMD->setPropertyValue("OutputWorkspace", "out");
294 binMD->executeAsChildAlg();
295 Workspace_sptr outputWS = binMD->getProperty("OutputWorkspace");
296 return std::dynamic_pointer_cast<MDHistoWorkspace>(outputWS);
297}
298
304MDHistoWorkspace_sptr IntegratePeaksMDHKL::cropHisto(int h, int k, int l, double box, const IMDWorkspace_sptr &ws) {
305 auto cropMD = createChildAlgorithm("IntegrateMDHistoWorkspace", 0.0, 0.3);
306 cropMD->setProperty("InputWorkspace", ws);
307
308 cropMD->setProperty("P1Bin",
309 boost::lexical_cast<std::string>(h - box) + ",0," + boost::lexical_cast<std::string>(h + box));
310 cropMD->setProperty("P2Bin",
311 boost::lexical_cast<std::string>(k - box) + ",0," + boost::lexical_cast<std::string>(k + box));
312 cropMD->setProperty("P3Bin",
313 boost::lexical_cast<std::string>(l - box) + ",0," + boost::lexical_cast<std::string>(l + box));
314
315 cropMD->setPropertyValue("OutputWorkspace", "out");
316 cropMD->executeAsChildAlg();
317 IMDHistoWorkspace_sptr outputWS = cropMD->getProperty("OutputWorkspace");
318 return std::dynamic_pointer_cast<MDHistoWorkspace>(outputWS);
319}
320
321} // namespace Mantid::MDAlgorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
#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.
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
A validator which provides a TENTATIVE check that a workspace contains common bins in each spectrum.
Kernel::IValidator_sptr clone() const override
Clone the current state.
A validator which checks that a workspace has a valid instrument.
A property class for workspaces.
A validator which checks that the unit of the workspace referred to by a WorkspaceProperty is the exp...
MDFrameFromMDWorkspace: Each dimension of the MDWorkspace contains an MDFrame.
Structure describing a single-crystal peak.
Definition: IPeak.h:26
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
The concrete, templated class for properties.
Integrate single-crystal peaks in reciprocal-space.
DataObjects::MDHistoWorkspace_sptr normalize(int h, int k, int l, double box, int gridPts, const API::MatrixWorkspace_sptr &flux, const API::MatrixWorkspace_sptr &sa, const API::IMDEventWorkspace_sptr &ws)
DataObjects::MDHistoWorkspace_sptr cropHisto(int h, int k, int l, double box, const API::IMDWorkspace_sptr &ws)
Runs the BinMD algorithm on the input to provide the output workspace All slicing algorithm propertie...
DataObjects::MDHistoWorkspace_sptr binEvent(int h, int k, int l, double box, int gridPts, const API::IMDWorkspace_sptr &ws)
Runs the BinMD algorithm on the input to provide the output workspace All slicing algorithm propertie...
void integratePeak(const int neighborPts, const DataObjects::MDHistoWorkspace_sptr &out, double &intensity, double &errorSquared)
void exec() override
Run the algorithm.
std::shared_ptr< IMDEventWorkspace > IMDEventWorkspace_sptr
Shared pointer to Mantid::API::IMDEventWorkspace.
std::shared_ptr< Workspace > Workspace_sptr
shared pointer to Mantid::API::Workspace
Definition: Workspace_fwd.h:20
std::shared_ptr< IMDHistoWorkspace > IMDHistoWorkspace_sptr
shared pointer to Mantid::API::IMDHistoWorkspace
std::shared_ptr< IMDWorkspace > IMDWorkspace_sptr
Shared pointer to the IMDWorkspace base class.
Definition: IMDWorkspace.h:146
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
std::shared_ptr< MDHistoWorkspace > MDHistoWorkspace_sptr
A shared pointer to a MDHistoWorkspace.
std::shared_ptr< PeaksWorkspace > PeaksWorkspace_sptr
Typedef for a shared pointer to a peaks workspace.
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
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
Definition: EmptyValues.h:43
std::string to_string(const wide_integer< Bits, Signed > &n)
Describes the direction (within an algorithm) of a Property.
Definition: Property.h:50
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54