Mantid
Loading...
Searching...
No Matches
EQSANSTofStructure.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"
16
17#include <vector>
18
19namespace Mantid::Algorithms {
20
21// Register the algorithm into the AlgorithmFactory
22DECLARE_ALGORITHM(EQSANSTofStructure)
23
24using namespace Kernel;
25using namespace API;
26using namespace DataObjects;
27using namespace Geometry;
28using Types::Event::TofEvent;
29
31 : API::Algorithm(), frame_tof0(0.), flight_path_correction(false), low_tof_cut(0.), high_tof_cut(0.) {}
32
34 declareProperty(std::make_unique<WorkspaceProperty<EventWorkspace>>("InputWorkspace", "", Direction::Input,
35 std::make_shared<WorkspaceUnitValidator>("TOF")),
36 "Workspace to apply the TOF correction to");
37 declareProperty("FlightPathCorrection", false, "If True, the neutron flight path correction will be applied",
39 declareProperty("LowTOFCut", 0.0,
40 "Width of the TOF margin to cut on the "
41 "lower end of the TOF distribution of each "
42 "frame",
44 declareProperty("HighTOFCut", 0.0,
45 "Width of the TOF margin to cut on the "
46 "upper end of the TOF distribution of "
47 "each frame",
49
50 // Output parameters
51 declareProperty("FrameSkipping", false, "If True, the data was taken in frame skipping mode",
53 declareProperty("TofOffset", 0.0, "TOF offset that was applied to the data", Kernel::Direction::Output);
54 declareProperty("WavelengthMin", 0.0, "Lower bound of the wavelength distribution of the first frame",
56 declareProperty("WavelengthMax", 0.0, "Upper bound of the wavelength distribution of the first frame",
58 declareProperty("WavelengthMinFrame2", 0.0, "Lower bound of the wavelength distribution of the second frame",
60 declareProperty("WavelengthMaxFrame2", 0.0, "Upper bound of the wavelength distribution of the second frame",
62}
63
65 EventWorkspace_sptr inputWS = getProperty("InputWorkspace");
66 flight_path_correction = getProperty("FlightPathCorrection");
67 low_tof_cut = getProperty("LowTOFCut");
68 high_tof_cut = getProperty("HighTOFCut");
69
70 // Calculate the frame width
71 auto frequencyLog = dynamic_cast<TimeSeriesProperty<double> *>(inputWS->run().getLogData("frequency"));
72 if (!frequencyLog) {
73 throw std::runtime_error("Frequency log not found.");
74 }
75 double frequency = frequencyLog->getStatistics().mean;
76 double tof_frame_width = 1.0e6 / frequency;
77
78 // Determine whether we need frame skipping or not by checking the chopper
79 // speed
80 bool frame_skipping = false;
81 auto chopper_speedLog = dynamic_cast<TimeSeriesProperty<double> *>(inputWS->run().getLogData("Speed1"));
82 if (!chopper_speedLog) {
83 throw std::runtime_error("Chopper speed log not found.");
84 }
85 const double chopper_speed = chopper_speedLog->getStatistics().mean;
86 if (std::fabs(chopper_speed - frequency / 2.0) < 1.0)
87 frame_skipping = true;
88
89 // Get TOF offset
90 frame_tof0 = getTofOffset(inputWS, frame_skipping);
91
92 // Calculate the frame width
93 double tmp_frame_width = frame_skipping ? tof_frame_width * 2.0 : tof_frame_width;
94 double frame_offset = 0.0;
95 if (frame_tof0 >= tmp_frame_width)
96 frame_offset = tmp_frame_width * (static_cast<int>(frame_tof0 / tmp_frame_width));
97
98 this->execEvent(inputWS, frame_tof0, frame_offset, tof_frame_width, tmp_frame_width, frame_skipping);
99}
100
102 double frame_offset, double tof_frame_width, double tmp_frame_width,
103 bool frame_skipping) {
104 const size_t numHists = inputWS->getNumberHistograms();
105 Progress progress(this, 0.0, 1.0, numHists);
106
107 // This now points to the correct distance and makes the naming clearer
108 // Get the nominal sample flange-to-detector distance (in mm)
109 Mantid::Kernel::Property *prop = inputWS->run().getProperty("sampleflange_detector_distance");
110 auto dp = dynamic_cast<Mantid::Kernel::PropertyWithValue<double> *>(prop);
111 if (!dp) {
112 throw std::runtime_error("sampleflange_detector_distance log not found.");
113 }
114 const double SFDD = *dp / 1000.0;
115
116 const auto &spectrumInfo = inputWS->spectrumInfo();
117 const auto l1 = spectrumInfo.l1();
118
119 // Loop through the spectra and apply correction
121 for (int64_t ispec = 0; ispec < int64_t(numHists); ++ispec) {
123
124 if (!spectrumInfo.hasDetectors(ispec)) {
125 g_log.warning() << "Workspace index " << ispec << " has no detector assigned to it - discarding\n";
126 continue;
127 }
128 const auto l2 = spectrumInfo.l2(ispec);
129 double tof_factor = (l1 + l2) / (l1 + SFDD);
130
131 // Get the pointer to the output event list
132 std::vector<TofEvent> &events = inputWS->getSpectrum(ispec).getEvents();
133 std::vector<TofEvent>::iterator it;
134 std::vector<TofEvent> clean_events;
135
136 for (it = events.begin(); it < events.end(); ++it) {
137 double newtof = it->tof();
138 newtof += frame_offset;
139 // Correct for the scattered neutron flight path
141 newtof /= tof_factor;
142
143 while (newtof < threshold)
144 newtof += tmp_frame_width;
145
146 // Remove events that don't fall within the accepted time window
147 double rel_tof = newtof - frame_tof0;
148 double x = (static_cast<int>(floor(rel_tof * 10)) % static_cast<int>(floor(tof_frame_width * 10))) * 0.1;
149 if (x < low_tof_cut || x > tof_frame_width - high_tof_cut) {
150 continue;
151 }
152 // At this point the events in the second frame are still off by a frame
153 if (frame_skipping && rel_tof > tof_frame_width)
154 newtof += tof_frame_width;
155 clean_events.emplace_back(newtof, it->pulseTime());
156 }
157 events.clear();
158 events.reserve(clean_events.size());
159 for (it = clean_events.begin(); it < clean_events.end(); ++it) {
160 events.emplace_back(*it);
161 }
162
163 progress.report("TOF structure");
165 }
167}
168
169double EQSANSTofStructure::getTofOffset(const EventWorkspace_const_sptr &inputWS, bool frame_skipping) {
170 //# Storage for chopper information read from the logs
171 double chopper_set_phase[4] = {0, 0, 0, 0};
172 double chopper_speed[4] = {0, 0, 0, 0};
173 double chopper_actual_phase[4] = {0, 0, 0, 0};
174 double chopper_wl_1[4] = {0, 0, 0, 0};
175 double chopper_wl_2[4] = {0, 0, 0, 0};
176 double frame_wl_1 = 0;
177 double frame_srcpulse_wl_1 = 0;
178 double frame_wl_2 = 0;
179 double chopper_srcpulse_wl_1[4] = {0, 0, 0, 0};
180 double chopper_frameskip_wl_1[4] = {0, 0, 0, 0};
181 double chopper_frameskip_wl_2[4] = {0, 0, 0, 0};
182 double chopper_frameskip_srcpulse_wl_1[4] = {0, 0, 0, 0};
183
184 // Calculate the frame width
185 auto frequencyLog = dynamic_cast<TimeSeriesProperty<double> *>(inputWS->run().getLogData("frequency"));
186 if (!frequencyLog) {
187 throw std::runtime_error("Frequency log not found.");
188 }
189 double frequency = frequencyLog->getStatistics().mean;
190 double tof_frame_width = 1.0e6 / frequency;
191
192 double tmp_frame_width = tof_frame_width;
193 if (frame_skipping)
194 tmp_frame_width *= 2.0;
195
196 // Choice of parameter set
197 int m_set = 0;
198 if (frame_skipping)
199 m_set = 1;
200
201 bool first = true;
202 bool first_skip = true;
203 double frameskip_wl_1 = 0;
204 double frameskip_srcpulse_wl_1 = 0;
205 double frameskip_wl_2 = 0;
206
207 for (int i = 0; i < 4; i++) {
208 // Read chopper information
209 std::ostringstream phase_str;
210 phase_str << "Phase" << i + 1;
211 auto log = dynamic_cast<TimeSeriesProperty<double> *>(inputWS->run().getLogData(phase_str.str()));
212 if (!log) {
213 throw std::runtime_error("Phase log not found.");
214 }
215 chopper_set_phase[i] = log->getStatistics().mean;
216 std::ostringstream speed_str;
217 speed_str << "Speed" << i + 1;
218 log = dynamic_cast<TimeSeriesProperty<double> *>(inputWS->run().getLogData(speed_str.str()));
219 if (!log) {
220 throw std::runtime_error("Speed log not found.");
221 }
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 if (x2 > 0)
250 chopper_wl_2[i] = 3.9560346 * x2 / CHOPPER_LOCATION[i];
251 else
252 chopper_wl_2[i] = 0.;
253
254 if (first) {
255 frame_wl_1 = chopper_wl_1[i];
256 frame_srcpulse_wl_1 = chopper_srcpulse_wl_1[i];
257 frame_wl_2 = chopper_wl_2[i];
258 first = false;
259 } else {
260 if (frame_skipping && i == 2) // ignore chopper 1 and 2 forthe shortest wl.
261 {
262 frame_wl_1 = chopper_wl_1[i];
263 frame_srcpulse_wl_1 = chopper_srcpulse_wl_1[i];
264 }
265 if (frame_wl_1 < chopper_wl_1[i])
266 frame_wl_1 = chopper_wl_1[i];
267 if (frame_wl_2 > chopper_wl_2[i])
268 frame_wl_2 = chopper_wl_2[i];
269 if (frame_srcpulse_wl_1 < chopper_srcpulse_wl_1[i])
270 frame_srcpulse_wl_1 = chopper_srcpulse_wl_1[i];
271 }
272
273 if (frame_skipping) {
274 if (x1 > 0) {
275 x1 += tof_frame_width; // skipped pulse
276 chopper_frameskip_wl_1[i] = 3.9560346 * x1 / CHOPPER_LOCATION[i];
277 chopper_frameskip_srcpulse_wl_1[i] = 3.9560346 * (x1 - chopper_wl_1[i] * PULSEWIDTH) / CHOPPER_LOCATION[i];
278 } else
279 chopper_wl_1[i] = chopper_srcpulse_wl_1[i] = 0.;
280
281 if (x2 > 0) {
282 x2 += tof_frame_width;
283 chopper_frameskip_wl_2[i] = 3.9560346 * x2 / CHOPPER_LOCATION[i];
284 } else
285 chopper_wl_2[i] = 0.;
286
287 if (i < 2 && chopper_frameskip_wl_1[i] > chopper_frameskip_wl_2[i])
288 continue;
289
290 if (first_skip) {
291 frameskip_wl_1 = chopper_frameskip_wl_1[i];
292 frameskip_srcpulse_wl_1 = chopper_frameskip_srcpulse_wl_1[i];
293 frameskip_wl_2 = chopper_frameskip_wl_2[i];
294 first_skip = false;
295 } else {
296 if (i == 2) // ignore chopper 1 and 2 forthe longest wl.
297 frameskip_wl_2 = chopper_frameskip_wl_2[i];
298
299 if (chopper_frameskip_wl_1[i] < chopper_frameskip_wl_2[i] && frameskip_wl_1 < chopper_frameskip_wl_1[i])
300 frameskip_wl_1 = chopper_frameskip_wl_1[i];
301
302 if (chopper_frameskip_wl_1[i] < chopper_frameskip_wl_2[i] &&
303 frameskip_srcpulse_wl_1 < chopper_frameskip_srcpulse_wl_1[i])
304 frameskip_srcpulse_wl_1 = chopper_frameskip_srcpulse_wl_1[i];
305
306 if (frameskip_wl_2 > chopper_frameskip_wl_2[i])
307 frameskip_wl_2 = chopper_frameskip_wl_2[i];
308 }
309 }
310 }
311
312 if (frame_wl_1 >= frame_wl_2) // too many frames later. So figure it out
313 {
314 double n_frame[4] = {0, 0, 0, 0};
315 double c_wl_1[4] = {0, 0, 0, 0};
316 double c_wl_2[4] = {0, 0, 0, 0};
317 bool passed = false;
318
319 do {
320 frame_wl_1 = c_wl_1[0] = chopper_wl_1[0] + 3.9560346 * n_frame[0] * tof_frame_width / CHOPPER_LOCATION[0];
321 frame_wl_2 = c_wl_2[0] = chopper_wl_2[0] + 3.9560346 * n_frame[0] * tof_frame_width / CHOPPER_LOCATION[0];
322
323 for (int i = 1; i < 4; i++) {
324 n_frame[i] = n_frame[i - 1] - 1;
325 passed = false;
326
327 do {
328 n_frame[i] += 1;
329 c_wl_1[i] = chopper_wl_1[i] + 3.9560346 * n_frame[i] * tof_frame_width / CHOPPER_LOCATION[i];
330 c_wl_2[i] = chopper_wl_2[i] + 3.9560346 * n_frame[i] * tof_frame_width / CHOPPER_LOCATION[i];
331
332 if (frame_wl_1 < c_wl_2[i] && frame_wl_2 > c_wl_1[i]) {
333 passed = true;
334 break;
335 }
336 if (frame_wl_2 < c_wl_1[i])
337 break; // over shot
338 } while (n_frame[i] - n_frame[i - 1] < 10);
339
340 if (!passed) {
341 n_frame[0] += 1;
342 break;
343 } else {
344 if (frame_wl_1 < c_wl_1[i])
345 frame_wl_1 = c_wl_1[i];
346 if (frame_wl_2 > c_wl_2[i])
347 frame_wl_2 = c_wl_2[i];
348 }
349 }
350 } while (!passed && n_frame[0] < 99);
351
352 if (frame_wl_2 > frame_wl_1) {
353 int n = 3;
354 if (c_wl_1[2] > c_wl_1[3])
355 n = 2;
356
357 frame_srcpulse_wl_1 = c_wl_1[n] - 3.9560346 * c_wl_1[n] * PULSEWIDTH / CHOPPER_LOCATION[n];
358
359 for (int i = 0; i < 4; i++) {
360 chopper_wl_1[i] = c_wl_1[i];
361 chopper_wl_2[i] = c_wl_2[i];
362 if (frame_skipping) {
363 chopper_frameskip_wl_1[i] = c_wl_1[i] + 3.9560346 * 2. * tof_frame_width / CHOPPER_LOCATION[i];
364 chopper_frameskip_wl_2[i] = c_wl_2[i] + 3.9560346 * 2. * tof_frame_width / CHOPPER_LOCATION[i];
365 if (i == 0) {
366 frameskip_wl_1 = chopper_frameskip_wl_1[i];
367 frameskip_wl_2 = chopper_frameskip_wl_2[i];
368 } else {
369 if (frameskip_wl_1 < chopper_frameskip_wl_1[i])
370 frameskip_wl_1 = chopper_frameskip_wl_1[i];
371 if (frameskip_wl_2 > chopper_frameskip_wl_2[i])
372 frameskip_wl_2 = chopper_frameskip_wl_2[i];
373 }
374 }
375 }
376 } else
377 frame_srcpulse_wl_1 = 0.0;
378 }
379 // Get source and detector locations
380 // get the name of the mapping file as set in the parameter files
381 std::vector<std::string> temp = inputWS->getInstrument()->getStringParameter("detector-name");
382 std::string det_name = "detector1";
383 if (temp.empty())
384 g_log.information() << "The instrument parameter file does not contain the "
385 "'detector-name' parameter: trying 'detector1'";
386 else
387 det_name = temp[0];
388
389 // Checked 8/11/2017 here detector_z is sfdd which has been updated
390 // in eqsansload.cpp
391 double source_z = inputWS->getInstrument()->getSource()->getPos().Z();
392 double detector_z = inputWS->getInstrument()->getComponentByName(det_name)->getPos().Z();
393
394 double source_to_detector = (detector_z - source_z) * 1000.0;
395 frame_tof0 = frame_srcpulse_wl_1 / 3.9560346 * source_to_detector;
396
397 g_log.information() << "Frame width " << tmp_frame_width << '\n';
398 g_log.information() << "TOF offset = " << frame_tof0 << " microseconds\n";
399 g_log.information() << "Band defined by T1-T4 " << frame_wl_1 << " " << frame_wl_2;
400 if (frame_skipping)
401 g_log.information() << " + " << frameskip_wl_1 << " " << frameskip_wl_2 << '\n';
402 else
403 g_log.information() << '\n';
404 g_log.information() << "Chopper Actual Phase Lambda1 Lambda2\n";
405 for (int i = 0; i < 4; i++)
406 g_log.information() << i << " " << chopper_actual_phase[i] << " " << chopper_wl_1[i] << " " << chopper_wl_2[i]
407 << '\n';
408
409 // Checked 8/10/2017
410 double low_wl_discard = 3.9560346 * low_tof_cut / source_to_detector;
411 double high_wl_discard = 3.9560346 * high_tof_cut / source_to_detector;
412
413 setProperty("FrameSkipping", frame_skipping);
414 setProperty("TofOffset", frame_tof0);
415 setProperty("WavelengthMin", frame_wl_1 + low_wl_discard);
416 setProperty("WavelengthMax", frame_wl_2 - high_wl_discard);
417 if (frame_skipping) {
418 setProperty("WavelengthMinFrame2", frameskip_wl_1 + low_wl_discard);
419 setProperty("WavelengthMaxFrame2", frameskip_wl_2 - high_wl_discard);
420 }
421
422 return frame_tof0;
423}
424
425} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
#define PARALLEL_START_INTERRUPT_REGION
Begins a block to skip processing is the algorithm has been interupted Note the end of the block if n...
Definition: MultiThreaded.h:94
#define PARALLEL_END_INTERRUPT_REGION
Ends a block to skip processing is the algorithm has been interupted Note the start of the block if n...
#define PARALLEL_FOR_IF(condition)
Empty definitions - to enable set your complier to enable openMP.
#define PARALLEL_CHECK_INTERRUPT_REGION
Adds a check after a Parallel region to see if it was interupted.
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
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
void progress(double p, const std::string &msg="", double estimatedTime=0.0, int progressPrecision=0)
Sends ProgressNotification.
Definition: Algorithm.cpp:231
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
A property class for workspaces.
void exec() override
Execution code.
void execEvent(const Mantid::DataObjects::EventWorkspace_sptr &inputWS, double threshold, double frame_offset, double tof_frame_width, double tmp_frame_width, bool frame_skipping)
void init() override
Initialisation code.
double getTofOffset(const DataObjects::EventWorkspace_const_sptr &inputWS, bool frame_skipping)
Compute TOF offset.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
The concrete, templated class for properties.
Base class for properties.
Definition: Property.h:94
A specialised Property class for holding a series of time-value pairs.
const double PULSEWIDTH
Apply correction to EQSANS data to account for its TOF structure.
const double CHOPPER_LOCATION[4]
const double CHOPPER_PHASE_OFFSET[2][4]
const double CHOPPER_ANGLE[4]
std::shared_ptr< const EventWorkspace > EventWorkspace_const_sptr
shared pointer to a const Workspace2D
std::shared_ptr< EventWorkspace > EventWorkspace_sptr
shared pointer to the EventWorkspace class
std::enable_if< std::is_pointer< Arg >::value, bool >::type threadSafe(Arg workspace)
Thread-safety check Checks the workspace to ensure it is suitable for multithreaded access.
Definition: MultiThreaded.h:22
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54