diff --git a/Framework/DataObjects/inc/MantidDataObjects/EventList.h b/Framework/DataObjects/inc/MantidDataObjects/EventList.h
index 8468f490acd264b76f304bfc51a624c197977e58..ddccb5d7f1319f39c76303c3b711ae32b3a8f4bf 100644
--- a/Framework/DataObjects/inc/MantidDataObjects/EventList.h
+++ b/Framework/DataObjects/inc/MantidDataObjects/EventList.h
@@ -222,6 +222,9 @@ public:
virtual size_t histogram_size() const;
void compressEvents(double tolerance, EventList *destination);
+ void compressFatEvents(const double tolerance,
+ const Types::Core::DateAndTime &timeStart,
+ const double seconds, EventList *destination);
// get EventType declaration
void generateHistogram(const MantidVec &X, MantidVec &Y, MantidVec &E,
bool skipError = false) const override;
@@ -436,6 +439,12 @@ private:
std::vector<WeightedEventNoTime> &out,
double tolerance);
template <class T>
+ static void compressFatEventsHelper(
+ const std::vector<T> &events, std::vector<WeightedEvent> &out,
+ const double tolerance, const Mantid::Types::Core::DateAndTime &timeStart,
+ const double seconds);
+
+ template <class T>
static void histogramForWeightsHelper(const std::vector<T> &events,
const MantidVec &X, MantidVec &Y,
MantidVec &E);
diff --git a/Framework/DataObjects/src/EventList.cpp b/Framework/DataObjects/src/EventList.cpp
index 917b57063dc0a12a74fd1c72b9ed9c7bd5f31a92..c615f3362ffaa59fcafe12ba86154263ecb76bde 100644
--- a/Framework/DataObjects/src/EventList.cpp
+++ b/Framework/DataObjects/src/EventList.cpp
@@ -2,6 +2,7 @@
#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"
@@ -28,7 +29,6 @@ using std::vector;
namespace Mantid {
namespace DataObjects {
-using Kernel::Exception::NotImplementedError;
using Types::Core::DateAndTime;
using Types::Event::TofEvent;
using namespace Mantid::API;
@@ -1567,6 +1567,87 @@ void EventList::compressEventsParallelHelper(
out.insert(out.end(), outputs[thread].begin(), outputs[thread].end());
}
+template <class T>
+inline void EventList::compressFatEventsHelper(
+ const std::vector<T> &events, std::vector<WeightedEvent> &out,
+ const double tolerance, const Types::Core::DateAndTime &timeStart,
+ const double seconds) {
+ // Clear the output. We can't know ahead of time how much space to reserve :(
+ out.clear();
+ // We will make a starting guess of 1/20th of the number of input events.
+ out.reserve(events.size() / 20);
+
+ // The last TOF to which we are comparing.
+ double lastTof = std::numeric_limits<double>::lowest();
+ // For getting an accurate average TOF
+ double totalTof = 0;
+
+ // pulsetime bin information - stored as int nanoseconds because it
+ // implementation type for DateAndTime object
+ const double SEC_TO_NANO = 1.e9;
+ const int64_t PULSETIME_TOL = static_cast<int64_t>(seconds * SEC_TO_NANO);
+
+ // pusletime information
+ DateAndTime timeRight = timeStart + PULSETIME_TOL; // right boundary <
+ std::vector<DateAndTime> pulsetimes; // all the times for new event
+
+ // Carrying weight and error
+ double weight = 0;
+ double errorSquared = 0;
+
+ // move up to first event that has a large enough pusletime
+ auto it = events.cbegin();
+ for (; it != events.cend(); ++it) {
+ if (it->m_pulsetime >= timeStart)
+ break;
+ }
+ // loop through events and accumulate weight
+ for (; it != events.cend(); ++it) {
+ if ((it->m_pulsetime < timeRight) && (it->m_tof - lastTof) <= tolerance) {
+ // Carry the error and weight
+ weight += it->weight();
+ errorSquared += it->errorSquared();
+ // Track the average tof
+ totalTof += it->m_tof;
+ // Accumulate the pulse times
+ pulsetimes.push_back(it->m_pulsetime);
+ } else {
+ // We exceeded the tolerance
+ if (!pulsetimes.empty()) {
+ // Create a new event with the average TOF and summed weights and
+ // squared errors.
+ out.emplace_back(totalTof / static_cast<double>(pulsetimes.size()),
+ Kernel::DateAndTimeHelpers::averageSorted(pulsetimes),
+ weight, errorSquared);
+ }
+ // Start a new combined object
+ totalTof = it->m_tof;
+ weight = it->weight();
+ errorSquared = it->errorSquared();
+ lastTof = it->m_tof;
+ timeRight += PULSETIME_TOL;
+ pulsetimes.clear();
+ pulsetimes.push_back(it->m_pulsetime);
+ }
+ }
+
+ // Put the last event in there too.
+ if (!pulsetimes.empty()) {
+ // Create a new event with the average TOF and summed weights and squared
+ // errors.
+ out.emplace_back(totalTof / static_cast<double>(pulsetimes.size()),
+ Kernel::DateAndTimeHelpers::averageSorted(pulsetimes),
+ weight, errorSquared);
+ }
+
+ // If you have over-allocated by more than 5%, reduce the size.
+ size_t excess_limit = out.size() / 20;
+ if ((out.capacity() - out.size()) > excess_limit) {
+ // Note: This forces a copy!
+ std::vector<WeightedEvent>(out).swap(out);
+ }
+}
+
// --------------------------------------------------------------------------
/** Compress the event list by grouping events with the same
* TOF (within a given tolerance). PulseTime is ignored.
@@ -1628,6 +1709,41 @@ void EventList::compressEvents(double tolerance, EventList *destination) {
destination->clearUnused();
}
+void EventList::compressFatEvents(
+ const double tolerance, const Mantid::Types::Core::DateAndTime &timeStart,
+ const double seconds, EventList *destination) {
+ switch (eventType) {
+ case WEIGHTED_NOTIME:
+ throw std::invalid_argument(
+ "Cannot compress events that do not have pulsetime");
+ case TOF:
+ this->sortPulseTimeTOF();
+ compressFatEventsHelper(this->events, destination->weightedEvents,
+ tolerance, timeStart, seconds);
+ break;
+ case WEIGHTED:
+ this->sortPulseTimeTOF();
+ if (destination == this) {
+ // Put results in a temp output
+ std::vector<WeightedEvent> out;
+ compressFatEventsHelper(this->weightedEvents, out, tolerance, timeStart,
+ seconds);
+ // Put it back
+ this->weightedEvents.swap(out);
+ } else {
+ compressFatEventsHelper(this->weightedEvents, destination->weightedEvents,
+ tolerance, timeStart, seconds);
+ }
+ break;
+ }
+ // In all cases, you end up WEIGHTED_NOTIME.
+ destination->eventType = WEIGHTED;
+ // The sort is still valid!
+ destination->order = PULSETIMETOF_SORT;
+ // Empty out storage for vectors that are now unused.
+ destination->clearUnused();
+}
+
// --------------------------------------------------------------------------
/** Utility function:
* Returns the iterator into events of the first TofEvent with
diff --git a/Framework/DataObjects/test/EventListTest.h b/Framework/DataObjects/test/EventListTest.h
index 93b909c7efff214949b1ce4af59d0d3736aa8550..233b33800764ee0d802f4578d2fe61361120acaf 100644
--- a/Framework/DataObjects/test/EventListTest.h
+++ b/Framework/DataObjects/test/EventListTest.h
@@ -1417,18 +1417,19 @@ public:
// Go through each possible EventType (except the no-time one) as the input
for (int this_type = 0; this_type < 3; this_type++) {
EventType curType = static_cast<EventType>(this_type);
+ // Since input data has varying pulse time, and same TOF on all events.
+ // Nothing interesting happens here.
+ if (curType == WEIGHTED_NOTIME) {
+ continue;
+ }
+
EventList el = this->fake_uniform_pulse_data(curType, 1);
el.switchTo(curType);
const double tofFactor = 1; // L1 / (L1 + L2)
const double tofShift = 0;
- if (curType == WEIGHTED_NOTIME) {
- TS_ASSERT_THROWS_NOTHING(el.sortTimeAtSample(tofFactor, tofShift));
- // Since input data has varying pulse time, and same TOF on all events.
- // Nothing interesting happens here.
- } else {
- TS_ASSERT_THROWS_NOTHING(el.sortTimeAtSample(tofFactor, tofShift));
-
+ TS_ASSERT_THROWS_NOTHING(el.sortTimeAtSample(tofFactor, tofShift));
+ if (curType != WEIGHTED_NOTIME) {
for (size_t i = 1; i < el.getNumberEvents(); i++) {
TSM_ASSERT_LESS_THAN_EQUALS(this_type, el.getEvent(i - 1).pulseTime(),
el.getEvent(i).pulseTime());
@@ -1440,13 +1441,12 @@ public:
void test_sortByTimeAtSample_random_tof_time() {
for (int this_type = 0; this_type < 3; this_type++) {
EventType curType = static_cast<EventType>(this_type);
+ if (curType == WEIGHTED_NOTIME)
+ continue;
+
EventList el = this->fake_random_tof_constant_pulse_data(curType, 10);
el.switchTo(curType);
- if (curType == WEIGHTED_NOTIME) {
- continue;
- }
-
const double tofFactor = 1; // L1 / (L1 + L2)
const double tofShift = 0;
@@ -2040,6 +2040,45 @@ public:
} // starting event type
}
+ void test_compressFatEvents() {
+ // no pulse time should throw an exception
+ EventList el_notime_output;
+ EventList el_notime = this->fake_data(WEIGHTED_NOTIME);
+ TS_ASSERT_THROWS(el_notime.compressFatEvents(10., DateAndTime(0), 10.,
+ &el_notime_output),
+ std::invalid_argument);
+
+ // integration range
+ const double XMIN = 0.;
+ const double XMAX = 1.e7;
+
+ // regular events should compress decently well
+ EventList el_output;
+ this->fake_uniform_data_weights(TOF);
+ TS_ASSERT_THROWS_NOTHING(
+ el.compressFatEvents(20000., el.getPulseTimeMin(), 5., &el_output));
+ TS_ASSERT_EQUALS(el_output.getNumberEvents(), 400);
+ TS_ASSERT_EQUALS(el_output.integrate(XMIN, XMAX, true),
+ el.integrate(XMIN, XMAX, true));
+
+ // weighted events with time is the main use case
+ EventList el_weight_output;
+ this->fake_uniform_data_weights(WEIGHTED);
+ TS_ASSERT_THROWS_NOTHING(el.compressFatEvents(20000., el.getPulseTimeMin(),
+ 5., &el_weight_output));
+ TS_ASSERT_EQUALS(el_weight_output.getNumberEvents(), 400);
+ TS_ASSERT_EQUALS(el_weight_output.integrate(XMIN, XMAX, true),
+ el.integrate(XMIN, XMAX, true));
+
+ // change the start time to see that events don't make it in
+ el_weight_output.clear();
+ // previous run's this->fake_uniform_data_weights(WEIGHTED);
+ TS_ASSERT_THROWS_NOTHING(el.compressFatEvents(20000., el.getPulseTimeMax(),
+ 5., &el_weight_output));
+ TS_ASSERT_EQUALS(el_weight_output.getNumberEvents(), 1);
+ TS_ASSERT_DELTA(el_weight_output.integrate(XMIN, XMAX, true), 2., .0001);
+ }
+
void test_getEventsFrom() {
std::vector<TofEvent> *rel;
TS_ASSERT_THROWS_NOTHING(getEventsFrom(el, rel));
@@ -2283,15 +2322,37 @@ public:
// Mocking functions
//==================================================================================
- EventList fake_data() {
+ EventList fake_data(EventType eventType = TOF) {
// Clear the list
el = EventList();
+ if (eventType != TOF)
+ el.switchTo(eventType);
+
// Create some mostly-reasonable fake data.
srand(1234); // Fixed random seed
- for (int i = 0; i < NUMEVENTS; i++) {
- // Random tof up to 10 ms
- // Random pulse time up to 1000
- el += TofEvent(1e7 * (rand() * 1.0 / RAND_MAX), rand() % 1000);
+ switch (eventType) {
+ case TOF:
+ for (int i = 0; i < NUMEVENTS; i++) {
+ // Random tof up to 10 ms
+ // Random pulse time up to 1000
+ el += TofEvent(1e7 * (rand() * 1.0 / RAND_MAX), rand() % 1000);
+ }
+ break;
+ case WEIGHTED:
+ for (int i = 0; i < NUMEVENTS; i++) {
+ // Random tof up to 10 ms
+ // Random pulse time up to 1000
+ el += WeightedEvent(
+ TofEvent(1e7 * (rand() * 1.0 / RAND_MAX), rand() % 1000));
+ }
+ break;
+ case WEIGHTED_NOTIME:
+ for (int i = 0; i < NUMEVENTS; i++) {
+ // Random tof up to 10 ms
+ // Random pulse time up to 1000
+ el += TofEvent(1e7 * (rand() * 1.0 / RAND_MAX), rand() % 1000);
+ }
+ break;
}
return el;
}
@@ -2313,6 +2374,8 @@ public:
pulse_time += BIN_DELTA / events_per_bin) {
el += WeightedEvent(TofEvent(100, static_cast<size_t>(pulse_time)));
}
+ } else {
+ throw std::runtime_error("not implemented: fake_uniform_pulse_data");
}
return el;
}
@@ -2330,6 +2393,9 @@ public:
// Random tof up to 10 ms
el += WeightedEvent(TofEvent(1e7 * (rand() * 1.0 / RAND_MAX), 0));
}
+ } else {
+ throw std::runtime_error(
+ "not implemented: fake_random_tof_constant_pulse_data");
}
return el;
}
@@ -2358,15 +2424,27 @@ public:
}
/** Create a uniform event list with each event weight of 2.0, error 2.5 */
- void fake_uniform_data_weights() {
+ void fake_uniform_data_weights(EventType eventType = WEIGHTED) {
// Clear the list
el = EventList();
- el.switchTo(WEIGHTED);
+ if (eventType != TOF)
+ el.switchTo(WEIGHTED);
+
// Create some mostly-reasonable fake data.
srand(1234); // Fixed random seed
+ int64_t pulsetime = 0;
+ // 10^6 nanoseconds for pulse time delta
+ const int64_t PULSETIME_DELTA = static_cast<int64_t>(BIN_DELTA / 1000);
for (double tof = 100; tof < MAX_TOF; tof += BIN_DELTA / 2) {
+ // add some jitter into the pulse time with rand
+ pulsetime = 10000 * (static_cast<int64_t>(tof) / PULSETIME_DELTA) +
+ (rand() % 1000);
// tof steps of 5 microseconds, starting at 100 ns, up to 20 msec
- el += WeightedEvent(tof, rand() % 1000, 2.0, 2.5 * 2.5);
+ if (eventType == TOF) {
+ el += TofEvent(tof, pulsetime);
+ } else if (eventType == WEIGHTED) {
+ el += WeightedEvent(tof, pulsetime, 2.0, 2.5 * 2.5);
+ }
}
}