21#include "Poco/String.h"
28using namespace Kernel;
30using namespace Geometry;
31using namespace DataObjects;
34 auto wsValidator = std::make_shared<WorkspaceUnitValidator>(
"Wavelength");
38 "The name of the input event Nexus file to load as dark current.");
42 "If true, the algorithm will be persistent and re-used when "
43 "other data sets are processed");
51 std::string output_message;
53 const std::string reductionManagerName =
getProperty(
"ReductionProperties");
54 std::shared_ptr<PropertyManager> reductionManager;
58 reductionManager = std::make_shared<PropertyManager>();
63 const bool persistent =
getProperty(
"PersistentCorrection");
64 if (!reductionManager->existsProperty(
"DarkCurrentAlgorithm") && persistent) {
65 auto algProp = std::make_unique<AlgorithmProperty>(
"DarkCurrentAlgorithm");
67 reductionManager->declareProperty(std::move(algProp));
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.");
88 progress.report(
"Subtracting dark current");
91 Poco::Path path(fileName);
92 const std::string entryName =
"DarkCurrent" + path.getBaseName();
93 std::string darkWSName =
"__dark_current_" + path.getBaseName();
95 if (reductionManager->existsProperty(entryName)) {
96 darkWS = reductionManager->getProperty(entryName);
97 darkWSName = reductionManager->getPropertyValue(entryName);
98 output_message += darkWSName +
'\n';
102 if (!reductionManager->existsProperty(
"LoadAlgorithm")) {
104 loadAlg->setProperty(
"Filename", fileName);
105 if (loadAlg->existsProperty(
"LoadMonitors"))
106 loadAlg->setProperty(
"LoadMonitors",
false);
107 loadAlg->executeAsChildAlg();
108 darkWS = loadAlg->getProperty(
"OutputWorkspace");
112 IAlgorithm_sptr loadAlg0 = reductionManager->getProperty(
"LoadAlgorithm");
113 const std::string loadString = loadAlg0->toString();
115 loadAlg->setChild(
true);
116 loadAlg->setProperty(
"Filename", fileName);
117 if (loadAlg->existsProperty(
"LoadMonitors"))
118 loadAlg->setProperty(
"LoadMonitors",
false);
119 loadAlg->setPropertyValue(
"OutputWorkspace", darkWSName);
121 darkWS = loadAlg->getProperty(
"OutputWorkspace");
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";
130 std::string darkWSOutputName =
getPropertyValue(
"OutputDarkCurrentWorkspace");
131 if (!darkWSOutputName.empty())
135 reductionManager->setPropertyValue(entryName, darkWSName);
136 reductionManager->setProperty(entryName, darkWS);
138 progress.report(3,
"Loaded dark current");
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");
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;
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");
158 scaling_factor = duration / dark_duration;
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";
169 progress.report(
"Scaling dark current");
173 rebinAlg->setProperty(
"InputWorkspace", darkWS);
174 rebinAlg->setProperty(
"OutputWorkspace", darkWS);
175 rebinAlg->executeAsChildAlg();
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");
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 "
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";
200 progress.report(
"Subtracting dark current");
201 const auto &spectrumInfo = inputWS->spectrumInfo();
203 for (
int i = 0; i < numberOfSpectra; i++) {
205 if (spectrumInfo.isMasked(i))
208 const MantidVec &YDarkValues = scaledDarkWS->readY(i);
209 const MantidVec &YDarkErrors = scaledDarkWS->readE(i);
210 const MantidVec &XValues = inputWS->readX(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);
220 setProperty(
"OutputMessage",
"Dark current subtracted: " + output_message);
222 progress.report(
"Subtracted dark current");
#define DECLARE_ALGORITHM(classname)
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
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.
std::string toString() const override
Serialize an object to a string.
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 IAlgorithm_sptr fromString(const std::string &input)
De-serialize an object from a string.
void progress(double p, const std::string &msg="", double estimatedTime=0.0, int progressPrecision=0)
Sends ProgressNotification.
@ Load
allowed here which will be passed to the algorithm
Helper class for reporting progress from algorithms.
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.
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
void exec() override
Execution code.
void init() override
Initialisation code.
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
@ Input
An input workspace.
@ Output
An output workspace.