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_dEmin(0.f), m_dEmax(0.f), m_Ei(0.), m_ki(0.), m_kfmin(0.), m_kfmax(0.), m_hIntegrated(true),
47 m_kIntegrated(true), m_lIntegrated(true), m_dEIntegrated(true), m_rubw(3, 3), m_hIdx(-1), m_kIdx(-1), m_lIdx(-1),
48 m_eIdx(-1), m_hX(), m_kX(), m_lX(), m_eX(), m_samplePos(), m_beamDir() {}
58 return "Calculate normalization for an MDEvent workspace for single crystal "
59 "direct geometry inelastic measurement.";
70 "An input MDWorkspace.");
75 for (
size_t i = 0; i < dimChars.size(); i++) {
78 std::string propName =
"AlignedDim" + dim;
82 "Enter it as a comma-separated list of values with the format: "
83 "'name,minimum,maximum,number_of_bins'. Leave blank for NONE.");
86 auto solidAngleValidator = std::make_shared<CompositeValidator>();
92 "An input workspace containing integrated vanadium (a measure of the "
96 "If set to true, the algorithm does "
97 "not check history if the workspace was modified since the"
98 "ConvertToMD algorithm was run, and assume that the direct "
99 "geometry inelastic mode is used.");
103 "An input MDHistoWorkspace used to accumulate normalization "
104 "from multiple MDEventWorkspaces. If unspecified a blank "
105 "MDHistoWorkspace will be created.");
109 "An input MDHistoWorkspace used to accumulate data from "
110 "multiple MDEventWorkspaces. If unspecified a blank "
111 "MDHistoWorkspace will be created.");
114 "A name for the output data MDHistoWorkspace.");
116 "A name for the output normalization MDHistoWorkspace.");
128 setProperty<Workspace_sptr>(
"OutputWorkspace", outputWS);
135 for (uint16_t expInfoIndex = 0; expInfoIndex <
m_numExptInfos; expInfoIndex++) {
138 bool skipNormalization =
false;
143 if (!skipNormalization) {
146 g_log.
warning(
"Binning limits are outside the limits of the MDWorkspace. "
147 "Not applying normalization.");
154 outputWS->setDisplayNormalization(
m_inputWS->displayNormalizationHisto());
164 throw std::invalid_argument(
"Invalid energy transfer mode. Algorithm only "
165 "supports direct geometry spectrometers.");
170 m_hmin = hdim->getMinimum();
171 m_kmin = kdim->getMinimum();
172 m_lmin = ldim->getMinimum();
174 m_hmax = hdim->getMaximum();
175 m_kmax = kdim->getMaximum();
176 m_lmax = ldim->getMaximum();
179 const auto &exptInfoZero = *(
m_inputWS->getExperimentInfo(0));
180 auto source = exptInfoZero.getInstrument()->getSource();
181 auto sample = exptInfoZero.getInstrument()->getSample();
182 if (source ==
nullptr || sample ==
nullptr) {
184 "Instrument not sufficiently defined: failed to get source and/or "
190 double originaldEmin = exptInfoZero.run().getBinBoundaries().front();
191 double originaldEmax = exptInfoZero.run().getBinBoundaries().back();
192 if (exptInfoZero.run().hasProperty(
"Ei")) {
193 m_Ei = exptInfoZero.run().getPropertyValueAsType<
double>(
"Ei");
195 throw std::invalid_argument(
"Ei stored in the workspace is not positive");
198 throw std::invalid_argument(
"Could not find Ei value in the workspace.");
201 if (
m_Ei - originaldEmin < eps) {
202 originaldEmin =
m_Ei - eps;
204 if (
m_Ei - originaldEmax < eps) {
205 originaldEmax =
m_Ei - 1e-7;
207 if (originaldEmin == originaldEmax) {
208 throw std::runtime_error(
"The limits of the original workspace used in "
209 "ConvertToMD are incorrect");
214 m_kfmin = std::sqrt(energyToK * (
m_Ei - originaldEmin));
215 m_kfmax = std::sqrt(energyToK * (
m_Ei - originaldEmax));
223 const auto &hist =
m_inputWS->getHistory();
224 const size_t nalgs = hist.size();
225 const auto &lastAlgHist = hist.getAlgorithmHistory(nalgs - 1);
226 const auto &penultimateAlgHist = hist.getAlgorithmHistory(nalgs - 2);
229 if (lastAlgHist->name() ==
"ConvertToMD") {
230 emode = lastAlgHist->getPropertyValue(
"dEAnalysisMode");
231 }
else if ((lastAlgHist->name() ==
"Load" || lastAlgHist->name() ==
"LoadMD") &&
232 penultimateAlgHist->name() ==
"ConvertToMD") {
234 emode = penultimateAlgHist->getPropertyValue(
"dEAnalysisMode");
236 throw std::invalid_argument(
"The last algorithm in the history of the "
237 "input workspace is not ConvertToMD");
250 binMD->setPropertyValue(
"AxisAligned",
"1");
251 for (
auto prop : props) {
252 const auto &propName = prop->name();
253 if (propName !=
"SolidAngleWorkspace" && propName !=
"TemporaryNormalizationWorkspace" &&
254 propName !=
"OutputNormalizationWorkspace" && propName !=
"SkipSafetyCheck") {
255 binMD->setPropertyValue(propName, prop->value());
258 binMD->executeAsChildAlg();
260 return std::dynamic_pointer_cast<MDHistoWorkspace>(outputWS);
269 std::shared_ptr<IMDHistoWorkspace>
tmp = this->
getProperty(
"TemporaryNormalizationWorkspace");
270 m_normWS = std::dynamic_pointer_cast<MDHistoWorkspace>(
tmp);
288 uint16_t expInfoIndex)
const {
289 const auto ¤tRun =
m_inputWS->getExperimentInfo(expInfoIndex)->run();
291 std::vector<coord_t> otherDimValues;
292 for (
size_t i = 4; i <
m_inputWS->getNumDims(); i++) {
293 const auto dimension =
m_inputWS->getDimension(i);
294 auto dimMin =
static_cast<float>(dimension->getMinimum());
295 auto dimMax =
static_cast<float>(dimension->getMaximum());
298 auto value =
static_cast<coord_t>(dimProp->firstValue());
299 otherDimValues.emplace_back(
value);
302 if (value < dimMin || value > dimMax) {
303 skipNormalization =
true;
307 return otherDimValues;
320 bool &skipNormalization) {
325 const size_t nrm1 = affineMat.
numRows() - 1;
326 const size_t ncm1 = affineMat.
numCols() - 1;
327 for (
size_t row = 0; row < nrm1; row++)
329 const auto dimen =
m_normWS->getDimension(row);
330 const auto dimMin(dimen->getMinimum()), dimMax(dimen->getMaximum());
331 if (affineMat[row][0] == 1.) {
337 skipNormalization =
true;
340 if (affineMat[row][1] == 1.) {
346 skipNormalization =
true;
349 if (affineMat[row][2] == 1.) {
355 skipNormalization =
true;
359 if (affineMat[row][3] == 1.) {
365 skipNormalization =
true;
368 for (
size_t col = 4; col < ncm1; col++)
370 if (affineMat[row][col] == 1.) {
371 double val = otherDimValues.at(col - 3);
372 if (val > dimMax || val < dimMin) {
373 skipNormalization =
true;
391 m_hX.resize(hDim.getNBoundaries());
392 for (
size_t i = 0; i <
m_hX.size(); ++i) {
393 m_hX[i] = hDim.getX(i);
398 m_kX.resize(kDim.getNBoundaries());
399 for (
size_t i = 0; i <
m_kX.size(); ++i) {
400 m_kX[i] = kDim.getX(i);
405 m_lX.resize(lDim.getNBoundaries());
406 for (
size_t i = 0; i <
m_lX.size(); ++i) {
407 m_lX[i] = lDim.getX(i);
413 m_eX.resize(eDim.getNBoundaries());
414 for (
size_t i = 0; i <
m_eX.size(); ++i) {
415 double temp =
m_Ei - eDim.getX(i);
416 temp = std::max(temp, 0.);
417 m_eX[i] = std::sqrt(energyToK * temp);
433 const auto ¤tExptInfo = *(
m_inputWS->getExperimentInfo(expInfoIndex));
437 throw std::runtime_error(
"Wokspace does not contain a log entry for the RUBW matrix."
441 m_rubw = currentExptInfo.run().getGoniometerMatrix() * rubwValue;
444 const double protonCharge = currentExptInfo.run().getProtonCharge();
446 const auto &spectrumInfo = currentExptInfo.spectrumInfo();
449 const auto ndets =
static_cast<int64_t
>(spectrumInfo.size());
453 if (solidAngleWS !=
nullptr) {
455 solidAngDetToIdx = solidAngleWS->getDetectorIDToWorkspaceIndexMap();
458 const size_t vmdDims = 4;
459 std::vector<std::atomic<signal_t>> signalArray(
m_normWS->getNPoints());
460 std::vector<std::array<double, 4>> intersections;
461 std::vector<coord_t> pos, posNew;
464 std::make_unique<API::Progress>(
this, 0.3 + progStep * expInfoIndex, 0.3 + progStep * (expInfoIndex + 1.), ndets);
466PRAGMA_OMP(parallel
for private(intersections, pos, posNew))
467for (int64_t i = 0; i < ndets; i++) {
470 if (!spectrumInfo.hasDetectors(i) || spectrumInfo.isMonitor(i) || spectrumInfo.isMasked(i)) {
473 const auto &detector = spectrumInfo.detector(i);
475 double phi = detector.getPhi();
477 const auto detID = detector.getID();
481 if (intersections.empty())
485 double solid = protonCharge;
487 solid = solidAngleWS->y(solidAngDetToIdx.find(detID)->second)[0] * protonCharge;
491 pos.resize(vmdDims + otherValues.size() + 1);
492 std::copy(otherValues.begin(), otherValues.end(), pos.begin() + vmdDims);
493 pos.emplace_back(1.f);
494 auto intersectionsBegin = intersections.begin();
495 for (
auto it = intersectionsBegin + 1; it != intersections.end(); ++it) {
496 const auto &curIntSec = *it;
497 const auto &prevIntSec = *(it - 1);
499 double delta = (curIntSec[3] * curIntSec[3] - prevIntSec[3] * prevIntSec[3]) / energyToK;
504 std::transform(curIntSec.data(), curIntSec.data() + vmdDims, prevIntSec.data(), pos.begin(),
505 [](
const double rhs,
const double lhs) { return static_cast<coord_t>(0.5 * (rhs + lhs)); });
508 pos[3] =
static_cast<coord_t>(
m_Ei - pos[3] * pos[3] / energyToK);
510 size_t linIndex =
m_normWS->getLinearIndexAtCoord(posNew.data());
511 if (linIndex ==
size_t(-1))
516 double signal = solid *
delta;
525 std::transform(signalArray.cbegin(), signalArray.cend(),
m_normWS->getSignalArray(),
m_normWS->mutableSignalArray(),
526 [](
const std::atomic<signal_t> &a,
const signal_t &b) { return a + b; });
528 std::copy(signalArray.cbegin(), signalArray.cend(),
m_normWS->mutableSignalArray());
542 V3D qout(sin(theta) * cos(phi), sin(theta) * sin(phi), cos(theta)), qin(0., 0.,
m_ki);
550 double hStart = qin.
X() - qout.X() *
m_kfmin, hEnd = qin.
X() - qout.X() *
m_kfmax;
551 double kStart = qin.
Y() - qout.Y() *
m_kfmin, kEnd = qin.
Y() - qout.Y() *
m_kfmax;
552 double lStart = qin.
Z() - qout.Z() *
m_kfmin, lEnd = qin.
Z() - qout.Z() *
m_kfmax;
554 auto hNBins =
m_hX.size();
555 auto kNBins =
m_kX.size();
556 auto lNBins =
m_lX.size();
557 auto eNBins =
m_eX.size();
558 intersections.clear();
559 intersections.reserve(hNBins + kNBins + lNBins + eNBins + 8);
562 if (
fabs(hStart - hEnd) > eps) {
564 double fk = (kEnd - kStart) / (hEnd - hStart);
565 double fl = (lEnd - lStart) / (hEnd - hStart);
567 for (
size_t i = 0; i < hNBins; i++) {
569 if ((hi >=
m_hmin) && (hi <=
m_hmax) && ((hStart - hi) * (hEnd - hi) < 0)) {
573 double ki = fk * (hi - hStart) + kStart;
574 double li = fl * (hi - hStart) + lStart;
576 double momi = fmom * (hi - hStart) +
m_kfmin;
577 intersections.push_back({{hi, ki, li, momi}});
586 double khmin = fk * (
m_hmin - hStart) + kStart;
587 double lhmin = fl * (
m_hmin - hStart) + lStart;
589 intersections.push_back({{
m_hmin, khmin, lhmin, momhMin}});
595 double khmax = fk * (
m_hmax - hStart) + kStart;
596 double lhmax = fl * (
m_hmax - hStart) + lStart;
598 intersections.push_back({{
m_hmax, khmax, lhmax, momhMax}});
604 if (
fabs(kStart - kEnd) > eps) {
606 double fh = (hEnd - hStart) / (kEnd - kStart);
607 double fl = (lEnd - lStart) / (kEnd - kStart);
609 for (
size_t i = 0; i < kNBins; i++) {
611 if ((ki >=
m_kmin) && (ki <=
m_kmax) && ((kStart - ki) * (kEnd - ki) < 0)) {
615 double hi = fh * (ki - kStart) + hStart;
616 double li = fl * (ki - kStart) + lStart;
618 double momi = fmom * (ki - kStart) +
m_kfmin;
619 intersections.push_back({{hi, ki, li, momi}});
627 double hkmin = fh * (
m_kmin - kStart) + hStart;
628 double lkmin = fl * (
m_kmin - kStart) + lStart;
630 intersections.push_back({{hkmin,
m_kmin, lkmin, momkMin}});
636 double hkmax = fh * (
m_kmax - kStart) + hStart;
637 double lkmax = fl * (
m_kmax - kStart) + lStart;
639 intersections.push_back({{hkmax,
m_kmax, lkmax, momkMax}});
645 if (
fabs(lStart - lEnd) > eps) {
647 double fh = (hEnd - hStart) / (lEnd - lStart);
648 double fk = (kEnd - kStart) / (lEnd - lStart);
650 for (
size_t i = 0; i < lNBins; i++) {
652 if ((li >=
m_lmin) && (li <=
m_lmax) && ((lStart - li) * (lEnd - li) < 0)) {
653 double hi = fh * (li - lStart) + hStart;
654 double ki = fk * (li - lStart) + kStart;
656 double momi = fmom * (li - lStart) +
m_kfmin;
657 intersections.push_back({{hi, ki, li, momi}});
665 double hlmin = fh * (
m_lmin - lStart) + hStart;
666 double klmin = fk * (
m_lmin - lStart) + kStart;
668 intersections.push_back({{hlmin, klmin,
m_lmin, momlMin}});
674 double hlmax = fh * (
m_lmax - lStart) + hStart;
675 double klmax = fk * (
m_lmax - lStart) + kStart;
677 intersections.push_back({{hlmax, klmax,
m_lmax, momlMax}});
684 for (
size_t i = 0; i < eNBins; i++) {
685 double kfi =
m_eX[i];
687 double h = qin.
X() - qout.X() * kfi;
688 double k = qin.
Y() - qout.Y() * kfi;
689 double l = qin.
Z() - qout.Z() * kfi;
691 intersections.push_back({{h, k, l, kfi}});
700 intersections.push_back({{hStart, kStart, lStart,
m_kfmin}});
704 intersections.push_back({{hEnd, kEnd, lEnd,
m_kfmax}});
708 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.
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.
A validator which checks that a workspace has a valid instrument.
A property class for workspaces.
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.
DataObjects::MDHistoWorkspace_sptr binInputWS()
Runs the BinMD algorithm on the input to provide the output workspace All slicing algorithm propertie...
Mantid::Kernel::DblMatrix m_rubw
(2*PiRUBW)^-1
API::IMDEventWorkspace_sptr m_inputWS
Input workspace.
double m_Ei
cached values for incident energy and momentum, final momentum min/max
int version() const override
Algorithm's version for identification.
Kernel::V3D m_samplePos
Sample position.
Kernel::V3D m_beamDir
Beam direction.
void calculateNormalization(const std::vector< coord_t > &otherValues, const Kernel::Matrix< coord_t > &affineTrans, uint16_t expInfoIndex)
Computed the normalization for the input workspace.
void createNormalizationWS(const DataObjects::MDHistoWorkspace &dataWS)
Create & cached the normalization workspace.
DataObjects::MDHistoWorkspace_sptr m_normWS
Normalization workspace.
std::vector< double > m_lX
bool m_hIntegrated
flag for integrated h,k,l, dE dimensions
Kernel::Matrix< coord_t > findIntergratedDimensions(const std::vector< coord_t > &otherDimValues, bool &skipNormalization)
Checks the normalization workspace against the indices of the original dimensions.
const std::string name() const override
Algorithm's name for use in the GUI and help.
coord_t m_hmin
limits for h,k,l, dE dimensions
std::vector< coord_t > getValuesFromOtherDimensions(bool &skipNormalization, uint16_t expInfoIndex=0) const
Retrieve logged values from non-HKL dimensions.
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...
void cacheInputs()
Set up starting values for cached variables.
void init() override
Initialize the algorithm's properties.
const std::string summary() const override
Algorithm's summary for use in the GUI and help.
uint16_t m_numExptInfos
number of experiment infos
std::vector< double > m_eX
std::vector< double > m_hX
cached X values along dimensions h,k,l. dE
const std::string category() const override
Algorithm's category for identification.
std::vector< double > m_kX
size_t m_hIdx
index of h,k,l, dE dimensions in the output workspaces
std::string convention
ki-kf for Inelastic convention; kf-ki for Crystallography convention
void cacheDimensionXValues()
Stores the X values from each H,K,L,E dimension as member variables Energy dimension is transformed t...
bool m_accumulate
internal flag to accumulate to an existing workspace
std::string inputEnergyMode() const
Currently looks for the ConvertToMD algorithm in the history.
void exec() override
Execute the algorithm.
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.
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...
static constexpr double NeutronMass
Mass of the neutron in kg.
static constexpr double h
Planck constant in J*s.
static constexpr double meV
1 meV in Joules.
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.
Describes the direction (within an algorithm) of a Property.
@ Input
An input workspace.
@ Output
An output workspace.