Mantid
Loading...
Searching...
No Matches
EQSANSLoad.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 +
9#include "MantidAPI/Axis.h"
11#include "MantidAPI/Run.h"
19
20#include <boost/regex.hpp>
21
22#include "Poco/DirectoryIterator.h"
23#include "Poco/NumberFormatter.h"
24#include "Poco/NumberParser.h"
25#include "Poco/String.h"
26
27#include <fstream>
28
30
31// Register the algorithm into the AlgorithmFactory
32DECLARE_ALGORITHM(EQSANSLoad)
33
34using namespace Kernel;
35using namespace API;
36using namespace Geometry;
37using namespace DataObjects;
38
40 declareProperty(std::make_unique<API::FileProperty>("Filename", "", API::FileProperty::OptionalLoad, "_event.nxs"),
41 "The name of the input event Nexus file to load");
42
43 auto wsValidator = std::make_shared<WorkspaceUnitValidator>("TOF");
44 declareProperty(std::make_unique<WorkspaceProperty<EventWorkspace>>("InputWorkspace", "", Direction::Input,
45 PropertyMode::Optional, wsValidator),
46 "Input event workspace. Assumed to be unmodified events "
47 "straight from LoadEventNexus");
48
49 declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
50 "Then name of the output EventWorkspace");
51 declareProperty("NoBeamCenter", false, "If true, the detector will not be moved according to the beam center");
52 declareProperty("UseConfigBeam", false,
53 "If true, the beam center defined in "
54 "the configuration file will be "
55 "used");
56 declareProperty("BeamCenterX", EMPTY_DBL(),
57 "Beam position in X pixel "
58 "coordinates (used only if "
59 "UseConfigBeam is false)");
60 declareProperty("BeamCenterY", EMPTY_DBL(),
61 "Beam position in Y pixel "
62 "coordinates (used only if "
63 "UseConfigBeam is false)");
64 declareProperty("UseConfigTOFCuts", false,
65 "If true, the edges of the TOF distribution will be cut "
66 "according to the configuration file");
67 declareProperty("LowTOFCut", 0.0,
68 "TOF value below which events will not be "
69 "loaded into the workspace at load-time");
70 declareProperty("HighTOFCut", 0.0,
71 "TOF value above which events will not be "
72 "loaded into the workspace at load-time");
73 declareProperty("SkipTOFCorrection", false, "IF true, the EQSANS TOF correction will be skipped");
74 declareProperty("WavelengthStep", 0.1,
75 "Wavelength steps to be used when "
76 "rebinning the data before performing "
77 "the reduction");
78 declareProperty("UseConfigMask", false,
79 "If true, the masking information "
80 "found in the configuration file "
81 "will be used");
82 declareProperty("UseConfig", true, "If true, the best configuration file found will be used");
83 declareProperty("CorrectForFlightPath", false,
84 "If true, the TOF will be modified for the true flight path "
85 "from the sample to the detector pixel");
86 declareProperty("PreserveEvents", true, "If true, the output workspace will be an event workspace");
87 declareProperty("SampleDetectorDistance", EMPTY_DBL(),
88 "Sample to detector distance to use (overrides meta data), in mm");
89 declareProperty("SampleDetectorDistanceOffset", EMPTY_DBL(),
90 "Offset to the sample to detector distance (use only when "
91 "using the distance found in the meta data), in mm");
92 declareProperty("SampleOffset", EMPTY_DBL(),
93 "Offset to be applied to the sample position (use only when "
94 "using the detector distance found in the meta data), in mm");
95 declareProperty("DetectorOffset", EMPTY_DBL(),
96 "Offset to be applied to the detector position (use only when "
97 "using the distance found in the meta data), in mm");
98 declareProperty("LoadMonitors", true, "If true, the monitor workspace will be loaded");
99 declareProperty("LoadNexusInstrumentXML", true,
100 "Reads the embedded Instrument XML from the NeXus file "
101 "(optional, default True). ");
102 declareProperty("OutputMessage", "", Direction::Output);
103 declareProperty("ReductionProperties", "__sans_reduction_properties", Direction::Input);
104}
105
109double getRunPropertyDbl(const MatrixWorkspace_sptr &inputWS, const std::string &pname) {
110 Mantid::Kernel::Property *prop = inputWS->run().getProperty(pname);
111 auto *dp = dynamic_cast<Mantid::Kernel::PropertyWithValue<double> *>(prop);
112 if (!dp)
113 throw std::runtime_error("Could not cast (interpret) the property " + pname +
114 " as a floating point numeric value.");
115 return *dp;
116}
117
120std::string EQSANSLoad::findConfigFile(const int &run) {
121 // Append the standard location of EQSANS config file to the data search
122 // directory list
123 std::string sns_folder = "/SNS/EQSANS/shared/instrument_configuration";
124 if (Poco::File(sns_folder).exists())
125 Kernel::ConfigService::Instance().appendDataSearchDir(sns_folder);
126
127 const std::vector<std::string> &searchPaths = Kernel::ConfigService::Instance().getDataSearchDirs();
128
129 int max_run_number = 0;
130 std::string config_file;
131 static boost::regex re1("eqsans_configuration\\.([0-9]+)$");
132 boost::smatch matches;
133 for (const auto &searchPath : searchPaths) {
134 if (Poco::File(searchPath).exists()) {
135 Poco::DirectoryIterator file_it(searchPath);
136 Poco::DirectoryIterator end;
137 for (; file_it != end; ++file_it) {
138 if (boost::regex_search(file_it.name(), matches, re1)) {
139 std::string s = matches[1];
140 int run_number = 0;
141 Poco::NumberParser::tryParse(s, run_number);
142 if (run_number > max_run_number && run_number <= run) {
143 max_run_number = run_number;
144 config_file = file_it.path().toString();
145 }
146 }
147 }
148 }
149 }
150 return config_file;
151}
152
155void EQSANSLoad::readRectangularMasks(const std::string &line) {
156 // Looking for rectangular mask
157 // Rectangular mask = 7, 0; 7, 255
158 boost::regex re_key("rectangular mask", boost::regex::icase);
159 boost::regex re_key_alt("elliptical mask", boost::regex::icase);
160 if (boost::regex_search(line, re_key) || boost::regex_search(line, re_key_alt)) {
161 boost::regex re_sig("=[ ]*([0-9]+)[ ]*[ ,][ ]*([0-9]+)[ ]*[ ;,][ "
162 "]*([0-9]+)[ ]*[ ,][ ]*([0-9]+)");
163 boost::smatch posVec;
164 if (boost::regex_search(line, posVec, re_sig)) {
165 if (posVec.size() == 5) {
166 for (int i = 0; i < 4; i++) {
167 std::string num_str = posVec[i + 1];
168 m_mask_as_string = m_mask_as_string + " " + num_str;
169 }
170 m_mask_as_string += ",";
171 }
172 }
173 }
174}
175
178void EQSANSLoad::readTOFcuts(const std::string &line) {
179 boost::regex re_key("tof edge discard", boost::regex::icase);
180 if (boost::regex_search(line, re_key)) {
181 boost::regex re_sig("=[ ]*([0-9]+)[ ]*[ ,][ ]*([0-9]+)");
182 boost::smatch posVec;
183 if (boost::regex_search(line, posVec, re_sig)) {
184 if (posVec.size() == 3) {
185 std::string num_str = posVec[1];
186 Poco::NumberParser::tryParseFloat(num_str, m_low_TOF_cut);
187 num_str = posVec[2];
188 Poco::NumberParser::tryParseFloat(num_str, m_high_TOF_cut);
189 }
190 }
191 }
192}
193
196void EQSANSLoad::readBeamCenter(const std::string &line) {
197 boost::regex re_key("spectrum center", boost::regex::icase);
198 if (boost::regex_search(line, re_key)) {
199 boost::regex re_sig("=[ ]*([0-9]+.[0-9]*)[ ]*[ ,][ ]*([0-9]+.[0-9]+)");
200 boost::smatch posVec;
201 if (boost::regex_search(line, posVec, re_sig)) {
202 if (posVec.size() == 3) {
203 std::string num_str = posVec[1];
204 Poco::NumberParser::tryParseFloat(num_str, m_center_x);
205 num_str = posVec[2];
206 Poco::NumberParser::tryParseFloat(num_str, m_center_y);
207 }
208 }
209 }
210}
211
214void EQSANSLoad::readModeratorPosition(const std::string &line) {
215 boost::regex re_key("sample location", boost::regex::icase);
216 if (boost::regex_search(line, re_key)) {
217 boost::regex re_sig("=[ ]*([0-9]+)");
218 boost::smatch posVec;
219 if (boost::regex_search(line, posVec, re_sig)) {
220 if (posVec.size() == 2) {
221 std::string num_str = posVec[1];
222 Poco::NumberParser::tryParseFloat(num_str, m_moderator_position);
224 }
225 }
226 }
227}
228
231void EQSANSLoad::readSourceSlitSize(const std::string &line) {
232 boost::regex re_key("wheel", boost::regex::icase);
233 if (boost::regex_search(line, re_key)) {
234 boost::regex re_sig(R"(([1-8]) wheel[ ]*([1-3])[ \t]*=[ \t]*(\w+))");
235 boost::smatch posVec;
236 if (boost::regex_search(line, posVec, re_sig)) {
237 if (posVec.size() == 4) {
238 std::string num_str = posVec[1];
239 int slit_number = 0;
240 Poco::NumberParser::tryParse(num_str, slit_number);
241 slit_number--;
242
243 num_str = posVec[2];
244 int wheel_number = 0;
245 Poco::NumberParser::tryParse(num_str, wheel_number);
246 wheel_number--;
247
248 num_str = posVec[3];
249 boost::regex re_size("\\w*?([0-9]+)mm");
250 int slit_size = 0;
251 if (boost::regex_search(num_str, posVec, re_size)) {
252 if (posVec.size() == 2) {
253 num_str = posVec[1];
254 Poco::NumberParser::tryParse(num_str, slit_size);
255 }
256 }
257 m_slit_positions[wheel_number][slit_number] = slit_size;
258 }
259 }
260 }
261}
262
265 if (!dataWS->run().hasProperty("vBeamSlit")) {
266 m_output_message += " Could not determine source aperture diameter: ";
267 m_output_message += "slit parameters were not found in the run log\n";
268 return;
269 }
270
271 const std::string slit1Name = "vBeamSlit";
272 Mantid::Kernel::Property *prop = dataWS->run().getProperty(slit1Name);
273 auto *dp = dynamic_cast<Mantid::Kernel::TimeSeriesProperty<double> *>(prop);
274 auto *ip = dynamic_cast<Mantid::Kernel::TimeSeriesProperty<int> *>(prop);
275 int slit1;
276 if (dp)
277 slit1 = static_cast<int>(dp->getStatistics().mean);
278 else if (ip)
279 slit1 = static_cast<int>(ip->getStatistics().mean);
280 else
281 throw std::runtime_error("Could not cast (interpret) the property " + slit1Name +
282 " as a time series property with "
283 "int or floating point values.");
284
285 const std::string slit2Name = "vBeamSlit2";
286 prop = dataWS->run().getProperty(slit2Name);
287 dp = dynamic_cast<Mantid::Kernel::TimeSeriesProperty<double> *>(prop);
288 ip = dynamic_cast<Mantid::Kernel::TimeSeriesProperty<int> *>(prop);
289 int slit2;
290 if (dp)
291 slit2 = static_cast<int>(dp->getStatistics().mean);
292 else if (ip)
293 slit2 = static_cast<int>(ip->getStatistics().mean);
294 else
295 throw std::runtime_error("Could not cast (interpret) the property " + slit2Name +
296 " as a time series property with "
297 "int or floating point values.");
298
299 const std::string slit3Name = "vBeamSlit3";
300 prop = dataWS->run().getProperty(slit3Name);
301 dp = dynamic_cast<Mantid::Kernel::TimeSeriesProperty<double> *>(prop);
302 ip = dynamic_cast<Mantid::Kernel::TimeSeriesProperty<int> *>(prop);
303 int slit3;
304 if (dp)
305 slit3 = static_cast<int>(dp->getStatistics().mean);
306 else if (ip)
307 slit3 = static_cast<int>(ip->getStatistics().mean);
308 else
309 throw std::runtime_error("Could not cast (interpret) the property " + slit3Name +
310 " as a time series property with "
311 "int or floating point values.");
312
313 if (slit1 < 0 && slit2 < 0 && slit3 < 0) {
314 m_output_message += " Could not determine source aperture diameter\n";
315 return;
316 }
317
318 // Default slit size
319 double S1 = 20.0;
320 double L1 = -1.0;
321 const double ssd = fabs(dataWS->getInstrument()->getSource()->getPos().Z()) * 1000.0;
322 int slits[3] = {slit1, slit2, slit3};
323 for (int i = 0; i < 3; i++) {
324 int m = slits[i] - 1;
325 if (m >= 0 && m < 6) {
326 double x = m_slit_positions[i][m];
327 double y = ssd - m_slit_to_source[i];
328 if (L1 < 0 || x / y < S1 / L1) {
329 L1 = y;
330 S1 = x;
331 }
332 }
333 }
334 dataWS->mutableRun().addProperty("source-aperture-diameter", S1, "mm", true);
335 m_output_message += " Source aperture diameter: ";
336 Poco::NumberFormatter::append(m_output_message, S1, 1);
337 m_output_message += " mm\n";
338}
339
342 // Check that we have a beam center defined, otherwise set the
343 // default beam center
346 g_log.information() << "Setting beam center to [" << m_center_x << ", " << m_center_y << "]\n";
347 return;
348 }
349
350 // Check that the center of the detector really is at (0,0)
351 int nx_pixels = static_cast<int>(dataWS->getInstrument()->getNumberParameter("number-of-x-pixels")[0]);
352 int ny_pixels = static_cast<int>(dataWS->getInstrument()->getNumberParameter("number-of-y-pixels")[0]);
353 const auto &detInfo = dataWS->detectorInfo();
354 const V3D pixel_first = detInfo.position(detInfo.indexOf(0));
355 int detIDx = EQSANSInstrument::getDetectorFromPixel(nx_pixels - 1, 0, dataWS);
356 int detIDy = EQSANSInstrument::getDetectorFromPixel(0, ny_pixels - 1, dataWS);
357
358 const V3D pixel_last_x = detInfo.position(detInfo.indexOf(detIDx));
359 const V3D pixel_last_y = detInfo.position(detInfo.indexOf(detIDy));
360 double x_offset = (pixel_first.X() + pixel_last_x.X()) / 2.0;
361 double y_offset = (pixel_first.Y() + pixel_last_y.Y()) / 2.0;
362 double beam_ctr_x = 0.0;
363 double beam_ctr_y = 0.0;
365
366 auto mvAlg = createChildAlgorithm("MoveInstrumentComponent", 0.5, 0.50);
367 mvAlg->setProperty<MatrixWorkspace_sptr>("Workspace", dataWS);
368 mvAlg->setProperty("ComponentName", "detector1");
369 mvAlg->setProperty("X", -x_offset - beam_ctr_x);
370 mvAlg->setProperty("Y", -y_offset - beam_ctr_y);
371 mvAlg->setProperty("RelativePosition", true);
372 mvAlg->executeAsChildAlg();
373 m_output_message += " Beam center offset: " + Poco::NumberFormatter::format(x_offset) + ", " +
374 Poco::NumberFormatter::format(y_offset) + " m\n";
375 // m_output_message += " Beam center in real-space: " +
376 // Poco::NumberFormatter::format(-x_offset-beam_ctr_x)
377 // + ", " + Poco::NumberFormatter::format(-y_offset-beam_ctr_y) + " m\n";
378 g_log.information() << "Moving beam center to " << m_center_x << " " << m_center_y << '\n';
379
380 dataWS->mutableRun().addProperty("beam_center_x", m_center_x, "pixel", true);
381 dataWS->mutableRun().addProperty("beam_center_y", m_center_y, "pixel", true);
382 m_output_message += " Beam center: " + Poco::NumberFormatter::format(m_center_x) + ", " +
383 Poco::NumberFormatter::format(m_center_y) + "\n";
384}
385
388void EQSANSLoad::readConfigFile(const std::string &filePath) {
389 // Initialize parameters
390 m_mask_as_string = "";
392
393 // The following should be properties
394 bool use_config_mask = getProperty("UseConfigMask");
395 bool use_config_cutoff = getProperty("UseConfigTOFCuts");
396 bool use_config_center = getProperty("UseConfigBeam");
397
398 std::ifstream file(filePath.c_str());
399 if (!file) {
400 g_log.error() << "Unable to open file: " << filePath << '\n';
401 throw Exception::FileError("Unable to open file: ", filePath);
402 }
403 g_log.information() << "Using config file: " << filePath << '\n';
404 m_output_message += " Using configuration file: " + filePath + "\n";
405
406 std::string line;
407 while (getline(file, line)) {
408 boost::trim(line);
409 std::string comment = line.substr(0, 1);
410 if (Poco::icompare(comment, "#") == 0)
411 continue;
412 if (use_config_mask)
414 if (use_config_cutoff)
415 readTOFcuts(line);
416 if (use_config_center)
417 readBeamCenter(line);
419 readSourceSlitSize(line);
420 }
421
422 if (use_config_mask)
423 dataWS->mutableRun().addProperty("rectangular_masks", m_mask_as_string, "pixels", true);
424
425 dataWS->mutableRun().addProperty("low_tof_cut", m_low_TOF_cut, "microsecond", true);
426 dataWS->mutableRun().addProperty("high_tof_cut", m_high_TOF_cut, "microsecond", true);
427 m_output_message += " Discarding lower " + Poco::NumberFormatter::format(m_low_TOF_cut) + " and upper " +
428 Poco::NumberFormatter::format(m_high_TOF_cut) + " microsec\n";
429
430 if (m_moderator_position != 0)
431 dataWS->mutableRun().addProperty("moderator_position", m_moderator_position, "mm", true);
432}
433
435 // Verify the validity of the inputs
436 // TODO: this should be done by the new data management algorithm used for
437 // live data reduction (when it's implemented...)
438 const std::string fileName = getPropertyValue("Filename");
439 EventWorkspace_sptr inputEventWS = getProperty("InputWorkspace");
440 if (fileName.empty() && !inputEventWS) {
441 g_log.error() << "EQSANSLoad input error: Either a valid file path or an "
442 "input workspace must be provided\n";
443 throw std::runtime_error("EQSANSLoad input error: Either a valid file path "
444 "or an input workspace must be provided");
445 } else if (!fileName.empty() && inputEventWS) {
446 g_log.error() << "EQSANSLoad input error: Either a valid file path or an "
447 "input workspace must be provided, but not both\n";
448 throw std::runtime_error("EQSANSLoad input error: Either a valid file path "
449 "or an input workspace must be provided, but not "
450 "both");
451 }
452
453 // Read in default TOF cuts
454 const bool skipTOFCorrection = getProperty("SkipTOFCorrection");
455 m_low_TOF_cut = getProperty("LowTOFCut");
456 m_high_TOF_cut = getProperty("HighTOFCut");
457
458 // Read in default beam center
459 m_center_x = getProperty("BeamCenterX");
460 m_center_y = getProperty("BeamCenterY");
461 const bool noBeamCenter = getProperty("NoBeamCenter");
462
463 // Reduction property manager
464 const std::string reductionManagerName = getProperty("ReductionProperties");
465 std::shared_ptr<PropertyManager> reductionManager;
466 if (PropertyManagerDataService::Instance().doesExist(reductionManagerName)) {
467 reductionManager = PropertyManagerDataService::Instance().retrieve(reductionManagerName);
468 } else {
469 reductionManager = std::make_shared<PropertyManager>();
470 PropertyManagerDataService::Instance().addOrReplace(reductionManagerName, reductionManager);
471 }
472
473 if (!reductionManager->existsProperty("LoadAlgorithm")) {
474 auto loadProp = std::make_unique<AlgorithmProperty>("LoadAlgorithm");
475 setPropertyValue("InputWorkspace", "");
476 setProperty("NoBeamCenter", false);
477 loadProp->setValue(toString());
478 reductionManager->declareProperty(std::move(loadProp));
479 }
480
481 if (!reductionManager->existsProperty("InstrumentName"))
482 reductionManager->declareProperty(std::make_unique<PropertyWithValue<std::string>>("InstrumentName", "EQSANS"));
483
484 // Output log
485 m_output_message = "";
486
487 // Check whether we need to load the data
488 if (!inputEventWS) {
489 const bool loadMonitors = getProperty("LoadMonitors");
490 const bool loadNexusInstrumentXML = getProperty("LoadNexusInstrumentXML");
491 auto loadAlg = createChildAlgorithm("LoadEventNexus", 0, 0.2);
492 loadAlg->setProperty("LoadMonitors", loadMonitors);
493 loadAlg->setProperty("Filename", fileName);
494 loadAlg->setProperty("LoadNexusInstrumentXML", loadNexusInstrumentXML);
495 if (skipTOFCorrection) {
496 if (m_low_TOF_cut > 0.0)
497 loadAlg->setProperty("FilterByTofMin", m_low_TOF_cut);
498 if (m_high_TOF_cut > 0.0)
499 loadAlg->setProperty("FilterByTofMax", m_high_TOF_cut);
500 }
501 loadAlg->execute();
502 Workspace_sptr dataWS_asWks = loadAlg->getProperty("OutputWorkspace");
503 dataWS = std::dynamic_pointer_cast<MatrixWorkspace>(dataWS_asWks);
504
505 // Get monitor workspace as necessary
506 std::string mon_wsname = getPropertyValue("OutputWorkspace") + "_monitors";
507 if (loadMonitors && loadAlg->existsProperty("MonitorWorkspace")) {
508 Workspace_sptr monWSOutput = loadAlg->getProperty("MonitorWorkspace");
509 MatrixWorkspace_sptr monWS = std::dynamic_pointer_cast<MatrixWorkspace>(monWSOutput);
510 if ((monWSOutput) && (!monWS))
511 // this was a group workspace - EQSansLoad does not support multi period
512 // data yet
513 throw Exception::NotImplementedError("The file contains multi period "
514 "data, support for this is not "
515 "implemented in EQSANSLoad yet");
516 declareProperty(std::make_unique<WorkspaceProperty<>>("MonitorWorkspace", mon_wsname, Direction::Output),
517 "Monitors from the Event NeXus file");
518 setProperty("MonitorWorkspace", monWS);
519 }
520 } else {
521 MatrixWorkspace_sptr outputWS = getProperty("OutputWorkspace");
522 EventWorkspace_sptr outputEventWS = std::dynamic_pointer_cast<EventWorkspace>(outputWS);
523 if (inputEventWS != outputEventWS) {
524 auto copyAlg = createChildAlgorithm("CloneWorkspace", 0, 0.2);
525 copyAlg->setProperty("InputWorkspace", inputEventWS);
526 copyAlg->executeAsChildAlg();
527 Workspace_sptr dataWS_asWks = copyAlg->getProperty("OutputWorkspace");
528 dataWS = std::dynamic_pointer_cast<MatrixWorkspace>(dataWS_asWks);
529 } else {
530 dataWS = std::dynamic_pointer_cast<MatrixWorkspace>(inputEventWS);
531 }
532 }
533
534 // Get the sample flange-to-detector distance
535 // We have to call it "SampleDetectorDistance" in the workspace
536 double sfdd = 0.0;
537 double s2d = 0.0;
538 const double sampleflange_det_dist = getProperty("SampleDetectorDistance");
539 if (!isEmpty(sampleflange_det_dist)) {
540 sfdd = sampleflange_det_dist;
541 } else {
542 if (!dataWS->run().hasProperty("detectorZ")) {
543 g_log.error() << "Could not determine Z position: the "
544 "SampleDetectorDistance property was not set "
545 "and the run logs do not contain the detectorZ property\n";
546 throw std::invalid_argument("Could not determine Z position: stopping execution");
547 }
548
549 const std::string dzName = "detectorZ";
550 Mantid::Kernel::Property *prop = dataWS->run().getProperty(dzName);
551 auto *dp = dynamic_cast<Mantid::Kernel::TimeSeriesProperty<double> *>(prop);
552 if (!dp)
553 throw std::runtime_error("Could not cast (interpret) the property " + dzName +
554 " as a time series property value.");
555 sfdd = dp->getStatistics().mean;
556
557 // Modify SDD according to the DetectorDistance offset if given
558 const double sampleflange_det_offset = getProperty("DetectorOffset");
559 if (!isEmpty(sampleflange_det_offset))
560 sfdd += sampleflange_det_offset;
561
562 // Modify SDD according to SampleDetectorDistanceOffset offset if given.
563 // This is here for backward compatibility.
564 const double sample_det_offset = getProperty("SampleDetectorDistanceOffset");
565 if (!isEmpty(sample_det_offset))
566 sfdd += sample_det_offset;
567 if (!isEmpty(sample_det_offset) && !isEmpty(sampleflange_det_offset))
568 g_log.error() << "Both DetectorOffset and SampleDetectorDistanceOffset "
569 "are set. Only one should be used.\n";
570
571 s2d = sfdd;
572 // Modify S2D according to the SampleDistance offset if given
573 // This assumes that a positive offset moves the sample toward the detector
574 const double sampleflange_sample_offset = getProperty("SampleOffset");
575 if (!isEmpty(sampleflange_sample_offset)) {
576 s2d -= sampleflange_sample_offset;
577
578 // Move the sample to its correct position
579 auto mvAlg = createChildAlgorithm("MoveInstrumentComponent", 0.2, 0.4);
580 mvAlg->setProperty<MatrixWorkspace_sptr>("Workspace", dataWS);
581 mvAlg->setProperty("ComponentName", "sample-position");
582 mvAlg->setProperty("Z", sampleflange_sample_offset / 1000.0);
583 mvAlg->setProperty("RelativePosition", false);
584 mvAlg->executeAsChildAlg();
585 g_log.information() << "Moving sample to " << sampleflange_sample_offset / 1000.0 << " meters\n";
587 " Sample position: " + Poco::NumberFormatter::format(sampleflange_sample_offset / 1000.0, 3) + " m\n";
588 }
589 }
590 dataWS->mutableRun().addProperty("sampleflange_detector_distance", sfdd, "mm", true);
591 dataWS->mutableRun().addProperty("sample_detector_distance", s2d, "mm", true);
592
593 // Move the detector to its correct position
594 auto mvAlg = createChildAlgorithm("MoveInstrumentComponent", 0.2, 0.4);
595 mvAlg->setProperty<MatrixWorkspace_sptr>("Workspace", dataWS);
596 mvAlg->setProperty("ComponentName", "detector1");
597 mvAlg->setProperty("Z", sfdd / 1000.0);
598 mvAlg->setProperty("RelativePosition", false);
599 mvAlg->executeAsChildAlg();
600 g_log.information() << "Moving detector to " << sfdd / 1000.0 << " meters\n";
601 m_output_message += " Detector position: " + Poco::NumberFormatter::format(sfdd / 1000.0, 3) + " m\n";
602
603 // Get the run number so we can find the proper config file
604 int run_number = 0;
605 std::string config_file;
606 if (dataWS->run().hasProperty("run_number")) {
607 const std::string run_str = dataWS->run().getPropertyValueAsType<std::string>("run_number");
608 Poco::NumberParser::tryParse(run_str, run_number);
609 // Find a proper config file
610 config_file = findConfigFile(run_number);
611 } else {
612 g_log.error() << "Could not find run number for workspace " << getPropertyValue("OutputWorkspace") << '\n';
613 m_output_message += " Could not find run number for data file\n";
614 }
615
616 // Process the config file
617 bool use_config = getProperty("UseConfig");
618 if (use_config && !config_file.empty()) {
619 // Special case to force reading the beam center from the config file
620 // We're adding this to be compatible with the original EQSANS load
621 // written in python
622 if (m_center_x == 0.0 && m_center_y == 0.0)
623 setProperty("UseConfigBeam", true);
624
625 readConfigFile(config_file);
626 } else if (use_config) {
627 use_config = false;
628 g_log.error() << "Cound not find config file for workspace " << getPropertyValue("OutputWorkspace") << '\n';
630 " Could not find configuration file for run " + Poco::NumberFormatter::format(run_number) + "\n";
631 }
632
633 // If we use the config file, move the moderator position
634 if (use_config) {
635 if (m_moderator_position > -13.0)
636 g_log.error() << "Moderator position seems close to the sample, please check\n";
637 g_log.information() << "Moving moderator to " << m_moderator_position << '\n';
638 m_output_message += " Moderator position: " + Poco::NumberFormatter::format(m_moderator_position) + " m\n";
639 mvAlg = createChildAlgorithm("MoveInstrumentComponent", 0.4, 0.45);
640 mvAlg->setProperty<MatrixWorkspace_sptr>("Workspace", dataWS);
641 mvAlg->setProperty("ComponentName", "moderator");
642 mvAlg->setProperty("Z", m_moderator_position);
643 mvAlg->setProperty("RelativePosition", false);
644 mvAlg->executeAsChildAlg();
645 }
646
647 // Get source aperture radius
649
650 // Move the beam center to its proper position
651 if (!noBeamCenter) {
653 if (reductionManager->existsProperty("LatestBeamCenterX") &&
654 reductionManager->existsProperty("LatestBeamCenterY")) {
655 m_center_x = reductionManager->getProperty("LatestBeamCenterX");
656 m_center_y = reductionManager->getProperty("LatestBeamCenterY");
657 }
658 }
660
661 // Add beam center to reduction properties, as the last beam center position
662 // that was used.
663 // This will give us our default position next time.
664 if (!reductionManager->existsProperty("LatestBeamCenterX"))
665 reductionManager->declareProperty(std::make_unique<PropertyWithValue<double>>("LatestBeamCenterX", m_center_x));
666 else
667 reductionManager->setProperty("LatestBeamCenterX", m_center_x);
668 if (!reductionManager->existsProperty("LatestBeamCenterY"))
669 reductionManager->declareProperty(std::make_unique<PropertyWithValue<double>>("LatestBeamCenterY", m_center_y));
670 else
671 reductionManager->setProperty("LatestBeamCenterY", m_center_y);
672 }
673
674 // Modify TOF
675 const bool correct_for_flight_path = getProperty("CorrectForFlightPath");
676 double wl_min = 0.0;
677 double wl_max;
678 double wl_combined_max = 0.0;
679 if (skipTOFCorrection) {
680 m_output_message += " Skipping EQSANS TOF correction: assuming a single frame\n";
681 dataWS->mutableRun().addProperty("is_frame_skipping", 0, true);
682 if (correct_for_flight_path) {
683 g_log.error() << "CorrectForFlightPath and SkipTOFCorrection can't be "
684 "set to true at the same time\n";
685 m_output_message += " Skipped flight path correction: see error log\n";
686 }
687 } else {
688 m_output_message += " Flight path correction ";
689 if (!correct_for_flight_path)
690 m_output_message += "NOT ";
691 m_output_message += "applied\n";
692 DataObjects::EventWorkspace_sptr dataWS_evt = std::dynamic_pointer_cast<EventWorkspace>(dataWS);
693 auto tofAlg = createChildAlgorithm("EQSANSTofStructure", 0.5, 0.7);
694 tofAlg->setProperty<EventWorkspace_sptr>("InputWorkspace", dataWS_evt);
695 tofAlg->setProperty("LowTOFCut", m_low_TOF_cut);
696 tofAlg->setProperty("HighTOFCut", m_high_TOF_cut);
697 tofAlg->setProperty("FlightPathCorrection", correct_for_flight_path);
698 tofAlg->executeAsChildAlg();
699 wl_min = tofAlg->getProperty("WavelengthMin");
700 wl_max = tofAlg->getProperty("WavelengthMax");
701 if (wl_min >= wl_max) {
702 g_log.error() << "Bad wavelength range\n";
703 g_log.error() << m_output_message << '\n';
704 }
705
706 const bool frame_skipping = tofAlg->getProperty("FrameSkipping");
707 dataWS->mutableRun().addProperty("wavelength_min", wl_min, "Angstrom", true);
708 dataWS->mutableRun().addProperty("wavelength_max", wl_max, "Angstrom", true);
709 dataWS->mutableRun().addProperty("is_frame_skipping", int(frame_skipping), true);
710 wl_combined_max = wl_max;
712 " Wavelength range: " + Poco::NumberFormatter::format(wl_min) + " - " + Poco::NumberFormatter::format(wl_max);
713 if (frame_skipping) {
714 const double wl_min2 = tofAlg->getProperty("WavelengthMinFrame2");
715 const double wl_max2 = tofAlg->getProperty("WavelengthMaxFrame2");
716 wl_combined_max = wl_max2;
717 dataWS->mutableRun().addProperty("wavelength_min_frame2", wl_min2, "Angstrom", true);
718 dataWS->mutableRun().addProperty("wavelength_max_frame2", wl_max2, "Angstrom", true);
719 m_output_message += " and " + Poco::NumberFormatter::format(wl_min2) + " - " +
720 Poco::NumberFormatter::format(wl_max2) + " Angstrom\n";
721 } else
722 m_output_message += " Angstrom\n";
723 }
724
725 // Convert to wavelength
726 // Checked on 8/10/17 - changed from "sdd" to "sfdd" as was done above
727 // sfdd + ssd gives total distance (corrected by offset) from the source
728 const double ssd = fabs(dataWS->getInstrument()->getSource()->getPos().Z()) * 1000.0;
729 const double conversion_factor = 3.9560346 / (sfdd + ssd);
731 " TOF to wavelength conversion factor: " + Poco::NumberFormatter::format(conversion_factor) + "\n";
732
733 if (skipTOFCorrection) {
734 DataObjects::EventWorkspace_sptr dataWS_evt = std::dynamic_pointer_cast<EventWorkspace>(dataWS);
735 if (dataWS_evt->getNumberEvents() == 0)
736 throw std::invalid_argument("No event to process: check your TOF cuts");
737 wl_min = dataWS_evt->getTofMin() * conversion_factor;
738 wl_max = dataWS_evt->getTofMax() * conversion_factor;
739 wl_combined_max = wl_max;
740 g_log.information() << "Wavelength range: " << wl_min << " to " << wl_max << '\n';
741 dataWS->mutableRun().addProperty("wavelength_min", wl_min, "Angstrom", true);
742 dataWS->mutableRun().addProperty("wavelength_max", wl_max, "Angstrom", true);
743 }
744
745 auto scAlg = createChildAlgorithm("ScaleX", 0.7, 0.71);
746 scAlg->setProperty<MatrixWorkspace_sptr>("InputWorkspace", dataWS);
747 scAlg->setProperty<MatrixWorkspace_sptr>("OutputWorkspace", dataWS);
748 scAlg->setProperty("Factor", conversion_factor);
749 scAlg->executeAsChildAlg();
750 dataWS->getAxis(0)->setUnit("Wavelength");
751
752 // Rebin so all the wavelength bins are aligned
753 const bool preserveEvents = getProperty("PreserveEvents");
754 const double wl_step = getProperty("WavelengthStep");
755
756 const double wl_min_rounded = round(wl_min * 100.0) / 100.0;
757 const double wl_max_rounded = round(wl_combined_max * 100.0) / 100.0;
758 std::string params = Poco::NumberFormatter::format(wl_min_rounded, 2) + "," + Poco::NumberFormatter::format(wl_step) +
759 "," + Poco::NumberFormatter::format(wl_max_rounded, 2);
760 g_log.information() << "Rebin parameters: " << params << '\n';
761 auto rebinAlg = createChildAlgorithm("Rebin", 0.71, 0.72);
762 rebinAlg->setProperty<MatrixWorkspace_sptr>("InputWorkspace", dataWS);
763 if (preserveEvents)
764 rebinAlg->setProperty<MatrixWorkspace_sptr>("OutputWorkspace", dataWS);
765 rebinAlg->setPropertyValue("Params", params);
766 rebinAlg->setProperty("PreserveEvents", preserveEvents);
767 rebinAlg->executeAsChildAlg();
768
769 if (!preserveEvents)
770 dataWS = rebinAlg->getProperty("OutputWorkspace");
771
772 dataWS->mutableRun().addProperty("event_ws", getPropertyValue("OutputWorkspace"), true);
773 setProperty<MatrixWorkspace_sptr>("OutputWorkspace", std::dynamic_pointer_cast<MatrixWorkspace>(dataWS));
774 setPropertyValue("OutputMessage", m_output_message);
775}
776
777} // namespace Mantid::WorkflowAlgorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
#define fabs(x)
Definition: Matrix.cpp:22
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
Kernel::Logger & g_log
Definition: Algorithm.h:451
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
A property class for workspaces.
Records the filename and the description of failure.
Definition: Exception.h:98
Marks code as not implemented yet.
Definition: Exception.h:138
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
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
The concrete, templated class for properties.
Base class for properties.
Definition: Property.h:94
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
A specialised Property class for holding a series of time-value pairs.
TimeSeriesPropertyStatistics getStatistics() const
Return a TimeSeriesPropertyStatistics object.
Class for 3D vectors.
Definition: V3D.h:34
constexpr double X() const noexcept
Get x.
Definition: V3D.h:232
constexpr double Y() const noexcept
Get y.
Definition: V3D.h:233
void readBeamCenter(const std::string &line)
Read the beam center from a config file string.
Definition: EQSANSLoad.cpp:196
void readTOFcuts(const std::string &line)
Read the TOF cuts from a config file string.
Definition: EQSANSLoad.cpp:178
void exec() override
Execution code.
Definition: EQSANSLoad.cpp:434
API::MatrixWorkspace_sptr dataWS
Definition: EQSANSLoad.h:70
void readModeratorPosition(const std::string &line)
Read the moderator position from a config file string.
Definition: EQSANSLoad.cpp:214
void init() override
Initialisation code.
Definition: EQSANSLoad.cpp:39
void getSourceSlitSize()
Get the source slit size from the slit information of the run properties.
Definition: EQSANSLoad.cpp:264
void readSourceSlitSize(const std::string &line)
Read the source slit sizes from a config file string.
Definition: EQSANSLoad.cpp:231
void moveToBeamCenter()
Move the detector according to the beam center.
Definition: EQSANSLoad.cpp:341
void readRectangularMasks(const std::string &line)
Read rectangular masks from a config file string.
Definition: EQSANSLoad.cpp:155
void readConfigFile(const std::string &filePath)
Read a config file.
Definition: EQSANSLoad.cpp:388
std::string findConfigFile(const int &run)
Find the most appropriate configuration file for a given run.
Definition: EQSANSLoad.cpp:120
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
bool exists(::NeXus::File &file, const std::string &name)
Based on the current group in the file, does the named sub-entry exist?
std::shared_ptr< EventWorkspace > EventWorkspace_sptr
shared pointer to the EventWorkspace class
void getDefaultBeamCenter(const API::MatrixWorkspace_sptr &dataWS, double &pixel_x, double &pixel_y)
int getDetectorFromPixel(const int &pixel_x, const int &pixel_y, const API::MatrixWorkspace_sptr &dataWS)
void getCoordinateFromPixel(const double &pixel_x, const double &pixel_y, const API::MatrixWorkspace_sptr &dataWS, double &x, double &y)
double getRunPropertyDbl(const MatrixWorkspace_sptr &inputWS, const std::string &pname)
Returns the value of a run property from a given workspace.
Definition: EQSANSLoad.cpp:109
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