23using namespace Kernel;
25using namespace Geometry;
29 std::make_shared<WorkspaceUnitValidator>(
"TOF")),
30 "Workspace to apply the TOF correction to");
34 "Workspace to store the corrected data in");
35 declareProperty(
"FrameSkipping",
false,
"True if the data was taken in frame-skipping mode",
44 if (outputWS != inputWS) {
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";
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";
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;
71 auto log = inputWS->run().getTimeSeriesProperty<
double>(
"frequency");
72 double frequency = log->getStatistics().mean;
73 double tof_frame_width = 1.0e6 / frequency;
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;
86 double frame_tof0 =
getTofOffset(inputWS, frame_skipping, source_to_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));
107 const auto nTOF =
static_cast<int>(XIn.size());
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)
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";
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];
153 for (
int i = tof_bin_range - 1; i < nTOF - 1; i++) {
154 XOut[i] = XOut[i - 1] + 10.0;
158 XOut[nTOF - 1] = XOut[nTOF - 2] + 10.0;
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];
167 XOut[tof_bin_range - 2 - cutoff] = XIn[tof_bin_range] + frame_offset;
171 YOut[tof_bin_range - 2 - cutoff] = 0.0;
172 EOut[tof_bin_range - 2 - cutoff] = 0.0;
178 double source_to_monitor) {
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};
194 auto log = inputWS->run().getTimeSeriesProperty<
double>(
"frequency");
195 double frequency = log->getStatistics().mean;
196 double tof_frame_width = 1.0e6 / frequency;
198 double tmp_frame_width = tof_frame_width;
200 tmp_frame_width *= 2.0;
208 bool first_skip =
true;
209 double frameskip_wl_1 = 0;
210 double frameskip_srcpulse_wl_1 = 0;
211 double frameskip_wl_2 = 0;
213 for (
int i = 0; i < 4; i++) {
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;
225 if (chopper_speed[i] <= 0)
230 while (chopper_actual_phase[i] < 0)
231 chopper_actual_phase[i] += tmp_frame_width;
233 double x1 = (chopper_actual_phase[i] - (tmp_frame_width * 0.5 *
CHOPPER_ANGLE[i] / 360.));
234 double x2 = (chopper_actual_phase[i] + (tmp_frame_width * 0.5 *
CHOPPER_ANGLE[i] / 360.));
238 x1 += tmp_frame_width;
239 x2 += tmp_frame_width;
247 chopper_wl_1[i] = chopper_srcpulse_wl_1[i] = 0.;
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];
257 if (frame_skipping && i == 2)
259 frame_wl_1 = chopper_wl_1[i];
260 frame_srcpulse_wl_1 = chopper_srcpulse_wl_1[i];
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];
270 if (frame_skipping) {
272 x1 += tof_frame_width;
276 chopper_wl_1[i] = chopper_srcpulse_wl_1[i] = 0.;
279 x2 += tof_frame_width;
282 chopper_wl_2[i] = 0.;
284 if (i < 2 && chopper_frameskip_wl_1[i] > chopper_frameskip_wl_2[i])
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];
294 frameskip_wl_2 = chopper_frameskip_wl_2[i];
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];
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];
303 if (frameskip_wl_2 > chopper_frameskip_wl_2[i])
304 frameskip_wl_2 = chopper_frameskip_wl_2[i];
309 if (frame_wl_1 >= frame_wl_2)
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};
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];
320 for (
int i = 1; i < 4; i++) {
321 n_frame[i] = n_frame[i - 1] - 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];
329 if (frame_wl_1 < c_wl_2[i] && frame_wl_2 > c_wl_1[i]) {
333 if (frame_wl_2 < c_wl_1[i])
335 }
while (n_frame[i] - n_frame[i - 1] < 10);
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];
347 }
while (!passed && n_frame[0] < 99);
349 if (frame_wl_2 > frame_wl_1) {
351 if (c_wl_1[2] > c_wl_1[3])
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];
363 frameskip_wl_1 = chopper_frameskip_wl_1[i];
364 frameskip_wl_2 = chopper_frameskip_wl_2[i];
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];
374 frame_srcpulse_wl_1 = 0.0;
377 double frame_tof0 = frame_srcpulse_wl_1 / 3.9560346 * source_to_monitor;
381 g_log.
information() <<
"Band defined by T1-T4 " << frame_wl_1 <<
" " << frame_wl_2;
383 g_log.
information() <<
" + " << frameskip_wl_1 <<
" " << frameskip_wl_2 <<
'\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]
#define DECLARE_ALGORITHM(classname)
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
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.
void information(const std::string &msg)
Logs at information level.
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
void exec() override
Execution code.
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]
const double CHOPPER_ANGLE[4]
std::vector< double > MantidVec
typedef for the data storage used in Mantid matrix workspaces
@ Input
An input workspace.
@ Output
An output workspace.