Mantid
Loading...
Searching...
No Matches
LoadIsawDetCal.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
12#include "MantidAPI/Run.h"
14
19
23
25#include "MantidKernel/V3D.h"
26
27#include <algorithm>
28#include <boost/algorithm/string/trim.hpp>
29#include <fstream>
30#include <numeric>
31#include <sstream>
32#include <utility>
33
34namespace Mantid::DataHandling {
35
36// Register the class into the algorithm factory
37DECLARE_ALGORITHM(LoadIsawDetCal)
38
39using namespace Kernel;
40using namespace API;
41using namespace Geometry;
42using namespace DataObjects;
43
47 declareProperty(std::make_unique<WorkspaceProperty<Workspace>>("InputWorkspace", "", Direction::InOut,
48 std::make_shared<InstrumentValidator>()),
49 "The workspace containing the geometry to be calibrated.");
50
51 const auto exts = std::vector<std::string>({".DetCal", ".detcal", ".peaks", ".integrate"});
52 declareProperty(std::make_unique<API::MultipleFileProperty>("Filename", exts),
53 "The input filename of the ISAW DetCal file (Two files "
54 "allowed for SNAP) ");
55
56 declareProperty(std::make_unique<API::FileProperty>("Filename2", "", API::FileProperty::OptionalLoad, exts),
57 "The input filename of the second ISAW DetCal file (West "
58 "banks for SNAP) ");
59
60 declareProperty("TimeOffset", 0.0, "Time Offset", Direction::Output);
61}
62
63namespace {
64const constexpr double DegreesPerRadian = 180.0 / M_PI;
65
66std::string getBankName(const std::string &bankPart, const int idnum) {
67 if (bankPart == "WISHpanel" && idnum < 10) {
68 return bankPart + "0" + std::to_string(idnum);
69 } else {
70 return bankPart + std::to_string(idnum);
71 }
72}
73
74std::string getInstName(const API::Workspace_const_sptr &wksp) {
75 MatrixWorkspace_const_sptr matrixWksp = std::dynamic_pointer_cast<const MatrixWorkspace>(wksp);
76 if (matrixWksp) {
77 return matrixWksp->getInstrument()->getName();
78 }
79
80 PeaksWorkspace_const_sptr peaksWksp = std::dynamic_pointer_cast<const PeaksWorkspace>(wksp);
81 if (peaksWksp) {
82 return peaksWksp->getInstrument()->getName();
83 }
84
85 throw std::runtime_error("Failed to determine instrument name");
86}
87} // namespace
88
89std::map<std::string, std::string> LoadIsawDetCal::validateInputs() {
90 std::map<std::string, std::string> result;
91
92 // two detcal files is only valid for snap
93 std::vector<std::string> filenames = getFilenames();
94 if (filenames.empty()) {
95 result["Filename"] = "Must supply .detcal file";
96 } else if (filenames.size() == 2) {
97 Workspace_const_sptr wksp = getProperty("InputWorkspace");
98 const auto instname = getInstName(wksp);
99
100 if (instname != "SNAP") {
101 result["Filename"] = "Two files is only valid for SNAP";
102 }
103 } else if (filenames.size() > 2) {
104 result["Filename"] = "Supply at most two .detcal files";
105 }
106
107 return result;
108}
109
115 // Get the input workspace
116 Workspace_sptr ws = getProperty("InputWorkspace");
117 MatrixWorkspace_sptr inputW = std::dynamic_pointer_cast<MatrixWorkspace>(ws);
118 PeaksWorkspace_sptr inputP = std::dynamic_pointer_cast<PeaksWorkspace>(ws);
119
120 Instrument_sptr inst = getCheckInst(ws);
121
122 std::string instname = inst->getName();
123
124 const auto filenames = getFilenames();
125
126 // Output summary to log file
127 int count, id, nrows, ncols;
128 double width, height, depth, detd, x, y, z, base_x, base_y, base_z, up_x, up_y, up_z;
129 std::ifstream input(filenames[0].c_str(), std::ios_base::in);
130 std::string line;
131 std::string detname;
132 // Build a list of Rectangular Detectors
133 std::vector<std::shared_ptr<RectangularDetector>> detList;
134 for (int i = 0; i < inst->nelements(); i++) {
135 std::shared_ptr<RectangularDetector> det;
136 std::shared_ptr<ICompAssembly> assem;
137 std::shared_ptr<ICompAssembly> assem2;
138
139 det = std::dynamic_pointer_cast<RectangularDetector>((*inst)[i]);
140 if (det) {
141 detList.emplace_back(det);
142 } else {
143 // Also, look in the first sub-level for RectangularDetectors (e.g. PG3).
144 // We are not doing a full recursive search since that will be very long
145 // for lots of pixels.
146 assem = std::dynamic_pointer_cast<ICompAssembly>((*inst)[i]);
147 if (assem) {
148 for (int j = 0; j < assem->nelements(); j++) {
149 det = std::dynamic_pointer_cast<RectangularDetector>((*assem)[j]);
150 if (det) {
151 detList.emplace_back(det);
152
153 } else {
154 // Also, look in the second sub-level for RectangularDetectors (e.g.
155 // PG3).
156 // We are not doing a full recursive search since that will be very
157 // long for lots of pixels.
158 assem2 = std::dynamic_pointer_cast<ICompAssembly>((*assem)[j]);
159 if (assem2) {
160 for (int k = 0; k < assem2->nelements(); k++) {
161 det = std::dynamic_pointer_cast<RectangularDetector>((*assem2)[k]);
162 if (det) {
163 detList.emplace_back(det);
164 }
165 }
166 }
167 }
168 }
169 }
170 }
171 }
172 std::unordered_set<int> uniqueBanks; // for CORELLI and WISH
173 std::string bankPart = "bank";
174 if (instname == "WISH")
175 bankPart = "WISHpanel";
176 if (detList.empty()) {
177 // Get all children
178 std::vector<IComponent_const_sptr> comps;
179 inst->getChildren(comps, true);
180
181 for (auto &comp : comps) {
182 std::string bankName = comp->getName();
183 boost::trim(bankName);
184 boost::erase_all(bankName, bankPart);
185 int bank = 0;
186 Strings::convert(bankName, bank);
187 if (bank == 0)
188 continue;
189 // Track unique bank numbers
190 uniqueBanks.insert(bank);
191 }
192 }
193
194 auto expInfoWS = std::dynamic_pointer_cast<ExperimentInfo>(ws);
195 auto &componentInfo = expInfoWS->mutableComponentInfo();
196 std::vector<ComponentScaling> rectangularDetectorScalings;
197
198 while (std::getline(input, line)) {
199 if (line[0] == '7') {
200 double mL1, mT0;
201 std::stringstream(line) >> count >> mL1 >> mT0;
202 setProperty("TimeOffset", mT0);
203 // Convert from cm to m
204 if (instname == "WISH")
205 center(0.0, 0.0, -mL1, "undulator", ws, componentInfo);
206 else
207 center(0.0, 0.0, -mL1, "moderator", ws, componentInfo);
208 // mT0 and time of flight are both in microsec
209 if (mT0 != 0.0) {
210 if (inputW) {
211 API::Run &run = inputW->mutableRun();
212 // Check to see if LoadEventNexus had T0 from TOPAZ Parameter file
213 auto alg1 = createChildAlgorithm("ChangeBinOffset");
214 alg1->setProperty<MatrixWorkspace_sptr>("InputWorkspace", inputW);
215 alg1->setProperty<MatrixWorkspace_sptr>("OutputWorkspace", inputW);
216 if (run.hasProperty("T0")) {
217 auto T0IDF = run.getPropertyValueAsType<double>("T0");
218 alg1->setProperty("Offset", mT0 - T0IDF);
219 } else {
220 alg1->setProperty("Offset", mT0);
221 }
222 alg1->executeAsChildAlg();
223 inputW = alg1->getProperty("OutputWorkspace");
224 // set T0 in the run parameters
225 run.addProperty<double>("T0", mT0, true);
226 } else if (inputP) {
227 // set T0 in the run parameters
228 API::Run &run = inputP->mutableRun();
229 run.addProperty<double>("T0", mT0, true);
230 }
231 }
232 }
233
234 if (line[0] != '5')
235 continue;
236
237 std::stringstream(line) >> count >> id >> nrows >> ncols >> width >> height >> depth >> detd >> x >> y >> z >>
238 base_x >> base_y >> base_z >> up_x >> up_y >> up_z;
239 if (id == 10 && filenames.size() == 2 && instname == "SNAP") {
240 input.close();
241 input.open(filenames[1].c_str());
242 while (std::getline(input, line)) {
243 if (line[0] != '5')
244 continue;
245
246 std::stringstream(line) >> count >> id >> nrows >> ncols >> width >> height >> depth >> detd >> x >> y >> z >>
247 base_x >> base_y >> base_z >> up_x >> up_y >> up_z;
248 if (id == 10)
249 break;
250 }
251 }
252 std::shared_ptr<RectangularDetector> det;
253 std::string bankName = getBankName(bankPart, id);
254 auto matchingDetector =
255 std::find_if(detList.begin(), detList.end(), [&bankName](const std::shared_ptr<RectangularDetector> &detector) {
256 return detector->getName() == bankName;
257 });
258 if (matchingDetector != detList.end()) {
259 det = *matchingDetector;
260 }
261
262 V3D rX(base_x, base_y, base_z);
263 V3D rY(up_x, up_y, up_z);
264
265 if (det) {
266 detname = det->getName();
267 center(x, y, z, detname, ws, componentInfo);
268
269 ComponentScaling detScaling;
270 detScaling.scaleX = CM_TO_M * width / det->xsize();
271 detScaling.scaleY = CM_TO_M * height / det->ysize();
272 detScaling.componentName = detname;
273 // Scaling will need both scale factors if LoadIsawPeaks or LoadIsawDetCal
274 // has already
275 // applied a calibration
276 if (inputW) {
277 Geometry::ParameterMap &pmap = inputW->instrumentParameters();
278 auto oldscalex = pmap.getDouble(detname, std::string("scalex"));
279 auto oldscaley = pmap.getDouble(detname, std::string("scaley"));
280 if (!oldscalex.empty())
281 detScaling.scaleX *= oldscalex[0];
282 if (!oldscaley.empty())
283 detScaling.scaleY *= oldscaley[0];
284 }
285 if (inputP) {
286 Geometry::ParameterMap &pmap = inputP->instrumentParameters();
287 auto oldscalex = pmap.getDouble(detname, std::string("scalex"));
288 auto oldscaley = pmap.getDouble(detname, std::string("scaley"));
289 if (!oldscalex.empty())
290 detScaling.scaleX *= oldscalex[0];
291 if (!oldscaley.empty())
292 detScaling.scaleY *= oldscaley[0];
293 }
294
295 rectangularDetectorScalings.emplace_back(detScaling);
296
297 doRotation(rX, rY, componentInfo, det);
298 }
299 auto bank = uniqueBanks.find(id);
300 if (bank == uniqueBanks.end())
301 continue;
302 int idnum = *bank;
303
304 bankName = getBankName(bankPart, idnum);
305 // Retrieve it
306 auto comp = inst->getComponentByName(bankName);
307 // for Corelli with sixteenpack under bank
308 if (instname == "CORELLI") {
309 std::vector<Geometry::IComponent_const_sptr> children;
310 std::shared_ptr<const Geometry::ICompAssembly> asmb =
311 std::dynamic_pointer_cast<const Geometry::ICompAssembly>(inst->getComponentByName(bankName));
312 asmb->getChildren(children, false);
313 comp = children[0];
314 }
315 if (comp) {
316 // Omitted scaling tubes
317 detname = comp->getFullName();
318 center(x, y, z, detname, ws, componentInfo);
319
320 bool doWishCorrection = (instname == "WISH"); // TODO: find out why this is needed for WISH
321 doRotation(rX, rY, componentInfo, comp, doWishCorrection);
322 }
323 }
324
325 // Do this last, to avoid the issue of invalidating DetectorInfo
326 applyScalings(ws, rectangularDetectorScalings);
327
328 setProperty("InputWorkspace", ws);
329}
330
341void LoadIsawDetCal::center(const double x, const double y, const double z, const std::string &detname,
342 const API::Workspace_sptr &ws, Geometry::ComponentInfo &componentInfo) {
343
344 Instrument_sptr inst = getCheckInst(ws);
345
346 IComponent_const_sptr comp = inst->getComponentByName(detname);
347 if (comp == nullptr) {
348 throw std::runtime_error("Component with name " + detname + " was not found.");
349 }
350
351 const V3D position(x * CM_TO_M, y * CM_TO_M, z * CM_TO_M);
352
353 const auto componentIndex = componentInfo.indexOf(comp->getComponentID());
354 componentInfo.setPosition(componentIndex, position);
355}
356
367 MatrixWorkspace_sptr inputW = std::dynamic_pointer_cast<MatrixWorkspace>(ws);
368 PeaksWorkspace_sptr inputP = std::dynamic_pointer_cast<PeaksWorkspace>(ws);
369
370 // Get some stuff from the input workspace
371 Instrument_sptr inst;
372 if (inputW) {
373 inst = std::const_pointer_cast<Instrument>(inputW->getInstrument());
374 if (!inst)
375 throw std::runtime_error("Could not get a valid instrument from the "
376 "MatrixWorkspace provided as input");
377 } else if (inputP) {
378 inst = std::const_pointer_cast<Instrument>(inputP->getInstrument());
379 if (!inst)
380 throw std::runtime_error("Could not get a valid instrument from the "
381 "PeaksWorkspace provided as input");
382 } else {
383 throw std::runtime_error("Could not get a valid instrument from the "
384 "workspace which does not seem to be valid as "
385 "input (must be either MatrixWorkspace or "
386 "PeaksWorkspace");
387 }
388
389 return inst;
390}
391
392std::vector<std::string> LoadIsawDetCal::getFilenames() {
393 std::vector<std::string> filenamesFromPropertyUnraveld;
394 std::vector<std::vector<std::string>> filenamesFromProperty = this->getProperty("Filename");
395 for (const auto &outer : filenamesFromProperty) {
396 std::copy(outer.begin(), outer.end(), std::back_inserter(filenamesFromPropertyUnraveld));
397 }
398
399 // shouldn't be used except for legacy cases
400 const std::string filename2 = this->getProperty("Filename2");
401 if (!filename2.empty())
402 filenamesFromPropertyUnraveld.emplace_back(filename2);
403
404 return filenamesFromPropertyUnraveld;
405}
406
417 const std::shared_ptr<const IComponent> &comp, bool doWishCorrection) {
418 // These are the ISAW axes
419 rX.normalize();
420 rY.normalize();
421
422 // These are the original axes
423 constexpr V3D oX(1., 0., 0.);
424 constexpr V3D oY(0., 1., 0.);
425
426 // Axis that rotates X
427 const V3D ax1 = oX.cross_prod(rX);
428 Quat Q1;
429 if (!ax1.nullVector(1e-12)) {
430 // Rotation angle from oX to rX
431 double angle1 = oX.angle(rX) * DegreesPerRadian;
432 if (doWishCorrection)
433 angle1 += 180.0;
434 // Create the first quaternion
435 Q1.setAngleAxis(angle1, ax1);
436 }
437
438 // Now we rotate the original Y using Q1
439 V3D roY = oY;
440 Q1.rotate(roY);
441 // Find the axis that rotates oYr onto rY
442 const V3D ax2 = roY.cross_prod(rY);
443 Quat Q2;
444 if (!ax2.nullVector(1e-12)) {
445 const double angle2 = roY.angle(rY) * DegreesPerRadian;
446 Q2.setAngleAxis(angle2, ax2);
447 }
448 // Final = those two rotations in succession; Q1 is done first.
449 const Quat Rot = Q2 * Q1;
450
451 // Then find the corresponding relative position
452 const auto componentIndex = componentInfo.indexOf(comp->getComponentID());
453
454 componentInfo.setRotation(componentIndex, Rot);
455}
456
467 const std::vector<ComponentScaling> &rectangularDetectorScalings) {
468
469 for (const auto &scaling : rectangularDetectorScalings) {
470 auto alg1 = createChildAlgorithm("ResizeRectangularDetector");
471 alg1->setProperty<Workspace_sptr>("Workspace", ws);
472 alg1->setProperty("ComponentName", scaling.componentName);
473 alg1->setProperty("ScaleX", scaling.scaleX);
474 alg1->setProperty("ScaleY", scaling.scaleY);
475 alg1->executeAsChildAlg();
476 }
477}
478
479} // namespace Mantid::DataHandling
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
double height
Definition: GetAllEi.cpp:155
double position
Definition: GetAllEi.cpp:154
int count
counter
Definition: Matrix.cpp:37
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
Definition: Algorithm.cpp:1913
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Definition: Algorithm.cpp:2076
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
@ OptionalLoad
to specify a file to read but the file doesn't have to exist
Definition: FileProperty.h:53
bool hasProperty(const std::string &name) const
Does the property exist on the object.
Definition: LogManager.cpp:265
void addProperty(Kernel::Property *prop, bool overwrite=false)
Add data to the object in the form of a property.
Definition: LogManager.h:79
HeldType getPropertyValueAsType(const std::string &name) const
Get the value of a property as the given TYPE.
Definition: LogManager.cpp:332
This class stores information regarding an experimental run as a series of log entries.
Definition: Run.h:38
A property class for workspaces.
std::vector< std::string > getFilenames()
Geometry::Instrument_sptr getCheckInst(const API::Workspace_sptr &ws)
Gets the instrument of the workspace, checking that the workspace and the instrument are as expected.
void init() override
Initialisation method.
void doRotation(Kernel::V3D rX, Kernel::V3D rY, Geometry::ComponentInfo &componentInfo, const std::shared_ptr< const Geometry::IComponent > &comp, bool doWishCorrection=false)
Perform the rotation for the calibration.
std::map< std::string, std::string > validateInputs() override
Perform validation of ALL the input properties of the algorithm.
void center(const double x, const double y, const double z, const std::string &detname, const API::Workspace_sptr &ws, Geometry::ComponentInfo &componentInfo)
Set the center of the supplied detector name.
void applyScalings(API::Workspace_sptr &ws, const std::vector< ComponentScaling > &rectangularDetectorScalings)
Apply the scalings from the calibration file.
void exec() override
Executes the algorithm.
ComponentInfo : Provides a component centric view on to the instrument.
Definition: ComponentInfo.h:40
void setRotation(size_t componentIndex, const Kernel::Quat &newRotation)
size_t indexOf(Geometry::IComponent *id) const
void setPosition(size_t componentIndex, const Kernel::V3D &newPosition)
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
Class for quaternions.
Definition: Quat.h:39
void rotate(V3D &) const
Rotate a vector.
Definition: Quat.cpp:397
void setAngleAxis(const double _deg, const V3D &_axis)
Constructor from an angle and axis.
Definition: Quat.cpp:114
Class for 3D vectors.
Definition: V3D.h:34
double normalize()
Make a normalized vector (return norm value)
Definition: V3D.cpp:130
constexpr V3D cross_prod(const V3D &v) const noexcept
Cross product (this * argument)
Definition: V3D.h:278
double angle(const V3D &) const
Angle between this and another vector.
Definition: V3D.cpp:165
bool nullVector(const double tolerance=1e-3) const noexcept
Determine if the point is null.
Definition: V3D.cpp:241
std::shared_ptr< Workspace > Workspace_sptr
shared pointer to Mantid::API::Workspace
Definition: Workspace_fwd.h:20
std::shared_ptr< const Workspace > Workspace_const_sptr
shared pointer to Mantid::API::Workspace (const version)
Definition: Workspace_fwd.h:22
std::shared_ptr< const MatrixWorkspace > MatrixWorkspace_const_sptr
shared pointer to the matrix workspace base class (const version)
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
std::shared_ptr< const PeaksWorkspace > PeaksWorkspace_const_sptr
Typedef for a shared pointer to a const peaks workspace.
std::shared_ptr< PeaksWorkspace > PeaksWorkspace_sptr
Typedef for a shared pointer to a peaks workspace.
std::shared_ptr< const IComponent > IComponent_const_sptr
Typdef of a shared pointer to a const IComponent.
Definition: IComponent.h:161
std::shared_ptr< Instrument > Instrument_sptr
Shared pointer to an instrument object.
int convert(const std::string &A, T &out)
Convert a string into a number.
Definition: Strings.cpp:665
Generate a tableworkspace to store the calibration results.
std::string to_string(const wide_integer< Bits, Signed > &n)
@ InOut
Both an input & output workspace.
Definition: Property.h:55
@ Output
An output workspace.
Definition: Property.h:54