Mantid
Loading...
Searching...
No Matches
EQSANSDarkCurrentSubtraction2.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
4// NScD Oak Ridge National Laboratory, European Spallation Source,
5// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6// SPDX - License - Identifier: GPL - 3.0 +
11#include "MantidAPI/Run.h"
19
20#include "Poco/Path.h"
21#include "Poco/String.h"
22
24
25// Register the algorithm into the AlgorithmFactory
26DECLARE_ALGORITHM(EQSANSDarkCurrentSubtraction2)
27
28using namespace Kernel;
29using namespace API;
30using namespace Geometry;
31using namespace DataObjects;
32
34 auto wsValidator = std::make_shared<WorkspaceUnitValidator>("Wavelength");
35 declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input, wsValidator));
36
37 declareProperty(std::make_unique<API::FileProperty>("Filename", "", API::FileProperty::Load, "_event.nxs"),
38 "The name of the input event Nexus file to load as dark current.");
39
40 declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output));
41 declareProperty("PersistentCorrection", true,
42 "If true, the algorithm will be persistent and re-used when "
43 "other data sets are processed");
44 declareProperty("ReductionProperties", "__sans_reduction_properties", Direction::Input);
45 declareProperty(std::make_unique<WorkspaceProperty<MatrixWorkspace>>("OutputDarkCurrentWorkspace", "",
47 declareProperty("OutputMessage", "", Direction::Output);
48}
49
51 std::string output_message;
52 // Reduction property manager
53 const std::string reductionManagerName = getProperty("ReductionProperties");
54 std::shared_ptr<PropertyManager> reductionManager;
55 if (PropertyManagerDataService::Instance().doesExist(reductionManagerName)) {
56 reductionManager = PropertyManagerDataService::Instance().retrieve(reductionManagerName);
57 } else {
58 reductionManager = std::make_shared<PropertyManager>();
59 PropertyManagerDataService::Instance().addOrReplace(reductionManagerName, reductionManager);
60 }
61
62 // If the load algorithm isn't in the reduction properties, add it
63 const bool persistent = getProperty("PersistentCorrection");
64 if (!reductionManager->existsProperty("DarkCurrentAlgorithm") && persistent) {
65 auto algProp = std::make_unique<AlgorithmProperty>("DarkCurrentAlgorithm");
66 algProp->setValue(toString());
67 reductionManager->declareProperty(std::move(algProp));
68 }
69
70 Progress progress(this, 0.0, 1.0, 10);
71
72 MatrixWorkspace_sptr inputWS = getProperty("InputWorkspace");
73
74 // This version of dark current subtraction only works on histograms.
75 // Users need to either make sure the EQSANSLoad algorithm produces
76 // histograms, or turn off the dark current subtraction.
77 EventWorkspace_const_sptr inputEventWS = std::dynamic_pointer_cast<const EventWorkspace>(inputWS);
78 if (inputEventWS) {
79 g_log.error() << "To use this version of EQSANSDarkCurrentSubtraction, "
80 << "you need to make sure EQSANSLoad produces histograms. "
81 << "You can also turn the dark current subtraction off.\n";
82 throw std::invalid_argument("EQSANSDarkCurrentSubtraction-v2 only works on histograms.");
83 }
84
85 const std::string fileName = getPropertyValue("Filename");
87
88 progress.report("Subtracting dark current");
89
90 // Look for an entry for the dark current in the reduction table
91 Poco::Path path(fileName);
92 const std::string entryName = "DarkCurrent" + path.getBaseName();
93 std::string darkWSName = "__dark_current_" + path.getBaseName();
94
95 if (reductionManager->existsProperty(entryName)) {
96 darkWS = reductionManager->getProperty(entryName);
97 darkWSName = reductionManager->getPropertyValue(entryName);
98 output_message += darkWSName + '\n';
99 } else {
100 // Load the dark current if we don't have it already
101 IAlgorithm_sptr loadAlg;
102 if (!reductionManager->existsProperty("LoadAlgorithm")) {
103 loadAlg = createChildAlgorithm("EQSANSLoad", 0.1, 0.3);
104 loadAlg->setProperty("Filename", fileName);
105 if (loadAlg->existsProperty("LoadMonitors"))
106 loadAlg->setProperty("LoadMonitors", false);
107 loadAlg->executeAsChildAlg();
108 darkWS = loadAlg->getProperty("OutputWorkspace");
109 } else {
110 // Get load algorithm as a string so that we can create a completely
111 // new proxy and ensure that we don't overwrite existing properties
112 IAlgorithm_sptr loadAlg0 = reductionManager->getProperty("LoadAlgorithm");
113 const std::string loadString = loadAlg0->toString();
114 loadAlg = Algorithm::fromString(loadString);
115 loadAlg->setChild(true);
116 loadAlg->setProperty("Filename", fileName);
117 if (loadAlg->existsProperty("LoadMonitors"))
118 loadAlg->setProperty("LoadMonitors", false);
119 loadAlg->setPropertyValue("OutputWorkspace", darkWSName);
120 loadAlg->execute();
121 darkWS = loadAlg->getProperty("OutputWorkspace");
122 }
123
124 output_message += "\n Loaded " + fileName + "\n";
125 if (loadAlg->existsProperty("OutputMessage")) {
126 std::string msg = loadAlg->getPropertyValue("OutputMessage");
127 output_message += " |" + Poco::replace(msg, "\n", "\n |") + "\n";
128 }
129
130 std::string darkWSOutputName = getPropertyValue("OutputDarkCurrentWorkspace");
131 if (!darkWSOutputName.empty())
132 setProperty("OutputDarkCurrentWorkspace", darkWS);
133 AnalysisDataService::Instance().addOrReplace(darkWSName, darkWS);
134 reductionManager->declareProperty(std::make_unique<WorkspaceProperty<>>(entryName, "", Direction::Output));
135 reductionManager->setPropertyValue(entryName, darkWSName);
136 reductionManager->setProperty(entryName, darkWS);
137 }
138 progress.report(3, "Loaded dark current");
139
140 // Normalize the dark current and data to counting time
141 double scaling_factor = 1.0;
142 if (inputWS->run().hasProperty("duration")) {
143 auto duration = inputWS->run().getPropertyValueAsType<double>("duration");
144 auto dark_duration = darkWS->run().getPropertyValueAsType<double>("duration");
145 ;
146 scaling_factor = duration / dark_duration;
147 } else if (inputWS->run().hasProperty("proton_charge")) {
148 auto dp = inputWS->run().getTimeSeriesProperty<double>("proton_charge");
149 double duration = dp->getStatistics().duration;
150
151 dp = darkWS->run().getTimeSeriesProperty<double>("proton_charge");
152 double dark_duration = dp->getStatistics().duration;
153 scaling_factor = duration / dark_duration;
154 } else if (inputWS->run().hasProperty("timer")) {
155 auto duration = inputWS->run().getPropertyValueAsType<double>("timer");
156 auto dark_duration = darkWS->run().getPropertyValueAsType<double>("timer");
157 ;
158 scaling_factor = duration / dark_duration;
159 } else {
160 output_message += "\n Could not find proton charge or duration in sample logs";
161 g_log.error() << "ERROR: Could not find proton charge or duration in sample logs\n";
162 };
163 // The scaling factor should account for the TOF cuts on each side of a frame
164 // The EQSANSLoad algorithm cuts the beginning and end of the TOF distribution
165 // so we don't need to correct the scaling factor here. When using
166 // LoadEventNexus
167 // we have to scale by (t_frame-t_low_cut-t_high_cut)/t_frame.
168
169 progress.report("Scaling dark current");
170
171 // Get the dark current counts per pixel
172 auto rebinAlg = createChildAlgorithm("Integration", 0.4, 0.5);
173 rebinAlg->setProperty("InputWorkspace", darkWS);
174 rebinAlg->setProperty("OutputWorkspace", darkWS);
175 rebinAlg->executeAsChildAlg();
176 MatrixWorkspace_sptr scaledDarkWS = rebinAlg->getProperty("OutputWorkspace");
177
178 // Scale the dark current
179 auto scaleAlg = createChildAlgorithm("Scale", 0.5, 0.6);
180 scaleAlg->setProperty("InputWorkspace", scaledDarkWS);
181 scaleAlg->setProperty("Factor", scaling_factor);
182 scaleAlg->setProperty("OutputWorkspace", scaledDarkWS);
183 scaleAlg->setProperty("Operation", "Multiply");
184 scaleAlg->executeAsChildAlg();
185 scaledDarkWS = rebinAlg->getProperty("OutputWorkspace");
186
187 // Scale the dark counts to the bin width and perform subtraction
188 const auto numberOfSpectra = static_cast<int>(inputWS->getNumberHistograms());
189 const auto numberOfDarkSpectra = static_cast<int>(scaledDarkWS->getNumberHistograms());
190 if (numberOfSpectra != numberOfDarkSpectra) {
191 g_log.error() << "Incompatible number of pixels between sample run and "
192 "dark current\n";
193 }
194 const auto nBins = static_cast<int>(inputWS->readY(0).size());
195 const auto xLength = static_cast<int>(inputWS->readX(0).size());
196 if (xLength != nBins + 1) {
197 g_log.error() << "The input workspaces are expected to be histograms\n";
198 }
199
200 progress.report("Subtracting dark current");
201 const auto &spectrumInfo = inputWS->spectrumInfo();
202 // Loop over all tubes and patch as necessary
203 for (int i = 0; i < numberOfSpectra; i++) {
204 // If this detector is a monitor, skip to the next one
205 if (spectrumInfo.isMasked(i))
206 continue;
207
208 const MantidVec &YDarkValues = scaledDarkWS->readY(i);
209 const MantidVec &YDarkErrors = scaledDarkWS->readE(i);
210 const MantidVec &XValues = inputWS->readX(i);
211 MantidVec &YValues = inputWS->dataY(i);
212 MantidVec &YErrors = inputWS->dataE(i);
213 for (int j = 0; j < nBins; j++) {
214 double bin_scale = (XValues[j + 1] - XValues[j]) / (XValues[nBins] - XValues[0]);
215 YValues[j] -= YDarkValues[0] * bin_scale;
216 YErrors[j] = sqrt(YErrors[j] * YErrors[j] + YDarkErrors[0] * YDarkErrors[0] * bin_scale * bin_scale);
217 }
218 }
219 setProperty("OutputWorkspace", inputWS);
220 setProperty("OutputMessage", "Dark current subtracted: " + output_message);
221
222 progress.report("Subtracted dark current");
223}
224
225} // namespace Mantid::WorkflowAlgorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
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
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
std::string toString() const override
Serialize an object to a string.
Definition: Algorithm.cpp:905
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
static IAlgorithm_sptr fromString(const std::string &input)
De-serialize an object from a string.
Definition: Algorithm.cpp:967
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
@ Load
allowed here which will be passed to the algorithm
Definition: FileProperty.h:52
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
A property class for workspaces.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void error(const std::string &msg)
Logs at error level.
Definition: Logger.cpp:77
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
std::shared_ptr< IAlgorithm > IAlgorithm_sptr
shared pointer to Mantid::API::IAlgorithm
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
std::shared_ptr< const EventWorkspace > EventWorkspace_const_sptr
shared pointer to a const Workspace2D
std::vector< double > MantidVec
typedef for the data storage used in Mantid matrix workspaces
Definition: cow_ptr.h:172
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54