Mantid
Loading...
Searching...
No Matches
LoadDNSSCD.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"
33#include <Poco/DateTime.h>
34#include <Poco/DateTimeFormat.h>
35#include <Poco/DateTimeFormatter.h>
36#include <Poco/DateTimeParser.h>
37#include <Poco/DirectoryIterator.h>
38#include <Poco/File.h>
39#include <Poco/Path.h>
40#include <boost/date_time/posix_time/posix_time.hpp>
41#include <boost/exception/diagnostic_information.hpp>
42#include <boost/exception_ptr.hpp>
43#include <boost/regex.hpp>
44
45#include <algorithm>
46#include <iomanip>
47#include <iterator>
48#include <map>
49
50//========================
51// helper functions
52namespace {
53void eraseSubStr(std::string &str, const std::string &toErase) {
54 // Search for the substring in string
55 size_t pos = str.find(toErase);
56 if (pos != std::string::npos) {
57 // If found then erase it from string
58 str.erase(pos, toErase.length());
59 }
60}
61
62std::string parseTime(std::string &str) {
63 // remove unnecessary symbols
64 eraseSubStr(str, "#");
65 eraseSubStr(str, "start");
66 eraseSubStr(str, "stopped");
67 eraseSubStr(str, "at");
68 auto it =
69 std::find_if(str.begin(), str.end(), [](char ch) { return !std::isspace<char>(ch, std::locale::classic()); });
70 str.erase(str.begin(), it);
71 using namespace boost::posix_time;
72 // try to parse as a posix time
73 try {
74 auto time = time_from_string(str);
75 return to_iso_extended_string(time);
76 } catch (std::exception &) {
77 int tzd;
78 Poco::DateTime dt;
79 bool ok = Poco::DateTimeParser::tryParse(str, dt, tzd);
80 if (ok) {
81 auto time = Poco::DateTimeFormatter::format(dt, "%Y-%m-%dT%H:%M:%S");
82 return time;
83 }
84 std::string result("");
85 return result;
86 }
87}
88
89} // anonymous namespace
90//============================
91
92using namespace Mantid::Kernel;
93using namespace Mantid::API;
94using namespace Mantid::DataObjects;
95using namespace Mantid::Geometry;
96
97namespace Mantid::MDAlgorithms {
98
100
101//----------------------------------------------------------------------------------------------
104LoadDNSSCD::LoadDNSSCD() : m_columnSep("\t, ;"), m_nDims(4), m_tof_max(20000.0) {}
105
113 // DNS data acquisition writes ascii files with .d_dat extension
114 int confidence(0);
115 if ((descriptor.extension() == ".d_dat") && descriptor.isAscii()) {
116 confidence = 80;
117 }
118 return confidence;
119}
120
121//----------------------------------------------------------------------------------------------
125 std::vector<std::string> exts(1, ".d_dat");
126 declareProperty(std::make_unique<MultipleFileProperty>("Filenames", exts),
127 "Select one or more DNS SCD .d_dat files to load."
128 "Files must be measured at the same conditions.");
129
130 declareProperty(std::make_unique<WorkspaceProperty<IMDEventWorkspace>>("OutputWorkspace", "", Direction::Output),
131 "An output MDEventWorkspace.");
132
134 std::make_unique<WorkspaceProperty<IMDEventWorkspace>>("NormalizationWorkspace", "", Direction::Output),
135 "An output normalization MDEventWorkspace.");
136
137 const std::vector<std::string> normOptions = {"monitor", "time"};
138 declareProperty("Normalization", "monitor", std::make_shared<StringListValidator>(normOptions),
139 "Algorithm will create a separate normalization workspace. "
140 "Choose whether it should contain monitor counts or time.");
141
142 const std::vector<std::string> wsOptions = {"raw", "HKL"};
143 declareProperty("LoadAs", "HKL", std::make_shared<StringListValidator>(wsOptions),
144 "Choose whether the algorithm should load raw data"
145 "or convert to H,K,L,dE space");
146
147 auto mustBePositive = std::make_shared<BoundedValidator<double>>();
148 mustBePositive->setLower(0.0);
149 auto reasonableAngle = std::make_shared<BoundedValidator<double>>();
150 reasonableAngle->setLower(5.0);
151 reasonableAngle->setUpper(175.0);
152 // clang-format off
153 auto mustBe3D = std::make_shared<ArrayLengthValidator<double> >(3);
154 auto mustBe2D = std::make_shared<ArrayLengthValidator<double> >(2);
155 // clang-format on
156 std::vector<double> u0(3, 0), v0(3, 0);
157 u0[0] = 1.;
158 u0[1] = 1.;
159 v0[2] = 1.;
160
161 declareProperty(std::make_unique<PropertyWithValue<double>>("a", 1.0, mustBePositive->clone(), Direction::Input),
162 "Lattice parameter a in Angstrom");
163 declareProperty(std::make_unique<PropertyWithValue<double>>("b", 1.0, mustBePositive->clone(), Direction::Input),
164 "Lattice parameter b in Angstrom");
165 declareProperty(std::make_unique<PropertyWithValue<double>>("c", 1.0, std::move(mustBePositive), Direction::Input),
166 "Lattice parameter c in Angstrom");
168 std::make_unique<PropertyWithValue<double>>("alpha", 90.0, reasonableAngle->clone(), Direction::Input),
169 "Angle between b and c in degrees");
170 declareProperty(std::make_unique<PropertyWithValue<double>>("beta", 90.0, reasonableAngle->clone(), Direction::Input),
171 "Angle between a and c in degrees");
173 std::make_unique<PropertyWithValue<double>>("gamma", 90.0, std::move(reasonableAngle), Direction::Input),
174 "Angle between a and b in degrees");
175
177 "OmegaOffset", 0.0, std::make_shared<BoundedValidator<double>>(), Direction::Input),
178 "Angle in degrees between (HKL1) and the beam axis"
179 "if the goniometer is at zero.");
180 declareProperty(std::make_unique<ArrayProperty<double>>("HKL1", std::move(u0), mustBe3D->clone()),
181 "Indices of the vector in reciprocal space in the horizontal plane at "
182 "angle Omegaoffset, "
183 "if the goniometer is at zero.");
184
185 declareProperty(std::make_unique<ArrayProperty<double>>("HKL2", std::move(v0), std::move(mustBe3D)),
186 "Indices of a second vector in reciprocal space in the horizontal plane "
187 "not parallel to HKL1");
188
189 std::vector<double> ttl(2, 0);
190 ttl[1] = 180.0;
191 declareProperty(std::make_unique<ArrayProperty<double>>("TwoThetaLimits", std::move(ttl), std::move(mustBe2D)),
192 "Range (min, max) of scattering angles (2theta, in degrees) to consider. "
193 "Everything out of this range will be cut.");
194
195 declareProperty(std::make_unique<WorkspaceProperty<API::ITableWorkspace>>("LoadHuberFrom", "", Direction::Input,
197 "A table workspace to load a list of raw sample rotation angles. "
198 "Huber angles given in the data files will be ignored.");
199
202 "A workspace name to save a list of raw sample rotation angles.");
203
204 auto mustBeIntPositive = std::make_shared<BoundedValidator<int>>();
205 mustBeIntPositive->setLower(0);
207 std::make_unique<PropertyWithValue<int>>("ElasticChannel", 0, std::move(mustBeIntPositive), Direction::Input),
208 "Elastic channel number. Only for TOF data.");
209
210 auto mustBeNegative = std::make_shared<BoundedValidator<double>>();
211 mustBeNegative->setUpper(0.0);
213 std::make_unique<PropertyWithValue<double>>("DeltaEmin", -10.0, std::move(mustBeNegative), Direction::Input),
214 "Minimal energy transfer to consider. Should be <=0. Only for TOF data.");
215}
216
217//----------------------------------------------------------------------------------------------
222 ColumnVector<double> huber = tws->getVector("Huber(degrees)");
223 // set huber[0] for each run in m_data
224 for (auto &ds : m_data) {
225 ds.huber = huber[0];
226 }
227 // dublicate runs for each huber in the table
228 std::vector<ExpData> old(m_data);
229 for (size_t i = 1; i < huber.size(); ++i) {
230 for (auto &ds : old) {
231 ds.huber = huber[i];
232 m_data.emplace_back(ds);
233 }
234 }
235}
236
237//----------------------------------------------------------------------------------------------
241 std::vector<double> huber;
242 huber.reserve(m_data.size());
243 std::transform(m_data.cbegin(), m_data.cend(), std::back_inserter(huber), [](const auto &ds) { return ds.huber; });
244 // remove dublicates
245 std::sort(huber.begin(), huber.end());
246 huber.erase(unique(huber.begin(), huber.end()), huber.end());
247
248 Mantid::API::ITableWorkspace_sptr huberWS = WorkspaceFactory::Instance().createTable("TableWorkspace");
249 huberWS->addColumn("double", "Huber(degrees)");
250 for (size_t i = 0; i < huber.size(); i++) {
251 huberWS->appendRow();
252 huberWS->cell<double>(i, 0) = huber[i];
253 }
254 return huberWS;
255}
256//----------------------------------------------------------------------------------------------
260 MultipleFileProperty *multiFileProp = dynamic_cast<MultipleFileProperty *>(getPointerToProperty("Filenames"));
261 if (!multiFileProp) {
262 throw std::logic_error("Filenames property must have MultipleFileProperty type.");
263 }
264 std::vector<std::string> filenames = VectorHelper::flattenVector(multiFileProp->operator()());
265 if (filenames.empty())
266 throw std::invalid_argument("Must specify at least one filename.");
267
268 // set type of normalization
269 std::string normtype = getProperty("Normalization");
270 if (normtype == "monitor") {
271 m_normtype = "Monitor";
272 m_normfactor = 1.0;
273 } else {
274 m_normtype = "Timer";
275 m_normfactor = 0.0; // error for time should be 0
276 }
277
278 g_log.notice() << "The normalization workspace will contain " << m_normtype << ".\n";
279
280 ExperimentInfo_sptr expinfo = std::make_shared<ExperimentInfo>();
281 API::Run &run = expinfo->mutableRun();
282 for (auto fname : filenames) {
283 std::map<std::string, std::string> str_metadata;
284 std::map<std::string, double> num_metadata;
285 try {
286 read_data(fname, str_metadata, num_metadata);
287 // if no stop_time, take file_save_time
288 std::string time(str_metadata["stop_time"]);
289 if (time.empty()) {
290 g_log.warning() << "stop_time is empty! File save time will be used instead." << std::endl;
291 time = str_metadata["file_save_time"];
292 }
293 updateProperties<std::string>(run, str_metadata, time);
294 updateProperties<double>(run, num_metadata, time);
295 } catch (...) {
296 g_log.warning() << "Failed to read file " << fname;
297 g_log.warning() << ". This file will be ignored. " << std::endl;
298 g_log.debug() << boost::current_exception_diagnostic_information() << std::endl;
299 }
300 }
301
302 if (m_data.empty())
303 throw std::runtime_error("No valid DNS files have been provided. Nothing to load.");
304
305 // merge data with different time channel number is not allowed
306 auto ch_n = m_data.front().nchannels;
307 bool same_channel_number =
308 std::all_of(m_data.cbegin(), m_data.cend(), [ch_n](const ExpData &d) { return (d.nchannels == ch_n); });
309 if (!same_channel_number)
310 throw std::runtime_error("Error: cannot merge data with different TOF channel numbers.");
311
312 std::string load_as = getProperty("LoadAs");
313 if (load_as == "raw")
314 m_nDims = 3;
315
317 m_OutWS->addExperimentInfo(expinfo);
318
319 // load huber angles from a table workspace if given
320 ITableWorkspace_sptr huberWS = getProperty("LoadHuberFrom");
321 if (huberWS) {
322 g_log.notice() << "Huber angles will be loaded from " << huberWS->getName() << std::endl;
323 loadHuber(huberWS);
324 }
325
326 // get wavelength
327 TimeSeriesProperty<double> *wlprop = dynamic_cast<TimeSeriesProperty<double> *>(expinfo->run().getProperty("Lambda"));
328 // assume, that lambda is in nm
329 double wavelength = wlprop->minValue() * 10.0; // needed to estimate extents => minValue
330 run.addProperty("wavelength", wavelength);
331 run.getProperty("wavelength")->setUnits("Angstrom");
332
333 if (load_as == "raw") {
334 fillOutputWorkspaceRaw(wavelength);
335 } else {
336 fillOutputWorkspace(wavelength);
337 }
338
339 std::string saveHuberTableWS = getProperty("SaveHuberTo");
340 if (!saveHuberTableWS.empty()) {
342 setProperty("SaveHuberTo", huber_table);
343 }
344 setProperty("OutputWorkspace", m_OutWS);
345}
346
347int LoadDNSSCD::splitIntoColumns(std::list<std::string> &columns, std::string &str) {
348 boost::split(columns, str, boost::is_any_of(m_columnSep), boost::token_compress_on);
349 return static_cast<int>(columns.size());
350}
351
352//----------------------------------------------------------------------------------------------
353
354template <class T>
355void LoadDNSSCD::updateProperties(API::Run &run, std::map<std::string, T> &metadata, std::string time) {
356 auto it = metadata.begin();
357 while (it != metadata.end()) {
358 TimeSeriesProperty<T> *timeSeries(nullptr);
359 std::string name(it->first);
360 std::string units;
361 // std::regex does not work for rhel7, thus boost
362 boost::regex reg("([-_a-zA-Z]+)\\[(.*)]");
363 boost::smatch match;
364 if (boost::regex_search(name, match, reg) && match.size() > 2) {
365 std::string new_name(match.str(1));
366 units.assign(match.str(2));
367 name = new_name;
368 }
369 if (run.hasProperty(name)) {
370 timeSeries = dynamic_cast<TimeSeriesProperty<T> *>(run.getLogData(name));
371 if (!timeSeries)
372 throw std::invalid_argument("Log '" + name + "' already exists but the values are a different type.");
373 } else {
374 timeSeries = new TimeSeriesProperty<T>(name);
375 if (!units.empty())
376 timeSeries->setUnits(units);
377 run.addProperty(timeSeries);
378 }
379 timeSeries->addValue(time, it->second);
380 ++it;
381 }
382}
383//----------------------------------------------------------------------------------------------
385void LoadDNSSCD::fillOutputWorkspace(double wavelength) {
386
387 // dimensions
388 std::vector<std::string> vec_ID(4);
389 vec_ID[0] = "H";
390 vec_ID[1] = "K";
391 vec_ID[2] = "L";
392 vec_ID[3] = "DeltaE";
393
394 std::vector<std::string> dimensionNames(4);
395 dimensionNames[0] = "H";
396 dimensionNames[1] = "K";
397 dimensionNames[2] = "L";
398 dimensionNames[3] = "DeltaE";
399
401
402 double a, b, c, alpha, beta, gamma;
403 a = getProperty("a");
404 b = getProperty("b");
405 c = getProperty("c");
406 alpha = getProperty("alpha");
407 beta = getProperty("beta");
408 gamma = getProperty("gamma");
409 std::vector<double> u = getProperty("HKL1");
410 std::vector<double> v = getProperty("HKL2");
411
412 // load empty DNS instrument to access L1 and L2
413 auto loadAlg = AlgorithmManager::Instance().create("LoadEmptyInstrument");
414 loadAlg->setChild(true);
415 loadAlg->setLogging(false);
416 loadAlg->initialize();
417 loadAlg->setProperty("InstrumentName", "DNS");
418 loadAlg->setProperty("OutputWorkspace", "__DNS_Inst");
419 loadAlg->execute();
420 MatrixWorkspace_sptr instWS = loadAlg->getProperty("OutputWorkspace");
421 const auto &instrument = instWS->getInstrument();
422 const auto &samplePosition = instrument->getSample()->getPos();
423 const auto &sourcePosition = instrument->getSource()->getPos();
424 const auto beamVector = samplePosition - sourcePosition;
425 const auto l1 = beamVector.norm();
426 // calculate tof1
427 auto velocity = PhysicalConstants::h / (PhysicalConstants::NeutronMass * wavelength * 1e-10); // m/s
428 auto tof1 = 1e+06 * l1 / velocity; // microseconds
429 g_log.debug() << "TOF1 = " << tof1 << std::endl;
430 // calculate incident energy
431 auto Ei = 0.5 * PhysicalConstants::NeutronMass * velocity * velocity / PhysicalConstants::meV;
432 g_log.debug() << "Ei = " << Ei << std::endl;
433
434 double dEmin = getProperty("DeltaEmin");
435 // estimate extents
436 double qmax = 4.0 * M_PI / wavelength;
437 std::vector<double> extentMins = {-qmax * a, -qmax * b, -qmax * c, dEmin};
438 std::vector<double> extentMaxs = {qmax * a, qmax * b, qmax * c, Ei};
439
440 // Get MDFrame of HKL type with RLU
441 auto unitFactory = makeMDUnitFactoryChain();
442 auto unit = unitFactory->create(Units::Symbol::RLU.ascii());
443 Mantid::Geometry::HKL frame(unit);
444
445 // add dimensions
446 for (size_t i = 0; i < m_nDims; ++i) {
447 std::string id = vec_ID[i];
448 std::string name = dimensionNames[i];
450 id, name, frame, static_cast<coord_t>(extentMins[i]), static_cast<coord_t>(extentMaxs[i]), 5)));
451 }
452
453 // Set coordinate system
454 m_OutWS->setCoordinateSystem(coordinateSystem);
455
456 // calculate RUB matrix
458 o = Mantid::Geometry::OrientedLattice(a, b, c, alpha, beta, gamma);
459 o.setUFromVectors(Mantid::Kernel::V3D(u[0], u[1], u[2]), Mantid::Kernel::V3D(v[0], v[1], v[2]));
460
461 double omega_offset = getProperty("OmegaOffset");
462 omega_offset *= -1.0 * deg2rad;
463 DblMatrix rotm(3, 3);
464 rotm[0][0] = std::cos(omega_offset);
465 rotm[0][1] = 0.0;
466 rotm[0][2] = std::sin(omega_offset);
467 rotm[1][0] = 0.0;
468 rotm[1][1] = 1.0;
469 rotm[1][2] = 0.0;
470 rotm[2][0] = -std::sin(omega_offset);
471 rotm[2][1] = 0.0;
472 rotm[2][2] = std::cos(omega_offset);
473
474 DblMatrix ub(o.getUB());
475 ub = rotm * ub;
476 o.setUB(ub);
477 DblMatrix ub_inv(ub);
478 // invert the UB matrix
479 ub_inv.Invert();
480
481 // Creates a new instance of the MDEventInserter to output workspace
482 MDEventWorkspace<MDEvent<4>, 4>::sptr mdws_mdevt_4 =
483 std::dynamic_pointer_cast<MDEventWorkspace<MDEvent<4>, 4>>(m_OutWS);
484 MDEventInserter<MDEventWorkspace<MDEvent<4>, 4>::sptr> inserter(mdws_mdevt_4);
485
486 // create a normalization workspace
487 IMDEventWorkspace_sptr normWS = m_OutWS->clone();
488
489 // Creates a new instance of the MDEventInserter to norm workspace
490 MDEventWorkspace<MDEvent<4>, 4>::sptr normws_mdevt_4 =
491 std::dynamic_pointer_cast<MDEventWorkspace<MDEvent<4>, 4>>(normWS);
492 MDEventInserter<MDEventWorkspace<MDEvent<4>, 4>::sptr> norm_inserter(normws_mdevt_4);
493
494 // scattering angle limits
495 std::vector<double> tth_limits = getProperty("TwoThetaLimits");
496 double theta_min = tth_limits[0] * deg2rad / 2.0;
497 double theta_max = tth_limits[1] * deg2rad / 2.0;
498
499 // get elastic channel from the user input
500 int echannel_user = getProperty("ElasticChannel");
501
502 // Go though each element of m_data to convert to MDEvent
503 for (ExpData ds : m_data) {
504 uint16_t expInfoIndex = 0;
505 signal_t norm_signal(ds.norm);
506 signal_t norm_error = std::sqrt(m_normfactor * norm_signal);
507 double ki = 2.0 * M_PI / ds.wavelength;
508 for (size_t i = 0; i < ds.detID.size(); i++) {
509 const auto &detector = instWS->getDetector(i);
510 const auto &detectorPosition = detector->getPos();
511 const auto detectorVector = detectorPosition - samplePosition;
512 const auto l2 = detectorVector.norm();
513 auto tof2_elastic = 1e+06 * l2 / velocity;
514 // geometric elastic channel
515 auto echannel_geom = static_cast<int>(std::ceil(tof2_elastic / ds.chwidth));
516 // rotate the signal array to get elastic peak at right position
517 int ch_diff = echannel_geom - echannel_user;
518 if ((echannel_user > 0) && (ch_diff < 0)) {
519 std::rotate(ds.signal[i].begin(), ds.signal[i].begin() - ch_diff, ds.signal[i].end());
520 } else if ((echannel_user > 0) && (ch_diff > 0)) {
521 std::rotate(ds.signal[i].rbegin(), ds.signal[i].rbegin() + ch_diff, ds.signal[i].rend());
522 }
523 detid_t detid(ds.detID[i]);
524 double theta = 0.5 * (ds.detID[i] * 5.0 - ds.deterota) * deg2rad;
525 auto nchannels = static_cast<int64_t>(ds.signal[i].size());
526 if ((theta > theta_min) && (theta < theta_max)) {
528 for (int64_t channel = 0; channel < nchannels; channel++) {
530 double signal = ds.signal[i][channel];
531 signal_t error = std::sqrt(signal);
532 double tof2 = static_cast<double>(channel) * ds.chwidth + 0.5 * ds.chwidth; // bin centers
533 double dE = 0.0;
534 if (nchannels > 1) {
535 double v2 = 1e+06 * l2 / tof2;
536 dE = Ei - 0.5 * PhysicalConstants::NeutronMass * v2 * v2 / PhysicalConstants::meV;
537 }
538 if (dE > dEmin) {
539 double kf = std::sqrt(ki * ki - 2.0e-20 * PhysicalConstants::NeutronMass * dE * PhysicalConstants::meV /
541 double tlab = std::atan2(ki - kf * cos(2.0 * theta), kf * sin(2.0 * theta));
542 double omega = (ds.huber - ds.deterota) * deg2rad - tlab;
543 V3D uphi(-cos(omega), 0, -sin(omega));
544 double qabs = 0.5 * std::sqrt(ki * ki + kf * kf - 2.0 * ki * kf * cos(2.0 * theta)) / M_PI;
545 V3D hphi = uphi * qabs; // qabs = ki * sin(theta), for elastic case;
546 V3D hkl = ub_inv * hphi;
547 std::vector<Mantid::coord_t> millerindex(4);
548 millerindex[0] = static_cast<float>(hkl.X());
549 millerindex[1] = static_cast<float>(hkl.Y());
550 millerindex[2] = static_cast<float>(hkl.Z());
551 millerindex[3] = static_cast<float>(dE);
552 PARALLEL_CRITICAL(addValues) {
553 inserter.insertMDEvent(static_cast<float>(signal), static_cast<float>(error * error),
554 static_cast<uint16_t>(expInfoIndex), 0, detid, millerindex.data());
555
556 norm_inserter.insertMDEvent(static_cast<float>(norm_signal), static_cast<float>(norm_error * norm_error),
557 static_cast<uint16_t>(expInfoIndex), 0, detid, millerindex.data());
558 }
559 }
561 }
563 }
564 }
565 }
566 setProperty("NormalizationWorkspace", normWS);
567}
568
569//----------------------------------------------------------------------------------------------
572void LoadDNSSCD::fillOutputWorkspaceRaw(double wavelength) {
573
574 // dimensions
575 std::vector<std::string> vec_ID(3);
576 vec_ID[0] = "Theta";
577 vec_ID[1] = "Omega";
578 vec_ID[2] = "TOF";
579
580 std::vector<std::string> dimensionNames(3);
581 dimensionNames[0] = "Scattering Angle";
582 dimensionNames[1] = "Omega";
583 dimensionNames[2] = "TOF";
584
586
587 // load empty DNS instrument to access L1 and L2
588 auto loadAlg = AlgorithmManager::Instance().create("LoadEmptyInstrument");
589 loadAlg->setChild(true);
590 loadAlg->setLogging(false);
591 loadAlg->initialize();
592 loadAlg->setProperty("InstrumentName", "DNS");
593 loadAlg->setProperty("OutputWorkspace", "__DNS_Inst");
594 loadAlg->execute();
595 MatrixWorkspace_sptr instWS = loadAlg->getProperty("OutputWorkspace");
596 const auto &instrument = instWS->getInstrument();
597 const auto &samplePosition = instrument->getSample()->getPos();
598 const auto &sourcePosition = instrument->getSource()->getPos();
599 const auto beamVector = samplePosition - sourcePosition;
600 const auto l1 = beamVector.norm();
601 // calculate tof1
602 auto velocity = PhysicalConstants::h / (PhysicalConstants::NeutronMass * wavelength * 1e-10); // m/s
603 auto tof1 = 1e+06 * l1 / velocity; // microseconds
604 g_log.debug() << "TOF1 = " << tof1 << std::endl;
605 // calculate incident energy
606 auto Ei = 0.5 * PhysicalConstants::NeutronMass * velocity * velocity / PhysicalConstants::meV;
607 g_log.debug() << "Ei = " << Ei << std::endl;
608
609 // estimate extents
610 // scattering angle limits
611 std::vector<double> tth_limits = getProperty("TwoThetaLimits");
612 double theta_min = tth_limits[0] / 2.0;
613 double theta_max = tth_limits[1] / 2.0;
614
615 std::vector<double> extentMins = {theta_min, 0.0, tof1};
616 std::vector<double> extentMaxs = {theta_max, 360.0, m_tof_max};
617
618 // Get MDFrame of HKL type with RLU
619 // auto unitFactory = makeMDUnitFactoryChain();
620 // auto unit = unitFactory->create(Units::Symbol::RLU.ascii());
621 // Mantid::Geometry::HKL frame(unit);
622 const Kernel::UnitLabel unitLabel("Degrees");
623 Mantid::Geometry::GeneralFrame frame("Scattering Angle", unitLabel);
624
625 // add dimensions
626 for (size_t i = 0; i < 3; ++i) {
627 std::string id = vec_ID[i];
628 std::string name = dimensionNames[i];
630 name, id, frame, static_cast<coord_t>(extentMins[i]), static_cast<coord_t>(extentMaxs[i]), 5)));
631 }
632
633 // Set coordinate system
634 m_OutWS->setCoordinateSystem(coordinateSystem);
635
636 // Creates a new instance of the MDEventInserter to output workspace
637 MDEventWorkspace<MDEvent<3>, 3>::sptr mdws_mdevt_3 =
638 std::dynamic_pointer_cast<MDEventWorkspace<MDEvent<3>, 3>>(m_OutWS);
639 MDEventInserter<MDEventWorkspace<MDEvent<3>, 3>::sptr> inserter(mdws_mdevt_3);
640
641 // create a normalization workspace
642 IMDEventWorkspace_sptr normWS = m_OutWS->clone();
643
644 // Creates a new instance of the MDEventInserter to norm workspace
645 MDEventWorkspace<MDEvent<3>, 3>::sptr normws_mdevt_3 =
646 std::dynamic_pointer_cast<MDEventWorkspace<MDEvent<3>, 3>>(normWS);
647 MDEventInserter<MDEventWorkspace<MDEvent<3>, 3>::sptr> norm_inserter(normws_mdevt_3);
648
649 // get elastic channel from the user input
650 int echannel_user = getProperty("ElasticChannel");
651
652 // Go though each element of m_data to convert to MDEvent
653 for (ExpData ds : m_data) {
654 uint16_t expInfoIndex = 0;
655 signal_t norm_signal(ds.norm);
656 signal_t norm_error = std::sqrt(m_normfactor * norm_signal);
657 for (size_t i = 0; i < ds.detID.size(); i++) {
658 const auto &detector = instWS->getDetector(i);
659 const auto &detectorPosition = detector->getPos();
660 const auto detectorVector = detectorPosition - samplePosition;
661 const auto l2 = detectorVector.norm();
662 auto tof2_elastic = 1e+06 * l2 / velocity;
663 // geometric elastic channel
664 auto echannel_geom = static_cast<int>(std::ceil(tof2_elastic / ds.chwidth));
665 // rotate the signal array to get elastic peak at right position
666 int ch_diff = echannel_geom - echannel_user;
667 if ((echannel_user > 0) && (ch_diff < 0)) {
668 std::rotate(ds.signal[i].begin(), ds.signal[i].begin() - ch_diff, ds.signal[i].end());
669 } else if ((echannel_user > 0) && (ch_diff > 0)) {
670 std::rotate(ds.signal[i].rbegin(), ds.signal[i].rbegin() + ch_diff, ds.signal[i].rend());
671 }
672
673 detid_t detid(ds.detID[i]);
674 double theta = 0.5 * (ds.detID[i] * 5.0 - ds.deterota);
675 auto nchannels = static_cast<int64_t>(ds.signal[i].size());
676 if ((theta > theta_min) && (theta < theta_max)) {
678 for (int64_t channel = 0; channel < nchannels; channel++) {
680 double signal = ds.signal[i][channel];
681 signal_t error = std::sqrt(signal);
682 double tof2(tof2_elastic);
683 if (nchannels > 1) {
684 tof2 = static_cast<double>(channel) * ds.chwidth + 0.5 * ds.chwidth; // bin centers
685 }
686 double omega = (ds.huber - ds.deterota);
687
688 std::vector<Mantid::coord_t> datapoint(3);
689 datapoint[0] = static_cast<float>(theta);
690 datapoint[1] = static_cast<float>(omega);
691 datapoint[2] = static_cast<float>(tof1 + tof2);
692 PARALLEL_CRITICAL(addValues) {
693 inserter.insertMDEvent(static_cast<float>(signal), static_cast<float>(error * error),
694 static_cast<uint16_t>(expInfoIndex), 0, detid, datapoint.data());
695
696 norm_inserter.insertMDEvent(static_cast<float>(norm_signal), static_cast<float>(norm_error * norm_error),
697 static_cast<uint16_t>(expInfoIndex), 0, detid, datapoint.data());
698 }
700 }
702 }
703 }
704 }
705 setProperty("NormalizationWorkspace", normWS);
706}
707
708void LoadDNSSCD::read_data(const std::string &fname, std::map<std::string, std::string> &str_metadata,
709 std::map<std::string, double> &num_metadata) {
710 std::ifstream file(fname);
711 std::string line;
712 std::string::size_type n;
713 std::string s;
714 boost::regex reg1("^#\\s+(\\w+):(.*)");
715 boost::regex reg2("^#\\s+((\\w+\\s)+)\\s+(-?\\d+(,\\d+)*(\\.\\d+(e\\d+)?)?)");
716 boost::smatch match;
717 getline(file, line);
718 n = line.find("DNS");
719 if (n == std::string::npos) {
720 throw std::invalid_argument("Not a DNS file");
721 }
722 // get file save time
723 Poco::File pfile(fname);
724 Poco::DateTime lastModified = pfile.getLastModified();
725 std::string wtime(Poco::DateTimeFormatter::format(lastModified, "%Y-%m-%dT%H:%M:%S"));
726 str_metadata.insert(std::make_pair("file_save_time", wtime));
727
728 // get file basename
729 Poco::Path p(fname);
730 str_metadata.insert(std::make_pair("run_number", p.getBaseName()));
731
732 // parse metadata
733 while (getline(file, line)) {
734 n = line.find("Lambda");
735 if (n != std::string::npos) {
736 boost::regex re("[\\s]+");
737 s = line.substr(5);
738 boost::sregex_token_iterator it(s.begin(), s.end(), re, -1);
739 boost::sregex_token_iterator reg_end;
740 getline(file, line);
741 std::string s2 = line.substr(2);
742 boost::sregex_token_iterator it2(s2.begin(), s2.end(), re, -1);
743 for (; (it != reg_end) && (it2 != reg_end); ++it) {
744 std::string token(it->str());
745 if (token.find_first_not_of(' ') == std::string::npos) {
746 ++it2;
747 continue;
748 }
749 if (token == "Mono") {
750 str_metadata.insert(std::make_pair(token, it2->str()));
751 } else {
752 num_metadata.insert(std::make_pair(token, std::stod(it2->str())));
753 }
754 ++it2;
755 }
756 }
757 // parse start and stop time
758 n = line.find("start");
759 if (n != std::string::npos) {
760 str_metadata.insert(std::make_pair("start_time", parseTime(line)));
761 getline(file, line);
762 str_metadata.insert(std::make_pair("stop_time", parseTime(line)));
763 getline(file, line);
764 }
765 if (boost::regex_search(line, match, reg1) && match.size() > 2) {
766 str_metadata.insert(std::make_pair(match.str(1), match.str(2)));
767 }
768 if (boost::regex_search(line, match, reg2) && match.size() > 2) {
769 s = match.str(1);
770 s.erase(std::find_if_not(s.rbegin(), s.rend(), ::isspace).base(), s.end());
771 num_metadata.insert(std::make_pair(s, std::stod(match.str(3))));
772 }
773 n = line.find("DATA");
774 if (n != std::string::npos) {
775 break;
776 }
777 }
778
779 std::map<std::string, double>::const_iterator m = num_metadata.lower_bound("TOF");
780 g_log.debug() << "TOF Channels number: " << m->second << std::endl;
781 std::map<std::string, double>::const_iterator w = num_metadata.lower_bound("Time");
782 g_log.debug() << "Channel width: " << w->second << std::endl;
783
784 ExpData ds;
785 ds.deterota = num_metadata["DeteRota"];
786 ds.huber = num_metadata["Huber"];
787 ds.wavelength = 10.0 * num_metadata["Lambda[nm]"];
788 ds.norm = num_metadata[m_normtype];
789 ds.chwidth = w->second;
790 ds.nchannels = static_cast<size_t>(std::ceil(m->second));
791
792 // read data array
793 getline(file, line);
794
795 std::list<std::string> columns;
796 while (getline(file, line)) {
797 boost::trim(line);
798 const int cols = splitIntoColumns(columns, line);
799 if (cols > 0) {
800 ds.detID.emplace_back(std::stoi(columns.front()));
801 columns.pop_front();
802 std::vector<double> signal;
803 std::transform(columns.begin(), columns.end(), std::back_inserter(signal),
804 [](const std::string &s) { return std::stod(s); });
805 ds.signal.emplace_back(signal);
806 }
807 }
808 // DNS PA detector bank has only 24 detectors
809 ds.detID.resize(24);
810 ds.signal.resize(24);
811 m_data.emplace_back(ds);
812}
813
814} // namespace Mantid::MDAlgorithms
double error
Definition: IndexPeaks.cpp:133
#define PARALLEL_START_INTERRUPT_REGION
Begins a block to skip processing is the algorithm has been interupted Note the end of the block if n...
Definition: MultiThreaded.h:94
#define PARALLEL_CRITICAL(name)
#define PARALLEL_END_INTERRUPT_REGION
Ends a block to skip processing is the algorithm has been interupted Note the start of the block if n...
#define PARALLEL_FOR_IF(condition)
Empty definitions - to enable set your complier to enable openMP.
#define PARALLEL_CHECK_INTERRUPT_REGION
Adds a check after a Parallel region to see if it was interupted.
#define DECLARE_FILELOADER_ALGORITHM(classname)
DECLARE_FILELOADER_ALGORITHM should be used in place of the standard DECLARE_ALGORITHM macro when wri...
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
Kernel::Property * getPointerToProperty(const std::string &name) const override
Get a property by name.
Definition: Algorithm.cpp:2033
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Definition: Algorithm.cpp:2076
ColumnVector gives access to the column elements without alowing its resizing.
size_t size()
Size of the vector.
bool hasProperty(const std::string &name) const
Does the property exist on the object.
Definition: LogManager.cpp:265
Kernel::Property * getLogData(const std::string &name) const
Access a single log entry.
Definition: LogManager.h:129
Kernel::Property * getProperty(const std::string &name) const
Returns the named property as a pointer.
Definition: LogManager.cpp:404
void addProperty(Kernel::Property *prop, bool overwrite=false)
Add data to the object in the form of a property.
Definition: LogManager.h:79
A property to allow a user to specify multiple files to load.
This class stores information regarding an experimental run as a series of log entries.
Definition: Run.h:38
A property class for workspaces.
static API::IMDEventWorkspace_sptr CreateMDWorkspace(size_t nd, const std::string &eventType="MDLeanEvent", const Mantid::API::MDNormalization &preferredNormalization=Mantid::API::MDNormalization::VolumeNormalization, const Mantid::API::MDNormalization &preferredNormalizationHisto=Mantid::API::MDNormalization::VolumeNormalization)
Create a MDEventWorkspace of the given type.
MDEventInserter : Helper class that provides a generic interface for adding events to an MDWorkspace ...
Templated class for the multi-dimensional event workspace.
GeneralFrame : Any MDFrame that isn't related to momemtum transfer.
Definition: GeneralFrame.h:21
HKL : HKL MDFrame.
Definition: HKL.h:21
Class to implement UB matrix.
void setUB(const Kernel::DblMatrix &newUB)
Sets the UB matrix and recalculates lattice parameters.
const Kernel::DblMatrix & setUFromVectors(const Kernel::V3D &u, const Kernel::V3D &v)
Create the U matrix from two vectors.
const Kernel::DblMatrix & getUB() const
Get the UB matrix.
Support for a property that holds an array of values.
Definition: ArrayProperty.h:28
BoundedValidator is a validator that requires the values to be between upper or lower bounds,...
Defines a wrapper around an open file.
static bool isAscii(const std::string &filename, const size_t nbytes=256)
Returns true if the file is considered ascii.
const std::string & extension() const
Access the file extension.
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 notice(const std::string &msg)
Logs at notice level.
Definition: Logger.cpp:95
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
T Invert()
LU inversion routine.
Definition: Matrix.cpp:924
The concrete, templated class for properties.
virtual void setUnits(const std::string &unit)
Sets the units of the property, as a string.
Definition: Property.cpp:186
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.
TYPE minValue() const
Returns the minimum value found in the series.
void addValue(const Types::Core::DateAndTime &time, const TYPE &value)
Add a value to the map using a DateAndTime object.
A base-class for the a class that is able to return unit labels in different representations.
Definition: UnitLabel.h:20
static const UnitLabel RLU
Reciprocal lattice units.
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
constexpr double Z() const noexcept
Get z.
Definition: V3D.h:234
LoadDNSSCD : Load a list of DNS .d_dat files into a MDEventWorkspace.
Definition: LoadDNSSCD.h:27
double m_normfactor
factor to multiply the error^2 for normalization
Definition: LoadDNSSCD.h:68
void exec() override
Run the algorithm.
Definition: LoadDNSSCD.cpp:259
size_t m_nDims
number of workspace dimensions
Definition: LoadDNSSCD.h:60
void loadHuber(const API::ITableWorkspace_sptr &tws)
Read Huber angles from a given table workspace.
Definition: LoadDNSSCD.cpp:221
std::string m_normtype
type of normalization;
Definition: LoadDNSSCD.h:66
void fillOutputWorkspaceRaw(double wavelength)
Fill output workspace with raw data theta, omega, tof space.
Definition: LoadDNSSCD.cpp:572
std::string m_columnSep
The column separator.
Definition: LoadDNSSCD.h:57
Mantid::API::IMDEventWorkspace_sptr m_OutWS
Output IMDEventWorkspace.
Definition: LoadDNSSCD.h:85
API::ITableWorkspace_sptr saveHuber()
Save Huber angles to a given table workspace.
Definition: LoadDNSSCD.cpp:240
void read_data(const std::string &fname, std::map< std::string, std::string > &str_metadata, std::map< std::string, double > &num_metadata)
Definition: LoadDNSSCD.cpp:708
void updateProperties(API::Run &run, std::map< std::string, T > &metadata, std::string time)
Definition: LoadDNSSCD.cpp:355
int splitIntoColumns(std::list< std::string > &columns, std::string &str)
Definition: LoadDNSSCD.cpp:347
void fillOutputWorkspace(double wavelength)
Fill output workspace with data converted to H, K, L, dE space.
Definition: LoadDNSSCD.cpp:385
void init() override
Initialise the properties.
Definition: LoadDNSSCD.cpp:124
std::vector< ExpData > m_data
Definition: LoadDNSSCD.h:82
int confidence(Kernel::FileDescriptor &descriptor) const override
Returns a confidence value that this algorithm can load a file.
Definition: LoadDNSSCD.cpp:112
const std::string name() const override
Algorithm's name for identification.
Definition: LoadDNSSCD.h:32
double m_tof_max
maximal TOF (for extends)
Definition: LoadDNSSCD.h:63
std::shared_ptr< IMDEventWorkspace > IMDEventWorkspace_sptr
Shared pointer to Mantid::API::IMDEventWorkspace.
std::shared_ptr< ITableWorkspace > ITableWorkspace_sptr
shared pointer to Mantid::API::ITableWorkspace
Kernel::Logger g_log("ExperimentInfo")
static logger object
std::shared_ptr< ExperimentInfo > ExperimentInfo_sptr
Shared pointer to ExperimentInfo.
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
constexpr double deg2rad
Defines units/enum for Crystal work.
Definition: AngleUnits.h:20
std::shared_ptr< MDHistoDimension > MDHistoDimension_sptr
Shared pointer to a MDHistoDimension.
std::vector< T > flattenVector(const std::vector< std::vector< T > > &v)
A convenience function to "flatten" the given vector of vectors into a single vector.
Definition: VectorHelper.h:69
SpecialCoordinateSystem
Special coordinate systems for Q3D.
MDUnitFactory_uptr MANTID_KERNEL_DLL makeMDUnitFactoryChain()
Convience method. Pre-constructed builder chain.
std::enable_if< std::is_pointer< Arg >::value, bool >::type threadSafe(Arg workspace)
Thread-safety check Checks the workspace to ensure it is suitable for multithreaded access.
Definition: MultiThreaded.h:22
static constexpr double NeutronMass
Mass of the neutron in kg.
static constexpr double h
Planck constant in J*s.
static constexpr double meV
1 meV in Joules.
static constexpr double h_bar
Planck constant in J*s, divided by 2 PI.
float coord_t
Typedef for the data type to use for coordinate axes in MD objects such as MDBox, MDEventWorkspace,...
Definition: MDTypes.h:27
int32_t detid_t
Typedef for a detector ID.
Definition: SpectrumInfo.h:21
double signal_t
Typedef for the signal recorded in a MDBox, etc.
Definition: MDTypes.h:36
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54
structure for experimental data
Definition: LoadDNSSCD.h:71
std::vector< std::vector< double > > signal
Definition: LoadDNSSCD.h:78