25#include "MantidIndexing/IndexInfo.h"
27#include "MantidTypes/SpectrumDefinition.h"
29#include <boost/math/special_functions/pow.hpp>
31using boost::math::pow;
76 double minTwoTheta{std::numeric_limits<double>::max()};
77 double maxTwoTheta{std::numeric_limits<double>::lowest()};
78 for (
int width = 0; width < 2; ++width) {
79 const auto offset =
right *
static_cast<double>(width);
80 for (
const auto &pointInRing : capRing) {
81 auto point = pointInRing;
85 const auto current = twoThetaFromLocalPoint(detInfo, detInfoIndex, samplePos, beamDir, point);
86 minTwoTheta = std::min(minTwoTheta, current);
87 maxTwoTheta = std::max(maxTwoTheta, current);
90 const auto beltOffset =
right * 0.5;
91 for (
size_t beltIndex = 0; beltIndex < capRing.size(); beltIndex += 2) {
92 const auto point = capRing[beltIndex] + beltOffset;
93 const auto current = twoThetaFromLocalPoint(detInfo, detInfoIndex, samplePos, beamDir, point);
94 minTwoTheta = std::min(minTwoTheta, current);
95 maxTwoTheta = std::max(maxTwoTheta, current);
97 return std::make_pair(minTwoTheta, maxTwoTheta);
116 if (geometry.
axis.
X() != 0. && geometry.
axis.
Z() != 0) {
117 const auto inverseXZSumSq = 1. / (pow<2>(geometry.
axis.
X()) + pow<2>(geometry.
axis.
Z()));
118 basis1.
setX(std::sqrt(1. - pow<2>(geometry.
axis.
X()) * inverseXZSumSq));
119 basis1.setY(geometry.
axis.
X() * std::sqrt(inverseXZSumSq));
122 const std::array<double, 8> angles{
123 {0., 0.25 * M_PI, 0.5 * M_PI, 0.75 * M_PI, M_PI, 1.25 * M_PI, 1.5 * M_PI, 1.75 * M_PI}};
124 double minTwoTheta{std::numeric_limits<double>::max()};
125 double maxTwoTheta{std::numeric_limits<double>::lowest()};
126 for (
const double &angle : angles) {
127 const auto basePoint =
129 for (
int i = 0; i < 3; ++i) {
130 const auto point = basePoint + geometry.
axis * (0.5 * geometry.
height *
static_cast<double>(i));
131 const auto current = twoThetaFromLocalPoint(detInfo, detInfoIndex, samplePos, beamDir, point);
132 minTwoTheta = std::min(minTwoTheta, current);
133 maxTwoTheta = std::max(maxTwoTheta, current);
136 return std::make_pair(minTwoTheta, maxTwoTheta);
153 const auto maxPoint = bbox.
maxPoint();
154 auto side = maxPoint * sideDir;
155 return twoThetaFromLocalPoint(detInfo, detInfoIndex, samplePos, beamDir, side);
170 double minTwoTheta{std::numeric_limits<double>::max()};
171 double maxTwoTheta{std::numeric_limits<double>::lowest()};
172 const std::array<Mantid::Kernel::V3D, 6> dirs{
173 {
Mantid::Kernel::V3D{1., 0., 0.}, {0., 1., 0.}, {0., 0., 1.}, {-1., 0., 0.}, {0., -1., 0.}, {0., 0., -1.}}};
174 for (
const auto &dir : dirs) {
175 const auto current = twoThetaFromBoundingBox(detInfo, detInfoIndex, samplePos, beamDir, dir);
176 minTwoTheta = std::min(minTwoTheta, current);
177 maxTwoTheta = std::max(maxTwoTheta, current);
179 return std::make_pair(minTwoTheta, maxTwoTheta);
195 switch (shape->shape()) {
205 return generalTwoThetaRange(detInfo, detInfoIndex, samplePos, beamDir);
216std::pair<double, double> twoThetasFromTable(
const Mantid::detid_t detID,
const std::vector<int> &detectorIDs,
217 const std::vector<double> &lowers,
const std::vector<double> &uppers) {
218 const auto range = std::equal_range(detectorIDs.cbegin(), detectorIDs.cend(),
static_cast<int>(detID));
219 if (std::distance(range.first, range.second) > 1) {
220 throw std::invalid_argument(
"Duplicate detector IDs in 'DetectorTwoThetaRanges'.");
222 if (range.first == detectorIDs.cend()) {
223 throw std::invalid_argument(
"No min/max 2thetas found for detector ID " +
std::to_string(detID));
225 const auto index = std::distance(detectorIDs.cbegin(), range.first);
226 const auto minmax = std::make_pair(lowers[
index], uppers[
index]);
227 if (minmax.first <= 0) {
228 throw std::invalid_argument(
"Non-positive min 2theta for detector ID " +
std::to_string(detID));
230 if (minmax.first > M_PI) {
231 throw std::invalid_argument(
"Min 2theta greater than pi for detector ID " +
std::to_string(detID));
233 if (minmax.second > M_PI) {
234 throw std::invalid_argument(
"Max 2theta greater than pi for detector ID" +
std::to_string(detID));
236 if (minmax.first >= minmax.second) {
237 throw std::invalid_argument(
"Min 2theta larger than max for detector ID " +
std::to_string(detID));
281 std::map<std::string, std::string> result;
284 result[
"InputWorkspace"] =
"InputWorkspace is of Incorrect type. Please "
285 "provide a MatrixWorkspace as the "
291 if (tableWS->columnCount() != 3) {
292 result[
"DetectorTwoThetaRanges"] =
"DetectorTwoThetaRanges requires 3 columns";
295 else if (!tableWS->getColumn(0)->isType<
int>()) {
296 result[
"DetectorTwoThetaRanges"] =
"The first column of DetectorTwoThetaRanges should be of type int";
299 else if (!tableWS->getColumn(1)->isType<
double>()) {
300 result[
"DetectorTwoThetaRanges"] =
"The second column of DetectorTwoThetaRanges should be of type "
304 else if (!tableWS->getColumn(2)->isType<
double>()) {
305 result[
"DetectorTwoThetaRanges"] =
"The third column of DetectorTwoThetaRanges should be of type double";
308 else if (tableWS->rowCount() != inputWS->getNumberHistograms()) {
309 result[
"DetectorTwoThetaRanges"] =
"The table and workspace do not have the same number of detectors";
325 g_log.
debug() <<
"Workspace type: " << outputWS->id() <<
'\n';
327 const size_t nEnergyBins = inputWS->blocksize();
328 const size_t nHistos = inputWS->getNumberHistograms();
331 std::vector<SpectrumDefinition> detIDMapping(outputWS->getNumberHistograms());
334 const size_t nreports(nHistos * nEnergyBins);
335 m_progress = std::make_unique<API::Progress>(
this, 0.0, 1.0, nreports);
342 std::vector<double> par = inputWS->getInstrument()->getNumberParameter(
"detector-neighbour-offset");
346 g_log.
debug() <<
"Offset: " << par[0] <<
'\n';
351 const auto &
X = inputWS->x(0);
353 const auto &inputIndices = inputWS->indexInfo();
354 const auto &spectrumInfo = inputWS->spectrumInfo();
357 for (int64_t i = 0; i < static_cast<int64_t>(nHistos); ++i) {
360 if (spectrumInfo.isMasked(i) || spectrumInfo.isMonitor(i)) {
368 const auto specNo =
static_cast<specnum_t>(inputIndices.spectrumNumber(i));
369 std::stringstream logStream;
370 for (
size_t j = 0; j < nEnergyBins; ++j) {
371 m_progress->report(
"Computing polygon intersections");
374 const double dE_j =
X[j];
375 const double dE_jp1 =
X[j + 1];
380 const V2D lr(dE_jp1, lrQ);
383 if (
g_log.
is(Logger::Priority::PRIO_DEBUG)) {
384 logStream <<
"Spectrum=" << specNo <<
", lower theta=" << thetaLower *
rad2deg
385 <<
", upper theta=" << thetaUpper *
rad2deg <<
". QE polygon: ll=" << ll <<
", lr=" << lr
386 <<
", ur=" << ur <<
", ul=" << ul <<
"\n";
393 const MantidVec::difference_type qIndex = std::upper_bound(
m_Qout.begin(),
m_Qout.end(), lrQ) -
m_Qout.begin();
394 if (qIndex != 0 && qIndex <
static_cast<int>(
m_Qout.size())) {
400 detIDMapping[qIndex - 1].add(spectrumInfo.spectrumDefinition(i)[0].first);
404 if (
g_log.
is(Logger::Priority::PRIO_DEBUG)) {
413 outputWS->finalize();
417 auto outputIndices = outputWS->indexInfo();
418 outputIndices.setSpectrumDefinitions(std::move(detIDMapping));
419 outputWS->setIndexInfo(outputIndices);
424 replaceNans->setChild(
true);
425 replaceNans->initialize();
426 replaceNans->setProperty(
"InputWorkspace", outputWS);
427 replaceNans->setProperty(
"OutputWorkspace", outputWS);
428 replaceNans->setProperty(
"NaNValue", 0.0);
429 replaceNans->setProperty(
"InfinityValue", 0.0);
430 replaceNans->setProperty(
"BigNumberThreshold", DBL_MAX);
431 replaceNans->execute();
444 const size_t nhist =
workspace.getNumberHistograms();
445 constexpr double skipDetector{-1.};
450 const auto referenceFrame = inst->getReferenceFrame();
451 const auto beamDir = referenceFrame->vecPointingAlongBeam();
453 const auto &detectorInfo =
workspace.detectorInfo();
454 const auto &spectrumInfo =
workspace.spectrumInfo();
455 const auto samplePos = spectrumInfo.samplePosition();
456 for (
size_t i = 0; i < nhist; ++i) {
457 m_progress->report(
"Calculating detector angles");
460 if (!spectrumInfo.hasDetectors(i) || spectrumInfo.isMonitor(i)) {
464 const auto &det = spectrumInfo.detector(i);
468 }
catch (std::runtime_error &) {
472 double minTwoTheta{spectrumInfo.twoTheta(i)};
473 double maxTwoTheta{minTwoTheta};
474 if (spectrumInfo.hasUniqueDetector(i)) {
475 const auto detInfoIndex = detectorInfo.indexOf(det.getID());
476 std::tie(minTwoTheta, maxTwoTheta) = minMaxTwoTheta(detectorInfo, detInfoIndex, samplePos, beamDir);
478 const auto &group =
dynamic_cast<const DetectorGroup &
>(det);
480 for (
const auto id : ids) {
481 const auto detInfoIndex = detectorInfo.indexOf(
id);
482 const auto current = minMaxTwoTheta(detectorInfo, detInfoIndex, samplePos, beamDir);
483 minTwoTheta = std::min(minTwoTheta, current.first);
484 maxTwoTheta = std::max(maxTwoTheta, current.second);
489 if (
g_log.
is(Logger::Priority::PRIO_DEBUG)) {
504 const size_t nHistos =
workspace.getNumberHistograms();
506 bool ignoreMasked =
true;
507 const int numNeighbours = 4;
513 const auto &spectrumInfo =
workspace.spectrumInfo();
515 for (
size_t i = 0; i < nHistos; ++i) {
516 m_progress->report(
"Calculating detector angular widths");
519 if (!spectrumInfo.hasDetectors(i) || spectrumInfo.isMonitor(i)) {
527 double thetaWidth = std::numeric_limits<double>::lowest();
530 const double theta = spectrumInfo.twoTheta(i);
533 const specnum_t deltaMinus1 = inSpec - 1;
537 for (
auto &neighbour : neighbours) {
539 if (spec == deltaPlus1 || spec == deltaMinus1 || spec == deltaPlusT || spec == deltaMinusT) {
540 const double theta_n = spectrumInfo.twoTheta(spec - 1) * 0.5;
542 const double dTheta = std::abs(theta - theta_n);
543 thetaWidth = std::max(thetaWidth, dTheta);
548 if (
g_log.
is(Logger::Priority::PRIO_DEBUG)) {
549 g_log.
debug() <<
"Detector at spectrum = " << inSpec <<
", width = " << thetaWidth *
rad2deg <<
" degrees\n";
556 constexpr double skipDetector{-1.};
557 const size_t nhist =
workspace.getNumberHistograms();
560 const auto &detIDs = angleTable.
getColVector<
int>(
"Detector ID");
561 const auto &lowers = angleTable.
getColVector<
double>(
"Min two theta");
562 const auto &uppers = angleTable.
getColVector<
double>(
"Max two theta");
563 const auto &spectrumInfo =
workspace.spectrumInfo();
564 for (
size_t i = 0; i < nhist; ++i) {
565 m_progress->report(
"Reading detector angles");
566 if (!spectrumInfo.hasDetectors(i) || spectrumInfo.isMonitor(i)) {
569 const auto &det = spectrumInfo.detector(i);
572 }
catch (std::runtime_error &) {
576 if (spectrumInfo.hasUniqueDetector(i)) {
577 const auto twoThetas = twoThetasFromTable(det.getID(), detIDs, lowers, uppers);
581 const auto &group =
dynamic_cast<const DetectorGroup &
>(det);
583 double minTwoTheta{spectrumInfo.twoTheta(i)};
584 double maxTwoTheta{minTwoTheta};
585 for (
const auto id : ids) {
586 const auto twoThetas = twoThetasFromTable(
id, detIDs, lowers, uppers);
587 minTwoTheta = std::min(minTwoTheta, twoThetas.first);
588 maxTwoTheta = std::max(maxTwoTheta, twoThetas.second);
593 if (
g_log.
is(Logger::Priority::PRIO_DEBUG)) {
#define DECLARE_ALGORITHM(classname)
IPeaksWorkspace_sptr workspace
std::map< DeltaEMode::Type, std::string > index
#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.
Mantid::Kernel::Quat(ComponentInfo::* rotation)(const size_t) const
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.
Base MatrixWorkspace Abstract Class.
WorkspaceNearestNeighbourInfo provides easy access to nearest-neighbour information for a workspace.
std::map< specnum_t, Kernel::V3D > getNeighboursExact(specnum_t spec) const
Queries the WorkspaceNearestNeighbours object for the selected spectrum number.
std::unique_ptr< API::Progress > m_progress
Progress reporter.
Converts a 2D workspace that has axes of energy transfer against spectrum number to one that gives in...
void exec() override
Run the algorithm.
std::vector< double > m_twoThetaLowers
Array for the lower 2theta angles, in radians.
int m_detNeighbourOffset
Offset for finding neighbor in nearest tube.
std::map< std::string, std::string > validateInputs() override
validate the inputs
void initAngularCachesNonPSD(const API::MatrixWorkspace &workspace)
Init the theta index.
std::vector< double > m_twoThetaUppers
Array for the upper 2theta angles, in radians.
void initAngularCachesPSD(const API::MatrixWorkspace &workspace)
Get angles and calculate angular widths.
int version() const override
Algorithm's version for identification.
SofQCommon m_EmodeProperties
void init() override
Initialize the algorithm.
const std::string category() const override
Algorithm's category for identification.
const std::string name() const override
Algorithm's name for identification.
void initAngularCachesTable(const API::MatrixWorkspace &workspace, const DataObjects::TableWorkspace &widthTable)
std::vector< double > m_Qout
Output Q axis.
static void createCommonInputProperties(API::Algorithm &alg)
Create the input properties on the given algorithm object.
TableWorkspace is an implementation of Workspace in which the data are organised in columns of same s...
std::vector< T > & getColVector(size_t index)
get access to column vector for index i.
A simple structure that defines an axis-aligned cuboid shaped bounding box for a geometrical object.
const Kernel::V3D & maxPoint() const
Returns the min point of the box.
Constructive Solid Geometry object.
const detail::ShapeInfo & shapeInfo() const override
Holds a collection of detectors.
std::vector< detid_t > getDetectorIDs() const
What detectors are contained in the group?
Geometry::DetectorInfo is an intermediate step towards a DetectorInfo that is part of Instrument-2....
Kernel::Quat rotation(const size_t index) const
Returns the rotation of the detector with given index.
Kernel::V3D position(const size_t index) const
Returns the position of the detector with given index.
const Geometry::IDetector & detector(const size_t index) const
Return a const reference to the detector with given index.
virtual const std::shared_ptr< const IObject > shape() const =0
Returns the shape of the Object.
A ConvexPolygon with only 4 vertices.
CylinderGeometry cylinderGeometry() const
CuboidGeometry cuboidGeometry() const
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.
bool is(int level) const
Returns true if at least the given log level is set.
void rotate(V3D &) const
Rotate a vector.
Implements a 2-dimensional vector embedded in a 3D space, i.e.
constexpr double X() const noexcept
Get x.
constexpr V3D cross_prod(const V3D &v) const noexcept
Cross product (this * argument)
double angle(const V3D &) const
Angle between this and another vector.
void setX(const double xx) noexcept
Set is x position.
constexpr double Z() const noexcept
Get z.
std::shared_ptr< const MatrixWorkspace > MatrixWorkspace_const_sptr
shared pointer to the matrix workspace base class (const version)
Geometry::IDetector_const_sptr DetConstPtr
std::map< specnum_t, Mantid::Kernel::V3D > SpectraDistanceMap
MANTID_DATAOBJECTS_DLL void normaliseOutput(const API::MatrixWorkspace_sptr &outputWS, const API::MatrixWorkspace_const_sptr &inputWS, API::Progress *progress=nullptr)
Compute sqrt of errors and put back in bin width division if necessary.
MANTID_DATAOBJECTS_DLL void rebinToFractionalOutput(const Geometry::Quadrilateral &inputQ, const API::MatrixWorkspace_const_sptr &inputWS, const size_t i, const size_t j, DataObjects::RebinnedOutput &outputWS, const std::vector< double > &verticalAxis, const DataObjects::RebinnedOutput_const_sptr &inputRB=nullptr)
Rebin the input quadrilateral to to output grid.
MANTID_DATAOBJECTS_DLL void finalizeFractionalRebin(DataObjects::RebinnedOutput &outputWS)
Set finalize flag after fractional rebinning loop.
std::shared_ptr< TableWorkspace > TableWorkspace_sptr
shared pointer to Mantid::DataObjects::TableWorkspace
std::shared_ptr< RebinnedOutput > RebinnedOutput_sptr
shared pointer to the RebinnedOutput class
constexpr double rad2deg
Radians to degrees conversion factor.
std::shared_ptr< const Mantid::Geometry::IDetector > IDetector_const_sptr
Shared pointer to IDetector (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.
int32_t detid_t
Typedef for a detector ID.
int32_t specnum_t
Typedef for a spectrum Number.
std::string to_string(const wide_integer< Bits, Signed > &n)
double getEFixed(const Geometry::IDetector &det) const
Get the efixed value for the given detector.
double q(const double deltaE, const double twoTheta, const Geometry::IDetector *det) const
Calculate the Q value.
void initCachedValues(const API::MatrixWorkspace &workspace, API::Algorithm *const hostAlgorithm)
The procedure analyses emode and efixed properties provided to the algorithm and identify the energy ...
const Kernel::V3D & leftFrontTop
const Kernel::V3D & rightFrontBottom
const Kernel::V3D & leftBackBottom
const Kernel::V3D & leftFrontBottom
const Kernel::V3D & centreOfBottomBase