Mantid
Loading...
Searching...
No Matches
EQSANSMonitorTOF.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 "MantidAPI/Run.h"
14
15using namespace Mantid::Kernel;
16using namespace Mantid::Geometry;
17
19
20// Register the algorithm into the AlgorithmFactory
21DECLARE_ALGORITHM(EQSANSMonitorTOF)
22
23using namespace Kernel;
24using namespace API;
25using namespace Geometry;
26
28 declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input,
29 std::make_shared<WorkspaceUnitValidator>("TOF")),
30 "Workspace to apply the TOF correction to");
31
32 // Output parameters
33 declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
34 "Workspace to store the corrected data in");
35 declareProperty("FrameSkipping", false, "True if the data was taken in frame-skipping mode",
37}
38
40 MatrixWorkspace_const_sptr inputWS = getProperty("InputWorkspace");
41
42 // Now create the output workspace
43 MatrixWorkspace_sptr outputWS = getProperty("OutputWorkspace");
44 if (outputWS != inputWS) {
45 outputWS = WorkspaceFactory::Instance().create(inputWS);
46 setProperty("OutputWorkspace", outputWS);
47 }
48
49 // Get the nominal sample-to-detector distance (in mm)
50 // const double MD = MONITORPOS/1000.0;
51
52 // Get the monitor
53 const std::vector<detid_t> monitor_list = inputWS->getInstrument()->getMonitors();
54 if (monitor_list.size() != 1)
55 g_log.error() << "EQSANS workspace does not have exactly ones monitor! "
56 "This should not happen\n";
57
58 const auto &detInfo = inputWS->detectorInfo();
59 const size_t monIndex0 = detInfo.indexOf(0);
60 if (!detInfo.isMonitor(monIndex0)) {
61 g_log.error() << "Spectrum number " << monIndex0 << " has no detector assigned to it - discarding\n";
62 return;
63 }
64
65 // Get the source to monitor distance in mm
66 double source_z = inputWS->getInstrument()->getSource()->getPos().Z();
67 double monitor_z = detInfo.position(monIndex0).Z();
68 double source_to_monitor = (monitor_z - source_z) * 1000.0;
69
70 // Calculate the frame width
71 auto log = inputWS->run().getTimeSeriesProperty<double>("frequency");
72 double frequency = log->getStatistics().mean;
73 double tof_frame_width = 1.0e6 / frequency;
74
75 // Determine whether we need frame skipping or not by checking the chopper
76 // speed
77 bool frame_skipping = false;
78 log = inputWS->run().getTimeSeriesProperty<double>("Speed1");
79 const double chopper_speed = log->getStatistics().mean;
80 if (std::fabs(chopper_speed - frequency / 2.0) < 1.0)
81 frame_skipping = true;
82
83 // Get TOF offset
84 // this is the call to the chopper code to say where
85 // the start of the data frame is relative to the native facility frame
86 double frame_tof0 = getTofOffset(inputWS, frame_skipping, source_to_monitor);
87
88 // Calculate the frame width
89 // none of this changes in response to just looking at the monitor
90 double tmp_frame_width = frame_skipping ? tof_frame_width * 2.0 : tof_frame_width;
91 double frame_offset = 0.0;
92 if (frame_tof0 >= tmp_frame_width)
93 frame_offset = tmp_frame_width * (static_cast<int>(frame_tof0 / tmp_frame_width));
94
95 // Find the new binning first
96 const MantidVec XIn = inputWS->readX(0); // Copy here to avoid holding on to
97 // reference for too long (problem
98 // with managed workspaces)
99
100 // Since we are swapping the low-TOF and high-TOF regions around the cutoff
101 // value,
102 // there is the potential for having an overlap between the two regions. We
103 // exclude
104 // the region beyond a single frame by considering only the first 1/60 sec of
105 // the
106 // TOF histogram. (Bins 1 to 1666, as opposed to 1 to 2000)
107 const auto nTOF = static_cast<int>(XIn.size());
108
109 // Loop through each bin to find the cutoff where the TOF distribution wraps
110 // around
111 int cutoff = 0;
112 double threshold = frame_tof0 - frame_offset;
113 int tof_bin_range = 0;
114 double frame = 1000000.0 / frequency;
115 for (int i = 0; i < nTOF; i++) {
116 if (XIn[i] < threshold)
117 cutoff = i;
118 if (XIn[i] < frame)
119 tof_bin_range = i;
120 }
121 g_log.information() << "Cutoff=" << cutoff << "; Threshold=" << threshold << '\n';
122 g_log.information() << "Low TOFs: old = [" << (cutoff + 1) << ", " << (tof_bin_range - 2) << "] -> new = [0, "
123 << (tof_bin_range - 3 - cutoff) << "]\n";
124 g_log.information() << "High bin boundary of the Low TOFs: old = " << tof_bin_range - 1
125 << "; new = " << (tof_bin_range - 2 - cutoff) << '\n';
126 g_log.information() << "High TOFs: old = [0, " << (cutoff - 1) << "] -> new = [" << (tof_bin_range - 1 - cutoff)
127 << ", " << (tof_bin_range - 2) << "]\n";
128 g_log.information() << "Overlap: new = [" << (tof_bin_range - 1) << ", " << (nTOF - 2) << "]\n";
129
130 // Keep a copy of the input data since we may end up overwriting it
131 // if the input workspace is equal to the output workspace.
132 // This is necessary since we are shuffling around the TOF bins.
133 const MantidVec YIn = MantidVec(inputWS->readY(0));
134 const MantidVec EIn = MantidVec(inputWS->readE(0));
135
136 MantidVec &XOut = outputWS->dataX(0);
137 MantidVec &YOut = outputWS->dataY(0);
138 MantidVec &EOut = outputWS->dataE(0);
139
140 // Here we modify the TOF according to the offset we calculated.
141 // Since this correction will change the order of the TOF bins,
142 // we do it in sequence so that we obtain a valid distribution
143 // as our result (with increasing TOF values).
144
145 // Move up the low TOFs
146 for (int i = 0; i < cutoff; i++) {
147 XOut[i + tof_bin_range - 1 - cutoff] = XIn[i] + frame_offset + tmp_frame_width;
148 YOut[i + tof_bin_range - 1 - cutoff] = YIn[i];
149 EOut[i + tof_bin_range - 1 - cutoff] = EIn[i];
150 }
151
152 // Get rid of extra bins
153 for (int i = tof_bin_range - 1; i < nTOF - 1; i++) {
154 XOut[i] = XOut[i - 1] + 10.0;
155 YOut[i] = 0.0;
156 EOut[i] = 0.0;
157 }
158 XOut[nTOF - 1] = XOut[nTOF - 2] + 10.0;
159
160 // Move down the high TOFs
161 for (int i = cutoff + 1; i < tof_bin_range - 1; i++) {
162 XOut[i - cutoff - 1] = XIn[i] + frame_offset;
163 YOut[i - cutoff - 1] = YIn[i];
164 EOut[i - cutoff - 1] = EIn[i];
165 }
166 // Don't forget the low boundary
167 XOut[tof_bin_range - 2 - cutoff] = XIn[tof_bin_range] + frame_offset;
168
169 // Zero out the cutoff bin, which no longer makes sense because
170 // len(x) = len(y)+1
171 YOut[tof_bin_range - 2 - cutoff] = 0.0;
172 EOut[tof_bin_range - 2 - cutoff] = 0.0;
173
174 setProperty("OutputWorkspace", outputWS);
175}
176
177double EQSANSMonitorTOF::getTofOffset(const MatrixWorkspace_const_sptr &inputWS, bool frame_skipping,
178 double source_to_monitor) {
179 //# Storage for chopper information read from the logs
180 double chopper_set_phase[4] = {0, 0, 0, 0};
181 double chopper_speed[4] = {0, 0, 0, 0};
182 double chopper_actual_phase[4] = {0, 0, 0, 0};
183 double chopper_wl_1[4] = {0, 0, 0, 0};
184 double chopper_wl_2[4] = {0, 0, 0, 0};
185 double frame_wl_1 = 0;
186 double frame_srcpulse_wl_1 = 0;
187 double frame_wl_2 = 0;
188 double chopper_srcpulse_wl_1[4] = {0, 0, 0, 0};
189 double chopper_frameskip_wl_1[4] = {0, 0, 0, 0};
190 double chopper_frameskip_wl_2[4] = {0, 0, 0, 0};
191 double chopper_frameskip_srcpulse_wl_1[4] = {0, 0, 0, 0};
192
193 // Calculate the frame width
194 auto log = inputWS->run().getTimeSeriesProperty<double>("frequency");
195 double frequency = log->getStatistics().mean;
196 double tof_frame_width = 1.0e6 / frequency;
197
198 double tmp_frame_width = tof_frame_width;
199 if (frame_skipping)
200 tmp_frame_width *= 2.0;
201
202 // Choice of parameter set
203 int m_set = 0;
204 if (frame_skipping)
205 m_set = 1;
206
207 bool first = true;
208 bool first_skip = true;
209 double frameskip_wl_1 = 0;
210 double frameskip_srcpulse_wl_1 = 0;
211 double frameskip_wl_2 = 0;
212
213 for (int i = 0; i < 4; i++) {
214 // Read chopper information
215 std::ostringstream phase_str;
216 phase_str << "Phase" << i + 1;
217 log = inputWS->run().getTimeSeriesProperty<double>(phase_str.str());
218 chopper_set_phase[i] = log->getStatistics().mean;
219 std::ostringstream speed_str;
220 speed_str << "Speed" << i + 1;
221 log = inputWS->run().getTimeSeriesProperty<double>(speed_str.str());
222 chopper_speed[i] = log->getStatistics().mean;
223
224 // Only process choppers with non-zero speed
225 if (chopper_speed[i] <= 0)
226 continue;
227
228 chopper_actual_phase[i] = chopper_set_phase[i] - CHOPPER_PHASE_OFFSET[m_set][i];
229
230 while (chopper_actual_phase[i] < 0)
231 chopper_actual_phase[i] += tmp_frame_width;
232
233 double x1 = (chopper_actual_phase[i] - (tmp_frame_width * 0.5 * CHOPPER_ANGLE[i] / 360.)); // opening edge
234 double x2 = (chopper_actual_phase[i] + (tmp_frame_width * 0.5 * CHOPPER_ANGLE[i] / 360.)); // closing edge
235 if (!frame_skipping) // not skipping
236 {
237 while (x1 < 0) {
238 x1 += tmp_frame_width;
239 x2 += tmp_frame_width;
240 }
241 }
242
243 if (x1 > 0) {
244 chopper_wl_1[i] = 3.9560346 * x1 / CHOPPER_LOCATION[i];
245 chopper_srcpulse_wl_1[i] = 3.9560346 * (x1 - chopper_wl_1[i] * PULSEWIDTH) / CHOPPER_LOCATION[i];
246 } else
247 chopper_wl_1[i] = chopper_srcpulse_wl_1[i] = 0.;
248
249 chopper_wl_2[i] = (x2 > 0) ? 3.9560346 * x2 / CHOPPER_LOCATION[i] : 0.;
250
251 if (first) {
252 frame_wl_1 = chopper_wl_1[i];
253 frame_srcpulse_wl_1 = chopper_srcpulse_wl_1[i];
254 frame_wl_2 = chopper_wl_2[i];
255 first = false;
256 } else {
257 if (frame_skipping && i == 2) // ignore chopper 1 and 2 forthe shortest wl.
258 {
259 frame_wl_1 = chopper_wl_1[i];
260 frame_srcpulse_wl_1 = chopper_srcpulse_wl_1[i];
261 }
262 if (frame_wl_1 < chopper_wl_1[i])
263 frame_wl_1 = chopper_wl_1[i];
264 if (frame_wl_2 > chopper_wl_2[i])
265 frame_wl_2 = chopper_wl_2[i];
266 if (frame_srcpulse_wl_1 < chopper_srcpulse_wl_1[i])
267 frame_srcpulse_wl_1 = chopper_srcpulse_wl_1[i];
268 }
269
270 if (frame_skipping) {
271 if (x1 > 0) {
272 x1 += tof_frame_width; // skipped pulse
273 chopper_frameskip_wl_1[i] = 3.9560346 * x1 / CHOPPER_LOCATION[i];
274 chopper_frameskip_srcpulse_wl_1[i] = 3.9560346 * (x1 - chopper_wl_1[i] * PULSEWIDTH) / CHOPPER_LOCATION[i];
275 } else
276 chopper_wl_1[i] = chopper_srcpulse_wl_1[i] = 0.;
277
278 if (x2 > 0) {
279 x2 += tof_frame_width;
280 chopper_frameskip_wl_2[i] = 3.9560346 * x2 / CHOPPER_LOCATION[i];
281 } else
282 chopper_wl_2[i] = 0.;
283
284 if (i < 2 && chopper_frameskip_wl_1[i] > chopper_frameskip_wl_2[i])
285 continue;
286
287 if (first_skip) {
288 frameskip_wl_1 = chopper_frameskip_wl_1[i];
289 frameskip_srcpulse_wl_1 = chopper_frameskip_srcpulse_wl_1[i];
290 frameskip_wl_2 = chopper_frameskip_wl_2[i];
291 first_skip = false;
292 } else {
293 if (i == 2) // ignore chopper 1 and 2 forthe longest wl.
294 frameskip_wl_2 = chopper_frameskip_wl_2[i];
295
296 if (chopper_frameskip_wl_1[i] < chopper_frameskip_wl_2[i] && frameskip_wl_1 < chopper_frameskip_wl_1[i])
297 frameskip_wl_1 = chopper_frameskip_wl_1[i];
298
299 if (chopper_frameskip_wl_1[i] < chopper_frameskip_wl_2[i] &&
300 frameskip_srcpulse_wl_1 < chopper_frameskip_srcpulse_wl_1[i])
301 frameskip_srcpulse_wl_1 = chopper_frameskip_srcpulse_wl_1[i];
302
303 if (frameskip_wl_2 > chopper_frameskip_wl_2[i])
304 frameskip_wl_2 = chopper_frameskip_wl_2[i];
305 }
306 }
307 }
308
309 if (frame_wl_1 >= frame_wl_2) // too many frames later. So figure it out
310 {
311 double n_frame[4] = {0, 0, 0, 0};
312 double c_wl_1[4] = {0, 0, 0, 0};
313 double c_wl_2[4] = {0, 0, 0, 0};
314 bool passed = false;
315
316 do {
317 frame_wl_1 = c_wl_1[0] = chopper_wl_1[0] + 3.9560346 * n_frame[0] * tof_frame_width / CHOPPER_LOCATION[0];
318 frame_wl_2 = c_wl_2[0] = chopper_wl_2[0] + 3.9560346 * n_frame[0] * tof_frame_width / CHOPPER_LOCATION[0];
319
320 for (int i = 1; i < 4; i++) {
321 n_frame[i] = n_frame[i - 1] - 1;
322 passed = false;
323
324 do {
325 n_frame[i] += 1;
326 c_wl_1[i] = chopper_wl_1[i] + 3.9560346 * n_frame[i] * tof_frame_width / CHOPPER_LOCATION[i];
327 c_wl_2[i] = chopper_wl_2[i] + 3.9560346 * n_frame[i] * tof_frame_width / CHOPPER_LOCATION[i];
328
329 if (frame_wl_1 < c_wl_2[i] && frame_wl_2 > c_wl_1[i]) {
330 passed = true;
331 break;
332 }
333 if (frame_wl_2 < c_wl_1[i])
334 break; // over shot
335 } while (n_frame[i] - n_frame[i - 1] < 10);
336
337 if (!passed) {
338 n_frame[0] += 1;
339 break;
340 } else {
341 if (frame_wl_1 < c_wl_1[i])
342 frame_wl_1 = c_wl_1[i];
343 if (frame_wl_2 > c_wl_2[i])
344 frame_wl_2 = c_wl_2[i];
345 }
346 }
347 } while (!passed && n_frame[0] < 99);
348
349 if (frame_wl_2 > frame_wl_1) {
350 int n = 3;
351 if (c_wl_1[2] > c_wl_1[3])
352 n = 2;
353
354 frame_srcpulse_wl_1 = c_wl_1[n] - 3.9560346 * c_wl_1[n] * PULSEWIDTH / CHOPPER_LOCATION[n];
355
356 for (int i = 0; i < 4; i++) {
357 chopper_wl_1[i] = c_wl_1[i];
358 chopper_wl_2[i] = c_wl_2[i];
359 if (frame_skipping) {
360 chopper_frameskip_wl_1[i] = c_wl_1[i] + 3.9560346 * 2. * tof_frame_width / CHOPPER_LOCATION[i];
361 chopper_frameskip_wl_2[i] = c_wl_2[i] + 3.9560346 * 2. * tof_frame_width / CHOPPER_LOCATION[i];
362 if (i == 0) {
363 frameskip_wl_1 = chopper_frameskip_wl_1[i];
364 frameskip_wl_2 = chopper_frameskip_wl_2[i];
365 } else {
366 if (frameskip_wl_1 < chopper_frameskip_wl_1[i])
367 frameskip_wl_1 = chopper_frameskip_wl_1[i];
368 if (frameskip_wl_2 > chopper_frameskip_wl_2[i])
369 frameskip_wl_2 = chopper_frameskip_wl_2[i];
370 }
371 }
372 }
373 } else
374 frame_srcpulse_wl_1 = 0.0;
375 }
376
377 double frame_tof0 = frame_srcpulse_wl_1 / 3.9560346 * source_to_monitor;
378
379 g_log.information() << "Frame width " << tmp_frame_width << '\n';
380 g_log.information() << "TOF offset = " << frame_tof0 << " microseconds\n";
381 g_log.information() << "Band defined by T1-T4 " << frame_wl_1 << " " << frame_wl_2;
382 if (frame_skipping)
383 g_log.information() << " + " << frameskip_wl_1 << " " << frameskip_wl_2 << '\n';
384 else
385 g_log.information() << '\n';
386 g_log.information() << "Chopper Actual Phase Lambda1 Lambda2\n";
387 for (int i = 0; i < 4; i++)
388 g_log.information() << i << " " << chopper_actual_phase[i] << " " << chopper_wl_1[i] << " " << chopper_wl_2[i]
389 << '\n';
390
391 setProperty("FrameSkipping", frame_skipping);
392
393 return frame_tof0;
394}
395
396} // namespace Mantid::WorkflowAlgorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
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
Kernel::Logger & g_log
Definition: Algorithm.h:451
A property class for workspaces.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
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
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
double getTofOffset(const API::MatrixWorkspace_const_sptr &inputWS, bool frame_skipping, double source_to_monitor)
Compute TOF offset.
void init() override
Initialisation code.
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
const double PULSEWIDTH
Determine the wavelength from the TOF in the beam monitor histogram.
const double CHOPPER_PHASE_OFFSET[2][4]
const double CHOPPER_LOCATION[4]
std::vector< double > MantidVec
typedef for the data storage used in Mantid matrix workspaces
Definition: cow_ptr.h:172
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54