Mantid
Loading...
Searching...
No Matches
SaveIsawDetCal.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 +
11#include "MantidAPI/Run.h"
12#include "MantidAPI/Workspace.h"
17#include "MantidKernel/System.h"
18
19#include <Poco/File.h>
20#include <boost/algorithm/string/trim.hpp>
21#include <fstream>
22
23using namespace Mantid::Kernel;
24using namespace Mantid::API;
25using namespace Mantid::Geometry;
26using Mantid::Types::Core::DateAndTime;
27using std::string;
28
29namespace Mantid::DataHandling {
30
31// Register the algorithm into the AlgorithmFactory
33
34
36void SaveIsawDetCal::init() {
37 declareProperty(std::make_unique<WorkspaceProperty<Workspace>>("InputWorkspace", "", Direction::Input),
38 "An input workspace.");
39
40 declareProperty(std::make_unique<FileProperty>("Filename", "", FileProperty::Save, ".DetCal"),
41 "Path to an ISAW-style .detcal file to save.");
42
43 declareProperty("TimeOffset", 0.0, "Offsets to be applied to times");
44 declareProperty(std::make_unique<ArrayProperty<string>>("BankNames", Direction::Input),
45 "Optional: Only select the specified banks");
46 declareProperty("AppendFile", false,
47 "Append to file if true.\n"
48 "If false, new file (default).");
49}
50
54 std::string filename = getPropertyValue("Filename");
55 // MatrixWorkspace_sptr ws = getProperty("InputWorkspace");
56
57 Workspace_sptr ws1 = getProperty("InputWorkspace");
58 ExperimentInfo_sptr ws = std::dynamic_pointer_cast<ExperimentInfo>(ws1);
59
60 double T0 = getProperty("TimeOffset");
61 const API::Run &run = ws->run();
62 // Use T0 from workspace if T0 not specified in input
63 if (T0 == 0.0 && run.hasProperty("T0")) {
64 T0 = run.getPropertyValueAsType<double>("T0");
65 if (T0 != 0) {
66 g_log.notice() << "T0 = " << T0 << '\n';
67 }
68 }
69
70 std::vector<std::string> bankNames = getProperty("BankNames");
71
72 inst = ws->getInstrument();
73 if (!inst)
74 throw std::runtime_error("No instrument in the Workspace. Cannot save DetCal file.");
75 // We cannot assume the peaks have bank type detector modules, so we have a
76 // string to check this
77 std::string bankPart = "bank";
78 if (inst->getName() == "WISH")
79 bankPart = "WISHpanel";
80
81 std::set<int> uniqueBanks;
82 if (bankNames.empty()) {
83 // Get all children
84 std::vector<IComponent_const_sptr> comps;
85 inst->getChildren(comps, true);
86
87 for (auto &comp : comps) {
88 std::string bankName = comp->getName();
89 boost::trim(bankName);
90 boost::erase_all(bankName, bankPart);
91 int bank = 0;
92 Strings::convert(bankName, bank);
93 if (bank == 0)
94 continue;
95 // Track unique bank numbers
96 uniqueBanks.insert(bank);
97 }
98 } else {
99 for (auto bankName : bankNames) {
100 boost::trim(bankName);
101 boost::erase_all(bankName, bankPart);
102 int bank = 0;
103 Strings::convert(bankName, bank);
104 if (bank == 0)
105 continue;
106 // Track unique bank numbers
107 uniqueBanks.insert(bank);
108 }
109 }
110 if (!inst)
111 throw std::runtime_error("No instrument in PeaksWorkspace. Cannot save peaks file.");
112
113 double l1;
114 V3D beamline;
115 double beamline_norm;
116 V3D samplePos;
117 inst->getInstrumentParameters(l1, beamline, beamline_norm, samplePos);
118
119 std::ofstream out;
120 bool append = getProperty("AppendFile");
121
122 // do not append if file does not exist
123 if (!Poco::File(filename.c_str()).exists())
124 append = false;
125
126 if (append) {
127 out.open(filename.c_str(), std::ios::app);
128 } else {
129 out.open(filename.c_str());
130 out << "# NEW CALIBRATION FILE FORMAT (in NeXus/SNS coordinates):\n";
131 out << "# Lengths are in centimeters.\n";
132 out << "# Base and up give directions of unit vectors for a local \n";
133 out << "# x,y coordinate system on the face of the detector.\n";
134 out << "#\n";
135 out << "#\n";
136 out << "# " << DateAndTime::getCurrentTime().toISO8601String() << '\n';
137
138 out << "6 L1 T0_SHIFT\n";
139 out << "7 " << std::setw(10);
140 out << std::setprecision(4) << std::fixed << (l1 * 100);
141 out << std::setw(13) << std::setprecision(3) << T0 << '\n';
142
143 out << "4 DETNUM NROWS NCOLS WIDTH HEIGHT DEPTH DETD CenterX "
144 " CenterY CenterZ BaseX BaseY BaseZ UpX UpY "
145 " UpZ\n";
146 }
147 // Here would save each detector...
148 std::set<int>::iterator it;
149 for (it = uniqueBanks.begin(); it != uniqueBanks.end(); ++it) {
150 // Build up the bank name
151 int bank = *it;
152 std::ostringstream mess;
153 if (bankPart == "WISHpanel" && bank < 10)
154 mess << bankPart << "0" << bank;
155 else
156 mess << bankPart << bank;
157
158 std::string bankName = mess.str();
159 // Retrieve it
160 std::shared_ptr<const IComponent> det = inst->getComponentByName(bankName);
161 if (inst->getName() == "CORELLI") // for Corelli with sixteenpack under bank
162 {
163 std::vector<Geometry::IComponent_const_sptr> children;
164 std::shared_ptr<const Geometry::ICompAssembly> asmb =
165 std::dynamic_pointer_cast<const Geometry::ICompAssembly>(inst->getComponentByName(bankName));
166 asmb->getChildren(children, false);
167 det = children[0];
168 }
169 if (det) {
170 // Center of the detector
171 V3D center = det->getPos();
172
173 // Distance to center of detector
174 double detd = (center - inst->getSample()->getPos()).norm();
175 int NCOLS, NROWS;
176 double xsize, ysize;
177 sizeBanks(bankName, NCOLS, NROWS, xsize, ysize);
178 // Base unit vector (along the horizontal, X axis)
179 // Since WISH banks are curved use center and increment 2 for tubedown
180 int midX = NCOLS / 2;
181 int midY = NROWS / 2;
182 V3D base = findPixelPos(bankName, midX + 2, midY) - findPixelPos(bankName, midX, midY);
183 base.normalize();
184
185 // Up unit vector (along the vertical, Y axis)
186 V3D up = findPixelPos(bankName, midX, midY + 1) - findPixelPos(bankName, midX, midY);
187 up.normalize();
188
189 // Write the line
190 out << "5 " << std::setw(6) << std::right << bank << " " << std::setw(6) << std::right << NROWS << " "
191 << std::setw(6) << std::right << NCOLS << " " << std::setw(7) << std::right << std::fixed
192 << std::setprecision(4) << 100.0 * xsize << " " << std::setw(7) << std::right << std::fixed
193 << std::setprecision(4) << 100.0 * ysize << " "
194 << " 0.2000 " << std::setw(6) << std::right << std::fixed << std::setprecision(2) << 100.0 * detd << " "
195 << std::setw(9) << std::right << std::fixed << std::setprecision(4) << 100.0 * center.X() << " "
196 << std::setw(9) << std::right << std::fixed << std::setprecision(4) << 100.0 * center.Y() << " "
197 << std::setw(9) << std::right << std::fixed << std::setprecision(4) << 100.0 * center.Z() << " "
198 << std::setw(8) << std::right << std::fixed << std::setprecision(5) << base.X() << " " << std::setw(8)
199 << std::right << std::fixed << std::setprecision(5) << base.Y() << " " << std::setw(8) << std::right
200 << std::fixed << std::setprecision(5) << base.Z() << " " << std::setw(8) << std::right << std::fixed
201 << std::setprecision(5) << up.X() << " " << std::setw(8) << std::right << std::fixed << std::setprecision(5)
202 << up.Y() << " " << std::setw(8) << std::right << std::fixed << std::setprecision(5) << up.Z() << " \n";
203
204 } else
205 g_log.warning() << "Information about detector module " << bankName << " not found and recognised\n";
206 }
207
208 out.close();
209}
210
211V3D SaveIsawDetCal::findPixelPos(const std::string &bankName, int col, int row) {
212 std::shared_ptr<const IComponent> parent = inst->getComponentByName(bankName);
213 if (parent->type() == "RectangularDetector") {
214 std::shared_ptr<const RectangularDetector> RDet = std::dynamic_pointer_cast<const RectangularDetector>(parent);
215
216 std::shared_ptr<Detector> pixel = RDet->getAtXY(col, row);
217 return pixel->getPos();
218 } else {
219 std::vector<Geometry::IComponent_const_sptr> children;
220 std::shared_ptr<const Geometry::ICompAssembly> asmb =
221 std::dynamic_pointer_cast<const Geometry::ICompAssembly>(parent);
222 asmb->getChildren(children, false);
223 if (children[0]->getName() == "sixteenpack") {
224 asmb = std::dynamic_pointer_cast<const Geometry::ICompAssembly>(children[0]);
225 children.clear();
226 asmb->getChildren(children, false);
227 }
228 int col0 = col - 1;
229 // WISH detectors are in bank in this order in instrument
230 if (inst->getName() == "WISH")
231 col0 = (col % 2 == 0 ? col / 2 + 75 : (col - 1) / 2);
232 std::shared_ptr<const Geometry::ICompAssembly> asmb2 =
233 std::dynamic_pointer_cast<const Geometry::ICompAssembly>(children[col0]);
234 std::vector<Geometry::IComponent_const_sptr> grandchildren;
235 asmb2->getChildren(grandchildren, false);
236 Geometry::IComponent_const_sptr first = grandchildren[row - 1];
237 return first->getPos();
238 }
239}
240
241void SaveIsawDetCal::sizeBanks(const std::string &bankName, int &NCOLS, int &NROWS, double &xsize, double &ysize) {
242 if (bankName == "None")
243 return;
244 std::shared_ptr<const IComponent> parent = inst->getComponentByName(bankName);
245 if (parent->type() == "RectangularDetector") {
246 std::shared_ptr<const RectangularDetector> RDet = std::dynamic_pointer_cast<const RectangularDetector>(parent);
247
248 NCOLS = RDet->xpixels();
249 NROWS = RDet->ypixels();
250 xsize = RDet->xsize();
251 ysize = RDet->ysize();
252 } else {
253 std::vector<Geometry::IComponent_const_sptr> children;
254 std::shared_ptr<const Geometry::ICompAssembly> asmb =
255 std::dynamic_pointer_cast<const Geometry::ICompAssembly>(parent);
256 asmb->getChildren(children, false);
257 if (children[0]->getName() == "sixteenpack") {
258 asmb = std::dynamic_pointer_cast<const Geometry::ICompAssembly>(children[0]);
259 children.clear();
260 asmb->getChildren(children, false);
261 }
262 std::shared_ptr<const Geometry::ICompAssembly> asmb2 =
263 std::dynamic_pointer_cast<const Geometry::ICompAssembly>(children[0]);
264 std::vector<Geometry::IComponent_const_sptr> grandchildren;
265 asmb2->getChildren(grandchildren, false);
266 NROWS = static_cast<int>(grandchildren.size());
267 NCOLS = static_cast<int>(children.size());
268 Geometry::IComponent_const_sptr first = children[0];
269 Geometry::IComponent_const_sptr last = children[NCOLS - 1];
270 xsize = first->getDistance(*last);
271 first = grandchildren[0];
272 last = grandchildren[NROWS - 1];
273 ysize = first->getDistance(*last);
274 }
275}
276} // namespace Mantid::DataHandling
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
std::string getName(const IMDDimension &self)
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
@ Save
to specify a file to write to, the file may or may not exist
Definition: FileProperty.h:49
bool hasProperty(const std::string &name) const
Does the property exist on the object.
Definition: LogManager.cpp:265
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.
Saves an instrument with RectangularDetectors to an ISAW .DetCal file.
void sizeBanks(const std::string &bankName, int &NCOLS, int &NROWS, double &xsize, double &ysize)
Kernel::V3D findPixelPos(const std::string &bankName, int col, int row)
find position for rectangular and non-rectangular
void exec() override
Run the algorithm.
Geometry::Instrument_const_sptr inst
Support for a property that holds an array of values.
Definition: ArrayProperty.h:28
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
Class for 3D vectors.
Definition: V3D.h:34
constexpr double X() const noexcept
Get x.
Definition: V3D.h:232
double normalize()
Make a normalized vector (return norm value)
Definition: V3D.cpp:130
constexpr double Y() const noexcept
Get y.
Definition: V3D.h:233
constexpr double Z() const noexcept
Get z.
Definition: V3D.h:234
std::shared_ptr< Workspace > Workspace_sptr
shared pointer to Mantid::API::Workspace
Definition: Workspace_fwd.h:20
std::shared_ptr< ExperimentInfo > ExperimentInfo_sptr
Shared pointer to ExperimentInfo.
std::shared_ptr< const IComponent > IComponent_const_sptr
Typdef of a shared pointer to a const IComponent.
Definition: IComponent.h:161
int convert(const std::string &A, T &out)
Convert a string into a number.
Definition: Strings.cpp:665
@ Input
An input workspace.
Definition: Property.h:53