39 return "Recalculates trajectory limits set by CropWorkspaceForMDNorm";
47 "An input MDEventWorkspace. Must be in Q_sample frame.");
49 "Copy of the input MDEventWorkspace with the corrected "
50 "trajectory extents.");
55 std::map<std::string, std::string> errorMessage;
59 if (inputWS->getNumDims() < 3) {
60 errorMessage.emplace(
"InputWorkspace",
"The input workspace must be at least 3D");
62 for (
size_t i = 0; i < 3; i++) {
64 errorMessage.emplace(
"InputWorkspace",
"The input workspace must be in Q_sample");
69 size_t nExperimentInfos = inputWS->getNumExperimentInfo();
70 if (nExperimentInfos == 0) {
71 errorMessage.emplace(
"InputWorkspace",
"There must be at least one experiment info");
73 for (
size_t iExpInfo = 0; iExpInfo < nExperimentInfos; iExpInfo++) {
74 auto ¤tExptInfo = *(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");
80 if (!currentExptInfo.run().hasProperty(
"MDNorm_high")) {
81 errorMessage.emplace(
"InputWorkspace",
"Missing MDNorm_high log. Please use "
82 "CropWorkspaceForMDNorm before converting to MD");
99 outWS = inWS->clone();
103 bool diffraction(
true);
105 if (outWS->getNumDims() > 3) {
106 if (outWS->getDimension(3)->getMDFrame().name() ==
"DeltaE") {
108 if (outWS->getExperimentInfo(0)->run().hasProperty(
"Ei")) {
109 Ei = outWS->getExperimentInfo(0)->run().getPropertyValueAsType<
double>(
"Ei");
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");
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);
133 dEmin = outWS->getDimension(3)->getMinimum();
134 dEmax = outWS->getDimension(3)->getMaximum();
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());
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");
151 for (
size_t iExpInfo = 0; iExpInfo < nExperimentInfos; iExpInfo++) {
152 auto ¤tExptInfo = *(outWS->getExperimentInfo(
static_cast<uint16_t
>(iExpInfo)));
154 const auto &spectrumInfo = currentExptInfo.spectrumInfo();
155 const auto nspectra =
static_cast<int64_t
>(spectrumInfo.size());
156 std::vector<double> lowValues, highValues;
158 auto *lowValuesLog =
dynamic_cast<VectorDoubleProperty *
>(currentExptInfo.getLog(
"MDNorm_low"));
159 lowValues = (*lowValuesLog)();
161 auto *highValuesLog =
dynamic_cast<VectorDoubleProperty *
>(currentExptInfo.getLog(
"MDNorm_high"));
162 highValues = (*highValuesLog)();
165 bool zeroWeights(
false);
166 for (
size_t iOtherDims = 0; iOtherDims < otherDimsNames.size(); iOtherDims++) {
170 currentExptInfo.run().getProperty(otherDimsNames[iOtherDims]));
171 if ((otherDimsLog->firstValue() < otherDimsMin[iOtherDims]) ||
172 (otherDimsLog->firstValue() > otherDimsMax[iOtherDims])) {
174 g_log.
warning() <<
"In experimentInfo " << iExpInfo <<
", log " << otherDimsNames[iOtherDims]
175 <<
" is outside limits\n";
180 highValues = lowValues;
182 auto source = currentExptInfo.getInstrument()->getSource()->getPos();
183 auto sample = currentExptInfo.getInstrument()->getSample()->getPos();
184 auto beamDir = sample - source;
186 auto gon = currentExptInfo.run().getGoniometerMatrix();
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];
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;
201 if (convention ==
"Crystallography") {
207 kfmin = lowValues[i];
208 kfmax = highValues[i];
209 qLabLow = (qin - qout) * kfmin;
210 qLabHigh = (qin - qout) * kfmax;
212 if (dEmin > lowValues[i]) {
213 lowValues[i] = dEmin;
215 if (dEmax < highValues[i]) {
216 highValues[i] = dEmax;
218 if (lowValues[i] > Ei) {
221 if (highValues[i] > Ei) {
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;
231 qSampleLow = gon * qLabLow;
232 qSampleHigh = gon * qLabHigh;
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];
245 if ((qxmin - qSampleLow.
X()) * (qxmin - qSampleHigh.
X()) < 0) {
246 double kfIntersection =
247 (qxmin - qSampleLow.
X()) * (kfmax - kfmin) / (qSampleHigh.
X() - qSampleLow.
X()) + kfmin;
249 kfIntersection = Ei - kfIntersection * kfIntersection / energyToK;
251 if ((qSampleLow.
X() < qxmin) && (lowValues[i] < kfIntersection)) {
252 lowValues[i] = kfIntersection;
254 if ((qSampleHigh.
X() < qxmin) && (highValues[i] > kfIntersection)) {
255 highValues[i] = kfIntersection;
259 if ((qxmax - qSampleLow.
X()) * (qxmax - qSampleHigh.
X()) < 0) {
260 double kfIntersection =
261 (qxmax - qSampleLow.
X()) * (kfmax - kfmin) / (qSampleHigh.
X() - qSampleLow.
X()) + kfmin;
263 kfIntersection = Ei - kfIntersection * kfIntersection / energyToK;
265 if ((qSampleLow.
X() > qxmax) && (lowValues[i] < kfIntersection)) {
266 lowValues[i] = kfIntersection;
268 if ((qSampleHigh.
X() > qxmax) && (highValues[i] > kfIntersection)) {
269 highValues[i] = kfIntersection;
273 if ((qymin - qSampleLow.
Y()) * (qymin - qSampleHigh.
Y()) < 0) {
274 double kfIntersection =
275 (qymin - qSampleLow.
Y()) * (kfmax - kfmin) / (qSampleHigh.
Y() - qSampleLow.
Y()) + kfmin;
277 kfIntersection = Ei - kfIntersection * kfIntersection / energyToK;
279 if ((qSampleLow.
Y() < qymin) && (lowValues[i] < kfIntersection)) {
280 lowValues[i] = kfIntersection;
282 if ((qSampleHigh.
Y() < qymin) && (highValues[i] > kfIntersection)) {
283 highValues[i] = kfIntersection;
287 if ((qymax - qSampleLow.
Y()) * (qymax - qSampleHigh.
Y()) < 0) {
288 double kfIntersection =
289 (qymax - qSampleLow.
Y()) * (kfmax - kfmin) / (qSampleHigh.
Y() - qSampleLow.
Y()) + kfmin;
291 kfIntersection = Ei - kfIntersection * kfIntersection / energyToK;
293 if ((qSampleLow.
Y() > qymax) && (lowValues[i] < kfIntersection)) {
294 lowValues[i] = kfIntersection;
296 if ((qSampleHigh.
Y() > qymax) && (highValues[i] > kfIntersection)) {
297 highValues[i] = kfIntersection;
301 if ((qzmin - qSampleLow.
Z()) * (qzmin - qSampleHigh.
Z()) < 0) {
302 double kfIntersection =
303 (qzmin - qSampleLow.
Z()) * (kfmax - kfmin) / (qSampleHigh.
Z() - qSampleLow.
Z()) + kfmin;
305 kfIntersection = Ei - kfIntersection * kfIntersection / energyToK;
307 if ((qSampleLow.
Z() < qzmin) && (lowValues[i] < kfIntersection)) {
308 lowValues[i] = kfIntersection;
310 if ((qSampleHigh.
Z() < qzmin) && (highValues[i] > kfIntersection)) {
311 highValues[i] = kfIntersection;
315 if ((qzmax - qSampleLow.
Z()) * (qzmax - qSampleHigh.
Z()) < 0) {
316 double kfIntersection =
317 (qzmax - qSampleLow.
Z()) * (kfmax - kfmin) / (qSampleHigh.
Z() - qSampleLow.
Z()) + kfmin;
319 kfIntersection = Ei - kfIntersection * kfIntersection / energyToK;
321 if ((qSampleLow.
Z() > qzmax) && (lowValues[i] < kfIntersection)) {
322 lowValues[i] = kfIntersection;
324 if ((qSampleHigh.
Z() > qzmax) && (highValues[i] > kfIntersection)) {
325 highValues[i] = kfIntersection;
332 currentExptInfo.mutableRun().addProperty(
"MDNorm_low", lowValues,
true);
333 currentExptInfo.mutableRun().addProperty(
"MDNorm_high", highValues,
true);
#define DECLARE_ALGORITHM(classname)
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.
A property class for workspaces.
static const std::string QSampleName
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.
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.
constexpr double X() const noexcept
Get x.
constexpr double Y() const noexcept
Get y.
constexpr double Z() const noexcept
Get z.
RecalculateTrajectoriesExtents :
void init() override
Initialize the algorithm's properties.
std::map< std::string, std::string > validateInputs() override final
Validate the input workspace.
void exec() override
Execute the algorithm.
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.
@ Input
An input workspace.
@ Output
An output workspace.