Mantid
Loading...
Searching...
No Matches
DeadTimeCorrection.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 +
8#include "MantidAPI/Axis.h"
11#include "MantidAPI/Progress.h"
14
15#include <limits>
16
17namespace Mantid::Algorithms {
18
19using API::FileProperty;
21using API::Progress;
22using API::WorkspaceProperty;
23using Kernel::BoundedValidator;
24using Kernel::Direction;
25using Kernel::PropertyWithValue;
26
27// Register the algorithm into the AlgorithmFactory
28DECLARE_ALGORITHM(DeadTimeCorrection)
29
30//----------------------------------------------------------------------------------------------
31
32
33const std::string DeadTimeCorrection::name() const { return "DeadTimeCorrection"; }
34
36int DeadTimeCorrection::version() const { return 1; }
37
39const std::string DeadTimeCorrection::category() const { return "CorrectionFunctions"; }
40
42const std::string DeadTimeCorrection::summary() const { return "Performs a dead time correction based on count rate."; }
43
44//----------------------------------------------------------------------------------------------
48 declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input), "An input workspace.");
49
50 declareProperty(std::make_unique<PropertyWithValue<std::string>>("GroupingPattern", "", Direction::Input),
51 "See the GroupingPattern documentation of GroupDetectors.");
52
53 auto positive = std::make_shared<BoundedValidator<double>>();
54 positive->setLower(0.);
55 declareProperty("Tau", 0., positive, "The count rate coefficient.");
56
57 declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
58 "An output workspace.");
59
60 declareProperty(std::make_unique<FileProperty>("MapFile", "", FileProperty::OptionalLoad,
61 std::vector<std::string>{".map", ".xml"}),
62 "A file that consists of lists of spectra numbers to group. See the "
63 "help of GroupDetectors for the file format");
64}
65
66//----------------------------------------------------------------------------------------------
70 MatrixWorkspace_sptr inputWorkspace = getProperty("InputWorkspace");
71 MatrixWorkspace_sptr outputWorkspace = getProperty("OutputWorkspace");
72 if (inputWorkspace != outputWorkspace) {
73 outputWorkspace = inputWorkspace->clone();
74 }
75 const auto map = outputWorkspace->getDetectorIDToWorkspaceIndexMap();
76 const double tau = getProperty("Tau");
77 MatrixWorkspace_sptr integrated = inputWorkspace;
78 auto xAxisUnitID = inputWorkspace->getAxis(0)->unit()->unitID();
79 if (outputWorkspace->blocksize() != 1 && xAxisUnitID != "Empty") {
80 auto integrator = createChildAlgorithm("Integration");
81 integrator->setProperty("InputWorkspace", inputWorkspace);
82 integrator->setPropertyValue("OutputWorkspace", "unused");
83 integrator->executeAsChildAlg();
84 integrated = integrator->getProperty("OutputWorkspace");
85 // after integration we end up with one bin
86 // however the bin edges might vary, which does not matter, we just need to
87 // group the counts, hence we need to do this before we can group the pixels
88 for (size_t index = 1; index < integrated->getNumberHistograms(); ++index) {
89 integrated->setSharedX(index, integrated->sharedX(0));
90 }
91 }
92 const std::string groupingPattern = getProperty("GroupingPattern");
93 MatrixWorkspace_sptr grouped = integrated;
94 if (!groupingPattern.empty()) {
95 auto grouper = createChildAlgorithm("GroupDetectors");
96 grouper->setProperty("InputWorkspace", integrated);
97 grouper->setPropertyValue("OutputWorkspace", "unused");
98 grouper->setPropertyValue("GroupingPattern", groupingPattern);
99 grouper->setPropertyValue("Behaviour", "Sum");
100 grouper->setProperty("KeepUngroupedSpectra", true);
101 grouper->executeAsChildAlg();
102 grouped = grouper->getProperty("OutputWorkspace");
103 } else {
104 const std::string filename = getProperty("MapFile");
105 if (!filename.empty()) {
106 auto grouper = createChildAlgorithm("GroupDetectors");
107 grouper->setProperty("InputWorkspace", integrated);
108 grouper->setPropertyValue("OutputWorkspace", "unused");
109 grouper->setPropertyValue("MapFile", filename);
110 grouper->setPropertyValue("Behaviour", "Sum");
111 grouper->setProperty("KeepUngroupedSpectra", true);
112 grouper->executeAsChildAlg();
113 grouped = grouper->getProperty("OutputWorkspace");
114 }
115 }
116 Progress progress(this, 0.0, 1.0, grouped->getNumberHistograms());
117 PARALLEL_FOR_IF(Kernel::threadSafe(*outputWorkspace))
118 for (int index = 0; index < static_cast<int>(grouped->getNumberHistograms()); ++index) {
120 progress.report("Performing the correction for the group at index " + std::to_string(index));
121 auto correction = grouped->readY(index);
122 for (size_t bin = 0; bin < correction.size(); bin++) {
123 auto &y = correction[bin];
124 if (y >= 1.0 / tau) {
125 g_log.warning() << "Saturation count rate reached for grouped detector at index " << index << ", in bin " << bin
126 << ". Correction will be infinity. Check your tau or input "
127 "workspace, make sure it is normalised by acquisition time.\n";
128 y = std::numeric_limits<double>::infinity();
129 }
130 }
131 std::for_each(correction.begin(), correction.end(), [&tau](double &y) { y = 1.0 / (1.0 - y * tau); });
132
133 const auto detIDs = grouped->getSpectrum(index).getDetectorIDs();
134 for (const auto id : detIDs) {
135 const size_t originalIndex = map.at(id);
136 auto &spectrum = outputWorkspace->mutableY(originalIndex);
137 auto &errors = outputWorkspace->mutableE(originalIndex);
138 if (correction.size() == 1) {
139 spectrum *= correction[0];
140 errors *= correction[0];
141 } else {
142 spectrum *= correction;
143 errors *= correction;
144 }
145 }
147 }
149 setProperty("OutputWorkspace", outputWorkspace);
150}
151
152} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
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_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.
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
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
void progress(double p, const std::string &msg="", double estimatedTime=0.0, int progressPrecision=0)
Sends ProgressNotification.
Definition: Algorithm.cpp:231
@ OptionalLoad
to specify a file to read but the file doesn't have to exist
Definition: FileProperty.h:53
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
A property class for workspaces.
DeadTimeCorrection : Performs dead time correction.
void exec() override
Execute the algorithm.
int version() const override
Algorithm's version for identification.
const std::string summary() const override
Algorithm's summary for use in the GUI and help.
void init() override
Initialize the algorithm's properties.
const std::string category() const override
Algorithm's category for identification.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
The concrete, templated class for properties.
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
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
STL namespace.
std::string to_string(const wide_integer< Bits, Signed > &n)
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54