25#include "MantidIndexing/IndexInfo.h"
32#include "MantidParallel/Communicator.h"
33#include "MantidTypes/SpectrumDefinition.h"
40using namespace Kernel;
42using namespace Geometry;
43using namespace DataObjects;
48 auto dataVal = std::make_shared<CompositeValidator>();
54 "Particle counts as a function of wavelength");
56 "Name of the workspace that will contain the result of the calculation");
58 "A comma separated list of first bin boundary, width, last bin boundary. "
60 "this can be followed by a comma and more widths and last boundary "
62 "Negative width values indicate logarithmic binning.");
64 "Scaling to apply to each spectrum. Must have\n"
65 "the same number of spectra as the DetBankWorkspace");
66 auto wavVal = std::make_shared<CompositeValidator>();
71 "Scaling to apply to each bin.\n"
72 "Must have the same number of bins as the DetBankWorkspace");
75 "Scaling that depends on both pixel and wavelength together.\n"
76 "Must have the same number of bins and spectra as the DetBankWorkspace.");
77 declareProperty(
"AccountForGravity",
false,
"Whether to correct for the effects of gravity");
78 declareProperty(
"SolidAngleWeighting",
true,
"If true, pixels will be weighted by their solid angle.",
80 auto mustBePositive = std::make_shared<BoundedValidator<double>>();
81 mustBePositive->setLower(0.0);
84 "To increase resolution some wavelengths are excluded within "
85 "this distance from the beam center (mm). Note that RadiusCut\n"
86 " and WaveCut both need to be larger than 0 to affect \n"
87 "the effective cutoff. See the algorithm description for\n"
88 " a detailed explanation of the cutoff.");
91 "To increase resolution by starting to remove some wavelengths below this"
92 " threshold (angstrom). Note that WaveCut\n"
93 " and RadiusCut both need to be larger than 0 to affect \n"
94 "on the effective cutoff. See the algorithm description for\n"
95 " a detailed explanation of the cutoff.");
98 "Set to true to output two additional workspaces which will "
99 "have the names OutputWorkspace_sumOfCounts "
100 "OutputWorkspace_sumOfNormFactors. The division of "
101 "_sumOfCounts and _sumOfNormFactors equals the workspace"
102 " returned by the property OutputWorkspace "
103 "(default is false).");
104 declareProperty(
"ExtraLength", 0.0, mustBePositive,
"Additional length for gravity correction.");
108 "Workspace to calculate the Q resolution.\n");
120 const bool doGravity =
getProperty(
"AccountForGravity");
127 g_log.
debug() <<
"All input workspaces were found to be valid\n";
129 double const *
const binNorms = waveAdj ? &(waveAdj->y(0)[0]) :
nullptr;
131 double const *
const binNormEs = waveAdj ? &(waveAdj->e(0)[0]) :
nullptr;
137 auto useQResolution =
static_cast<bool>(qResolution);
142 auto &QOut = outputWS->x(0);
143 auto &YOut = outputWS->mutableY(0);
144 auto &EOutTo2 = outputWS->mutableE(0);
146 HistogramData::HistogramY normSum(YOut.size(), 0.0);
148 HistogramData::HistogramE normError2(EOutTo2.size(), 0.0);
150 HistogramData::HistogramDx qResolutionOut(YOut.size(), 0.0);
152 const auto numSpec =
static_cast<int>(
m_dataWS->getNumberHistograms());
155 const auto &spectrumInfo =
m_dataWS->spectrumInfo();
157 for (
int i = 0; i < numSpec; ++i) {
159 if (!spectrumInfo.hasDetectors(i)) {
160 g_log.
warning() <<
"Workspace index " << i <<
" (SpectrumIndex = " <<
m_dataWS->getSpectrum(i).getSpectrumNo()
161 <<
") has no detector assigned to it - discarding\n";
165 if (spectrumInfo.isMonitor(i) || spectrumInfo.isMasked(i))
171 const size_t wavStart =
173 if (wavStart >=
m_dataWS->y(i).size()) {
178 const size_t numWavbins =
m_dataWS->y(i).size() - wavStart;
182 HistogramData::HistogramY _noDirectUseStorage_(3 * numWavbins);
184 auto norms = _noDirectUseStorage_.begin();
187 auto normETo2s = norms + numWavbins;
189 auto QIn = normETo2s + numWavbins;
198 auto YIn =
m_dataWS->y(i).cbegin() + wavStart;
199 auto EIn =
m_dataWS->e(i).cbegin() + wavStart;
207 auto QResIn = useQResolution ? (qResolution->y(i).cbegin() + wavStart) : YIn;
211 auto loc = QOut.cend();
214 const auto end =
m_dataWS->y(i).cend();
215 for (; YIn != end; ++YIn, ++EIn, ++QIn, ++norms, ++normETo2s) {
220 if ((loc != QOut.begin()) && (loc != QOut.end())) {
222 const size_t bin = loc - QOut.begin() - 1;
225 normSum[bin] += *norms;
228 EOutTo2[bin] += (*EIn) * (*EIn);
229 normError2[bin] += *normETo2s;
230 if (useQResolution) {
231 auto QBin = (QOut[bin + 1] - QOut[bin]);
235 qResolutionOut[bin] += (*YIn) * std::sqrt((*QResIn) * (*QResIn) + QBin * QBin / 12.0);
241 if (useQResolution) {
249 const auto &inSpec =
m_dataWS->getSpectrum(i);
250 auto &outSpec = outputWS->getSpectrum(0);
251 outSpec.addDetectorIDs(inSpec.getDetectorIDs());
260 auto size =
static_cast<int>(YOut.size());
262 for (
int rank = 1; rank <
communicator().size(); ++rank) {
263 HistogramData::HistogramY
y(YOut.size());
264 HistogramData::HistogramE e2(YOut.size());
275 std::vector<detid_t> detIds(detCount);
276 communicator().recv(rank, tag, detIds.data(), detCount);
277 outputWS->getSpectrum(0).addDetectorIDs(detIds);
280 communicator().send(0, tag, YOut.rawData().data(), size);
281 communicator().send(0, tag, EOutTo2.rawData().data(), size);
282 communicator().send(0, tag, normSum.rawData().data(), size);
283 communicator().send(0, tag, normError2.rawData().data(), size);
284 const auto detIdSet = outputWS->getSpectrum(0).getDetectorIDs();
285 std::vector<detid_t> detIds(detIdSet.begin(), detIdSet.end());
286 const auto nDets =
static_cast<int>(detIds.size());
292 if (useQResolution) {
296 auto countsIterator = YOut.cbegin();
297 auto qResolutionIterator = qResolutionOut.begin();
298 for (; countsIterator != YOut.end(); ++countsIterator, ++qResolutionIterator) {
301 if ((*countsIterator) > 0.0) {
302 *qResolutionIterator = (*qResolutionIterator) / (*countsIterator);
305 outputWS->setPointStandardDeviations(0, std::move(qResolutionOut));
311 ws_sumOfCounts->setSharedX(0, outputWS->sharedX(0));
313 ws_sumOfCounts->mutableY(0) = YOut;
314 ws_sumOfCounts->setSharedDx(0, outputWS->sharedDx(0));
315 ws_sumOfCounts->setFrequencyVariances(0, outputWS->e(0));
318 ws_sumOfNormFactors->setSharedX(0, outputWS->sharedX(0));
319 ws_sumOfNormFactors->mutableY(0) = normSum;
320 ws_sumOfNormFactors->setSharedDx(0, outputWS->sharedDx(0));
321 ws_sumOfNormFactors->setFrequencyVariances(0, normError2);
323 helper.
outputParts(
this, ws_sumOfCounts, ws_sumOfNormFactors);
324 }
else if (doOutputParts) {
328 progress.report(
"Normalizing I(Q)");
330 normalize(normSum, normError2, YOut, EOutTo2);
344 HistogramData::BinEdges XOut(0);
348 Indexing::IndexInfo indexInfo(
349 1,
communicator().rank() == 0 ? Parallel::StorageMode::MasterOnly : Parallel::StorageMode::Cloned,
351 indexInfo.setSpectrumDefinitions(std::vector<SpectrumDefinition>(1));
352 auto outputWS = create<MatrixWorkspace>(*
m_dataWS, indexInfo, XOut);
354 outputWS->setYUnitLabel(
"1/cm");
356 outputWS->setDistribution(
true);
358 outputWS->getSpectrum(0).clearDetectorIDs();
359 outputWS->getSpectrum(0).setSpectrumNo(1);
384 double const *
const binNormEs, HistogramData::HistogramY::iterator norm,
385 HistogramData::HistogramY::iterator normETo2)
const {
386 double detectorAdj, detAdjErr;
387 pixelWeight(pixelAdj, wsIndex, detectorAdj, detAdjErr);
389 for (
auto n = norm, e = normETo2;
n != normETo2; ++
n, ++e) {
391 *e = detAdjErr * detAdjErr;
394 if (binNorms && binNormEs) {
397 addWaveAdj(binNorms + wavStart, binNormEs + wavStart, norm, normETo2, wavePixelAdj->y(wsIndex).begin() + wavStart,
398 wavePixelAdj->e(wsIndex).begin() + wavStart);
400 addWaveAdj(binNorms + wavStart, binNormEs + wavStart, norm, normETo2);
402 normToMask(wavStart, wsIndex, norm, normETo2);
417 double &
error)
const {
418 const auto &detectorInfo =
m_dataWS->detectorInfo();
419 const V3D samplePos = detectorInfo.samplePosition();
423 for (
const auto detID :
m_dataWS->getSpectrum(wsIndex).getDetectorIDs()) {
424 const auto index = detectorInfo.indexOf(detID);
425 if (!detectorInfo.isMasked(
index))
426 weight += detectorInfo.detector(
index).solidAngle(samplePos);
431 if (weight < 1e-200) {
432 throw std::logic_error(
"Invalid (zero or negative) solid angle for one detector");
436 weight *= pixelAdj->readY(wsIndex)[0];
437 error = weight * pixelAdj->readE(wsIndex)[0];
455void Q1D2::addWaveAdj(
const double *c,
const double *Dc, HistogramData::HistogramY::iterator bInOut,
456 HistogramData::HistogramY::iterator e2InOut)
const {
468 const auto end = e2InOut;
469 for (; bInOut != end; ++e2InOut, ++c, ++Dc, ++bInOut) {
471 *e2InOut = ((*e2InOut) * (*c) * (*c)) + ((*Dc) * (*Dc) * (*bInOut) * (*bInOut));
473 *bInOut = (*bInOut) * (*c);
492void Q1D2::addWaveAdj(
const double *c,
const double *Dc, HistogramData::HistogramY::iterator bInOut,
493 HistogramData::HistogramY::iterator e2InOut,
494 HistogramData::HistogramY::const_iterator wavePixelAdjData,
495 HistogramData::HistogramE::const_iterator wavePixelAdjError)
const {
517 const auto end = e2InOut;
518 for (; bInOut != end; ++e2InOut, ++c, ++Dc, ++bInOut, ++wavePixelAdjData, ++wavePixelAdjError) {
520 *e2InOut = ((*e2InOut) * (*c) * (*c) * (*wavePixelAdjData) * (*wavePixelAdjData)) +
521 ((*Dc) * (*Dc) * (*bInOut) * (*bInOut) * (*wavePixelAdjData) * (*wavePixelAdjData)) +
522 ((*wavePixelAdjError) * (*wavePixelAdjError) * (*c) * (*c) * (*bInOut) * (*bInOut));
524 *bInOut = (*bInOut) * (*c) * (*wavePixelAdjData);
537void Q1D2::normToMask(
const size_t offSet,
const size_t wsIndex,
const HistogramData::HistogramY::iterator theNorms,
538 const HistogramData::HistogramY::iterator errorSquared)
const {
540 if (
m_dataWS->hasMaskedBins(wsIndex)) {
544 MatrixWorkspace::MaskList::const_iterator it;
545 for (it = mask.begin(); it != mask.end(); ++it) {
546 size_t outBin = it->first;
547 if (outBin < offSet) {
554 const double factor = 1.0 - (it->second);
555 *(theNorms + outBin) *= factor;
556 *(errorSquared + outBin) *= factor * factor;
575 const size_t offset, HistogramData::HistogramY::iterator Qs,
const double extraLength)
const {
576 static const double FOUR_PI = 4.0 * M_PI;
579 auto waves =
m_dataWS->x(wsInd).cbegin() + offset;
581 const auto end =
m_dataWS->x(wsInd).end() - 1;
584 for (; waves != end; ++Qs, ++waves) {
587 const double lambda = 0.5 * (*(waves + 1) + (*waves));
592 *Qs = FOUR_PI * sinTheta /
lambda;
597 const double factor = 2.0 * FOUR_PI * sin(spectrumInfo.
twoTheta(wsInd) * 0.5);
598 for (; waves != end; ++Qs, ++waves) {
601 *Qs = factor / (*(waves + 1) + (*waves));
619 HistogramData::HistogramX::const_iterator &loc)
const {
620 if (loc != OutQs.end()) {
621 while (loc != OutQs.begin()) {
622 if ((QToFind >= *(loc - 1)) && (QToFind < *loc)) {
627 if (QToFind < *loc) {
633 if (OutQs.empty() || QToFind > *(loc - 1)) {
641 loc = std::lower_bound(OutQs.begin(), OutQs.end(), QToFind);
654void Q1D2::normalize(
const HistogramData::HistogramY &normSum,
const HistogramData::HistogramE &normError2,
655 HistogramData::HistogramY &counts, HistogramData::HistogramE &errors)
const {
656 for (
size_t k = 0; k < counts.size(); ++k) {
658 const double c = normSum[k];
659 const double a = counts[k] /= c;
666 const double aOverc = a / c;
667 errors[k] = std::sqrt(errors[k] / (c * c) + normError2[k] * aOverc * aOverc);
672void checkStorageMode(
const std::map<std::string, Parallel::StorageMode> &storageModes,
const std::string &name) {
673 if (storageModes.count(name) && storageModes.at(name) != Parallel::StorageMode::Cloned)
674 throw std::runtime_error(name +
" must have " + Parallel::toString(Parallel::StorageMode::Cloned));
678Parallel::ExecutionMode
680 if (storageModes.count(
"PixelAdj") || storageModes.count(
"WavePixelAdj") || storageModes.count(
"QResolution"))
681 throw std::runtime_error(
"Using in PixelAdj, WavePixelAdj, or QResolution in an MPI run of " +
name() +
682 " is currently not supported.");
683 checkStorageMode(storageModes,
"WavelengthAdj");
684 return Parallel::getCorrespondingExecutionMode(storageModes.at(
"DetBankWorkspace"));
#define DECLARE_ALGORITHM(classname)
const std::vector< double > * lambda
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.
Base class from which all concrete algorithm classes should be derived.
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.
void progress(double p, const std::string &msg="", double estimatedTime=0.0, int progressPrecision=0)
Sends ProgressNotification.
const Parallel::Communicator & communicator() const
Returns a const reference to the (MPI) communicator of the algorithm.
A validator which provides a TENTATIVE check that a workspace contains common bins in each spectrum.
A validator which checks that a workspace contains histogram data (the default) or point data as requ...
A validator which checks that a workspace has a valid instrument.
std::map< size_t, double > MaskList
Masked bins for each spectrum are stored as a set of pairs containing <bin index, weight>
Helper class for reporting progress from algorithms.
API::SpectrumInfo is an intermediate step towards a SpectrumInfo that is part of Instrument-2....
double twoTheta(const size_t index) const
Returns the scattering angle 2 theta in radians (angle w.r.t.
A property class for workspaces.
A validator which checks that the unit of the workspace referred to by a WorkspaceProperty is the exp...
A helper class for calculating neutron's gravitional drop.
double calcSinTheta(const double wavAngstroms) const
Caclulates the sin of the that the neutron left the sample at, before the effect of gravity.
Parallel::ExecutionMode getParallelExecutionMode(const std::map< std::string, Parallel::StorageMode > &storageModes) const override
Get correct execution mode based on input storage modes for an MPI run.
void normToMask(const size_t offSet, const size_t wsIndex, const HistogramData::HistogramY::iterator theNorms, const HistogramData::HistogramY::iterator errorSquared) const
Scaled to bin masking, to the normalization.
void calculateNormalization(const size_t wavStart, const size_t wsIndex, const API::MatrixWorkspace_const_sptr &pixelAdj, const API::MatrixWorkspace_const_sptr &wavePixelAdj, double const *const binNorms, double const *const binNormEs, HistogramData::HistogramY::iterator norm, HistogramData::HistogramY::iterator normETo2) const
Calculate the normalization term for each output bin.
void normalize(const HistogramData::HistogramY &normSum, const HistogramData::HistogramE &normError2, HistogramData::HistogramY &counts, HistogramData::HistogramE &errors) const
Divides the number of counts in each output Q bin by the wrighting ("number that would expected to ar...
void pixelWeight(const API::MatrixWorkspace_const_sptr &pixelAdj, const size_t wsIndex, double &weight, double &error) const
Calculates the normalisation for the spectrum specified by the index number that was passed as the so...
Q1D2()
Default constructor.
const std::string name() const override
Algorithm's name.
void convertWavetoQ(const API::SpectrumInfo &spectrumInfo, const size_t wsInd, const bool doGravity, const size_t offset, HistogramData::HistogramY::iterator Qs, const double extraLength) const
Fills a vector with the Q values calculated from the wavelength bin centers from the input workspace ...
API::MatrixWorkspace_sptr setUpOutputWorkspace(const std::vector< double > &binParams) const
Creates the output workspace, its size, units, etc.
API::MatrixWorkspace_const_sptr m_dataWS
the experimental workspace with counts across the detector
void init() override
Initialisation code.
void getQBinPlus1(const HistogramData::HistogramX &OutQs, const double QToFind, HistogramData::HistogramY::const_iterator &loc) const
This is a slightly "clever" method as it makes some guesses about where is best to look for the right...
void exec() override
Execution code.
void addWaveAdj(const double *c, const double *Dc, HistogramData::HistogramY::iterator bInOut, HistogramData::HistogramY::iterator e2InOut) const
Calculates the contribution to the normalization terms from each bin in a spectrum.
Helper class for the Q1D and Qxy algorithms.
size_t waveLengthCutOff(const API::MatrixWorkspace_const_sptr &dataWS, const API::SpectrumInfo &spectrumInfo, const double RCut, const double WCut, const size_t wsInd) const
Finds the first index number of the first wavelength bin that should included based on the the calcul...
void outputParts(API::Algorithm *alg, const API::MatrixWorkspace_sptr &sumOfCounts, const API::MatrixWorkspace_sptr &sumOfNormFactors)
This method performs the common work between Qxy and Q1D2 if algorihtm parameter OutputParts=True.
void examineInput(const API::MatrixWorkspace_const_sptr &dataWS, const API::MatrixWorkspace_const_sptr &binAdj, const API::MatrixWorkspace_const_sptr &detectAdj, const API::MatrixWorkspace_const_sptr &qResolution)
Checks if workspaces input to Q1D or Qxy are reasonable.
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 warning(const std::string &msg)
Logs at warning level.
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
std::shared_ptr< const MatrixWorkspace > MatrixWorkspace_const_sptr
shared pointer to the matrix workspace base class (const version)
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
int MANTID_KERNEL_DLL createAxisFromRebinParams(const std::vector< double > ¶ms, std::vector< double > &xnew, const bool resize_xnew=true, const bool full_bins_only=false, const double xMinHint=std::nan(""), const double xMaxHint=std::nan(""), const bool useReverseLogarithmic=false, const double power=-1)
Creates a new output X array given a 'standard' set of rebinning parameters.
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.
@ Input
An input workspace.
@ Output
An output workspace.