20#include "MantidHistogramData/Histogram.h"
26#include "MantidTypes/SpectrumDefinition.h"
78 throw(std::invalid_argument(
"modify allowed value has been called on wrong property"));
87 spectra_max =
static_cast<int>(inputWS->getNumberHistograms()) + 1;
108 std::vector<detid_t> monitorIDList = pInstr->getMonitors();
110 const auto &specInfo = inputWS->spectrumInfo();
111 std::set<detid_t> idsInWorkspace;
114 while (i < specInfo.size() && idsInWorkspace.size() < monitorIDList.size()) {
115 if (specInfo.isMonitor(i))
116 idsInWorkspace.insert(specInfo.detector(i).getID());
119 monitorIDList = std::vector<detid_t>(idsInWorkspace.begin(), idsInWorkspace.end());
122 if (monitorIDList.empty()) {
139 bool values_redefined =
false;
140 for (
size_t i = 0; i < monitorIDList.size(); i++) {
142 values_redefined =
true;
146 return values_redefined;
150 return std::none_of(specDef.cbegin(), specDef.cend(),
151 [timeIndex](
const auto &spec) { return spec.second != timeIndex; });
157using namespace Kernel;
164 auto validatorHistSingle = std::make_shared<CompositeValidator>(CompositeRelation::OR);
167 auto validator = std::make_shared<CompositeValidator>();
168 validator->add(validatorHistSingle);
172 "Name of the input workspace. Must be a non-distribution histogram.");
175 "Name to use for the output workspace");
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)",
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.");
197 std::make_unique<MonIDPropChanger>(
"InputWorkspace",
"MonitorSpectrum",
"MonitorWorkspace"));
202 "A workspace containing one or more spectra to normalize the "
203 "InputWorkspace by.");
207 "The index of the spectrum within the MonitorWorkspace(2 "
208 "(0<=ind<=nHistograms in MonitorWorkspace) you want to "
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)",
216 std::make_unique<Kernel::EnabledWhenProperty>(
"MonitorSpectrum",
IS_DEFAULT));
223 "If set, normalization will be by integrated count from this "
226 "If set, normalization will be by integrated count up to "
227 "this maximum x value");
229 "If true and an integration range is set then partial bins at either \n"
230 "end of the integration range are also included");
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");
246 bool isSingleCountWorkspace =
false;
248 isSingleCountWorkspace = (!inputWS->isHistogramData()) && (inputWS->blocksize() == 1);
249 }
catch (std::length_error &) {
251 isSingleCountWorkspace =
false;
264 if (!norm_ws_name.empty()) {
266 if (out_name == norm_ws_name) {
269 norm_ws_name = norm_ws_name +
"_normFactor";
271 pProp->setValue(norm_ws_name);
286 const std::vector<std::size_t> &workspaceIndexes) {
289 childAlg->setProperty(
"WorkspaceIndexList", workspaceIndexes);
290 childAlg->executeAsChildAlg();
299 std::map<std::string, std::string> issues;
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;
313 const double intMin =
getProperty(
"IntegrationRangeMin");
314 const double intMax =
getProperty(
"IntegrationRangeMax");
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.";
323 const int monIndex =
getProperty(
"MonitorWorkspaceIndex");
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.";
330 if (monWS->getInstrument()->getName() != inWS->getInstrument()->getName()) {
331 issues[
"MonitorWorkspace"] =
"The Input and Monitor workspaces must come "
332 "from the same instrument.";
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";
353 m_scanInput = inputWorkspace->detectorInfo().isScanning();
357 throw std::runtime_error(
"Can not currently use a separate monitor "
358 "workspace with a detector scan input workspace.");
365 g_log.
information(
"Both input workspace MonitorSpectrum number and monitor "
366 "workspace are specified. Ignoring Monitor Workspace");
370 if (inWS && monIDs) {
372 "detector ID are specified. Ignoring Detector ID");
377 if (sepWS && monIDs) {
379 "specified. Ignoring Detector ID");
391 const auto &monitorSpecInfo =
m_monitor->spectrumInfo();
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.");
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.");
421 if (monitorSpec < 0) {
425 throw std::runtime_error(
"Both MonitorSpectrum and MonitorID can not be negative");
428 std::vector<detid_t> detID(1, monitorID);
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");
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.");
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.");
447 throw std::runtime_error(
"Cannot retrieve monitor spectrum - spectrum "
448 "numbers not attached to workspace");
451 if (!specs.count(monitorSpec)) {
452 throw std::runtime_error(
"Input workspace does not contain spectrum "
453 "number given for MonitorSpectrum");
457 return inputWorkspace;
466 const int wsID =
getProperty(
"MonitorWorkspaceIndex");
524 const bool isSingleCountWorkspace) {
529 if (!isSingleCountWorkspace) {
536 integrate->setProperty<
bool>(
"IncludePartialBins",
getProperty(
"IncludePartialBins"));
537 integrate->executeAsChildAlg();
538 m_monitor = integrate->getProperty(
"OutputWorkspace");
541 EventWorkspace_sptr inputEvent = std::dynamic_pointer_cast<EventWorkspace>(inputWorkspace);
549 divide->executeAsChildAlg();
552 outputWorkspace = divide->getProperty(
"OutputWorkspace");
569 if (outputWorkspace != inputWorkspace)
570 outputWorkspace = inputWorkspace->clone();
572 size_t monitorWorkspaceIndex = 0;
575 const auto &specInfo = inputWorkspace->spectrumInfo();
581 prog.
report(
"Performing normalisation");
583 size_t timeIndex = 0;
585 timeIndex = specInfo.spectrumDefinition(workspaceIndex)[0].second;
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++;
593 for (int64_t i = 0; i < int64_t(outputWorkspace->getNumberHistograms()); ++i) {
595 const auto &specDef = specInfo.spectrumDefinition(i);
600 auto hist = outputWorkspace->histogram(i);
601 auto &yValues = hist.mutableY();
602 auto &eValues = hist.mutableE();
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;
609 outputWorkspace->setHistogram(i, hist);
622 EventWorkspace_sptr inputEvent = std::dynamic_pointer_cast<EventWorkspace>(inputWorkspace);
625 if (outputWorkspace != inputWorkspace) {
627 outputWorkspace = inputWorkspace->clone();
629 outputWorkspace = create<MatrixWorkspace>(*inputWorkspace);
631 auto outputEvent = std::dynamic_pointer_cast<EventWorkspace>(outputWorkspace);
633 const auto &inputSpecInfo = inputWorkspace->spectrumInfo();
634 const auto &monitorSpecInfo =
m_monitor->spectrumInfo();
636 const auto specLength = inputWorkspace->blocksize();
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;
644 timeIndex = monitorSpecInfo.spectrumDefinition(workspaceIndex)[0].second;
650 const size_t numHists = inputWorkspace->getNumberHistograms();
652 bool hasZeroDivision =
false;
653 Progress prog(
this, 0.0, 1.0, numHists);
656 for (int64_t i = 0; i < int64_t(numHists); ++i) {
660 const auto &specDef = inputSpecInfo.spectrumDefinition(i);
664 const auto &
X = inputWorkspace->binEdges(i);
669 auto E = (
m_commonBins ? monE : CountStandardDeviations(specLength));
675 if (
X.back() == 0.0 &&
X.front() == 0.0)
680 Y.mutableRawData(), E.mutableRawData(),
false);
687 EventList &outEL = outputEvent->getSpectrum(i);
688 outEL.
divide(
X.rawData(),
Y.mutableRawData(), E.mutableRawData());
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));
698 for (
size_t k = 0; k < specLength; ++k) {
700 const double leftY = inY[k];
701 const double rightY =
Y[k];
704 hasZeroDivision =
true;
710 const double newY = leftY / rightY;
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);
727 if (hasZeroDivision) {
728 g_log.
warning() <<
"Division by zero in some of the bins.\n";
731 outputEvent->clearMRU();
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();
748 auto &yNew =
Y.mutableRawData();
749 auto &eNew = E.mutableRawData();
751 for (
size_t j = 0; j < specLength; ++j) {
752 const double factor = range / ((
X[j + 1] -
X[j]) * monitorSum);
#define DECLARE_ALGORITHM(classname)
#define PARALLEL_START_INTERRUPT_REGION
Begins a block to skip processing is the algorithm has been interupted Note the end of the block if n...
#define PARALLEL_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 ...
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
Kernel::Property * getPointerToProperty(const std::string &name) const override
Get a property by name.
std::string getPropertyValue(const std::string &name) const override
Get the value of a property as a string.
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
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.
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.
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.
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.
std::string MonitorWorkspaceProp
std::vector< int > iExistingAllowedValues
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.
std::vector< size_t > m_workspaceIndexes
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.
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.
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.
const char * what() const noexcept override
Writes out the range and limits.
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...
void warning(const std::string &msg)
Logs at warning level.
void information(const std::string &msg)
Logs at information level.
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
The concrete, templated class for properties.
Base class for properties.
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.
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
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.
@ Input
An input workspace.
@ Output
An output workspace.