Mantid
Loading...
Searching...
No Matches
SANSBeamFinder.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 +
18
19#include "Poco/NumberFormatter.h"
20#include "Poco/Path.h"
21#include "Poco/String.h"
22
24
25// Register the algorithm into the AlgorithmFactory
26DECLARE_ALGORITHM(SANSBeamFinder)
27
28using namespace Kernel;
29using namespace API;
30using namespace Geometry;
31using namespace DataObjects;
32
34 const std::vector<std::string> fileExts{"_event.nxs", ".xml"};
35 declareProperty(std::make_unique<API::FileProperty>("Filename", "", API::FileProperty::Load, fileExts),
36 "Data filed used to find beam center");
37
38 declareProperty("BeamCenterX", EMPTY_DBL(), "Beam position in X pixel coordinates");
39 declareProperty("BeamCenterY", EMPTY_DBL(), "Beam position in Y pixel coordinates");
40
41 declareProperty("UseDirectBeamMethod", true, "If true, the direct beam method will be used");
42 declareProperty("BeamRadius", 3.0, "Beam radius in pixels, used with the scattered beam method");
43
44 declareProperty("FoundBeamCenterX", EMPTY_DBL(), Direction::Output);
45 declareProperty("FoundBeamCenterY", EMPTY_DBL(), Direction::Output);
46
47 declareProperty("PersistentCorrection", true,
48 "If true, the algorithm will be persistent and re-used when "
49 "other data sets are processed");
50 declareProperty("ReductionProperties", "__sans_reduction_properties", Direction::Input);
51 declareProperty("OutputMessage", "", Direction::Output);
52}
53
54MatrixWorkspace_sptr SANSBeamFinder::loadBeamFinderFile(const std::string &beamCenterFile) {
55 Poco::Path path(beamCenterFile);
56 const std::string entryName = "SANSBeamFinder" + path.getBaseName();
57 const std::string reductionManagerName = getProperty("ReductionProperties");
58
59 MatrixWorkspace_sptr finderWS;
60
61 if (m_reductionManager->existsProperty(entryName)) {
62 finderWS = m_reductionManager->getProperty(entryName);
63 m_output_message += " |Using existing workspace: " + finderWS->getName() + '\n';
64 } else {
65 // Load the dark current if we don't have it already
66 std::string finderWSName = "__beam_finder_" + path.getBaseName();
67
68 if (!m_reductionManager->existsProperty("LoadAlgorithm")) {
69 auto loadAlg = createChildAlgorithm("EQSANSLoad", 0.1, 0.3);
70 loadAlg->setProperty("Filename", beamCenterFile);
71 loadAlg->setProperty("NoBeamCenter", true);
72 loadAlg->setProperty("BeamCenterX", EMPTY_DBL());
73 loadAlg->setProperty("BeamCenterY", EMPTY_DBL());
74 loadAlg->setProperty("ReductionProperties", reductionManagerName);
75 loadAlg->executeAsChildAlg();
76 finderWS = loadAlg->getProperty("OutputWorkspace");
77 m_output_message += " |Loaded " + beamCenterFile + "\n";
78 std::string msg = loadAlg->getPropertyValue("OutputMessage");
79 m_output_message += " |" + Poco::replace(msg, "\n", "\n |") + "\n";
80 } else {
81 // Get load algorithm as a string so that we can create a completely
82 // new proxy and ensure that we don't overwrite existing properties
83 IAlgorithm_sptr loadAlg0 = m_reductionManager->getProperty("LoadAlgorithm");
84 const std::string loadString = loadAlg0->toString();
85 auto loadAlg = Algorithm::fromString(loadString);
86
87 loadAlg->setProperty("Filename", beamCenterFile);
88 if (loadAlg->existsProperty("NoBeamCenter"))
89 loadAlg->setProperty("NoBeamCenter", true);
90 if (loadAlg->existsProperty("BeamCenterX"))
91 loadAlg->setProperty("BeamCenterX", EMPTY_DBL());
92 if (loadAlg->existsProperty("BeamCenterY"))
93 loadAlg->setProperty("BeamCenterY", EMPTY_DBL());
94 if (loadAlg->existsProperty("ReductionProperties"))
95 loadAlg->setProperty("ReductionProperties", reductionManagerName);
96 loadAlg->setPropertyValue("OutputWorkspace", finderWSName);
97 loadAlg->execute();
98 std::shared_ptr<Workspace> wks = AnalysisDataService::Instance().retrieve(finderWSName);
99 finderWS = std::dynamic_pointer_cast<MatrixWorkspace>(wks);
100
101 m_output_message += " |Loaded " + beamCenterFile + "\n";
102 if (loadAlg->existsProperty("OutputMessage")) {
103 std::string msg = loadAlg->getPropertyValue("OutputMessage");
104 m_output_message += " |" + Poco::replace(msg, "\n", "\n |") + "\n";
105 }
106 }
107 m_reductionManager->declareProperty(std::make_unique<WorkspaceProperty<>>(entryName, "", Direction::Output));
108 m_reductionManager->setPropertyValue(entryName, finderWSName);
109 m_reductionManager->setProperty(entryName, finderWS);
110 }
111 return finderWS;
112}
113
115 // Reduction property manager
116 const std::string reductionManagerName = getProperty("ReductionProperties");
117 if (PropertyManagerDataService::Instance().doesExist(reductionManagerName)) {
118 m_reductionManager = PropertyManagerDataService::Instance().retrieve(reductionManagerName);
119 } else {
120 m_reductionManager = std::make_shared<PropertyManager>();
121 PropertyManagerDataService::Instance().addOrReplace(reductionManagerName, m_reductionManager);
122 }
123
124 const bool persistent = getProperty("PersistentCorrection");
125 if (!m_reductionManager->existsProperty("SANSBeamFinderAlgorithm") && persistent) {
126 auto algProp = std::make_unique<AlgorithmProperty>("SANSBeamFinderAlgorithm");
127 algProp->setValue(toString());
128 m_reductionManager->declareProperty(std::move(algProp));
129 }
130
131 m_output_message = "Beam center determination\n";
132
133 // Pixel coordinate to real-space coordinate mapping scheme
134 bool specialMapping = false;
135 if (m_reductionManager->existsProperty("InstrumentName")) {
136 const std::string instrumentName = m_reductionManager->getPropertyValue("InstrumentName");
137 specialMapping = instrumentName == "HFIRSANS";
138 }
139
140 // Find beam center
141 double center_x = getProperty("BeamCenterX");
142 double center_y = getProperty("BeamCenterY");
143
144 // Check whether we already know the position
145 const std::string beamCenterFile = getProperty("Filename");
146 Poco::Path path(beamCenterFile);
147 const std::string entryNameX = "SANSBeamFinder_X_" + path.getBaseName();
148 const std::string entryNameY = "SANSBeamFinder_Y_" + path.getBaseName();
149
150 // If the beam enter was provided, simply add it to the reduction
151 // property manager for other algorithms to find it
152 if (!isEmpty(center_x) && !isEmpty(center_y)) {
153 m_output_message += " |Using supplied beam center: ";
154 } else if (m_reductionManager->existsProperty(entryNameX) && m_reductionManager->existsProperty(entryNameY)) {
155 center_x = m_reductionManager->getProperty(entryNameX);
156 center_y = m_reductionManager->getProperty(entryNameY);
157 } else {
158 // Load the beam center file
159 MatrixWorkspace_sptr beamCenterWS = loadBeamFinderFile(beamCenterFile);
160
161 // HFIR reduction masks the first pixels on each edge of the detector
162 if (specialMapping)
163 // int high, int low, int left, int right
164 maskEdges(beamCenterWS, 1, 1, 1, 1);
165
166 auto ctrAlg = createChildAlgorithm("FindCenterOfMassPosition");
167 ctrAlg->setProperty("InputWorkspace", beamCenterWS);
168
169 const bool directBeam = getProperty("UseDirectBeamMethod");
170 ctrAlg->setProperty("DirectBeam", directBeam);
171
172 double beamRadius = getProperty("BeamRadius");
173 if (!directBeam && !isEmpty(beamRadius)) {
174 std::vector<double> pars = beamCenterWS->getInstrument()->getNumberParameter("x-pixel-size");
175 if (pars.empty()) {
176 g_log.error() << "Could not read pixel size from instrument "
177 "parameters: using default\n";
178 } else {
179 ctrAlg->setProperty("BeamRadius", beamRadius * pars[0] / 1000.0);
180 }
181 }
182 ctrAlg->execute();
183 std::vector<double> centerOfMass = ctrAlg->getProperty("CenterOfMass");
184
185 if (specialMapping) {
186 HFIRInstrument::getPixelFromCoordinate(centerOfMass[0], centerOfMass[1], beamCenterWS, center_x, center_y);
187 } else {
188 EQSANSInstrument::getPixelFromCoordinate(centerOfMass[0], centerOfMass[1], beamCenterWS, center_x, center_y);
189 }
190
191 m_output_message += " |Found beam center: ";
192 }
193
194 // Store for later use
195 if (persistent) {
196 if (!m_reductionManager->existsProperty("LatestBeamCenterX"))
197 m_reductionManager->declareProperty(std::make_unique<PropertyWithValue<double>>("LatestBeamCenterX", center_x));
198 if (!m_reductionManager->existsProperty("LatestBeamCenterY"))
199 m_reductionManager->declareProperty(std::make_unique<PropertyWithValue<double>>("LatestBeamCenterY", center_y));
200
201 m_reductionManager->setProperty("LatestBeamCenterX", center_x);
202 m_reductionManager->setProperty("LatestBeamCenterY", center_y);
203 }
204
206 "[" + Poco::NumberFormatter::format(center_x, 3) + ", " + Poco::NumberFormatter::format(center_y, 3) + "]\n";
207
208 // Workflow algorithms can use the LatestBeamCenterX/Y entries, but to be
209 // compatible with the old ReductionSteps we also set output properties
210 // with the beam center position
211 setProperty("FoundBeamCenterX", center_x);
212 setProperty("FoundBeamCenterY", center_y);
213 setProperty("OutputMessage", m_output_message);
214}
215
216/*
217 * The standard HFIR reduction masks the edges of the detector
218 * This is here mostly to allow a direct comparison with old HFIR code
219 * and ensure that we reproduce the same results
220 *
221 * 2016/05/06 : this only works for RectangularDetector
222 *
223 */
224void SANSBeamFinder::maskEdges(const MatrixWorkspace_sptr &beamCenterWS, int high, int low, int left, int right,
225 const std::string &componentName) {
226
227 auto instrument = beamCenterWS->getInstrument();
228
229 std::shared_ptr<Mantid::Geometry::RectangularDetector> component;
230 try {
231 component = std::const_pointer_cast<Mantid::Geometry::RectangularDetector>(
232 std::dynamic_pointer_cast<const Mantid::Geometry::RectangularDetector>(
233 instrument->getComponentByName(componentName)));
234 } catch (std::exception &) {
235 g_log.warning("Expecting the component " + componentName + " to be a RectangularDetector. maskEdges not executed.");
236 return;
237 }
238
239 if (!component) {
240 g_log.warning("Expecting the component " + componentName + " to be a RectangularDetector. maskEdges not executed.");
241 return;
242 }
243
244 std::vector<int> IDs;
245
246 // right
247 for (int i = 0; i < right * component->idstep(); i++) {
248 IDs.emplace_back(component->idstart() + i);
249 }
250 // left
251 for (int i = component->maxDetectorID(); i > (component->maxDetectorID() - left * component->idstep()); i--) {
252 IDs.emplace_back(i);
253 }
254 // low
255 // 0,256,512,768,..,1,257,513
256 for (int row = 0; row < low; row++) { // 0,1
257 for (int i = row + component->idstart();
258 i < component->nelements() * component->idstep() - component->idstep() + low + component->idstart();
259 i += component->idstep()) {
260 IDs.emplace_back(i);
261 }
262 }
263 // high
264 // 255, 511, 767..
265 for (int row = 0; row < high; row++) {
266 for (int i = component->idstep() + component->idstart() - row - 1;
267 i < component->nelements() * component->idstep() + component->idstart(); i += component->idstep()) {
268 IDs.emplace_back(i);
269 }
270 }
271
272 g_log.debug() << "SANSBeamFinder::maskEdges Detector Ids to Mask:" << std::endl;
273 for (auto id : IDs) {
274 g_log.debug() << id << " ";
275 }
276 g_log.debug() << std::endl;
277
278 auto maskAlg = createChildAlgorithm("MaskDetectors");
279 maskAlg->setChild(true);
280 maskAlg->setProperty("Workspace", beamCenterWS);
281 maskAlg->setProperty("DetectorList", IDs);
282 maskAlg->execute();
283
284 // HAcking to mask Wing Detector when finding the beam center
285 if (instrument->getComponentByName("wing_detector")) {
286 maskAlg = createChildAlgorithm("SANSMask");
287 maskAlg->setProperty("Workspace", beamCenterWS);
288 maskAlg->setPropertyValue("Facility", "HFIR");
289 maskAlg->setProperty("MaskedFullComponent", "wing_detector");
290 maskAlg->execute();
291 }
292}
293
294} // namespace Mantid::WorkflowAlgorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
double left
Definition: LineProfile.cpp:80
double right
Definition: LineProfile.cpp:81
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
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
static bool isEmpty(const NumT toCheck)
checks that the value was not set by users, uses the value in empty double/int.
@ Load
allowed here which will be passed to the algorithm
Definition: FileProperty.h:52
A property class for workspaces.
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
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
The concrete, templated class for properties.
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.
std::shared_ptr< Kernel::PropertyManager > m_reductionManager
API::MatrixWorkspace_sptr loadBeamFinderFile(const std::string &beamCenterFile)
void maskEdges(const API::MatrixWorkspace_sptr &beamCenterWS, int high, int low, int left, int right, const std::string &componentName="detector1")
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
void getPixelFromCoordinate(const double &x, const double &y, const API::MatrixWorkspace_sptr &dataWS, double &pixel_x, double &pixel_y)
void getPixelFromCoordinate(const double &x, const double &y, const API::MatrixWorkspace_sptr &dataWS, double &pixel_x, double &pixel_y)
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
Definition: EmptyValues.h:43
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54