25#include <boost/lexical_cast.hpp>
44 : m_hmin(0.f), m_hmax(0.f), m_kmin(0.f), m_kmax(0.f), m_lmin(0.f), m_lmax(0.f), m_dEmin(0.f), m_dEmax(0.f),
45 m_Ei(0.), m_ki(0.), m_kfmin(0.), m_kfmax(0.), m_hIntegrated(false), m_kIntegrated(false), m_lIntegrated(false),
46 m_dEIntegrated(false), m_hX(), m_kX(), m_lX(), m_eX(), m_hIdx(-1), m_kIdx(-1), m_lIdx(-1), m_eIdx(-1),
47 m_rubw(3, 3), m_normWS() {}
60 return "Calculate the reciprocal space coverage for direct geometry "
72 m_hX.resize(hDim.getNBins());
73 for (
size_t i = 0; i <
m_hX.size(); ++i) {
74 m_hX[i] = hDim.getX(i);
77 m_kX.resize(kDim.getNBins());
78 for (
size_t i = 0; i <
m_kX.size(); ++i) {
79 m_kX[i] = kDim.getX(i);
82 m_lX.resize(lDim.getNBins());
83 for (
size_t i = 0; i <
m_lX.size(); ++i) {
84 m_lX[i] = lDim.getX(i);
88 m_eX.resize(eDim.getNBins());
89 for (
size_t i = 0; i <
m_eX.size(); ++i) {
90 m_eX[i] = std::sqrt(energyToK * (
m_Ei - eDim.getX(i)));
98 std::make_shared<InstrumentValidator>()),
99 "An input workspace.");
102 auto mustBe3D = std::make_shared<ArrayLengthValidator<double> >(3);
104 auto mustBePositive = std::make_shared<BoundedValidator<double>>();
105 mustBePositive->setLower(0.0);
107 std::vector<double> Q1(3, 0.), Q2(3, 0), Q3(3, 0);
112 "Q1 projection direction in the x,y,z format. Q1, Q2, Q3 "
113 "must not be coplanar");
115 "Q2 projection direction in the x,y,z format. Q1, Q2, Q3 "
116 "must not be coplanar");
118 "Q3 projection direction in the x,y,z format. Q1, Q2, Q3 "
119 "must not be coplanar");
122 "Incident energy. If set, will override Ei in the input workspace");
124 std::vector<std::string> options{
"Q1",
"Q2",
"Q3",
"DeltaE"};
126 for (
int i = 1; i <= 4; i++) {
127 std::string dim(
"Dimension");
128 dim += boost::lexical_cast<std::string>(i);
129 declareProperty(dim, options[i - 1], std::make_shared<StringListValidator>(options),
130 "Dimension to bin or integrate");
134 dim +
" step size. If empty the dimension will be "
135 "integrated between minimum and maximum values");
139 "A name for the output data MDHistoWorkspace.");
151 const auto &detectorInfo = inputWS->detectorInfo();
152 std::vector<double> tt, phi;
153 for (
size_t i = 0; i < detectorInfo.size(); ++i) {
154 if (!detectorInfo.isMasked(i) && !detectorInfo.isMonitor(i)) {
155 const auto &detector = detectorInfo.detector(i);
156 tt.emplace_back(detector.getTwoTheta(
V3D(0, 0, 0),
V3D(0, 0, 1)));
157 phi.emplace_back(detector.getPhi());
161 double ttmax = *(std::max_element(tt.begin(), tt.end()));
164 if (inputWS->run().hasProperty(
"Ei")) {
166 m_Ei = boost::lexical_cast<double>(eiprop->
value());
168 throw std::invalid_argument(
"Ei stored in the workspace is not positive");
171 throw std::invalid_argument(
"Could not find Ei in the workspace. Please enter a positive value");
178 double q1min(0.), q1max(0.), q2min(0.), q2max(0.), q3min(0.), q3max(0.);
179 double q1step(0.), q2step(0.), q3step(0.), dEstep(0.);
180 size_t q1NumBins = 1, q2NumBins = 1, q3NumBins = 1, dENumBins = 1;
181 for (
int i = 1; i <= 4; i++) {
182 std::string dim(
"Dimension");
183 dim += boost::lexical_cast<std::string>(i);
185 if (dimensioni ==
"Q1") {
186 affineMat[i - 1][0] = 1.;
192 if (dimensioni ==
"Q2") {
193 affineMat[i - 1][1] = 1.;
199 if (dimensioni ==
"Q3") {
200 affineMat[i - 1][2] = 1.;
206 if (dimensioni ==
"DeltaE") {
207 affineMat[i - 1][3] = 1.;
226 if (dEstep *
static_cast<double>(dENumBins) +
m_dEmin <
m_dEmax) {
236 throw std::invalid_argument(
"Please make sure each dimension is selected only once.");
253 double Qmax = QmaxTemp;
260 if (inputWS->run().getGoniometer().isDefined()) {
261 gon = inputWS->run().getGoniometerMatrix();
266 if (inputWS->sample().hasOrientedLattice()) {
267 UB = inputWS->sample().getOrientedLattice().getUB();
272 std::vector<double> Q1Basis =
getProperty(
"Q1Basis");
273 std::vector<double> Q2Basis =
getProperty(
"Q2Basis");
274 std::vector<double> Q3Basis =
getProperty(
"Q3Basis");
279 m_rubw = gon * UB * W * (2.0 * M_PI);
298 if (q1min >= q1max) {
299 throw std::invalid_argument(
"Q1max has to be greater than Q1min");
304 q1NumBins =
static_cast<size_t>((q1max - q1min) / q1step);
305 if (q1step *
static_cast<double>(q1NumBins) + q1min < q1max) {
307 q1max = q1step *
static_cast<double>(q1NumBins) + q1min;
317 if (q2min >= q2max) {
318 throw std::invalid_argument(
"Q2max has to be greater than Q2min");
323 q2NumBins =
static_cast<size_t>((q2max - q2min) / q2step);
324 if (q2step *
static_cast<double>(q2NumBins) + q2min < q2max) {
326 q2max = q2step *
static_cast<double>(q2NumBins) + q2min;
336 if (q3min >= q3max) {
337 throw std::invalid_argument(
"Q3max has to be greater than Q3min");
342 q3NumBins =
static_cast<size_t>((q3max - q3min) / q3step);
343 if (q3step *
static_cast<double>(q3NumBins) + q3min < q3max) {
345 q3max = q3step *
static_cast<double>(q3NumBins) + q3min;
350 std::vector<Mantid::Geometry::MDHistoDimension_sptr> binDimensions;
357 auto out1 = std::make_shared<MDHistoDimension>(
"Q1",
"Q1", frame1,
static_cast<coord_t>(q1min),
358 static_cast<coord_t>(q1max), q1NumBins);
359 auto out2 = std::make_shared<MDHistoDimension>(
"Q2",
"Q2", frame2,
static_cast<coord_t>(q2min),
360 static_cast<coord_t>(q2max), q2NumBins);
361 auto out3 = std::make_shared<MDHistoDimension>(
"Q3",
"Q3", frame3,
static_cast<coord_t>(q3min),
362 static_cast<coord_t>(q3max), q3NumBins);
363 auto out4 = std::make_shared<MDHistoDimension>(
"DeltaE",
"DeltaE", frame4,
static_cast<coord_t>(
m_dEmin),
366 for (
size_t row = 0; row <= 3; row++) {
367 if (affineMat[row][0] == 1.) {
368 binDimensions.emplace_back(out1);
370 if (affineMat[row][1] == 1.) {
371 binDimensions.emplace_back(out2);
373 if (affineMat[row][2] == 1.) {
374 binDimensions.emplace_back(out3);
376 if (affineMat[row][3] == 1.) {
377 binDimensions.emplace_back(out4);
381 m_normWS = std::make_shared<MDHistoWorkspace>(binDimensions);
387 const auto ndets =
static_cast<int64_t
>(tt.size());
390 for (int64_t i = 0; i < ndets; i++) {
393 if (intersections.empty())
395 auto intersectionsBegin = intersections.begin();
396 for (
auto it = intersectionsBegin + 1; it != intersections.end(); ++it) {
397 const auto &curIntSec = *it;
398 const auto &prevIntSec = *(it - 1);
400 double delta = curIntSec[3] - prevIntSec[3];
404 std::vector<coord_t> pos(4);
405 std::transform(curIntSec.getBareArray(), curIntSec.getBareArray() + 4, prevIntSec.getBareArray(), pos.begin(),
406 [](
const double lhs,
const double rhs) { return static_cast<coord_t>(0.5 * (lhs + rhs)); });
408 pos[3] =
static_cast<coord_t>(
m_Ei - pos[3] * pos[3] / energyToK);
410 std::vector<coord_t> posNew = affineMat * pos;
411 size_t linIndex =
m_normWS->getLinearIndexAtCoord(posNew.data());
412 if (linIndex ==
size_t(-1))
430 V3D qout(sin(theta) * cos(phi), sin(theta) * sin(phi), cos(theta)), qin(0., 0.,
m_ki);
437 double hStart = qin.
X() - qout.X() *
m_kfmin, hEnd = qin.
X() - qout.X() *
m_kfmax;
438 double kStart = qin.
Y() - qout.Y() *
m_kfmin, kEnd = qin.
Y() - qout.Y() *
m_kfmax;
439 double lStart = qin.
Z() - qout.Z() *
m_kfmin, lEnd = qin.
Z() - qout.Z() *
m_kfmax;
441 auto hNBins =
m_hX.size();
442 auto kNBins =
m_kX.size();
443 auto lNBins =
m_lX.size();
444 auto eNBins =
m_eX.size();
445 std::vector<Kernel::VMD> intersections;
446 intersections.reserve(hNBins + kNBins + lNBins + eNBins + 8);
449 if (
fabs(hStart - hEnd) > eps) {
451 double fk = (kEnd - kStart) / (hEnd - hStart);
452 double fl = (lEnd - lStart) / (hEnd - hStart);
454 for (
size_t i = 0; i < hNBins; i++) {
456 if ((hi >=
m_hmin) && (hi <=
m_hmax) && ((hStart - hi) * (hEnd - hi) < 0)) {
460 double ki = fk * (hi - hStart) + kStart;
461 double li = fl * (hi - hStart) + lStart;
463 double momi = fmom * (hi - hStart) +
m_kfmin;
464 intersections.emplace_back(hi, ki, li, momi);
473 double khmin = fk * (
m_hmin - hStart) + kStart;
474 double lhmin = fl * (
m_hmin - hStart) + lStart;
476 intersections.emplace_back(
m_hmin, khmin, lhmin, momhMin);
482 double khmax = fk * (
m_hmax - hStart) + kStart;
483 double lhmax = fl * (
m_hmax - hStart) + lStart;
485 intersections.emplace_back(
m_hmax, khmax, lhmax, momhMax);
491 if (
fabs(kStart - kEnd) > eps) {
493 double fh = (hEnd - hStart) / (kEnd - kStart);
494 double fl = (lEnd - lStart) / (kEnd - kStart);
496 for (
size_t i = 0; i < kNBins; i++) {
498 if ((ki >=
m_kmin) && (ki <=
m_kmax) && ((kStart - ki) * (kEnd - ki) < 0)) {
502 double hi = fh * (ki - kStart) + hStart;
503 double li = fl * (ki - kStart) + lStart;
505 double momi = fmom * (ki - kStart) +
m_kfmin;
506 intersections.emplace_back(hi, ki, li, momi);
514 double hkmin = fh * (
m_kmin - kStart) + hStart;
515 double lkmin = fl * (
m_kmin - kStart) + lStart;
517 intersections.emplace_back(hkmin,
m_kmin, lkmin, momkMin);
523 double hkmax = fh * (
m_kmax - kStart) + hStart;
524 double lkmax = fl * (
m_kmax - kStart) + lStart;
526 intersections.emplace_back(hkmax,
m_kmax, lkmax, momkMax);
532 if (
fabs(lStart - lEnd) > eps) {
534 double fh = (hEnd - hStart) / (lEnd - lStart);
535 double fk = (kEnd - kStart) / (lEnd - lStart);
537 for (
size_t i = 0; i < lNBins; i++) {
539 if ((li >=
m_lmin) && (li <=
m_lmax) && ((lStart - li) * (lEnd - li) < 0)) {
540 double hi = fh * (li - lStart) + hStart;
541 double ki = fk * (li - lStart) + kStart;
543 double momi = fmom * (li - lStart) +
m_kfmin;
544 intersections.emplace_back(hi, ki, li, momi);
552 double hlmin = fh * (
m_lmin - lStart) + hStart;
553 double klmin = fk * (
m_lmin - lStart) + kStart;
555 intersections.emplace_back(hlmin, klmin,
m_lmin, momlMin);
561 double hlmax = fh * (
m_lmax - lStart) + hStart;
562 double klmax = fk * (
m_lmax - lStart) + kStart;
564 intersections.emplace_back(hlmax, klmax,
m_lmax, momlMax);
571 for (
size_t i = 0; i < eNBins; i++) {
572 double kfi =
m_eX[i];
574 double h = qin.
X() - qout.X() * kfi;
575 double k = qin.
Y() - qout.Y() * kfi;
576 double l = qin.
Z() - qout.Z() * kfi;
578 intersections.emplace_back(h, k, l, kfi);
587 intersections.emplace_back(hStart, kStart, lStart,
m_kfmin);
591 intersections.emplace_back(hEnd, kEnd, lEnd,
m_kfmax);
595 using IterType = std::vector<Mantid::Kernel::VMD>::iterator;
596 std::stable_sort<IterType, bool (*)(const Mantid::Kernel::VMD &, const Mantid::Kernel::VMD &)>(
597 intersections.begin(), intersections.end(), compareMomentum);
599 return intersections;
#define DECLARE_ALGORITHM(classname)
const std::vector< double > & rhs
#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_CRITICAL(name)
#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 PARALLEL_FOR_IF(condition)
Empty definitions - to enable set your complier to enable openMP.
#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.
A property class for workspaces.
GeneralFrame : Any MDFrame that isn't related to momemtum transfer.
Class to implement UB matrix.
void setUB(const Kernel::DblMatrix &newUB)
Sets the UB matrix and recalculates lattice parameters.
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.
T determinant() const
Calculate the determinant.
T Invert()
LU inversion routine.
void setColumn(const size_t nCol, const std::vector< T > &newCol)
The concrete, templated class for properties.
Base class for properties.
virtual std::string value() const =0
Returns the value of the property as a string.
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.
constexpr double Z() const noexcept
Get z.
CalculateCoverageDGS : Calculate coverage for single crystal direct geometry scattering.
bool m_hIntegrated
flag for integrated h,k,l,dE dimensions
const std::string category() const override
Algorithm's category for identification.
Mantid::DataObjects::MDHistoWorkspace_sptr m_normWS
Normalization workspace (this is the coverage workspace)
const std::string name() const override
Algorithms name for identification.
coord_t m_hmin
limits for h,k,l,dE dimensions
void cacheDimensionXValues()
Stores the X values from each H,K,L dimension as member variables.
double m_Ei
cached values for incident energy and momentum, final momentum min/max
Mantid::Kernel::DblMatrix m_rubw
(2*PiRUBW)^-1
std::vector< double > m_lX
int version() const override
Algorithm's version for identification.
const std::string summary() const override
Algorithm's summary for use in the GUI and help.
std::vector< double > m_kX
std::vector< double > m_hX
cached X values along dimensions h,k,l, dE
size_t m_hIdx
index of h,k,l,dE dimensions in the output workspaces
void exec() override
Execute the algorithm.
void init() override
Initialize the algorithm's properties.
std::vector< double > m_eX
std::vector< Kernel::VMD > calculateIntersections(const double theta, const double phi)
Calculate the points of intersection for the given detector with cuboid surrounding the detector posi...
std::shared_ptr< const MatrixWorkspace > MatrixWorkspace_const_sptr
shared pointer to the matrix workspace base class (const version)
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::Matrix< double > DblMatrix
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,...
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
@ Input
An input workspace.
@ Output
An output workspace.