Mantid
Loading...
Searching...
No Matches
NormaliseToMonitor.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 +
20#include "MantidHistogramData/Histogram.h"
26#include "MantidTypes/SpectrumDefinition.h"
27
28#include <cfloat>
29#include <numeric>
30
31using namespace Mantid::API;
32using namespace Mantid::DataObjects;
33using namespace Mantid::HistogramData;
35
36namespace Mantid::Algorithms {
37
38// Method of complex validator class
39// method checks if the property is enabled
41 int sp_id = algo->getProperty(SpectraNum);
42 // if there is spectra number set to normalize by, nothing else can be
43 // selected;
44 if (sp_id > 0) {
45 is_enabled = false;
46 return false;
47 } else {
48 is_enabled = true;
49 }
50
51 // is there the ws property, which describes monitors ws. It also disables the
52 // monitors ID property
54 if (monitorsWS) {
55 is_enabled = false;
56 } else {
57 is_enabled = true;
58 }
59 return is_enabled;
60}
61
62// method checks if other properties have changed and these changes affected
63// MonID property
64bool MonIDPropChanger::isConditionChanged(const IPropertyManager *algo, const std::string &changedPropName) const {
65 UNUSED_ARG(changedPropName);
66 // is enabled is based on other properties:
67 if (!is_enabled)
68 return false;
69 // read monitors list from the input workspace
71 return monitorIdReader(inputWS);
72}
73
74// function which modifies the allowed values for the list of monitors.
76 auto *piProp = dynamic_cast<Kernel::PropertyWithValue<int> *>(pProp);
77 if (!piProp) {
78 throw(std::invalid_argument("modify allowed value has been called on wrong property"));
79 }
80
81 if (iExistingAllowedValues.empty()) {
83 int spectra_max(-1);
84 if (inputWS) {
85 // let's assume that detectors IDs correspond to spectraID --
86 // not always the case but often. TODO: should be fixed
87 spectra_max = static_cast<int>(inputWS->getNumberHistograms()) + 1;
88 }
89 piProp->replaceValidator(std::make_shared<Kernel::BoundedValidator<int>>(-1, spectra_max));
90 } else {
91 piProp->replaceValidator(std::make_shared<Kernel::ListValidator<int>>(iExistingAllowedValues));
92 }
93}
94
95// read the monitors list from the workspace and try to do it once for any
96// particular ws;
98 // no workspace
99 if (!inputWS)
100 return false;
101
102 // no instrument
103 Geometry::Instrument_const_sptr pInstr = inputWS->getInstrument();
104 if (!pInstr)
105 return false;
106
107 // are these monitors really there?
108 std::vector<detid_t> monitorIDList = pInstr->getMonitors();
109 {
110 const auto &specInfo = inputWS->spectrumInfo();
111 std::set<detid_t> idsInWorkspace;
112 size_t i = 0;
113 // Loop over spectra, but finish early if we find everything
114 while (i < specInfo.size() && idsInWorkspace.size() < monitorIDList.size()) {
115 if (specInfo.isMonitor(i))
116 idsInWorkspace.insert(specInfo.detector(i).getID());
117 ++i;
118 }
119 monitorIDList = std::vector<detid_t>(idsInWorkspace.begin(), idsInWorkspace.end());
120 }
121
122 if (monitorIDList.empty()) {
123 if (iExistingAllowedValues.empty()) {
124 return false;
125 } else {
127 return true;
128 }
129 }
130
131 // are known values the same as the values we have just identified?
132 if (iExistingAllowedValues.size() != monitorIDList.size()) {
134 iExistingAllowedValues.assign(monitorIDList.begin(), monitorIDList.end());
135 return true;
136 }
137 // the monitor list has the same size as before. Is it equivalent to the
138 // existing one?
139 bool values_redefined = false;
140 for (size_t i = 0; i < monitorIDList.size(); i++) {
141 if (iExistingAllowedValues[i] != monitorIDList[i]) {
142 values_redefined = true;
143 iExistingAllowedValues[i] = monitorIDList[i];
144 }
145 }
146 return values_redefined;
147}
148
149bool spectrumDefinitionsMatchTimeIndex(const SpectrumDefinition &specDef, const size_t timeIndex) {
150 return std::none_of(specDef.cbegin(), specDef.cend(),
151 [timeIndex](const auto &spec) { return spec.second != timeIndex; });
152}
153
154// Register with the algorithm factory
155DECLARE_ALGORITHM(NormaliseToMonitor)
156
157using namespace Kernel;
158using namespace API;
159using std::size_t;
160
162 // Must be histograms OR one count per bin
163 // Must be raw counts
164 auto validatorHistSingle = std::make_shared<CompositeValidator>(CompositeRelation::OR);
165 validatorHistSingle->add<HistogramValidator>();
166 validatorHistSingle->add<SingleCountValidator>();
167 auto validator = std::make_shared<CompositeValidator>();
168 validator->add(validatorHistSingle);
169 validator->add<RawCountValidator>();
170
171 declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input, validator),
172 "Name of the input workspace. Must be a non-distribution histogram.");
173
174 declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
175 "Name to use for the output workspace");
176 // should be any spectrum number, but named this property MonitorSpectrum to
177 // keep compatibility with previous scripts
178 // Can either set a spectrum within the workspace to be the monitor
179 // spectrum.....
180 declareProperty("MonitorSpectrum", -1,
181 "The spectrum number within the InputWorkspace you want to "
182 "normalize by (It can be a monitor spectrum or a spectrum "
183 "responsible for a group of detectors or monitors)",
185
186 // Or take monitor ID to identify the spectrum one wish to use or
187 declareProperty("MonitorID", -1,
188 "The MonitorID (detector ID), which defines the monitor's data "
189 "within the InputWorkspace. Will be overridden by the values "
190 "correspondent to MonitorSpectrum field if one is provided "
191 "in the field above.\n"
192 "If workspace do not have monitors, the MonitorID can refer "
193 "to empty data and the field then can accepts any MonitorID "
194 "within the InputWorkspace.");
195 // set up the validator, which would verify if spectrum is correct
196 setPropertySettings("MonitorID",
197 std::make_unique<MonIDPropChanger>("InputWorkspace", "MonitorSpectrum", "MonitorWorkspace"));
198
199 // ...or provide it in a separate workspace (note: optional WorkspaceProperty)
200 declareProperty(std::make_unique<WorkspaceProperty<>>("MonitorWorkspace", "", Direction::Input,
201 PropertyMode::Optional, validator),
202 "A workspace containing one or more spectra to normalize the "
203 "InputWorkspace by.");
204 setPropertySettings("MonitorWorkspace", std::make_unique<Kernel::EnabledWhenProperty>("MonitorSpectrum", IS_DEFAULT));
205
206 declareProperty("MonitorWorkspaceIndex", 0,
207 "The index of the spectrum within the MonitorWorkspace(2 "
208 "(0<=ind<=nHistograms in MonitorWorkspace) you want to "
209 "normalize by\n"
210 "(usually related to the index, responsible for the "
211 "monitor's data but can be any).\n"
212 "If no value is provided in this field, '''InputWorkspace''' "
213 "will be normalized by first spectra (with index 0)",
215 setPropertySettings("MonitorWorkspaceIndex",
216 std::make_unique<Kernel::EnabledWhenProperty>("MonitorSpectrum", IS_DEFAULT));
217
218 // If users set either of these optional properties two things happen
219 // 1) normalization is by an integrated count instead of bin-by-bin
220 // 2) if the value is within the range of X's in the spectrum it crops the
221 // spectrum
222 declareProperty("IntegrationRangeMin", EMPTY_DBL(),
223 "If set, normalization will be by integrated count from this "
224 "minimum x value");
225 declareProperty("IntegrationRangeMax", EMPTY_DBL(),
226 "If set, normalization will be by integrated count up to "
227 "this maximum x value");
228 declareProperty("IncludePartialBins", false,
229 "If true and an integration range is set then partial bins at either \n"
230 "end of the integration range are also included");
231
232 declareProperty(std::make_unique<WorkspaceProperty<>>("NormFactorWS", "", Direction::Output, PropertyMode::Optional),
233 "Name of the workspace, containing the normalization factor.\n"
234 "If this name is empty, normalization workspace is not returned. If the "
235 "name coincides with the output workspace name, _normFactor suffix is "
236 "added to this name");
237}
238
240 // Get the input workspace
241 const MatrixWorkspace_sptr inputWS = getProperty("InputWorkspace");
242 MatrixWorkspace_sptr outputWS = getProperty("OutputWorkspace");
243 // First check the inputs
244 checkProperties(inputWS);
245
246 bool isSingleCountWorkspace = false;
247 try {
248 isSingleCountWorkspace = (!inputWS->isHistogramData()) && (inputWS->blocksize() == 1);
249 } catch (std::length_error &) {
250 // inconsistent bin size, not a single count workspace
251 isSingleCountWorkspace = false;
252 }
253
254 // See if the normalization with integration properties are set.
255 const bool integrate = setIntegrationProps(isSingleCountWorkspace);
256
257 if (integrate)
258 normaliseByIntegratedCount(inputWS, outputWS, isSingleCountWorkspace);
259 else
260 normaliseBinByBin(inputWS, outputWS);
261
262 setProperty("OutputWorkspace", outputWS);
263 std::string norm_ws_name = getPropertyValue("NormFactorWS");
264 if (!norm_ws_name.empty()) {
265 std::string out_name = getPropertyValue("OutputWorkspace");
266 if (out_name == norm_ws_name) {
267 // if the normalization factor workspace name coincides with output
268 // workspace name, add _normFactor suffix to this name
269 norm_ws_name = norm_ws_name + "_normFactor";
270 auto pProp = (this->getPointerToProperty("NormFactorWS"));
271 pProp->setValue(norm_ws_name);
272 }
273 if (!integrate)
275 setProperty("NormFactorWS", m_monitor);
276 }
277}
278
286 const std::vector<std::size_t> &workspaceIndexes) {
287 auto childAlg = createChildAlgorithm("ExtractSpectra");
288 childAlg->setProperty<MatrixWorkspace_sptr>("InputWorkspace", ws);
289 childAlg->setProperty("WorkspaceIndexList", workspaceIndexes);
290 childAlg->executeAsChildAlg();
291 MatrixWorkspace_sptr outWS = childAlg->getProperty("OutputWorkspace");
292 return outWS;
293}
294
298std::map<std::string, std::string> NormaliseToMonitor::validateInputs() {
299 std::map<std::string, std::string> issues;
300 // Check where the monitor spectrum should come from
301 const Property *monSpecProp = getProperty("MonitorSpectrum");
302 const Property *monIDProp = getProperty("MonitorID");
303 MatrixWorkspace_const_sptr monWS = getProperty("MonitorWorkspace");
304 // something has to be set
305 if (monSpecProp->isDefault() && !monWS && monIDProp->isDefault()) {
306 const std::string mess("Either MonitorSpectrum, MonitorID or "
307 "MonitorWorkspace has to be provided.");
308 issues["MonitorSpectrum"] = mess;
309 issues["MonitorID"] = mess;
310 issues["MonitorWorkspace"] = mess;
311 }
312
313 const double intMin = getProperty("IntegrationRangeMin");
314 const double intMax = getProperty("IntegrationRangeMax");
315 if (!isEmpty(intMin) && !isEmpty(intMax)) {
316 if (intMin > intMax) {
317 issues["IntegrationRangeMin"] = "Range minimum set to a larger value than maximum.";
318 issues["IntegrationRangeMax"] = "Range maximum set to a smaller value than minimum.";
319 }
320 }
321
322 if (monWS && monSpecProp->isDefault()) {
323 const int monIndex = getProperty("MonitorWorkspaceIndex");
324 if (monIndex < 0) {
325 issues["MonitorWorkspaceIndex"] = "A workspace index cannot be negative.";
326 } else if (monWS->getNumberHistograms() <= static_cast<size_t>(monIndex)) {
327 issues["MonitorWorkspaceIndex"] = "The MonitorWorkspace must contain the MonitorWorkspaceIndex.";
328 }
329 MatrixWorkspace_const_sptr inWS = getProperty("InputWorkspace");
330 if (monWS->getInstrument()->getName() != inWS->getInstrument()->getName()) {
331 issues["MonitorWorkspace"] = "The Input and Monitor workspaces must come "
332 "from the same instrument.";
333 }
334 if (monWS->getAxis(0)->unit()->unitID() != inWS->getAxis(0)->unit()->unitID()) {
335 issues["MonitorWorkspace"] = "The Input and Monitor workspaces must have the same unit";
336 }
337 }
338
339 return issues;
340}
341
346
347 // Check where the monitor spectrum should come from
348 Property *monSpec = getProperty("MonitorSpectrum");
349 MatrixWorkspace_sptr monWS = getProperty("MonitorWorkspace");
350 Property *monID = getProperty("MonitorID");
351 // Is the monitor spectrum within the main input workspace
352 const bool inWS = !monSpec->isDefault();
353 m_scanInput = inputWorkspace->detectorInfo().isScanning();
354 // Or is it in a separate workspace
355 bool sepWS{monWS};
356 if (m_scanInput && sepWS)
357 throw std::runtime_error("Can not currently use a separate monitor "
358 "workspace with a detector scan input workspace.");
359 // or monitor ID
360 bool monIDs = !monID->isDefault();
361 // something has to be set
362 // One and only one of these properties should have been set
363 // input from separate workspace is overwritten by monitor spectrum
364 if (inWS && sepWS) {
365 g_log.information("Both input workspace MonitorSpectrum number and monitor "
366 "workspace are specified. Ignoring Monitor Workspace");
367 sepWS = false;
368 }
369 // input from detector ID is rejected in favor of monitor sp
370 if (inWS && monIDs) {
371 g_log.information("Both input workspace MonitorSpectrum number and "
372 "detector ID are specified. Ignoring Detector ID");
373 monIDs = false;
374 }
375 // separate ws takes over detectorID (this logic is duplicated within
376 // getInWSMonitorSpectrum)
377 if (sepWS && monIDs) {
378 g_log.information("Both input MonitorWorkspace and detector ID are "
379 "specified. Ignoring Detector ID");
380 }
381
382 // Do a check for common binning and store
383 m_commonBins = inputWorkspace->isCommonBins();
384
385 // Check the monitor spectrum or workspace and extract into new workspace
386 m_monitor = sepWS ? getMonitorWorkspace(inputWorkspace) : getInWSMonitorSpectrum(inputWorkspace);
387
388 // Check that the 'monitor' spectrum actually relates to a monitor - warn if
389 // not
390 try {
391 const auto &monitorSpecInfo = m_monitor->spectrumInfo();
392 for (const auto workspaceIndex : m_workspaceIndexes)
393 if (!monitorSpecInfo.isMonitor(workspaceIndex))
394 g_log.warning() << "The spectrum N: " << workspaceIndex << " in MonitorWorkspace does not refer to a monitor.\n"
395 << "Continuing with normalization regardless.";
397 g_log.warning("Unable to check if the spectrum provided relates to a "
398 "monitor - the instrument is not fully specified.\n "
399 "Continuing with normalization regardless.");
400 g_log.warning() << "Error was: " << e.what() << "\n";
401 if (m_scanInput)
402 throw std::runtime_error("Can not continue, spectrum can not be obtained "
403 "for monitor workspace, but the input workspace "
404 "has a detector scan.");
405 }
406}
407
416 // this is the index of the spectra within the workspace and we need to
417 // identify it either from DetID or from SpecID
418 // size_t spectra_num(-1);
419 // try monitor spectrum. If it is specified, it overrides everything
420 int monitorSpec = getProperty("MonitorSpectrum");
421 if (monitorSpec < 0) {
422 // Get hold of the monitor spectrum through detector ID
423 int monitorID = getProperty("MonitorID");
424 if (monitorID < 0) {
425 throw std::runtime_error("Both MonitorSpectrum and MonitorID can not be negative");
426 }
427 // set spectra of detector's ID of one selected monitor ID
428 std::vector<detid_t> detID(1, monitorID);
429 // got the index of correspondent spectra (should be only one).
430 auto indexList = inputWorkspace->getIndicesFromDetectorIDs(detID);
431 if (indexList.empty()) {
432 throw std::runtime_error("Can not find spectra, corresponding to the requested monitor ID");
433 }
434 if (indexList.size() > 1 && !m_scanInput) {
435 throw std::runtime_error("More then one spectrum corresponds to the "
436 "requested monitor ID. This is unexpected in a "
437 "non-scanning workspace.");
438 }
439 m_workspaceIndexes = indexList;
440 } else { // monitor spectrum is specified.
441 if (m_scanInput)
442 throw std::runtime_error("For a scanning input workspace the monitor ID "
443 "must be provided. Normalisation can not be "
444 "performed to a spectrum.");
445 const SpectraAxis *axis = dynamic_cast<const SpectraAxis *>(inputWorkspace->getAxis(1));
446 if (!axis) {
447 throw std::runtime_error("Cannot retrieve monitor spectrum - spectrum "
448 "numbers not attached to workspace");
449 }
450 auto specs = axis->getSpectraIndexMap();
451 if (!specs.count(monitorSpec)) {
452 throw std::runtime_error("Input workspace does not contain spectrum "
453 "number given for MonitorSpectrum");
454 }
455 m_workspaceIndexes = std::vector<size_t>(1, specs[monitorSpec]);
456 }
457 return inputWorkspace;
458}
459
465 MatrixWorkspace_sptr monitorWS = getProperty("MonitorWorkspace");
466 const int wsID = getProperty("MonitorWorkspaceIndex");
467 m_workspaceIndexes = std::vector<size_t>(1, wsID);
468 // In this case we need to test whether the bins in the monitor workspace
469 // match
470 m_commonBins = (m_commonBins && WorkspaceHelpers::matchingBins(*inputWorkspace, *monitorWS, true));
471 // Copy the monitor spectrum because it will get changed
472 return monitorWS;
473}
474
488bool NormaliseToMonitor::setIntegrationProps(const bool isSingleCountWorkspace) {
489 m_integrationMin = getProperty("IntegrationRangeMin");
490 m_integrationMax = getProperty("IntegrationRangeMax");
491
492 // Check if neither of these have been changed from their defaults
493 // (EMPTY_DBL())
494 if ((isEmpty(m_integrationMin) && isEmpty(m_integrationMax)) && !isSingleCountWorkspace) {
495 // Nothing has been set so the user doesn't want to use integration so let's
496 // move on
497 return false;
498 }
499 // Yes integration is going to be used...
500
501 // Now check the end X values are within the X value range of the workspace
502 if ((isEmpty(m_integrationMin) || m_integrationMin < m_monitor->x(0).front()) && !isSingleCountWorkspace) {
503 g_log.warning() << "Integration range minimum set to workspace min: " << m_integrationMin << '\n';
504 m_integrationMin = m_monitor->x(0).front();
505 }
506 if ((isEmpty(m_integrationMax) || m_integrationMax > m_monitor->x(0).back()) && !isSingleCountWorkspace) {
507 g_log.warning() << "Integration range maximum set to workspace max: " << m_integrationMax << '\n';
508 m_integrationMax = m_monitor->x(0).back();
509 }
510
511 // Return indicating that these properties should be used
512 return true;
513}
514
523 MatrixWorkspace_sptr &outputWorkspace,
524 const bool isSingleCountWorkspace) {
526
527 // If single counting no need to integrate, monitor already guaranteed to be a
528 // single count
529 if (!isSingleCountWorkspace) {
530 // Add up all the bins so it's just effectively a series of values with
531 // errors
532 auto integrate = createChildAlgorithm("Integration");
533 integrate->setProperty<MatrixWorkspace_sptr>("InputWorkspace", m_monitor);
534 integrate->setProperty("RangeLower", m_integrationMin);
535 integrate->setProperty("RangeUpper", m_integrationMax);
536 integrate->setProperty<bool>("IncludePartialBins", getProperty("IncludePartialBins"));
537 integrate->executeAsChildAlg();
538 m_monitor = integrate->getProperty("OutputWorkspace");
539 }
540
541 EventWorkspace_sptr inputEvent = std::dynamic_pointer_cast<EventWorkspace>(inputWorkspace);
542
543 if (inputEvent) {
544 // Run the divide algorithm explicitly to enable progress reporting
545 auto divide = createChildAlgorithm("Divide", 0.0, 1.0);
546 divide->setProperty<MatrixWorkspace_sptr>("LHSWorkspace", inputWorkspace);
547 divide->setProperty<MatrixWorkspace_sptr>("RHSWorkspace", m_monitor);
548 divide->setProperty<MatrixWorkspace_sptr>("OutputWorkspace", outputWorkspace);
549 divide->executeAsChildAlg();
550
551 // Get back the result
552 outputWorkspace = divide->getProperty("OutputWorkspace");
553 } else {
554 performHistogramDivision(inputWorkspace, outputWorkspace);
555 }
556}
557
568 MatrixWorkspace_sptr &outputWorkspace) {
569 if (outputWorkspace != inputWorkspace)
570 outputWorkspace = inputWorkspace->clone();
571
572 size_t monitorWorkspaceIndex = 0;
573
574 Progress prog(this, 0.0, 1.0, m_workspaceIndexes.size());
575 const auto &specInfo = inputWorkspace->spectrumInfo();
576 for (const auto workspaceIndex : m_workspaceIndexes) {
577 // Errors propagated according to
578 // http://docs.mantidproject.org/nightly/concepts/ErrorPropagation.html#error-propagation
579 // This is similar to that in MantidAlgorithms::Divide
580
581 prog.report("Performing normalisation");
582
583 size_t timeIndex = 0;
584 if (m_scanInput)
585 timeIndex = specInfo.spectrumDefinition(workspaceIndex)[0].second;
586
587 const auto newYFactor = 1.0 / m_monitor->histogram(monitorWorkspaceIndex).y()[0];
588 const auto divisorError = m_monitor->histogram(monitorWorkspaceIndex).e()[0];
589 const double yErrorFactor = pow(divisorError * newYFactor, 2);
590 monitorWorkspaceIndex++;
591
592 PARALLEL_FOR_IF(Kernel::threadSafe(*outputWorkspace))
593 for (int64_t i = 0; i < int64_t(outputWorkspace->getNumberHistograms()); ++i) {
595 const auto &specDef = specInfo.spectrumDefinition(i);
596
597 if (!spectrumDefinitionsMatchTimeIndex(specDef, timeIndex))
598 continue;
599
600 auto hist = outputWorkspace->histogram(i);
601 auto &yValues = hist.mutableY();
602 auto &eValues = hist.mutableE();
603
604 for (size_t j = 0; j < yValues.size(); ++j) {
605 eValues[j] = newYFactor * sqrt(eValues[j] * eValues[j] + yValues[j] * yValues[j] * yErrorFactor);
606 yValues[j] *= newYFactor;
607 }
608
609 outputWorkspace->setHistogram(i, hist);
611 }
613 }
614}
615
621 MatrixWorkspace_sptr &outputWorkspace) {
622 EventWorkspace_sptr inputEvent = std::dynamic_pointer_cast<EventWorkspace>(inputWorkspace);
623
624 // Only create output workspace if different to input one
625 if (outputWorkspace != inputWorkspace) {
626 if (inputEvent) {
627 outputWorkspace = inputWorkspace->clone();
628 } else
629 outputWorkspace = create<MatrixWorkspace>(*inputWorkspace);
630 }
631 auto outputEvent = std::dynamic_pointer_cast<EventWorkspace>(outputWorkspace);
632
633 const auto &inputSpecInfo = inputWorkspace->spectrumInfo();
634 const auto &monitorSpecInfo = m_monitor->spectrumInfo();
635
636 const auto specLength = inputWorkspace->blocksize();
637 for (auto &workspaceIndex : m_workspaceIndexes) {
638 // Get hold of the monitor spectrum
639 const auto &monX = m_monitor->binEdges(workspaceIndex);
640 auto monY = m_monitor->counts(workspaceIndex);
641 auto monE = m_monitor->countStandardDeviations(workspaceIndex);
642 size_t timeIndex = 0;
643 if (m_scanInput)
644 timeIndex = monitorSpecInfo.spectrumDefinition(workspaceIndex)[0].second;
645 // Calculate the overall normalization just the once if bins are all
646 // matching
647 if (m_commonBins)
648 this->normalisationFactor(monX, monY, monE);
649
650 const size_t numHists = inputWorkspace->getNumberHistograms();
651 // Flag set when a division by 0 is found
652 bool hasZeroDivision = false;
653 Progress prog(this, 0.0, 1.0, numHists);
654 // Loop over spectra
655 PARALLEL_FOR_IF(Kernel::threadSafe(*inputWorkspace, *outputWorkspace, *m_monitor))
656 for (int64_t i = 0; i < int64_t(numHists); ++i) {
658 prog.report();
659
660 const auto &specDef = inputSpecInfo.spectrumDefinition(i);
661 if (!spectrumDefinitionsMatchTimeIndex(specDef, timeIndex))
662 continue;
663
664 const auto &X = inputWorkspace->binEdges(i);
665 // If not rebinning, just point to our monitor spectra, otherwise create
666 // new vectors
667
668 auto Y = (m_commonBins ? monY : Counts(specLength));
669 auto E = (m_commonBins ? monE : CountStandardDeviations(specLength));
670
671 if (!m_commonBins) {
672 // ConvertUnits can give X vectors of all zeros - skip these, they
673 // cause
674 // problems
675 if (X.back() == 0.0 && X.front() == 0.0)
676 continue;
677 // Rebin the monitor spectrum to match the binning of the current data
678 // spectrum
679 VectorHelper::rebinHistogram(monX.rawData(), monY.mutableRawData(), monE.mutableRawData(), X.rawData(),
680 Y.mutableRawData(), E.mutableRawData(), false);
681 // Recalculate the overall normalization factor
682 this->normalisationFactor(X, Y, E);
683 }
684
685 if (inputEvent) {
686 // --- EventWorkspace ---
687 EventList &outEL = outputEvent->getSpectrum(i);
688 outEL.divide(X.rawData(), Y.mutableRawData(), E.mutableRawData());
689 } else {
690 // --- Workspace2D ---
691 auto &YOut = outputWorkspace->mutableY(i);
692 auto &EOut = outputWorkspace->mutableE(i);
693 const auto &inY = inputWorkspace->y(i);
694 const auto &inE = inputWorkspace->e(i);
695 outputWorkspace->setSharedX(i, inputWorkspace->sharedX(i));
696
697 // The code below comes more or less straight out of Divide.cpp
698 for (size_t k = 0; k < specLength; ++k) {
699 // Get the input Y's
700 const double leftY = inY[k];
701 const double rightY = Y[k];
702
703 if (rightY == 0.0) {
704 hasZeroDivision = true;
705 }
706
707 // Calculate result and store in local variable to avoid overwriting
708 // original data if output workspace is same as one of the input
709 // ones
710 const double newY = leftY / rightY;
711
712 if (fabs(rightY) > 1.0e-12 && fabs(newY) > 1.0e-12) {
713 const double lhsFactor = (inE[k] < 1.0e-12 || fabs(leftY) < 1.0e-12) ? 0.0 : pow((inE[k] / leftY), 2);
714 const double rhsFactor = E[k] < 1.0e-12 ? 0.0 : pow((E[k] / rightY), 2);
715 EOut[k] = std::abs(newY) * sqrt(lhsFactor + rhsFactor);
716 }
717
718 // Now store the result
719 YOut[k] = newY;
720 } // end Workspace2D case
721 } // end loop over current spectrum
722
724 } // end loop over spectra
726
727 if (hasZeroDivision) {
728 g_log.warning() << "Division by zero in some of the bins.\n";
729 }
730 if (inputEvent)
731 outputEvent->clearMRU();
732 }
733}
734
743void NormaliseToMonitor::normalisationFactor(const BinEdges &X, Counts &Y, CountStandardDeviations &E) {
744 const double monitorSum = std::accumulate(Y.begin(), Y.end(), 0.0);
745 const double range = X.back() - X.front();
746 auto specLength = Y.size();
747
748 auto &yNew = Y.mutableRawData();
749 auto &eNew = E.mutableRawData();
750
751 for (size_t j = 0; j < specLength; ++j) {
752 const double factor = range / ((X[j + 1] - X[j]) * monitorSum);
753 yNew[j] *= factor;
754 eNew[j] *= factor;
755 }
756}
757
758} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
#define fabs(x)
Definition: Matrix.cpp:22
#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_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.
#define UNUSED_ARG(x)
Function arguments are sometimes unused in certain implmentations but are required for documentation ...
Definition: System.h:64
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
Kernel::Property * getPointerToProperty(const std::string &name) const override
Get a property by name.
Definition: Algorithm.cpp:2033
std::string getPropertyValue(const std::string &name) const override
Get the value of a property as a string.
Definition: Algorithm.cpp:2026
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Definition: Algorithm.cpp:2076
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.
Definition: Algorithm.cpp:842
Kernel::Logger & g_log
Definition: Algorithm.h:451
static bool isEmpty(const NumT toCheck)
checks that the value was not set by users, uses the value in empty double/int.
A validator which checks that a workspace contains histogram data (the default) or point data as requ...
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
A validator which checks that a workspace contains raw counts in its bins.
SingleCountValidator : This validator checks that there is only a single entry per spectrum,...
Class to represent the spectra axis of a workspace.
Definition: SpectraAxis.h:31
spec2index_map getSpectraIndexMap() const
Returns a map where spectra is the key and index is the value This is used for efficient search of sp...
A property class for workspaces.
void applyChanges(const Mantid::Kernel::IPropertyManager *algo, Kernel::Property *const pProp) override
The function user have to overload it in their custom code to modify the property according to the ch...
bool monitorIdReader(const API::MatrixWorkspace_const_sptr &inputWS) const
bool isConditionChanged(const Mantid::Kernel::IPropertyManager *algo, const std::string &changedPropName="") const override
to verify if the properties, this one depends on have changed or other special condition occurs which...
bool isEnabled(const Mantid::Kernel::IPropertyManager *algo) const override
Is the property to be shown as "enabled" in the GUI.
void performHistogramDivision(const API::MatrixWorkspace_sptr &inputWorkspace, API::MatrixWorkspace_sptr &outputWorkspace)
This performs a similar operation to divide, but is a separate algorithm so that the correct spectra ...
bool m_commonBins
Whether the input workspace has common bins.
double m_integrationMax
The upper bound of the integration range.
void normalisationFactor(const HistogramData::BinEdges &X, HistogramData::Counts &Y, HistogramData::CountStandardDeviations &E)
Calculates the overall normalization factor.
void normaliseBinByBin(const API::MatrixWorkspace_sptr &inputWorkspace, API::MatrixWorkspace_sptr &outputWorkspace)
Carries out the bin-by-bin normalization.
void exec() override
Virtual method - must be overridden by concrete algorithm.
void normaliseByIntegratedCount(const API::MatrixWorkspace_sptr &inputWorkspace, API::MatrixWorkspace_sptr &outputWorkspace, const bool isSingleCountWorkspace)
Carries out a normalization based on the integrated count of the monitor over a range.
void init() override
Virtual method - must be overridden by concrete algorithm.
API::MatrixWorkspace_sptr getInWSMonitorSpectrum(const API::MatrixWorkspace_sptr &inputWorkspace)
Checks and retrieves the requested spectrum out of the input workspace.
double m_integrationMin
The lower bound of the integration range.
API::MatrixWorkspace_sptr m_monitor
A single spectrum workspace containing the monitor.
std::map< std::string, std::string > validateInputs() override
Validates input properties.
API::MatrixWorkspace_sptr getMonitorWorkspace(const API::MatrixWorkspace_sptr &inputWorkspace)
Checks and retrieves the monitor spectrum out of the input workspace.
void checkProperties(const API::MatrixWorkspace_sptr &inputWorkspace)
Makes sure that the input properties are set correctly.
bool setIntegrationProps(const bool isSingleCountWorkspace)
Sets the maximum and minimum X values of the monitor spectrum to use for integration.
API::MatrixWorkspace_sptr extractMonitorSpectra(const API::MatrixWorkspace_sptr &ws, const std::vector< size_t > &workspaceIndexes)
Pulls the monitor spectra out of a larger workspace.
A class for holding :
Definition: EventList.h:56
void divide(const double value, const double error=0.0) override
Divide the weights in this event list by a scalar with an (optional) error.
Definition: EventList.cpp:3633
BoundedValidator is a validator that requires the values to be between upper or lower bounds,...
Exception for when an item is not found in a collection.
Definition: Exception.h:145
const char * what() const noexcept override
Writes out the range and limits.
Definition: Exception.cpp:105
Interface to PropertyManager.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void setPropertySettings(const std::string &name, std::unique_ptr< IPropertySettings > settings)
virtual TypedValue getProperty(const std::string &name) const =0
Get the value of a property.
ListValidator is a validator that requires the value of a property to be one of a defined list of pos...
Definition: ListValidator.h:29
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
Definition: ProgressBase.h:51
The concrete, templated class for properties.
Base class for properties.
Definition: Property.h:94
virtual bool isDefault() const =0
Overriden function that returns if property has the same value that it was initialised with,...
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
bool spectrumDefinitionsMatchTimeIndex(const SpectrumDefinition &specDef, const size_t timeIndex)
std::shared_ptr< EventWorkspace > EventWorkspace_sptr
shared pointer to the EventWorkspace class
std::shared_ptr< const Instrument > Instrument_const_sptr
Shared pointer to an const instrument object.
void MANTID_KERNEL_DLL rebinHistogram(const std::vector< double > &xold, const std::vector< double > &yold, const std::vector< double > &eold, const std::vector< double > &xnew, std::vector< double > &ynew, std::vector< double > &enew, bool addition)
Rebins histogram data according to a new output X array.
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
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
Definition: EmptyValues.h:43
static bool matchingBins(const MatrixWorkspace &ws1, const MatrixWorkspace &ws2, const bool firstOnly=false)
Checks whether the bins (X values) of two workspace are the same.
@ InOut
Both an input & output workspace.
Definition: Property.h:55
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54