Mantid
Loading...
Searching...
No Matches
LoadDetectorInfo.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"
13
14#include <Poco/Path.h>
15#include <boost/algorithm/string/compare.hpp>
16#include <boost/algorithm/string/predicate.hpp>
17#include <nexus/NeXusFile.hpp>
18
19#include <fstream>
20
21namespace Mantid::DataHandling {
22using namespace Kernel;
23using namespace API;
24using namespace Geometry;
25
26namespace {
27// Name of the offset parameter
28const char *DELAY_PARAM = "DelayTime";
29// Name of pressure parameter
30const char *PRESSURE_PARAM = "TubePressure";
31// Name of wall thickness parameter
32const char *THICKNESS_PARAM = "TubeThickness";
33} // namespace
34
35// Register the algorithm into the algorithm factory
36DECLARE_ALGORITHM(LoadDetectorInfo)
37
38
40 : Algorithm(), m_baseInstrument(), m_samplePos(), m_moveDets(false), m_workspace(), m_instDelta(-1.0),
41 m_instPressure(-1.0), m_instThickness(-1.0) {}
42
44
45 declareProperty(std::make_unique<WorkspaceProperty<>>("Workspace", "", Direction::InOut),
46 "The name of the workspace to that the detector information "
47 "will be loaded into.");
48
49 const std::vector<std::string> exts{".dat", ".raw", ".sca", ".nxs"};
50 declareProperty(std::make_unique<FileProperty>("DataFilename", "", FileProperty::Load, exts),
51 "A **raw, dat, nxs** or **sca** file that contains information about the "
52 "detectors in the "
53 "workspace. The description of **dat** and **nxs** file format is "
54 "provided below.");
55
56 declareProperty("RelocateDets", false,
57 "If true, the detectors are moved to "
58 "the positions specified in the file "
59 "defined by the field above.",
61}
62
65 std::string filename = getPropertyValue("DataFilename");
66 std::string ext = Poco::Path(filename).getExtension();
67 if (boost::iequals(ext, "dat") || boost::iequals(ext, "sca")) {
68 loadFromDAT(filename);
69 } else if (boost::iequals(ext, "raw")) {
70 loadFromRAW(filename);
71 } else if (boost::iequals(ext, "nxs")) {
72 loadFromIsisNXS(filename);
73 } else {
74 throw std::invalid_argument("Unknown file type with extension=." + ext);
75 }
76}
77
82 m_workspace = getProperty("Workspace");
83 m_moveDets = getProperty("RelocateDets");
84
85 // Cache base instrument
86 m_baseInstrument = m_workspace->getInstrument()->baseInstrument();
87 Geometry::IComponent_const_sptr sample = m_workspace->getInstrument()->getSample();
88 if (sample)
89 m_samplePos = sample->getPos();
90
91 // cache values of instrument level parameters so we only change then if they
92 // are different
93 const auto &pmap = m_workspace->constInstrumentParameters();
94 // delay
95 auto param = pmap.get(m_baseInstrument->getComponentID(), DELAY_PARAM);
96 if (param)
97 m_instDelta = param->value<double>();
98 // pressure
99 param = pmap.get(m_baseInstrument->getComponentID(), PRESSURE_PARAM);
100 if (param)
101 m_instPressure = param->value<double>();
102 // thickness
103 param = pmap.get(m_baseInstrument->getComponentID(), THICKNESS_PARAM);
104 if (param)
105 m_instThickness = param->value<double>();
106}
107
112void LoadDetectorInfo::loadFromDAT(const std::string &filename) {
113 std::ifstream datFile(filename.c_str());
114 if (!datFile) {
115 throw Exception::FileError("Unable to access dat file", filename);
116 }
117
118 std::string line;
119 // skip 3 lines of header info
120 for (int i = 0; i < 3; ++i)
121 getline(datFile, line);
122
123 // start loop over file
124 auto &pmap = m_workspace->instrumentParameters();
125 auto &wsDetInfo = m_workspace->mutableDetectorInfo();
126 while (getline(datFile, line)) {
127 if (line.empty() || line[0] == '#')
128 continue;
129
130 std::istringstream is(line);
131 detid_t detID(0);
132 int code(0);
133 float delta(0.0f), l2(0.0f), theta(0.0f), phi(0.0f);
134 is >> detID >> delta >> l2 >> code >> theta >> phi;
135 // offset value is be subtracted so store negative
136 delta *= -1.0f;
137
138 try {
139 size_t index = wsDetInfo.indexOf(detID);
140 if (wsDetInfo.isMonitor(index) || code == 1)
141 continue;
142
143 // drop 10 float columns
144 for (int i = 0; i < 10; ++i) {
145 float droppedFloat(0.0f);
146 is >> droppedFloat;
147 }
148
149 // pressure, wall thickness
150 float pressure(0.0), thickness(0.0);
151 is >> pressure >> thickness;
152
153 updateParameterMap(wsDetInfo, index, pmap, l2, theta, phi, delta, pressure, thickness);
154 } catch (std::out_of_range &) {
155 continue;
156 }
157 }
158}
159
163void LoadDetectorInfo::loadFromRAW(const std::string &filename) {
164 ISISRAW2 iraw;
165 if (iraw.readFromFile(filename.c_str(), false) != 0) {
166 throw Exception::FileError("Unable to access raw file:", filename);
167 }
168
169 const int numDets = iraw.i_det;
170 const int numUserTables = iraw.i_use;
171 int pressureTabNum(0), thicknessTabNum(0);
172 if (numUserTables == 10) {
173 pressureTabNum = 7;
174 thicknessTabNum = 8;
175 } else if (numUserTables == 14) {
176 pressureTabNum = 11;
177 thicknessTabNum = 12;
178 } else {
179 throw std::invalid_argument("RAW file contains unexpected number of user tables=" + std::to_string(numUserTables) +
180 ". Expected 10 or 14.");
181 }
182
183 // Is ut01 (=phi) present? Sometimes an array is present but has wrong values
184 // e.g.all 1.0 or all 2.0
185 bool phiPresent = (iraw.ut[0] != 1.0 && iraw.ut[0] != 2.0);
186
187 // Start loop over detectors
188 auto &pmap = m_workspace->instrumentParameters();
189 auto &wsDetInfo = m_workspace->mutableDetectorInfo();
190 for (int i = 0; i < numDets; ++i) {
191 auto detID = static_cast<detid_t>(iraw.udet[i]);
192 int code = iraw.code[i];
193 try {
194 size_t index = wsDetInfo.indexOf(detID);
195 if (wsDetInfo.isMonitor(index) || code == 1)
196 continue;
197
198 // Positions
199 float l2 = iraw.len2[i];
200 float theta = iraw.tthe[i];
201 float phi = (phiPresent ? iraw.ut[i] : 0.0f);
202
203 // Parameters
204 float delta = iraw.delt[i];
205 // offset value is be subtracted so store negative
206 delta *= -1.0f;
207 // pressure, wall thickness
208 float pressure = iraw.ut[i + pressureTabNum * numDets];
209 float thickness = iraw.ut[i + thicknessTabNum * numDets];
210
211 updateParameterMap(wsDetInfo, index, pmap, l2, theta, phi, delta, pressure, thickness);
212 } catch (std::out_of_range &) {
213 continue;
214 }
215 }
216}
217
222void LoadDetectorInfo::loadFromIsisNXS(const std::string &filename) {
223 ::NeXus::File nxsFile(filename,
224 NXACC_READ); // will throw if file can't be opened
225
226 // two types of file:
227 // - new type entry per detector
228 // - old libisis with single pressure, thickness entry for whole file
229
230 // hold data read from file
231 DetectorInfo detInfo;
232
233 std::map<std::string, std::string> entries;
234 nxsFile.getEntries(entries);
235 if (entries.find("full_reference_detector") != entries.end()) {
236 nxsFile.openGroup("full_reference_detector", "NXIXTdetector");
237 readLibisisNxs(nxsFile, detInfo);
238 nxsFile.closeGroup();
239 } else if (entries.find("detectors.dat") != entries.end()) {
240 nxsFile.openGroup("detectors.dat", "NXEntry");
241 readNXSDotDat(nxsFile, detInfo);
242 nxsFile.closeGroup();
243 } else {
244 throw std::invalid_argument("Unknown NeXus file type");
245 }
246 nxsFile.close();
247
248 // Start loop over detectors
249 auto &pmap = m_workspace->instrumentParameters();
250 auto &wsDetInfo = m_workspace->mutableDetectorInfo();
251 auto numDets = static_cast<int>(detInfo.ids.size());
252 for (int i = 0; i < numDets; ++i) {
253 detid_t detID = detInfo.ids[i];
254 int code = detInfo.codes[i];
255 try {
256 size_t index = wsDetInfo.indexOf(detID);
257 if (wsDetInfo.isMonitor(index) || code == 1)
258 continue;
259
260 // Positions
261 double l2 = detInfo.l2[i];
262 double theta = detInfo.theta[i];
263 double phi = detInfo.phi[i];
264
265 // Parameters
266 double delta = detInfo.delays[i];
267 // offset value is be subtracted so store negative
268 delta *= -1.0;
269 // pressure, wall thickness
270 double pressure = detInfo.pressures[i];
271 double thickness = detInfo.thicknesses[i];
272
273 updateParameterMap(wsDetInfo, index, pmap, l2, theta, phi, delta, pressure, thickness);
274 } catch (std::out_of_range &) {
275 continue;
276 }
277 }
278}
279
287void LoadDetectorInfo::readLibisisNxs(::NeXus::File &nxsFile, DetectorInfo &detInfo) const {
288 nxsFile.readData<int32_t>("det_no", detInfo.ids);
289 nxsFile.readData<int32_t>("det_type", detInfo.codes);
290 nxsFile.readData<double>("delay_time", detInfo.delays);
291 const size_t numDets = detInfo.ids.size();
292
293 if (m_moveDets) {
294 nxsFile.readData<double>("L2", detInfo.l2);
295 nxsFile.readData<double>("theta", detInfo.theta);
296 nxsFile.readData<double>("phi", detInfo.phi);
297 } else {
298 // these will get ignored
299 detInfo.l2.resize(numDets, -1.0);
300 detInfo.theta.resize(numDets, -1.0);
301 detInfo.phi.resize(numDets, -1.0);
302 }
303
304 // pressure & wall thickness are global here
305 double pressure = -1.0;
306 double thickness = -1.0;
307 nxsFile.openGroup("det_he3", "NXIXTdet_he3");
308 nxsFile.readData<double>("gas_pressure", pressure);
309 nxsFile.readData<double>("wall_thickness", thickness);
310 nxsFile.closeGroup();
311 if (pressure <= 0.0) {
312 g_log.warning("The data file does not contain correct He3 pressure, "
313 "default value of 10 bar is used instead");
314 pressure = 10.0;
315 }
316 if (thickness <= 0.0) {
317 g_log.warning("The data file does not contain correct detector's wall "
318 "thickness, default value of 0.8mm is used instead");
319 thickness = 0.0008;
320 }
321 detInfo.pressures.resize(numDets, pressure);
322 detInfo.thicknesses.resize(numDets, thickness);
323}
324
332void LoadDetectorInfo::readNXSDotDat(::NeXus::File &nxsFile, DetectorInfo &detInfo) const {
333 std::vector<int32_t> fileIDs;
334 nxsFile.readData<int32_t>("detID", fileIDs); // containts both ids & codes
335 std::vector<float> fileOffsets;
336 nxsFile.readData<float>("timeOffsets", fileOffsets);
337 const size_t numDets = fileOffsets.size() / 2;
338
339 std::vector<float> detCoords;
340 if (m_moveDets) {
341 nxsFile.readData<float>("detSphericalCoord", detCoords);
342 } else {
343 detCoords.resize(3 * numDets, -1.0f);
344 }
345
346 // pressure & wall thickness
347 std::vector<float> pressureAndWall;
348 nxsFile.readData<float>("detPressureAndWall", pressureAndWall);
349
350 if (numDets != fileIDs.size() / 2 || numDets != detCoords.size() / 3 || numDets != pressureAndWall.size() / 2) {
351 std::ostringstream os;
352 os << "The sizes of NeXus data columns are inconsistent in detectors.dat.\n"
353 << "detIDs=" << fileIDs.size() << ", offsets=" << fileOffsets.size()
354 << ", pressure & thickness=" << pressureAndWall.size() << "\n";
355 throw std::runtime_error(os.str());
356 }
357
358 detInfo.ids.resize(numDets);
359 detInfo.codes.resize(numDets);
360 detInfo.delays.resize(numDets);
361 detInfo.l2.resize(numDets);
362 detInfo.theta.resize(numDets);
363 detInfo.phi.resize(numDets);
364 detInfo.pressures.resize(numDets);
365 detInfo.thicknesses.resize(numDets);
366
368 for (int i = 0; i < static_cast<int>(numDets); i++) {
369 detInfo.ids[i] = fileIDs[2 * i];
370 detInfo.codes[i] = fileIDs[2 * i + 1];
371 detInfo.delays[i] = fileOffsets[2 * i];
372
373 detInfo.l2[i] = detCoords[3 * i + 0]; // L2,
374 detInfo.theta[i] = detCoords[3 * i + 1]; // Theta
375 detInfo.phi[i] = detCoords[3 * i + 2]; // Phi
376
377 detInfo.pressures[i] = pressureAndWall[2 * i]; // pressure;
378 detInfo.thicknesses[i] = pressureAndWall[2 * i + 1]; // wallThickness;
379 }
380}
381
394void LoadDetectorInfo::updateParameterMap(Geometry::DetectorInfo &detectorInfo, const size_t detIndex,
395 Geometry::ParameterMap &pmap, const double l2, const double theta,
396 const double phi, const double delay, const double pressure,
397 const double thickness) const {
398
399 const auto detCompID = detectorInfo.detector(detIndex).getComponentID();
400
401 // store detector params that are different to instrument level
402 if (fabs(delay - m_instDelta) > 1e-06)
403 pmap.addDouble(detCompID, DELAY_PARAM, delay);
404 if (fabs(pressure - m_instPressure) > 1e-06)
405 pmap.addDouble(detCompID, PRESSURE_PARAM, pressure);
406 if (fabs(thickness - m_instThickness) > 1e-06)
407 pmap.addDouble(detCompID, THICKNESS_PARAM, thickness);
408
409 // move
410 if (m_moveDets) {
411 V3D newPos;
412 newPos.spherical(l2, theta, phi);
413 // The sample position may not be at 0,0,0
414 newPos += m_samplePos;
415 detectorInfo.setPosition(detIndex, newPos);
416 }
417}
418} // namespace Mantid::DataHandling
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
std::map< DeltaEMode::Type, std::string > index
Definition: DeltaEMode.cpp:19
#define fabs(x)
Definition: Matrix.cpp:22
#define PARALLEL_FOR_NO_WSP_CHECK()
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
float * delt
hold off table (size NDET)
Definition: isisraw.h:303
int i_det
number of detectors NDET
Definition: isisraw.h:296
int * code
code for UTn tables (size NDET)
Definition: isisraw.h:305
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
Base class from which all concrete algorithm classes should be derived.
Definition: Algorithm.h:85
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
@ Load
allowed here which will be passed to the algorithm
Definition: FileProperty.h:52
A property class for workspaces.
double m_instThickness
Wall thickness value set on the instrument level.
void readNXSDotDat(::NeXus::File &nxsFile, DetectorInfo &detInfo) const
Read data from old-style libisis NeXus file.
double m_instDelta
Delta value that has been set on the instrument.
Geometry::Instrument_const_sptr m_baseInstrument
Cached instrument for this workspace.
void loadFromRAW(const std::string &filename)
Use a .raw file as input.
void cacheInputs()
Cache the user input that will be frequently accessed.
bool m_moveDets
If set to true then update the detector positions base on the information in the given file.
void updateParameterMap(Geometry::DetectorInfo &detectorInfo, const size_t detIndex, Geometry::ParameterMap &pmap, const double l2, const double theta, const double phi, const double delay, const double pressure, const double thickness) const
Update the parameter map with the new values for the given detector.
void readLibisisNxs(::NeXus::File &nxsFile, DetectorInfo &detInfo) const
Read data from old-style libisis NeXus file.
void exec() override
Virtual method - must be overridden by concrete algorithm.
Kernel::V3D m_samplePos
Cached sample position for this workspace.
API::MatrixWorkspace_sptr m_workspace
store a pointer to the user selected workspace
void loadFromIsisNXS(const std::string &filename)
Use a isis raw nexus or event file as input.
void init() override
Virtual method - must be overridden by concrete algorithm.
double m_instPressure
Pressure value set on the instrument level.
void loadFromDAT(const std::string &filename)
Use a .dat or .sca file as input.
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.
const Geometry::IDetector & detector(const size_t index) const
Return a const reference to the detector with given index.
virtual ComponentID getComponentID() const =0
Returns the ComponentID - a unique identifier of the component.
Records the filename and the description of failure.
Definition: Exception.h:98
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
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< const IComponent > IComponent_const_sptr
Typdef of a shared pointer to a const IComponent.
Definition: IComponent.h:161
int32_t detid_t
Typedef for a detector ID.
Definition: SpectrumInfo.h:21
Generate a tableworkspace to store the calibration results.
std::string to_string(const wide_integer< Bits, Signed > &n)
Simple data holder for passing the detector info around when dealing with the NeXus data.
@ InOut
Both an input & output workspace.
Definition: Property.h:55
@ Input
An input workspace.
Definition: Property.h:53