28using HistogramData::BinEdges;
29using HistogramData::Counts;
30using HistogramData::CountStandardDeviations;
39 "Name of the input MDEventWorkspace that stores detectors "
40 "counts from a constant-wave powder diffraction experiment.");
43 "Name of the input MDEventWorkspace that stores monitor "
44 "counts from a constant-wave powder diffraciton experiment.");
47 "A comma separated list of first bin boundary, width, last bin boundary. "
49 "this can be followed by a comma and more widths and last boundary "
51 "Negative width values indicate logarithmic binning.");
54 "Name of the output workspace for reduced data.");
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.");
60 declareProperty(
"NeutronWaveLength",
EMPTY_DBL(),
"Constant wavelength of the neutrons from reactor source.");
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.");
68 declareProperty(
"ScaleFactor", 1.0,
"Scaling factor on the normalized counts.");
71 "A comma separated list of integers to indicate the IDs of "
72 "the detectors that will be excluded from binning.");
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. ");
87 const std::vector<double> binParams =
getProperty(
"BinningParams");
91 bool doLinearInterpolation =
getProperty(
"LinearInterpolateZeroCounts");
94 double wavelength =
getProperty(
"NeutronWaveLength");
96 std::vector<detid_t> excluded_detids =
getProperty(
"ExcludedDetectorIDs");
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.");
107 std::map<int, double> map_runWavelength;
108 if (outputunit !=
"2theta") {
110 std::string wavelengthpropertyname =
getProperty(
"NeutornWaveLengthPropertyName");
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());
118 double thislambda = wavelength;
119 if (inputDataWS->getExperimentInfo(iexp)->run().hasProperty(wavelengthpropertyname))
121 std::stod(inputDataWS->getExperimentInfo(iexp)->run().getProperty(wavelengthpropertyname)->value());
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());
130 map_runWavelength.emplace(runid, thislambda);
135 double xmin, xmax, binsize;
136 xmin = xmax = binsize = -1;
137 if (binParams.size() == 1) {
138 binsize = binParams[0];
140 <<
" is specified. Xmin and Xmax "
141 " will be calcualted from motor positions and wavelength. "
142 "More CPU time will be used."
144 }
else if (binParams.size() == 3) {
146 binsize = binParams[1];
149 throw std::runtime_error(
"Min value of the bin must be smaller than maximum value.");
152 throw std::runtime_error(
"Binning parameters must have either 1 or 3 items.");
156 std::sort(excluded_detids.begin(), excluded_detids.end());
158 xmax, binsize, doLinearInterpolation, excluded_detids);
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) {
195 int64_t numevents = dataws->getNEvents();
198 double lowerboundary, upperboundary;
199 if (xmin < 0 || xmax < 0) {
202 findXBoundary(dataws, targetunit, map_runwavelength, lowerboundary, upperboundary);
204 lowerboundary = xmin;
205 upperboundary = xmax;
208 g_log.
debug() <<
"Binning from " << lowerboundary <<
" to " << upperboundary <<
"\n";
212 sizex =
static_cast<size_t>((upperboundary - lowerboundary) / binsize + 1);
213 if (lowerboundary +
static_cast<double>(sizex) * binsize < upperboundary)
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);
222 for (
size_t i = 0; i < sizex; ++i) {
223 vecx[i] = lowerboundary +
static_cast<double>(i) * binsize;
228 if (targetunit ==
"dSpacing")
230 else if (targetunit ==
"Momentum Transfer (Q)")
233 binMD(dataws, unitchar, map_runwavelength, vecx, vecy, vec_excludeddets);
234 binMD(monitorws, unitchar, map_runwavelength, vecx, vecm, vec_excludeddets);
237 double maxmonitorcounts = 0;
238 for (
size_t i = 0; i < vecm.size(); ++i) {
246 vece[i] = vecy[i] * sqrt((ey /
y) * (ey /
y) + (em /
m) * (em /
m));
248 if (
m > maxmonitorcounts)
249 maxmonitorcounts =
m;
259 pdws->setYUnitLabel(
"Intensity");
261 pdws->getAxis(0)->setUnit(
"dSpacing");
262 else if (unitchar ==
'q')
263 pdws->getAxis(0)->setUnit(
"MomentumTransfer");
266 pdws->getAxis(0)->setUnit(
"Degrees");
269 pdws->setHistogram(0, BinEdges(vecx), Counts(vecy), CountStandardDeviations(vece));
274 if (dolinearinterpolation)
290 const std::string &targetunit,
291 const std::map<int, double> &map_runwavelength,
double &xmin,
double &xmax) {
293 uint16_t numruns = dataws->getNumExperimentInfo();
298 for (uint16_t irun = 0; irun < numruns; ++irun) {
300 if (!dataws->getExperimentInfo(irun)->getInstrument()) {
301 g_log.
warning() <<
"iRun = " << irun <<
" of total " << numruns <<
" does not have instrument associated"
307 int runnumber = dataws->getExperimentInfo(irun)->getRunNumber();
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";
320 std::vector<detid_t> vec_detid = dataws->getExperimentInfo(irun)->getInstrument()->getDetectorIDs(
true);
321 if (vec_detid.empty()) {
326 const V3D samplepos = dataws->getExperimentInfo(irun)->getInstrument()->getSample()->getPos();
327 const V3D sourcepos = dataws->getExperimentInfo(irun)->getInstrument()->getSource()->getPos();
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";
338 for (
size_t idet = 0; idet < numdets; ++idet) {
340 const V3D detpos = tmpdet->getPos();
342 double R, theta, phi;
349 double twotheta =
calculate2Theta(v_det_sample, v_sample_src) / M_PI * 180.;
353 if (targetunit ==
"2theta")
357 throw std::runtime_error(
"Wavelength is not defined!");
359 if (targetunit ==
"dSpacing")
361 else if (targetunit ==
"Momentum Transfer (Q)")
364 throw std::runtime_error(
"Unrecognized unit.");
375 g_log.
debug() <<
"Find boundary for unit " << targetunit <<
": [" << xmin <<
", " << xmax <<
"]"
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) {
393 if (mdws->getNumExperimentInfo() == 0)
394 throw std::runtime_error(
"There is no ExperimentInfo object that has been set to "
395 "input MDEventWorkspace!");
397 g_log.
information() <<
"Number of ExperimentInfo objects of MDEventWrokspace is " << mdws->getNumExperimentInfo()
403 const V3D samplepos = sample->getPos();
404 g_log.
debug() <<
"Sample position is " << samplepos.
X() <<
", " << samplepos.
Y() <<
", " << samplepos.
Z() <<
"\n";
407 const V3D sourcepos = source->getPos();
408 g_log.
debug() <<
"Source position is " << sourcepos.
X() <<
"," << sourcepos.
Y() <<
", " << sourcepos.
Z() <<
"\n";
411 auto mditer = mdws->createIterator();
412 bool scancell =
true;
413 size_t nextindex = 1;
414 int currExpInfoIndex = -1;
415 double currWavelength = -1;
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";
423 for (
size_t iev = 0; iev < numev2; ++iev) {
425 detid_t detid = mditer->getInnerDetectorID(iev);
427 g_log.
debug() <<
"Detector " << detid <<
" is excluded. Signal = " << mditer->getInnerSignal(iev) <<
"\n";
431 double tempx = mditer->getInnerPosition(iev, 0);
432 double tempy = mditer->getInnerPosition(iev, 1);
433 double tempz = mditer->getInnerPosition(iev, 2);
437 double twotheta =
calculate2Theta(v_det_sample, v_sample_src) / M_PI * 180.;
440 auto temprun =
static_cast<int>(mditer->getInnerExpInfoIndex(iev));
445 if (temprun != currExpInfoIndex) {
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());
454 currWavelength = miter->second;
464 const double SMALL = 1.0E-5;
465 if (outx + SMALL < vecx.front()) {
468 }
else if (
fabs(outx - vecx.front()) < SMALL) {
471 }
else if (outx - SMALL > vecx.back()) {
473 xindex =
static_cast<int>(vecx.size());
474 }
else if (
fabs(outx - vecx.back()) < SMALL) {
476 xindex =
static_cast<int>(vecy.size()) - 1;
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])) {
486 g_log.
debug() <<
"Case for almost same. Event X = " << outx <<
", Boundary = " << vecx[xindex] <<
"\n";
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";
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";
505 }
else if (xindex >=
static_cast<int>(vecy.size())) {
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";
516 double signal = mditer->getInnerSignal(iev);
517 vecy[xindex] += signal;
522 if (mditer->next()) {
524 mditer->jumpTo(nextindex);
542 const double &infinitesimal) {
543 g_log.
debug() <<
"Number of spectrum = " << matrixws->getNumberHistograms() <<
" Infinitesimal = " << infinitesimal
545 size_t numspec = matrixws->getNumberHistograms();
546 for (
size_t i = 0; i < numspec; ++i) {
548 bool onsearch =
true;
549 size_t minNonZeroIndex = 0;
551 if (matrixws->readY(i)[minNonZeroIndex] > infinitesimal)
556 if (minNonZeroIndex == matrixws->readY(i).size())
559 size_t maxNonZeroIndex = matrixws->readY(i).size() - 1;
562 if (matrixws->readY(i)[maxNonZeroIndex] > infinitesimal)
564 else if (maxNonZeroIndex == 0)
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!");
575 for (
size_t j = minNonZeroIndex + 1; j < maxNonZeroIndex; ++j) {
576 if (matrixws->readY(i)[j] < infinitesimal) {
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;
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);
612 auto lastindex =
static_cast<uint16_t
>(inputmdws->getNumExperimentInfo() - 1);
616 Run &targetrun = matrixws->mutableRun();
617 const Run &srcrun = lastexpinfo->run();
619 const std::vector<Kernel::Property *> &vec_srcprop = srcrun.
getProperties();
620 for (
auto p : vec_srcprop) {
622 g_log.
debug() <<
"Cloned property " << p->name() <<
"\n";
634 const double &infinitesimal) {
635 size_t numspec = matrixws->getNumberHistograms();
636 for (
size_t iws = 0; iws < numspec; ++iws) {
639 size_t numelements = datay.size();
640 for (
size_t i = 0; i < numelements; ++i) {
642 if (datay[i] >= infinitesimal) {
643 datay[i] *= scalefactor;
644 datae[i] *= scalefactor;
659 return std::find(vec_excludedet.begin(), vec_excludedet.end(), detid) != vec_excludedet.end();
#define DECLARE_ALGORITHM(classname)
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
const std::vector< Kernel::Property * > & getProperties() const
Return all of the current properties.
void addProperty(Kernel::Property *prop, bool overwrite=false)
Add data to the object in the form of a property.
This class stores information regarding an experimental run as a series of log entries.
A property class for workspaces.
Support for a property that holds an array of values.
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.
void error(const std::string &msg)
Logs at error level.
void warning(const std::string &msg)
Logs at warning level.
void information(const std::string &msg)
Logs at information level.
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
constexpr double X() const noexcept
Get x.
constexpr double Y() const noexcept
Get y.
void getSpherical(double &R, double &theta, double &phi) const noexcept
Return the vector's position in spherical coordinates.
constexpr double Z() const noexcept
Get z.
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.
double m_infitesimal
Infinitesimal value.
void exec() override
Execution code.
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.
std::shared_ptr< const Mantid::Geometry::IDetector > IDetector_const_sptr
Shared pointer to IDetector (const version)
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.
std::vector< double > MantidVec
typedef for the data storage used in Mantid matrix workspaces
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
@ Input
An input workspace.
@ Output
An output workspace.