EventList.cpp

#include "MantidDataObjects/EventList.h"
#include "MantidDataObjects/Histogram1D.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidDataObjects/EventWorkspaceMRU.h"
#include "MantidKernel/DateAndTime.h"
#include "MantidKernel/DateAndTimeHelpers.h"
#include "MantidKernel/Exception.h"
#include "MantidKernel/Logger.h"
#include "MantidKernel/Unit.h"

#ifdef _MSC_VER
// qualifier applied to function type has no meaning; ignored
#pragma warning(disable : 4180)
#endif
#include "tbb/parallel_sort.h"
#ifdef _MSC_VER
#pragma warning(default : 4180)
#endif

#include <cfloat>
#include <cmath>
#include <functional>
#include <limits>
#include <stdexcept>

using std::ostream;
using std::runtime_error;
using std::size_t;
using std::vector;

namespace Mantid {
namespace DataObjects {
using Types::Core::DateAndTime;
using Types::Event::TofEvent;
using namespace Mantid::API;

namespace {

const double SEC_TO_NANO = 1.e9;

/**
 * Calculate the corrected full time in nanoseconds
 * @param event : The event with pulse time and time-of-flight
 * @param tofFactor : Time of flight coefficient factor
 * @param tofShift : Tof shift in seconds
 * @return Corrected full time at sample in Nanoseconds.
 */
template <typename EventType>
int64_t calculateCorrectedFullTime(const EventType &event,
                                   const double tofFactor,
                                   const double tofShift) {
  return event.pulseTime().totalNanoseconds() +
         static_cast<int64_t>(tofFactor * (event.tof() * 1.0E3) +
                              (tofShift * 1.0E9));
}

/**
 * Type for comparing events in terms of time at sample
 */
template <typename EventType>
class CompareTimeAtSample
    : public std::binary_function<EventType, EventType, bool> {
private:
  const double m_tofFactor;
  const double m_tofShift;

public:
  CompareTimeAtSample(const double tofFactor, const double tofShift)
      : m_tofFactor(tofFactor), m_tofShift(tofShift) {}

  /**
   * Compare two events based on the time they arrived at the sample.
   * Coefficient is used to provide scaling.
   * For elastic scattering coefficient is L1 / (L1 + L2)
   * @param e1 :: first event to compare
   * @param e2 :: second event to compare
   * @param coefficient :: scaling coefficient
   * @return True if first event evaluates to be < second event, otherwise false
   */
  bool operator()(const EventType &e1, const EventType &e2) const {
    const auto tAtSample1 =
        calculateCorrectedFullTime(e1, m_tofFactor, m_tofShift);
    const auto tAtSample2 =
        calculateCorrectedFullTime(e2, m_tofFactor, m_tofShift);
    return (tAtSample1 < tAtSample2);
  }
};
}
//==========================================================================
/// --------------------- TofEvent Comparators
/// ----------------------------------
//==========================================================================
/** Compare two events' TOF, return true if e1 should be before e2.
 * @param e1 :: first event
 * @param e2 :: second event
 *  */
template <typename T> bool compareEventTof(const T &e1, const T &e2) {
  return (e1.tof() < e2.tof());
}

/** Compare two events' FRAME id, return true if e1 should be before e2.
 * @param e1 :: first event
 * @param e2 :: second event
 *  */
bool compareEventPulseTime(const TofEvent &e1, const TofEvent &e2) {
  return (e1.pulseTime() < e2.pulseTime());
}

/** Compare two events' FRAME id, return true if e1 should be before e2.
 *  Assuming that if e1's pulse time is earlier than e2's, then e1 must be
 * earlier regardless TOF value
 * @param e1 :: first event
 * @param e2 :: second event
 *  */
bool compareEventPulseTimeTOF(const TofEvent &e1, const TofEvent &e2) {

  if (e1.pulseTime() < e2.pulseTime()) {
    return true;
  } else if ((e1.pulseTime() == e2.pulseTime()) && (e1.tof() < e2.tof())) {
    return true;
  }

  return false;
}

// comparator for pulse time with tolerance
struct comparePulseTimeTOFDelta {
  explicit comparePulseTimeTOFDelta(const Types::Core::DateAndTime &start,
                                    const double seconds)
      : startNano(start.totalNanoseconds()),
        deltaNano(static_cast<int64_t>(seconds * SEC_TO_NANO)) {}

