24using namespace Kernel;
26using namespace DataObjects;
27using namespace Geometry;
28using Types::Event::TofEvent;
31 :
API::
Algorithm(), frame_tof0(0.), flight_path_correction(false), low_tof_cut(0.), high_tof_cut(0.) {}
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",
40 "Width of the TOF margin to cut on the "
41 "lower end of the TOF distribution of each "
45 "Width of the TOF margin to cut on the "
46 "upper end of the TOF distribution of "
51 declareProperty(
"FrameSkipping",
false,
"If True, the data was taken in frame skipping mode",
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",
73 throw std::runtime_error(
"Frequency log not found.");
75 double frequency = frequencyLog->getStatistics().mean;
76 double tof_frame_width = 1.0e6 / frequency;
80 bool frame_skipping =
false;
82 if (!chopper_speedLog) {
83 throw std::runtime_error(
"Chopper speed log not found.");
85 const double chopper_speed = chopper_speedLog->getStatistics().mean;
86 if (std::fabs(chopper_speed - frequency / 2.0) < 1.0)
87 frame_skipping =
true;
93 double tmp_frame_width = frame_skipping ? tof_frame_width * 2.0 : tof_frame_width;
94 double frame_offset = 0.0;
96 frame_offset = tmp_frame_width * (
static_cast<int>(
frame_tof0 / tmp_frame_width));
98 this->
execEvent(inputWS,
frame_tof0, frame_offset, tof_frame_width, tmp_frame_width, frame_skipping);
102 double frame_offset,
double tof_frame_width,
double tmp_frame_width,
103 bool frame_skipping) {
104 const size_t numHists = inputWS->getNumberHistograms();
112 throw std::runtime_error(
"sampleflange_detector_distance log not found.");
114 const double SFDD = *dp / 1000.0;
116 const auto &spectrumInfo = inputWS->spectrumInfo();
117 const auto l1 = spectrumInfo.l1();
121 for (int64_t ispec = 0; ispec < int64_t(numHists); ++ispec) {
124 if (!spectrumInfo.hasDetectors(ispec)) {
125 g_log.
warning() <<
"Workspace index " << ispec <<
" has no detector assigned to it - discarding\n";
128 const auto l2 = spectrumInfo.l2(ispec);
129 double tof_factor = (l1 +
l2) / (l1 + SFDD);
132 std::vector<TofEvent> &events = inputWS->getSpectrum(ispec).getEvents();
133 std::vector<TofEvent>::iterator it;
134 std::vector<TofEvent> clean_events;
136 for (it = events.begin(); it < events.end(); ++it) {
137 double newtof = it->tof();
138 newtof += frame_offset;
141 newtof /= tof_factor;
143 while (newtof < threshold)
144 newtof += tmp_frame_width;
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) {
153 if (frame_skipping && rel_tof > tof_frame_width)
154 newtof += tof_frame_width;
155 clean_events.emplace_back(newtof, it->pulseTime());
158 events.reserve(clean_events.size());
159 for (it = clean_events.begin(); it < clean_events.end(); ++it) {
160 events.emplace_back(*it);
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};
187 throw std::runtime_error(
"Frequency log not found.");
189 double frequency = frequencyLog->getStatistics().mean;
190 double tof_frame_width = 1.0e6 / frequency;
192 double tmp_frame_width = tof_frame_width;
194 tmp_frame_width *= 2.0;
202 bool first_skip =
true;
203 double frameskip_wl_1 = 0;
204 double frameskip_srcpulse_wl_1 = 0;
205 double frameskip_wl_2 = 0;
207 for (
int i = 0; i < 4; i++) {
209 std::ostringstream phase_str;
210 phase_str <<
"Phase" << i + 1;
213 throw std::runtime_error(
"Phase log not found.");
215 chopper_set_phase[i] = log->getStatistics().mean;
216 std::ostringstream speed_str;
217 speed_str <<
"Speed" << i + 1;
220 throw std::runtime_error(
"Speed log not found.");
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 chopper_wl_2[i] = 0.;
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];
260 if (frame_skipping && i == 2)
262 frame_wl_1 = chopper_wl_1[i];
263 frame_srcpulse_wl_1 = chopper_srcpulse_wl_1[i];
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];
273 if (frame_skipping) {
275 x1 += tof_frame_width;
279 chopper_wl_1[i] = chopper_srcpulse_wl_1[i] = 0.;
282 x2 += tof_frame_width;
285 chopper_wl_2[i] = 0.;
287 if (i < 2 && chopper_frameskip_wl_1[i] > chopper_frameskip_wl_2[i])
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];
297 frameskip_wl_2 = chopper_frameskip_wl_2[i];
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];
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];
306 if (frameskip_wl_2 > chopper_frameskip_wl_2[i])
307 frameskip_wl_2 = chopper_frameskip_wl_2[i];
312 if (frame_wl_1 >= frame_wl_2)
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};
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];
323 for (
int i = 1; i < 4; i++) {
324 n_frame[i] = n_frame[i - 1] - 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];
332 if (frame_wl_1 < c_wl_2[i] && frame_wl_2 > c_wl_1[i]) {
336 if (frame_wl_2 < c_wl_1[i])
338 }
while (n_frame[i] - n_frame[i - 1] < 10);
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];
350 }
while (!passed && n_frame[0] < 99);
352 if (frame_wl_2 > frame_wl_1) {
354 if (c_wl_1[2] > c_wl_1[3])
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];
366 frameskip_wl_1 = chopper_frameskip_wl_1[i];
367 frameskip_wl_2 = chopper_frameskip_wl_2[i];
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];
377 frame_srcpulse_wl_1 = 0.0;
381 std::vector<std::string> temp = inputWS->getInstrument()->getStringParameter(
"detector-name");
382 std::string det_name =
"detector1";
384 g_log.
information() <<
"The instrument parameter file does not contain the "
385 "'detector-name' parameter: trying 'detector1'";
391 double source_z = inputWS->getInstrument()->getSource()->getPos().Z();
392 double detector_z = inputWS->getInstrument()->getComponentByName(det_name)->getPos().Z();
394 double source_to_detector = (detector_z - source_z) * 1000.0;
395 frame_tof0 = frame_srcpulse_wl_1 / 3.9560346 * source_to_detector;
399 g_log.
information() <<
"Band defined by T1-T4 " << frame_wl_1 <<
" " << frame_wl_2;
401 g_log.
information() <<
" + " << frameskip_wl_1 <<
" " << frameskip_wl_2 <<
'\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]
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;
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);
#define DECLARE_ALGORITHM(classname)
#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...
#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.
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.
void progress(double p, const std::string &msg="", double estimatedTime=0.0, int progressPrecision=0)
Sends ProgressNotification.
Helper class for reporting progress from algorithms.
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.
bool flight_path_correction
EQSANSTofStructure()
Default constructor.
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.
void information(const std::string &msg)
Logs at information level.
The concrete, templated class for properties.
Base class for properties.
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.
@ Input
An input workspace.
@ Output
An output workspace.