35bool compareMomentum(
const std::array<double, 4> &v1,
const std::array<double, 4> &v2) {
return (v1[3] < v2[3]); }
45 : m_normWS(), m_inputWS(), m_hmin(0.0f), m_hmax(0.0f), m_kmin(0.0f), m_kmax(0.0f), m_lmin(0.0f), m_lmax(0.0f),
46 m_hIntegrated(true), m_kIntegrated(true), m_lIntegrated(true), m_rubw(3, 3), m_kiMin(0.0), m_kiMax(
EMPTY_DBL()),
47 m_hIdx(-1), m_kIdx(-1), m_lIdx(-1), m_hX(), m_kX(), m_lX(), m_samplePos(), m_beamDir() {}
57 return "Calculate normalization for an MDEvent workspace for single crystal "
69 "An input MDWorkspace.");
74 for (
size_t i = 0; i < dimChars.size(); i++) {
77 std::string propName =
"AlignedDim" + dim;
81 "Enter it as a comma-separated list of values with the format: "
82 "'name,minimum,maximum,number_of_bins'. Leave blank for NONE.");
85 auto fluxValidator = std::make_shared<CompositeValidator>();
89 auto solidAngleValidator = fluxValidator->
clone();
92 "An input workspace containing momentum dependent flux.");
95 "An input workspace containing momentum integrated vanadium "
96 "(a measure of the solid angle).");
99 "If set to true, the algorithm does "
100 "not check history if the workspace was modified since the"
101 "ConvertToMD algorithm was run, and assume that the elastic "
106 "An input MDHistoWorkspace used to accumulate normalization "
107 "from multiple MDEventWorkspaces. "
108 "If unspecified a blank MDHistoWorkspace will be created.");
112 "An input MDHistoWorkspace used to accumulate data from "
113 "multiple MDEventWorkspaces. If "
114 "unspecified a blank MDHistoWorkspace will be created.");
117 "A name for the output data MDHistoWorkspace.");
119 "A name for the output normalization MDHistoWorkspace.");
130 setProperty<Workspace_sptr>(
"OutputWorkspace", outputWS);
137 for (uint16_t expInfoIndex = 0; expInfoIndex <
m_numExptInfos; expInfoIndex++) {
140 bool skipNormalization =
false;
145 if (!skipNormalization) {
148 g_log.
warning(
"Binning limits are outside the limits of the MDWorkspace. "
149 "Not applying normalization.");
162 throw std::invalid_argument(
"Invalid energy transfer mode. Algorithm "
163 "currently only supports elastic data.");
167 m_hmin = hdim->getMinimum();
168 m_kmin = kdim->getMinimum();
169 m_lmin = ldim->getMinimum();
170 m_hmax = hdim->getMaximum();
171 m_kmax = kdim->getMaximum();
172 m_lmax = ldim->getMaximum();
174 const auto &exptInfoZero = *(
m_inputWS->getExperimentInfo(0));
175 auto source = exptInfoZero.getInstrument()->getSource();
176 auto sample = exptInfoZero.getInstrument()->getSample();
177 if (source ==
nullptr || sample ==
nullptr) {
179 "Instrument not sufficiently defined: failed to get source and/or "
192 const size_t nalgs =
history.size();
193 const auto &lastAlgorithm =
history.lastAlgorithm();
196 if (lastAlgorithm->name() ==
"ConvertToMD") {
197 emode = lastAlgorithm->getPropertyValue(
"dEAnalysisMode");
198 }
else if ((lastAlgorithm->name() ==
"Load" ||
history.lastAlgorithm()->name() ==
"LoadMD") &&
199 history.getAlgorithmHistory(nalgs - 2)->name() ==
"ConvertToMD") {
202 for (
auto &hist : histvec) {
203 if (hist->name() ==
"dEAnalysisMode") {
204 emode = hist->value();
209 throw std::invalid_argument(
"The last algorithm in the history of the "
210 "input workspace is not ConvertToMD");
223 binMD->setPropertyValue(
"AxisAligned",
"1");
224 for (
auto prop : props) {
225 const auto &propName = prop->name();
226 if (propName !=
"FluxWorkspace" && propName !=
"SolidAngleWorkspace" &&
227 propName !=
"TemporaryNormalizationWorkspace" && propName !=
"OutputNormalizationWorkspace" &&
228 propName !=
"SkipSafetyCheck") {
229 binMD->setPropertyValue(propName, prop->value());
232 binMD->executeAsChildAlg();
234 return std::dynamic_pointer_cast<MDHistoWorkspace>(outputWS);
243 std::shared_ptr<IMDHistoWorkspace>
tmp = this->
getProperty(
"TemporaryNormalizationWorkspace");
244 m_normWS = std::dynamic_pointer_cast<MDHistoWorkspace>(
tmp);
262 const auto ¤tRun =
m_inputWS->getExperimentInfo(expInfoIndex)->run();
264 std::vector<coord_t> otherDimValues;
265 for (
size_t i = 3; i <
m_inputWS->getNumDims(); i++) {
266 const auto dimension =
m_inputWS->getDimension(i);
267 auto dimMin =
static_cast<float>(dimension->getMinimum());
268 auto dimMax =
static_cast<float>(dimension->getMaximum());
271 auto value =
static_cast<coord_t>(dimProp->firstValue());
272 otherDimValues.emplace_back(
value);
275 if (value < dimMin || value > dimMax) {
276 skipNormalization =
true;
280 return otherDimValues;
293 bool &skipNormalization) {
298 const size_t nrm1 = affineMat.
numRows() - 1;
299 const size_t ncm1 = affineMat.
numCols() - 1;
300 for (
size_t row = 0; row < nrm1; row++)
302 const auto dimen =
m_normWS->getDimension(row);
303 const auto dimMin(dimen->getMinimum()), dimMax(dimen->getMaximum());
304 if (affineMat[row][0] == 1.) {
310 skipNormalization =
true;
313 if (affineMat[row][1] == 1.) {
319 skipNormalization =
true;
322 if (affineMat[row][2] == 1.) {
328 skipNormalization =
true;
332 for (
size_t col = 3; col < ncm1; col++)
334 if (affineMat[row][col] == 1.) {
335 double val = otherDimValues.at(col - 3);
336 if (val > dimMax || val < dimMin) {
337 skipNormalization =
true;
352 m_hX.resize(hDim.getNBoundaries());
353 for (
size_t i = 0; i <
m_hX.size(); ++i) {
354 m_hX[i] = hDim.getX(i);
359 m_kX.resize(kDim.getNBoundaries());
360 for (
size_t i = 0; i <
m_kX.size(); ++i) {
361 m_kX[i] = kDim.getX(i);
366 m_lX.resize(lDim.getNBoundaries());
367 for (
size_t i = 0; i <
m_lX.size(); ++i) {
368 m_lX[i] = lDim.getX(i);
386 const auto ¤tExptInfo = *(
m_inputWS->getExperimentInfo(expInfoIndex));
390 throw std::runtime_error(
"Wokspace does not contain a log entry for the RUBW matrix."
394 m_rubw = currentExptInfo.run().getGoniometerMatrix() * rubwValue;
397 const double protonCharge = currentExptInfo.run().getProtonCharge();
399 const auto &spectrumInfo = currentExptInfo.spectrumInfo();
402 const auto ndets =
static_cast<int64_t
>(spectrumInfo.size());
403 const detid2index_map fluxDetToIdx = integrFlux->getDetectorIDToWorkspaceIndexMap();
404 const detid2index_map solidAngDetToIdx = solidAngleWS->getDetectorIDToWorkspaceIndexMap();
406 const size_t vmdDims = 4;
407 std::vector<std::atomic<signal_t>> signalArray(
m_normWS->getNPoints());
408 std::vector<std::array<double, 4>> intersections;
409 std::vector<double> xValues, yValues;
410 std::vector<coord_t> pos, posNew;
413 std::make_unique<API::Progress>(
this, 0.3 + progStep * expInfoIndex, 0.3 + progStep * (expInfoIndex + 1.), ndets);
416for (int64_t i = 0; i < ndets; i++) {
419 if (!spectrumInfo.hasDetectors(i) || spectrumInfo.isMonitor(i) || spectrumInfo.isMasked(i)) {
422 const auto &detector = spectrumInfo.detector(i);
424 double phi = detector.getPhi();
426 const auto detID = detector.getID();
430 if (intersections.empty())
434 size_t wsIdx = fluxDetToIdx.find(detID)->second;
436 double solid = solidAngleWS->y(solidAngDetToIdx.find(detID)->second)[0] * protonCharge;
440 auto intersectionsBegin = intersections.begin();
442 xValues.resize(intersections.size());
443 yValues.resize(intersections.size());
444 auto x = xValues.begin();
445 for (
auto it = intersectionsBegin; it != intersections.end(); ++it, ++
x) {
455 pos.resize(vmdDims + otherValues.size());
456 std::copy(otherValues.begin(), otherValues.end(), pos.begin() + vmdDims - 1);
457 pos.emplace_back(1.f);
459 for (
auto it = intersectionsBegin + 1; it != intersections.end(); ++it) {
460 const auto &curIntSec = *it;
461 const auto &prevIntSec = *(it - 1);
463 double delta = curIntSec[3] - prevIntSec[3];
468 std::transform(curIntSec.begin(), curIntSec.begin() + vmdDims - 1, prevIntSec.begin(), pos.begin(),
469 [](
const double lhs,
const double rhs) { return static_cast<coord_t>(0.5 * (lhs + rhs)); });
471 size_t linIndex =
m_normWS->getLinearIndexAtCoord(posNew.data());
472 if (linIndex ==
size_t(-1))
476 auto k =
static_cast<size_t>(std::distance(intersectionsBegin, it));
478 signal_t signal = (yValues[k] - yValues[k - 1]) * solid;
487 std::transform(signalArray.cbegin(), signalArray.cend(),
m_normWS->getSignalArray(),
m_normWS->mutableSignalArray(),
488 [](
const std::atomic<signal_t> &a,
const signal_t &b) { return a + b; });
490 std::copy(signalArray.cbegin(), signalArray.cend(),
m_normWS->mutableSignalArray());
504 std::vector<double> &yValues)
const {
505 assert(xValues.size() == yValues.size());
508 const auto &xData = integrFlux.
x(sp);
509 const double xStart = xData.front();
510 const double xEnd = xData.back();
515 const auto &yData = integrFlux.
y(sp);
516 size_t spSize = yData.size();
518 const double yMin = 0.0;
519 const double yMax = yData.back();
521 size_t nData = xValues.size();
523 if (xValues[nData - 1] < xStart) {
524 std::fill(yValues.begin(), yValues.end(), yMin);
529 if (xValues[0] > xEnd) {
530 std::fill(yValues.begin(), yValues.end(), yMax);
536 while (i < nData - 1 && xValues[i] < xStart) {
541 for (; i < nData; i++) {
543 if (j >= spSize - 1) {
546 double xi = xValues[i];
547 while (j < spSize - 1 && xi > xData[j])
550 if (xi == xData[j]) {
551 yValues[i] = yData[j];
552 }
else if (j == spSize - 1) {
557 double x0 = xData[j - 1];
558 double x1 = xData[j];
559 double y0 = yData[j - 1];
560 double y1 = yData[j];
561 yValues[i] = y0 + (y1 - y0) * (xi - x0) / (x1 - x0);
580 V3D q(-sin(theta) * cos(phi), -sin(theta) * sin(phi), 1. - cos(theta));
592 auto hNBins =
m_hX.size();
593 auto kNBins =
m_kX.size();
594 auto lNBins =
m_lX.size();
595 intersections.clear();
596 intersections.reserve(hNBins + kNBins + lNBins + 8);
599 if (
fabs(hStart - hEnd) > eps) {
601 double fk = (kEnd - kStart) / (hEnd - hStart);
602 double fl = (lEnd - lStart) / (hEnd - hStart);
604 for (
size_t i = 0; i < hNBins; i++) {
606 if ((hi >=
m_hmin) && (hi <=
m_hmax) && ((hStart - hi) * (hEnd - hi) < 0)) {
610 double ki = fk * (hi - hStart) + kStart;
611 double li = fl * (hi - hStart) + lStart;
613 double momi = fmom * (hi - hStart) +
m_kiMin;
614 intersections.push_back({{hi, ki, li, momi}});
623 double khmin = fk * (
m_hmin - hStart) + kStart;
624 double lhmin = fl * (
m_hmin - hStart) + lStart;
626 intersections.push_back({{
m_hmin, khmin, lhmin, momhMin}});
632 double khmax = fk * (
m_hmax - hStart) + kStart;
633 double lhmax = fl * (
m_hmax - hStart) + lStart;
635 intersections.push_back({{
m_hmax, khmax, lhmax, momhMax}});
641 if (
fabs(kStart - kEnd) > eps) {
643 double fh = (hEnd - hStart) / (kEnd - kStart);
644 double fl = (lEnd - lStart) / (kEnd - kStart);
646 for (
size_t i = 0; i < kNBins; i++) {
648 if ((ki >=
m_kmin) && (ki <=
m_kmax) && ((kStart - ki) * (kEnd - ki) < 0)) {
651 double hi = fh * (ki - kStart) + hStart;
652 double li = fl * (ki - kStart) + lStart;
654 double momi = fmom * (ki - kStart) +
m_kiMin;
655 intersections.push_back({{hi, ki, li, momi}});
664 double hkmin = fh * (
m_kmin - kStart) + hStart;
665 double lkmin = fl * (
m_kmin - kStart) + lStart;
667 intersections.push_back({{hkmin,
m_kmin, lkmin, momkMin}});
673 double hkmax = fh * (
m_kmax - kStart) + hStart;
674 double lkmax = fl * (
m_kmax - kStart) + lStart;
676 intersections.push_back({{hkmax,
m_kmax, lkmax, momkMax}});
682 if (
fabs(lStart - lEnd) > eps) {
684 double fh = (hEnd - hStart) / (lEnd - lStart);
685 double fk = (kEnd - kStart) / (lEnd - lStart);
687 for (
size_t i = 0; i < lNBins; i++) {
689 if ((li >=
m_lmin) && (li <=
m_lmax) && ((lStart - li) * (lEnd - li) < 0)) {
692 double hi = fh * (li - lStart) + hStart;
693 double ki = fk * (li - lStart) + kStart;
695 double momi = fmom * (li - lStart) +
m_kiMin;
696 intersections.push_back({{hi, ki, li, momi}});
705 double hlmin = fh * (
m_lmin - lStart) + hStart;
706 double klmin = fk * (
m_lmin - lStart) + kStart;
708 intersections.push_back({{hlmin, klmin,
m_lmin, momlMin}});
714 double hlmax = fh * (
m_lmax - lStart) + hStart;
715 double klmax = fk * (
m_lmax - lStart) + kStart;
717 intersections.push_back({{hlmax, klmax,
m_lmax, momlMax}});
725 intersections.push_back({{hStart, kStart, lStart,
m_kiMin}});
729 intersections.push_back({{hEnd, kEnd, lEnd,
m_kiMax}});
733 std::stable_sort(intersections.begin(), intersections.end(), compareMomentum);
#define DECLARE_ALGORITHM(classname)
const std::vector< double > & rhs
double value
The value of the point.
#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...
#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::vector< history_type > history
history information
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
virtual std::shared_ptr< Algorithm > createChildAlgorithm(const std::string &name, const double startProgress=-1., const double endProgress=-1., const bool enableLogging=true, const int &version=-1)
Create a Child Algorithm.
const std::vector< Kernel::Property * > & getProperties() const override
Get the list of managed properties.
A validator which provides a TENTATIVE check that a workspace contains common bins in each spectrum.
Kernel::IValidator_sptr clone() const override
Clone the current state.
A validator which checks that a workspace has a valid instrument.
Base MatrixWorkspace Abstract Class.
const HistogramData::HistogramX & x(const size_t index) const
const HistogramData::HistogramY & y(const size_t index) const
A property class for workspaces.
A validator which checks that the unit of the workspace referred to by a WorkspaceProperty is the exp...
std::unique_ptr< MDHistoWorkspace > clone() const
Returns a clone of the workspace.
Exception for errors associated with the instrument definition.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void warning(const std::string &msg)
Logs at warning level.
T Invert()
LU inversion routine.
void multiplyPoint(const std::vector< T > &in, std::vector< T > &out) const
Multiply M*Vec.
size_t numRows() const
Return the number of rows in the matrix.
size_t numCols() const
Return the number of columns in the matrix.
The concrete, templated class for properties.
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
A specialised Property class for holding a series of time-value pairs.
constexpr double X() const noexcept
Get x.
constexpr double Y() const noexcept
Get y.
constexpr double Z() const noexcept
Get z.
MDNormSCD : Generate MD normalization for single crystal diffraction.
const std::string category() const override
Algorithm's category for identification.
bool m_accumulate
internal flag to accumulate to an existing workspace
DataObjects::MDHistoWorkspace_sptr m_normWS
Normalization workspace.
std::vector< double > m_kX
bool m_hIntegrated
flag for integrated h,k,l dimensions
void cacheDimensionXValues()
Stores the X values from each H,K,L dimension as member variables.
coord_t m_hmin
limits for h,k,l dimensions
void exec() override
Execute the algorithm.
std::string inputEnergyMode() const
Currently looks for the ConvertToMD algorithm in the history.
Mantid::Kernel::DblMatrix m_rubw
(2*PiRUBW)^-1
DataObjects::MDHistoWorkspace_sptr binInputWS()
Runs the BinMD algorithm on the input to provide the output workspace All slicing algorithm propertie...
std::vector< double > m_hX
cached X values along dimensions h,k,l
const std::string name() const override
Algorithm's name for use in the GUI and help.
int version() const override
Algorithm's version for identification.
void calculateNormalization(const std::vector< coord_t > &otherValues, const Kernel::Matrix< coord_t > &affineTrans, uint16_t expInfoIndex)
Computed the normalization for the input workspace.
double m_kiMin
limits for momentum
std::vector< double > m_lX
uint16_t m_numExptInfos
number of experiment infos
size_t m_hIdx
index of h,k,l dimensions in the output workspaces
Kernel::V3D m_samplePos
Sample position.
void init() override
Initialize the algorithm's properties.
Kernel::Matrix< coord_t > findIntergratedDimensions(const std::vector< coord_t > &otherDimValues, bool &skipNormalization)
Checks the normalization workspace against the indices of the original dimensions.
API::IMDEventWorkspace_sptr m_inputWS
Input workspace.
void calculateIntersections(std::vector< std::array< double, 4 > > &intersections, const double theta, const double phi)
Calculate the points of intersection for the given detector with cuboid surrounding the detector posi...
const std::string summary() const override
Algorithm's summary for use in the GUI and help.
std::vector< coord_t > getValuesFromOtherDimensions(bool &skipNormalization, uint16_t expInfoIndex=0) const
Retrieve logged values from non-HKL dimensions.
Kernel::V3D m_beamDir
Beam direction.
void createNormalizationWS(const DataObjects::MDHistoWorkspace &dataWS)
Create & cached the normalization workspace.
void cacheInputs()
Set up starting values for cached variables.
std::string convention
ki-kf for Inelastic convention; kf-ki for Crystallography convention
void calcIntegralsForIntersections(const std::vector< double > &xValues, const API::MatrixWorkspace &integrFlux, size_t sp, std::vector< double > &yValues) const
Linearly interpolate between the points in integrFlux at xValues and save the results in yValues.
static std::string getDimensionChars()
std::shared_ptr< Workspace > Workspace_sptr
shared pointer to Mantid::API::Workspace
std::shared_ptr< const MatrixWorkspace > MatrixWorkspace_const_sptr
shared pointer to the matrix workspace base class (const version)
@ NoNormalization
Don't normalize = return raw counts.
std::shared_ptr< MDHistoWorkspace > MDHistoWorkspace_sptr
A shared pointer to a MDHistoWorkspace.
std::string toString(const T &value)
Convert a number to a string.
std::vector< PropertyHistory_sptr > PropertyHistories
std::enable_if< std::is_pointer< Arg >::value, bool >::type threadSafe(Arg workspace)
Thread-safety check Checks the workspace to ensure it is suitable for multithreaded access.
MANTID_KERNEL_DLL V3D normalize(V3D v)
Normalizes a V3D.
void AtomicOp(std::atomic< T > &f, T d, BinaryOp op)
Uses std::compare_exchange_weak to update the atomic value f = op(f, d) Used to improve parallel scal...
float coord_t
Typedef for the data type to use for coordinate axes in MD objects such as MDBox, MDEventWorkspace,...
std::unordered_map< detid_t, size_t > detid2index_map
Map with key = detector ID, value = workspace index.
double signal_t
Typedef for the signal recorded in a MDBox, etc.
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
Describes the direction (within an algorithm) of a Property.
@ Input
An input workspace.
@ Output
An output workspace.