  bool operator()(const TofEvent &e1, const TofEvent &e2) {
    // get the pulse times converted into bin number from start time
    const int64_t e1Pulse =
        (e1.pulseTime().totalNanoseconds() - startNano) / deltaNano;
    const int64_t e2Pulse =
        (e2.pulseTime().totalNanoseconds() - startNano) / deltaNano;

    // compare with the calculated bin information
    if (e1Pulse < e2Pulse) {
      return true;
    } else if ((e1Pulse == e2Pulse) && (e1.tof() < e2.tof())) {
      return true;
    }

    return false;
  }

  int64_t startNano;
  int64_t deltaNano;
};

/// Constructor (empty)
// EventWorkspace is always histogram data and so is thus EventList
EventList::EventList()
    : m_histogram(HistogramData::Histogram::XMode::BinEdges,
                  HistogramData::Histogram::YMode::Counts),
      eventType(TOF), order(UNSORTED), mru(nullptr) {}

/** Constructor with a MRU list
 * @param mru :: pointer to the MRU of the parent EventWorkspace
 * @param specNo :: the spectrum number for the event list
 */
EventList::EventList(EventWorkspaceMRU *mru, specnum_t specNo)
    : IEventList(specNo), m_histogram(HistogramData::Histogram::XMode::BinEdges,
                                      HistogramData::Histogram::YMode::Counts),
      eventType(TOF), order(UNSORTED), mru(mru) {}

/** Constructor copying from an existing event list
 * @param rhs :: EventList object to copy*/
EventList::EventList(const EventList &rhs)
    : IEventList(rhs), m_histogram(rhs.m_histogram), mru{nullptr} {
  // Note that operator= also assigns m_histogram, but the above use of the copy
  // constructor avoid a memory allocation and is thus faster.
  this->operator=(rhs);
}

/** Constructor, taking a vector of events.
 * @param events :: Vector of TofEvent's */
EventList::EventList(const std::vector<TofEvent> &events)
    : m_histogram(HistogramData::Histogram::XMode::BinEdges,
                  HistogramData::Histogram::YMode::Counts),
      eventType(TOF), mru(nullptr) {
  this->events.assign(events.begin(), events.end());
  this->eventType = TOF;
  this->order = UNSORTED;
}

/** Constructor, taking a vector of events.
 * @param events :: Vector of WeightedEvent's */
EventList::EventList(const std::vector<WeightedEvent> &events)
    : m_histogram(HistogramData::Histogram::XMode::BinEdges,
                  HistogramData::Histogram::YMode::Counts),
      mru(nullptr) {
  this->weightedEvents.assign(events.begin(), events.end());
  this->eventType = WEIGHTED;
  this->order = UNSORTED;
}

/** Constructor, taking a vector of events.
 * @param events :: Vector of WeightedEventNoTime's */
EventList::EventList(const std::vector<WeightedEventNoTime> &events)
    : m_histogram(HistogramData::Histogram::XMode::BinEdges,
                  HistogramData::Histogram::YMode::Counts),
      mru(nullptr) {
  this->weightedEventsNoTime.assign(events.begin(), events.end());
  this->eventType = WEIGHTED_NOTIME;
  this->order = UNSORTED;
}

/// Destructor
EventList::~EventList() {
  // Note: These two lines do not seem to have an effect on releasing memory
  //  at least on Linux. (Memory usage seems to increase event after deleting
  //  EventWorkspaces.
  //  Therefore, for performance, they are kept commented:
  clear();

  // this->events.clear();
  // std::vector<TofEvent>().swap(events); //Trick to release the vector memory.
}

/// Copy data from another EventList, via ISpectrum reference.
void EventList::copyDataFrom(const ISpectrum &source) {
  source.copyDataInto(*this);
}

/// Used by copyDataFrom for dynamic dispatch for its `source`.
void EventList::copyDataInto(EventList &sink) const {
  sink.m_histogram = m_histogram;
  sink.events = events;
  sink.weightedEvents = weightedEvents;
  sink.weightedEventsNoTime = weightedEventsNoTime;
  sink.eventType = eventType;
  sink.order = order;
}

/// Used by Histogram1D::copyDataFrom for dynamic dispatch for `other`.
void EventList::copyDataInto(Histogram1D &sink) const {
  sink.setHistogram(histogram());
}

// --------------------------------------------------------------------------
/** Create an EventList from a histogram. This converts bins to weighted
 * events.
 * Any existing events are cleared.
 *
 * @param inSpec :: ISpectrum ptr to histogram data.
 * @param GenerateZeros :: if true, generate event(s) for empty bins
 * @param GenerateMultipleEvents :: if true, create several evenly-spaced fake
 *events inside the bin
 * @param MaxEventsPerBin :: max number of events to generate in one bin, if
 *GenerateMultipleEvents
 */
void EventList::createFromHistogram(const ISpectrum *inSpec, bool GenerateZeros,
                                    bool GenerateMultipleEvents,
                                    int MaxEventsPerBin) {
  // Fresh start
  this->clear(true);

  // Get the input histogram
  const MantidVec &X = inSpec->readX();
  const MantidVec &Y = inSpec->readY();
  const MantidVec &E = inSpec->readE();
  if (Y.size() + 1 != X.size())
    throw std::runtime_error(
        "Expected a histogram (X vector should be 1 longer than the Y vector)");

  // Copy detector IDs and spectra
  this->copyInfoFrom(*inSpec);
  // We need weights but have no way to set the time. So use weighted, no time
  this->switchTo(WEIGHTED_NOTIME);
  if (GenerateZeros)
    this->weightedEventsNoTime.reserve(Y.size());

  for (size_t i = 0; i < X.size() - 1; i++) {
    double weight = Y[i];
    if ((weight != 0.0 || GenerateZeros) && std::isfinite(weight)) {
      double error = E[i];
      // Also check that the error is not a bad number
      if (std::isfinite(error)) {
        if (GenerateMultipleEvents) {
          // --------- Multiple events per bin ----------
          double errorSquared = error * error;
          // Find how many events to fake
          double val = weight / E[i];
          val *= val;
          // Convert to int with slight rounding up. This is to avoid rounding
          // errors
          int numEvents = int(val + 0.2);
          if (numEvents < 1)
            numEvents = 1;
          if (numEvents > MaxEventsPerBin)
            numEvents = MaxEventsPerBin;
          // Scale the weight and error for each
          weight /= numEvents;
          errorSquared /= numEvents;

          // Spread the TOF. e.g. 2 events = 0.25, 0.75.
          double tofStep = (X[i + 1] - X[i]) / (numEvents);
          for (size_t j = 0; j < size_t(numEvents); j++) {
            double tof = X[i] + tofStep * (0.5 + double(j));
            // Create and add the event
            // TODO: try emplace_back() here.
            weightedEventsNoTime.emplace_back(tof, weight, errorSquared);
          }
        } else {
          // --------- Single event per bin ----------
          // TOF = midpoint of the bin
          double tof = (X[i] + X[i + 1]) / 2.0;
          // Error squared is carried in the event
          double errorSquared = E[i];
          errorSquared *= errorSquared;
          // Create and add the event
          weightedEventsNoTime.emplace_back(tof, weight, errorSquared);
        }
      } // error is nont NAN or infinite
    }   // weight is non-zero, not NAN, and non-infinite
  }     // (each bin)

  // Set the X binning parameters
  this->setX(inSpec->ptrX());

  // Manually set that this is sorted by TOF, since it is. This will make it
  // "threadSafe" in other algos.
  this->setSortOrder(TOF_SORT);
}

// --------------------------------------------------------------------------
// --- Operators
// -------------------------------------------------------------------

/** Copy into this event list from another
 * @param rhs :: We will copy all the events from that into this object.
 * @return reference to this
 * */
EventList &EventList::operator=(const EventList &rhs) {
  // Note that we are NOT copying the MRU pointer.
  IEventList::operator=(rhs);
  m_histogram = rhs.m_histogram;
  events = rhs.events;
  weightedEvents = rhs.weightedEvents;
  weightedEventsNoTime = rhs.weightedEventsNoTime;
  eventType = rhs.eventType;
  order = rhs.order;
  return *this;
}

// --------------------------------------------------------------------------
/** Append an event to the histogram.
 * @param event :: TofEvent to add at the end of the list.
 * @return reference to this
 * */
EventList &EventList::operator+=(const TofEvent &event) {

  switch (this->eventType) {
  case TOF:
    // Simply push the events
    this->events.push_back(event);
    break;

  case WEIGHTED:
    this->weightedEvents.emplace_back(event);
    break;

  case WEIGHTED_NOTIME:
    this->weightedEventsNoTime.emplace_back(event);
    break;
  }

  this->order = UNSORTED;
  return *this;
}

// --------------------------------------------------------------------------
/** Append a list of events to the histogram.
 * The internal event list will switch to the required type.
 *
 * @param more_events :: A vector of events to append.
 * @return reference to this
 * */
EventList &EventList::operator+=(const std::vector<TofEvent> &more_events) {
  switch (this->eventType) {
  case TOF:
    // Simply push the events
    this->events.insert(this->events.end(), more_events.begin(),
                        more_events.end());
    break;

  case WEIGHTED:
    // Add default weights to all the un-weighted incoming events from the list.
    // and append to the list
    this->weightedEvents.reserve(this->weightedEvents.size() +
                                 more_events.size());
    for (const auto &event : more_events) {
      this->weightedEvents.emplace_back(event);
    }
    break;

  case WEIGHTED_NOTIME:
    // Add default weights to all the un-weighted incoming events from the list.
    // and append to the list
    this->weightedEventsNoTime.reserve(this->weightedEventsNoTime.size() +
                                       more_events.size());
    for (const auto &more_event : more_events)
      this->weightedEventsNoTime.emplace_back(more_event);
    break;
  }

  this->order = UNSORTED;
  return *this;
}

// --------------------------------------------------------------------------
/** Append a WeightedEvent to the histogram.
 * Note: The whole list will switch to weights (a possibly lengthy operation)
 *  if it did not have weights before.
 *
 * @param event :: WeightedEvent to add at the end of the list.
 * @return reference to this
 * */
EventList &EventList::operator+=(const WeightedEvent &event) {
  this->switchTo(WEIGHTED);
  this->weightedEvents.push_back(event);
  this->order = UNSORTED;
  return *this;
}

// --------------------------------------------------------------------------
/** Append a list of events to the histogram.
 * Note: The whole list will switch to weights (a possibly lengthy operation)
 *  if it did not have weights before.
 *
 * @param more_events :: A vector of events to append.
 * @return reference to this
 * */
EventList &EventList::
operator+=(const std::vector<WeightedEvent> &more_events) {
  switch (this->eventType) {
  case TOF:
    // Need to switch to weighted
    this->switchTo(WEIGHTED);
  // Fall through to the insertion!

  case WEIGHTED:
    // Append the two lists
    this->weightedEvents.insert(weightedEvents.end(), more_events.begin(),
                                more_events.end());
    break;

  case WEIGHTED_NOTIME:
    // Add default weights to all the un-weighted incoming events from the list.
    // and append to the list
    this->weightedEventsNoTime.reserve(this->weightedEventsNoTime.size() +
                                       more_events.size());
    for (const auto &event : more_events) {
      this->weightedEventsNoTime.emplace_back(event);
    }
    break;
  }

  this->order = UNSORTED;
  return *this;
}

// --------------------------------------------------------------------------
/** Append a list of events to the histogram.
 * Note: The whole list will switch to weights (a possibly lengthy operation)
 *  if it did not have weights before.
 *
 * @param more_events :: A vector of events to append.
 * @return reference to this
 * */
EventList &EventList::
operator+=(const std::vector<WeightedEventNoTime> &more_events) {
  switch (this->eventType) {
  case TOF:
  case WEIGHTED:
    // Need to switch to weighted with no time
    this->switchTo(WEIGHTED_NOTIME);
  // Fall through to the insertion!

  case WEIGHTED_NOTIME:
    // Simple appending of the two lists
    this->weightedEventsNoTime.insert(weightedEventsNoTime.end(),
                                      more_events.begin(), more_events.end());
    break;
  }

  this->order = UNSORTED;
  return *this;
}

// --------------------------------------------------------------------------
/** Append another EventList to this event list.
 * The event lists are concatenated, and a union of the sets of detector ID's is
 *done.
 * Switching of event types may occur if the two are different.
 *
 * @param more_events :: Another EventList.
 * @return reference to this
 * */
EventList &EventList::operator+=(const EventList &more_events) {
  // We'll let the += operator for the given vector of event lists handle it
  switch (more_events.getEventType()) {
  case TOF:
    this->operator+=(more_events.events);
    break;

  case WEIGHTED:
    this->operator+=(more_events.weightedEvents);
    break;

  case WEIGHTED_NOTIME:
    this->operator+=(more_events.weightedEventsNoTime);
    break;
  }

  // No guaranteed order
  this->order = UNSORTED;
  // Do a union between the detector IDs of both lists
  addDetectorIDs(more_events.getDetectorIDs());

  return *this;
}

// --------------------------------------------------------------------------
/** SUBTRACT another EventList from this event list.
 * The event lists are concatenated, but the weights of the incoming
 *    list are multiplied by -1.0.
 *
 * @tparam T1, T2 :: TofEvent, WeightedEvent or WeightedEventNoTime
 * @param events :: The event vector being changed.
 * @param more_events :: Another event vector being subtracted from this.
 * @return reference to this
 * */
template <class T1, class T2>
void EventList::minusHelper(std::vector<T1> &events,
                            const std::vector<T2> &more_events) {
  // Make the end vector big enough in one go (avoids repeated re-allocations).
  events.reserve(events.size() + more_events.size());
  /* In the event of subtracting in place, calling the end() vector would make
   * it point
   * at the wrong place
   * Using it caused a segault, Ticket #2306.
   * So we cache the end (this speeds up too).
   */
  for (const auto &ev : more_events) {
    // We call the constructor for T1. In the case of WeightedEventNoTime, the
    // pulse time will just be ignored.
    events.emplace_back(ev.tof(), ev.pulseTime(), ev.weight() * (-1.0),
                        ev.errorSquared());
  }
}

// --------------------------------------------------------------------------
/** SUBTRACT another EventList from this event list.
 * The event lists are concatenated, but the weights of the incoming
 *    list are multiplied by -1.0.
 *
 * @param more_events :: Another EventList.
 * @return reference to this
 * */
EventList &EventList::operator-=(const EventList &more_events) {
  if (this == &more_events) {
    // Special case, ticket #3844 part 2.
    // When doing this = this - this,
    // simply clear the input event list. Saves memory!
    this->clearData();
    return *this;
  }

  // We'll let the -= operator for the given vector of event lists handle it
  switch (this->getEventType()) {
  case TOF:
    this->switchTo(WEIGHTED);
  // Fall through

  case WEIGHTED:
    switch (more_events.getEventType()) {
    case TOF:
      minusHelper(this->weightedEvents, more_events.events);
      break;
    case WEIGHTED:
      minusHelper(this->weightedEvents, more_events.weightedEvents);
      break;
    case WEIGHTED_NOTIME:
      // TODO: Should this throw?
      minusHelper(this->weightedEvents, more_events.weightedEventsNoTime);
      break;
    }
    break;

  case WEIGHTED_NOTIME:
    switch (more_events.getEventType()) {
    case TOF:
      minusHelper(this->weightedEventsNoTime, more_events.events);
      break;
    case WEIGHTED:
      minusHelper(this->weightedEventsNoTime, more_events.weightedEvents);
      break;
    case WEIGHTED_NOTIME:
      minusHelper(this->weightedEventsNoTime, more_events.weightedEventsNoTime);
      break;
    }
    break;
  }

  // No guaranteed order
  this->order = UNSORTED;

  // NOTE: What to do about detector ID's?
  return *this;
}

// --------------------------------------------------------------------------
/** Equality operator between EventList's
 * @param rhs :: other EventList to compare
 * @return :: true if equal.
 */
bool EventList::operator==(const EventList &rhs) const {
  if (this->getNumberEvents() != rhs.getNumberEvents())
    return false;
  if (this->eventType != rhs.eventType)
    return false;
  // Check all event lists; The empty ones will compare equal
  if (events != rhs.events)
    return false;
  if (weightedEvents != rhs.weightedEvents)
    return false;
  if (weightedEventsNoTime != rhs.weightedEventsNoTime)
    return false;
  return true;
}

/** Inequality comparator
 * @param rhs :: other EventList to compare
 * @return :: true if not equal.
 */
bool EventList::operator!=(const EventList &rhs) const {
  return (!this->operator==(rhs));
}

bool EventList::equals(const EventList &rhs, const double tolTof,
                       const double tolWeight, const int64_t tolPulse) const {
  // generic checks
  if (this->getNumberEvents() != rhs.getNumberEvents())
    return false;
  if (this->eventType != rhs.eventType)
    return false;

  // loop over the events
  size_t numEvents = this->getNumberEvents();
  switch (this->eventType) {
  case TOF:
    for (size_t i = 0; i < numEvents; ++i) {
      if (!this->events[i].equals(rhs.events[i], tolTof, tolPulse))
        return false;
    }
    break;
  case WEIGHTED:
    for (size_t i = 0; i < numEvents; ++i) {
      if (!this->weightedEvents[i].equals(rhs.weightedEvents[i], tolTof,
                                          tolWeight, tolPulse))
        return false;
    }
    break;
  case WEIGHTED_NOTIME:
    for (size_t i = 0; i < numEvents; ++i) {
      if (!this->weightedEventsNoTime[i].equals(rhs.weightedEventsNoTime[i],
                                                tolTof, tolWeight))
        return false;
    }
    break;
  default:
    break;
  }

  // anything that gets this far is equal within tolerances
  return true;
}

// -----------------------------------------------------------------------------------------------
/** Return the type of Event vector contained within.
 * @return :: a EventType value.
 */
EventType EventList::getEventType() const { return eventType; }

// -----------------------------------------------------------------------------------------------
/** Switch the EventList to use the given EventType (TOF, WEIGHTED, or
 * WEIGHTED_NOTIME)
 */
void EventList::switchTo(EventType newType) {
  switch (newType) {
  case TOF:
    if (eventType != TOF)
      throw std::runtime_error("EventList::switchTo() called on an EventList "
                               "with weights to go down to TofEvent's. This "
                               "would remove weight information and therefore "
                               "is not possible.");
    break;

  case WEIGHTED:
    switchToWeightedEvents();
    break;

  case WEIGHTED_NOTIME:
    switchToWeightedEventsNoTime();
    break;
  }
  // Make sure to free memory
  this->clearUnused();
}

// -----------------------------------------------------------------------------------------------
/** Switch the EventList to use WeightedEvents instead
 * of TofEvent.
 */
void EventList::switchToWeightedEvents() {
  switch (eventType) {
  case WEIGHTED:
    // Do nothing; it already is weighted
    return;

  case WEIGHTED_NOTIME:
    throw std::runtime_error("EventList::switchToWeightedEvents() called on an "
                             "EventList with WeightedEventNoTime's. It has "
                             "lost the pulse time information and can't go "
                             "back to WeightedEvent's.");
    break;

  case TOF:
    weightedEventsNoTime.clear();
    // Convert and copy all TofEvents to the weightedEvents list.
    this->weightedEvents.assign(events.cbegin(), events.cend());
    // Get rid of the old events
    events.clear();
    eventType = WEIGHTED;
    break;
  }
}

// -----------------------------------------------------------------------------------------------
/** Switch the EventList to use WeightedEventNoTime's instead
 * of TofEvent.
 */
void EventList::switchToWeightedEventsNoTime() {
  switch (eventType) {
  case WEIGHTED_NOTIME:
    // Do nothing if already there
    return;

  case TOF: {
    // Convert and copy all TofEvents to the weightedEvents list.
    this->weightedEventsNoTime.assign(events.cbegin(), events.cend());
    // Get rid of the old events
    events.clear();
    weightedEvents.clear();
    eventType = WEIGHTED_NOTIME;
  } break;

  case WEIGHTED: {
    // Convert and copy all TofEvents to the weightedEvents list.
    this->weightedEventsNoTime.assign(weightedEvents.cbegin(),
                                      weightedEvents.cend());
    // Get rid of the old events
    events.clear();
    weightedEvents.clear();
    eventType = WEIGHTED_NOTIME;
  } break;
  }
}

// ==============================================================================================
// --- Testing functions (mostly)
// ---------------------------------------------------------------
// ==============================================================================================

/** Return the given event in the list.
 * Handles the different types of events by converting to WeightedEvent (the
 * most general type).
 * @param event_number :: the index of the event to retrieve
 * @return a WeightedEvent
 */
WeightedEvent EventList::getEvent(size_t event_number) {
  switch (eventType) {
  case TOF:
    return WeightedEvent(events[event_number]);
  case WEIGHTED:
    return weightedEvents[event_number];
  case WEIGHTED_NOTIME:
    return WeightedEvent(weightedEventsNoTime[event_number].tof(), 0,
                         weightedEventsNoTime[event_number].weight(),
                         weightedEventsNoTime[event_number].errorSquared());
  }
  throw std::runtime_error("EventList: invalid event type value was found.");
}

// ==============================================================================================
// --- Handling the event list
// -------------------------------------------------------------------
// ==============================================================================================

/** Return the const list of TofEvents contained.
 * NOTE! This should be used for testing purposes only, as much as possible. The
 *EventList
 * may contain weighted events, requiring use of getWeightedEvents() instead.
 *
 * @return a const reference to the list of non-weighted events
 * */
const std::vector<TofEvent> &EventList::getEvents() const {
  if (eventType != TOF)
    throw std::runtime_error("EventList::getEvents() called for an EventList "
                             "that has weights. Use getWeightedEvents() or "
                             "getWeightedEventsNoTime().");
  return this->events;
}

/** Return the list of TofEvents contained.
 * NOTE! This should be used for testing purposes only, as much as possible. The
 *EventList
 * may contain weighted events, requiring use of getWeightedEvents() instead.
 *
 * @return a reference to the list of non-weighted events
 * */
std::vector<TofEvent> &EventList::getEvents() {
  if (eventType != TOF)
    throw std::runtime_error("EventList::getEvents() called for an EventList "
                             "that has weights. Use getWeightedEvents() or "
                             "getWeightedEventsNoTime().");
  return this->events;
}

/** Return the list of WeightedEvent contained.
 * NOTE! This should be used for testing purposes only, as much as possible. The
 *EventList
 * may contain un-weighted events, requiring use of getEvents() instead.
 *
 * @return a reference to the list of weighted events
 * */
std::vector<WeightedEvent> &EventList::getWeightedEvents() {
  if (eventType != WEIGHTED)
    throw std::runtime_error("EventList::getWeightedEvents() called for an "
                             "EventList not of type WeightedEvent. Use "
                             "getEvents() or getWeightedEventsNoTime().");
  return this->weightedEvents;
}

/** Return the list of WeightedEvent contained.
 * NOTE! This should be used for testing purposes only, as much as possible. The
 *EventList
 * may contain un-weighted events, requiring use of getEvents() instead.
 *
 * @return a const reference to the list of weighted events
 * */
const std::vector<WeightedEvent> &EventList::getWeightedEvents() const {
  if (eventType != WEIGHTED)
    throw std::runtime_error("EventList::getWeightedEvents() called for an "
                             "EventList not of type WeightedEvent. Use "
                             "getEvents() or getWeightedEventsNoTime().");
  return this->weightedEvents;
}

/** Return the list of WeightedEvent contained.
 * NOTE! This should be used for testing purposes only, as much as possible.
 *
 * @return a reference to the list of weighted events
 * */
std::vector<WeightedEventNoTime> &EventList::getWeightedEventsNoTime() {
  if (eventType != WEIGHTED_NOTIME)
    throw std::runtime_error("EventList::getWeightedEvents() called for an "
                             "EventList not of type WeightedEventNoTime. Use "
                             "getEvents() or getWeightedEvents().");
  return this->weightedEventsNoTime;
}

/** Return the list of WeightedEventNoTime contained.
 * NOTE! This should be used for testing purposes only, as much as possible.
 *
 * @return a const reference to the list of weighted events
 * */
const std::vector<WeightedEventNoTime> &
EventList::getWeightedEventsNoTime() const {
  if (eventType != WEIGHTED_NOTIME)
    throw std::runtime_error("EventList::getWeightedEventsNoTime() called for "
                             "an EventList not of type WeightedEventNoTime. "
                             "Use getEvents() or getWeightedEvents().");
  return this->weightedEventsNoTime;
}

/** Clear the list of events and any
 * associated detector ID's.
 * */
void EventList::clear(const bool removeDetIDs) {
  if (mru)
    mru->deleteIndex(this);
  this->events.clear();
  std::vector<TofEvent>().swap(this->events); // STL Trick to release memory
  this->weightedEvents.clear();
  std::vector<WeightedEvent>().swap(
      this->weightedEvents); // STL Trick to release memory
  this->weightedEventsNoTime.clear();
  std::vector<WeightedEventNoTime>().swap(
      this->weightedEventsNoTime); // STL Trick to release memory
  if (removeDetIDs)
    this->clearDetectorIDs();
}

/** Clear any unused event lists (the ones that do not
 * match the currently used type).
 * Memory is freed.
 * */
void EventList::clearUnused() {
  if (eventType != TOF) {
    this->events.clear();
    std::vector<TofEvent>().swap(this->events); // STL Trick to release memory
  }
  if (eventType != WEIGHTED) {
    this->weightedEvents.clear();
    std::vector<WeightedEvent>().swap(
        this->weightedEvents); // STL Trick to release memory
  }
  if (eventType != WEIGHTED_NOTIME) {
    this->weightedEventsNoTime.clear();
    std::vector<WeightedEventNoTime>().swap(
        this->weightedEventsNoTime); // STL Trick to release memory
  }
}

/// Mask the spectrum to this value. Removes all events.
void EventList::clearData() { this->clear(false); }

/** Sets the MRU list for this event list
 *
 * @param newMRU :: new MRU for the workspace containing this EventList
 */
void EventList::setMRU(EventWorkspaceMRU *newMRU) { mru = newMRU; }

/** Reserve a certain number of entries in the (NOT-WEIGHTED) event list. Do NOT
 *call
 * on weighted events!
 *
 * Calls std::vector<>::reserve() in order to pre-allocate the length of the
 *event list vector.
 *
 * @param num :: number of events that will be in this EventList
 */
void EventList::reserve(size_t num) { this->events.reserve(num); }

// ==============================================================================================
// --- Sorting functions -----------------------------------------------------
// ==============================================================================================

// --------------------------------------------------------------------------
/** Sort events by TOF or Frame
 * @param order :: Order by which to sort.
 * */
void EventList::sort(const EventSortType order) const {
  if (order == UNSORTED) {
    return; // don't bother doing anything. Why did you ask to unsort?
  } else if (order == TOF_SORT) {
    this->sortTof();
  } else if (order == PULSETIME_SORT) {
    this->sortPulseTime();
  } else if (order == PULSETIMETOF_SORT) {
    this->sortPulseTimeTOF();
  } else if (order == PULSETIMETOF_DELTA_SORT) {
    throw std::invalid_argument("sorting by pulse time with delta requires "
                                "extra parameters. Use sortPulseTimeTOFDelta "
                                "instead.");
  } else if (order == TIMEATSAMPLE_SORT) {
    throw std::invalid_argument("sorting by time at sample requires extra "
                                "parameters. Use sortTimeAtSample instead.");
  } else {
    throw runtime_error("Invalid sort type in EventList::sort(EventSortType)");
  }
}

// --------------------------------------------------------------------------
/** Manually set the event list sort order value. No actual sorting takes place.
 * SHOULD ONLY BE USED IN TESTS or if you know what you are doing.
 * @param order :: sort order to set.
 */
void EventList::setSortOrder(const EventSortType order) const {
  this->order = order;
}

//  // MergeSort from:
//  http://en.literateprograms.org/Merge_sort_%28C_Plus_Plus%29#chunk%20def:merge
//  template<typename IT, typename VT> void insert(IT begin, IT end, const VT
//  &v)
//  {
//    while(begin+1!=end && *(begin+1)<v) {
//      std::swap(*begin, *(begin+1));
//      ++begin;
//    }
//    *begin=v;
//  }
//