Mantid
Loading...
Searching...
No Matches
ApplyDetailedBalanceMD.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 +
10// #include "MantidAPI/IMDEventWorkspace.h"
14#include "MantidAPI/Run.h"
22#include <math.h>
23
24using namespace Mantid::Kernel;
25using namespace Mantid::API;
26using namespace Mantid::Geometry;
27using namespace Mantid::DataObjects;
28
29namespace Mantid::MDAlgorithms {
30
31// Register the algorithm into the AlgorithmFactory
32DECLARE_ALGORITHM(ApplyDetailedBalanceMD)
33
34
35const std::string ApplyDetailedBalanceMD::name() const { return "ApplyDetailedBalanceMD"; }
36
38int ApplyDetailedBalanceMD::version() const { return 1; }
39
41const std::string ApplyDetailedBalanceMD::summary() const { return "Apply detailed balance to MDEventWorkspace"; }
42
44const std::string ApplyDetailedBalanceMD::category() const { return "MDAlgorithms"; }
45
46//---------------------------------------------------------------------------------------------------------
51
53 std::make_unique<WorkspaceProperty<API::IMDEventWorkspace>>("InputWorkspace", "", Kernel::Direction::Input),
54 "An input MDEventWorkspace. Must be in Q_sample/Q_lab frame. Must have an axis as DeltaE");
55
56 declareProperty(std::make_unique<PropertyWithValue<std::string>>("Temperature", "", Direction::Input),
57 "SampleLog variable name that contains the temperature or a number");
59 std::make_unique<WorkspaceProperty<API::IMDEventWorkspace>>("OutputWorkspace", "", Kernel::Direction::Output),
60 "The output MDEventWorkspace with detailed balance applied");
61}
62
63//---------------------------------------------------------------------------------------------------------
68 // Get input workspace
69 API::IMDEventWorkspace_sptr input_ws = getProperty("InputWorkspace");
70
71 // Process input workspace and create output workspace
72 std::string output_ws_name = getPropertyValue("OutputWorkspace");
73
74 API::IMDEventWorkspace_sptr output_ws(0);
75 if (input_ws->getName() == output_ws_name) {
76 // Calcualte in-place
77 output_ws = input_ws;
78 } else {
79 // Clone input workace to output workspace
80 output_ws = input_ws->clone();
81 }
82
83 // Apply detailed balance to MDEvents
85
86 // refresh cache for MDBoxes: set correct Box signal
87 output_ws->refreshCache();
88
89 // Clear masking (box flags) from the output workspace
90 output_ws->clearMDMasking();
91
92 // Set output
93 setProperty("OutputWorkspace", output_ws);
94}
95
96//---------------------------------------------------------------------------------------------------------
107std::map<std::string, std::string> ApplyDetailedBalanceMD::validateInputs() {
108 std::map<std::string, std::string> output;
109
110 // Get input workspace
111 API::IMDEventWorkspace_sptr input_ws = getProperty("InputWorkspace");
112
113 // check input dimension
114 std::string dim_error = checkInputMDDimension();
115 if (dim_error != "") {
116 output["InputWorkspace"] = dim_error;
117 }
118
119 std::string kerror = getTemperature(input_ws);
120 if (kerror != "")
121 output["Temperature"] = kerror;
122
123 return output;
124}
125
126//---------------------------------------------------------------------------------------------------------
134 std::string errormsg("");
135
136 API::IMDEventWorkspace_sptr inputws = getProperty("InputWorkspace");
137 size_t numdims = inputws->getNumDims();
138
139 // Get and check the dimensions: Q3D or Q1D
140 const Mantid::Kernel::SpecialCoordinateSystem coordsys = inputws->getSpecialCoordinateSystem();
141 size_t qdim(0);
142 std::string qdimstr("Not Q3D or |Q|");
145 // q3d
146 qdim = 3;
147 qdimstr = "Q3D";
148 } else {
149 // search Q1D: at any place
150 for (size_t i = 0; i < numdims; ++i) {
151 if (inputws->getDimension(i)->getName() == "|Q|") {
152 qdim = 1;
153 qdimstr = "|Q|";
154 break;
155 }
156 }
157 }
158
159 // Check DeltaE
160 if (qdim == 1 && inputws->getDimension(1)->getName() == "DeltaE") {
161 // 2nd dimension is DeltaE
162 mDeltaEIndex = 1;
163 } else if (qdim == 3 && inputws->getDimension(3)->getName() == "DeltaE") {
164 // 4th dimension is DeltaE
165 mDeltaEIndex = 3;
166 } else {
167 // Error
168 g_log.error() << "Coordiate system = " << coordsys << " does not meet requirement: \n";
169 for (size_t i = 0; i < numdims; ++i) {
170 g_log.error() << i << "-th dim: " << inputws->getDimension(i)->getName() << "\n";
171 }
172 errormsg += "Q Dimension (" + qdimstr +
173 ") is neither Q3D nor |Q|. Or DeltaE is found in proper place (2nd or 4th dimension).";
174 }
175
176 return errormsg;
177}
178
179//---------------------------------------------------------------------------------------------
183template <typename MDE, size_t nd>
185 // Get Box from MDEventWorkspace
186 MDBoxBase<MDE, nd> *box1 = ws->getBox();
187 std::vector<API::IMDNode *> boxes;
188 box1->getBoxes(boxes, 1000, true);
189 auto numBoxes = int(boxes.size());
190
191 // Add the boxes in parallel. They should be spread out enough on each
192 // core to avoid stepping on each other.
193
194 PRAGMA_OMP( parallel for if (!ws->isFileBacked()))
195 for (int i = 0; i < numBoxes; ++i) {
197 auto *box = dynamic_cast<MDBox<MDE, nd> *>(boxes[i]);
198 if (box && !box->getIsMasked()) {
199 // get the MEEvents from box
200 std::vector<MDE> &events = box->getEvents();
201 // Add events, with bounds checking
202 for (auto it = events.begin(); it != events.end(); ++it) {
203 // Create the event
204 // do calculattion
205 float temperatue(static_cast<float>(mExpinfoTemperatureMean[it->getExpInfoIndex()]));
206
207 // delta_e = it->getCenter(mDeltaEIndex);
208 // factor = pi * (1 - exp(-deltaE/(kb*T)))
209 float factor = static_cast<float>(M_PI) *
210 (static_cast<float>(1.) - exp(-it->getCenter(mDeltaEIndex) *
211 static_cast<float>(PhysicalConstants::meVtoKelvin / temperatue)));
212
213 // calcalate and set intesity
214 auto intensity = it->getSignal() * factor;
215 it->setSignal(intensity);
216
217 // calculate and set error
218 auto error2 = it->getErrorSquared() * factor * factor;
219 // error2 *= factor * factor;
220 it->setErrorSquared(error2);
221 }
222 }
223 if (box) {
224 box->releaseEvents();
225 }
227 }
229
230 return;
231}
232
233//---------------------------------------------------------------------------------------------------------
240 // Get temperture sample log name
241 std::string Tstring = getProperty("Temperature");
242 std::string temperature_error("");
243
244 // Try to convert Tstring to a float
245 float temperature;
246 try {
247 temperature = boost::lexical_cast<float>(Tstring);
248 } catch (...) {
249 // set to a unphysical value
250 temperature = -10;
251 }
252
253 // the input property could be a valid float; if not must search the experiment info
255 uint16_t numexpinfo = mdws->getNumExperimentInfo();
256
257 for (uint16_t i = 0; i < numexpinfo; ++i) {
258 if (temperature < 0) {
259 // if user specified is not a valid float
260 ExperimentInfo_const_sptr expinfo = mdws->getExperimentInfo(i);
261 if (expinfo->run().hasProperty(Tstring)) {
262 auto log = dynamic_cast<Kernel::TimeSeriesProperty<double> *>(expinfo->run().getProperty(Tstring));
263 if (!log) {
264 // wrong type of sample log: must be TimeSeriesProperty<double>
265 std::stringstream errss;
266 errss << "ExperimentInfo" << i << " has " << Tstring << ", which is not a valid double-valuesd log";
267 temperature_error += errss.str() + "\n";
268 } else {
269 mExpinfoTemperatureMean[i] = log->getStatistics().mean;
270 }
271 } else {
272 // specified sample log does not exist
273 std::stringstream errss;
274 errss << "ExperimentInfo " << i << " does not have tempertaure log " << Tstring;
275 temperature_error += errss.str() + "\n";
276 }
277 } else {
278 // set user specified temperature to map
279 mExpinfoTemperatureMean[i] = temperature;
280 }
281 }
282
283 return temperature_error;
284}
285
286} // 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.
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
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.
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.
A specialised Property class for holding a series of time-value pairs.
ApplyDetailedBalance : Perform the And boolean operation on two MDHistoWorkspaces.
std::string checkInputMDDimension()
Check input workspace dimension.
std::string getTemperature(const API::IMDEventWorkspace_sptr &mdws)
Get temperature.
const std::string summary() const override
Summary.
void applyDetailedBalance(typename Mantid::DataObjects::MDEventWorkspace< MDE, nd >::sptr ws)
Apply detailed balance to each MDEvent.
void init() override
Initialize the proeprties.
std::map< uint16_t, double > mExpinfoTemperatureMean
map of temperature retrieved from sample logs
size_t mDeltaEIndex
index of the MD dimension index for DeltaE
const std::string category() const override
category
int version() const override
Algorithm's version for identification.
std::map< std::string, std::string > validateInputs() override
Validate inputs.
std::shared_ptr< IMDEventWorkspace > IMDEventWorkspace_sptr
Shared pointer to Mantid::API::IMDEventWorkspace.
std::shared_ptr< const ExperimentInfo > ExperimentInfo_const_sptr
Shared pointer to const ExperimentInfo.
SpecialCoordinateSystem
Special coordinate systems for Q3D.
static constexpr double meVtoKelvin
1 meV in Kelvin.
STL namespace.
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54