Mantid
Loading...
Searching...
No Matches
SANSSensitivityCorrection.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 +
13#include "MantidAPI/Run.h"
14#include "MantidAPI/TableRow.h"
23#include "Poco/File.h"
24#include "Poco/NumberFormatter.h"
25#include "Poco/Path.h"
26#include "Poco/String.h"
28
29// Register the algorithm into the AlgorithmFactory
30DECLARE_ALGORITHM(SANSSensitivityCorrection)
31
32using namespace Kernel;
33using namespace API;
34using namespace Geometry;
35using namespace DataObjects;
36
39 std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input, PropertyMode::Optional));
40 const std::vector<std::string> fileExts{"_event.nxs", ".xml"};
41 declareProperty(std::make_unique<API::FileProperty>("Filename", "", API::FileProperty::Load, fileExts),
42 "Flood field or sensitivity file.");
43 declareProperty("UseSampleDC", true,
44 "If true, the dark current subtracted "
45 "from the sample data will also be "
46 "subtracted from the flood field.");
47 declareProperty(std::make_unique<API::FileProperty>("DarkCurrentFile", "", API::FileProperty::OptionalLoad, fileExts),
48 "The name of the input file to load as dark current.");
49
50 auto positiveDouble = std::make_shared<BoundedValidator<double>>();
51 positiveDouble->setLower(0);
52 declareProperty("MinEfficiency", EMPTY_DBL(), positiveDouble,
53 "Minimum efficiency for a pixel to be considered (default: no minimum).");
54 declareProperty("MaxEfficiency", EMPTY_DBL(), positiveDouble,
55 "Maximum efficiency for a pixel to be considered (default: no maximum).");
56
57 declareProperty("FloodTransmissionValue", EMPTY_DBL(), positiveDouble,
58 "Transmission value for the flood field material "
59 "(default: no transmission).");
60 declareProperty("FloodTransmissionError", 0.0, positiveDouble,
61 "Transmission error for the flood field material "
62 "(default: no transmission).");
63
64 declareProperty("BeamCenterX", EMPTY_DBL(),
65 "Beam position in X pixel coordinates (optional: otherwise "
66 "sample beam center is used)");
67 declareProperty("BeamCenterY", EMPTY_DBL(),
68 "Beam position in Y pixel coordinates (optional: otherwise "
69 "sample beam center is used)");
70 declareProperty("MaskedFullComponent", "", "Component Name to fully mask according to the IDF file.");
71 declareProperty(std::make_unique<ArrayProperty<int>>("MaskedEdges"),
72 "Number of pixels to mask on the edges: X-low, X-high, Y-low, Y-high");
73 declareProperty("MaskedComponent", "", "Component Name to mask the edges according to the IDF file.");
74
76 std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output, PropertyMode::Optional));
77 declareProperty("ReductionProperties", "__sans_reduction_properties");
78 declareProperty(std::make_unique<WorkspaceProperty<MatrixWorkspace>>("OutputSensitivityWorkspace", "",
80 declareProperty("OutputMessage", "", Direction::Output);
81}
82
87bool SANSSensitivityCorrection::fileCheck(const std::string &filePath) {
88 // Check the file extension
89 Poco::Path path(filePath);
90 const std::string extn = path.getExtension();
91 const std::string nxs("nxs");
92 const std::string nx5("nx5");
93 if (!(Poco::icompare(nxs, extn) == 0 || Poco::icompare(nx5, extn) == 0))
94 return false;
95
96 // If we have a Nexus file, check that is comes from Mantid
97 std::vector<std::string> entryName, definition;
98 int count = NeXus::getNexusEntryTypes(filePath, entryName, definition);
99 if (count <= -1) {
100 g_log.error("Error reading file " + filePath);
101 throw Exception::FileError("Unable to read data in File:", filePath);
102 } else if (count == 0) {
103 g_log.error("Error no entries found in " + filePath);
104 return false;
105 }
106
107 return entryName[0] == "mantid_workspace_1";
108}
109
111 // Output log
112 m_output_message = "";
113
114 Progress progress(this, 0.0, 1.0, 10);
115
116 // Reduction property manager
117 const std::string reductionManagerName = getProperty("ReductionProperties");
118 std::shared_ptr<PropertyManager> reductionManager;
119 if (PropertyManagerDataService::Instance().doesExist(reductionManagerName)) {
120 reductionManager = PropertyManagerDataService::Instance().retrieve(reductionManagerName);
121 } else {
122 reductionManager = std::make_shared<PropertyManager>();
123 PropertyManagerDataService::Instance().addOrReplace(reductionManagerName, reductionManager);
124 }
125
126 if (!reductionManager->existsProperty("SensitivityAlgorithm")) {
127 auto algProp = std::make_unique<AlgorithmProperty>("SensitivityAlgorithm");
128 algProp->setValue(toString());
129 reductionManager->declareProperty(std::move(algProp));
130 }
131
132 progress.report("Loading sensitivity file");
133 const std::string fileName = getPropertyValue("Filename");
134
135 // Look for an entry for the dark current in the reduction table
136 Poco::Path path(fileName);
137 const std::string entryName = "Sensitivity" + path.getBaseName();
138 MatrixWorkspace_sptr floodWS;
139 std::string floodWSName = "__sensitivity_" + path.getBaseName();
140
141 if (reductionManager->existsProperty(entryName)) {
142 std::string wsName = reductionManager->getPropertyValue(entryName);
143 floodWS =
144 std::dynamic_pointer_cast<MatrixWorkspace>(AnalysisDataService::Instance().retrieveWS<MatrixWorkspace>(wsName));
145 m_output_message += " |Using " + wsName + "\n";
146 g_log.debug() << "SANSSensitivityCorrection :: Using sensitivity workspace: " << wsName << "\n";
147 } else {
148 // Load the flood field if we don't have it already
149 // First, try to determine whether we need to load data or a sensitivity
150 // workspace...
151 if (fileCheck(fileName)) {
152 g_log.debug() << "SANSSensitivityCorrection :: Loading sensitivity file: " << fileName << "\n";
153 auto loadAlg = createChildAlgorithm("Load", 0.1, 0.3);
154 loadAlg->setProperty("Filename", fileName);
155 loadAlg->executeAsChildAlg();
156 Workspace_sptr floodWS_ws = loadAlg->getProperty("OutputWorkspace");
157 floodWS = std::dynamic_pointer_cast<MatrixWorkspace>(floodWS_ws);
158
159 // Check that it's really a sensitivity file
160 if (!floodWS->run().hasProperty("is_sensitivity")) {
161 // Reset pointer
162 floodWS.reset();
163 g_log.error() << "A processed Mantid workspace was loaded but it "
164 "wasn't a sensitivity file!\n";
165 }
166 }
167
168 // ... if we don't, just load the data and process it
169 if (!floodWS) {
170 // Read in default beam center
171 double center_x = getProperty("BeamCenterX");
172 double center_y = getProperty("BeamCenterY");
173 if (isEmpty(center_x) || isEmpty(center_y)) {
174 if (reductionManager->existsProperty("LatestBeamCenterX") &&
175 reductionManager->existsProperty("LatestBeamCenterY")) {
176 center_x = reductionManager->getProperty("LatestBeamCenterX");
177 center_y = reductionManager->getProperty("LatestBeamCenterY");
178 m_output_message += " |Setting beam center to [" + Poco::NumberFormatter::format(center_x, 1) + ", " +
179 Poco::NumberFormatter::format(center_y, 1) + "]\n";
180 } else
181 m_output_message += " |No beam center provided: skipping!\n";
182 }
183
184 const std::string rawFloodWSName = "__flood_data_" + path.getBaseName();
185 MatrixWorkspace_sptr rawFloodWS;
186 if (!reductionManager->existsProperty("LoadAlgorithm")) {
187 auto loadAlg = createChildAlgorithm("Load", 0.1, 0.3);
188 loadAlg->setProperty("Filename", fileName);
189 if (!isEmpty(center_x) && loadAlg->existsProperty("BeamCenterX"))
190 loadAlg->setProperty("BeamCenterX", center_x);
191 if (!isEmpty(center_y) && loadAlg->existsProperty("BeamCenterY"))
192 loadAlg->setProperty("BeamCenterY", center_y);
193 loadAlg->setPropertyValue("OutputWorkspace", rawFloodWSName);
194 loadAlg->executeAsChildAlg();
195 Workspace_sptr tmpWS = loadAlg->getProperty("OutputWorkspace");
196 rawFloodWS = std::dynamic_pointer_cast<MatrixWorkspace>(tmpWS);
197 m_output_message += " | Loaded " + fileName + " (Load algorithm)\n";
198 } else {
199 // Get load algorithm as a string so that we can create a completely
200 // new proxy and ensure that we don't overwrite existing properties
201 IAlgorithm_sptr loadAlg0 = reductionManager->getProperty("LoadAlgorithm");
202 const std::string loadString = loadAlg0->toString();
203 auto loadAlg = Algorithm::fromString(loadString);
204 loadAlg->setChild(true);
205 loadAlg->setProperty("Filename", fileName);
206 loadAlg->setPropertyValue("OutputWorkspace", rawFloodWSName);
207 if (!isEmpty(center_x) && loadAlg->existsProperty("BeamCenterX"))
208 loadAlg->setProperty("BeamCenterX", center_x);
209 if (!isEmpty(center_y) && loadAlg->existsProperty("BeamCenterY"))
210 loadAlg->setProperty("BeamCenterY", center_y);
211 loadAlg->execute();
212 rawFloodWS = loadAlg->getProperty("OutputWorkspace");
213 m_output_message += " |Loaded " + fileName + "\n";
214 if (loadAlg->existsProperty("OutputMessage")) {
215 std::string msg = loadAlg->getPropertyValue("OutputMessage");
216 m_output_message += " |" + Poco::replace(msg, "\n", "\n |") + "\n";
217 }
218 }
219
220 // Check whether we just loaded a flood field data set, or the actual
221 // sensitivity
222 if (!rawFloodWS->run().hasProperty("is_sensitivity")) {
223 const std::string darkCurrentFile = getPropertyValue("DarkCurrentFile");
224
225 // Look for a dark current subtraction algorithm
226 std::string dark_result;
227 if (reductionManager->existsProperty("DarkCurrentAlgorithm")) {
228 IAlgorithm_sptr darkAlg = reductionManager->getProperty("DarkCurrentAlgorithm");
229 darkAlg->setChild(true);
230 darkAlg->setProperty("InputWorkspace", rawFloodWS);
231 darkAlg->setProperty("OutputWorkspace", rawFloodWS);
232
233 // Execute as-is if we use the sample dark current, otherwise check
234 // whether a dark current file was provided.
235 // Otherwise do nothing
236 if (getProperty("UseSampleDC")) {
237 darkAlg->execute();
238 if (darkAlg->existsProperty("OutputMessage"))
239 dark_result = darkAlg->getPropertyValue("OutputMessage");
240 } else if (!darkCurrentFile.empty()) {
241 darkAlg->setProperty("Filename", darkCurrentFile);
242 darkAlg->setProperty("PersistentCorrection", false);
243 darkAlg->execute();
244 if (darkAlg->existsProperty("OutputMessage"))
245 dark_result = darkAlg->getPropertyValue("OutputMessage");
246 else
247 dark_result = " Dark current subtracted\n";
248 }
249 } else if (!darkCurrentFile.empty()) {
250 // We need to subtract the dark current for the flood field but no
251 // dark
252 // current subtraction was set for the sample! Use the default dark
253 // current algorithm if we can find it.
254 if (reductionManager->existsProperty("DefaultDarkCurrentAlgorithm")) {
255 IAlgorithm_sptr darkAlg = reductionManager->getProperty("DefaultDarkCurrentAlgorithm");
256 darkAlg->setChild(true);
257 darkAlg->setProperty("InputWorkspace", rawFloodWS);
258 darkAlg->setProperty("OutputWorkspace", rawFloodWS);
259 darkAlg->setProperty("Filename", darkCurrentFile);
260 darkAlg->setProperty("PersistentCorrection", false);
261 darkAlg->execute();
262 if (darkAlg->existsProperty("OutputMessage"))
263 dark_result = darkAlg->getPropertyValue("OutputMessage");
264 } else {
265 // We are running out of options
266 g_log.error() << "No dark current algorithm provided to load [" << getPropertyValue("DarkCurrentFile")
267 << "]: skipped!\n";
268 dark_result = " No dark current algorithm provided: skipped\n";
269 }
270 }
271 m_output_message += " |" + Poco::replace(dark_result, "\n", "\n |") + "\n";
272
273 // Look for solid angle correction algorithm
274 if (reductionManager->existsProperty("SANSSolidAngleCorrection")) {
275 IAlgorithm_sptr solidAlg = reductionManager->getProperty("SANSSolidAngleCorrection");
276 solidAlg->setChild(true);
277 solidAlg->setProperty("InputWorkspace", rawFloodWS);
278 solidAlg->setProperty("OutputWorkspace", rawFloodWS);
279 solidAlg->execute();
280 std::string msg = "Solid angle correction applied\n";
281 if (solidAlg->existsProperty("OutputMessage"))
282 msg = solidAlg->getPropertyValue("OutputMessage");
283 m_output_message += " |" + Poco::replace(msg, "\n", "\n |") + "\n";
284 }
285
286 // Apply transmission correction as needed
287 double floodTransmissionValue = getProperty("FloodTransmissionValue");
288 double floodTransmissionError = getProperty("FloodTransmissionError");
289
290 if (!isEmpty(floodTransmissionValue)) {
291 g_log.debug() << "SANSSensitivityCorrection :: Applying transmission "
292 "to flood field\n";
293 auto transAlg = createChildAlgorithm("ApplyTransmissionCorrection");
294 transAlg->setProperty("InputWorkspace", rawFloodWS);
295 transAlg->setProperty("OutputWorkspace", rawFloodWS);
296 transAlg->setProperty("TransmissionValue", floodTransmissionValue);
297 transAlg->setProperty("TransmissionError", floodTransmissionError);
298 transAlg->setProperty("ThetaDependent", true);
299 transAlg->execute();
300 rawFloodWS = transAlg->getProperty("OutputWorkspace");
301 m_output_message += " |Applied transmission to flood field\n";
302 }
303
304 // Calculate detector sensitivity
305 auto effAlg = createChildAlgorithm("CalculateEfficiency", -1, -1, true, 1);
306 effAlg->setProperty("InputWorkspace", rawFloodWS);
307
308 const double minEff = getProperty("MinEfficiency");
309 const double maxEff = getProperty("MaxEfficiency");
310 const std::string maskFullComponent = getPropertyValue("MaskedFullComponent");
311 const std::string maskEdges = getPropertyValue("MaskedEdges");
312 const std::string maskComponent = getPropertyValue("MaskedComponent");
313
314 effAlg->setProperty("MinEfficiency", minEff);
315 effAlg->setProperty("MaxEfficiency", maxEff);
316 effAlg->setProperty("MaskedFullComponent", maskFullComponent);
317 effAlg->setProperty("MaskedEdges", maskEdges);
318 effAlg->setProperty("MaskedComponent", maskComponent);
319 effAlg->execute();
320 floodWS = effAlg->getProperty("OutputWorkspace");
321 } else {
322 floodWS = rawFloodWS;
323 }
324 // Patch as needed
325 if (reductionManager->existsProperty("SensitivityPatchAlgorithm")) {
326 IAlgorithm_sptr patchAlg = reductionManager->getProperty("SensitivityPatchAlgorithm");
327 patchAlg->setChild(true);
328 patchAlg->setProperty("Workspace", floodWS);
329 patchAlg->execute();
330 m_output_message += " |Sensitivity patch applied\n";
331 }
332
333 floodWS->mutableRun().addProperty("is_sensitivity", 1, "", true);
334 }
335 std::string floodWSOutputName = getPropertyValue("OutputSensitivityWorkspace");
336 if (floodWSOutputName.empty()) {
337 setPropertyValue("OutputSensitivityWorkspace", floodWSName);
338 AnalysisDataService::Instance().addOrReplace(floodWSName, floodWS);
339 reductionManager->declareProperty(
340 std::make_unique<WorkspaceProperty<>>(entryName, floodWSName, Direction::InOut));
341 reductionManager->setPropertyValue(entryName, floodWSName);
342 reductionManager->setProperty(entryName, floodWS);
343 }
344 setProperty("OutputSensitivityWorkspace", floodWS);
345 }
346
347 progress.report(3, "Loaded flood field");
348
349 // Check whether we need to apply the correction to a workspace
350 MatrixWorkspace_sptr inputWS = getProperty("InputWorkspace");
351 if (inputWS) {
352 // Divide sample data by detector efficiency
353 auto divideAlg = createChildAlgorithm("Divide", 0.6, 0.7);
354 divideAlg->setProperty("LHSWorkspace", inputWS);
355 divideAlg->setProperty("RHSWorkspace", floodWS);
356 divideAlg->executeAsChildAlg();
357 MatrixWorkspace_sptr outputWS = divideAlg->getProperty("OutputWorkspace");
358
359 // Copy over the efficiency's masked pixels to the reduced workspace
360 auto maskAlg = createChildAlgorithm("MaskDetectors", 0.75, 0.85);
361 maskAlg->setProperty("Workspace", outputWS);
362 maskAlg->setProperty("MaskedWorkspace", floodWS);
363 maskAlg->executeAsChildAlg();
364
365 setProperty("OutputWorkspace", outputWS);
366 }
367 setProperty("OutputMessage", "Sensitivity correction computed\n" + m_output_message);
368
369 progress.report("Performed sensitivity correction");
370}
371
372} // namespace Mantid::WorkflowAlgorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
int count
counter
Definition: Matrix.cpp:37
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
void setPropertyValue(const std::string &name, const std::string &value) override
Set the value of a property by string N.B.
Definition: Algorithm.cpp:1975
static bool isEmpty(const NumT toCheck)
checks that the value was not set by users, uses the value in empty double/int.
@ OptionalLoad
to specify a file to read but the file doesn't have to exist
Definition: FileProperty.h:53
@ 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.
Support for a property that holds an array of values.
Definition: ArrayProperty.h:28
Records the filename and the description of failure.
Definition: Exception.h:98
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void debug(const std::string &msg)
Logs at debug level.
Definition: Logger.cpp:114
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...
bool fileCheck(const std::string &filePath)
Check whether we have a processed file of not.
std::shared_ptr< IAlgorithm > IAlgorithm_sptr
shared pointer to Mantid::API::IAlgorithm
std::shared_ptr< Workspace > Workspace_sptr
shared pointer to Mantid::API::Workspace
Definition: Workspace_fwd.h:20
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
MANTID_NEXUS_DLL int getNexusEntryTypes(const std::string &fileName, std::vector< std::string > &entryName, std::vector< std::string > &definition)
Get all the Nexus entry types for a file.
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
Definition: EmptyValues.h:43
@ 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