Mantid
Loading...
Searching...
No Matches
DgsScatteredTransmissionCorrectionMD.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 +
7
8// local
10
11// 3rd party
18
19// standard
20#include <limits>
21#include <math.h>
22
23using namespace Mantid::API;
24using namespace Mantid::DataObjects;
25using namespace Mantid::Geometry;
26using namespace Mantid::Kernel;
27
28namespace Mantid::MDAlgorithms {
29
30constexpr double EMPTY_FLT() noexcept { return std::numeric_limits<float>::max() / 2; }
31
32// Register the algorithm into the AlgorithmFactory
33DECLARE_ALGORITHM(DgsScatteredTransmissionCorrectionMD)
34
35//---------------------------------------------------------------------------------------------------------
36
38 declareProperty(std::make_unique<WorkspaceProperty<IMDEventWorkspace>>("InputWorkspace", "", Direction::Input),
39 "Input MDEventWorkspace. Either QSample (or QLab) frame plus DeltaE, or just Qmod plus DeltaE");
40
41 auto mustBePositive = std::make_shared<BoundedValidator<double>>();
42 mustBePositive->setLower(0.0);
43 mustBePositive->setLowerExclusive(true);
44 declareProperty(
45 std::make_unique<PropertyWithValue<double>>("ExponentFactor", EMPTY_DBL(), mustBePositive, Direction::Input),
46 "Depletion rate exponent");
47
48 declareProperty(std::make_unique<WorkspaceProperty<IMDEventWorkspace>>("OutputWorkspace", "", Direction::Output),
49 "The output MDEventWorkspace with the correction applied");
50}
51
52//---------------------------------------------------------------------------------------------------------
53
54std::map<std::string, std::string> DgsScatteredTransmissionCorrectionMD::validateInputs() {
55 std::map<std::string, std::string> output;
56 // validate input workspace
57 std::string workspace_error = checkInputWorkspace();
58 if (workspace_error != "")
59 output["InputWorkspace"] = workspace_error;
60 // ExponentFactor should be positive
61 double c = getProperty("ExponentFactor");
62 if (c <= 0.0)
63 output["ExponentFactor"] = "ExponentFactor should be positive";
64 return output;
65}
66
67//---------------------------------------------------------------------------------------------------------
68
70
71 // Verify the dimension of the input workspace
72 std::string dimensionError = checkInputMDDimensions();
73 if (dimensionError != "")
74 return dimensionError;
75
76 // Verify input workspace has an Efixed metadata
77 IMDEventWorkspace_sptr inputws = getProperty("InputWorkspace");
78 uint16_t nRuns(inputws->getNumExperimentInfo());
79 mEfixedValues.resize(nRuns);
80 for (uint16_t i = 0; i < inputws->getNumExperimentInfo(); i++) {
81 const auto expInfo = inputws->getExperimentInfo(i);
82 try {
83 mEfixedValues[i] = static_cast<float>(expInfo->getEFixed());
84 } catch (std::runtime_error &e) {
85 return e.what();
86 }
87 }
88 return ""; // all is well
89}
90
91//---------------------------------------------------------------------------------------------------------
92
94 std::string errormsg("");
95 IMDEventWorkspace_sptr inputws = getProperty("InputWorkspace");
96 size_t numdims = inputws->getNumDims();
97
98 // Get and check the dimensions: Q3D or Q1D
99 const SpecialCoordinateSystem coordsys = inputws->getSpecialCoordinateSystem();
100 size_t qdim(0);
101 std::string qdimstr("Not Q3D or |Q|");
102 if (coordsys == SpecialCoordinateSystem::QLab || coordsys == SpecialCoordinateSystem::QSample) {
103 // q3d
104 qdim = 3;
105 qdimstr = "Q3D";
106 } else {
107 // search Q1D: at any place
108 for (size_t i = 0; i < numdims; ++i) {
109 if (inputws->getDimension(i)->getName() == "|Q|") {
110 qdim = 1;
111 qdimstr = "|Q|";
112 break;
113 }
114 }
115 }
116
117 // Check dimension index for dimension "DeltaE"
118 bool qModCase = (qdim == 1) && (inputws->getDimension(1)->getName() == "DeltaE");
119 bool qVecCase = (qdim == 3) && (inputws->getDimension(3)->getName() == "DeltaE");
120 if (!qModCase && !qVecCase) { // both allowed q-cases failed
121 g_log.error() << "Coordinate system = " << coordsys << " does not meet requirement: \n";
122 for (size_t i = 0; i < numdims; ++i) {
123 g_log.error() << i << "-th dim: " << inputws->getDimension(i)->getName() << "\n";
124 }
125 errormsg += "Q Dimension (" + qdimstr +
126 ") is neither Q3D nor |Q|. Or DeltaE is found in an improper place (2nd or 4th dimension).";
127 }
128
129 return errormsg;
130}
131
132//---------------------------------------------------------------------------------------------------------
133
134template <typename MDE, size_t nd>
136 double cDouble = getProperty("ExponentFactor");
137 float c = static_cast<float>(cDouble);
138 size_t deltaEIndex = ws->getNumDims() - 1;
139
140 // Get Box from MDEventWorkspace
141 MDBoxBase<MDE, nd> *box1 = ws->getBox();
142 std::vector<API::IMDNode *> boxes;
143 box1->getBoxes(boxes, 1000, true);
144 auto numBoxes = int(boxes.size());
145
146 // Add the boxes in parallel. They should be spread out enough on each core to avoid stepping on each other.
147
148 PRAGMA_OMP( parallel for if (!ws->isFileBacked()))
149 for (int i = 0; i < numBoxes; ++i) {
151 auto *box = dynamic_cast<MDBox<MDE, nd> *>(boxes[i]);
152 if (box && !box->getIsMasked()) {
153 std::vector<MDE> &events = box->getEvents();
154 for (auto it = events.begin(); it != events.end(); ++it) {
155 float Ei = mEfixedValues[it->getExpInfoIndex()];
156 float Ef = Ei - it->getCenter(deltaEIndex); // Ei - deltaE
157 float correction(exp(c * Ef));
158
159 it->setSignal(it->getSignal() * correction);
160 it->setErrorSquared(it->getErrorSquared() * correction * correction); // no error from `correction`
161 }
162 }
163 if (box) {
164 box->releaseEvents();
165 }
167 }
169}
170
171//---------------------------------------------------------------------------------------------------------
172
174 // Get input workspace
175 IMDEventWorkspace_sptr inputWs = getProperty("InputWorkspace");
176
177 // Initialize the output workspace
178 std::string outputWsName = getPropertyValue("OutputWorkspace");
179 IMDEventWorkspace_sptr outputWs(nullptr);
180 if (inputWs->getName() == outputWsName)
181 outputWs = inputWs;
182 else
183 outputWs = inputWs->clone();
184
185 // Apply detailed balance to MDEvents
187
188 // refresh cache for MDBoxes: set correct Box signal
189 outputWs->refreshCache();
190
191 // Clear masking (box flags) from the output workspace
192 outputWs->clearMDMasking();
193
194 // Set output
195 setProperty("OutputWorkspace", outputWs);
196}
197
198} // namespace Mantid::MDAlgorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
#define CALL_MDEVENT_FUNCTION(funcname, workspace)
Macro that makes it possible to call a templated method for a MDEventWorkspace using a IMDEventWorksp...
#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 PRAGMA_OMP(expression)
#define PARALLEL_CHECK_INTERRUPT_REGION
Adds a check after a Parallel region to see if it was interupted.
std::string getPropertyValue(const std::string &name) const override
Get the value of a property as a string.
Definition: Algorithm.cpp:2026
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
virtual void getBoxes(std::vector< IMDNode * > &boxes, size_t maxDepth, bool leafOnly)=0
Fill a vector with all the boxes who are the childred of this one up to a certain depth.
A property class for workspaces.
Templated super-class of a multi-dimensional event "box".
Definition: MDBoxBase.h:50
Templated class for a multi-dimensional event "box".
Definition: MDBox.h:45
std::vector< MDE > & getEvents()
Get vector of events to change.
std::shared_ptr< MDEventWorkspace< MDE, nd > > sptr
Typedef for a shared pointer of this kind of event workspace.
size_t getNumDims() const override
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.
Definition: Logger.cpp:77
The concrete, templated class for properties.
std::map< std::string, std::string > validateInputs() override
Verify the input properties meets certain requirements.
std::string checkInputMDDimensions()
Verify the input Workspace dimensions are either QSample (or QLab) frame plus DeltaE,...
void correctForTransmission(typename MDEventWorkspace< MDE, nd >::sptr ws)
Apply the transmission correction to each event.
std::string checkInputWorkspace()
Verify the input workspace meets certain requirements.
std::shared_ptr< IMDEventWorkspace > IMDEventWorkspace_sptr
Shared pointer to Mantid::API::IMDEventWorkspace.
SpecialCoordinateSystem
Special coordinate systems for Q3D.
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
Definition: EmptyValues.h:43
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54