Mantid
Loading...
Searching...
No Matches
GetSpiceDataRawCountsFromMD.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 +
8
9#include "MantidAPI/Axis.h"
13#include "MantidAPI/Run.h"
18
19#include <algorithm>
20#include <utility>
21
22namespace Mantid::MDAlgorithms {
23
24using namespace Mantid::HistogramData;
25using namespace Mantid::API;
26using namespace Mantid::Kernel;
27
28DECLARE_ALGORITHM(GetSpiceDataRawCountsFromMD)
29
30
34 declareProperty(std::make_unique<WorkspaceProperty<IMDEventWorkspace>>("InputWorkspace", "", Direction::Input),
35 "Name of the input data MDEventWorkspace from which the raw "
36 "values are retrieved.");
37
38 declareProperty(std::make_unique<WorkspaceProperty<IMDEventWorkspace>>("MonitorWorkspace", "", Direction::Input),
39 "Name of the input monitor MDEventWorkspace paired with "
40 "input data workspace. ");
41
42 declareProperty(std::make_unique<WorkspaceProperty<MatrixWorkspace>>("OutputWorkspace", "", Direction::Output),
43 "Name of the output MatrixWorkspace containing the raw data required.");
44
45 std::array<std::string, 3> vecmode = {{"Pt.", "Detector", "Sample Log"}};
46 auto modevalidator = std::make_shared<ListValidator<std::string>>(vecmode);
47 declareProperty("Mode", "Detector", modevalidator,
48 "Mode selector. (1) Pt.: get the raw detectors' signal of the "
49 "specified Pt.; "
50 "(2) Detector: get one detector's signals along all experiment points; "
51 "(3) Sample Log: get the value of a sample log along all experiment "
52 "points.");
53
54 declareProperty("XLabel", "",
55 "Label for the X-values of the output MatrixWorkspace. "
56 "For mode 'Pt.', it won't take value. The output of X-axis is always "
57 "2-theta. "
58 "For mode 'Detector', if it is left blank, "
59 "the default then will be 2-theta value of detector's position. "
60 "For mode 'Sample Log', the default value is 'Pt.'. "
61 "In the later 2 modes, XLabel can be any supported sample log's name.");
62
63 declareProperty("Pt", EMPTY_INT(),
64 "Experiment point number (i.e., run_number in MDEventWorkspace "
65 "of the detectors' counts to be exported. "
66 "It is used in mode 'Pt.' only. ");
67
68 declareProperty("DetectorID", EMPTY_INT(),
69 "Detector ID of the detector whose raw counts "
70 "will be exported for all experiment points."
71 "It is used in mode 'Detector' only. ");
72
73 declareProperty("SampleLogName", "",
74 "Name of the sample log whose value to be exported. "
75 "It is used in mode 'Sample Log' only. ");
76
77 declareProperty("NormalizeByMonitorCounts", true,
78 "If specified as true, all the output detectors' counts will "
79 "be normalized by their monitor counts. ");
80}
81
86 // Process input parameters
87 IMDEventWorkspace_sptr datamdws = getProperty("InputWorkspace");
88 IMDEventWorkspace_sptr monitormdws = getProperty("MonitorWorkspace");
89 std::string mode = getProperty("Mode");
90 bool donormalize = getProperty("NormalizeByMonitorCounts");
91 std::string xlabel = getProperty("XLabel");
92
93 // Branch out by mode
94 std::vector<double> vecX;
95 std::vector<double> vecY;
96 std::string ylabel;
97 if (mode == "Pt.") {
98 // export detector counts for one specific Pt./run number
99 int runnumber = getProperty("Pt");
100 if (isEmpty(runnumber))
101 throw std::runtime_error("For 'Pt.', value of 'Pt.' must be specified.");
102 exportDetCountsOfRun(datamdws, monitormdws, runnumber, vecX, vecY, xlabel, ylabel, donormalize);
103 } else if (mode == "Detector") {
104 int detid = getProperty("DetectorID");
105 if (isEmpty(detid))
106 throw std::runtime_error("For mode 'Detector', value of 'DetectorID' must be specified.");
107 exportIndividualDetCounts(datamdws, monitormdws, detid, vecX, vecY, xlabel, ylabel, donormalize);
108 } else if (mode == "Sample Log") {
109 std::string samplelogname = getProperty("SampleLogName");
110 if (samplelogname.empty())
111 throw std::runtime_error("For mode 'Sample Log', value of 'SampleLogName' must be specified.");
112 exportSampleLogValue(datamdws, samplelogname, vecX, vecY, xlabel, ylabel);
113 } else {
114 // Raise exception
115 std::stringstream ess;
116 ess << "Mode " << mode << " is not supported.";
117 throw std::runtime_error(ess.str());
118 }
119
120 // Create and export output workspace
121 MatrixWorkspace_sptr outws = createOutputWorkspace(vecX, vecY, xlabel, ylabel);
122 setProperty("OutputWorkspace", outws);
123}
124
125//----------------------------------------------------------------------------------------------
138 const API::IMDEventWorkspace_const_sptr &monitormdws,
139 const int runnumber, std::vector<double> &vecX,
140 std::vector<double> &vecY, std::string &xlabel,
141 std::string &ylabel, bool donormalize) {
142 // Get detector counts
143 std::vector<double> vec2theta;
144 std::vector<double> vecDetCounts;
145 int detid = -1;
146 getDetCounts(datamdws, runnumber, detid, vec2theta, vecDetCounts, true);
147 if (vec2theta.size() != vecDetCounts.size())
148 throw std::runtime_error("Logic error! Vector of 2theta must have same size as "
149 "vector of detectors' counts.");
150
151 // Get monitor counts
152 std::vector<double> vec2thetaMon;
153 std::vector<double> vecMonitorCounts;
154 // Normalize if required
155 if (donormalize) {
156 getDetCounts(monitormdws, runnumber, detid, vec2thetaMon, vecMonitorCounts, false);
157 // check
158 if (vecDetCounts.size() != vecMonitorCounts.size())
159 throw std::runtime_error("Number of detectors' counts' is different from that of "
160 "monitor counts.");
161
162 for (size_t i = 0; i < vecDetCounts.size(); ++i)
163 if (vecMonitorCounts[i] > 1.0E-9)
164 vecDetCounts[i] = vecDetCounts[i] / vecMonitorCounts[i];
165 else
166 vecDetCounts[i] = 0.;
167 }
168
169 // Sort and output
170 const size_t numdets = vecDetCounts.size();
171 std::vector<std::pair<double, double>> vecPair(numdets);
172 for (size_t i = 0; i < numdets; ++i) {
173 vecPair[i] = std::make_pair(vec2theta[i], vecDetCounts[i]);
174 }
175 std::sort(vecPair.begin(), vecPair.end());
176
177 // Apply to output
178 vecX.resize(numdets);
179 vecY.resize(numdets);
180 for (size_t i = 0; i < numdets; ++i) {
181 vecX[i] = vecPair[i].first;
182 vecY[i] = vecPair[i].second;
183 }
184
185 // Set up label
186 xlabel = "Degrees";
187 if (donormalize)
188 ylabel = "Noramlized Intensity";
189 else
190 ylabel = "Counts";
191}
192
193//----------------------------------------------------------------------------------------------
206 const API::IMDEventWorkspace_const_sptr &monitormdws,
207 const int detid, std::vector<double> &vecX,
208 std::vector<double> &vecY, std::string &xlabel,
209 std::string &ylabel, const bool &donormalize) {
210 // Get detector counts
211 std::vector<double> vecSampleLog;
212 std::vector<double> vecDetCounts;
213 int runnumber = -1;
214 bool get2theta = false;
215 if (xlabel.empty()) {
216 // xlabel is in default and thus use 2-theta for X
217 get2theta = true;
218 }
219 getDetCounts(datamdws, runnumber, detid, vecSampleLog, vecDetCounts, get2theta);
220 if (!get2theta) {
221 getSampleLogValues(datamdws, xlabel, runnumber, vecSampleLog);
222 }
223 if (vecSampleLog.size() != vecDetCounts.size())
224 throw std::runtime_error("Logic error! number of sample logs should be "
225 "same as number of detector's counts.");
226
227 // Get monitor counts
228 std::vector<double> vec2thetaMon;
229 std::vector<double> vecMonitorCounts;
230 // FIXME - Consider refactoring in future
231 // Normalize if required
232 if (donormalize) {
233 getDetCounts(monitormdws, runnumber, detid, vec2thetaMon, vecMonitorCounts, false);
234 if (vecDetCounts.size() != vecMonitorCounts.size())
235 throw std::runtime_error("Number of detectors' counts' is different from that of "
236 "monitor counts.");
237
238 for (size_t i = 0; i < vecDetCounts.size(); ++i)
239 if (vecMonitorCounts[i] > 1.0E-9)
240 vecDetCounts[i] = vecDetCounts[i] / vecMonitorCounts[i];
241 else
242 vecDetCounts[i] = 0.;
243 }
244
245 // Sort and output
246 const size_t numpts = vecDetCounts.size();
247 std::vector<std::pair<double, double>> vecPair(numpts);
248 for (size_t i = 0; i < numpts; ++i) {
249 vecPair[i] = std::make_pair(vecSampleLog[i], vecDetCounts[i]);
250 }
251 std::sort(vecPair.begin(), vecPair.end());
252
253 // Apply to output
254 vecX.resize(numpts);
255 vecY.resize(numpts);
256 for (size_t i = 0; i < numpts; ++i) {
257 vecX[i] = vecPair[i].first;
258 vecY[i] = vecPair[i].second;
259 }
260
261 // Set up label
262 if (get2theta)
263 xlabel = "Degrees";
264 if (donormalize)
265 ylabel = "Noramlized Intensity";
266 else
267 ylabel = "Counts";
268}
269
270//----------------------------------------------------------------------------------------------
281 const std::string &samplelogname, std::vector<double> &vecX,
282 std::vector<double> &vecY, std::string &xlabel,
283 std::string &ylabel) {
284 // prepare
285 vecX.clear();
286 vecY.clear();
287
288 // Y values
289 int runnumber = -1;
290 getSampleLogValues(datamdws, samplelogname, runnumber, vecY);
291 ylabel = samplelogname;
292
293 // X values
294 if (xlabel.empty()) {
295 // default
296 xlabel = "Pt.";
297 }
298 getSampleLogValues(datamdws, xlabel, runnumber, vecX);
299
300 if (vecX.size() != vecY.size())
301 throw std::runtime_error("It is not correct to have two different sizes vectors.");
302
303 // Sort
304 const size_t numpts = vecX.size();
305 std::vector<std::pair<double, double>> vecPair(numpts);
306 for (size_t i = 0; i < numpts; ++i) {
307 vecPair[i] = std::make_pair(vecX[i], vecY[i]);
308 }
309 std::sort(vecPair.begin(), vecPair.end());
310
311 // Apply to output
312 vecX.resize(numpts);
313 vecY.resize(numpts);
314 for (size_t i = 0; i < numpts; ++i) {
315 vecX[i] = vecPair[i].first;
316 vecY[i] = vecPair[i].second;
317 }
318}
319
320//----------------------------------------------------------------------------------------------
333 const int &detid, std::vector<double> &vecX, std::vector<double> &vecY,
334 bool formX) {
335 // Get sample and source position
336 if (mdws->getNumExperimentInfo() == 0)
337 throw std::runtime_error("There is no ExperimentInfo object that has been set to "
338 "input MDEventWorkspace!");
339
340 V3D samplepos;
341 V3D sourcepos;
342
343 if (formX) {
344 ExperimentInfo_const_sptr expinfo = mdws->getExperimentInfo(0);
345 Geometry::IComponent_const_sptr sample = expinfo->getInstrument()->getSample();
346 samplepos = sample->getPos();
347 g_log.debug() << "Sample position is " << samplepos.X() << ", " << samplepos.Y() << ", " << samplepos.Z() << "\n";
348
349 Geometry::IComponent_const_sptr source = expinfo->getInstrument()->getSource();
350 sourcepos = source->getPos();
351 g_log.debug() << "Source position is " << sourcepos.X() << "," << sourcepos.Y() << ", " << sourcepos.Z() << "\n";
352 vecX.clear();
353 }
354 vecY.clear();
355
356 // Go through all events to find out their positions
357 auto mditer = mdws->createIterator();
358
359 bool scancell = true;
360 size_t nextindex = 1;
361 while (scancell) {
362 // get the number of events of this cell
363 size_t numev2 = mditer->getNumEvents();
364 g_log.debug() << "MDWorkspace " << mdws->getName() << " Cell " << nextindex - 1 << ": Number of events = " << numev2
365 << " Does NEXT cell exist = " << mditer->next() << "\n";
366
367 // loop over all the events in current cell
368 for (size_t iev = 0; iev < numev2; ++iev) {
369 // filter out the events with uninterrested run numbers and detid
370 // runnumber/detid < 0 indicates that all run number or all detectors will
371 // be taken
372 int thisrunnumber = mditer->getInnerExpInfoIndex(iev);
373 if (runnumber >= 0 && thisrunnumber != runnumber)
374 continue;
375
376 int thisdetid = mditer->getInnerDetectorID(iev);
377 if (detid >= 0 && thisdetid != detid)
378 continue;
379
380 // get detector position for 2theta
381 if (formX) {
382 double tempx = mditer->getInnerPosition(iev, 0);
383 double tempy = mditer->getInnerPosition(iev, 1);
384 double tempz = mditer->getInnerPosition(iev, 2);
385 Kernel::V3D detpos(tempx, tempy, tempz);
386 Kernel::V3D v_det_sample = detpos - samplepos;
387 Kernel::V3D v_sample_src = samplepos - sourcepos;
388 double twotheta = v_det_sample.angle(v_sample_src) / M_PI * 180.;
389 vecX.emplace_back(twotheta);
390 }
391
392 // add new value to vecPair
393 double signal = mditer->getInnerSignal(iev);
394 vecY.emplace_back(signal);
395 } // ENDFOR (iev)
396
397 // Advance to next cell
398 if (mditer->next()) {
399 // advance to next cell
400 mditer->jumpTo(nextindex);
401 ++nextindex;
402 } else {
403 // break the loop
404 scancell = false;
405 }
406 } // ENDOF(while)
407}
408
409//----------------------------------------------------------------------------------------------
418 const std::string &samplelogname, const int runnumber,
419 std::vector<double> &vecSampleLog) {
420 // Clear input
421 vecSampleLog.clear();
422
423 // Loop over all experiment info objects
424 uint16_t numexpinfo = mdws->getNumExperimentInfo();
425 for (uint16_t iexp = 0; iexp < numexpinfo; ++iexp) {
426 // Get handler to experiment info object
427 ExperimentInfo_const_sptr expinfo = mdws->getExperimentInfo(iexp);
428 // Skip invalid run
429 int thisrunnumber = expinfo->getRunNumber();
430 if (thisrunnumber < 0)
431 continue;
432 // Skip unmatched run: input runnumber < 0 means that all run will be
433 // accepted
434 if (runnumber >= 0 && thisrunnumber != runnumber)
435 continue;
436 // Check property exists
437 if (!expinfo->run().hasProperty(samplelogname)) {
438 std::stringstream ess;
439 ess << "Workspace " << mdws->getName() << "'s " << iexp
440 << "-th ExperimentInfo with "
441 "run number "
442 << thisrunnumber << " does not have specified property " << samplelogname;
443 throw std::runtime_error(ess.str());
444 }
445 // Get experiment value
446 double logvalue = expinfo->run().getPropertyAsSingleValue(samplelogname);
447 vecSampleLog.emplace_back(logvalue);
448 g_log.debug() << "Add sample log (" << samplelogname << ") " << logvalue << " of " << iexp << "-th ExperimentInfo "
449 << "\n";
450 }
451}
452
453//----------------------------------------------------------------------------------------------
463 const std::vector<double> &vecY,
464 const std::string &xlabel,
465 const std::string &ylabel) {
466 // Create MatrixWorkspace
467 size_t sizex = vecX.size();
468 size_t sizey = vecY.size();
469 if (sizex != sizey || sizex == 0)
470 throw std::runtime_error("Unable to create output matrix workspace.");
471
473 std::dynamic_pointer_cast<MatrixWorkspace>(WorkspaceFactory::Instance().create("Workspace2D", 1, sizex, sizey));
474 if (!outws)
475 throw std::runtime_error("Failed to create output matrix workspace.");
476
477 // Set data
478 outws->setHistogram(0, Points(vecX), Counts(vecY));
479 auto &dataE = outws->mutableE(0);
480 std::replace_if(
481 dataE.begin(), dataE.end(), [](double val) { return val < 1.0; }, 1.0);
482
483 // Set label
484 outws->setYUnitLabel(ylabel);
485 if (!xlabel.empty()) {
486 try {
487 outws->getAxis(0)->setUnit(xlabel);
488 } catch (...) {
489 g_log.information() << "Label " << xlabel
490 << " for X-axis is not a unit "
491 "registered."
492 << "\n";
493 }
494 }
495
496 return outws;
497}
498
499} // namespace Mantid::MDAlgorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
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
static bool isEmpty(const NumT toCheck)
checks that the value was not set by users, uses the value in empty double/int.
A property class for workspaces.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void debug(const std::string &msg)
Logs at debug level.
Definition: Logger.cpp:114
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
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
double angle(const V3D &) const
Angle between this and another vector.
Definition: V3D.cpp:165
constexpr double Z() const noexcept
Get z.
Definition: V3D.h:234
GetSpiceDataRawCountsFromMD : Export raw detectors' counts or sample log values from IMDEventWorkspac...
void getSampleLogValues(const API::IMDEventWorkspace_const_sptr &mdws, const std::string &samplelogname, const int runnumber, std::vector< double > &vecSampleLog)
Get sample log values.
void getDetCounts(const API::IMDEventWorkspace_const_sptr &mdws, const int &runnumber, const int &detid, std::vector< double > &vecX, std::vector< double > &vecY, bool formX)
Get detectors' counts.
API::MatrixWorkspace_sptr createOutputWorkspace(const std::vector< double > &vecX, const std::vector< double > &vecY, const std::string &xlabel, const std::string &ylabel)
Create output workspace.
void exportDetCountsOfRun(const API::IMDEventWorkspace_const_sptr &datamdws, const API::IMDEventWorkspace_const_sptr &monitormdws, const int runnumber, std::vector< double > &vecX, std::vector< double > &vecY, std::string &xlabel, std::string &ylabel, bool donormalize)
Export all detectors' counts for a run.
void exportIndividualDetCounts(const API::IMDEventWorkspace_const_sptr &datamdws, const API::IMDEventWorkspace_const_sptr &monitormdws, const int detid, std::vector< double > &vecX, std::vector< double > &vecY, std::string &xlabel, std::string &ylabel, const bool &donormalize)
Export a detector's counts accross all runs.
void exportSampleLogValue(const API::IMDEventWorkspace_const_sptr &datamdws, const std::string &samplelogname, std::vector< double > &vecX, std::vector< double > &vecY, std::string &xlabel, std::string &ylabel)
Export sample log values accross all runs.
std::shared_ptr< IMDEventWorkspace > IMDEventWorkspace_sptr
Shared pointer to Mantid::API::IMDEventWorkspace.
std::shared_ptr< const IMDEventWorkspace > IMDEventWorkspace_const_sptr
Shared pointer to Mantid::API::IMDEventWorkspace (const version)
std::shared_ptr< const ExperimentInfo > ExperimentInfo_const_sptr
Shared pointer to const ExperimentInfo.
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
std::unique_ptr< T > create(const P &parent, const IndexArg &indexArg, const HistArg &histArg)
This is the create() method that all the other create() methods call.
std::shared_ptr< const IComponent > IComponent_const_sptr
Typdef of a shared pointer to a const IComponent.
Definition: IComponent.h:161
constexpr int EMPTY_INT() noexcept
Returns what we consider an "empty" integer within a property.
Definition: EmptyValues.h:25
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54