Mantid
Loading...
Searching...
No Matches
UpdateInstrumentFromFile.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 +
8#include "LoadRaw/isisraw2.h"
20
21#include <boost/algorithm/string/predicate.hpp>
22#include <boost/scoped_ptr.hpp>
23// clang-format off
24#include <nexus/NeXusFile.hpp>
25#include <nexus/NeXusException.hpp>
26// clang-format on
27
28#include <algorithm>
29#include <fstream>
30
31namespace Mantid::DataHandling {
32
33DECLARE_ALGORITHM(UpdateInstrumentFromFile)
34
35using namespace Kernel;
36using namespace API;
39using Kernel::V3D;
40
42UpdateInstrumentFromFile::UpdateInstrumentFromFile() : m_workspace(), m_ignorePhi(false), m_ignoreMonitors(true) {}
43
46 // When used as a Child Algorithm the workspace name is not used - hence the
47 // "Anonymous" to satisfy the validator
48 declareProperty(std::make_unique<WorkspaceProperty<MatrixWorkspace>>("Workspace", "Anonymous", Direction::InOut),
49 "The name of the workspace in which to store the imported instrument");
50 declareProperty(std::make_unique<FileProperty>("Filename", "", FileProperty::Load,
51 std::vector<std::string>{".raw", ".nxs", ".dat", ".s*"}),
52 "The filename of the input file.\n"
53 "Currently supports RAW, ISIS NeXus, DAT & multi-column (at "
54 "least 2) ascii files");
55 declareProperty("MoveMonitors", (!m_ignoreMonitors),
56 "If true the positions of any detectors marked as monitors "
57 "in the IDF will be moved also");
58 declareProperty("IgnorePhi", m_ignorePhi, "If true the phi values form the file will be ignored ");
59 declareProperty("AsciiHeader", "",
60 "If the file is a simple text file, then this property is used to"
61 "define the values in each column of the file. For example: "
62 "spectrum,theta,t0,-,R"
63 "Keywords=spectrum,ID,R,theta,phi. A dash means skip column. Keywords "
64 "are recognised"
65 "as identifying components to move to new positions. Any other names in "
66 "the list"
67 "are added as instrument parameters.");
68 declareProperty("SkipFirstNLines", 0,
69 "If the file is ASCII, then skip this "
70 "number of lines at the start of the "
71 "file");
72}
73
80 // Retrieve the filename from the properties
81 const std::string filename = getPropertyValue("Filename");
82 m_workspace = getProperty("Workspace");
83
84 if (!m_workspace->getInstrument()) {
85 throw std::runtime_error("Input workspace has no defined instrument");
86 }
87
88 m_ignorePhi = getProperty("IgnorePhi");
89 const bool moveMonitors = getProperty("MoveMonitors");
90 m_ignoreMonitors = (!moveMonitors);
91
92 // Check file type
93 if (NexusDescriptor::isReadable(filename)) {
94 LoadISISNexus2 isisNexus;
95 LoadEventNexus eventNexus;
96
97 // we open and close the HDF5 file.
98 // there is an issue with how HDF5 files are opened (only one at a time)
99 // swap the order of descriptors
100 boost::scoped_ptr<Kernel::NexusHDF5Descriptor> descriptorNexusHDF5(new Kernel::NexusHDF5Descriptor(filename));
101
102 boost::scoped_ptr<Kernel::NexusDescriptor> descriptor(new Kernel::NexusDescriptor(filename));
103
104 if (isisNexus.confidence(*descriptor) > 0 || eventNexus.confidence(*descriptorNexusHDF5) > 0) {
105 auto &nxFile = descriptor->data();
106 const auto &rootEntry = descriptor->firstEntryNameType();
107 nxFile.openGroup(rootEntry.first, rootEntry.second);
108 updateFromNeXus(nxFile);
109 return;
110 }
111 }
112
113 if (FileDescriptor::isAscii(filename)) {
114 // If no header specified & the extension is .dat or .sca, then assume ISIS
115 // DAT file structure
116 if (getPropertyValue("AsciiHeader").empty() &&
117 (boost::iends_with(filename, ".dat") || boost::iends_with(filename, ".sca"))) {
118 this->setPropertyValue("AsciiHeader", "ID,-,R,-,theta,phi,-,-,-,-,-,-,-,-,-,-,-,-,-");
119 this->setProperty("SkipFirstNLines", 2);
120 }
121 updateFromAscii(filename);
122 return;
123 }
124
125 LoadRawHelper isisRAW;
126 auto descriptor = std::make_unique<Kernel::FileDescriptor>(filename);
127 if (isisRAW.confidence(*descriptor) > 0) {
128 updateFromRaw(filename);
129 } else {
130 throw std::invalid_argument("File \"" + filename + "\" is not a valid input file.");
131 }
132}
133
138void UpdateInstrumentFromFile::updateFromRaw(const std::string &filename) {
139 ISISRAW2 iraw;
140 if (iraw.readFromFile(filename.c_str(), false) != 0) {
141 g_log.error("Unable to open file " + filename);
142 throw Exception::FileError("Unable to open File:", filename);
143 }
144
145 const int32_t numDetector = iraw.i_det;
146 std::vector<int32_t> detID(iraw.udet, iraw.udet + numDetector);
147 std::vector<float> l2(iraw.len2, iraw.len2 + numDetector);
148 std::vector<float> theta(iraw.tthe, iraw.tthe + numDetector);
149 // Is ut01 (=phi) present? Sometimes an array is present but has wrong values
150 // e.g.all 1.0 or all 2.0
151 bool phiPresent = iraw.i_use > 0 && iraw.ut[0] != 1.0 && iraw.ut[0] != 2.0;
152 std::vector<float> phi(0);
153 if (phiPresent) {
154 phi = std::vector<float>(iraw.ut, iraw.ut + numDetector);
155 } else {
156 phi = std::vector<float>(numDetector, 0.0);
157 }
158 g_log.information() << "Setting detector postions from RAW file.\n";
159 setDetectorPositions(detID, l2, theta, phi);
160}
161
166void UpdateInstrumentFromFile::updateFromNeXus(::NeXus::File &nxFile) {
167 try {
168 nxFile.openGroup("isis_vms_compat", "IXvms");
169 } catch (::NeXus::Exception &) {
170 throw std::runtime_error("Unknown NeXus flavour. Cannot update instrument "
171 "positions using this type of file");
172 }
173 // Det ID
174 std::vector<int32_t> detID;
175 nxFile.openData("UDET");
176 nxFile.getData(detID);
177 nxFile.closeData();
178 // Position information
179 std::vector<float> l2, theta, phi;
180 nxFile.openData("LEN2");
181 nxFile.getData(l2);
182 nxFile.closeData();
183 nxFile.openData("TTHE");
184 nxFile.getData(theta);
185 nxFile.closeData();
186 nxFile.openData("UT01");
187 nxFile.getData(phi);
188 nxFile.closeData();
189
190 g_log.information() << "Setting detector postions from NeXus file.\n";
191 setDetectorPositions(detID, l2, theta, phi);
192}
193
198void UpdateInstrumentFromFile::updateFromAscii(const std::string &filename) {
199 AsciiFileHeader header;
200 const bool isSpectrum = parseAsciiHeader(header);
201
202 // Throws for multiple detectors
203 const spec2index_map specToIndex(m_workspace->getSpectrumToWorkspaceIndexMap());
204
205 std::ifstream datfile(filename.c_str(), std::ios_base::in);
206 const int skipNLines = getProperty("SkipFirstNLines");
207 std::string line;
208 int lineCount(0);
209 while (lineCount < skipNLines) {
210 std::getline(datfile, line);
211 ++lineCount;
212 }
213
214 Geometry::ParameterMap &pmap = m_workspace->instrumentParameters();
215 auto &detectorInfo = m_workspace->mutableDetectorInfo();
216 const auto &spectrumInfo = m_workspace->spectrumInfo();
217
218 std::vector<double> colValues(header.colCount - 1, 0.0);
219 while (std::getline(datfile, line)) {
220 boost::trim(line);
221 std::istringstream is(line);
222 // Column 0 should be ID/spectrum number
223 int32_t detOrSpec(-1000);
224 is >> detOrSpec;
225 // If first thing read is not a number then skip the line
226 if (is.fail()) {
227 g_log.debug() << "Skipping \"" << line << "\". Cannot interpret as list of numbers.\n";
228 continue;
229 }
230
231 bool skip{false};
232 auto index = static_cast<size_t>(-1);
233 const Geometry::IDetector *det{nullptr};
234 if (isSpectrum) {
235 auto it = specToIndex.find(detOrSpec);
236 if (it != specToIndex.end()) {
237 index = it->second;
238 if (spectrumInfo.hasDetectors(index)) {
239 det = &spectrumInfo.detector(index);
240 } else {
241 skip = true;
242 }
243 } else {
244 g_log.debug() << "Skipping \"" << line << "\". Spectrum is not in workspace.\n";
245 continue;
246 }
247 } else {
248 try {
249 index = detectorInfo.indexOf(detOrSpec);
250 det = &detectorInfo.detector(index);
251 } catch (std::out_of_range &) {
252 skip = true;
253 }
254 }
255 if (skip || index == static_cast<size_t>(-1)) {
256 g_log.debug() << "Skipping \"" << line << "\". Spectrum in workspace but cannot find associated detector.\n";
257 continue;
258 }
259
260 std::vector<size_t> indices;
261 if (isSpectrum) {
262 if (auto group = dynamic_cast<const Geometry::DetectorGroup *>(det)) {
263 const auto detIDs = group->getDetectorIDs();
264 std::transform(detIDs.cbegin(), detIDs.cend(), std::back_inserter(indices),
265 [detectorInfo](const auto detID) { return detectorInfo.indexOf(detID); });
266 } else {
267 indices.emplace_back(detectorInfo.indexOf(det->getID()));
268 }
269 } else {
270 indices.emplace_back(index);
271 }
272
273 // Special cases for detector r,t,p. Everything else is
274 // attached as an detector parameter
275 double R(0.0), theta(0.0), phi(0.0);
276 for (size_t i = 1; i < header.colCount; ++i) {
277 double value(0.0);
278 is >> value;
279 if (i < header.colCount - 1 && is.eof()) {
280 // If stringstream is at EOF & we are not at the last column then
281 // there aren't enought columns in the file
282 throw std::runtime_error("UpdateInstrumentFromFile::updateFromAscii - "
283 "File contains fewer than expected number of "
284 "columns, check AsciiHeader property.");
285 }
286
287 if (i == header.rColIdx)
288 R = value;
289 else if (i == header.thetaColIdx)
290 theta = value;
291 else if (i == header.phiColIdx)
292 phi = value;
293 else if (header.detParCols.count(i) == 1) {
294 for (const auto detIndex : indices) {
295 auto id = detectorInfo.detector(detIndex).getComponentID();
296 pmap.addDouble(id, header.colToName[i], value);
297 }
298 }
299 }
300 // Check stream state. stringstream::EOF should have been reached, if not
301 // then there is still more to
302 // read and the file has more columns than the header indicated
303 if (!is.eof()) {
304 throw std::runtime_error("UpdateInstrumentFromFile::updateFromAscii - "
305 "File contains more than expected number of "
306 "columns, check AsciiHeader property.");
307 }
308
309 // If not supplied use current values
310 double r, t, p;
311 if (isSpectrum)
312 spectrumInfo.position(index).getSpherical(r, t, p);
313 else
314 detectorInfo.position(index).getSpherical(r, t, p);
315 if (header.rColIdx == 0)
316 R = r;
317 if (header.thetaColIdx == 0)
318 theta = t;
319 if (header.phiColIdx == 0 || m_ignorePhi)
320 phi = p;
321
322 for (const auto detIndex : indices)
323 setDetectorPosition(detectorInfo, detIndex, static_cast<float>(R), static_cast<float>(theta),
324 static_cast<float>(phi));
325 }
326}
327
336 const std::string header = getProperty("AsciiHeader");
337 if (header.empty()) {
338 throw std::invalid_argument("Ascii file provided but the AsciiHeader "
339 "property is empty, cannot interpret columns");
340 }
341
343 headerInfo.colCount = splitter.count();
344 auto it = splitter.begin(); // First column must be spectrum number or detector ID
345 const std::string &col0 = *it;
346 bool isSpectrum(false);
347 if (boost::iequals("spectrum", col0))
348 isSpectrum = true;
349 if (!isSpectrum && !boost::iequals("id", col0)) {
350 throw std::invalid_argument("Invalid AsciiHeader, first column name must "
351 "be either 'spectrum' or 'id'");
352 }
353
354 ++it;
355 size_t counter(1);
356 for (; it != splitter.end(); ++it) {
357 const std::string &colName = *it;
358 if (boost::iequals("R", colName))
359 headerInfo.rColIdx = counter;
360 else if (boost::iequals("theta", colName))
361 headerInfo.thetaColIdx = counter;
362 else if (boost::iequals("phi", colName))
363 headerInfo.phiColIdx = counter;
364 else if (boost::iequals("-", colName)) // Skip dashed
365 {
366 ++counter;
367 continue;
368 } else {
369 headerInfo.detParCols.insert(counter);
370 headerInfo.colToName.emplace(counter, colName);
371 }
372 ++counter;
373 }
374
375 return isSpectrum;
376}
377
385void UpdateInstrumentFromFile::setDetectorPositions(const std::vector<int32_t> &detID, const std::vector<float> &l2,
386 const std::vector<float> &theta, const std::vector<float> &phi) {
387 const auto numDetector = static_cast<int>(detID.size());
388 g_log.information() << "Setting new positions for " << numDetector << " detectors\n";
389
390 auto &detectorInfo = m_workspace->mutableDetectorInfo();
391 for (int i = 0; i < numDetector; ++i) {
392 try {
393 auto index = detectorInfo.indexOf(detID[i]);
394 double p{phi[i]};
395 if (m_ignorePhi) {
396 double r, t;
397 detectorInfo.position(index).getSpherical(r, t, p);
398 }
399 setDetectorPosition(detectorInfo, index, l2[i], theta[i], static_cast<float>(p));
400 } catch (std::out_of_range &) {
401 // Invalid detID[i]
402 continue;
403 }
404 progress(static_cast<double>(i) / numDetector, "Updating Detector Positions from File");
405 }
406}
407
417 const float l2, const float theta, const float phi) {
418 if (m_ignoreMonitors && detectorInfo.isMonitor(index))
419 return;
420
421 Kernel::V3D pos;
422 pos.spherical(l2, theta, phi);
423 detectorInfo.setPosition(index, pos);
424}
425
426} // namespace Mantid::DataHandling
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
double value
The value of the point.
Definition: FitMW.cpp:51
std::map< DeltaEMode::Type, std::string > index
Definition: DeltaEMode.cpp:19
isis raw file.
Definition: isisraw2.h:13
float * ut
nuse UT* user tables (total size NUSE*NDET) ut01=phi
Definition: isisraw.h:307
int i_use
number of user defined UTn tables NUSE
Definition: isisraw.h:298
int i_det
number of detectors NDET
Definition: isisraw.h:296
int readFromFile(const char *filename, bool read_data=true)
stuff
Definition: isisraw.cpp:923
int * udet
user detector number for each detector (size NDET)
Definition: isisraw.h:322
float * tthe
2theta scattering angle (size NDET)
Definition: isisraw.h:306
float * len2
L2 table (size NDET)
Definition: isisraw.h:304
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
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
@ Load
allowed here which will be passed to the algorithm
Definition: FileProperty.h:52
A property class for workspaces.
int confidence(Kernel::NexusHDF5Descriptor &descriptor) const override
Return the confidence with with this algorithm can load the file.
Loads a file in a NeXus format and stores it in a 2D workspace.
int confidence(Kernel::NexusDescriptor &descriptor) const override
Returns a confidence value that this algorithm can load a file.
Helper class for LoadRaw algorithms.
Definition: LoadRawHelper.h:40
int confidence(Kernel::FileDescriptor &descriptor) const override
Returns a confidence value that this algorithm can load a file.
void setDetectorPositions(const std::vector< int32_t > &detID, const std::vector< float > &l2, const std::vector< float > &theta, const std::vector< float > &phi)
Set the new detector positions.
void updateFromRaw(const std::string &filename)
Assumes the file is a raw file.
void updateFromAscii(const std::string &filename)
Updates from a more generic ascii file.
void setDetectorPosition(Geometry::DetectorInfo &detectorInfo, const size_t index, const float l2, const float theta, const float phi)
Set the new detector position for a single det ID.
API::MatrixWorkspace_sptr m_workspace
The input workspace to modify.
void updateFromNeXus(::NeXus::File &nxFile)
Assumes the file is an ISIS NeXus file.
bool parseAsciiHeader(AsciiFileHeader &headerInfo)
Parse the header and fill the headerInfo struct.
void init() override
Overwrites Algorithm method. Does nothing at present.
void exec() override
Overwrites Algorithm method.
Holds a collection of detectors.
Definition: DetectorGroup.h:28
Geometry::DetectorInfo is an intermediate step towards a DetectorInfo that is part of Instrument-2....
Definition: DetectorInfo.h:49
void setPosition(const size_t index, const Kernel::V3D &position)
Set the absolute position of the detector with given index. Not thread safe.
bool isMonitor(const size_t index) const
Returns true if the detector is a monitor.
Interface class for detector objects.
Definition: IDetector.h:43
Records the filename and the description of failure.
Definition: Exception.h:98
bool isAscii() const
Returns true if the descriptor is looking at an ascii file.
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 information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
Defines a wrapper around a file whose internal structure can be accessed using the NeXus API.
static bool isReadable(const std::string &filename, const Version version=AnyVersion)
Returns true if the file is considered to store data in a hierarchy.
Iterator begin()
Iterator referring to first element in the container.
@ TOK_TRIM
remove leading and trailing whitespace from tokens
Iterator end()
Iterator referring to the past-the-end element in the container.
std::size_t count() const
Get the total number of tokens.
Class for 3D vectors.
Definition: V3D.h:34
void spherical(const double R, const double theta, const double phi) noexcept
Sets the vector position based on spherical coordinates.
Definition: V3D.cpp:57
std::shared_ptr< Instrument > Instrument_sptr
Shared pointer to an instrument object.
std::shared_ptr< Mantid::Geometry::IDetector > IDetector_sptr
Shared pointer to IDetector.
Definition: IDetector.h:100
std::unordered_map< specnum_t, size_t > spec2index_map
Map with key = spectrum number, value = workspace index.
Generate a tableworkspace to store the calibration results.
Simple structure to store information about the ASCII file header.
@ InOut
Both an input & output workspace.
Definition: Property.h:55