Mantid
Loading...
Searching...
No Matches
ConvertCWPDMDToSpectra.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#include <utility>
8
10
11#include "MantidAPI/Axis.h"
16#include "MantidAPI/Run.h"
22
23namespace Mantid::MDAlgorithms {
24
25using namespace Mantid::API;
26using namespace Mantid::Kernel;
27using namespace Mantid::MDAlgorithms;
28using HistogramData::BinEdges;
29using HistogramData::Counts;
30using HistogramData::CountStandardDeviations;
31
33
34const double BIGNUMBER = 1.0E100;
35
37
38 declareProperty(std::make_unique<WorkspaceProperty<IMDEventWorkspace>>("InputWorkspace", "", Direction::Input),
39 "Name of the input MDEventWorkspace that stores detectors "
40 "counts from a constant-wave powder diffraction experiment.");
41
42 declareProperty(std::make_unique<WorkspaceProperty<IMDEventWorkspace>>("InputMonitorWorkspace", "", Direction::Input),
43 "Name of the input MDEventWorkspace that stores monitor "
44 "counts from a constant-wave powder diffraciton experiment.");
45
46 declareProperty(std::make_unique<ArrayProperty<double>>("BinningParams"),
47 "A comma separated list of first bin boundary, width, last bin boundary. "
48 "Optionally\n"
49 "this can be followed by a comma and more widths and last boundary "
50 "pairs.\n"
51 "Negative width values indicate logarithmic binning.");
52
53 declareProperty(std::make_unique<WorkspaceProperty<MatrixWorkspace>>("OutputWorkspace", "", Direction::Output),
54 "Name of the output workspace for reduced data.");
55
56 std::array<std::string, 3> vecunits = {{"2theta", "dSpacing", "Momentum Transfer (Q)"}};
57 auto unitval = std::make_shared<ListValidator<std::string>>(vecunits);
58 declareProperty("UnitOutput", "2theta", unitval, "Unit of the output workspace.");
59
60 declareProperty("NeutronWaveLength", EMPTY_DBL(), "Constant wavelength of the neutrons from reactor source.");
61
62 declareProperty("NeutornWaveLengthPropertyName", "wavelength",
63 "Property name of the neutron wavelength in the sample log."
64 "If output unit is other than 2theta and NeutronWaveLength is not given,"
65 "then the neutron wavelength will be searched in sample logs by "
66 "name specified by this property.");
67
68 declareProperty("ScaleFactor", 1.0, "Scaling factor on the normalized counts.");
69
70 declareProperty(std::make_unique<ArrayProperty<int>>("ExcludedDetectorIDs"),
71 "A comma separated list of integers to indicate the IDs of "
72 "the detectors that will be excluded from binning.");
73
74 declareProperty("LinearInterpolateZeroCounts", true,
75 "If set to true and if a bin has zero count, a linear "
76 "interpolation will be made to set the value of this bin. It "
77 "is applied to the case that the bin size is small. ");
78}
79
81 // Process input workspaces
82 // input data workspace
83 IMDEventWorkspace_sptr inputDataWS = getProperty("InputWorkspace");
84 // input monitor workspace
85 IMDEventWorkspace_sptr inputMonitorWS = getProperty("InputMonitorWorkspace");
86 // input binning parameters
87 const std::vector<double> binParams = getProperty("BinningParams");
88 // scale factor
89 double scaleFactor = getProperty("ScaleFactor");
90 // do linear interpolation
91 bool doLinearInterpolation = getProperty("LinearInterpolateZeroCounts");
92 // unit
93 std::string outputunit = getProperty("UnitOutput");
94 double wavelength = getProperty("NeutronWaveLength");
95
96 std::vector<detid_t> excluded_detids = getProperty("ExcludedDetectorIDs");
97
98 // Validate inputs
99 // input data workspace and monitor workspace should match
100 size_t numdataevents = inputDataWS->getNEvents();
101 size_t nummonitorevents = inputMonitorWS->getNEvents();
102 if (numdataevents != nummonitorevents)
103 throw std::runtime_error("Input data workspace and monitor workspace have "
104 "different number of MDEvents.");
105
106 // output unit: make a map for wavelength
107 std::map<int, double> map_runWavelength;
108 if (outputunit != "2theta") {
109 // set up runid and wavelength map
110 std::string wavelengthpropertyname = getProperty("NeutornWaveLengthPropertyName");
111
112 uint16_t numexpinfo = inputDataWS->getNumExperimentInfo();
113 for (uint16_t iexp = 0; iexp < numexpinfo; ++iexp) {
114 int runid = std::stoi(inputDataWS->getExperimentInfo(iexp)->run().getProperty("run_number")->value());
115 // skip if run id is not a valid one
116 if (runid < 0)
117 continue;
118 double thislambda = wavelength;
119 if (inputDataWS->getExperimentInfo(iexp)->run().hasProperty(wavelengthpropertyname))
120 thislambda =
121 std::stod(inputDataWS->getExperimentInfo(iexp)->run().getProperty(wavelengthpropertyname)->value());
122 else if (wavelength == EMPTY_DBL()) {
123 std::stringstream errss;
124 errss << "In order to convert unit to " << outputunit
125 << ", either NeutronWaveLength "
126 " is to be specified or property "
127 << wavelengthpropertyname << " must exist for run " << runid << ".";
128 throw std::runtime_error(errss.str());
129 }
130 map_runWavelength.emplace(runid, thislambda);
131 }
132 }
133
134 // bin parameters
135 double xmin, xmax, binsize;
136 xmin = xmax = binsize = -1;
137 if (binParams.size() == 1) {
138 binsize = binParams[0];
139 g_log.warning() << "Only bin size " << binParams[0]
140 << " is specified. Xmin and Xmax "
141 " will be calcualted from motor positions and wavelength. "
142 "More CPU time will be used."
143 << "\n";
144 } else if (binParams.size() == 3) {
145 xmin = binParams[0];
146 binsize = binParams[1];
147 xmax = binParams[2];
148 if (xmin >= xmax)
149 throw std::runtime_error("Min value of the bin must be smaller than maximum value.");
150 } else {
151 // Either 1 or 3 parameters. Throw exception
152 throw std::runtime_error("Binning parameters must have either 1 or 3 items.");
153 }
154
155 // Rebin
156 std::sort(excluded_detids.begin(), excluded_detids.end());
157 API::MatrixWorkspace_sptr outws = reducePowderData(inputDataWS, inputMonitorWS, outputunit, map_runWavelength, xmin,
158 xmax, binsize, doLinearInterpolation, excluded_detids);
159
160 // Scale
161 scaleMatrixWorkspace(outws, scaleFactor, m_infitesimal);
162
163 // Set up the sample logs
164 setupSampleLogs(outws, inputDataWS);
165
166 // Return
167 setProperty("OutputWorkspace", outws);
168}
169
170//----------------------------------------------------------------------------------------------
192 const std::string &targetunit, const std::map<int, double> &map_runwavelength, const double xmin, const double xmax,
193 const double binsize, bool dolinearinterpolation, const std::vector<detid_t> &vec_excludeddets) {
194 // Get some information
195 int64_t numevents = dataws->getNEvents();
196
197 // check xmin and xmax
198 double lowerboundary, upperboundary;
199 if (xmin < 0 || xmax < 0) {
200 // xmin or xmax cannot be negative (2theta, dspace and q are always
201 // positive)
202 findXBoundary(dataws, targetunit, map_runwavelength, lowerboundary, upperboundary);
203 } else {
204 lowerboundary = xmin;
205 upperboundary = xmax;
206 }
207
208 g_log.debug() << "Binning from " << lowerboundary << " to " << upperboundary << "\n";
209
210 // Create bins in 2theta (degree)
211 size_t sizex, sizey;
212 sizex = static_cast<size_t>((upperboundary - lowerboundary) / binsize + 1);
213 if (lowerboundary + static_cast<double>(sizex) * binsize < upperboundary)
214 ++lowerboundary;
215
216 sizey = sizex - 1;
217 g_log.debug() << "Number of events = " << numevents << ", bin size = " << binsize << ", SizeX = " << sizex << ", "
218 << ", SizeY = " << sizey << ", Delta = " << upperboundary - lowerboundary << ", Bin size = " << binsize
219 << ", sizex_d = " << (upperboundary - lowerboundary) / binsize + 2 << "\n";
220 std::vector<double> vecx(sizex), vecy(sizex - 1, 0), vecm(sizex - 1, 0), vece(sizex - 1, 0);
221
222 for (size_t i = 0; i < sizex; ++i) {
223 vecx[i] = lowerboundary + static_cast<double>(i) * binsize;
224 }
225
226 // Convert unit to unit char bit
227 char unitchar = 't'; // default 2theta
228 if (targetunit == "dSpacing")
229 unitchar = 'd';
230 else if (targetunit == "Momentum Transfer (Q)")
231 unitchar = 'q';
232
233 binMD(dataws, unitchar, map_runwavelength, vecx, vecy, vec_excludeddets);
234 binMD(monitorws, unitchar, map_runwavelength, vecx, vecm, vec_excludeddets);
235
236 // Normalize by division
237 double maxmonitorcounts = 0;
238 for (size_t i = 0; i < vecm.size(); ++i) {
239 if (vecm[i] >= 1.) {
240 double y = vecy[i];
241 double ey = sqrt(y);
242 double m = vecm[i];
243 double em = sqrt(m);
244 vecy[i] = y / m;
245 // using standard deviation's error propagation
246 vece[i] = vecy[i] * sqrt((ey / y) * (ey / y) + (em / m) * (em / m));
247 // maximum monitor counts
248 if (m > maxmonitorcounts)
249 maxmonitorcounts = m;
250 } else {
251 vecy[i] = 0.0;
252 vece[i] = 1.0;
253 }
254 }
255
256 // Create workspace and set values
257 API::MatrixWorkspace_sptr pdws = WorkspaceFactory::Instance().create("Workspace2D", 1, sizex, sizey);
258 // Set unit
259 pdws->setYUnitLabel("Intensity");
260 if (unitchar == 'd')
261 pdws->getAxis(0)->setUnit("dSpacing");
262 else if (unitchar == 'q')
263 pdws->getAxis(0)->setUnit("MomentumTransfer");
264 else {
265 // Twotheta
266 pdws->getAxis(0)->setUnit("Degrees");
267 }
268
269 pdws->setHistogram(0, BinEdges(vecx), Counts(vecy), CountStandardDeviations(vece));
270
271 // Interpolation
272 m_infitesimal = 0.1 / (maxmonitorcounts);
273
274 if (dolinearinterpolation)
276
277 return pdws;
278}
279
280//----------------------------------------------------------------------------------------------
290 const std::string &targetunit,
291 const std::map<int, double> &map_runwavelength, double &xmin, double &xmax) {
292 // Go through all instruments
293 uint16_t numruns = dataws->getNumExperimentInfo();
294
295 xmin = BIGNUMBER;
296 xmax = -1;
297
298 for (uint16_t irun = 0; irun < numruns; ++irun) {
299 // Skip the Experiment Information does not have run
300 if (!dataws->getExperimentInfo(irun)->getInstrument()) {
301 g_log.warning() << "iRun = " << irun << " of total " << numruns << " does not have instrument associated"
302 << "\n";
303 continue;
304 }
305
306 // Get run number
307 int runnumber = dataws->getExperimentInfo(irun)->getRunNumber();
308 g_log.debug() << "Run " << runnumber << ": ";
309 auto miter = map_runwavelength.find(runnumber);
310 double wavelength = -1;
311 if (miter != map_runwavelength.end()) {
312 wavelength = miter->second;
313 g_log.debug() << " wavelength = " << wavelength << "\n";
314 } else {
315 g_log.debug() << " no matched wavelength."
316 << "\n";
317 }
318
319 // Get source and sample position
320 std::vector<detid_t> vec_detid = dataws->getExperimentInfo(irun)->getInstrument()->getDetectorIDs(true);
321 if (vec_detid.empty()) {
322 g_log.information() << "Run " << runnumber << " has no detectors."
323 << "\n";
324 continue;
325 }
326 const V3D samplepos = dataws->getExperimentInfo(irun)->getInstrument()->getSample()->getPos();
327 const V3D sourcepos = dataws->getExperimentInfo(irun)->getInstrument()->getSource()->getPos();
328
329 // Get all detectors
330 // std::vector<detid_t> vec_detid =
331 // dataws->getExperimentInfo(irun)->getInstrument()->getDetectorIDs(true);
332 std::vector<Geometry::IDetector_const_sptr> vec_det =
333 dataws->getExperimentInfo(irun)->getInstrument()->getDetectors(vec_detid);
334 size_t numdets = vec_det.size();
335 g_log.debug() << "Run = " << runnumber << ": Number of detectors = " << numdets << "\n";
336
337 // Scan all the detectors to get Xmin and Xmax
338 for (size_t idet = 0; idet < numdets; ++idet) {
339 Geometry::IDetector_const_sptr tmpdet = vec_det[idet];
340 const V3D detpos = tmpdet->getPos();
341
342 double R, theta, phi;
343 detpos.getSpherical(R, theta, phi);
344 if (R < 0.0001)
345 g_log.error("Invalid detector position");
346
347 Kernel::V3D v_det_sample = detpos - samplepos;
348 Kernel::V3D v_sample_src = samplepos - sourcepos;
349 double twotheta = calculate2Theta(v_det_sample, v_sample_src) / M_PI * 180.;
350
351 // convert unit optionally
352 double outx = -1;
353 if (targetunit == "2theta")
354 outx = twotheta;
355 else {
356 if (wavelength <= 0)
357 throw std::runtime_error("Wavelength is not defined!");
358
359 if (targetunit == "dSpacing")
360 outx = calculateDspaceFrom2Theta(twotheta, wavelength);
361 else if (targetunit == "Momentum Transfer (Q)")
362 outx = calculateQFrom2Theta(twotheta, wavelength);
363 else
364 throw std::runtime_error("Unrecognized unit.");
365 }
366
367 // Compare with xmin and xmax
368 if (outx < xmin)
369 xmin = outx;
370 if (outx > xmax)
371 xmax = outx;
372 }
373 }
374
375 g_log.debug() << "Find boundary for unit " << targetunit << ": [" << xmin << ", " << xmax << "]"
376 << "\n";
377}
378
379//----------------------------------------------------------------------------------------------
390 const std::map<int, double> &map_runlambda, const std::vector<double> &vecx,
391 std::vector<double> &vecy, const std::vector<detid_t> &vec_excludedet) {
392 // Check whether MD workspace has proper instrument and experiment Info
393 if (mdws->getNumExperimentInfo() == 0)
394 throw std::runtime_error("There is no ExperimentInfo object that has been set to "
395 "input MDEventWorkspace!");
396 else
397 g_log.information() << "Number of ExperimentInfo objects of MDEventWrokspace is " << mdws->getNumExperimentInfo()
398 << "\n";
399
400 // Get sample position
401 ExperimentInfo_const_sptr expinfo = mdws->getExperimentInfo(0);
402 Geometry::IComponent_const_sptr sample = expinfo->getInstrument()->getSample();
403 const V3D samplepos = sample->getPos();
404 g_log.debug() << "Sample position is " << samplepos.X() << ", " << samplepos.Y() << ", " << samplepos.Z() << "\n";
405
406 Geometry::IComponent_const_sptr source = expinfo->getInstrument()->getSource();
407 const V3D sourcepos = source->getPos();
408 g_log.debug() << "Source position is " << sourcepos.X() << "," << sourcepos.Y() << ", " << sourcepos.Z() << "\n";
409
410 // Go through all events to find out their positions
411 auto mditer = mdws->createIterator();
412 bool scancell = true;
413 size_t nextindex = 1;
414 int currExpInfoIndex = -1;
415 double currWavelength = -1;
416 while (scancell) {
417 // get the number of events of this cell
418 size_t numev2 = mditer->getNumEvents();
419 g_log.debug() << "MDWorkspace " << mdws->getName() << " Cell " << nextindex - 1 << ": Number of events = " << numev2
420 << " Does NEXT cell exist = " << mditer->next() << "\n";
421
422 // loop over all the events in current cell
423 for (size_t iev = 0; iev < numev2; ++iev) {
424 // get detector position for 2theta
425 detid_t detid = mditer->getInnerDetectorID(iev);
426 if (isExcluded(vec_excludedet, detid)) {
427 g_log.debug() << "Detector " << detid << " is excluded. Signal = " << mditer->getInnerSignal(iev) << "\n";
428 continue;
429 }
430
431 double tempx = mditer->getInnerPosition(iev, 0);
432 double tempy = mditer->getInnerPosition(iev, 1);
433 double tempz = mditer->getInnerPosition(iev, 2);
434 Kernel::V3D detpos(tempx, tempy, tempz);
435 Kernel::V3D v_det_sample = detpos - samplepos;
436 Kernel::V3D v_sample_src = samplepos - sourcepos;
437 double twotheta = calculate2Theta(v_det_sample, v_sample_src) / M_PI * 180.;
438
439 // convert unit optionally
440 auto temprun = static_cast<int>(mditer->getInnerExpInfoIndex(iev));
441 double outx;
442 if (unitbit == 't')
443 outx = twotheta;
444 else {
445 if (temprun != currExpInfoIndex) {
446 // use map to find a new wavelength
447 auto miter = map_runlambda.find(temprun);
448 if (miter == map_runlambda.end()) {
449 std::stringstream errss;
450 errss << "Event " << iev << " has run ID as " << temprun << ". "
451 << "It has no corresponding ExperimentInfo in MDWorkspace " << mdws->getName() << ".";
452 throw std::runtime_error(errss.str());
453 }
454 currWavelength = miter->second;
455 }
456 if (unitbit == 'd')
457 outx = calculateDspaceFrom2Theta(twotheta, currWavelength);
458 else
459 outx = calculateQFrom2Theta(twotheta, currWavelength);
460 }
461
462 // get signal and assign signal to bin
463 int xindex;
464 const double SMALL = 1.0E-5;
465 if (outx + SMALL < vecx.front()) {
466 // Significantly out of left boundary
467 xindex = -1;
468 } else if (fabs(outx - vecx.front()) < SMALL) {
469 // Almost on the left boundary
470 xindex = 0;
471 } else if (outx - SMALL > vecx.back()) {
472 // Significantly out of right boundary
473 xindex = static_cast<int>(vecx.size());
474 } else if (fabs(outx - vecx.back()) < SMALL) {
475 // Right on the right boundary
476 xindex = static_cast<int>(vecy.size()) - 1;
477 } else {
478 // Other situation
479 auto vfiter = std::lower_bound(vecx.cbegin(), vecx.cend(), outx);
480 xindex = static_cast<int>(vfiter - vecx.cbegin());
481 if ((xindex < static_cast<int>(vecx.size())) && (outx + 1.0E-5 < vecx[xindex])) {
482 // assume the bin's boundaries are of [...) and consider numerical
483 // error
484 xindex -= 1;
485 } else {
486 g_log.debug() << "Case for almost same. Event X = " << outx << ", Boundary = " << vecx[xindex] << "\n";
487 }
488 if (xindex >= static_cast<int>(vecy.size())) {
489 g_log.warning() << "Case unexpected: Event X = " << outx << ", Boundary = " << vecx[xindex] << "\n";
490 } else if (xindex < 0) {
491 g_log.warning() << "Case unexpected: Event X = " << outx << ", Boundary index is out of vector range.\n";
492 }
493 }
494
495 // add signal
496 if (xindex < 0) {
497 // Out of left boundary
498 int32_t innerDetid = mditer->getInnerDetectorID(iev);
499 uint16_t runid = mditer->getInnerExpInfoIndex(iev);
500 g_log.debug() << "Event is out of user-specified range by " << (outx - vecx.front()) << ", xindex = " << xindex
501 << ", " << unitbit << " = " << outx << " out of left boundeary [" << vecx.front() << ", "
502 << vecx.back() << "]. dep pos = " << detpos.X() << ", " << detpos.Y() << ", " << detpos.Z()
503 << ", Run = " << runid << ", DetectorID = " << innerDetid << "\n";
504 continue;
505 } else if (xindex >= static_cast<int>(vecy.size())) {
506 // Out of right boundary
507 int32_t innerDetid = mditer->getInnerDetectorID(iev);
508 uint16_t runid = mditer->getInnerExpInfoIndex(iev);
509 g_log.debug() << "Event is out of user-specified range "
510 << "xindex = " << xindex << ", " << unitbit << " = " << outx << " out of [" << vecx.front()
511 << ", " << vecx.back() << "]. dep pos = " << detpos.X() << ", " << detpos.Y() << ", "
512 << detpos.Z() << "; sample pos = " << samplepos.X() << ", " << samplepos.Y() << ", "
513 << samplepos.Z() << ", Run = " << runid << ", DetectorID = " << innerDetid << "\n";
514 continue;
515 } else {
516 double signal = mditer->getInnerSignal(iev);
517 vecy[xindex] += signal;
518 }
519 }
520
521 // Advance to next cell
522 if (mditer->next()) {
523 // advance to next cell
524 mditer->jumpTo(nextindex);
525 ++nextindex;
526 } else {
527 // break the loop
528 scancell = false;
529 }
530 } // ENDOF(while)
531}
532
533//----------------------------------------------------------------------------------------------
542 const double &infinitesimal) {
543 g_log.debug() << "Number of spectrum = " << matrixws->getNumberHistograms() << " Infinitesimal = " << infinitesimal
544 << "\n";
545 size_t numspec = matrixws->getNumberHistograms();
546 for (size_t i = 0; i < numspec; ++i) {
547 // search for the first nonzero value and last nonzero value
548 bool onsearch = true;
549 size_t minNonZeroIndex = 0;
550 while (onsearch) {
551 if (matrixws->readY(i)[minNonZeroIndex] > infinitesimal)
552 onsearch = false;
553 else
554 ++minNonZeroIndex;
555
556 if (minNonZeroIndex == matrixws->readY(i).size())
557 onsearch = false;
558 }
559 size_t maxNonZeroIndex = matrixws->readY(i).size() - 1;
560 onsearch = true;
561 while (onsearch) {
562 if (matrixws->readY(i)[maxNonZeroIndex] > infinitesimal)
563 onsearch = false;
564 else if (maxNonZeroIndex == 0)
565 onsearch = false;
566 else
567 --maxNonZeroIndex;
568 }
569 g_log.debug() << "iMinNonZero = " << minNonZeroIndex << ", iMaxNonZero = " << maxNonZeroIndex
570 << " Workspace index = " << i << ", Y size = " << matrixws->readY(i).size() << "\n";
571 if (minNonZeroIndex >= maxNonZeroIndex)
572 throw std::runtime_error("It is not right!");
573
574 // Do linear interpolation for zero count values
575 for (size_t j = minNonZeroIndex + 1; j < maxNonZeroIndex; ++j) {
576 if (matrixws->readY(i)[j] < infinitesimal) {
577 // Do interpolation
578 // gives y = y_0 + (y_1-y_0)\frac{x - x_0}{x_1-x_0}
579
580 double leftx = matrixws->readX(i)[j - 1];
581 double lefty = matrixws->readY(i)[j - 1];
582 bool findnonzeroy = true;
583 size_t iright = j + 1;
584 while (findnonzeroy) {
585 if (matrixws->readY(i)[iright] > infinitesimal)
586 findnonzeroy = false;
587 else
588 ++iright;
589 }
590 double rightx = matrixws->readX(i)[iright];
591 double righty = matrixws->readY(i)[iright];
592 double curx = matrixws->readX(i)[j];
593 double curinterpoy = lefty + (righty - lefty) * (curx - leftx) / (rightx - leftx);
594 matrixws->dataY(i)[j] = curinterpoy;
595 matrixws->dataE(i)[j] = sqrt(curinterpoy);
596 }
597 }
598
599 return;
600 }
601}
602
603//----------------------------------------------------------------------------------------------
610 const API::IMDEventWorkspace_const_sptr &inputmdws) {
611 // get hold of the last experiment info from md workspace to copy over
612 auto lastindex = static_cast<uint16_t>(inputmdws->getNumExperimentInfo() - 1);
613 ExperimentInfo_const_sptr lastexpinfo = inputmdws->getExperimentInfo(lastindex);
614
615 // get hold of experiment info from matrix ws
616 Run &targetrun = matrixws->mutableRun();
617 const Run &srcrun = lastexpinfo->run();
618
619 const std::vector<Kernel::Property *> &vec_srcprop = srcrun.getProperties();
620 for (auto p : vec_srcprop) {
621 targetrun.addProperty(p->clone());
622 g_log.debug() << "Cloned property " << p->name() << "\n";
623 }
624}
625
626//----------------------------------------------------------------------------------------------
633void ConvertCWPDMDToSpectra::scaleMatrixWorkspace(const API::MatrixWorkspace_sptr &matrixws, const double &scalefactor,
634 const double &infinitesimal) {
635 size_t numspec = matrixws->getNumberHistograms();
636 for (size_t iws = 0; iws < numspec; ++iws) {
637 MantidVec &datay = matrixws->dataY(iws);
638 MantidVec &datae = matrixws->dataE(iws);
639 size_t numelements = datay.size();
640 for (size_t i = 0; i < numelements; ++i) {
641 // bin with zero counts is not scaled up
642 if (datay[i] >= infinitesimal) {
643 datay[i] *= scalefactor;
644 datae[i] *= scalefactor;
645 }
646 }
647 } // FOR(iws)
648}
649
650//----------------------------------------------------------------------------------------------
657bool ConvertCWPDMDToSpectra::isExcluded(const std::vector<detid_t> &vec_excludedet, const detid_t detid) {
658
659 return std::find(vec_excludedet.begin(), vec_excludedet.end(), detid) != vec_excludedet.end();
660}
661
662} // namespace Mantid::MDAlgorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
#define fabs(x)
Definition: Matrix.cpp:22
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
const std::vector< Kernel::Property * > & getProperties() const
Return all of the current properties.
Definition: LogManager.cpp:287
void addProperty(Kernel::Property *prop, bool overwrite=false)
Add data to the object in the form of a property.
Definition: LogManager.h:79
This class stores information regarding an experimental run as a series of log entries.
Definition: Run.h:38
A property class for workspaces.
Support for a property that holds an array of values.
Definition: ArrayProperty.h:28
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 error(const std::string &msg)
Logs at error level.
Definition: Logger.cpp:77
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
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
void getSpherical(double &R, double &theta, double &phi) const noexcept
Return the vector's position in spherical coordinates.
Definition: V3D.cpp:117
constexpr double Z() const noexcept
Get z.
Definition: V3D.h:234
ConvertCWPDMDToSpectra : Convert one MDWorkspaces containing reactor-source powder diffractometer's d...
void linearInterpolation(const API::MatrixWorkspace_sptr &matrixws, const double &infinitesimal)
Do linear interpolation to zero counts if bin is too small.
void setupSampleLogs(const API::MatrixWorkspace_sptr &matrixws, const API::IMDEventWorkspace_const_sptr &inputmdws)
Set up sample logs.
void findXBoundary(const API::IMDEventWorkspace_const_sptr &dataws, const std::string &targetunit, const std::map< int, double > &map_runwavelength, double &xmin, double &xmax)
Find the binning boundary according to detectors' positions.
API::MatrixWorkspace_sptr reducePowderData(const API::IMDEventWorkspace_const_sptr &dataws, const API::IMDEventWorkspace_const_sptr &monitorws, const std::string &targetunit, const std::map< int, double > &map_runwavelength, const double xmin, const double xmax, const double binsize, bool dolinearinterpolation, const std::vector< detid_t > &vec_excludeddets)
Main algorithm to reduce powder diffraction data.
void binMD(const API::IMDEventWorkspace_const_sptr &mdws, const char &unitbit, const std::map< int, double > &map_runlambda, const std::vector< double > &vecx, std::vector< double > &vecy, const std::vector< detid_t > &vec_excludedet)
Bin signals to its 2theta position.
bool isExcluded(const std::vector< detid_t > &vec_excludedet, const detid_t detid)
Check whether a detector is excluded.
void scaleMatrixWorkspace(const API::MatrixWorkspace_sptr &matrixws, const double &scalefactor, const double &infinitesimal)
Scale reduced data.
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::shared_ptr< const IComponent > IComponent_const_sptr
Typdef of a shared pointer to a const IComponent.
Definition: IComponent.h:161
std::shared_ptr< const Mantid::Geometry::IDetector > IDetector_const_sptr
Shared pointer to IDetector (const version)
Definition: IDetector.h:102
double calculateDspaceFrom2Theta(const double &twotheta, const double &wavelength)
Calculate d-space value from detector's position (2theta/theta) and wavelength.
double calculateQFrom2Theta(const double &twotheta, const double &wavelength)
Calculate Q value from detector's positin (2theta/theta) and wavelength q = 2 k sin(theta) = (4 pi)/(...
double calculate2Theta(const Kernel::V3D &detpos, const Kernel::V3D &samplepos)
Calculate 2theta value of detector postion from sample position.
int32_t detid_t
Typedef for a detector ID.
Definition: SpectrumInfo.h:21
std::vector< double > MantidVec
typedef for the data storage used in Mantid matrix workspaces
Definition: cow_ptr.h:172
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