Mantid
Loading...
Searching...
No Matches
SetScalingPSD.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 +
7// @author Ronald Fowler
9#include "LoadRaw/isisraw.h"
10#include "MantidAPI/Axis.h"
18#include "MantidKernel/V3D.h"
19
20#include <cmath>
21#include <fstream>
22
23namespace Mantid::DataHandling {
24
25// Register the algorithm into the algorithm factory
26DECLARE_ALGORITHM(SetScalingPSD)
27
28using namespace Kernel;
29using namespace API;
30using namespace Geometry;
31
33SetScalingPSD::SetScalingPSD() : Algorithm(), m_scalingOption(0) {}
34
39 // Declare required input parameters for algorithm
40 declareProperty(std::make_unique<FileProperty>("ScalingFilename", "", FileProperty::Load,
41 std::vector<std::string>{".sca", ".raw"}),
42 "The name of the scaling calibrations file to read, including its\n"
43 "full or relative path. The file extension must be either .sca or\n"
44 ".raw (filenames are case sensitive on linux)");
45 declareProperty(std::make_unique<WorkspaceProperty<>>("Workspace", "", Direction::InOut),
46 "The name of the workspace to apply the scaling to. This must be\n"
47 "associated with an instrument appropriate for the scaling file");
48
49 auto mustBePositive = std::make_shared<BoundedValidator<int>>();
50 mustBePositive->setLower(0);
51 declareProperty("ScalingOption", 0, mustBePositive,
52 "Control scaling calculation - 0 => use average of left and right\n"
53 "scaling (default). 1 => use maximum scaling. 2 => maximum + 5%");
54}
55
61 // Retrieve the filename and output workspace name from the properties
62 m_filename = getPropertyValue("ScalingFilename");
63 m_workspace = getProperty("Workspace");
64 m_scalingOption = getProperty("ScalingOption");
65 std::vector<Kernel::V3D> truepos;
67 // calculateDetectorShifts(truepos);
68}
69
76bool SetScalingPSD::processScalingFile(const std::string &scalingFile, std::vector<Kernel::V3D> &truePos) {
77 // Read the scaling information from a text file (.sca extension) or from a
78 // raw file (.raw)
79 // This is really corrected positions as (r,theta,phi) for each detector
80 // Compare these with the instrument values to determine the change in
81 // position and the scaling
82 // which may be necessary for each pixel if in a tube.
83 // movePos is used to updated positions
84 std::map<int, Kernel::V3D> posMap;
85 std::map<int, double> scaleMap;
86 std::map<int, double>::iterator its;
87
88 const auto &detectorInfo = m_workspace->detectorInfo();
89 if (scalingFile.find(".sca") != std::string::npos || scalingFile.find(".SCA") != std::string::npos) {
90 // read a .sca text format file
91 // format consists of a short header followed by one line per detector
92
93 std::ifstream sFile(scalingFile.c_str());
94 if (!sFile) {
95 g_log.error() << "Unable to open scaling file " << scalingFile << '\n';
96 return false;
97 }
98 std::string str;
99 // skip header line should be <filename> generated by <prog>
100 getline(sFile, str);
101 int detectorCount;
102 getline(sFile, str); // get detector count line
103 std::istringstream detCountStream(str);
104 detCountStream >> detectorCount;
105 if (detectorCount < 1) {
106 g_log.error("Bad detector count in scaling file");
107 throw std::runtime_error("Bad detector count in scaling file");
108 }
109 truePos.reserve(detectorCount);
110 getline(sFile, str); // skip title line
111 int detIdLast = -10;
112 Kernel::V3D truPosLast, detPosLast;
113
114 Progress prog(this, 0.0, 0.5, detectorCount);
115 // Now loop through lines, one for each detector/monitor. The latter are
116 // ignored.
117
118 while (getline(sFile, str)) {
119 if (str.empty() || str[0] == '#')
120 continue;
121 std::istringstream istr(str);
122
123 // read 6 values from the line to get the 3 (l2,theta,phi) of interest
124 int detIndex, code;
125 double l2, theta, phi, offset;
126 istr >> detIndex >> offset >> l2 >> code >> theta >> phi;
127
128 // sanity check on angles - l2 should be +ve but sample file has a few -ve
129 // values
130 // on monitors
131 if (theta > 181.0 || theta < -1 || phi < -181 || phi > 181) {
132 g_log.error("Position angle data out of range in .sca file");
133 throw std::runtime_error("Position angle data out of range in .sca file");
134 }
135 Kernel::V3D truPos;
136 // use abs as correction file has -ve l2 for first few detectors
137 truPos.spherical(fabs(l2), theta, phi);
138 truePos.emplace_back(truPos);
139 try {
140 // detIndex is what Mantid usually calls detectorID
141 size_t index = detectorInfo.indexOf(detIndex);
142 Kernel::V3D detPos = detectorInfo.position(index);
143 Kernel::V3D shift = truPos - detPos;
144
145 // scaling applied to dets that are not monitors and have sequential IDs
146 if (detIdLast == detIndex - 1 && !detectorInfo.isMonitor(index)) {
147 Kernel::V3D diffI = detPos - detPosLast;
148 Kernel::V3D diffT = truPos - truPosLast;
149 double scale = diffT.norm() / diffI.norm();
150 // Wish to store the scaling in a map, if we already have a scaling
151 // for this detector (i.e. from the other side) we average the two
152 // values. End of tube detectors only have one scaling estimate.
153 scaleMap[detIndex] = scale;
154 its = scaleMap.find(detIndex - 1);
155 if (its == scaleMap.end())
156 scaleMap[detIndex - 1] = scale;
157 else
158 its->second = 0.5 * (its->second + scale);
159 }
160 detIdLast = detIndex;
161 detPosLast = detPos;
162 truPosLast = truPos;
163 posMap[detIndex] = shift;
164 prog.report();
165 } catch (std::out_of_range &) {
166 continue;
167 }
168 }
169 } else if (scalingFile.find(".raw") != std::string::npos || scalingFile.find(".RAW") != std::string::npos) {
170 std::vector<int> detID;
171 std::vector<Kernel::V3D> truepos;
172 getDetPositionsFromRaw(scalingFile, detID, truepos);
173 //
174 auto detectorCount = static_cast<int>(detID.size());
175 if (detectorCount < 1) {
176 g_log.error("Failed to read any detectors from RAW file");
177 throw std::runtime_error("Failed to read any detectors from RAW file");
178 }
179 int detIdLast = -10;
180 Kernel::V3D truPosLast, detPosLast;
181 Progress prog(this, 0.0, 0.5, detectorCount);
182 for (int i = 0; i < detectorCount; i++) {
183 int detIndex = detID[i];
184 try {
185 // detIndex is what Mantid usually calls detectorID
186 size_t index = detectorInfo.indexOf(detIndex);
187 Kernel::V3D detPos = detectorInfo.position(index);
188 Kernel::V3D shift = truepos[i] - detPos;
189
190 if (detIdLast == detIndex - 1 && !detectorInfo.isMonitor(index)) {
191 Kernel::V3D diffI = detPos - detPosLast;
192 Kernel::V3D diffT = truepos[i] - truPosLast;
193 double scale = diffT.norm() / diffI.norm();
194 scaleMap[detIndex] = scale;
195 its = scaleMap.find(detIndex - 1);
196 if (its == scaleMap.end()) {
197 scaleMap[detIndex - 1] = scale;
198 } else {
199 if (m_scalingOption == 0)
200 its->second = 0.5 * (its->second + scale); // average of two
201 else if (m_scalingOption == 1) {
202 if (its->second < scale)
203 its->second = scale; // max
204 } else if (m_scalingOption == 2) {
205 if (its->second < scale)
206 its->second = scale;
207 its->second *= 1.05; // max+5%
208 } else
209 its->second = 3.0; // crazy test value
210 }
211 }
212 detIdLast = detID[i];
213 detPosLast = detPos;
214 truPosLast = truepos[i];
215 posMap[detIndex] = shift;
216 prog.report();
217 } catch (std::out_of_range &) {
218 continue;
219 }
220 }
221 }
222 movePos(m_workspace, posMap, scaleMap);
223 return true;
224}
225
233void SetScalingPSD::movePos(API::MatrixWorkspace_sptr &WS, std::map<int, Kernel::V3D> &posMap,
234 std::map<int, double> &scaleMap) {
235
236 auto &detectorInfo = WS->mutableDetectorInfo();
237 auto &componentInfo = WS->mutableComponentInfo();
238
239 double maxScale = -1e6, minScale = 1e6, aveScale = 0.0;
240 int scaleCount = 0;
241 Progress prog(this, 0.5, 1.0, static_cast<int>(detectorInfo.size()));
242
243 for (size_t i = 0; i < detectorInfo.size(); ++i) {
244 int idet = detectorInfo.detectorIDs()[i];
245
246 // Check if we have a shift, else do nothing.
247 auto itPos = posMap.find(idet);
248 if (itPos == posMap.end())
249 continue;
250
251 // Do a relative move
252 const auto newPosition = detectorInfo.position(i) + itPos->second;
253 detectorInfo.setPosition(i, newPosition);
254
255 // Set the "sca" instrument parameter
256 auto itScale = scaleMap.find(idet);
257 if (itScale != scaleMap.end()) {
258 const double scale = itScale->second;
259 if (minScale > scale)
260 minScale = scale;
261 if (maxScale < scale)
262 maxScale = scale;
263 aveScale += fabs(1.0 - scale);
264 componentInfo.setScaleFactor(i, V3D(1.0, scale, 1.0));
265 ++scaleCount;
266 }
267 prog.report();
268 }
269 g_log.debug() << "Range of scaling factors is " << minScale << " to " << maxScale << "\n";
270 if (0 != scaleCount) {
271 g_log.debug() << "Average abs scaling fraction is " << aveScale / scaleCount << "\n";
272 } else {
273 g_log.debug() << "Average abs scaling fraction cannot ba calculated "
274 "because the scale count is 0! Its value before dividing "
275 "by the count is "
276 << aveScale << "\n";
277 }
278}
279
285void SetScalingPSD::getDetPositionsFromRaw(const std::string &rawfile, std::vector<int> &detID,
286 std::vector<Kernel::V3D> &pos) {
287 (void)rawfile; // Avoid compiler warning
288
289 // open raw file
290 ISISRAW iraw(nullptr);
291 if (iraw.readFromFile(m_filename.c_str(), false) != 0) {
292 g_log.error("Unable to open file " + m_filename);
293 throw Exception::FileError("Unable to open File:", m_filename);
294 }
295 // get detector information
296 const int numDetector = iraw.i_det; // number of detectors
297 const int *const rawDetID = iraw.udet; // detector IDs
298 const float *const r = iraw.len2; // distance from sample
299 const float *const angle = iraw.tthe; // angle between indicent beam and
300 // direction from sample to detector
301 // (two-theta)
302 const float *const phi = iraw.ut;
303 // Is ut01 (=phi) present? Sometimes an array is present but has wrong values
304 // e.g.all 1.0 or all 2.0
305 bool phiPresent = iraw.i_use > 0 && phi[0] != 1.0 && phi[0] != 2.0;
306 if (!phiPresent) {
307 g_log.error("Unable to get Phi values from the raw file");
308 }
309 detID.reserve(numDetector);
310 pos.reserve(numDetector);
311 Kernel::V3D point;
312 for (int i = 0; i < numDetector; ++i) {
313 point.spherical(r[i], angle[i], phi[i]);
314 pos.emplace_back(point);
315 detID.emplace_back(rawDetID[i]);
316 }
317}
318
319} // namespace Mantid::DataHandling
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
std::map< DeltaEMode::Type, std::string > index
Definition: DeltaEMode.cpp:19
IntArray detectorCount
#define fabs(x)
Definition: Matrix.cpp:22
isis raw file.
Definition: isisraw.h:272
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
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
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
A property class for workspaces.
void exec() override
Overwrites Algorithm method.
void getDetPositionsFromRaw(const std::string &rawfile, std::vector< int > &detID, std::vector< Kernel::V3D > &pos)
read the positions of detectors defined in the raw file
bool processScalingFile(const std::string &scalingFile, std::vector< Kernel::V3D > &truePos)
Read the scaling information from a file (e.g.
std::string m_filename
The name and path of the input file.
Definition: SetScalingPSD.h:75
SetScalingPSD()
Default constructor.
void movePos(API::MatrixWorkspace_sptr &WS, std::map< int, Kernel::V3D > &posMap, std::map< int, double > &scaleMap)
apply the shifts in posMap to the detectors in WS
API::MatrixWorkspace_sptr m_workspace
Pointer to the workspace.
Definition: SetScalingPSD.h:79
int m_scalingOption
An integer option controlling the scaling method.
Definition: SetScalingPSD.h:77
void init() override
Overwrites Algorithm method.
Records the filename and the description of failure.
Definition: Exception.h:98
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 report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
Definition: ProgressBase.h:51
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
double norm() const noexcept
Definition: V3D.h:263
ISIS VMS raw file definitions.
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
@ InOut
Both an input & output workspace.
Definition: Property.h:55