Mantid
Loading...
Searching...
No Matches
Q1D2.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 +
7#include <utility>
8
9#include "MantidAPI/Axis.h"
12#include "MantidAPI/ISpectrum.h"
25#include "MantidIndexing/IndexInfo.h"
32#include "MantidParallel/Communicator.h"
33#include "MantidTypes/SpectrumDefinition.h"
34
35namespace Mantid::Algorithms {
36
37// Register the algorithm into the AlgorithmFactory
39
40using namespace Kernel;
41using namespace API;
42using namespace Geometry;
43using namespace DataObjects;
44
45Q1D2::Q1D2() : API::Algorithm(), m_dataWS(), m_doSolidAngle(false) {}
46
47void Q1D2::init() {
48 auto dataVal = std::make_shared<CompositeValidator>();
49 dataVal->add<WorkspaceUnitValidator>("Wavelength");
50 dataVal->add<HistogramValidator>();
51 dataVal->add<InstrumentValidator>();
52 dataVal->add<CommonBinsValidator>();
53 declareProperty(std::make_unique<WorkspaceProperty<>>("DetBankWorkspace", "", Direction::Input, dataVal),
54 "Particle counts as a function of wavelength");
55 declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
56 "Name of the workspace that will contain the result of the calculation");
57 declareProperty(std::make_unique<ArrayProperty<double>>("OutputBinning", std::make_shared<RebinParamsValidator>()),
58 "A comma separated list of first bin boundary, width, last bin boundary. "
59 "Optionally\n"
60 "this can be followed by a comma and more widths and last boundary "
61 "pairs.\n"
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>();
67 wavVal->add<WorkspaceUnitValidator>("Wavelength");
68 wavVal->add<HistogramValidator>();
70 std::make_unique<WorkspaceProperty<>>("WavelengthAdj", "", Direction::Input, PropertyMode::Optional, wavVal),
71 "Scaling to apply to each bin.\n"
72 "Must have the same number of bins as the DetBankWorkspace");
74 std::make_unique<WorkspaceProperty<>>("WavePixelAdj", "", Direction::Input, PropertyMode::Optional, dataVal),
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);
82
83 declareProperty("RadiusCut", 0.0, mustBePositive,
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.");
89
90 declareProperty("WaveCut", 0.0, mustBePositive,
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.");
96
97 declareProperty("OutputParts", false,
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.");
105
107 std::make_unique<WorkspaceProperty<>>("QResolution", "", Direction::Input, PropertyMode::Optional, dataVal),
108 "Workspace to calculate the Q resolution.\n");
109}
114 m_dataWS = getProperty("DetBankWorkspace");
115 MatrixWorkspace_const_sptr waveAdj = getProperty("WavelengthAdj");
116 MatrixWorkspace_const_sptr pixelAdj = getProperty("PixelAdj");
117 MatrixWorkspace_const_sptr wavePixelAdj = getProperty("WavePixelAdj");
118 MatrixWorkspace_const_sptr qResolution = getProperty("QResolution");
119
120 const bool doGravity = getProperty("AccountForGravity");
121 m_doSolidAngle = getProperty("SolidAngleWeighting");
122
123 // throws if we don't have common binning or another incompatibility
124 Qhelper helper;
125 helper.examineInput(m_dataWS, waveAdj, pixelAdj, qResolution);
126 // FIXME: how to examine the wavePixelAdj?
127 g_log.debug() << "All input workspaces were found to be valid\n";
128 // normalization as a function of wavelength (i.e. centers of x-value bins)
129 double const *const binNorms = waveAdj ? &(waveAdj->y(0)[0]) : nullptr;
130 // error on the wavelength normalization
131 double const *const binNormEs = waveAdj ? &(waveAdj->e(0)[0]) : nullptr;
132
133 // define the (large number of) data objects that are going to be used in all
134 // iterations of the loop below
135
136 // Flag to decide if Q Resolution is to be used
137 auto useQResolution = static_cast<bool>(qResolution);
138
139 // this will become the output workspace from this algorithm
140 MatrixWorkspace_sptr outputWS = setUpOutputWorkspace(getProperty("OutputBinning"));
141
142 auto &QOut = outputWS->x(0);
143 auto &YOut = outputWS->mutableY(0);
144 auto &EOutTo2 = outputWS->mutableE(0);
145 // normalisation that is applied to counts in each Q bin
146 HistogramData::HistogramY normSum(YOut.size(), 0.0);
147 // the error on the normalisation
148 HistogramData::HistogramE normError2(EOutTo2.size(), 0.0);
149 // the averaged Q resolution.
150 HistogramData::HistogramDx qResolutionOut(YOut.size(), 0.0);
151
152 const auto numSpec = static_cast<int>(m_dataWS->getNumberHistograms());
153 Progress progress(this, 0.05, 1.0, numSpec + 1);
154
155 const auto &spectrumInfo = m_dataWS->spectrumInfo();
156 PARALLEL_FOR_IF(Kernel::threadSafe(*m_dataWS, *outputWS, pixelAdj.get()))
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";
162 continue;
163 }
164 // Skip if we have a monitor or if the detector is masked.
165 if (spectrumInfo.isMonitor(i) || spectrumInfo.isMasked(i))
166 continue;
167
168 // get the bins that are included inside the RadiusCut/WaveCutcut off, those
169 // to calculate for
170 // const size_t wavStart = waveLengthCutOff(i);
171 const size_t wavStart =
172 helper.waveLengthCutOff(m_dataWS, spectrumInfo, getProperty("RadiusCut"), getProperty("WaveCut"), i);
173 if (wavStart >= m_dataWS->y(i).size()) {
174 // all the spectra in this detector are out of range
175 continue;
176 }
177
178 const size_t numWavbins = m_dataWS->y(i).size() - wavStart;
179 // make just one call to new to reduce CPU overhead on each thread, access
180 // to these
181 // three "arrays" is via iterators
182 HistogramData::HistogramY _noDirectUseStorage_(3 * numWavbins);
183 // normalization term
184 auto norms = _noDirectUseStorage_.begin();
185 // the error on these weights, it contributes to the error calculation on
186 // the output workspace
187 auto normETo2s = norms + numWavbins;
188 // the Q values calculated from input wavelength workspace
189 auto QIn = normETo2s + numWavbins;
190
191 // the weighting for this input spectrum that is added to the normalization
192 calculateNormalization(wavStart, i, pixelAdj, wavePixelAdj, binNorms, binNormEs, norms, normETo2s);
193
194 // now read the data from the input workspace, calculate Q for each bin
195 convertWavetoQ(spectrumInfo, i, doGravity, wavStart, QIn, getProperty("ExtraLength"));
196
197 // Pointers to the counts data and it's error
198 auto YIn = m_dataWS->y(i).cbegin() + wavStart;
199 auto EIn = m_dataWS->e(i).cbegin() + wavStart;
200
201 // Pointers to the QResolution data. Note that the xdata was initially the
202 // same, hence
203 // the same indexing applies to the y values of m_dataWS and qResolution
204 // If we want to use it set it to the correct value, else to YIN, although
205 // that does not matter, as
206 // we won't use it
207 auto QResIn = useQResolution ? (qResolution->y(i).cbegin() + wavStart) : YIn;
208
209 // when finding the output Q bin remember that the input Q bins (from the
210 // convert to wavelength) start high and reduce
211 auto loc = QOut.cend();
212 // sum the Q contributions from each individual spectrum into the output
213 // array
214 const auto end = m_dataWS->y(i).cend();
215 for (; YIn != end; ++YIn, ++EIn, ++QIn, ++norms, ++normETo2s) {
216 // find the output bin that each input y-value will fall into, remembering
217 // there is one more bin boundary than bins
218 getQBinPlus1(QOut, *QIn, loc);
219 // ignore counts that are out of the output range
220 if ((loc != QOut.begin()) && (loc != QOut.end())) {
221 // the actual Q-bin to add something to
222 const size_t bin = loc - QOut.begin() - 1;
223 PARALLEL_CRITICAL(q1d_counts_sum) {
224 YOut[bin] += *YIn;
225 normSum[bin] += *norms;
226 // these are the errors squared which will be summed and square rooted
227 // at the end
228 EOutTo2[bin] += (*EIn) * (*EIn);
229 normError2[bin] += *normETo2s;
230 if (useQResolution) {
231 auto QBin = (QOut[bin + 1] - QOut[bin]);
232 // Here we need to take into account the Bin width and the count
233 // weigthing. The
234 // formula should be YIN* sqrt(QResIn^2 + (QBin/sqrt(12))^2)
235 qResolutionOut[bin] += (*YIn) * std::sqrt((*QResIn) * (*QResIn) + QBin * QBin / 12.0);
236 }
237 }
238 }
239
240 // Increment the QResolution iterator
241 if (useQResolution) {
242 ++QResIn;
243 }
244 }
245
246 PARALLEL_CRITICAL(q1d_spectra_map) {
247 progress.report("Computing I(Q)");
248 // Add up the detector IDs in the output spectrum at workspace index 0
249 const auto &inSpec = m_dataWS->getSpectrum(i);
250 auto &outSpec = outputWS->getSpectrum(0);
251 outSpec.addDetectorIDs(inSpec.getDetectorIDs());
252 }
253
255 }
257
258 if (communicator().size() > 1) {
259 int tag = 0;
260 auto size = static_cast<int>(YOut.size());
261 if (communicator().rank() == 0) {
262 for (int rank = 1; rank < communicator().size(); ++rank) {
263 HistogramData::HistogramY y(YOut.size());
264 HistogramData::HistogramE e2(YOut.size());
265 communicator().recv(rank, tag, &y[0], size);
266 YOut += y;
267 communicator().recv(rank, tag, &e2[0], size);
268 EOutTo2 += e2;
269 communicator().recv(rank, tag, &y[0], size);
270 normSum += y;
271 communicator().recv(rank, tag, &e2[0], size);
272 normError2 += e2;
273 int detCount;
274 communicator().recv(rank, tag, detCount);
275 std::vector<detid_t> detIds(detCount);
276 communicator().recv(rank, tag, detIds.data(), detCount);
277 outputWS->getSpectrum(0).addDetectorIDs(detIds);
278 }
279 } else {
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());
287 communicator().send(0, tag, nDets);
288 communicator().send(0, tag, detIds.data(), nDets);
289 }
290 }
291
292 if (useQResolution) {
293 // The number of Q (x)_ values is N, while the number of DeltaQ values is
294 // N-1,
295 // Richard Heenan suggested to duplicate the last entry of DeltaQ
296 auto countsIterator = YOut.cbegin();
297 auto qResolutionIterator = qResolutionOut.begin();
298 for (; countsIterator != YOut.end(); ++countsIterator, ++qResolutionIterator) {
299 // Divide by the counts of the Qbin, if the counts are 0, the the
300 // qresolution will also be 0
301 if ((*countsIterator) > 0.0) {
302 *qResolutionIterator = (*qResolutionIterator) / (*countsIterator);
303 }
304 }
305 outputWS->setPointStandardDeviations(0, std::move(qResolutionOut));
306 }
307
308 bool doOutputParts = getProperty("OutputParts");
309 if (doOutputParts && communicator().rank() == 0) {
310 MatrixWorkspace_sptr ws_sumOfCounts = WorkspaceFactory::Instance().create(outputWS);
311 ws_sumOfCounts->setSharedX(0, outputWS->sharedX(0));
312 // Copy now as YOut is modified in normalize
313 ws_sumOfCounts->mutableY(0) = YOut;
314 ws_sumOfCounts->setSharedDx(0, outputWS->sharedDx(0));
315 ws_sumOfCounts->setFrequencyVariances(0, outputWS->e(0));
316
317 MatrixWorkspace_sptr ws_sumOfNormFactors = WorkspaceFactory::Instance().create(outputWS);
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);
322
323 helper.outputParts(this, ws_sumOfCounts, ws_sumOfNormFactors);
324 } else if (doOutputParts) {
325 helper.outputParts(this, nullptr, nullptr);
326 }
327
328 progress.report("Normalizing I(Q)");
329 // finally divide the number of counts in each output Q bin by its weighting
330 normalize(normSum, normError2, YOut, EOutTo2);
331
332 if (communicator().rank() == 0) {
333 setProperty("OutputWorkspace", outputWS);
334 }
335}
336
342API::MatrixWorkspace_sptr Q1D2::setUpOutputWorkspace(const std::vector<double> &binParams) const {
343 // Calculate the output binning
344 HistogramData::BinEdges XOut(0);
345 static_cast<void>(VectorHelper::createAxisFromRebinParams(binParams, XOut.mutableRawData()));
346
347 // Create output workspace. On all but rank 0 this is a temporary workspace.
348 Indexing::IndexInfo indexInfo(
349 1, communicator().rank() == 0 ? Parallel::StorageMode::MasterOnly : Parallel::StorageMode::Cloned,
350 communicator());
351 indexInfo.setSpectrumDefinitions(std::vector<SpectrumDefinition>(1));
352 auto outputWS = create<MatrixWorkspace>(*m_dataWS, indexInfo, XOut);
353 outputWS->getAxis(0)->unit() = UnitFactory::Instance().create("MomentumTransfer");
354 outputWS->setYUnitLabel("1/cm");
355
356 outputWS->setDistribution(true);
357
358 outputWS->getSpectrum(0).clearDetectorIDs();
359 outputWS->getSpectrum(0).setSpectrumNo(1);
360
361 return outputWS;
362}
363
381void Q1D2::calculateNormalization(const size_t wavStart, const size_t wsIndex,
382 const API::MatrixWorkspace_const_sptr &pixelAdj,
383 const API::MatrixWorkspace_const_sptr &wavePixelAdj, double const *const binNorms,
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);
388 // use that the normalization array ends at the start of the error array
389 for (auto n = norm, e = normETo2; n != normETo2; ++n, ++e) {
390 *n = detectorAdj;
391 *e = detAdjErr * detAdjErr;
392 }
393
394 if (binNorms && binNormEs) {
395 if (wavePixelAdj)
396 // pass the iterator for the wave pixel Adj dependent
397 addWaveAdj(binNorms + wavStart, binNormEs + wavStart, norm, normETo2, wavePixelAdj->y(wsIndex).begin() + wavStart,
398 wavePixelAdj->e(wsIndex).begin() + wavStart);
399 else
400 addWaveAdj(binNorms + wavStart, binNormEs + wavStart, norm, normETo2);
401 }
402 normToMask(wavStart, wsIndex, norm, normETo2);
403}
404
416void Q1D2::pixelWeight(const API::MatrixWorkspace_const_sptr &pixelAdj, const size_t wsIndex, double &weight,
417 double &error) const {
418 const auto &detectorInfo = m_dataWS->detectorInfo();
419 const V3D samplePos = detectorInfo.samplePosition();
420
421 if (m_doSolidAngle) {
422 weight = 0.0;
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);
427 }
428 } else
429 weight = 1.0;
430
431 if (weight < 1e-200) {
432 throw std::logic_error("Invalid (zero or negative) solid angle for one detector");
433 }
434 // this input multiplies up the adjustment if it exists
435 if (pixelAdj) {
436 weight *= pixelAdj->readY(wsIndex)[0];
437 error = weight * pixelAdj->readE(wsIndex)[0];
438 } else {
439 error = 0.0;
440 }
441}
442
455void Q1D2::addWaveAdj(const double *c, const double *Dc, HistogramData::HistogramY::iterator bInOut,
456 HistogramData::HistogramY::iterator e2InOut) const {
457 // normalize by the wavelength dependent correction, keeping the percentage
458 // errors the same
459 // the error when a = b*c, the formula for Da, the error on a, in terms of Db,
460 // etc. is (Da/a)^2 = (Db/b)^2 + (Dc/c)^2
461 //(Da)^2 = ((Db*a/b)^2 + (Dc*a/c)^2) = (Db*c)^2 + (Dc*b)^2
462 // the variable names relate to those above as: existing values (b=bInOut)
463 // multiplied by the additional errors (Dc=binNormEs), existing errors
464 // (Db=sqrt(e2InOut)) times new factor (c=binNorms)
465
466 // use the fact that error array follows straight after the normalization
467 // array
468 const auto end = e2InOut;
469 for (; bInOut != end; ++e2InOut, ++c, ++Dc, ++bInOut) {
470 // first the error
471 *e2InOut = ((*e2InOut) * (*c) * (*c)) + ((*Dc) * (*Dc) * (*bInOut) * (*bInOut));
472 // now the actual calculation a = b*c
473 *bInOut = (*bInOut) * (*c);
474 }
475}
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 {
496 // normalize by the wavelength dependent correction, keeping the percentage
497 // errors the same
498 // the error when a = b*c*e, the formula for Da, the error on a, in terms of
499 // Db, etc. is
500 // (Da/a)^2 = (Db/b)^2 + (Dc/c)^2 + (De/e)^2
501 //(Da)^2 = ((Db*a/b)^2 + (Dc*a/c)^2) + (De * a/e)^2
502 // But: a/b = c*e; a/c = b*e; a/e = b*c;
503 // So:
504 // (Da)^2 = (c*e*Db)^2 + (b*e*Dc)^2 + (b*c*De)^2
505 //
506 // Consider:
507 // Da = Error (e2InOut)
508 // Db^2 = PixelDependentError (e2InOut)
509 // b = PixelDependentValue (bInOut)
510 // c = WaveDependentValue (c)
511 // Dc = WaveDependentError (Dc)
512 // e = PixelWaveDependentValue (wavePixelAdjData)
513 // De = PiexlWaveDependentError (wavePixelAdjError)
514
515 // use the fact that error array follows straight after the normalization
516 // array
517 const auto end = e2InOut;
518 for (; bInOut != end; ++e2InOut, ++c, ++Dc, ++bInOut, ++wavePixelAdjData, ++wavePixelAdjError) {
519 // first the error
520 *e2InOut = ((*e2InOut) * (*c) * (*c) * (*wavePixelAdjData) * (*wavePixelAdjData)) +
521 ((*Dc) * (*Dc) * (*bInOut) * (*bInOut) * (*wavePixelAdjData) * (*wavePixelAdjData)) +
522 ((*wavePixelAdjError) * (*wavePixelAdjError) * (*c) * (*c) * (*bInOut) * (*bInOut));
523 // now the actual calculation a = b*c*e : Pixel * Wave * PixelWave
524 *bInOut = (*bInOut) * (*c) * (*wavePixelAdjData);
525 }
526}
527
537void Q1D2::normToMask(const size_t offSet, const size_t wsIndex, const HistogramData::HistogramY::iterator theNorms,
538 const HistogramData::HistogramY::iterator errorSquared) const {
539 // if any bins are masked it is normally a small proportion
540 if (m_dataWS->hasMaskedBins(wsIndex)) {
541 // Get a reference to the list of masked bins
542 const MatrixWorkspace::MaskList &mask = m_dataWS->maskedBins(wsIndex);
543 // Now iterate over the list, adjusting the weights for the affected bins
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) {
548 // this masked bin wasn't in the range being delt with anyway
549 continue;
550 }
551 outBin -= offSet;
552 // The weight for this masked bin is 1 - the degree to which this bin is
553 // masked
554 const double factor = 1.0 - (it->second);
555 *(theNorms + outBin) *= factor;
556 *(errorSquared + outBin) *= factor * factor;
557 }
558 }
559}
560
574void Q1D2::convertWavetoQ(const SpectrumInfo &spectrumInfo, const size_t wsInd, const bool doGravity,
575 const size_t offset, HistogramData::HistogramY::iterator Qs, const double extraLength) const {
576 static const double FOUR_PI = 4.0 * M_PI;
577
578 // wavelengths (lamda) to be converted to Q
579 auto waves = m_dataWS->x(wsInd).cbegin() + offset;
580 // going from bin boundaries to bin centered x-values the size goes down one
581 const auto end = m_dataWS->x(wsInd).end() - 1;
582 if (doGravity) {
583 GravitySANSHelper grav(spectrumInfo, wsInd, extraLength);
584 for (; waves != end; ++Qs, ++waves) {
585 // the HistogramValidator at the start should ensure that we have one more
586 // bin on the input wavelengths
587 const double lambda = 0.5 * (*(waves + 1) + (*waves));
588 // as the fall under gravity is wavelength dependent sin theta is now
589 // different for each bin with each detector
590 const double sinTheta = grav.calcSinTheta(lambda);
591 // Now we're ready to go to Q
592 *Qs = FOUR_PI * sinTheta / lambda;
593 }
594 } else {
595 // Calculate the Q values for the current spectrum, using Q =
596 // 4*pi*sin(theta)/lambda
597 const double factor = 2.0 * FOUR_PI * sin(spectrumInfo.twoTheta(wsInd) * 0.5);
598 for (; waves != end; ++Qs, ++waves) {
599 // the HistogramValidator at the start should ensure that we have one more
600 // bin on the input wavelengths
601 *Qs = factor / (*(waves + 1) + (*waves));
602 }
603 }
604}
618void Q1D2::getQBinPlus1(const HistogramData::HistogramX &OutQs, const double QToFind,
619 HistogramData::HistogramX::const_iterator &loc) const {
620 if (loc != OutQs.end()) {
621 while (loc != OutQs.begin()) {
622 if ((QToFind >= *(loc - 1)) && (QToFind < *loc)) {
623 return;
624 }
625 --loc;
626 }
627 if (QToFind < *loc) {
628 // QToFind is outside the array leave loc == OutQs.begin()
629 return;
630 }
631 } else // loc == OutQs.end()
632 {
633 if (OutQs.empty() || QToFind > *(loc - 1)) {
634 // outside the array leave loc == OutQs.end()
635 return;
636 }
637 }
638
639 // we are lost, normally the order of the Q values means we only get here on
640 // the first iteration. It's slow
641 loc = std::lower_bound(OutQs.begin(), OutQs.end(), QToFind);
642}
643
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) {
657 // the normalisation is a = b/c where b = counts c =normalistion term
658 const double c = normSum[k];
659 const double a = counts[k] /= c;
660 // when a = b/c, the formula for Da, the error on a, in terms of Db, etc. is
661 // (Da/a)^2 = (Db/b)^2 + (Dc/c)^2
662 //(Da)^2 = ((Db/b)^2 + (Dc/c)^2)*(b^2/c^2) = ((Db/c)^2 + (b*Dc/c^2)^2) =
663 //(Db^2 + (b*Dc/c)^2)/c^2 = (Db^2 + (Dc*a)^2)/c^2
664 // this will work as long as c>0, but then the above formula above can't
665 // deal with 0 either
666 const double aOverc = a / c;
667 errors[k] = std::sqrt(errors[k] / (c * c) + normError2[k] * aOverc * aOverc);
668 }
669}
670
671namespace {
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));
675}
676} // namespace
677
678Parallel::ExecutionMode
679Q1D2::getParallelExecutionMode(const std::map<std::string, Parallel::StorageMode> &storageModes) const {
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"));
685}
686
687} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
const std::vector< double > * lambda
double error
Definition: IndexPeaks.cpp:133
std::map< DeltaEMode::Type, std::string > index
Definition: DeltaEMode.cpp:19
#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...
Definition: MultiThreaded.h:94
#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.
Definition: Algorithm.h:85
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
void progress(double p, const std::string &msg="", double estimatedTime=0.0, int progressPrecision=0)
Sends ProgressNotification.
Definition: Algorithm.cpp:231
const Parallel::Communicator & communicator() const
Returns a const reference to the (MPI) communicator of the algorithm.
Definition: Algorithm.cpp:1870
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.
Definition: Progress.h:25
API::SpectrumInfo is an intermediate step towards a SpectrumInfo that is part of Instrument-2....
Definition: SpectrumInfo.h:53
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.
Definition: Q1D2.cpp:679
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.
Definition: Q1D2.cpp:537
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.
Definition: Q1D2.cpp:381
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...
Definition: Q1D2.cpp:654
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...
Definition: Q1D2.cpp:416
Q1D2()
Default constructor.
Definition: Q1D2.cpp:45
const std::string name() const override
Algorithm's name.
Definition: Q1D2.h:31
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 ...
Definition: Q1D2.cpp:574
API::MatrixWorkspace_sptr setUpOutputWorkspace(const std::vector< double > &binParams) const
Creates the output workspace, its size, units, etc.
Definition: Q1D2.cpp:342
API::MatrixWorkspace_const_sptr m_dataWS
the experimental workspace with counts across the detector
Definition: Q1D2.h:50
void init() override
Initialisation code.
Definition: Q1D2.cpp:47
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...
Definition: Q1D2.cpp:618
void exec() override
Execution code.
Definition: Q1D2.cpp:113
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.
Definition: Q1D2.cpp:455
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
Support for a property that holds an array of values.
Definition: ArrayProperty.h:28
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 warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
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
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 > &params, 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.
Definition: MultiThreaded.h:22
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54