Mantid
Loading...
Searching...
No Matches
Qxy.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
4// NScD Oak Ridge National Laboratory, European Spallation Source,
5// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6// SPDX - License - Identifier: GPL - 3.0 +
11#include "MantidAPI/Run.h"
19#include "MantidHistogramData/LinearGenerator.h"
24
25#include <algorithm>
26#include <cmath>
27
28namespace Mantid::Algorithms {
29
30// Register the algorithm into the AlgorithmFactory
32
33using namespace Kernel;
34using namespace API;
35using namespace Geometry;
36
37void Qxy::init() {
38 auto wsValidator = std::make_shared<CompositeValidator>();
39 wsValidator->add<WorkspaceUnitValidator>("Wavelength");
40 wsValidator->add<HistogramValidator>();
41 wsValidator->add<InstrumentValidator>();
42
43 declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input, wsValidator),
44 "The corrected data in units of wavelength.");
45 declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
46 "The name to use for the corrected workspace.");
47
48 auto mustBePositive = std::make_shared<BoundedValidator<double>>();
49 mustBePositive->setLower(1.0e-12);
50
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>();
60 wavVal->add<WorkspaceUnitValidator>("Wavelength");
61 wavVal->add<HistogramValidator>();
63 std::make_unique<WorkspaceProperty<>>("WavelengthAdj", "", Direction::Input, PropertyMode::Optional, wavVal),
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 "
68 "status]].");
69 declareProperty("AccountForGravity", false, "Whether to correct for the effects of gravity.", Direction::Input);
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);
74 declareProperty("RadiusCut", 0.0, mustBePositive2,
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]])");
78 declareProperty("WaveCut", 0.0, mustBePositive2,
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]])");
82 declareProperty("OutputParts", false,
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.");
89}
90
91void Qxy::exec() {
92 MatrixWorkspace_const_sptr inputWorkspace = getProperty("InputWorkspace");
93 MatrixWorkspace_const_sptr waveAdj = getProperty("WavelengthAdj");
94 MatrixWorkspace_const_sptr pixelAdj = getProperty("PixelAdj");
95 const bool doGravity = getProperty("AccountForGravity");
96 const bool doSolidAngle = getProperty("SolidAngleWeighting");
97
98 // throws if we don't have common binning or another incompatibility
99 Qhelper helper;
100 helper.examineInput(inputWorkspace, waveAdj, pixelAdj);
101 g_log.debug() << "All input workspaces were found to be valid\n";
102
103 // Create the output Qx-Qy grid
104 MatrixWorkspace_sptr outputWorkspace = this->setUpOutputWorkspace(inputWorkspace);
105
106 // Will also need an identically-sized workspace to hold the solid angle/time
107 // bin masked weight
108 MatrixWorkspace_sptr weights = WorkspaceFactory::Instance().create(outputWorkspace);
109 // Copy the X values from the output workspace to the solidAngles one
110 for (size_t i = 0; i < weights->getNumberHistograms(); ++i)
111 weights->setSharedX(i, outputWorkspace->sharedX(0));
112
113 const size_t numSpec = inputWorkspace->getNumberHistograms();
114 const size_t numBins = inputWorkspace->blocksize();
115
116 // Set the progress bar (1 update for every one percent increase in progress)
117 Progress prog(this, 0.05, 1.0, numSpec);
118
119 const auto &spectrumInfo = inputWorkspace->spectrumInfo();
120 const auto &detectorInfo = inputWorkspace->detectorInfo();
121
122 // the samplePos is often not (0, 0, 0) because the instruments components are
123 // moved to account for the beam centre
124 const V3D samplePos = spectrumInfo.samplePosition();
125
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";
129 continue;
130 }
131 // If no detector found or if it's masked or a monitor, skip onto the next
132 // spectrum
133 if (spectrumInfo.isMonitor(i) || spectrumInfo.isMasked(i))
134 continue;
135
136 // get the bins that are included inside the RadiusCut/WaveCutcut off, those
137 // to calculate for
138 const size_t wavStart =
139 helper.waveLengthCutOff(inputWorkspace, spectrumInfo, getProperty("RadiusCut"), getProperty("WaveCut"), i);
140 if (wavStart >= inputWorkspace->y(i).size()) {
141 // all the spectra in this detector are out of range
142 continue;
143 }
144
145 V3D detPos = spectrumInfo.position(i) - samplePos;
146
147 // these will be re-calculated if gravity is on but without gravity there is
148 // no need
149 double phi = atan2(detPos.Y(), detPos.X());
150 double a = cos(phi);
151 double b = sin(phi);
152 double sinTheta = sin(spectrumInfo.twoTheta(i) * 0.5);
153
154 // Get references to the data for this spectrum
155 const auto &X = inputWorkspace->x(i);
156 const auto &Y = inputWorkspace->y(i);
157 const auto &E = inputWorkspace->e(i);
158
159 const auto &axis = outputWorkspace->x(0);
160
161 // the solid angle of the detector as seen by the sample is used for
162 // normalisation later on
163 double angle = 0.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);
168 }
169
170 // some bins are masked completely or partially, the following vector will
171 // contain the fractions
172 std::vector<double> maskFractions;
173 if (inputWorkspace->hasMaskedBins(i)) {
174 // go through the set and convert it to a vector
175 const MatrixWorkspace::MaskList &mask = inputWorkspace->maskedBins(i);
176 maskFractions.resize(numBins, 1.0);
177 MatrixWorkspace::MaskList::const_iterator it, itEnd(mask.end());
178 for (it = mask.begin(); it != itEnd; ++it) {
179 // The weight for this masked bin is 1 minus the degree to which this
180 // bin is masked
181 maskFractions[it->first] -= it->second;
182 }
183 }
184 double maskFraction(1);
185
186 // this object is not used if gravity correction is off, but it is only
187 // constructed once per spectrum
189 if (doGravity) {
190 grav = GravitySANSHelper(spectrumInfo, i, getProperty("ExtraLength"));
191 }
192
193 for (int j = static_cast<int>(numBins) - 1; j >= static_cast<int>(wavStart); --j) {
194 if (j < 0)
195 break; // Be careful with counting down. Need a better fix but this will
196 // work for now
197 const double binWidth = X[j + 1] - X[j];
198 // Calculate the wavelength at the mid-point of this bin
199 const double wavLength = X[j] + (binWidth) / 2.0;
200
201 if (doGravity) {
202 // SANS instruments must have their y-axis pointing up, show the
203 // detector position as where the neutron would be without gravity
204 sinTheta = grav.calcComponents(wavLength, a, b);
205 }
206
207 // Calculate |Q| for this bin
208 const double Q = 4.0 * M_PI * sinTheta / wavLength;
209
210 // Now get the x & y components of Q.
211 const double Qx = a * Q;
212 // Test whether they're in range, if not go to next spectrum.
213 if (Qx < axis.front() || Qx >= axis.back())
214 break;
215 const double Qy = b * Q;
216 if (Qy < axis.front() || Qy >= axis.back())
217 break;
218 // Find the indices pointing to the place in the 2D array where this bin's
219 // contents should go
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);
222 // PARALLEL_CRITICAL(qxy) /* Write to shared memory - must protect
223 // */
224 {
225 // the data will be copied to this bin in the output array
226 double &outputBinY = outputWorkspace->mutableY(yIndex)[xIndex];
227 double &outputBinE = outputWorkspace->mutableE(yIndex)[xIndex];
228
229 if (std::isnan(outputBinY)) {
230 outputBinY = outputBinE = 0;
231 }
232 // Add the contents of the current bin to the 2D array.
233 outputBinY += Y[j];
234 // add the errors in quadranture
235 outputBinE = std::sqrt((outputBinE * outputBinE) + (E[j] * E[j]));
236
237 // account for masked bins
238 if (!maskFractions.empty()) {
239 maskFraction = maskFractions[j];
240 }
241 // add the total weight for this bin in the weights workspace,
242 // in an equivalent bin to where the data was stored
243
244 // first take into account the product of contributions to the weight
245 // which have no errors
246 double weight = 0.0;
247 if (doSolidAngle)
248 weight = maskFraction * angle;
249 else
250 weight = maskFraction;
251
252 // then the product of contributions which have errors, i.e. optional
253 // pixelAdj and waveAdj contributions
254 auto &outWeightsY = weights->mutableY(yIndex);
255 auto &outWeightsE = weights->mutableE(yIndex);
256
257 if (pixelAdj && waveAdj) {
258 auto pixelY = pixelAdj->y(i)[0];
259 auto pixelE = pixelAdj->e(i)[0];
260
261 auto waveY = waveAdj->y(0)[j];
262 auto waveE = waveAdj->e(0)[j];
263
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;
269 // add product of errors from pixelAdj and waveAdj (note no error on
270 // weight is assumed)
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];
275
276 outWeightsY[xIndex] += weight * pixelY;
277 const double pixelESq = weight * pixelE;
278 // add error from pixelAdj
279 outWeightsE[xIndex] += pixelESq * pixelESq;
280 } else if (waveAdj) {
281 auto waveY = waveAdj->y(0)[j];
282 auto waveE = waveAdj->e(0)[j];
283
284 outWeightsY[xIndex] += weight * waveY;
285 const double waveESq = weight * waveE;
286 // add error from waveAdj
287 outWeightsE[xIndex] += waveESq * waveESq;
288 } else
289 outWeightsY[xIndex] += weight;
290 }
291 } // loop over single spectrum
292
293 prog.report("Calculating Q");
294
295 } // loop over all spectra
296
297 // take sqrt of error weight values
298 // left to be executed here for computational efficiency
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); });
303 }
304
305 bool doOutputParts = getProperty("OutputParts");
306 if (doOutputParts) {
307 // copy outputworkspace before it gets further modified
308 MatrixWorkspace_sptr ws_sumOfCounts = WorkspaceFactory::Instance().create(outputWorkspace);
309 for (size_t i = 0; i < ws_sumOfCounts->getNumberHistograms(); i++) {
310 ws_sumOfCounts->setHistogram(i, outputWorkspace->histogram(i));
311 }
312
313 helper.outputParts(this, ws_sumOfCounts, weights);
314 }
315
316 // Divide the output data by the solid angles
317 outputWorkspace /= weights;
318 outputWorkspace->setDistribution(true);
319
320 // Count of the number of empty cells
321 const size_t nhist = outputWorkspace->getNumberHistograms();
322 const size_t nbins = outputWorkspace->blocksize();
323 int emptyBins = 0;
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)
328 ++emptyBins;
329 }
330 }
331
332 // Log the number of empty bins
333 g_log.notice() << "There are a total of " << emptyBins << " (" << (100 * emptyBins) / (outputWorkspace->size())
334 << "%) empty Q bins.\n";
335}
336
349std::vector<double> Qxy::logBinning(double min, double max, int num) {
350 if (min < 0 || max < 0)
351 std::cerr << "Only positive numbers allowed\n";
352 if (min == 0)
353 min = 1e-3; // This is Qmin default! Might have to change this
354 std::vector<double> outBins(num);
355 min = log10(min);
356 max = log10(max);
357 double binWidth = fabs((max - min) / (num - 1));
358 for (int i = 0; i < num; ++i) {
359 outBins[i] = pow(10, min + i * binWidth);
360 }
361 return outBins;
362}
363
364double Qxy::getQminFromWs(const API::MatrixWorkspace &inputWorkspace) {
365 // get qmin from the run properties
366 double qmin = 0;
367 const API::Run &run = inputWorkspace.run();
368 if (run.hasProperty("qmin")) {
369 qmin = run.getPropertyValueAsType<double>("Qmin");
370 } else {
371 g_log.warning() << "Could not retrieve Qmin from run object.\n";
372 }
373 g_log.notice() << "QxQy: Using logarithm binning with qmin=" << qmin << ".\n";
374 return qmin;
375}
376
382 const double max = getProperty("MaxQxy");
383 const double delta = getProperty("DeltaQ");
384 const bool log_binning = getProperty("IQxQyLogBinning");
385
386 // number of bins
387 auto nBins = static_cast<int>(max / delta);
388
389 HistogramData::BinEdges axis;
390 if (log_binning) {
391 // get qmin from the run properties
392 double qmin = getQminFromWs(*inputWorkspace);
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;
401 axis = totalBinning;
402
403 } else {
404 if (nBins * delta != max)
405 ++nBins; // Stop at first boundary past MaxQxy if max is not a multiple of
406 // delta
407 double startVal = -1.0 * delta * nBins;
408 nBins *= 2; // go from -max to +max
409 nBins += 1; // Add 1 - this is a histogram
410
411 // Build up the X values
412 HistogramData::BinEdges linearAxis(nBins, HistogramData::LinearGenerator(startVal, delta));
413 axis = linearAxis;
414 }
415
416 // Create an output workspace with the same meta-data as the input
417 MatrixWorkspace_sptr outputWorkspace =
418 WorkspaceFactory::Instance().create(inputWorkspace, nBins - 1, nBins, nBins - 1);
419
420 // Create a numeric axis to replace the vertical one
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];
426 // Set the Y value on the axis
427 verticalAxisRaw->setValue(i, currentVal);
428 }
429
430 // Fill the X vectors in the output workspace
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);
435
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();
439 }
440 }
441
442 // Set the axis units
443 outputWorkspace->getAxis(1)->unit() = outputWorkspace->getAxis(0)->unit() =
444 UnitFactory::Instance().create("MomentumTransfer");
445 // Set the 'Y' unit (gets confusing here...this is probably a Z axis in this
446 // case)
447 outputWorkspace->setYUnitLabel("Cross Section (1/cm)");
448
449 setProperty("OutputWorkspace", outputWorkspace);
450 return outputWorkspace;
451}
452
453} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
std::map< DeltaEMode::Type, std::string > index
Definition: DeltaEMode.cpp:19
#define fabs(x)
Definition: Matrix.cpp:22
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
Definition: Algorithm.cpp:1913
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Definition: Algorithm.cpp:2076
Kernel::Logger & g_log
Definition: Algorithm.h:451
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.
Definition: LogManager.cpp:265
HeldType getPropertyValueAsType(const std::string &name) const
Get the value of a property as the given TYPE.
Definition: LogManager.cpp:332
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.
Definition: Progress.h:25
This class stores information regarding an experimental run as a series of log entries.
Definition: Run.h:38
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.
Definition: Qhelper.h:20
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...
Definition: Qhelper.cpp:150
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.
Definition: Qhelper.cpp:183
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.
Definition: Qhelper.cpp:29
double getQminFromWs(const API::MatrixWorkspace &inputWorkspace)
Definition: Qxy.cpp:364
API::MatrixWorkspace_sptr setUpOutputWorkspace(const API::MatrixWorkspace_const_sptr &inputWorkspace)
Creates the output workspace, setting the X vector to the bins boundaries in Qx.
Definition: Qxy.cpp:381
void init() override
Initialisation code.
Definition: Qxy.cpp:37
std::vector< double > logBinning(double min, double max, int num)
This function calculates a logarithm binning It gives the same result as scipy:
Definition: Qxy.cpp:349
void exec() override
Execution code.
Definition: Qxy.cpp:91
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.
Definition: Logger.cpp:114
void notice(const std::string &msg)
Logs at notice level.
Definition: Logger.cpp:95
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
Definition: ProgressBase.h:51
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
Class for 3D vectors.
Definition: V3D.h:34
constexpr double X() const noexcept
Get x.
Definition: V3D.h:232
constexpr double Y() const noexcept
Get y.
Definition: V3D.h:233
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.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54