EQSANSTofStructure.cpp

//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidAlgorithms/EQSANSTofStructure.h"
#include "MantidAPI/WorkspaceUnitValidator.h"
#include "MantidDataObjects/Events.h"
#include "MantidDataObjects/EventList.h"
#include "MantidDataObjects/EventWorkspace.h"
#include "MantidGeometry/Instrument.h"
#include "MantidKernel/TimeSeriesProperty.h"

#include <vector>

using namespace Mantid::Kernel;
using namespace Mantid::DataObjects;
using namespace Mantid::Geometry;

namespace Mantid {
namespace Algorithms {

// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(EQSANSTofStructure)

using namespace Kernel;
using namespace API;
using namespace Geometry;

EQSANSTofStructure::EQSANSTofStructure()
    : API::Algorithm(), frame_tof0(0.), flight_path_correction(false),
      low_tof_cut(0.), high_tof_cut(0.) {}

void EQSANSTofStructure::init() {
  declareProperty(make_unique<WorkspaceProperty<EventWorkspace>>(
                      "InputWorkspace", "", Direction::Input,
                      boost::make_shared<WorkspaceUnitValidator>("TOF")),
                  "Workspace to apply the TOF correction to");
  declareProperty("FlightPathCorrection", false,
                  "If True, the neutron flight path correction will be applied",
                  Kernel::Direction::Input);
  declareProperty("LowTOFCut", 0.0, "Width of the TOF margin to cut on the "
                                    "lower end of the TOF distribution of each "
                                    "frame",
                  Kernel::Direction::Input);
  declareProperty("HighTOFCut", 0.0, "Width of the TOF margin to cut on the "
                                     "upper end of the TOF distribution of "
                                     "each frame",
                  Kernel::Direction::Input);

  // Output parameters
  declareProperty("FrameSkipping", false,
                  "If True, the data was taken in frame skipping mode",
                  Kernel::Direction::Output);
  declareProperty("TofOffset", 0.0, "TOF offset that was applied to the data",
                  Kernel::Direction::Output);
  declareProperty(
      "WavelengthMin", 0.0,
      "Lower bound of the wavelength distribution of the first frame",
      Kernel::Direction::Output);
  declareProperty(
      "WavelengthMax", 0.0,
      "Upper bound of the wavelength distribution of the first frame",
      Kernel::Direction::Output);
  declareProperty(
      "WavelengthMinFrame2", 0.0,
      "Lower bound of the wavelength distribution of the second frame",
      Kernel::Direction::Output);
  declareProperty(
      "WavelengthMaxFrame2", 0.0,
      "Upper bound of the wavelength distribution of the second frame",
      Kernel::Direction::Output);
}

void EQSANSTofStructure::exec() {
  EventWorkspace_sptr inputWS = getProperty("InputWorkspace");
  flight_path_correction = getProperty("FlightPathCorrection");
  low_tof_cut = getProperty("LowTOFCut");
  high_tof_cut = getProperty("HighTOFCut");

  // Calculate the frame width
  auto frequencyLog = dynamic_cast<TimeSeriesProperty<double> *>(
      inputWS->run().getLogData("frequency"));
  if (!frequencyLog) {
    throw std::runtime_error("Frequency log not found.");
  }
  double frequency = frequencyLog->getStatistics().mean;
  double tof_frame_width = 1.0e6 / frequency;

  // Determine whether we need frame skipping or not by checking the chopper
  // speed
  bool frame_skipping = false;
  auto chopper_speedLog = dynamic_cast<TimeSeriesProperty<double> *>(
      inputWS->run().getLogData("Speed1"));
  if (!chopper_speedLog) {
    throw std::runtime_error("Chopper speed log not found.");
  }
  const double chopper_speed = chopper_speedLog->getStatistics().mean;
  if (std::fabs(chopper_speed - frequency / 2.0) < 1.0)
    frame_skipping = true;

  // Get TOF offset
  frame_tof0 = getTofOffset(inputWS, frame_skipping);

  // Calculate the frame width
  double tmp_frame_width =
      frame_skipping ? tof_frame_width * 2.0 : tof_frame_width;
  double frame_offset = 0.0;
  if (frame_tof0 >= tmp_frame_width)
    frame_offset =
        tmp_frame_width * (static_cast<int>(frame_tof0 / tmp_frame_width));

  this->execEvent(inputWS, frame_tof0, frame_offset, tof_frame_width,
                  tmp_frame_width, frame_skipping);
}

void EQSANSTofStructure::execEvent(
    Mantid::DataObjects::EventWorkspace_sptr inputWS, double threshold,
    double frame_offset, double tof_frame_width, double tmp_frame_width,
    bool frame_skipping) {
  const size_t numHists = inputWS->getNumberHistograms();
  Progress progress(this, 0.0, 1.0, numHists);

  // Get the nominal sample-to-detector distance (in mm)
  Mantid::Kernel::Property *prop =
      inputWS->run().getProperty("sample_detector_distance");
  auto dp = dynamic_cast<Mantid::Kernel::PropertyWithValue<double> *>(prop);
  if (!dp) {
    throw std::runtime_error("sample_detector_distance log not found.");
  }
  const double SDD = *dp / 1000.0;

  // Loop through the spectra and apply correction
  PARALLEL_FOR1(inputWS)
  for (int64_t ispec = 0; ispec < int64_t(numHists); ++ispec) {

    IDetector_const_sptr det;
    try {
      det = inputWS->getDetector(ispec);
    } catch (Exception::NotFoundError &) {
      g_log.warning() << "Workspace index " << ispec
                      << " has no detector assigned to it - discarding\n";
      // 'continue' statement moved outside catch block because Mac Intel
      // compiler has a problem with it being here in an openmp block.
    }
    if (!det)
      continue;

    // Get the flight path from the sample to the detector pixel
    const V3D samplePos = inputWS->getInstrument()->getSample()->getPos();
    const V3D scattered_flight_path = det->getPos() - samplePos;

    // Sample-to-source distance
    const V3D sourcePos = inputWS->getInstrument()->getSource()->getPos();
    const V3D SSD = samplePos - sourcePos;
    double tof_factor =
        (SSD.norm() + scattered_flight_path.norm()) / (SSD.norm() + SDD);

    PARALLEL_START_INTERUPT_REGION

    // Get the pointer to the output event list
    std::vector<TofEvent> &events = inputWS->getSpectrum(ispec).getEvents();
    std::vector<TofEvent>::iterator it;
    std::vector<TofEvent> clean_events;

    for (it = events.begin(); it < events.end(); ++it) {
      double newtof = it->tof();
      newtof += frame_offset;
      // Correct for the scattered neutron flight path
      if (flight_path_correction)
        newtof /= tof_factor;

      while (newtof < threshold)
        newtof += tmp_frame_width;

      // Remove events that don't fall within the accepted time window
      double rel_tof = newtof - frame_tof0;
      double x = (static_cast<int>(floor(rel_tof * 10)) %
                  static_cast<int>(floor(tof_frame_width * 10))) *
                 0.1;
      if (x < low_tof_cut || x > tof_frame_width - high_tof_cut) {
        continue;
      }
      // At this point the events in the second frame are still off by a frame
      if (frame_skipping && rel_tof > tof_frame_width)
        newtof += tof_frame_width;
      clean_events.emplace_back(newtof, it->pulseTime());
    }
    events.clear();
    events.reserve(clean_events.size());
    for (it = clean_events.begin(); it < clean_events.end(); ++it) {
      events.push_back(*it);
    }

    progress.report("TOF structure");
    PARALLEL_END_INTERUPT_REGION
  }
  PARALLEL_CHECK_INTERUPT_REGION
}

double EQSANSTofStructure::getTofOffset(EventWorkspace_const_sptr inputWS,
                                        bool frame_skipping) {
  //# Storage for chopper information read from the logs
  double chopper_set_phase[4] = {0, 0, 0, 0};
  double chopper_speed[4] = {0, 0, 0, 0};
  double chopper_actual_phase[4] = {0, 0, 0, 0};
  double chopper_wl_1[4] = {0, 0, 0, 0};
  double chopper_wl_2[4] = {0, 0, 0, 0};
  double frame_wl_1 = 0;
  double frame_srcpulse_wl_1 = 0;
  double frame_wl_2 = 0;
  double chopper_srcpulse_wl_1[4] = {0, 0, 0, 0};
  double chopper_frameskip_wl_1[4] = {0, 0, 0, 0};
  double chopper_frameskip_wl_2[4] = {0, 0, 0, 0};
  double chopper_frameskip_srcpulse_wl_1[4] = {0, 0, 0, 0};

  // Calculate the frame width
  auto frequencyLog = dynamic_cast<TimeSeriesProperty<double> *>(
      inputWS->run().getLogData("frequency"));
  if (!frequencyLog) {
    throw std::runtime_error("Frequency log not found.");
  }
  double frequency = frequencyLog->getStatistics().mean;
  double tof_frame_width = 1.0e6 / frequency;

  double tmp_frame_width = tof_frame_width;
  if (frame_skipping)
    tmp_frame_width *= 2.0;

  // Choice of parameter set
  int m_set = 0;
  if (frame_skipping)
    m_set = 1;

  bool first = true;
  bool first_skip = true;
  double frameskip_wl_1 = 0;
  double frameskip_srcpulse_wl_1 = 0;
  double frameskip_wl_2 = 0;

  for (int i = 0; i < 4; i++) {
    // Read chopper information
    std::ostringstream phase_str;
    phase_str << "Phase" << i + 1;
    auto log = dynamic_cast<TimeSeriesProperty<double> *>(
        inputWS->run().getLogData(phase_str.str()));
    if (!log) {
      throw std::runtime_error("Phase log not found.");
    }
    chopper_set_phase[i] = log->getStatistics().mean;
    std::ostringstream speed_str;
    speed_str << "Speed" << i + 1;
    log = dynamic_cast<TimeSeriesProperty<double> *>(
        inputWS->run().getLogData(speed_str.str()));
    if (!log) {
      throw std::runtime_error("Speed log not found.");
    }
    chopper_speed[i] = log->getStatistics().mean;

    // Only process choppers with non-zero speed
    if (chopper_speed[i] <= 0)
      continue;

    chopper_actual_phase[i] =
        chopper_set_phase[i] - CHOPPER_PHASE_OFFSET[m_set][i];

    while (chopper_actual_phase[i] < 0)
      chopper_actual_phase[i] += tmp_frame_width;

    double x1 =
        (chopper_actual_phase[i] -
         (tmp_frame_width * 0.5 * CHOPPER_ANGLE[i] / 360.)); // opening edge
    double x2 =
        (chopper_actual_phase[i] +
         (tmp_frame_width * 0.5 * CHOPPER_ANGLE[i] / 360.)); // closing edge
    if (!frame_skipping)                                     // not skipping
    {
      while (x1 < 0) {
        x1 += tmp_frame_width;
        x2 += tmp_frame_width;
      }
    }

    if (x1 > 0) {
      chopper_wl_1[i] = 3.9560346 * x1 / CHOPPER_LOCATION[i];
      chopper_srcpulse_wl_1[i] =
          3.9560346 * (x1 - chopper_wl_1[i] * PULSEWIDTH) / CHOPPER_LOCATION[i];
    } else
      chopper_wl_1[i] = chopper_srcpulse_wl_1[i] = 0.;

    if (x2 > 0)
      chopper_wl_2[i] = 3.9560346 * x2 / CHOPPER_LOCATION[i];
    else
      chopper_wl_2[i] = 0.;

    if (first) {
      frame_wl_1 = chopper_wl_1[i];
      frame_srcpulse_wl_1 = chopper_srcpulse_wl_1[i];
      frame_wl_2 = chopper_wl_2[i];
      first = false;
    } else {
      if (frame_skipping &&
          i == 2) // ignore chopper 1 and 2 forthe shortest wl.
      {
        frame_wl_1 = chopper_wl_1[i];
        frame_srcpulse_wl_1 = chopper_srcpulse_wl_1[i];
      }
      if (frame_wl_1 < chopper_wl_1[i])
        frame_wl_1 = chopper_wl_1[i];
      if (frame_wl_2 > chopper_wl_2[i])
        frame_wl_2 = chopper_wl_2[i];
      if (frame_srcpulse_wl_1 < chopper_srcpulse_wl_1[i])
        frame_srcpulse_wl_1 = chopper_srcpulse_wl_1[i];
    }

    if (frame_skipping) {
      if (x1 > 0) {
        x1 += tof_frame_width; // skipped pulse
        chopper_frameskip_wl_1[i] = 3.9560346 * x1 / CHOPPER_LOCATION[i];
        chopper_frameskip_srcpulse_wl_1[i] =
            3.9560346 * (x1 - chopper_wl_1[i] * PULSEWIDTH) /
            CHOPPER_LOCATION[i];
      } else
        chopper_wl_1[i] = chopper_srcpulse_wl_1[i] = 0.;

      if (x2 > 0) {
        x2 += tof_frame_width;
        chopper_frameskip_wl_2[i] = 3.9560346 * x2 / CHOPPER_LOCATION[i];
      } else
        chopper_wl_2[i] = 0.;

      if (i < 2 && chopper_frameskip_wl_1[i] > chopper_frameskip_wl_2[i])
        continue;

      if (first_skip) {
        frameskip_wl_1 = chopper_frameskip_wl_1[i];
        frameskip_srcpulse_wl_1 = chopper_frameskip_srcpulse_wl_1[i];
        frameskip_wl_2 = chopper_frameskip_wl_2[i];
        first_skip = false;
      } else {
        if (i == 2) // ignore chopper 1 and 2 forthe longest wl.
          frameskip_wl_2 = chopper_frameskip_wl_2[i];

        if (chopper_frameskip_wl_1[i] < chopper_frameskip_wl_2[i] &&
            frameskip_wl_1 < chopper_frameskip_wl_1[i])
          frameskip_wl_1 = chopper_frameskip_wl_1[i];

        if (chopper_frameskip_wl_1[i] < chopper_frameskip_wl_2[i] &&
            frameskip_srcpulse_wl_1 < chopper_frameskip_srcpulse_wl_1[i])
          frameskip_srcpulse_wl_1 = chopper_frameskip_srcpulse_wl_1[i];

        if (frameskip_wl_2 > chopper_frameskip_wl_2[i])
          frameskip_wl_2 = chopper_frameskip_wl_2[i];
      }
    }
  }

  if (frame_wl_1 >= frame_wl_2) // too many frames later. So figure it out
  {
    double n_frame[4] = {0, 0, 0, 0};
    double c_wl_1[4] = {0, 0, 0, 0};
    double c_wl_2[4] = {0, 0, 0, 0};
    bool passed = false;

    do {
      frame_wl_1 = c_wl_1[0] =
          chopper_wl_1[0] +
          3.9560346 * n_frame[0] * tof_frame_width / CHOPPER_LOCATION[0];
      frame_wl_2 = c_wl_2[0] =
          chopper_wl_2[0] +
          3.9560346 * n_frame[0] * tof_frame_width / CHOPPER_LOCATION[0];

      for (int i = 1; i < 4; i++) {
        n_frame[i] = n_frame[i - 1] - 1;
        passed = false;

        do {
          n_frame[i] += 1;
          c_wl_1[i] =
              chopper_wl_1[i] +
              3.9560346 * n_frame[i] * tof_frame_width / CHOPPER_LOCATION[i];
          c_wl_2[i] =
              chopper_wl_2[i] +
              3.9560346 * n_frame[i] * tof_frame_width / CHOPPER_LOCATION[i];

          if (frame_wl_1 < c_wl_2[i] && frame_wl_2 > c_wl_1[i]) {
            passed = true;
            break;
          }
          if (frame_wl_2 < c_wl_1[i])
            break; // over shot
        } while (n_frame[i] - n_frame[i - 1] < 10);

        if (!passed) {
          n_frame[0] += 1;
          break;
        } else {
          if (frame_wl_1 < c_wl_1[i])
            frame_wl_1 = c_wl_1[i];
          if (frame_wl_2 > c_wl_2[i])
            frame_wl_2 = c_wl_2[i];
        }
      }
    } while (!passed && n_frame[0] < 99);

    if (frame_wl_2 > frame_wl_1) {
      int n = 3;
      if (c_wl_1[2] > c_wl_1[3])
        n = 2;

      frame_srcpulse_wl_1 =
          c_wl_1[n] - 3.9560346 * c_wl_1[n] * PULSEWIDTH / CHOPPER_LOCATION[n];

      for (int i = 0; i < 4; i++) {
        chopper_wl_1[i] = c_wl_1[i];
        chopper_wl_2[i] = c_wl_2[i];
        if (frame_skipping) {
          chopper_frameskip_wl_1[i] =
              c_wl_1[i] +
              3.9560346 * 2. * tof_frame_width / CHOPPER_LOCATION[i];
          chopper_frameskip_wl_2[i] =
              c_wl_2[i] +
              3.9560346 * 2. * tof_frame_width / CHOPPER_LOCATION[i];
          if (i == 0) {
            frameskip_wl_1 = chopper_frameskip_wl_1[i];
            frameskip_wl_2 = chopper_frameskip_wl_2[i];
          } else {
            if (frameskip_wl_1 < chopper_frameskip_wl_1[i])
              frameskip_wl_1 = chopper_frameskip_wl_1[i];
            if (frameskip_wl_2 > chopper_frameskip_wl_2[i])
              frameskip_wl_2 = chopper_frameskip_wl_2[i];
          }
        }
      }
    } else
      frame_srcpulse_wl_1 = 0.0;
  }
  // Get source and detector locations
  // get the name of the mapping file as set in the parameter files
  std::vector<std::string> temp =
      inputWS->getInstrument()->getStringParameter("detector-name");
  std::string det_name = "detector1";
  if (temp.empty())
    g_log.information() << "The instrument parameter file does not contain the "
                           "'detector-name' parameter: trying 'detector1'";
  else
    det_name = temp[0];

  double source_z = inputWS->getInstrument()->getSource()->getPos().Z();
  double detector_z =
      inputWS->getInstrument()->getComponentByName(det_name)->getPos().Z();

  double source_to_detector = (detector_z - source_z) * 1000.0;
  frame_tof0 = frame_srcpulse_wl_1 / 3.9560346 * source_to_detector;

  g_log.information() << "Frame width " << tmp_frame_width << '\n';
  g_log.information() << "TOF offset = " << frame_tof0 << " microseconds\n";
  g_log.information() << "Band defined by T1-T4 " << frame_wl_1 << " "
                      << frame_wl_2;
  if (frame_skipping)
    g_log.information() << " + " << frameskip_wl_1 << " " << frameskip_wl_2
                        << '\n';
  else
    g_log.information() << '\n';
  g_log.information() << "Chopper    Actual Phase    Lambda1    Lambda2\n";
  for (int i = 0; i < 4; i++)
    g_log.information() << i << "    " << chopper_actual_phase[i] << "  "
                        << chopper_wl_1[i] << "  " << chopper_wl_2[i] << '\n';

  double low_wl_discard = 3.9560346 * low_tof_cut / source_to_detector;
  double high_wl_discard = 3.9560346 * high_tof_cut / source_to_detector;

  setProperty("FrameSkipping", frame_skipping);
  setProperty("TofOffset", frame_tof0);
  setProperty("WavelengthMin", frame_wl_1 + low_wl_discard);
  setProperty("WavelengthMax", frame_wl_2 - high_wl_discard);
  if (frame_skipping) {
    setProperty("WavelengthMinFrame2", frameskip_wl_1 + low_wl_discard);
    setProperty("WavelengthMaxFrame2", frameskip_wl_2 - high_wl_discard);
  }

  return frame_tof0;
}

} // namespace Algorithms
} // namespace Mantid