Mantid
Loading...
Searching...
No Matches
RecalculateTrajectoriesExtents.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/Run.h"
17
18namespace Mantid::MDAlgorithms {
19using namespace Mantid::Kernel;
20using namespace Mantid::API;
22
23// Register the algorithm into the AlgorithmFactory
25
26//----------------------------------------------------------------------------------------------
27
28
29const std::string RecalculateTrajectoriesExtents::name() const { return "RecalculateTrajectoriesExtents"; }
30
33
35const std::string RecalculateTrajectoriesExtents::category() const { return "MDAlgorithms\\Normalisation"; }
36
38const std::string RecalculateTrajectoriesExtents::summary() const {
39 return "Recalculates trajectory limits set by CropWorkspaceForMDNorm";
40}
41
42//----------------------------------------------------------------------------------------------
46 declareProperty(std::make_unique<WorkspaceProperty<IMDEventWorkspace>>("InputWorkspace", "", Direction::Input),
47 "An input MDEventWorkspace. Must be in Q_sample frame.");
48 declareProperty(std::make_unique<WorkspaceProperty<IMDEventWorkspace>>("OutputWorkspace", "", Direction::Output),
49 "Copy of the input MDEventWorkspace with the corrected "
50 "trajectory extents.");
51}
52
54std::map<std::string, std::string> RecalculateTrajectoriesExtents::validateInputs() {
55 std::map<std::string, std::string> errorMessage;
56
57 // Check for input workspace frame
58 Mantid::API::IMDEventWorkspace_sptr inputWS = this->getProperty("InputWorkspace");
59 if (inputWS->getNumDims() < 3) {
60 errorMessage.emplace("InputWorkspace", "The input workspace must be at least 3D");
61 } else {
62 for (size_t i = 0; i < 3; i++) {
63 if (inputWS->getDimension(i)->getMDFrame().name() != Mantid::Geometry::QSample::QSampleName) {
64 errorMessage.emplace("InputWorkspace", "The input workspace must be in Q_sample");
65 }
66 }
67 }
68 // Check for property MDNorm_low and MDNorm_high
69 size_t nExperimentInfos = inputWS->getNumExperimentInfo();
70 if (nExperimentInfos == 0) {
71 errorMessage.emplace("InputWorkspace", "There must be at least one experiment info");
72 } else {
73 for (size_t iExpInfo = 0; iExpInfo < nExperimentInfos; iExpInfo++) {
74 auto &currentExptInfo = *(inputWS->getExperimentInfo(static_cast<uint16_t>(iExpInfo)));
75 if (!currentExptInfo.run().hasProperty("MDNorm_low")) {
76 errorMessage.emplace("InputWorkspace", "Missing MDNorm_low log. Please "
77 "use CropWorkspaceForMDNorm "
78 "before converting to MD");
79 }
80 if (!currentExptInfo.run().hasProperty("MDNorm_high")) {
81 errorMessage.emplace("InputWorkspace", "Missing MDNorm_high log. Please use "
82 "CropWorkspaceForMDNorm before converting to MD");
83 }
84 }
85 }
86 return errorMessage;
87}
88
89//----------------------------------------------------------------------------------------------
93 IMDEventWorkspace_sptr inWS = getProperty("InputWorkspace");
94 IMDEventWorkspace_sptr outWS = getProperty("OutputWorkspace");
95
96 // If input and output workspaces are not the same, create a new workspace for
97 // the output
98 if (outWS != inWS) {
99 outWS = inWS->clone();
100 }
101
102 // check if using diffraction or direct inelastic
103 bool diffraction(true);
104 double Ei(0.0);
105 if (outWS->getNumDims() > 3) {
106 if (outWS->getDimension(3)->getMDFrame().name() == "DeltaE") {
107 diffraction = false;
108 if (outWS->getExperimentInfo(0)->run().hasProperty("Ei")) {
109 Ei = outWS->getExperimentInfo(0)->run().getPropertyValueAsType<double>("Ei");
110 } else {
111 throw std::runtime_error("Workspace contains energy transfer axis, but no Ei."
112 "The MD normalization workflow is not implemented for "
113 "indirect geometry");
114 }
115 }
116 }
117
118 const double energyToK = 8.0 * M_PI * M_PI * PhysicalConstants::NeutronMass * PhysicalConstants::meV * 1e-20 /
120
121 auto convention = Kernel::ConfigService::Instance().getString("Q.convention");
122 // get limits for all dimensions
123 double qxmin = outWS->getDimension(0)->getMinimum();
124 double qxmax = outWS->getDimension(0)->getMaximum();
125 double qymin = outWS->getDimension(1)->getMinimum();
126 double qymax = outWS->getDimension(1)->getMaximum();
127 double qzmin = outWS->getDimension(2)->getMinimum();
128 double qzmax = outWS->getDimension(2)->getMaximum();
129 double dEmin(0.0), dEmax(0.0);
130 size_t nqedims = 3;
131 if (!diffraction) {
132 nqedims = 4;
133 dEmin = outWS->getDimension(3)->getMinimum();
134 dEmax = outWS->getDimension(3)->getMaximum();
135 }
136 std::vector<double> otherDimsMin, otherDimsMax;
137 std::vector<std::string> otherDimsNames;
138 for (size_t i = nqedims; i < outWS->getNumDims(); i++) {
139 otherDimsMin.emplace_back(outWS->getDimension(i)->getMinimum());
140 otherDimsMax.emplace_back(outWS->getDimension(i)->getMaximum());
141 otherDimsNames.emplace_back(outWS->getDimension(i)->getName());
142 }
143
144 // Loop over experiment infos
145 size_t nExperimentInfos = outWS->getNumExperimentInfo();
146 if (nExperimentInfos > 1) {
147 g_log.warning("More than one experiment info. On merged workspaces, the "
148 "limits recalculations might be wrong");
149 }
150
151 for (size_t iExpInfo = 0; iExpInfo < nExperimentInfos; iExpInfo++) {
152 auto &currentExptInfo = *(outWS->getExperimentInfo(static_cast<uint16_t>(iExpInfo)));
153
154 const auto &spectrumInfo = currentExptInfo.spectrumInfo();
155 const auto nspectra = static_cast<int64_t>(spectrumInfo.size());
156 std::vector<double> lowValues, highValues;
157
158 auto *lowValuesLog = dynamic_cast<VectorDoubleProperty *>(currentExptInfo.getLog("MDNorm_low"));
159 lowValues = (*lowValuesLog)();
160
161 auto *highValuesLog = dynamic_cast<VectorDoubleProperty *>(currentExptInfo.getLog("MDNorm_high"));
162 highValues = (*highValuesLog)();
163
164 // deal with other dimensions first
165 bool zeroWeights(false);
166 for (size_t iOtherDims = 0; iOtherDims < otherDimsNames.size(); iOtherDims++) {
167 // check other dimensions. there might be no events, but if the first log
168 // value is not within limits, the weight should be zero as well
169 auto *otherDimsLog = dynamic_cast<Kernel::TimeSeriesProperty<double> *>(
170 currentExptInfo.run().getProperty(otherDimsNames[iOtherDims]));
171 if ((otherDimsLog->firstValue() < otherDimsMin[iOtherDims]) ||
172 (otherDimsLog->firstValue() > otherDimsMax[iOtherDims])) {
173 zeroWeights = true;
174 g_log.warning() << "In experimentInfo " << iExpInfo << ", log " << otherDimsNames[iOtherDims]
175 << " is outside limits\n";
176 continue;
177 }
178 }
179 if (zeroWeights) {
180 highValues = lowValues;
181 } else {
182 auto source = currentExptInfo.getInstrument()->getSource()->getPos();
183 auto sample = currentExptInfo.getInstrument()->getSample()->getPos();
184 auto beamDir = sample - source;
185 beamDir.normalize();
186 auto gon = currentExptInfo.run().getGoniometerMatrix();
187 gon.Invert();
188
189 // calculate limits in Q_sample
190 for (int64_t i = 0; i < nspectra; i++) {
191 if (!spectrumInfo.hasDetectors(i) || spectrumInfo.isMonitor(i) || spectrumInfo.isMasked(i)) {
192 highValues[i] = lowValues[i];
193 continue;
194 }
195 const auto &detector = spectrumInfo.detector(i);
196 double theta = detector.getTwoTheta(sample, beamDir);
197 double phi = detector.getPhi();
198 V3D qout(sin(theta) * cos(phi), sin(theta) * sin(phi), cos(theta)), qin(0., 0., 1.), qLabLow, qLabHigh,
199 qSampleLow, qSampleHigh;
200 double kfmin, kfmax;
201 if (convention == "Crystallography") {
202 qout *= -1;
203 qin *= -1;
204 }
205 if (diffraction) {
206 // units of limits are momentum
207 kfmin = lowValues[i];
208 kfmax = highValues[i];
209 qLabLow = (qin - qout) * kfmin;
210 qLabHigh = (qin - qout) * kfmax;
211 } else {
212 if (dEmin > lowValues[i]) {
213 lowValues[i] = dEmin;
214 }
215 if (dEmax < highValues[i]) {
216 highValues[i] = dEmax;
217 }
218 if (lowValues[i] > Ei) {
219 lowValues[i] = Ei;
220 }
221 if (highValues[i] > Ei) {
222 highValues[i] = Ei;
223 }
224
225 double ki = std::sqrt(energyToK * Ei);
226 kfmin = std::sqrt(energyToK * (Ei - lowValues[i]));
227 kfmax = std::sqrt(energyToK * (Ei - highValues[i]));
228 qLabLow = qin * ki - qout * kfmin;
229 qLabHigh = qin * ki - qout * kfmax;
230 }
231 qSampleLow = gon * qLabLow;
232 qSampleHigh = gon * qLabHigh;
233 // check intersection with the box
234 // completely outside the box -> no weight
235 if (((qSampleLow.X() < qxmin) && (qSampleHigh.X() < qxmin)) ||
236 ((qSampleLow.X() > qxmax) && (qSampleHigh.X() > qxmax)) ||
237 ((qSampleLow.Y() < qymin) && (qSampleHigh.Y() < qymin)) ||
238 ((qSampleLow.Y() > qymax) && (qSampleHigh.Y() > qymax)) ||
239 ((qSampleLow.Z() < qzmin) && (qSampleHigh.Z() < qzmin)) ||
240 ((qSampleLow.Z() > qzmax) && (qSampleHigh.Z() > qzmax))) {
241 highValues[i] = lowValues[i];
242 continue;
243 }
244 // either intersection or completely indide the box
245 if ((qxmin - qSampleLow.X()) * (qxmin - qSampleHigh.X()) < 0) {
246 double kfIntersection =
247 (qxmin - qSampleLow.X()) * (kfmax - kfmin) / (qSampleHigh.X() - qSampleLow.X()) + kfmin;
248 if (!diffraction) {
249 kfIntersection = Ei - kfIntersection * kfIntersection / energyToK;
250 }
251 if ((qSampleLow.X() < qxmin) && (lowValues[i] < kfIntersection)) {
252 lowValues[i] = kfIntersection;
253 }
254 if ((qSampleHigh.X() < qxmin) && (highValues[i] > kfIntersection)) {
255 highValues[i] = kfIntersection;
256 }
257 }
258
259 if ((qxmax - qSampleLow.X()) * (qxmax - qSampleHigh.X()) < 0) {
260 double kfIntersection =
261 (qxmax - qSampleLow.X()) * (kfmax - kfmin) / (qSampleHigh.X() - qSampleLow.X()) + kfmin;
262 if (!diffraction) {
263 kfIntersection = Ei - kfIntersection * kfIntersection / energyToK;
264 }
265 if ((qSampleLow.X() > qxmax) && (lowValues[i] < kfIntersection)) {
266 lowValues[i] = kfIntersection;
267 }
268 if ((qSampleHigh.X() > qxmax) && (highValues[i] > kfIntersection)) {
269 highValues[i] = kfIntersection;
270 }
271 }
272
273 if ((qymin - qSampleLow.Y()) * (qymin - qSampleHigh.Y()) < 0) {
274 double kfIntersection =
275 (qymin - qSampleLow.Y()) * (kfmax - kfmin) / (qSampleHigh.Y() - qSampleLow.Y()) + kfmin;
276 if (!diffraction) {
277 kfIntersection = Ei - kfIntersection * kfIntersection / energyToK;
278 }
279 if ((qSampleLow.Y() < qymin) && (lowValues[i] < kfIntersection)) {
280 lowValues[i] = kfIntersection;
281 }
282 if ((qSampleHigh.Y() < qymin) && (highValues[i] > kfIntersection)) {
283 highValues[i] = kfIntersection;
284 }
285 }
286
287 if ((qymax - qSampleLow.Y()) * (qymax - qSampleHigh.Y()) < 0) {
288 double kfIntersection =
289 (qymax - qSampleLow.Y()) * (kfmax - kfmin) / (qSampleHigh.Y() - qSampleLow.Y()) + kfmin;
290 if (!diffraction) {
291 kfIntersection = Ei - kfIntersection * kfIntersection / energyToK;
292 }
293 if ((qSampleLow.Y() > qymax) && (lowValues[i] < kfIntersection)) {
294 lowValues[i] = kfIntersection;
295 }
296 if ((qSampleHigh.Y() > qymax) && (highValues[i] > kfIntersection)) {
297 highValues[i] = kfIntersection;
298 }
299 }
300
301 if ((qzmin - qSampleLow.Z()) * (qzmin - qSampleHigh.Z()) < 0) {
302 double kfIntersection =
303 (qzmin - qSampleLow.Z()) * (kfmax - kfmin) / (qSampleHigh.Z() - qSampleLow.Z()) + kfmin;
304 if (!diffraction) {
305 kfIntersection = Ei - kfIntersection * kfIntersection / energyToK;
306 }
307 if ((qSampleLow.Z() < qzmin) && (lowValues[i] < kfIntersection)) {
308 lowValues[i] = kfIntersection;
309 }
310 if ((qSampleHigh.Z() < qzmin) && (highValues[i] > kfIntersection)) {
311 highValues[i] = kfIntersection;
312 }
313 }
314
315 if ((qzmax - qSampleLow.Z()) * (qzmax - qSampleHigh.Z()) < 0) {
316 double kfIntersection =
317 (qzmax - qSampleLow.Z()) * (kfmax - kfmin) / (qSampleHigh.Z() - qSampleLow.Z()) + kfmin;
318 if (!diffraction) {
319 kfIntersection = Ei - kfIntersection * kfIntersection / energyToK;
320 }
321 if ((qSampleLow.Z() > qzmax) && (lowValues[i] < kfIntersection)) {
322 lowValues[i] = kfIntersection;
323 }
324 if ((qSampleHigh.Z() > qzmax) && (highValues[i] > kfIntersection)) {
325 highValues[i] = kfIntersection;
326 }
327 }
328 } // end loop over spectra
329 }
330
331 // set the new values for the MDNorm_low and MDNorm_high
332 currentExptInfo.mutableRun().addProperty("MDNorm_low", lowValues, true);
333 currentExptInfo.mutableRun().addProperty("MDNorm_high", highValues, true);
334 }
335
336 setProperty("OutputWorkspace", outWS);
337}
338
339} // namespace Mantid::MDAlgorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
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
Kernel::Logger & g_log
Definition: Algorithm.h:451
A property class for workspaces.
static const std::string QSampleName
Definition: QSample.h:23
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
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...
A specialised Property class for holding a series of time-value pairs.
Class for 3D vectors.
Definition: V3D.h:34
constexpr double X() const noexcept
Get x.
Definition: V3D.h:232
constexpr double Y() const noexcept
Get y.
Definition: V3D.h:233
constexpr double Z() const noexcept
Get z.
Definition: V3D.h:234
void init() override
Initialize the algorithm's properties.
std::map< std::string, std::string > validateInputs() override final
Validate the input workspace.
const std::string category() const override
Algorithm's category for identification.
int version() const override
Algorithm's version for identification.
const std::string summary() const override
Algorithm's summary for use in the GUI and help.
std::shared_ptr< IMDEventWorkspace > IMDEventWorkspace_sptr
Shared pointer to Mantid::API::IMDEventWorkspace.
static constexpr double NeutronMass
Mass of the neutron in kg.
static constexpr double h
Planck constant in J*s.
static constexpr double meV
1 meV in Joules.
STL namespace.
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54