19#include "MantidHistogramData/LinearGenerator.h"
33using namespace Kernel;
35using namespace Geometry;
38 auto wsValidator = std::make_shared<CompositeValidator>();
44 "The corrected data in units of wavelength.");
46 "The name to use for the corrected workspace.");
48 auto mustBePositive = std::make_shared<BoundedValidator<double>>();
49 mustBePositive->setLower(1.0e-12);
51 declareProperty(
"MaxQxy", -1.0, mustBePositive,
"The upper limit of the Qx-Qy grid (goes from -MaxQxy to +MaxQxy).");
52 declareProperty(
"DeltaQ", -1.0, mustBePositive,
"The dimension of a Qx-Qy cell.");
53 declareProperty(
"IQxQyLogBinning",
false,
"I(qx,qy) log binning when binning is not specified.",
56 "The scaling to apply to each spectrum e.g. for detector "
57 "efficiency, must have just one bin per spectrum and the "
58 "same number of spectra as DetBankWorkspace.");
59 auto wavVal = std::make_shared<CompositeValidator>();
64 "The scaling to apply to each bin to account for monitor "
65 "counts, transmission fraction, etc. Must be one spectrum "
66 "with the same binning as the InputWorkspace, the same units "
67 "(counts) and the same [[ConvertToDistribution|distribution "
70 declareProperty(
"SolidAngleWeighting",
true,
"If true, pixels will be weighted by their solid angle.",
72 auto mustBePositive2 = std::make_shared<BoundedValidator<double>>();
73 mustBePositive2->setLower(0.0);
75 "The minimum distance in mm from the beam center at which "
76 "all wavelengths are used in the calculation (see section "
77 "[[Q1D#Resolution and Cutoffs|Resolution and Cutoffs]])");
79 "The shortest wavelength in angstrom at which counts should "
80 "be summed from all detector pixels (see section "
81 "[[Q1D#Resolution and Cutoffs|Resolution and Cutoffs]])");
83 "Set to true to output two additional workspaces which will "
84 "have the names OutputWorkspace_sumOfCounts "
85 "OutputWorkspace_sumOfNormFactors. The division of "
86 "_sumOfCounts and _sumOfNormFactors equals the workspace "
87 "returned by the property OutputWorkspace");
88 declareProperty(
"ExtraLength", 0.0, mustBePositive2,
"Additional length for gravity correction.");
95 const bool doGravity =
getProperty(
"AccountForGravity");
96 const bool doSolidAngle =
getProperty(
"SolidAngleWeighting");
101 g_log.
debug() <<
"All input workspaces were found to be valid\n";
110 for (
size_t i = 0; i < weights->getNumberHistograms(); ++i)
111 weights->setSharedX(i, outputWorkspace->sharedX(0));
113 const size_t numSpec = inputWorkspace->getNumberHistograms();
114 const size_t numBins = inputWorkspace->blocksize();
117 Progress prog(
this, 0.05, 1.0, numSpec);
119 const auto &spectrumInfo = inputWorkspace->spectrumInfo();
120 const auto &detectorInfo = inputWorkspace->detectorInfo();
124 const V3D samplePos = spectrumInfo.samplePosition();
126 for (int64_t i = 0; i < int64_t(numSpec); ++i) {
127 if (!spectrumInfo.hasDetectors(i)) {
128 g_log.
warning() <<
"Workspace index " << i <<
" has no detector assigned to it - discarding\n";
133 if (spectrumInfo.isMonitor(i) || spectrumInfo.isMasked(i))
138 const size_t wavStart =
140 if (wavStart >= inputWorkspace->y(i).size()) {
145 V3D detPos = spectrumInfo.position(i) - samplePos;
149 double phi = atan2(detPos.
Y(), detPos.
X());
152 double sinTheta = sin(spectrumInfo.twoTheta(i) * 0.5);
155 const auto &
X = inputWorkspace->x(i);
156 const auto &
Y = inputWorkspace->y(i);
157 const auto &E = inputWorkspace->e(i);
159 const auto &axis = outputWorkspace->x(0);
164 for (
const auto detID : inputWorkspace->getSpectrum(i).getDetectorIDs()) {
165 const auto index = detectorInfo.indexOf(detID);
166 if (!detectorInfo.isMasked(
index))
167 angle += detectorInfo.detector(
index).solidAngle(samplePos);
172 std::vector<double> maskFractions;
173 if (inputWorkspace->hasMaskedBins(i)) {
176 maskFractions.resize(numBins, 1.0);
177 MatrixWorkspace::MaskList::const_iterator it, itEnd(mask.end());
178 for (it = mask.begin(); it != itEnd; ++it) {
181 maskFractions[it->first] -= it->second;
184 double maskFraction(1);
193 for (
int j =
static_cast<int>(numBins) - 1; j >=
static_cast<int>(wavStart); --j) {
197 const double binWidth =
X[j + 1] -
X[j];
199 const double wavLength =
X[j] + (binWidth) / 2.0;
208 const double Q = 4.0 * M_PI * sinTheta / wavLength;
211 const double Qx = a * Q;
213 if (Qx < axis.front() || Qx >= axis.back())
215 const double Qy = b * Q;
216 if (Qy < axis.front() || Qy >= axis.back())
220 const auto xIndex = std::upper_bound(axis.begin(), axis.end(), Qx) - axis.begin() - 1;
221 const int yIndex =
static_cast<int>(std::upper_bound(axis.begin(), axis.end(), Qy) - axis.begin() - 1);
226 double &outputBinY = outputWorkspace->mutableY(yIndex)[xIndex];
227 double &outputBinE = outputWorkspace->mutableE(yIndex)[xIndex];
229 if (std::isnan(outputBinY)) {
230 outputBinY = outputBinE = 0;
235 outputBinE = std::sqrt((outputBinE * outputBinE) + (E[j] * E[j]));
238 if (!maskFractions.empty()) {
239 maskFraction = maskFractions[j];
248 weight = maskFraction * angle;
250 weight = maskFraction;
254 auto &outWeightsY = weights->mutableY(yIndex);
255 auto &outWeightsE = weights->mutableE(yIndex);
257 if (pixelAdj && waveAdj) {
258 auto pixelY = pixelAdj->y(i)[0];
259 auto pixelE = pixelAdj->e(i)[0];
261 auto waveY = waveAdj->y(0)[j];
262 auto waveE = waveAdj->e(0)[j];
264 outWeightsY[xIndex] += weight * pixelY * waveY;
265 const double pixelYSq = pixelY * pixelY;
266 const double pixelESq = pixelE * pixelE;
267 const double waveYSq = waveY * waveY;
268 const double waveESq = waveE * waveE;
271 outWeightsE[xIndex] += weight * weight * (waveESq * pixelYSq + pixelESq * waveYSq);
272 }
else if (pixelAdj) {
273 auto pixelY = pixelAdj->y(i)[0];
274 auto pixelE = pixelAdj->e(i)[0];
276 outWeightsY[xIndex] += weight * pixelY;
277 const double pixelESq = weight * pixelE;
279 outWeightsE[xIndex] += pixelESq * pixelESq;
280 }
else if (waveAdj) {
281 auto waveY = waveAdj->y(0)[j];
282 auto waveE = waveAdj->e(0)[j];
284 outWeightsY[xIndex] += weight * waveY;
285 const double waveESq = weight * waveE;
287 outWeightsE[xIndex] += waveESq * waveESq;
289 outWeightsY[xIndex] += weight;
293 prog.
report(
"Calculating Q");
299 size_t numHist = weights->getNumberHistograms();
300 for (
size_t i = 0; i < numHist; i++) {
301 auto &weightsE = weights->mutableE(i);
302 std::transform(weightsE.cbegin(), weightsE.cend(), weightsE.begin(), [&](
double val) { return sqrt(val); });
309 for (
size_t i = 0; i < ws_sumOfCounts->getNumberHistograms(); i++) {
310 ws_sumOfCounts->setHistogram(i, outputWorkspace->histogram(i));
317 outputWorkspace /= weights;
318 outputWorkspace->setDistribution(
true);
321 const size_t nhist = outputWorkspace->getNumberHistograms();
322 const size_t nbins = outputWorkspace->blocksize();
324 for (
size_t i = 0; i < nhist; ++i) {
325 auto &yOut = outputWorkspace->y(i);
326 for (
size_t j = 0; j < nbins; ++j) {
327 if (yOut[j] < 1.0e-12)
333 g_log.
notice() <<
"There are a total of " << emptyBins <<
" (" << (100 * emptyBins) / (outputWorkspace->size())
334 <<
"%) empty Q bins.\n";
350 if (min < 0 || max < 0)
351 std::cerr <<
"Only positive numbers allowed\n";
354 std::vector<double> outBins(num);
357 double binWidth =
fabs((max - min) / (num - 1));
358 for (
int i = 0; i < num; ++i) {
359 outBins[i] = pow(10, min + i * binWidth);
371 g_log.
warning() <<
"Could not retrieve Qmin from run object.\n";
373 g_log.
notice() <<
"QxQy: Using logarithm binning with qmin=" << qmin <<
".\n";
384 const bool log_binning =
getProperty(
"IQxQyLogBinning");
387 auto nBins =
static_cast<int>(max /
delta);
389 HistogramData::BinEdges axis;
393 std::vector<double> positiveBinning =
logBinning(qmin, max, nBins);
394 std::reverse(std::begin(positiveBinning), std::end(positiveBinning));
395 std::vector<double> totalBinning = positiveBinning;
396 std::for_each(std::begin(totalBinning), std::end(totalBinning), [](
double &
n) {
n = -1 *
n; });
397 totalBinning.emplace_back(0.0);
398 std::reverse(std::begin(positiveBinning), std::end(positiveBinning));
399 totalBinning.insert(std::end(totalBinning), std::begin(positiveBinning), std::end(positiveBinning));
400 nBins = nBins * 2 + 1;
404 if (nBins *
delta != max)
407 double startVal = -1.0 *
delta * nBins;
412 HistogramData::BinEdges linearAxis(nBins, HistogramData::LinearGenerator(startVal,
delta));
421 auto verticalAxis = std::make_unique<BinEdgeAxis>(nBins);
422 auto verticalAxisRaw = verticalAxis.get();
423 outputWorkspace->replaceAxis(1, std::move(verticalAxis));
424 for (
int i = 0; i < nBins; ++i) {
425 const double currentVal = axis[i];
427 verticalAxisRaw->setValue(i, currentVal);
431 for (
int i = 0; i < nBins - 1; ++i) {
432 outputWorkspace->setBinEdges(i, axis);
433 auto &
y = outputWorkspace->mutableY(i);
434 auto &e = outputWorkspace->mutableE(i);
436 for (
int j = 0; j < nBins - j; ++j) {
437 y[j] = std::numeric_limits<double>::quiet_NaN();
438 e[j] = std::numeric_limits<double>::quiet_NaN();
443 outputWorkspace->getAxis(1)->unit() = outputWorkspace->getAxis(0)->unit() =
447 outputWorkspace->setYUnitLabel(
"Cross Section (1/cm)");
450 return outputWorkspace;
#define DECLARE_ALGORITHM(classname)
std::map< DeltaEMode::Type, std::string > index
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.
const Run & run() const
Run details object access.
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.
bool hasProperty(const std::string &name) const
Does the property exist on the object.
HeldType getPropertyValueAsType(const std::string &name) const
Get the value of a property as the given TYPE.
Base MatrixWorkspace Abstract Class.
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.
This class stores information regarding an experimental run as a series of log entries.
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 calcComponents(const double wavAngstroms, double &xFrac, double &yFrac) const
Calculate the sins and cosins of angles as required to calculate Q is 2 dimensions.
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.
double getQminFromWs(const API::MatrixWorkspace &inputWorkspace)
API::MatrixWorkspace_sptr setUpOutputWorkspace(const API::MatrixWorkspace_const_sptr &inputWorkspace)
Creates the output workspace, setting the X vector to the bins boundaries in Qx.
void init() override
Initialisation code.
std::vector< double > logBinning(double min, double max, int num)
This function calculates a logarithm binning It gives the same result as scipy:
void exec() override
Execution code.
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 notice(const std::string &msg)
Logs at notice level.
void warning(const std::string &msg)
Logs at warning level.
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
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.
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
@ Input
An input workspace.
@ Output
An output workspace.