Mantid
Loading...
Searching...
No Matches
EventList.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
4// NScD Oak Ridge National Laboratory, European Spallation Source,
5// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6// SPDX - License - Identifier: GPL - 3.0 +
13#include "MantidKernel/Logger.h"
15#include "MantidKernel/Unit.h"
17
18#ifdef _MSC_VER
19// qualifier applied to function type has no meaning; ignored
20#pragma warning(disable : 4180)
21#endif
22#include "tbb/parallel_sort.h"
23#ifdef _MSC_VER
24#pragma warning(default : 4180)
25#endif
26
27#include <algorithm>
28#include <cfloat>
29#include <cmath>
30#include <functional>
31#include <limits>
32#include <stdexcept>
33
34using std::ostream;
35using std::runtime_error;
36using std::size_t;
37using std::vector;
38
39namespace Mantid::DataObjects {
40using Types::Core::DateAndTime;
41using Types::Event::TofEvent;
42using namespace Mantid::API;
43
44namespace {
45
46constexpr double SEC_TO_NANO{1.e9};
47
48// minimum event vector length to use tbb::parallel_sort
49// this is 4x what parallel_sort uses in the indidividual blocks
50constexpr size_t MIN_VEC_LENGTH_PARALLEL_SORT{2000};
51
59template <typename EventType>
60int64_t calculateCorrectedFullTime(const EventType &event, const double tofFactor, const double tofShift) {
61 return event.pulseTime().totalNanoseconds() +
62 static_cast<int64_t>(tofFactor * (event.tof() * 1.0E3) + (tofShift * 1.0E9));
63}
64
68template <typename EventType> class CompareTimeAtSample {
69private:
70 const double m_tofFactor;
71 const double m_tofShift;
72
73public:
74 CompareTimeAtSample(const double tofFactor, const double tofShift) : m_tofFactor(tofFactor), m_tofShift(tofShift) {}
75
85 bool operator()(const EventType &e1, const EventType &e2) const {
86 const auto tAtSample1 = calculateCorrectedFullTime(e1, m_tofFactor, m_tofShift);
87 const auto tAtSample2 = calculateCorrectedFullTime(e2, m_tofFactor, m_tofShift);
88 return (tAtSample1 < tAtSample2);
89 }
90};
91} // namespace
92//==========================================================================
95//==========================================================================
100bool compareEventPulseTime(const TofEvent &e1, const TofEvent &e2) { return (e1.pulseTime() < e2.pulseTime()); }
101
108bool compareEventPulseTimeTOF(const TofEvent &e1, const TofEvent &e2) {
109
110 if (e1.pulseTime() < e2.pulseTime()) {
111 return true;
112 } else if ((e1.pulseTime() == e2.pulseTime()) && (e1.tof() < e2.tof())) {
113 return true;
114 }
115
116 return false;
117}
118
119// comparator for pulse time with tolerance
121 explicit comparePulseTimeTOFDelta(const Types::Core::DateAndTime &start, const double seconds)
122 : startNano(start.totalNanoseconds()), deltaNano(static_cast<int64_t>(seconds * SEC_TO_NANO)) {}
123
124 bool operator()(const TofEvent &e1, const TofEvent &e2) {
125 // get the pulse times converted into bin number from start time
126 const int64_t e1Pulse = (e1.pulseTime().totalNanoseconds() - startNano) / deltaNano;
127 const int64_t e2Pulse = (e2.pulseTime().totalNanoseconds() - startNano) / deltaNano;
128
129 // compare with the calculated bin information
130 if (e1Pulse < e2Pulse) {
131 return true;
132 } else if ((e1Pulse == e2Pulse) && (e1.tof() < e2.tof())) {
133 return true;
134 }
135
136 return false;
137 }
138
139 int64_t startNano;
140 int64_t deltaNano;
141};
142
143struct FindBin {
144 double divisor;
145 double offset;
146 std::optional<size_t> (*findBin)(std::span<double const>, const double, const double, const double, const bool);
147 FindBin(double step, double xmin) {
148 if (step < 0) {
150 divisor = 1. / log1p(abs(step)); // use this to do change of base
151 offset = log(xmin) * divisor;
152 } else {
154 divisor = 1. / step;
155 offset = xmin * divisor;
156 }
157 }
158
159 std::optional<size_t> operator()(std::span<double const> X, const double tof, const bool findExact) {
160 return findBin(X, tof, divisor, offset, findExact);
161 }
162};
163
165// EventWorkspace is always histogram data and so is thus EventList
167 : m_histogram(HistogramData::Histogram::XMode::BinEdges, HistogramData::Histogram::YMode::Counts),
168 eventType(event_type), order(UNSORTED), mru(nullptr) {
169 switch (eventType) {
170 case TOF:
171 this->events = std::make_unique<std::vector<Mantid::Types::Event::TofEvent>>();
172 this->weightedEvents = nullptr;
173 this->weightedEventsNoTime = nullptr;
174 break;
175
176 case WEIGHTED:
177 this->events = nullptr;
178 this->weightedEvents = std::make_unique<std::vector<WeightedEvent>>();
179 this->weightedEventsNoTime = nullptr;
180 break;
181
182 case WEIGHTED_NOTIME:
183 this->events = nullptr;
184 this->weightedEventsNoTime = std::make_unique<std::vector<WeightedEventNoTime>>();
185 this->weightedEventsNoTime = nullptr;
186 break;
187 }
188}
189
195 : IEventList(specNo),
196 m_histogram(HistogramData::Histogram::XMode::BinEdges, HistogramData::Histogram::YMode::Counts),
197 weightedEvents(nullptr), weightedEventsNoTime(nullptr), eventType(TOF), order(UNSORTED), mru(mru) {
198 this->events = std::make_unique<std::vector<Mantid::Types::Event::TofEvent>>();
199}
200
203EventList::EventList(const EventList &rhs) : IEventList(rhs), m_histogram(rhs.m_histogram), mru{nullptr} {
204 // Note that operator= also assigns m_histogram, but the above use of the copy
205 // constructor avoid a memory allocation and is thus faster.
206 this->operator=(rhs);
207}
208
211EventList::EventList(const std::vector<Types::Event::TofEvent> &events)
212 : m_histogram(HistogramData::Histogram::XMode::BinEdges, HistogramData::Histogram::YMode::Counts),
213 weightedEvents(nullptr), weightedEventsNoTime(nullptr), eventType(TOF), mru(nullptr) {
214 this->events = std::make_unique<std::vector<Mantid::Types::Event::TofEvent>>(events.cbegin(), events.cend());
215 this->eventType = TOF;
216 this->order = UNSORTED;
217}
218
221EventList::EventList(const std::vector<WeightedEvent> &events)
222 : m_histogram(HistogramData::Histogram::XMode::BinEdges, HistogramData::Histogram::YMode::Counts), events(nullptr),
223 weightedEventsNoTime(nullptr), mru(nullptr) {
224 this->weightedEvents = std::make_unique<std::vector<WeightedEvent>>(events.cbegin(), events.cend());
225 this->eventType = WEIGHTED;
226 this->order = UNSORTED;
227}
228
231EventList::EventList(const std::vector<WeightedEventNoTime> &events)
232 : m_histogram(HistogramData::Histogram::XMode::BinEdges, HistogramData::Histogram::YMode::Counts), events(nullptr),
233 weightedEvents(nullptr), mru(nullptr) {
234 this->weightedEventsNoTime = std::make_unique<std::vector<WeightedEventNoTime>>(events.cbegin(), events.cend());
236 this->order = UNSORTED;
237}
238
241 // clear this out of the MRU (copy of code from EventList::clear()
242 if (mru) {
243 try {
244 mru->deleteIndex(this);
245 } catch (const std::runtime_error &) {
246 // this is an ignorable error
247 }
248 }
249
250 // set all member vectors to nullptr
251 this->events.reset();
252 this->weightedEvents.reset();
253 this->weightedEventsNoTime.reset();
254}
255
257void EventList::copyDataFrom(const ISpectrum &source) { source.copyDataInto(*this); }
258
262 if (events)
263 sink.events = std::make_unique<std::vector<Types::Event::TofEvent>>(events->cbegin(), events->cend());
264 else if (sink.events)
265 sink.events = std::make_unique<std::vector<Types::Event::TofEvent>>();
266 if (weightedEvents)
267 sink.weightedEvents =
268 std::make_unique<std::vector<WeightedEvent>>(weightedEvents->cbegin(), weightedEvents->cend());
269 else if (sink.weightedEvents)
270 sink.weightedEvents = std::make_unique<std::vector<WeightedEvent>>();
272 sink.weightedEventsNoTime = std::make_unique<std::vector<WeightedEventNoTime>>(weightedEventsNoTime->cbegin(),
273 weightedEventsNoTime->cend());
274 else if (sink.weightedEventsNoTime)
275 sink.weightedEventsNoTime = std::make_unique<std::vector<WeightedEventNoTime>>();
276
277 sink.eventType = eventType;
278 sink.order = order;
279}
280
283
284// --------------------------------------------------------------------------
296void EventList::createFromHistogram(const ISpectrum *inSpec, bool GenerateZeros, bool GenerateMultipleEvents,
297 int MaxEventsPerBin) {
298 // Fresh start
299 this->clear(true);
300
301 // Get the input histogram
302 Mantid::HistogramData::HistogramX const &X = inSpec->x();
303 Mantid::HistogramData::HistogramY const &Y = inSpec->y();
304 Mantid::HistogramData::HistogramE const &E = inSpec->e();
305 if (Y.size() + 1 != X.size()) {
306 throw std::runtime_error("Expected a histogram (X vector should be 1 longer than the Y vector)");
307 }
308
309 // Copy detector IDs and spectra
310 this->copyInfoFrom(*inSpec);
311 // We need weights but have no way to set the time. So use weighted, no time
313 if (GenerateZeros)
314 this->weightedEventsNoTime->reserve(Y.size());
315
316 for (size_t i = 0; i < X.size() - 1; i++) {
317 double weight = Y[i];
318 if ((weight != 0.0 || GenerateZeros) && std::isfinite(weight)) {
319 double error = E[i];
320 // Also check that the error is not a bad number
321 if (std::isfinite(error)) {
322 if (GenerateMultipleEvents) {
323 // --------- Multiple events per bin ----------
324 double errorSquared = error * error;
325 // Find how many events to fake
326 double val = weight / E[i];
327 val *= val;
328 // Convert to int with slight rounding up. This is to avoid rounding
329 // errors
330 auto numEvents = int(val + 0.2);
331 if (numEvents < 1)
332 numEvents = 1;
333 if (numEvents > MaxEventsPerBin)
334 numEvents = MaxEventsPerBin;
335 // Scale the weight and error for each
336 weight /= numEvents;
337 errorSquared /= numEvents;
338
339 // Spread the TOF. e.g. 2 events = 0.25, 0.75.
340 double tofStep = (X[i + 1] - X[i]) / (numEvents);
341 for (size_t j = 0; j < size_t(numEvents); j++) {
342 double tof = X[i] + tofStep * (0.5 + double(j));
343 // Create and add the event
344 // TODO: try emplace_back() here.
345 weightedEventsNoTime->emplace_back(tof, weight, errorSquared);
346 }
347 } else {
348 // --------- Single event per bin ----------
349 // TOF = midpoint of the bin
350 double tof = (X[i] + X[i + 1]) / 2.0;
351 // Error squared is carried in the event
352 double errorSquared = E[i];
353 errorSquared *= errorSquared;
354 // Create and add the event
355 weightedEventsNoTime->emplace_back(tof, weight, errorSquared);
356 }
357 } // error is nont NAN or infinite
358 } // weight is non-zero, not NAN, and non-infinite
359 } // (each bin)
360
361 // Set the X binning parameters
362 this->setSharedX(inSpec->sharedX());
363
364 // Manually set that this is sorted by TOF, since it is. This will make it
365 // "threadSafe" in other algos.
366 this->setSortOrder(TOF_SORT);
367}
368
369// --------------------------------------------------------------------------
370// --- Operators
371// -------------------------------------------------------------------
372
378 // Note that we are NOT copying the MRU pointer
379 // the EventWorkspace that possesses the EventList has already configured the mru
380 IEventList::operator=(rhs);
381 m_histogram = rhs.m_histogram;
382 rhs.copyDataInto(*this);
383 return *this;
384}
385
386// --------------------------------------------------------------------------
391EventList &EventList::operator+=(const Types::Event::TofEvent &event) {
392
393 switch (this->eventType) {
394 case TOF:
395 // Simply push the events
396 this->events->emplace_back(event);
397 break;
398
399 case WEIGHTED:
400 this->weightedEvents->emplace_back(event);
401 break;
402
403 case WEIGHTED_NOTIME:
404 this->weightedEventsNoTime->emplace_back(event);
405 break;
406 }
407
408 this->order = UNSORTED;
409 return *this;
410}
411
412// --------------------------------------------------------------------------
419EventList &EventList::operator+=(const std::vector<Types::Event::TofEvent> &more_events) {
420 switch (this->eventType) {
421 case TOF:
422 // Simply push the events
423 this->events->insert(this->events->end(), more_events.cbegin(), more_events.cend());
424 break;
425
426 case WEIGHTED:
427 // Add default weights to all the un-weighted incoming events from the list.
428 // and append to the list
429 this->weightedEvents->reserve(this->weightedEvents->size() + more_events.size());
430 std::copy(more_events.cbegin(), more_events.cend(), std::back_inserter(*this->weightedEvents));
431 break;
432
433 case WEIGHTED_NOTIME:
434 // Add default weights to all the un-weighted incoming events from the list.
435 // and append to the list
436 this->weightedEventsNoTime->reserve(this->weightedEventsNoTime->size() + more_events.size());
437 std::copy(more_events.cbegin(), more_events.cend(), std::back_inserter(*this->weightedEventsNoTime));
438 break;
439 }
440
441 this->order = UNSORTED;
442 return *this;
443}
444
445// --------------------------------------------------------------------------
454 this->switchTo(WEIGHTED);
455 this->weightedEvents->emplace_back(event);
456 this->order = UNSORTED;
457 return *this;
458}
459
460// --------------------------------------------------------------------------
468EventList &EventList::operator+=(const std::vector<WeightedEvent> &more_events) {
469 switch (this->eventType) {
470 case TOF:
471 // Need to switch to weighted
472 this->switchTo(WEIGHTED);
473 // Fall through to the insertion!
474
475 case WEIGHTED:
476 // Append the two lists
477 this->weightedEvents->insert(weightedEvents->end(), more_events.cbegin(), more_events.cend());
478 break;
479
480 case WEIGHTED_NOTIME:
481 // Add default weights to all the un-weighted incoming events from the list.
482 // and append to the list
483 this->weightedEventsNoTime->reserve(this->weightedEventsNoTime->size() + more_events.size());
484 std::copy(more_events.cbegin(), more_events.cend(), std::back_inserter(*this->weightedEventsNoTime));
485 break;
486 }
487
488 this->order = UNSORTED;
489 return *this;
490}
491
492// --------------------------------------------------------------------------
500EventList &EventList::operator+=(const std::vector<WeightedEventNoTime> &more_events) {
501 switch (this->eventType) {
502 case TOF:
503 case WEIGHTED:
504 // Need to switch to weighted with no time
506 // Fall through to the insertion!
507
508 case WEIGHTED_NOTIME:
509 // Simple appending of the two lists
510 this->weightedEventsNoTime->insert(weightedEventsNoTime->end(), more_events.cbegin(), more_events.cend());
511 break;
512 }
513
514 this->order = UNSORTED;
515 return *this;
516}
517
518// --------------------------------------------------------------------------
528 if (!more_events.empty()) {
529 // We'll let the += operator for the given vector of event lists handle it
530 switch (more_events.getEventType()) {
531 case TOF:
532 this->operator+=(*more_events.events);
533 break;
534
535 case WEIGHTED:
536 this->operator+=(*more_events.weightedEvents);
537 break;
538
539 case WEIGHTED_NOTIME:
540 this->operator+=(*more_events.weightedEventsNoTime);
541 break;
542 }
543
544 // No guaranteed order
545 if (this->empty()) {
546 this->order = more_events.order;
547 } else {
548 this->order = UNSORTED;
549 }
550 }
551
552 // Do a union between the detector IDs of both lists
553 addDetectorIDs(more_events.getDetectorIDs());
554
555 return *this;
556}
557
558// --------------------------------------------------------------------------
567template <class T1, class T2> void EventList::minusHelper(std::vector<T1> &events, const std::vector<T2> &more_events) {
568 // Make the end vector big enough in one go (avoids repeated re-allocations).
569 events.reserve(events.size() + more_events.size());
570 /* In the event of subtracting in place, calling the end() vector would make
571 * it point at the wrong place
572 * Using it caused a segault, Ticket #2306.
573 * So we cache the end (this speeds up too).
574 */
575 // We call the constructor for T1. In the case of WeightedEventNoTime, the pulse time will just be ignored.
576 std::transform(more_events.cbegin(), more_events.cend(), std::back_inserter(events),
577 [](const auto &ev) { return T1(ev.tof(), ev.pulseTime(), ev.weight() * (-1.0), ev.errorSquared()); });
578}
579
580// --------------------------------------------------------------------------
589 if (this == &more_events) {
590 // Special case, ticket #3844 part 2.
591 // When doing this = this - this,
592 // simply clear the input event list. Saves memory!
593 this->clearData();
594 return *this;
595 }
596
597 // We'll let the -= operator for the given vector of event lists handle it
598 switch (this->getEventType()) {
599 case TOF:
600 this->switchTo(WEIGHTED);
601 // Fall through
602
603 case WEIGHTED:
604 switch (more_events.getEventType()) {
605 case TOF:
606 minusHelper(*this->weightedEvents, *more_events.events);
607 break;
608 case WEIGHTED:
609 minusHelper(*this->weightedEvents, *more_events.weightedEvents);
610 break;
611 case WEIGHTED_NOTIME:
612 // TODO: Should this throw?
613 minusHelper(*this->weightedEvents, *more_events.weightedEventsNoTime);
614 break;
615 }
616 break;
617
618 case WEIGHTED_NOTIME:
619 switch (more_events.getEventType()) {
620 case TOF:
621 minusHelper(*this->weightedEventsNoTime, *more_events.events);
622 break;
623 case WEIGHTED:
624 minusHelper(*this->weightedEventsNoTime, *more_events.weightedEvents);
625 break;
626 case WEIGHTED_NOTIME:
628 break;
629 }
630 break;
631 }
632
633 // No guaranteed order
634 this->order = UNSORTED;
635
636 // NOTE: What to do about detector ID's?
637 return *this;
638}
639
640namespace {
641/*
642 * Both can be nullptr, or the values can be equal, but do not have one nullptr
643 */
644template <typename T>
645bool vectorPtrEquals(const std::unique_ptr<std::vector<T>> &left, const std::unique_ptr<std::vector<T>> &right) {
646 if (left && right) {
647 return (*left == *right);
648 ;
649 } else if ((left && !right) || (right && !left)) {
650 return false;
651 }
652 return true;
653}
654} // anonymous namespace
655
656// --------------------------------------------------------------------------
662 if (this->getNumberEvents() != rhs.getNumberEvents())
663 return false;
664 if (this->eventType != rhs.eventType)
665 return false;
666 if (this->empty())
667 return true;
668 // Check all event lists; The empty ones will compare equal
669 if (!vectorPtrEquals(events, rhs.events))
670 return false;
671 if (!vectorPtrEquals(weightedEvents, rhs.weightedEvents))
672 return false;
673 if (!vectorPtrEquals(weightedEventsNoTime, rhs.weightedEventsNoTime))
674 return false;
675
676 // nothing wasn't equal, so they are equal
677 return true;
678}
679
684bool EventList::operator!=(const EventList &rhs) const { return (!this->operator==(rhs)); }
685
686bool EventList::equals(const EventList &rhs, const double tolTof, const double tolWeight,
687 const int64_t tolPulse) const {
688 // generic checks
689 if (this->getNumberEvents() != rhs.getNumberEvents())
690 return false;
691 if (this->eventType != rhs.eventType)
692 return false;
693 if (this->empty())
694 return true;
695
696 // loop over the events
697 switch (this->eventType) {
698 case TOF: {
699 auto leftIter = this->events->cbegin();
700 auto leftEnd = this->events->cend();
701 auto rightIter = rhs.events->cbegin();
702 while (leftIter != leftEnd) {
703 if (!leftIter->equals(*rightIter, tolTof, tolPulse))
704 return false;
705 leftIter = std::next(leftIter);
706 rightIter = std::next(rightIter);
707 }
708 break;
709 }
710 case WEIGHTED: {
711 auto leftIter = this->weightedEvents->cbegin();
712 auto leftEnd = this->weightedEvents->cend();
713 auto rightIter = rhs.weightedEvents->cbegin();
714 while (leftIter != leftEnd) {
715 if (!leftIter->equals(*rightIter, tolTof, tolWeight, tolPulse))
716 return false;
717 leftIter = std::next(leftIter);
718 rightIter = std::next(rightIter);
719 }
720 break;
721 }
722 case WEIGHTED_NOTIME: {
723 auto leftIter = this->weightedEventsNoTime->cbegin();
724 auto leftEnd = this->weightedEventsNoTime->cend();
725 auto rightIter = rhs.weightedEventsNoTime->cbegin();
726 while (leftIter != leftEnd) {
727 if (!leftIter->equals(*rightIter, tolTof, tolWeight))
728 return false;
729 leftIter = std::next(leftIter);
730 rightIter = std::next(rightIter);
731 }
732 break;
733 }
734 default:
735 break;
736 }
737
738 // anything that gets this far is equal within tolerances
739 return true;
740}
741
742// -----------------------------------------------------------------------------------------------
747
748// -----------------------------------------------------------------------------------------------
753 switch (newType) {
754 case TOF:
755 if (eventType != TOF)
756 throw std::runtime_error("EventList::switchTo() called on an EventList with weights to go down to TofEvent's. "
757 "This would remove weight information and therefore is not possible.");
758 break;
759
760 case WEIGHTED:
762 break;
763
764 case WEIGHTED_NOTIME:
766 break;
767 }
768 // Make sure to free memory
769 this->clearUnused();
770}
771
772// -----------------------------------------------------------------------------------------------
777 switch (eventType) {
778 case WEIGHTED:
779 // Do nothing; it already is weighted
780 return;
781
782 case WEIGHTED_NOTIME:
783 throw std::runtime_error("EventList::switchToWeightedEvents() called on an EventList with WeightedEventNoTime's. "
784 "It has lost the pulse time information and can't go back to WeightedEvent's.");
785 break;
786
787 case TOF:
788 if (events && !events->empty()) {
789 // Convert and copy all TofEvents to the weightedEvents list.
790 weightedEvents = std::make_unique<std::vector<WeightedEvent>>(events->cbegin(), events->cend());
791 // Get rid of the old events
792 events.reset();
793 } else {
794 weightedEvents = std::make_unique<std::vector<WeightedEvent>>();
795 }
797 break;
798 }
799}
800
801// -----------------------------------------------------------------------------------------------
806 switch (eventType) {
807 case WEIGHTED_NOTIME:
808 // Do nothing if already there
809 return;
810
811 case TOF: {
812 if (events && !events->empty()) {
813 // Convert and copy all TofEvents to the weightedEvents list.
814 this->weightedEventsNoTime = std::make_unique<std::vector<WeightedEventNoTime>>(events->cbegin(), events->cend());
815 // Get rid of the old events
816 events.reset();
817 } else {
818 this->weightedEventsNoTime = std::make_unique<std::vector<WeightedEventNoTime>>();
819 }
820 break;
821 }
822
823 case WEIGHTED: {
824 // Convert and copy all TofEvents to the weightedEvents list.
825 if (weightedEvents && !weightedEvents->empty()) {
827 std::make_unique<std::vector<WeightedEventNoTime>>(weightedEvents->cbegin(), weightedEvents->cend());
828 // Get rid of the old events
829 weightedEvents.reset();
830 } else {
831 this->weightedEventsNoTime = std::make_unique<std::vector<WeightedEventNoTime>>();
832 }
833 break;
834 }
835 }
837}
838
839// ==============================================================================================
840// --- Testing functions (mostly)
841// ---------------------------------------------------------------
842// ==============================================================================================
843
850WeightedEvent EventList::getEvent(size_t event_number) {
851 switch (eventType) {
852 case TOF:
853 return WeightedEvent(events->at(event_number));
854 case WEIGHTED:
855 return weightedEvents->at(event_number);
856 case WEIGHTED_NOTIME: {
857 const auto event = weightedEventsNoTime->at(event_number);
858 return WeightedEvent(event.tof(), 0, event.weight(), event.errorSquared());
859 }
860 }
861 throw std::runtime_error("EventList: invalid event type value was found.");
862}
863
864// ==============================================================================================
865// --- Handling the event list
866// -------------------------------------------------------------------
867// ==============================================================================================
868
876const std::vector<TofEvent> &EventList::getEvents() const {
877 if (eventType != TOF)
878 throw std::runtime_error("EventList::getEvents() called for an EventList that has weights. Use getWeightedEvents() "
879 "or getWeightedEventsNoTime().");
880 if (this->events)
881 return *this->events;
882 else
883 throw std::runtime_error("unweighted event vector is not initialized");
884}
885
892std::vector<TofEvent> &EventList::getEvents() {
893 if (eventType != TOF)
894 throw std::runtime_error("EventList::getEvents() called for an EventList that has weights. Use getWeightedEvents() "
895 "or getWeightedEventsNoTime().");
896 if (this->events)
897 return *this->events;
898 else
899 throw std::runtime_error("unweighted event vector is not initialized");
900}
901
909std::vector<WeightedEvent> &EventList::getWeightedEvents() {
910 if (eventType != WEIGHTED)
911 throw std::runtime_error("EventList::getWeightedEvents() called for an EventList not of type WeightedEvent. Use "
912 "getEvents() or getWeightedEventsNoTime().");
913 if (this->weightedEvents)
914 return *this->weightedEvents;
915 else
916 throw std::runtime_error("weighted event vector is not initialized");
917}
918
926const std::vector<WeightedEvent> &EventList::getWeightedEvents() const {
927 if (eventType != WEIGHTED)
928 throw std::runtime_error("EventList::getWeightedEvents() called for an EventList not of type WeightedEvent. Use "
929 "getEvents() or getWeightedEventsNoTime().");
930 if (this->weightedEvents)
931 return *this->weightedEvents;
932 else
933 throw std::runtime_error("weighted event vector is not initialed");
934}
935
941std::vector<WeightedEventNoTime> &EventList::getWeightedEventsNoTime() {
943 throw std::runtime_error("EventList::getWeightedEventsNoTime() called for an EventList not of type "
944 "WeightedEventNoTime. Use getEvents() or getWeightedEvents().");
945 if (this->weightedEventsNoTime)
946 return *this->weightedEventsNoTime;
947 else
948 throw std::runtime_error("weighted event no time vector is not initialed");
949}
950
956const std::vector<WeightedEventNoTime> &EventList::getWeightedEventsNoTime() const {
958 throw std::runtime_error("EventList::getWeightedEventsNoTime() called for an EventList not of type "
959 "WeightedEventNoTime. Use getEvents() or getWeightedEvents().");
960 if (this->weightedEventsNoTime)
961 return *this->weightedEventsNoTime;
962 else
963 throw std::runtime_error("weighted event no time vector is not initialed");
964}
965
969void EventList::clear(const bool removeDetIDs) {
970 if (mru) {
971 try {
972 mru->deleteIndex(this);
973 } catch (const std::runtime_error &) {
974 // this is an ignorable error
975 }
976 }
977 // clear representations that aren't for the current type
978 this->clearUnused();
979
980 // release unused memory or allocate new vector
981 // rather than creating a new object, reset existing pointer
982 if (!this->empty()) {
983 if (this->events && eventType == TOF) {
984 this->events->clear();
985 std::vector<TofEvent>().swap(*this->events); // STL Trick to release memory
986 }
987 if (this->weightedEvents && eventType == WEIGHTED) {
988 this->weightedEvents->clear();
989 std::vector<WeightedEvent>().swap(*this->weightedEvents); // STL Trick to release memory
990 }
992 this->weightedEventsNoTime->clear();
993 std::vector<WeightedEventNoTime>().swap(*this->weightedEventsNoTime); // STL Trick to release memory
994 }
995 }
996 if (removeDetIDs)
997 this->clearDetectorIDs();
998}
999
1005 if (eventType != TOF && (this->events)) {
1006 this->events.reset();
1007 }
1008 if (eventType != WEIGHTED && (this->weightedEvents)) {
1009 this->weightedEvents.reset();
1010 }
1012 this->weightedEventsNoTime.reset();
1013 }
1014}
1015
1017void EventList::clearData() { this->clear(false); }
1018
1023void EventList::setMRU(EventWorkspaceMRU *newMRU) { mru = newMRU; }
1024
1032void EventList::reserve(size_t num) {
1033 switch (this->eventType) {
1034 case TOF:
1035 this->events->reserve(num);
1036 break;
1037 case WEIGHTED:
1038 this->weightedEvents->reserve(num);
1039 break;
1040 case WEIGHTED_NOTIME:
1041 this->weightedEventsNoTime->reserve(num);
1042 break;
1043 }
1044}
1045
1046// ==============================================================================================
1047// --- Sorting functions -----------------------------------------------------
1048// ==============================================================================================
1049
1050// --------------------------------------------------------------------------
1054void EventList::sort(const EventSortType order) const {
1055 if (order == UNSORTED) {
1056 return; // don't bother doing anything. Why did you ask to unsort?
1057 } else if (order == TOF_SORT) {
1058 this->sortTof();
1059 } else if (order == PULSETIME_SORT) {
1060 this->sortPulseTime();
1061 } else if (order == PULSETIMETOF_SORT) {
1062 this->sortPulseTimeTOF();
1063 } else if (order == PULSETIMETOF_DELTA_SORT) {
1064 throw std::invalid_argument("sorting by pulse time with delta requires "
1065 "extra parameters. Use sortPulseTimeTOFDelta "
1066 "instead.");
1067 } else if (order == TIMEATSAMPLE_SORT) {
1068 throw std::invalid_argument("sorting by time at sample requires extra "
1069 "parameters. Use sortTimeAtSample instead.");
1070 } else {
1071 throw runtime_error("Invalid sort type in EventList::sort(EventSortType)");
1072 }
1073}
1074
1075// --------------------------------------------------------------------------
1080void EventList::setSortOrder(const EventSortType order) const { this->order = order; }
1081
1082namespace {
1083// these are abstractions
1084template <class RandomIt> void switchable_sort(RandomIt first, RandomIt last) {
1085 const auto vec_size = static_cast<size_t>(std::distance(first, last));
1086 if (vec_size < 2)
1087 return;
1088 else if (vec_size < MIN_VEC_LENGTH_PARALLEL_SORT)
1089 std::sort(first, last);
1090 else
1091 tbb::parallel_sort(first, last);
1092}
1093
1094template <class RandomIt, class Compare> void switchable_sort(RandomIt first, RandomIt last, Compare comp) {
1095 const auto vec_size = static_cast<size_t>(std::distance(first, last));
1096 if (vec_size < 2)
1097 return;
1098 else if (vec_size < MIN_VEC_LENGTH_PARALLEL_SORT)
1099 std::sort(first, last, std::move(comp));
1100 else
1101 tbb::parallel_sort(first, last, comp);
1102}
1103} // anonymous namespace
1104
1105// --------------------------------------------------------------------------
1108 // nothing to do
1109 if (this->order == TOF_SORT)
1110 return;
1111
1112 // Avoid sorting from multiple threads
1113 std::lock_guard<std::mutex> _lock(m_sortMutex);
1114 // If the list was sorted while waiting for the lock, return.
1115 if (this->order == TOF_SORT) // cppcheck-suppress identicalConditionAfterEarlyExit
1116 return;
1117
1118 switch (eventType) {
1119 case TOF:
1120 switchable_sort(events->begin(), events->end());
1121 break;
1122 case WEIGHTED:
1123 switchable_sort(weightedEvents->begin(), weightedEvents->end());
1124 break;
1125 case WEIGHTED_NOTIME:
1126 switchable_sort(weightedEventsNoTime->begin(), weightedEventsNoTime->end());
1127 break;
1128 }
1129 // Save the order to avoid unnecessary re-sorting.
1130 this->order = TOF_SORT;
1131}
1132
1133// --------------------------------------------------------------------------
1142void EventList::sortTimeAtSample(const double &tofFactor, const double &tofShift, bool forceResort) const {
1143 // Check pre-cached sort flag.
1144 if (this->order == TIMEATSAMPLE_SORT && !forceResort)
1145 return;
1146
1147 // Avoid sorting from multiple threads
1148 std::lock_guard<std::mutex> _lock(m_sortMutex);
1149 // If the list was sorted while waiting for the lock, return.
1150 if (this->order == TIMEATSAMPLE_SORT && !forceResort)
1151 return;
1152
1153 // Perform sort.
1154 switch (eventType) {
1155 case TOF: {
1156 CompareTimeAtSample<TofEvent> comparitor(tofFactor, tofShift);
1157 switchable_sort(events->begin(), events->end(), comparitor);
1158 } break;
1159 case WEIGHTED: {
1160 CompareTimeAtSample<WeightedEvent> comparitor(tofFactor, tofShift);
1161 switchable_sort(weightedEvents->begin(), weightedEvents->end(), comparitor);
1162 } break;
1163 case WEIGHTED_NOTIME: {
1164 CompareTimeAtSample<WeightedEventNoTime> comparitor(tofFactor, tofShift);
1165 switchable_sort(weightedEventsNoTime->begin(), weightedEventsNoTime->end(), comparitor);
1166 } break;
1167 }
1168 // Save the order to avoid unnecessary re-sorting.
1169 this->order = TIMEATSAMPLE_SORT;
1170}
1171
1172// --------------------------------------------------------------------------
1175 if (this->order == PULSETIME_SORT || this->order == PULSETIMETOF_SORT)
1176 return; // nothing to do
1177
1178 // Avoid sorting from multiple threads
1179 std::lock_guard<std::mutex> _lock(m_sortMutex);
1180 // If the list was sorted while waiting for the lock, return.
1181 if (this->order == PULSETIME_SORT)
1182 return;
1183
1184 // Perform sort.
1185 switch (eventType) {
1186 case TOF:
1187 switchable_sort(events->begin(), events->end(), compareEventPulseTime);
1188 break;
1189 case WEIGHTED:
1190 switchable_sort(weightedEvents->begin(), weightedEvents->end(), compareEventPulseTime);
1191 break;
1192 case WEIGHTED_NOTIME:
1193 // Do nothing; there is no time to sort
1194 break;
1195 }
1196 // Save the order to avoid unnecessary re-sorting.
1197 this->order = PULSETIME_SORT;
1198}
1199
1200/*
1201 * Sort events by pulse time + TOF
1202 * (the absolute time)
1203 */
1205 if (this->order == PULSETIMETOF_SORT)
1206 return; // already ordered
1207
1208 // Avoid sorting from multiple threads
1209 std::lock_guard<std::mutex> _lock(m_sortMutex);
1210 // If the list was sorted while waiting for the lock, return.
1211 if (this->order == PULSETIMETOF_SORT) // cppcheck-suppress identicalConditionAfterEarlyExit
1212 return;
1213
1214 switch (eventType) {
1215 case TOF:
1216 switchable_sort(events->begin(), events->end(), compareEventPulseTimeTOF);
1217 break;
1218 case WEIGHTED:
1219 switchable_sort(weightedEvents->begin(), weightedEvents->end(), compareEventPulseTimeTOF);
1220 break;
1221 case WEIGHTED_NOTIME:
1222 // Do nothing; there is no time to sort
1223 break;
1224 }
1225
1226 // Save
1227 this->order = PULSETIMETOF_SORT;
1228}
1229
1237void EventList::sortPulseTimeTOFDelta(const Types::Core::DateAndTime &start, const double seconds) const {
1238 // Avoid sorting from multiple threads
1239 std::lock_guard<std::mutex> _lock(m_sortMutex);
1240
1241 std::function<bool(const TofEvent &, const TofEvent &)> comparator = comparePulseTimeTOFDelta(start, seconds);
1242
1243 switch (eventType) {
1244 case TOF:
1245 switchable_sort(events->begin(), events->end(), std::move(comparator));
1246 break;
1247 case WEIGHTED:
1248 switchable_sort(weightedEvents->begin(), weightedEvents->end(), std::move(comparator));
1249 break;
1250 case WEIGHTED_NOTIME:
1251 // Do nothing; there is no time to sort
1252 break;
1253 }
1254
1255 this->order = UNSORTED; // so the function always re-runs
1256}
1257
1258// --------------------------------------------------------------------------
1260bool EventList::isSortedByTof() const { return (this->order == TOF_SORT); }
1261
1262// --------------------------------------------------------------------------
1265
1266// --------------------------------------------------------------------------
1273 // reverse the histogram bin parameters. Changing X invalidates any histogram cached in the MRU,
1274 // so drop it first -- this is what the deprecated dataX() used to do implicitly.
1275 if (mru)
1276 mru->deleteIndex(this);
1277 auto &x = mutableX();
1278 std::reverse(x.begin(), x.end());
1279
1280 // flip the events if they are tof sorted
1281 if (this->isSortedByTof()) {
1282 switch (eventType) {
1283 case TOF:
1284 std::reverse(this->events->begin(), this->events->end());
1285 break;
1286 case WEIGHTED:
1287 std::reverse(this->weightedEvents->begin(), this->weightedEvents->end());
1288 break;
1289 case WEIGHTED_NOTIME:
1290 std::reverse(this->weightedEventsNoTime->begin(), this->weightedEventsNoTime->end());
1291 break;
1292 }
1293 // And we are still sorted! :)
1294 }
1295 // Otherwise, do nothing. If it was sorted by pulse time, then it still is
1296}
1297
1298// --------------------------------------------------------------------------
1307 switch (eventType) {
1308 case TOF:
1309 return (this->events) ? this->events->size() : 0;
1310 case WEIGHTED:
1311 return (this->weightedEvents) ? this->weightedEvents->size() : 0;
1312 case WEIGHTED_NOTIME:
1313 return (this->weightedEventsNoTime) ? this->weightedEventsNoTime->size() : 0;
1314 }
1315 throw std::runtime_error("EventList: invalid event type value was found.");
1316}
1317
1321bool EventList::empty() const {
1322 switch (eventType) {
1323 case TOF:
1324 if (this->events)
1325 return this->events->empty();
1326 else
1327 throw std::runtime_error("TOF events is nullptr");
1328 case WEIGHTED:
1329 if (this->weightedEvents)
1330 return this->weightedEvents->empty();
1331 else
1332 throw std::runtime_error("WEIGHTED events is nullptr");
1333 case WEIGHTED_NOTIME:
1334 if (this->weightedEventsNoTime)
1335 return this->weightedEventsNoTime->empty();
1336 else
1337 throw std::runtime_error("WEIGHTED_NOTIME events is nullptr");
1338 }
1339 throw std::runtime_error("EventList: invalid event type value was found.");
1340}
1341
1342// --------------------------------------------------------------------------
1350 switch (eventType) {
1351 case TOF:
1352 return this->events->capacity() * sizeof(TofEvent) + sizeof(EventList);
1353 case WEIGHTED:
1354 return this->weightedEvents->capacity() * sizeof(WeightedEvent) + sizeof(EventList);
1355 case WEIGHTED_NOTIME:
1356 return this->weightedEventsNoTime->capacity() * sizeof(WeightedEventNoTime) + sizeof(EventList);
1357 }
1358 throw std::runtime_error("EventList: invalid event type value was found.");
1359}
1360
1361// --------------------------------------------------------------------------
1365 size_t x_size = x().size();
1366 if (x_size > 1)
1367 return x_size - 1;
1368 else
1369 return 0;
1370}
1371
1372// ==============================================================================================
1373// --- Setting the Histogram X axis, without recalculating the histogram
1374// -----------------------
1375// ==============================================================================================
1376
1382 m_histogram.setSharedX(X);
1383 if (mru)
1384 mru->deleteIndex(this);
1385}
1386
1391 if (mru)
1392 mru->deleteIndex(this);
1393 // Handing out a mutable MantidVec would allow the length to be changed, so
1394 // FixedLengthVector::mutableRawData() is protected and Histogram is its only friend. The
1395 // deprecated legacy accessor is therefore the only way to express this.
1396 GNU_DIAG_OFF("deprecated-declarations")
1397 MSVC_DIAG_OFF(4996)
1398 return m_histogram.dataX();
1399 MSVC_DIAG_ON(4996)
1400 GNU_DIAG_ON("deprecated-declarations")
1401}
1402
1405const MantidVec &EventList::dataX() const { return m_histogram.x().rawData(); }
1406
1408const MantidVec &EventList::readX() const { return m_histogram.x().rawData(); }
1409
1412
1414GNU_DIAG_OFF("deprecated-declarations")
1415MSVC_DIAG_OFF(4996)
1416MantidVec &EventList::dataDx() { return m_histogram.dataDx(); }
1417MSVC_DIAG_ON(4996)
1418GNU_DIAG_ON("deprecated-declarations")
1420const MantidVec &EventList::dataDx() const { return m_histogram.dx().rawData(); }
1422const MantidVec &EventList::readDx() const { return m_histogram.dx().rawData(); }
1423
1424// ==============================================================================================
1425// --- Return Data Vectors --------------------------------------------------
1426// ==============================================================================================
1427
1434 auto Y = new MantidVec();
1435 MantidVec E;
1436 // Generate the Y histogram while skipping the E if possible.
1437 generateHistogram(x().rawData(), *Y, E, true);
1438 return Y;
1439}
1440
1447 MantidVec Y;
1448 auto E = new MantidVec();
1449 generateHistogram(x().rawData(), Y, *E);
1450 // Y is unused.
1451 return E;
1452}
1453
1455HistogramData::Histogram EventList::getHistogram() const { return m_histogram; }
1456
1457HistogramData::Histogram EventList::histogram() const {
1458 HistogramData::Histogram ret(m_histogram);
1459 ret.setSharedY(sharedY());
1460 ret.setSharedE(sharedE());
1461 return ret;
1462}
1463
1464HistogramData::Counts EventList::counts() const { return histogram().counts(); }
1465
1466HistogramData::CountVariances EventList::countVariances() const { return histogram().countVariances(); }
1467
1468HistogramData::CountStandardDeviations EventList::countStandardDeviations() const {
1469 return histogram().countStandardDeviations();
1470}
1471
1472HistogramData::Frequencies EventList::frequencies() const { return histogram().frequencies(); }
1473
1474HistogramData::FrequencyVariances EventList::frequencyVariances() const { return histogram().frequencyVariances(); }
1475
1476HistogramData::FrequencyStandardDeviations EventList::frequencyStandardDeviations() const {
1477 return histogram().frequencyStandardDeviations();
1478}
1479
1480const HistogramData::HistogramY &EventList::y() const {
1481 if (!mru)
1482 throw std::runtime_error("'EventList::y()' called with no MRU set. This is not allowed.");
1483
1484 return *sharedY();
1485}
1486const HistogramData::HistogramE &EventList::e() const {
1487 if (!mru)
1488 throw std::runtime_error("'EventList::e()' called with no MRU set. This is not allowed.");
1489
1490 return *sharedE();
1491}
1493 // This is the thread number from which this function was called.
1494 const int thread = PARALLEL_THREAD_NUMBER;
1495
1497
1498 // Is the data in the mrulist?
1499 if (mru) {
1500 mru->ensureEnoughBuffersY(static_cast<size_t>(thread));
1501 yData = mru->findY(static_cast<size_t>(thread), this);
1502 }
1503
1504 if (!yData) {
1505 MantidVec Y;
1506 MantidVec E;
1507 this->generateHistogram(x().rawData(), Y, E);
1508
1509 // Create the MRU object
1510 yData = Kernel::make_cow<HistogramData::HistogramY>(std::move(Y));
1511
1512 // Lets save it in the MRU
1513 if (mru) {
1514 mru->insertY(thread, yData, this);
1515 auto eData = Kernel::make_cow<HistogramData::HistogramE>(std::move(E));
1516 mru->ensureEnoughBuffersE(thread);
1517 mru->insertE(thread, eData, this);
1518 }
1519 }
1520 return yData;
1521}
1523 // This is the thread number from which this function was called.
1524 const auto thread = static_cast<size_t>(PARALLEL_THREAD_NUMBER);
1525
1527
1528 // Is the data in the mrulist?
1529 if (mru) {
1530 mru->ensureEnoughBuffersE(thread);
1531 eData = mru->findE(thread, this);
1532 }
1533
1534 if (!eData) {
1535 // Now use that to get E -- Y values are generated from another function
1536 MantidVec Y_ignored;
1537 MantidVec E;
1538 this->generateHistogram(x().rawData(), Y_ignored, E);
1539 eData = Kernel::make_cow<HistogramData::HistogramE>(std::move(E));
1540
1541 // Lets save it in the MRU
1542 if (mru)
1543 mru->insertE(thread, eData, this);
1544 }
1545 return eData;
1546}
1553 if (!mru)
1554 throw std::runtime_error("'EventList::dataY()' called with no MRU set. This is not allowed.");
1555
1556 // WARNING: The Y data of sharedY() is stored in MRU, returning reference fine
1557 // as long as it stays there.
1558 return sharedY()->rawData();
1559}
1560
1567 if (!mru)
1568 throw std::runtime_error("'EventList::dataE()' called with no MRU set. This is not allowed.");
1569
1570 // WARNING: The E data of sharedE() is stored in MRU, returning reference fine
1571 // as long as it stays there.
1572 return sharedE()->rawData();
1573}
1574
1575namespace {
1576inline double calcNorm(const double errorSquared) {
1577 if (errorSquared == 0.)
1578 return 0;
1579 else if (errorSquared == 1.)
1580 return 1.;
1581 else
1582 return 1. / std::sqrt(errorSquared);
1583}
1584} // namespace
1585
1586// --------------------------------------------------------------------------
1595template <class T>
1596inline void EventList::compressEventsHelper(const std::vector<T> &events, std::vector<WeightedEventNoTime> &out,
1597 double tolerance) {
1598 // Clear the output. We can't know ahead of time how much space to reserve :(
1599 out.clear();
1600 // We will make a starting guess of 1/20th of the number of input events.
1601 out.reserve(events.size() / 20);
1602
1603 // The last TOF to which we are comparing.
1604 double lastTof = events.front().m_tof;
1605 // For getting an accurate average TOF
1606 double totalTof = 0;
1607 int num = 0;
1608 // Carrying weight, error, and normalization
1609 double weight = 0;
1610 double errorSquared = 0;
1611 double normalization = 0.;
1612
1613 double bin_end = lastTof;
1614 std::function<bool(const double, const double)> compareTof;
1615 std::function<double(const double, double)> next_bin;
1616
1617 if (tolerance < 0) { // log
1618 if (lastTof < 0)
1619 throw std::runtime_error("compressEvents with log binning doesn't work with negative TOF");
1620
1621 if (lastTof == 0)
1622 bin_end = fabs(tolerance);
1623
1624 // for log we do "less than" so that is matches the log binning of the Rebin algorithm
1625 compareTof = [](const double lhs, const double rhs) { return lhs < rhs; };
1626 next_bin = [tolerance](const double lastTof, double bin_end) {
1627 // advance the bin_end until we find the one that this next event falls into
1628 while (lastTof >= bin_end)
1629 bin_end = bin_end * (1 - tolerance);
1630 return bin_end;
1631 };
1632 } else { // linear
1633 // for linear we do "less than or equals" because that is how it was originally implemented
1634 compareTof = [](const double lhs, const double rhs) { return lhs <= rhs; };
1635 next_bin = [tolerance](const double lastTof, double) { return lastTof + tolerance; };
1636 }
1637
1638 // get first bin_end
1639 bin_end = next_bin(lastTof, bin_end);
1640
1641 for (auto it = events.cbegin(); it != events.cend(); it++) {
1642 if (compareTof(it->m_tof, bin_end)) {
1643 // Carry the error and weight
1644 weight += it->weight();
1645 errorSquared += it->errorSquared();
1646 // Track the average tof
1647 num++;
1648 const double norm = calcNorm(it->errorSquared());
1649 normalization += norm;
1650 totalTof += it->m_tof * norm;
1651 } else {
1652 // We exceeded the tolerance
1653 // Create a new event with the average TOF and summed weights and
1654 // squared errors.
1655 if (num == 1) {
1656 // last time-of-flight is the only one contributing
1657 out.emplace_back(lastTof, weight, errorSquared);
1658 } else if (num > 1) {
1659 out.emplace_back(totalTof / normalization, weight, errorSquared);
1660 }
1661 // Start a new combined object
1662 num = 1;
1663 const double norm = calcNorm(it->errorSquared());
1664 normalization = norm;
1665 totalTof = it->m_tof * norm;
1666 weight = it->weight();
1667 errorSquared = it->errorSquared();
1668 lastTof = it->m_tof;
1669
1670 bin_end = next_bin(lastTof, bin_end);
1671 }
1672 }
1673
1674 // Put the last event in there too with the average TOF and summed weights and
1675 // squared errors.
1676 if (num == 1) {
1677 // last time-of-flight is the only one contributing
1678 out.emplace_back(lastTof, weight, errorSquared);
1679 } else if (num > 1) {
1680 out.emplace_back(totalTof / normalization, weight, errorSquared);
1681 }
1682
1683 // If you have over-allocated by more than 5%, reduce the size.
1684 size_t excess_limit = out.size() / 20;
1685 if ((out.capacity() - out.size()) > excess_limit) {
1686 out.shrink_to_fit();
1687 }
1688}
1689
1690template <class T>
1691inline void EventList::compressFatEventsHelper(const std::vector<T> &events, std::vector<WeightedEvent> &out,
1692 const double tolerance, const Types::Core::DateAndTime &timeStart,
1693 const double seconds) {
1694 // Clear the output. We can't know ahead of time how much space to reserve :(
1695 out.clear();
1696 // We will make a starting guess of 1/20th of the number of input events.
1697 out.reserve(events.size() / 20);
1698
1699 // The last TOF to which we are comparing.
1700 double lastTof = events.front().m_tof;
1701 // For getting an accurate average TOF
1702 double totalTof = 0;
1703
1704 // pulsetime bin information - stored as int nanoseconds because it
1705 // is the implementation type for DateAndTime object
1706 const int64_t pulsetimeStart = timeStart.totalNanoseconds();
1707 const auto pulsetimeDelta = static_cast<int64_t>(seconds * SEC_TO_NANO);
1708
1709 // pulsetime information
1710 std::vector<DateAndTime> pulsetimes; // all the times for new event
1711 std::vector<double> pulsetimeWeights;
1712
1713 // Carrying weight and error
1714 double weight = 0.;
1715 double errorSquared = 0.;
1716 double tofNormalization = 0.;
1717
1718 // Move up to first event that has a large enough pulsetime. This is just in case someone starts from after the
1719 // starttime of the run. It is expected that users will normally use the default which means this will only check the
1720 // first event.
1721 auto it = events.cbegin();
1722 for (; it != events.cend(); ++it) {
1723 if (it->m_pulsetime >= timeStart)
1724 break;
1725 }
1726
1727 if (it == events.cend())
1728 throw std::runtime_error("failed to find first pulse time in the events");
1729
1730 // bin if the pulses are histogrammed
1731 int64_t lastPulseBin = (it->m_pulsetime.totalNanoseconds() - pulsetimeStart) / pulsetimeDelta;
1732
1733 double bin_end = lastTof;
1734 double tof_min{0};
1735 std::function<bool(const double, const double)> compareTof;
1736 std::function<double(const double, double)> next_bin;
1737
1738 if (tolerance < 0) { // log
1739 // for log we do "less than" so that is matches the log binning of the Rebin algorithm
1740 compareTof = [](const double lhs, const double rhs) { return lhs < rhs; };
1741 next_bin = [tolerance](const double lastTof, double bin_end) {
1742 // advance the bin_end until we find the one that this next event falls into
1743 while (lastTof >= bin_end)
1744 bin_end = bin_end * (1 - tolerance);
1745 return bin_end;
1746 };
1747
1748 // get minimum Tof so that binning is consistent across all pulses
1749 const auto event_min = std::min_element(
1750 events.cbegin(), events.cend(), [](const auto &left, const auto &right) { return left.tof() < right.tof(); });
1751 bin_end = tof_min = event_min->tof();
1752
1753 if (tof_min < 0)
1754 throw std::runtime_error("compressEvents with log binning doesn't work with negative TOF");
1755
1756 // can't start at 0 as this will create an infinite loop
1757 if (tof_min == 0)
1758 bin_end = tof_min = fabs(tolerance);
1759
1760 } else { // linear
1761 // for linear we do "less than or equals" because that is how it was originally implemented
1762 compareTof = [](const double lhs, const double rhs) { return lhs <= rhs; };
1763 next_bin = [tolerance](const double lastTof, double) { return lastTof + tolerance; };
1764 }
1765
1766 // get first bin_end
1767 bin_end = next_bin(lastTof, bin_end);
1768
1769 // loop through events and accumulate weight
1770 for (; it != events.cend(); ++it) {
1771 const int64_t eventPulseBin = (it->m_pulsetime.totalNanoseconds() - pulsetimeStart) / pulsetimeDelta;
1772 if ((eventPulseBin <= lastPulseBin) && compareTof(it->m_tof, bin_end)) {
1773 // Carry the error and weight
1774 weight += it->weight();
1775 errorSquared += it->errorSquared();
1776 double norm = calcNorm(it->errorSquared());
1777 tofNormalization += norm;
1778 // Track the average tof
1779 totalTof += it->m_tof * norm;
1780 // Accumulate the pulse times
1781 pulsetimes.emplace_back(it->m_pulsetime);
1782 pulsetimeWeights.emplace_back(norm);
1783 } else {
1784 // We exceeded the tolerance
1785 if (!pulsetimes.empty()) {
1786 // Create a new event with the average TOF and summed weights and
1787 // squared errors. 1 event used doesn't need to average
1788 if (pulsetimes.size() == 1) {
1789 out.emplace_back(lastTof, pulsetimes.front(), weight, errorSquared);
1790 } else {
1791 out.emplace_back(totalTof / tofNormalization,
1792 Kernel::DateAndTimeHelpers::averageSorted(pulsetimes, pulsetimeWeights), weight,
1793 errorSquared);
1794 }
1795 }
1796 if (tolerance < 0 && eventPulseBin != lastPulseBin)
1797 // reset the bin_end for the new pulse bin
1798 bin_end = tof_min;
1799
1800 // Start a new combined object
1801 double norm = calcNorm(it->errorSquared());
1802 totalTof = it->m_tof * norm;
1803 weight = it->weight();
1804 errorSquared = it->errorSquared();
1805 tofNormalization = norm;
1806 lastTof = it->m_tof;
1807 lastPulseBin = eventPulseBin;
1808 pulsetimes.clear();
1809 pulsetimes.emplace_back(it->m_pulsetime);
1810 pulsetimeWeights.clear();
1811 pulsetimeWeights.emplace_back(norm);
1812
1813 bin_end = next_bin(lastTof, bin_end);
1814 }
1815 }
1816
1817 // Put the last event in there too.
1818 if (!pulsetimes.empty()) {
1819 // Create a new event with the average TOF and summed weights and
1820 // squared errors. 1 event used doesn't need to average
1821 if (pulsetimes.size() == 1) {
1822 out.emplace_back(lastTof, pulsetimes.front(), weight, errorSquared);
1823 } else {
1824 out.emplace_back(totalTof / tofNormalization,
1825 Kernel::DateAndTimeHelpers::averageSorted(pulsetimes, pulsetimeWeights), weight, errorSquared);
1826 }
1827 }
1828
1829 // If you have over-allocated by more than 5%, reduce the size.
1830 size_t excess_limit = out.size() / 20;
1831 if ((out.capacity() - out.size()) > excess_limit) {
1832 out.shrink_to_fit();
1833 }
1834}
1835
1836// --------------------------------------------------------------------------
1847 if (this->empty()) {
1848 // allocate memory in correct vector
1850 destination->weightedEventsNoTime = std::make_unique<std::vector<WeightedEventNoTime>>();
1851 } else {
1852 this->sortTof();
1853 switch (eventType) {
1854 case TOF:
1855 // if (parallel)
1856 // compressEventsParallelHelper(this->events,
1857 // destination->weightedEventsNoTime, tolerance);
1858 // else
1859 destination->weightedEventsNoTime = std::make_unique<std::vector<WeightedEventNoTime>>();
1861 break;
1862
1863 case WEIGHTED:
1864 // if (parallel)
1865 // compressEventsParallelHelper(this->weightedEvents,
1866 // destination->weightedEventsNoTime, tolerance);
1867 // else
1868 destination->weightedEventsNoTime = std::make_unique<std::vector<WeightedEventNoTime>>();
1870
1871 break;
1872
1873 case WEIGHTED_NOTIME:
1874 if (destination == this) {
1875 // Put results in a temp output
1876 auto out = std::make_unique<std::vector<WeightedEventNoTime>>();
1877 // if (parallel)
1878 // compressEventsParallelHelper(this->weightedEventsNoTime,
1879 // out,
1880 // tolerance);
1881 // else
1883 // Put it back
1884 this->weightedEventsNoTime.swap(out);
1885 } else {
1886 // if (parallel)
1887 // compressEventsParallelHelper(this->weightedEventsNoTime,
1888 // destination->weightedEventsNoTime, tolerance);
1889 // else
1890 destination->weightedEventsNoTime = std::make_unique<std::vector<WeightedEventNoTime>>();
1892 }
1893 break;
1894 }
1895 }
1896 // In all cases, you end up WEIGHTED_NOTIME.
1897 destination->eventType = WEIGHTED_NOTIME;
1898 // The sort is still valid!
1899 destination->order = TOF_SORT;
1900 // Empty out storage for vectors that are now unused.
1901 destination->clearUnused();
1902}
1903
1904template <class T>
1905inline void EventList::createWeightedEvents(std::vector<WeightedEventNoTime> &out, const std::vector<double> &tof,
1906 const std::vector<T> &weight, const std::vector<T> &error) {
1907 out.clear();
1908 for (size_t i = 0; i < weight.size(); ++i) {
1909 const auto errors = static_cast<float>(error[i]);
1910 if (errors > 0)
1911 out.emplace_back(tof[i], static_cast<float>(weight[i]), errors);
1912 }
1913}
1914
1915template <class T>
1916inline void EventList::processWeightedEvents(const std::vector<T> &events, std::vector<WeightedEventNoTime> &out,
1917 const std::shared_ptr<std::vector<double>> histogram_bin_edges,
1918 struct FindBin findBin) {
1919 const auto NUM_BINS = histogram_bin_edges->size() - 1;
1920 std::vector<double> tof(NUM_BINS, 0.);
1921 std::vector<double> normalization(NUM_BINS, 0.);
1922 std::vector<float> weight(NUM_BINS, 0.);
1923 std::vector<float> error(NUM_BINS, 0.);
1924 for (const auto &ev : events) {
1925 const auto &bin_optional = findBin(*histogram_bin_edges.get(), ev.m_tof, false);
1926 if (bin_optional) {
1927 const auto bin = bin_optional.value();
1928 const double norm = calcNorm(ev.m_errorSquared);
1929 tof[bin] += ev.m_tof * norm;
1930 normalization[bin] += norm;
1931 weight[bin] += ev.m_weight;
1932 error[bin] += ev.m_errorSquared;
1933 }
1934 }
1935
1936 // normalize TOFs
1937 std::transform(tof.begin(), tof.end(), normalization.begin(), tof.begin(), std::divides<double>());
1938
1939 createWeightedEvents(out, tof, weight, error);
1940}
1941
1943 const std::shared_ptr<std::vector<double>> histogram_bin_edges) {
1944 if (this->empty()) {
1945 // allocate memory in correct vector
1947 destination->weightedEventsNoTime = std::make_unique<std::vector<WeightedEventNoTime>>();
1948 } else {
1949 const auto NUM_BINS = histogram_bin_edges->size() - 1;
1950 const auto xmin = static_cast<double>(histogram_bin_edges->front());
1951
1952 auto findBin = FindBin(tolerance, xmin);
1953
1954 switch (eventType) {
1955 case TOF: {
1956 std::vector<double> tof(NUM_BINS, 0);
1957 std::vector<uint32_t> count(NUM_BINS, 0);
1958 for (const auto &ev : *this->events) {
1959 const auto &bin_optional = findBin(*histogram_bin_edges.get(), ev.m_tof, false);
1960 if (bin_optional) {
1961 const auto bin = bin_optional.value();
1962 count[bin]++;
1963 tof[bin] += ev.m_tof;
1964 }
1965 }
1966
1967 // average TOFs
1968 std::transform(tof.begin(), tof.end(), count.begin(), tof.begin(), std::divides<double>());
1969
1970 destination->weightedEventsNoTime = std::make_unique<std::vector<WeightedEventNoTime>>();
1972 break;
1973 }
1974
1975 case WEIGHTED: {
1976 destination->weightedEventsNoTime = std::make_unique<std::vector<WeightedEventNoTime>>();
1977 processWeightedEvents(*this->weightedEvents, *destination->weightedEventsNoTime, histogram_bin_edges, findBin);
1978 break;
1979 }
1980 case WEIGHTED_NOTIME:
1981 if (destination == this) {
1982 // Put results in a temp output
1983 auto out = std::make_unique<std::vector<WeightedEventNoTime>>();
1984 processWeightedEvents(*this->weightedEventsNoTime, *out, histogram_bin_edges, findBin);
1985 // Put it back
1986 this->weightedEventsNoTime.swap(out);
1987 } else {
1988 destination->weightedEventsNoTime = std::make_unique<std::vector<WeightedEventNoTime>>();
1989 processWeightedEvents(*this->weightedEventsNoTime, *destination->weightedEventsNoTime, histogram_bin_edges,
1990 findBin);
1991 }
1992 break;
1993 }
1994 }
1995
1996 // In all cases, you end up WEIGHTED_NOTIME.
1997 destination->eventType = WEIGHTED_NOTIME;
1998 // The result will be sorted
1999 destination->order = TOF_SORT;
2000 // Empty out storage for vectors that are now unused.
2001 destination->clearUnused();
2002}
2003
2004void EventList::compressFatEvents(const double tolerance, const Mantid::Types::Core::DateAndTime &timeStart,
2005 const double seconds, EventList *destination) {
2006 if (this->empty()) {
2007 // allocate memory in correct vector
2008 if (eventType != WEIGHTED)
2009 destination->weightedEvents = std::make_unique<std::vector<WeightedEvent>>();
2010 } else {
2011 switch (eventType) {
2012 case WEIGHTED_NOTIME:
2013 throw std::invalid_argument("Cannot compress events that do not have pulsetime");
2014 case TOF:
2015 this->sortPulseTimeTOFDelta(timeStart, seconds);
2016 destination->weightedEvents = std::make_unique<std::vector<WeightedEvent>>();
2017 compressFatEventsHelper(*this->events, *destination->weightedEvents, tolerance, timeStart, seconds);
2018 break;
2019 case WEIGHTED:
2020 this->sortPulseTimeTOFDelta(timeStart, seconds);
2021 if (destination == this) {
2022 // Put results in a temp output
2023 auto out = std::make_unique<std::vector<WeightedEvent>>();
2024 compressFatEventsHelper(*this->weightedEvents, *out, tolerance, timeStart, seconds);
2025 // Put it back
2026 this->weightedEvents.swap(out);
2027 } else {
2028 destination->weightedEvents = std::make_unique<std::vector<WeightedEvent>>();
2029 compressFatEventsHelper(*this->weightedEvents, *destination->weightedEvents, tolerance, timeStart, seconds);
2030 }
2031 break;
2032 }
2033 }
2034 // In all cases, you end up WEIGHTED_NOTIME.
2035 destination->eventType = WEIGHTED;
2036 // The sort order is pulsetimetof as we've compressed out the tolerance
2037 destination->order = PULSETIMETOF_SORT;
2038 // Empty out storage for vectors that are now unused.
2039 destination->clearUnused();
2040}
2041
2042// --------------------------------------------------------------------------
2052template <class T>
2053typename std::vector<T>::const_iterator static findFirstEvent(const std::vector<T> &events, T seek_tof) {
2054 return std::find_if_not(events.cbegin(), events.cend(), [seek_tof](const T &x) { return x < seek_tof; });
2055}
2056
2057// --------------------------------------------------------------------------
2067template <class T>
2068typename std::vector<T>::const_iterator EventList::findFirstPulseEvent(const std::vector<T> &events,
2069 const double seek_pulsetime) {
2070 auto itev = events.cbegin();
2071 auto itev_end = events.cend(); // cache for speed
2072
2073 // if tof < X[0], that means that you need to skip some events
2074 while ((itev != itev_end) && (static_cast<double>(itev->pulseTime().totalNanoseconds()) < seek_pulsetime))
2075 itev++;
2076 // Better fix would be to use a binary search instead of the linear one used
2077 // here.
2078 return itev;
2079}
2080
2081// --------------------------------------------------------------------------
2094template <class T>
2095typename std::vector<T>::const_iterator
2096EventList::findFirstTimeAtSampleEvent(const std::vector<T> &events, const double seek_time, const double &tofFactor,
2097 const double &tofOffset) const {
2098 auto itev = events.cbegin();
2099 auto itev_end = events.cend(); // cache for speed
2100
2101 // if tof < X[0], that means that you need to skip some events
2102 while ((itev != itev_end) &&
2103 (static_cast<double>(calculateCorrectedFullTime(*itev, tofFactor, tofOffset)) < seek_time))
2104 itev++;
2105 // Better fix would be to use a binary search instead of the linear one used
2106 // here.
2107 return itev;
2108}
2109
2110// --------------------------------------------------------------------------
2120template <class T> typename std::vector<T>::iterator static findFirstEvent(std::vector<T> &events, T seek_tof) {
2121 return std::find_if_not(events.begin(), events.end(), [seek_tof](const T &x) { return x < seek_tof; });
2122}
2123
2124// --------------------------------------------------------------------------
2134template <class T>
2135void EventList::histogramForWeightsHelper(const std::vector<T> &events, std::span<double const> X, MantidVec &Y,
2136 MantidVec &E) {
2137 // For slight speed=up.
2138 size_t x_size = X.size();
2139
2140 if (x_size <= 1) {
2141 // X was not set. Return an empty array.
2142 Y.resize(0, 0);
2143 return;
2144 }
2145
2146 // If the sizes are the same, then the "resize" command will NOT clear the original values.
2147 bool mustFill = (Y.size() == x_size - 1);
2148 // Clear the Y data, assign all to 0.
2149 Y.resize(x_size - 1, 0.0);
2150 // Clear the Error data, assign all to 0.
2151 // Note: Errors will be squared until the last step.
2152 E.resize(x_size - 1, 0.0);
2153
2154 if (mustFill) {
2155 // We must make sure the starting point is 0.0
2156 std::fill(Y.begin(), Y.end(), 0.0);
2157 std::fill(E.begin(), E.end(), 0.0);
2158 }
2159
2160 //---------------------- Histogram without weights
2161 //---------------------------------
2162
2163 // Do we even have any events to do?
2164 if (!events.empty()) {
2165 // Iterate through all events (sorted by tof)
2166 auto itev = findFirstEvent(events, T(X[0]));
2167 auto itev_end = events.cend();
2168 // The above can still take you to end() if no events above X[0], so check
2169 // again.
2170 if (itev == itev_end)
2171 return;
2172
2173 // Find the first bin
2174 size_t bin = 0;
2175 // The tof is greater the first bin boundary, so we need to find the first
2176 // bin
2177 double tof = itev->tof();
2178 while (bin < x_size - 1) {
2179 // Within range?
2180 if ((tof >= X[bin]) && (tof < X[bin + 1])) {
2181 // Add up the weight (convert to double before adding, to preserve
2182 // precision)
2183 Y[bin] += double(itev->m_weight);
2184 E[bin] += double(itev->m_errorSquared); // square of error
2185 break;
2186 }
2187 ++bin;
2188 }
2189 // Go to the next event, we've already binned this first one.
2190 ++itev;
2191
2192 // Keep going through all the events
2193 while ((itev != itev_end) && (bin < x_size - 1)) {
2194 tof = itev->tof();
2195 while (bin < x_size - 1) {
2196 // Within range? Since both events and X are sorted, they are going to
2197 // have
2198 // tof >= X[bin] because the previous event was.
2199 if (tof < X[bin + 1]) {
2200 // Add up the weight (convert to double before adding, to preserve
2201 // precision)
2202 Y[bin] += double(itev->m_weight);
2203 E[bin] += double(itev->m_errorSquared); // square of error
2204 break;
2205 }
2206 ++bin;
2207 }
2208 ++itev;
2209 }
2210 } // end if (there are any events to histogram)
2211
2212 // Now do the sqrt of all errors
2213 std::transform(E.cbegin(), E.cend(), E.begin(), static_cast<double (*)(double)>(sqrt));
2214}
2215
2216// --------------------------------------------------------------------------
2230template <class T>
2231void EventList::histogramForWeightsHelper(const std::vector<T> &events, const double step, std::span<double const> X,
2232 MantidVec &Y, MantidVec &E) {
2233 size_t x_size = X.size();
2234
2235 if (x_size <= 1) {
2236 // X was not set. Return an empty array.
2237 Y.resize(0, 0);
2238 return;
2239 }
2240
2241 // If the sizes are the same, then the "resize" command will NOT clear the original values.
2242 bool mustFill = (Y.size() == x_size - 1);
2243 Y.resize(x_size - 1, 0.0);
2244 E.resize(x_size - 1, 0.0);
2245 if (mustFill) {
2246 // We must make sure the starting point is 0.0
2247 std::fill(Y.begin(), Y.end(), 0.0);
2248 std::fill(E.begin(), E.end(), 0.0);
2249 }
2250
2251 if (events.empty())
2252 return;
2253
2254 const auto xmin = X.front();
2255 const auto xmax = X.back();
2256
2257 auto findBin = FindBin(step, xmin);
2258
2259 for (const T &ev : events) {
2260 const double tof = ev.tof();
2261 if (tof < xmin || tof >= xmax)
2262 continue;
2263
2264 std::optional<size_t> n_bin = findBin(X, tof, true);
2265
2266 if (n_bin) {
2267 Y[n_bin.value()] += ev.weight();
2268 E[n_bin.value()] += ev.errorSquared();
2269 }
2270 }
2271
2272 // Now do the sqrt of all errors
2273 std::transform(E.cbegin(), E.cend(), E.begin(), static_cast<double (*)(double)>(sqrt));
2274}
2275
2276// --------------------------------------------------------------------------
2287 bool skipError) const {
2288 // All types of weights need to be sorted by Pulse Time
2289 this->sortPulseTime();
2290
2291 switch (eventType) {
2292 case TOF:
2293 // Make the single ones
2295 if (!skipError)
2296 this->generateErrorsHistogram(Y, E);
2297 break;
2298
2299 case WEIGHTED:
2300 throw std::runtime_error("Cannot histogram by pulse time on Weighted "
2301 "Events currently"); // This could be supported.
2302
2303 case WEIGHTED_NOTIME:
2304 throw std::runtime_error("Cannot histogram by pulse time on Weighted Events NoTime");
2305 }
2306}
2307
2320 const double &tofFactor, const double &tofOffset, bool skipError) const {
2321 // All types of weights need to be sorted by time at sample
2322 this->sortTimeAtSample(tofFactor, tofOffset);
2323
2324 switch (eventType) {
2325 case TOF:
2326 // Make the single ones
2327 this->generateCountsHistogramTimeAtSample(X, Y, tofFactor, tofOffset);
2328 if (!skipError)
2329 this->generateErrorsHistogram(Y, E);
2330 break;
2331
2332 case WEIGHTED:
2333 throw std::runtime_error("Cannot histogram by time at sample on Weighted "
2334 "Events currently"); // This could be supported.
2335
2336 case WEIGHTED_NOTIME:
2337 throw std::runtime_error("Cannot histogram by time at sample on Weighted Events NoTime");
2338 }
2339}
2340
2341// --------------------------------------------------------------------------
2351void EventList::generateHistogram(std::span<double const> X, MantidVec &Y, MantidVec &E, bool skipError) const {
2352 // All types of weights need to be sorted by TOF
2353
2354 this->sortTof();
2355
2356 switch (eventType) {
2357 case TOF:
2358 // Make the single ones
2359 this->generateCountsHistogram(X, Y);
2360 if (!skipError)
2361 this->generateErrorsHistogram(Y, E);
2362 break;
2363
2364 case WEIGHTED:
2366 break;
2367
2368 case WEIGHTED_NOTIME:
2370 break;
2371 }
2372}
2373
2374// --------------------------------------------------------------------------
2389void EventList::generateHistogram(const double step, std::span<double const> X, MantidVec &Y, MantidVec &E,
2390 bool skipError) const {
2391 // if events are already sorted, use faster sorted histogram method
2392 if (isSortedByTof() || empty())
2393 return generateHistogram(X, Y, E, skipError);
2394
2395 switch (eventType) {
2396 case TOF:
2397 this->generateCountsHistogram(step, X, Y);
2398 if (!skipError)
2399 this->generateErrorsHistogram(Y, E);
2400 break;
2401
2402 case WEIGHTED:
2403 histogramForWeightsHelper(*this->weightedEvents, step, X, Y, E);
2404 break;
2405
2406 case WEIGHTED_NOTIME:
2408 break;
2409 }
2410}
2411
2412// --------------------------------------------------------------------------
2419void EventList::generateCountsHistogramPulseTime(std::span<double const> X, MantidVec &Y) const {
2420 // For slight speed=up.
2421 size_t x_size = X.size();
2422
2423 if (x_size <= 1) {
2424 // X was not set. Return an empty array.
2425 Y.resize(0, 0);
2426 return;
2427 }
2428
2429 // Sort the events by pulsetime
2430 this->sortPulseTime();
2431 // Clear the Y data, assign all to 0.
2432 Y.resize(x_size - 1, 0);
2433
2434 //---------------------- Histogram without weights
2435 //---------------------------------
2436
2437 if (!this->events->empty()) {
2438 // Iterate through all events (sorted by pulse time)
2439 auto itev = findFirstPulseEvent(*this->events, X[0]);
2440 auto itev_end = events->cend(); // cache for speed
2441 // The above can still take you to end() if no events above X[0], so check
2442 // again.
2443 if (itev == itev_end)
2444 return;
2445
2446 // Find the first bin
2447 size_t bin = 0;
2448
2449 // The tof is greater the first bin boundary, so we need to find the first
2450 // bin
2451 double pulsetime = static_cast<double>(itev->pulseTime().totalNanoseconds());
2452 while (bin < x_size - 1) {
2453 // Within range?
2454 if ((pulsetime >= X[bin]) && (pulsetime < X[bin + 1])) {
2455 Y[bin]++;
2456 break;
2457 }
2458 ++bin;
2459 }
2460 // Go to the next event, we've already binned this first one.
2461 ++itev;
2462
2463 // Keep going through all the events
2464 while ((itev != itev_end) && (bin < x_size - 1)) {
2465 pulsetime = static_cast<double>(itev->pulseTime().totalNanoseconds());
2466 while (bin < x_size - 1) {
2467 // Within range?
2468 if ((pulsetime >= X[bin]) && (pulsetime < X[bin + 1])) {
2469 Y[bin]++;
2470 break;
2471 }
2472 ++bin;
2473 }
2474 ++itev;
2475 }
2476 } // end if (there are any events to histogram)
2477}
2478
2493void EventList::generateCountsHistogramPulseTime(const double &xMin, const double &xMax, MantidVec &Y,
2494 const double TOF_min, const double TOF_max) const {
2495
2496 if (this->events->empty())
2497 return;
2498
2499 size_t nBins = Y.size();
2500
2501 if (nBins == 0)
2502 return;
2503
2504 double step = (xMax - xMin) / static_cast<double>(nBins);
2505
2506 for (const TofEvent &ev : *this->events) {
2507 double pulsetime = static_cast<double>(ev.pulseTime().totalNanoseconds());
2508 if (pulsetime < xMin || pulsetime >= xMax)
2509 continue;
2510 if (ev.tof() < TOF_min || ev.tof() >= TOF_max)
2511 continue;
2512
2513 auto n_bin = static_cast<size_t>((pulsetime - xMin) / step);
2514 Y[n_bin]++;
2515 }
2516}
2517
2518// --------------------------------------------------------------------------
2527void EventList::generateCountsHistogramTimeAtSample(std::span<double const> X, MantidVec &Y, const double &tofFactor,
2528 const double &tofOffset) const {
2529 // For slight speed=up.
2530 const size_t x_size = X.size();
2531
2532 if (x_size <= 1) {
2533 // X was not set. Return an empty array.
2534 Y.resize(0, 0);
2535 return;
2536 }
2537
2538 // Sort the events by pulsetime
2539 this->sortTimeAtSample(tofFactor, tofOffset);
2540 // Clear the Y data, assign all to 0.
2541 Y.resize(x_size - 1, 0);
2542
2543 //---------------------- Histogram without weights
2544 //---------------------------------
2545
2546 if (!this->events->empty()) {
2547 // Iterate through all events (sorted by pulse time)
2548 auto itev = findFirstTimeAtSampleEvent(*this->events, X[0], tofFactor, tofOffset);
2549 std::vector<TofEvent>::const_iterator itev_end = events->end(); // cache for speed
2550 // The above can still take you to end() if no events above X[0], so check
2551 // again.
2552 if (itev == itev_end)
2553 return;
2554
2555 // Find the first bin
2556 size_t bin = 0;
2557
2558 auto tAtSample = static_cast<double>(calculateCorrectedFullTime(*itev, tofFactor, tofOffset));
2559 while (bin < x_size - 1) {
2560 // Within range?
2561 if ((tAtSample >= X[bin]) && (tAtSample < X[bin + 1])) {
2562 Y[bin]++;
2563 break;
2564 }
2565 ++bin;
2566 }
2567 // Go to the next event, we've already binned this first one.
2568 ++itev;
2569
2570 // Keep going through all the events
2571 while ((itev != itev_end) && (bin < x_size - 1)) {
2572 tAtSample = static_cast<double>(calculateCorrectedFullTime(*itev, tofFactor, tofOffset));
2573 while (bin < x_size - 1) {
2574 // Within range?
2575 if ((tAtSample >= X[bin]) && (tAtSample < X[bin + 1])) {
2576 Y[bin]++;
2577 break;
2578 }
2579 ++bin;
2580 }
2581 ++itev;
2582 }
2583 } // end if (there are any events to histogram)
2584}
2585
2586// --------------------------------------------------------------------------
2592void EventList::generateCountsHistogram(std::span<double const> X, MantidVec &Y) const {
2593 // For slight speed=up.
2594 size_t x_size = X.size();
2595
2596 if (x_size <= 1) {
2597 // X was not set. Return an empty array.
2598 Y.resize(0, 0);
2599 return;
2600 }
2601
2602 // Sort the events by tof
2603 this->sortTof();
2604 // Clear the Y data, assign all to 0.
2605 Y.resize(x_size - 1, 0);
2606
2607 //---------------------- Histogram without weights
2608 //---------------------------------
2609
2610 // Do we even have any events to do?
2611 if (!this->events->empty()) {
2612 // Iterate through all events (sorted by tof) placing them in the correct
2613 // bin.
2614 auto itev = findFirstEvent(*this->events, TofEvent(X[0]));
2615 const auto itend = this->events->end();
2616 // Go through all the events,
2617 for (auto itx = X.begin(); itev != itend; ++itev) {
2618 const double tof = itev->tof();
2619 itx = std::find_if(itx, X.end(), [tof](const double x) { return tof < x; });
2620 if (itx == X.end()) {
2621 break;
2622 }
2623 const auto bin = static_cast<size_t>(std::max(std::distance(X.begin(), itx) - 1, std::ptrdiff_t{0}));
2624 ++Y[bin];
2625 }
2626 } // end if (there are any events to histogram)
2627}
2628
2639std::optional<size_t> EventList::findLinearBin(std::span<double const> X, const double tof, const double divisor,
2640 const double offset, const bool findExact) {
2641 const auto bin = static_cast<size_t>(tof * divisor - offset);
2642 if (bin >= X.size())
2643 return std::nullopt;
2644 else if (findExact)
2645 return findExactBin(X, tof, bin);
2646 else
2647 return bin;
2648}
2649
2668std::optional<size_t> EventList::findLogBin(std::span<double const> X, const double tof, const double divisor,
2669 const double offset, const bool findExact) {
2670 const auto bin = static_cast<size_t>(log(tof) * divisor - offset);
2671 if (bin >= X.size())
2672 return std::nullopt;
2673 else if (findExact)
2674 return findExactBin(X, tof, bin);
2675 else
2676 return bin;
2677}
2678
2686size_t EventList::findExactBin(std::span<double const> X, const double tof, const size_t n_bin) {
2687 // is tof slower than suggested bin
2688 auto tof_of_bin = X.begin() + n_bin; // boundary suggested
2689 if (tof < *tof_of_bin)
2690 return std::move(n_bin - 1);
2691
2692 // is tof higher than suggested bin
2693 ++tof_of_bin; // move to next boundary
2694 if (tof >= *tof_of_bin)
2695 return std::move(n_bin + 1);
2696
2697 // tof is in the bin
2698 return std::move(n_bin);
2699}
2700
2701// --------------------------------------------------------------------------
2712void EventList::generateCountsHistogram(const double step, std::span<double const> X, MantidVec &Y) const {
2713 // For slight speed=up.
2714 size_t x_size = X.size();
2715
2716 if (x_size <= 1) {
2717 // X was not set. Return an empty array.
2718 Y.resize(0, 0);
2719 return;
2720 }
2721
2722 // If the sizes are the same, then the "resize" command will NOT clear the original values.
2723 bool mustFill = (Y.size() == x_size - 1);
2724 // Clear the Y data, assign all to 0.
2725 Y.resize(x_size - 1, 0);
2726 if (mustFill) // starting point is no counts
2727 std::fill(Y.begin(), Y.end(), 0.0);
2728
2729 // Do we even have any events to do?
2730 if (this->events->empty())
2731 return;
2732
2733 const auto xmin = X.front();
2734 const auto xmax = X.back();
2735
2736 auto findBin = FindBin(step, xmin);
2737
2738 for (const TofEvent &ev : *this->events) {
2739 const double tof = ev.tof();
2740 if (tof < xmin || tof >= xmax)
2741 continue;
2742
2743 const std::optional<size_t> n_bin = findBin(X, tof, true);
2744
2745 if (n_bin)
2746 Y[n_bin.value()]++;
2747 }
2748}
2749
2750// --------------------------------------------------------------------------
2758void EventList::generateErrorsHistogram(std::span<double const> Y, MantidVec &E) const {
2759 // Fill the vector for the errors, containing sqrt(count)
2760 E.resize(Y.size(), 0);
2761
2762 // windows can get confused about std::sqrt
2763 std::transform(Y.begin(), Y.end(), E.begin(), static_cast<double (*)(double)>(sqrt));
2764
2765} //----------------------------------------------------------------------------------
2766
2777template <class T>
2778void EventList::integrateHelper(std::vector<T> &events, const double minX, const double maxX, const bool entireRange,
2779 double &sum, double &error) {
2780 sum = 0;
2781 error = 0;
2782 // Nothing in the list?
2783 if (events.empty())
2784 return;
2785
2786 // Iterators for limits - whole range by default
2787 auto lowit = events.cbegin();
2788 auto highit = events.cend();
2789
2790 // But maybe we don't want the entire range?
2791 if (!entireRange) {
2792 // If a silly range was given, return 0.
2793 if (maxX < minX)
2794 return;
2795
2796 // If the first element is lower that the xmin then search for new lowit
2797 if (lowit->tof() < minX)
2798 lowit = std::lower_bound(events.cbegin(), events.cend(), minX);
2799 // If the last element is higher that the xmax then search for new lowit
2800 if ((highit - 1)->tof() > maxX) {
2801 highit = std::upper_bound(lowit, events.cend(), T(maxX));
2802 }
2803 }
2804
2805 // Sum up all the weights
2806 for (auto it = lowit; it != highit; ++it) {
2807 sum += it->weight();
2808 error += it->errorSquared();
2809 }
2810 error = std::sqrt(error);
2811}
2812
2813// --------------------------------------------------------------------------
2822double EventList::integrate(const double minX, const double maxX, const bool entireRange) const {
2823 double sum(0), error(0);
2824 integrate(minX, maxX, entireRange, sum, error);
2825 return sum;
2826}
2827
2837void EventList::integrate(const double minX, const double maxX, const bool entireRange, double &sum,
2838 double &error) const {
2839 sum = 0;
2840 error = 0;
2841 if (!entireRange) {
2842 // The event list must be sorted by TOF!
2843 this->sortTof();
2844 }
2845
2846 // Convert the list
2847 switch (eventType) {
2848 case TOF:
2849 integrateHelper(*this->events, minX, maxX, entireRange, sum, error);
2850 break;
2851 case WEIGHTED:
2852 integrateHelper(*this->weightedEvents, minX, maxX, entireRange, sum, error);
2853 break;
2854 case WEIGHTED_NOTIME:
2855 integrateHelper(*this->weightedEventsNoTime, minX, maxX, entireRange, sum, error);
2856 break;
2857 default:
2858 throw std::runtime_error("EventList: invalid event type value was found.");
2859 }
2860}
2861
2862// ==============================================================================================
2863// ----------- Conversion Functions (changing tof values)
2864// ---------------------------------------
2865// ==============================================================================================
2866
2873void EventList::convertTof(std::function<double(double)> func, const int sorting) {
2874 // fix the histogram parameter. Changing X invalidates any histogram cached in the MRU, so drop
2875 // it first -- this is what the deprecated dataX() used to do implicitly.
2876 if (mru)
2877 mru->deleteIndex(this);
2878 auto &x = mutableX();
2879 transform(x.cbegin(), x.cend(), x.begin(), func);
2880
2881 // do nothing if sorting > 0
2882 if (sorting == 0) {
2883 this->setSortOrder(UNSORTED);
2884 } else if ((sorting < 0) && (this->getSortType() == TOF_SORT)) {
2885 this->reverse();
2886 }
2887
2888 if (this->getNumberEvents() == 0)
2889 return;
2890
2891 // Convert the list
2892 switch (eventType) {
2893 case TOF:
2894 this->convertTofHelper(*this->events, func);
2895 break;
2896 case WEIGHTED:
2897 this->convertTofHelper(*this->weightedEvents, func);
2898 break;
2899 case WEIGHTED_NOTIME:
2900 this->convertTofHelper(*this->weightedEventsNoTime, func);
2901 break;
2902 }
2903}
2904
2909template <class T> void EventList::convertTofHelper(std::vector<T> &events, const std::function<double(double)> &func) {
2910 // iterate through all events
2911 for (auto &ev : events)
2912 ev.m_tof = func(ev.m_tof);
2913}
2914
2915// --------------------------------------------------------------------------
2921void EventList::convertTof(const double factor, const double offset) {
2922 // fix the histogram parameter
2923 auto &x = mutableX();
2924 x *= factor;
2925 x += offset;
2926
2927 if ((factor < 0.) && (this->getSortType() == TOF_SORT))
2928 this->reverse();
2929
2930 if (this->getNumberEvents() == 0)
2931 return;
2932
2933 // Convert the list
2934 switch (eventType) {
2935 case TOF:
2936 this->convertTofHelper(*this->events, factor, offset);
2937 break;
2938 case WEIGHTED:
2939 this->convertTofHelper(*this->weightedEvents, factor, offset);
2940 break;
2941 case WEIGHTED_NOTIME:
2942 this->convertTofHelper(*this->weightedEventsNoTime, factor, offset);
2943 break;
2944 }
2945}
2946
2947// --------------------------------------------------------------------------
2956template <class T> void EventList::convertTofHelper(std::vector<T> &events, const double factor, const double offset) {
2957 // iterate through all events
2958 for (auto &event : events) {
2959 event.m_tof = event.m_tof * factor + offset;
2960 }
2961}
2962
2963// --------------------------------------------------------------------------
2970void EventList::scaleTof(const double factor) { this->convertTof(factor, 0.0); }
2971
2972// --------------------------------------------------------------------------
2977void EventList::addTof(const double offset) { this->convertTof(1.0, offset); }
2978
2979// --------------------------------------------------------------------------
2985template <class T> void EventList::addPulsetimeHelper(std::vector<T> &events, const double seconds) {
2986 // iterate through all events
2987 for (auto &event : events) {
2988 event.m_pulsetime += seconds;
2989 }
2990}
2991
2998template <class T> void EventList::addPulsetimesHelper(std::vector<T> &events, const std::vector<double> &seconds) {
2999 auto eventIterEnd{events.end()};
3000 auto secondsIter{seconds.cbegin()};
3001 for (auto eventIter = events.begin(); eventIter < eventIterEnd; ++eventIter, ++secondsIter) {
3002 eventIter->m_pulsetime += *secondsIter;
3003 }
3004}
3005
3006// --------------------------------------------------------------------------
3011void EventList::addPulsetime(const double seconds) {
3012 if (this->getNumberEvents() == 0)
3013 return;
3014
3015 // Convert the list
3016 switch (eventType) {
3017 case TOF:
3018 this->addPulsetimeHelper(*this->events, seconds);
3019 break;
3020 case WEIGHTED:
3021 this->addPulsetimeHelper(*this->weightedEvents, seconds);
3022 break;
3023 case WEIGHTED_NOTIME:
3024 throw std::runtime_error("EventList::addPulsetime() called on an event "
3025 "list with no pulse times. You must call this "
3026 "algorithm BEFORE CompressEvents.");
3027 break;
3028 }
3029}
3030
3031// --------------------------------------------------------------------------
3036void EventList::addPulsetimes(const std::vector<double> &seconds) {
3037 if (this->getNumberEvents() == 0)
3038 return;
3039 if (this->getNumberEvents() != seconds.size()) {
3040 throw std::runtime_error("");
3041 }
3042
3043 // Convert the list
3044 switch (eventType) {
3045 case TOF:
3046 this->addPulsetimesHelper(*this->events, seconds);
3047 break;
3048 case WEIGHTED:
3049 this->addPulsetimesHelper(*this->weightedEvents, seconds);
3050 break;
3051 case WEIGHTED_NOTIME:
3052 throw std::runtime_error("EventList::addPulsetime() called on an event "
3053 "list with no pulse times. You must call this "
3054 "algorithm BEFORE CompressEvents.");
3055 break;
3056 }
3057}
3058
3059// --------------------------------------------------------------------------
3067template <class T>
3068std::size_t EventList::maskTofHelper(std::vector<T> &events, const double tofMin, const double tofMax) {
3069 // quick checks to make sure that the masking range is even in the data
3070 if (tofMin > events.crbegin()->tof())
3071 return 0;
3072 if (tofMax < events.cbegin()->tof())
3073 return 0;
3074
3075 // Find the index of the first tofMin
3076 auto it_first = std::lower_bound(events.begin(), events.end(), tofMin);
3077 if ((it_first != events.end()) && (it_first->tof() < tofMax)) {
3078 // Something was found
3079 // Look for the first one > tofMax
3080 auto it_last = std::upper_bound(it_first, events.end(), T(tofMax));
3081
3082 if (it_first >= it_last) {
3083 throw std::runtime_error("Event filter is all messed up"); // TODO
3084 }
3085
3086 size_t tmp = std::size_t(std::distance(it_first, it_last));
3087 // it_last will either be at the end (if not found) or before it.
3088 // Erase this range from the vector
3089 events.erase(it_first, it_last);
3090
3091 // Done! Sorting is still valid, no need to redo.
3092 return tmp; //(it_last - it_first); the iterators get invalid after erase
3093 }
3094 return 0; // didn't remove any events
3095}
3096
3097// --------------------------------------------------------------------------
3104void EventList::maskTof(const double tofMin, const double tofMax) {
3105 if (tofMax <= tofMin)
3106 throw std::runtime_error("EventList::maskTof: tofMax must be > tofMin");
3107
3108 // don't do anything with an emply list
3109 if (this->getNumberEvents() == 0)
3110 return;
3111
3112 // Start by sorting by tof
3113 this->sortTof();
3114
3115 // Convert the list
3116 size_t numOrig = 0;
3117 size_t numDel = 0;
3118 switch (eventType) {
3119 case TOF:
3120 numOrig = this->events->size();
3121 numDel = this->maskTofHelper(*this->events, tofMin, tofMax);
3122 break;
3123 case WEIGHTED:
3124 numOrig = this->weightedEvents->size();
3125 numDel = this->maskTofHelper(*this->weightedEvents, tofMin, tofMax);
3126 break;
3127 case WEIGHTED_NOTIME:
3128 numOrig = this->weightedEventsNoTime->size();
3129 numDel = this->maskTofHelper(*this->weightedEventsNoTime, tofMin, tofMax);
3130 break;
3131 }
3132
3133 if (numDel >= numOrig)
3134 this->clear(false);
3135}
3136
3137// --------------------------------------------------------------------------
3144template <class T> std::size_t EventList::maskConditionHelper(std::vector<T> &events, const std::vector<bool> &mask) {
3145
3146 // runs through the two synchronized vectors and delete elements
3147 // for condition false
3148 auto itm = std::find(mask.begin(), mask.end(), false);
3149 auto first = events.begin() + (itm - mask.begin());
3150
3151 if (itm != mask.end()) {
3152 for (auto ite = first; ++ite != events.end() && ++itm != mask.end();) {
3153 if (*itm != false) {
3154 *first++ = std::move(*ite);
3155 }
3156 }
3157 }
3158
3159 const auto n = static_cast<size_t>(events.end() - first);
3160 if (n != 0)
3161 events.erase(first, events.end());
3162
3163 return n;
3164}
3165
3166// --------------------------------------------------------------------------
3172void EventList::maskCondition(const std::vector<bool> &mask) {
3173
3174 // mask size must match the number of events
3175 if (this->getNumberEvents() != mask.size())
3176 throw std::runtime_error("EventList::maskTof: tofMax must be > tofMin");
3177
3178 // don't do anything with an emply list
3179 if (this->getNumberEvents() == 0)
3180 return;
3181
3182 // Convert the list
3183 size_t numOrig = 0;
3184 size_t numDel = 0;
3185 switch (eventType) {
3186 case TOF:
3187 numOrig = this->events->size();
3188 numDel = this->maskConditionHelper(*this->events, mask);
3189 break;
3190 case WEIGHTED:
3191 numOrig = this->weightedEvents->size();
3192 numDel = this->maskConditionHelper(*this->weightedEvents, mask);
3193 break;
3194 case WEIGHTED_NOTIME:
3195 numOrig = this->weightedEventsNoTime->size();
3196 numDel = this->maskConditionHelper(*this->weightedEventsNoTime, mask);
3197 break;
3198 }
3199
3200 if (numDel >= numOrig)
3201 this->clear(false);
3202}
3203
3204// --------------------------------------------------------------------------
3210template <class T> void EventList::getTofsHelper(const std::vector<T> &events, std::vector<double> &tofs) {
3211 tofs.clear();
3212 for (auto itev = events.cbegin(); itev != events.cend(); ++itev)
3213 tofs.emplace_back(itev->m_tof);
3214}
3215
3219void EventList::getTofs(std::vector<double> &tofs) const {
3220 // Set the capacity of the vector to avoid multiple resizes
3221 tofs.reserve(this->getNumberEvents());
3222
3223 // Convert the list
3224 switch (eventType) {
3225 case TOF:
3226 this->getTofsHelper(*this->events, tofs);
3227 break;
3228 case WEIGHTED:
3229 this->getTofsHelper(*this->weightedEvents, tofs);
3230 break;
3231 case WEIGHTED_NOTIME:
3232 this->getTofsHelper(*this->weightedEventsNoTime, tofs);
3233 break;
3234 }
3235}
3236
3241std::vector<double> EventList::getTofs() const {
3242 std::vector<double> tofs;
3243 this->getTofs(tofs);
3244 return tofs;
3245}
3246
3247// --------------------------------------------------------------------------
3253template <class T> void EventList::getWeightsHelper(const std::vector<T> &events, std::vector<double> &weights) {
3254 weights.clear();
3255 weights.reserve(events.size());
3256 std::transform(events.cbegin(), events.cend(), std::back_inserter(weights),
3257 [](const auto &event) { return event.weight(); });
3258}
3259
3263void EventList::getWeights(std::vector<double> &weights) const {
3264 // Set the capacity of the vector to avoid multiple resizes
3265 weights.reserve(this->getNumberEvents());
3266
3267 // Convert the list
3268 switch (eventType) {
3269 case WEIGHTED:
3270 this->getWeightsHelper(*this->weightedEvents, weights);
3271 break;
3272 case WEIGHTED_NOTIME:
3273 this->getWeightsHelper(*this->weightedEventsNoTime, weights);
3274 break;
3275 default:
3276 // not a weighted event type, return 1.0 for all.
3277 weights.assign(this->getNumberEvents(), 1.0);
3278 break;
3279 }
3280}
3281
3286std::vector<double> EventList::getWeights() const {
3287 std::vector<double> weights;
3288 this->getWeights(weights);
3289 return weights;
3290}
3291
3292// --------------------------------------------------------------------------
3298template <class T>
3299void EventList::getWeightErrorsHelper(const std::vector<T> &events, std::vector<double> &weightErrors) {
3300 weightErrors.clear();
3301 weightErrors.reserve(events.size());
3302 std::transform(events.cbegin(), events.cend(), std::back_inserter(weightErrors),
3303 [](const auto &event) { return event.error(); });
3304}
3305
3309void EventList::getWeightErrors(std::vector<double> &weightErrors) const {
3310 // Set the capacity of the vector to avoid multiple resizes
3311 weightErrors.reserve(this->getNumberEvents());
3312
3313 // Convert the list
3314 switch (eventType) {
3315 case WEIGHTED:
3316 this->getWeightErrorsHelper(*this->weightedEvents, weightErrors);
3317 break;
3318 case WEIGHTED_NOTIME:
3319 this->getWeightErrorsHelper(*this->weightedEventsNoTime, weightErrors);
3320 break;
3321 default:
3322 // not a weighted event type, return 1.0 for all.
3323 weightErrors.assign(this->getNumberEvents(), 1.0);
3324 break;
3325 }
3326}
3327
3332std::vector<double> EventList::getWeightErrors() const {
3333 std::vector<double> weightErrors;
3334 this->getWeightErrors(weightErrors);
3335 return weightErrors;
3336}
3337
3343template <typename UnaryOperation>
3344std::vector<DateAndTime> EventList::eventTimesCalculator(const UnaryOperation &timesCalc) const {
3345 std::vector<DateAndTime> times;
3346 switch (eventType) {
3347 case TOF:
3348 times.reserve(events->size());
3349 std::transform(events->cbegin(), events->cend(), std::back_inserter(times), timesCalc);
3350 break;
3351 case WEIGHTED:
3352 times.reserve(weightedEvents->size());
3353 std::transform(weightedEvents->cbegin(), weightedEvents->cend(), std::back_inserter(times), timesCalc);
3354 break;
3355 case WEIGHTED_NOTIME:
3356 times.reserve(weightedEventsNoTime->size());
3357 std::transform(weightedEventsNoTime->cbegin(), weightedEventsNoTime->cend(), std::back_inserter(times), timesCalc);
3358 break;
3359 }
3360 return times;
3361}
3362
3367std::vector<Mantid::Types::Core::DateAndTime> EventList::getPulseTimes() const {
3368 auto timeCalc = [](const auto &event) { return event.pulseTime(); };
3369 return eventTimesCalculator(timeCalc);
3370}
3371
3373std::vector<DateAndTime> EventList::getPulseTOFTimes() const {
3374 auto timeCalc = [](const auto &event) { return event.pulseTOFTime(); };
3375 return eventTimesCalculator(timeCalc);
3376}
3377
3382std::vector<DateAndTime> EventList::getPulseTOFTimesAtSample(const double &factor, const double &shift) const {
3383 auto timeCalc = [factor, shift](const auto &event) { return event.pulseTOFTimeAtSample(factor, shift); };
3384 return eventTimesCalculator(timeCalc);
3385}
3386
3387// --------------------------------------------------------------------------
3388
3389namespace { // anonymous namespace
3390template <class T> double getTofMinimumHelper(const std::vector<T> &events) {
3391 const auto result = std::min_element(events.cbegin(), events.cend(),
3392 [](const auto &left, const auto &right) { return left.tof() < right.tof(); });
3393 return result->tof();
3394}
3395
3396template <class T> double getTofMaximumHelper(const std::vector<T> &events) {
3397 const auto result = std::max_element(events.cbegin(), events.cend(),
3398 [](const auto &left, const auto &right) { return left.tof() < right.tof(); });
3399 return result->tof();
3400}
3401} // anonymous namespace
3402
3406double EventList::getTofMin() const {
3407 // set up as the maximum available double
3408 double tMin = std::numeric_limits<double>::max();
3409
3410 // no events is a soft error
3411 if (this->empty())
3412 return tMin;
3413
3414 // when events are ordered by tof just need the first value
3415 if (this->order == TOF_SORT) {
3416 switch (eventType) {
3417 case TOF:
3418 return this->events->front().tof();
3419 case WEIGHTED:
3420 return this->weightedEvents->front().tof();
3421 case WEIGHTED_NOTIME:
3422 return this->weightedEventsNoTime->front().tof();
3423 }
3424 }
3425
3426 // now we are stuck with a linear search
3427 switch (eventType) {
3428 case TOF: {
3429 tMin = getTofMinimumHelper(*this->events);
3430 break;
3431 }
3432 case WEIGHTED: {
3433 tMin = getTofMinimumHelper(*this->weightedEvents);
3434 break;
3435 }
3436 case WEIGHTED_NOTIME: {
3437 tMin = getTofMinimumHelper(*this->weightedEventsNoTime);
3438 break;
3439 }
3440 }
3441
3442 return tMin;
3443}
3444
3448double EventList::getTofMax() const {
3449 // set up as the minimum available double
3450 double tMax = std::numeric_limits<double>::lowest();
3451
3452 // no events is a soft error
3453 if (this->empty())
3454 return tMax;
3455
3456 // when events are ordered by tof just need the first value
3457 if (this->order == TOF_SORT) {
3458 switch (eventType) {
3459 case TOF:
3460 return this->events->back().tof();
3461 case WEIGHTED:
3462 return this->weightedEvents->back().tof();
3463 case WEIGHTED_NOTIME:
3464 return this->weightedEventsNoTime->back().tof();
3465 }
3466 }
3467
3468 // now we are stuck with a linear search
3469 switch (eventType) {
3470 case TOF: {
3471 tMax = getTofMaximumHelper(*this->events);
3472 break;
3473 }
3474 case WEIGHTED: {
3475 tMax = getTofMaximumHelper(*this->weightedEvents);
3476 break;
3477 }
3478 case WEIGHTED_NOTIME: {
3479 tMax = getTofMaximumHelper(*this->weightedEventsNoTime);
3480 break;
3481 }
3482 }
3483
3484 return tMax;
3485}
3486
3487// --------------------------------------------------------------------------
3488namespace { // anonymous namespace
3489template <class T> DateAndTime getPulseMinimumHelper(const std::vector<T> &events) {
3490 const auto result = std::min_element(events.cbegin(), events.cend(), [](const auto &left, const auto &right) {
3491 return left.pulseTime() < right.pulseTime();
3492 });
3493 return result->pulseTime();
3494}
3495
3496template <class T> DateAndTime getPulseMaximumHelper(const std::vector<T> &events) {
3497 const auto result = std::max_element(events.cbegin(), events.cend(), [](const auto &left, const auto &right) {
3498 return left.pulseTime() < right.pulseTime();
3499 });
3500 return result->pulseTime();
3501}
3502} // anonymous namespace
3503
3507DateAndTime EventList::getPulseTimeMin() const {
3508 // no events is a soft error
3509 if (this->empty())
3510 return DateAndTime::maximum();
3511
3512 // when events are ordered by pulse time just need the first value
3513 if (this->order == PULSETIME_SORT) {
3514 switch (eventType) {
3515 case TOF:
3516 return this->events->front().pulseTime();
3517 case WEIGHTED:
3518 return this->weightedEvents->front().pulseTime();
3519 case WEIGHTED_NOTIME:
3520 return this->weightedEventsNoTime->front().pulseTime();
3521 }
3522 }
3523
3524 // now we are stuck with a linear search
3525 switch (eventType) {
3526 case TOF:
3527 return getPulseMinimumHelper(*this->events);
3528 case WEIGHTED:
3529 return getPulseMinimumHelper(*this->weightedEvents);
3530 case WEIGHTED_NOTIME:
3531 return getPulseMinimumHelper(*this->weightedEventsNoTime);
3532 }
3533
3534 return DateAndTime::maximum();
3535}
3536
3540DateAndTime EventList::getPulseTimeMax() const {
3541 // no events is a soft error
3542 if (this->empty())
3543 return DateAndTime::minimum();
3544
3545 // when events are ordered by pulse time just need the first value
3546 if (this->order == PULSETIME_SORT) {
3547 switch (eventType) {
3548 case TOF:
3549 return this->events->back().pulseTime();
3550 case WEIGHTED:
3551 return this->weightedEvents->back().pulseTime();
3552 case WEIGHTED_NOTIME:
3553 return this->weightedEventsNoTime->back().pulseTime();
3554 }
3555 }
3556
3557 // now we are stuck with a linear search
3558 switch (eventType) {
3559 case TOF:
3560 return getPulseMaximumHelper(*this->events);
3561 case WEIGHTED:
3562 return getPulseMaximumHelper(*this->weightedEvents);
3563 case WEIGHTED_NOTIME:
3564 return getPulseMaximumHelper(*this->weightedEventsNoTime);
3565 }
3566
3567 return DateAndTime::minimum();
3568}
3569
3570void EventList::getPulseTimeMinMax(Mantid::Types::Core::DateAndTime &tMin,
3571 Mantid::Types::Core::DateAndTime &tMax) const {
3572 // set up as the minimum available date time.
3573 tMax = DateAndTime::minimum();
3574 tMin = DateAndTime::maximum();
3575
3576 // no events is a soft error
3577 if (this->empty())
3578 return;
3579
3580 // when events are ordered by pulse time just need the first/last values
3581 if (this->order == PULSETIME_SORT) {
3582 switch (eventType) {
3583 case TOF:
3584 tMin = this->events->front().pulseTime();
3585 tMax = this->events->back().pulseTime();
3586 return;
3587 case WEIGHTED:
3588 tMin = this->weightedEvents->front().pulseTime();
3589 tMax = this->weightedEvents->back().pulseTime();
3590 return;
3591 case WEIGHTED_NOTIME:
3592 tMin = this->weightedEventsNoTime->front().pulseTime();
3593 tMax = this->weightedEventsNoTime->back().pulseTime();
3594 return;
3595 }
3596 }
3597
3598 // now we are stuck with a linear search
3599 // could this be done more efficiently than using ->at?
3600 size_t numEvents = this->getNumberEvents();
3601 DateAndTime temp = tMax; // start with the smallest possible value
3602 for (size_t i = 0; i < numEvents; i++) {
3603 switch (eventType) {
3604 case TOF:
3605 temp = this->events->at(i).pulseTime();
3606 break;
3607 case WEIGHTED:
3608 temp = this->weightedEvents->at(i).pulseTime();
3609 break;
3610 case WEIGHTED_NOTIME:
3611 temp = this->weightedEventsNoTime->at(i).pulseTime();
3612 break;
3613 }
3614 if (temp > tMax)
3615 tMax = temp;
3616 if (temp < tMin)
3617 tMin = temp;
3618 }
3619}
3620
3621DateAndTime EventList::getTimeAtSampleMax(const double &tofFactor, const double &tofOffset) const {
3622 // set up as the minimum available date time.
3623 DateAndTime tMax = DateAndTime::minimum();
3624
3625 // no events is a soft error
3626 if (this->empty())
3627 return tMax;
3628
3629 // when events are ordered by time at sample just need the first value
3630 if (this->order == TIMEATSAMPLE_SORT) {
3631 switch (eventType) {
3632 case TOF:
3633 return calculateCorrectedFullTime(this->events->back(), tofFactor, tofOffset);
3634 case WEIGHTED:
3635 return calculateCorrectedFullTime(this->weightedEvents->back(), tofFactor, tofOffset);
3636 case WEIGHTED_NOTIME:
3637 return calculateCorrectedFullTime(this->weightedEventsNoTime->back(), tofFactor, tofOffset);
3638 }
3639 }
3640
3641 // now we are stuck with a linear search
3642 size_t numEvents = this->getNumberEvents();
3643 DateAndTime temp = tMax; // start with the smallest possible value
3644 for (size_t i = 0; i < numEvents; i++) {
3645 switch (eventType) {
3646 case TOF:
3647 temp = calculateCorrectedFullTime(this->events->at(i), tofFactor, tofOffset);
3648 break;
3649 case WEIGHTED:
3650 temp = calculateCorrectedFullTime(this->weightedEvents->at(i), tofFactor, tofOffset);
3651 break;
3652 case WEIGHTED_NOTIME:
3653 temp = calculateCorrectedFullTime(this->weightedEventsNoTime->at(i), tofFactor, tofOffset);
3654 break;
3655 }
3656 if (temp > tMax)
3657 tMax = temp;
3658 }
3659 return tMax;
3660}
3661
3662DateAndTime EventList::getTimeAtSampleMin(const double &tofFactor, const double &tofOffset) const {
3663 // set up as the minimum available date time.
3664 DateAndTime tMin = DateAndTime::maximum();
3665
3666 // no events is a soft error
3667 if (this->empty())
3668 return tMin;
3669
3670 // when events are ordered by time at sample just need the first value
3671 if (this->order == TIMEATSAMPLE_SORT) {
3672 switch (eventType) {
3673 case TOF:
3674 return calculateCorrectedFullTime(this->events->front(), tofFactor, tofOffset);
3675 case WEIGHTED:
3676 return calculateCorrectedFullTime(this->weightedEvents->front(), tofFactor, tofOffset);
3677 case WEIGHTED_NOTIME:
3678 return calculateCorrectedFullTime(this->weightedEventsNoTime->front(), tofFactor, tofOffset);
3679 }
3680 }
3681
3682 // now we are stuck with a linear search
3683 size_t numEvents = this->getNumberEvents();
3684 DateAndTime temp = tMin; // start with the smallest possible value
3685 for (size_t i = 0; i < numEvents; i++) {
3686 switch (eventType) {
3687 case TOF:
3688 temp = calculateCorrectedFullTime(this->events->at(i), tofFactor, tofOffset);
3689 break;
3690 case WEIGHTED:
3691 temp = calculateCorrectedFullTime(this->weightedEvents->at(i), tofFactor, tofOffset);
3692 break;
3693 case WEIGHTED_NOTIME:
3694 temp = calculateCorrectedFullTime(this->weightedEventsNoTime->at(i), tofFactor, tofOffset);
3695 break;
3696 }
3697 if (temp < tMin)
3698 tMin = temp;
3699 }
3700 return tMin;
3701}
3702
3703// --------------------------------------------------------------------------
3709template <class T> void EventList::setTofsHelper(std::vector<T> &events, const std::vector<double> &tofs) {
3710 if (tofs.empty())
3711 return;
3712
3713 size_t x_size = tofs.size();
3714 if (events.size() != x_size)
3715 return; // should this throw an exception?
3716
3717 for (size_t i = 0; i < x_size; ++i)
3718 events[i].m_tof = tofs[i];
3719}
3720
3721// --------------------------------------------------------------------------
3728 this->order = UNSORTED;
3729
3730 // Convert the list
3731 switch (eventType) {
3732 case TOF:
3733 this->setTofsHelper(*this->events, tofs);
3734 break;
3735 case WEIGHTED:
3736 this->setTofsHelper(*this->weightedEvents, tofs);
3737 break;
3738 case WEIGHTED_NOTIME:
3739 this->setTofsHelper(*this->weightedEventsNoTime, tofs);
3740 break;
3741 }
3742}
3743
3744// ==============================================================================================
3745// ----------- MULTIPLY AND DIVIDE ---------------------------------------
3746// ==============================================================================================
3747
3748//------------------------------------------------------------------------------------------------
3756template <class T> void EventList::multiplyHelper(std::vector<T> &events, const double value, const double error) {
3757 // Square of the value
3758 const double valueSquared = value * value;
3759
3760 auto itev_end = events.end();
3761
3762 if (error == 0) {
3763 // Error-less calculation
3764 for (auto itev = events.begin(); itev != itev_end; itev++) {
3765 itev->m_errorSquared = static_cast<float>(itev->m_errorSquared * valueSquared);
3766 itev->m_weight *= static_cast<float>(value);
3767 }
3768 } else {
3769 // Carry the scalar error
3770 const double errorSquared = error * error; // Square of the value's error
3771 for (auto itev = events.begin(); itev != itev_end; itev++) {
3772 itev->m_errorSquared =
3773 static_cast<float>(itev->m_errorSquared * valueSquared + errorSquared * itev->m_weight * itev->m_weight);
3774 itev->m_weight *= static_cast<float>(value);
3775 }
3776 }
3777}
3778
3779//------------------------------------------------------------------------------------------------
3792 this->multiply(value);
3793 return *this;
3794}
3795
3796//------------------------------------------------------------------------------------------------
3825void EventList::multiply(const double value, const double error) {
3826 // Do nothing if multiplying by exactly one and there is no error
3827 if ((value == 1.0) && (error == 0.0))
3828 return;
3829
3830 switch (eventType) {
3831 case TOF:
3832 // Switch to weights if needed.
3833 this->switchTo(WEIGHTED);
3834 // Fall through
3835
3836 case WEIGHTED:
3837 multiplyHelper(*this->weightedEvents, value, error);
3838 break;
3839
3840 case WEIGHTED_NOTIME:
3842 break;
3843 }
3844}
3845
3846//------------------------------------------------------------------------------------------------
3855template <class T>
3856void EventList::multiplyHistogramHelper(std::vector<T> &events, std::span<double const> X, std::span<double const> Y,
3857 std::span<double const> E) {
3858 // Validate inputs
3859 if ((X.size() < 2) || (Y.size() != E.size()) || (X.size() != 1 + Y.size())) {
3860 std::stringstream msg;
3861 msg << "EventList::multiply() was given invalid size or "
3862 "inconsistent histogram arrays: X["
3863 << X.size() << "] "
3864 << "Y[" << Y.size() << " E[" << E.size() << "]";
3865 throw std::invalid_argument(msg.str());
3866 }
3867
3868 size_t x_size = X.size();
3869
3870 // Iterate through all events (sorted by tof)
3871 auto itev = findFirstEvent(events, T(X[0]));
3872 auto itev_end = events.end();
3873 // The above can still take you to end() if no events above X[0], so check
3874 // again.
3875 if (itev == itev_end)
3876 return;
3877
3878 // Find the first bin
3879 size_t bin = 0;
3880
3881 // Multiplier values
3882 double value;
3883 double error;
3884 double valueSquared;
3885 double errorSquared;
3886
3887 // If the tof is greater the first bin boundary, so we need to find the first
3888 // bin
3889 double tof = itev->tof();
3890 while (bin < x_size - 1) {
3891 // Within range?
3892 if ((tof >= X[bin]) && (tof < X[bin + 1]))
3893 break; // Stop increasing bin
3894 ++bin;
3895 }
3896
3897 // New bin! Find what you are multiplying!
3898 value = Y[bin];
3899 error = E[bin];
3900 valueSquared = value * value;
3901 errorSquared = error * error;
3902
3903 // Keep going through all the events
3904 while ((itev != itev_end) && (bin < x_size - 1)) {
3905 tof = itev->tof();
3906 while (bin < x_size - 1) {
3907 // Event is Within range?
3908 if ((tof >= X[bin]) && (tof < X[bin + 1])) {
3909 // Process this event. Multiply and calculate error.
3910 itev->m_errorSquared =
3911 static_cast<float>(itev->m_errorSquared * valueSquared + errorSquared * itev->m_weight * itev->m_weight);
3912 itev->m_weight *= static_cast<float>(value);
3913 break; // out of the bin-searching-while-loop
3914 }
3915 ++bin;
3916 if (bin >= x_size - 1)
3917 break;
3918
3919 // New bin! Find what you are multiplying!
3920 value = Y[bin];
3921 error = E[bin];
3922 valueSquared = value * value;
3923 errorSquared = error * error;
3924 }
3925 ++itev;
3926 }
3927}
3928
3929//------------------------------------------------------------------------------------------------
3950void EventList::multiply(std::span<double const> X, std::span<double const> Y, std::span<double const> E) {
3951 switch (eventType) {
3952 case TOF:
3953 // Switch to weights if needed.
3954 this->switchTo(WEIGHTED);
3955 // Fall through
3956
3957 case WEIGHTED:
3958 // Sorting by tof is necessary for the algorithm
3959 this->sortTof();
3961 break;
3962
3963 case WEIGHTED_NOTIME:
3964 // Sorting by tof is necessary for the algorithm
3965 this->sortTof();
3967 break;
3968 }
3969}
3970
3971//------------------------------------------------------------------------------------------------
3980template <class T>
3981void EventList::divideHistogramHelper(std::vector<T> &events, std::span<double const> X, std::span<double const> Y,
3982 std::span<double const> E) {
3983 // Validate inputs
3984 if ((X.size() < 2) || (Y.size() != E.size()) || (X.size() != 1 + Y.size())) {
3985 std::stringstream msg;
3986 msg << "EventList::divide() was given invalid size or "
3987 "inconsistent histogram arrays: X["
3988 << X.size() << "] "
3989 << "Y[" << Y.size() << " E[" << E.size() << "]";
3990 throw std::invalid_argument(msg.str());
3991 }
3992
3993 size_t x_size = X.size();
3994
3995 // Iterate through all events (sorted by tof)
3996 auto itev = findFirstEvent(events, T(X[0]));
3997 auto itev_end = events.end();
3998 // The above can still take you to end() if no events above X[0], so check
3999 // again.
4000 if (itev == itev_end)
4001 return;
4002
4003 // Find the first bin
4004 size_t bin = 0;
4005
4006 // Multiplier values
4007 double value;
4008 double error;
4009 double valError_over_value_squared;
4010
4011 // If the tof is greater the first bin boundary, so we need to find the first
4012 // bin
4013 double tof = itev->tof();
4014 while (bin < x_size - 1) {
4015 // Within range?
4016 if ((tof >= X[bin]) && (tof < X[bin + 1]))
4017 break; // Stop increasing bin
4018 ++bin;
4019 }
4020
4021 // New bin! Find what you are multiplying!
4022 value = Y[bin];
4023 error = E[bin];
4024
4025 // --- Division case ---
4026 if (value == 0) {
4027 value = std::numeric_limits<float>::quiet_NaN(); // Avoid divide by zero
4028 valError_over_value_squared = 0;
4029 } else
4030 valError_over_value_squared = error * error / (value * value);
4031
4032 // Keep going through all the events
4033 while ((itev != events.end()) && (bin < x_size - 1)) {
4034 tof = itev->tof();
4035 while (bin < x_size - 1) {
4036 // Event is Within range?
4037 if ((tof >= X[bin]) && (tof < X[bin + 1])) {
4038 // Process this event. Divide and calculate error.
4039 double newWeight = itev->m_weight / value;
4040 itev->m_errorSquared = static_cast<float>(
4041 newWeight * newWeight *
4042 ((itev->m_errorSquared / (itev->m_weight * itev->m_weight)) + valError_over_value_squared));
4043 itev->m_weight = static_cast<float>(newWeight);
4044 break; // out of the bin-searching-while-loop
4045 }
4046 ++bin;
4047 if (bin >= x_size - 1)
4048 break;
4049
4050 // New bin! Find what you are multiplying!
4051 value = Y[bin];
4052 error = E[bin];
4053
4054 // --- Division case ---
4055 if (value == 0) {
4056 value = std::numeric_limits<float>::quiet_NaN(); // Avoid divide by zero
4057 valError_over_value_squared = 0;
4058 } else
4059 valError_over_value_squared = error * error / (value * value);
4060 }
4061 ++itev;
4062 }
4063}
4064
4065//------------------------------------------------------------------------------------------------
4087void EventList::divide(std::span<double const> X, std::span<double const> Y, std::span<double const> E) {
4088 switch (eventType) {
4089 case TOF:
4090 // Switch to weights if needed.
4091 this->switchTo(WEIGHTED);
4092 // Fall through
4093
4094 case WEIGHTED:
4095 // Sorting by tof is necessary for the algorithm
4096 this->sortTof();
4097 divideHistogramHelper(*this->weightedEvents, X, Y, E);
4098 break;
4099
4100 case WEIGHTED_NOTIME:
4101 // Sorting by tof is necessary for the algorithm
4102 this->sortTof();
4104 break;
4105 }
4106}
4107
4108//------------------------------------------------------------------------------------------------
4118 if (value == 0.0)
4119 throw std::invalid_argument("EventList::divide() called with value of 0.0. Cannot divide by zero.");
4120 this->multiply(1.0 / value, 0.0);
4121 return *this;
4122}
4123
4124//------------------------------------------------------------------------------------------------
4134void EventList::divide(const double value, const double error) {
4135 if (value == 0.0)
4136 throw std::invalid_argument("EventList::divide() called with value of 0.0. Cannot divide by zero.");
4137 // Do nothing if dividing by exactly 1.0, no error
4138 else if (value == 1.0 && error == 0.0)
4139 return;
4140
4141 // We'll multiply by 1/value
4142 double invValue = 1.0 / value;
4143 // Relative error remains the same
4144 double invError = (error / value) * invValue;
4145
4146 this->multiply(invValue, invError);
4147}
4148
4149// ==============================================================================================
4150// ----------- SPLITTING AND FILTERING ---------------------------------------
4151// ==============================================================================================
4152//------------------------------------------------------------------------------------------------
4162void EventList::filterByPulseTime(Types::Core::DateAndTime start, Types::Core::DateAndTime stop,
4163 EventList &output) const {
4164 if (this == &output) {
4165 throw std::invalid_argument("In-place filtering is not allowed");
4166 }
4167
4168 // Start by sorting the event list by pulse time.
4169 this->sortPulseTime();
4170 // Clear the output
4171 output.clear();
4172 // Has to match the given type
4173 output.switchTo(eventType);
4174 output.setDetectorIDs(this->getDetectorIDs());
4175 output.setHistogram(m_histogram);
4176 output.setSortOrder(this->order);
4177
4178 // Iterate through all events (sorted by pulse time)
4179 switch (eventType) {
4180 case TOF:
4181 filterByPulseTimeHelper(*this->events, start, stop, *output.events);
4182 break;
4183 case WEIGHTED:
4184 filterByPulseTimeHelper(*this->weightedEvents, start, stop, *output.weightedEvents);
4185 break;
4186 case WEIGHTED_NOTIME:
4187 throw std::runtime_error("EventList::filterByPulseTime() called on an "
4188 "EventList that no longer has time information.");
4189 break;
4190 }
4191}
4192
4205void EventList::filterByPulseTime(Kernel::TimeROI const *timeRoi, EventList *output) const {
4206
4207 this->sortPulseTime();
4208 // Clear the output
4209
4210 output->clear();
4211 output->setDetectorIDs(this->getDetectorIDs());
4212 output->setHistogram(m_histogram);
4213 // Has to match the given type
4214 output->switchTo(eventType);
4215
4216 if ((timeRoi == nullptr) || (timeRoi->useAll())) {
4217 throw std::invalid_argument("TimeROI can not use all time");
4218 }
4219 const auto &intervals = timeRoi->toTimeIntervals();
4220 if (intervals.empty())
4221 return;
4222
4223 switch (eventType) {
4224 case TOF:
4225 filterByTimeROIHelper(*this->events, intervals, output);
4226 break;
4227 case WEIGHTED:
4228 filterByTimeROIHelper(*this->weightedEvents, intervals, output);
4229 break;
4230 case WEIGHTED_NOTIME:
4231 throw std::runtime_error("EventList::filterByPulseTime() called on an "
4232 "EventList that no longer has time information.");
4233 break;
4234 }
4235}
4236
4244template <class T>
4245void EventList::filterByPulseTimeHelper(std::vector<T> &events, DateAndTime start, DateAndTime stop,
4246 std::vector<T> &output) {
4247 std::copy_if(events.begin(), events.end(), std::back_inserter(output),
4248 [start, stop](const T &t) { return (t.m_pulsetime >= start) && (t.m_pulsetime < stop); });
4249}
4250
4257template <class T>
4258void EventList::filterByTimeROIHelper(std::vector<T> &events, const std::vector<Kernel::TimeInterval> &intervals,
4259 EventList *output) {
4260 // Iterate through the splitter at the same time
4261 auto itspl = intervals.cbegin();
4262 auto itspl_end = intervals.cend();
4263 // Iterate through all events (sorted by tof)
4264 auto itev = events.cbegin();
4265 auto itev_end = events.cend();
4266
4267 // This is the time of the first section. Anything before is thrown out.
4268 while (itspl != itspl_end) {
4269 // Get the splitting interval times and destination
4270 DateAndTime start = itspl->start();
4271 DateAndTime stop = itspl->stop();
4272 // Skip the events before the start of the time
4273 while ((itev != itev_end) && (itev->m_pulsetime < start))
4274 itev++;
4275
4276 // Go through all the events that are in the interval (if any)
4277 while ((itev != itev_end) && (itev->m_pulsetime < stop)) {
4278 // Copy the event into another
4279 const T eventCopy(*itev);
4280 output->addEventQuickly(eventCopy);
4281 ++itev;
4282 }
4283
4284 // Go to the next interval
4285 ++itspl;
4286 // But if we reached the end, then we are done.
4287 if (itspl == itspl_end)
4288 break;
4289
4290 // No need to keep looping through the filter if we are out of events
4291 if (itev == itev_end)
4292 break;
4293 }
4294 // Done!
4295}
4296
4302 if (timeRoi == nullptr) {
4303 throw std::runtime_error("TimeROI can not be a nullptr\n");
4304 }
4305 if (timeRoi->useAll()) {
4306 throw std::invalid_argument("TimeROI can not be empty\n");
4307 }
4308 // Start by sorting the event list by pulse time.
4309 this->sortPulseTime();
4310
4311 // Iterate through all events (sorted by pulse time)
4312 switch (eventType) {
4313 case TOF:
4314 filterInPlaceHelper(timeRoi, *this->events);
4315 break;
4316 case WEIGHTED:
4317 filterInPlaceHelper(timeRoi, *this->weightedEvents);
4318 break;
4319 case WEIGHTED_NOTIME:
4320 throw std::runtime_error("EventList::filterInPlace() called on an "
4321 "EventList that no longer has time information.");
4322 break;
4323 }
4324}
4325
4352template <class T>
4353void EventList::filterInPlaceHelper(Kernel::TimeROI const *timeRoi, typename std::vector<T> &events) {
4354
4355 const auto splitter = timeRoi->toTimeIntervals();
4356 // Iterate through the splitter at the same time
4357 auto itspl = splitter.cbegin();
4358 auto itspl_end = splitter.cend();
4359 DateAndTime start, stop;
4360
4361 // Iterate for the input
4362 auto itev = events.begin();
4363 auto itev_end = events.end();
4364
4365 // Iterator for the outputted list; will follow the input except when events
4366 // are dropped.
4367 auto itOut = events.begin();
4368
4369 // This is the time of the first section. Anything before is thrown out.
4370 while (itspl != itspl_end) {
4371 // Get the splitting interval times and destination
4372 start = itspl->start();
4373 stop = itspl->stop();
4374 // Skip the events before the start of the time
4375 while ((itev != itev_end) && (itev->m_pulsetime < start))
4376 itev++;
4377
4378 // Are we aligned in the input vs output?
4379 bool copyingInPlace = (itOut == itev);
4380 if (copyingInPlace) {
4381 while ((itev != itev_end) && (itev->m_pulsetime < stop))
4382 ++itev;
4383 // Make sure the iterators still match
4384 itOut = itev;
4385 } else {
4386 // Go through all the events that are in the interval (if any)
4387 while ((itev != itev_end) && (itev->m_pulsetime < stop)) {
4388 *itOut = *itev;
4389 ++itOut;
4390 ++itev;
4391 }
4392 }
4393
4394 // Go to the next interval
4395 ++itspl;
4396 // But if we reached the end, then we are done.
4397 if (itspl == itspl_end)
4398 break;
4399
4400 // No need to keep looping through the filter if we are out of events
4401 if (itev == itev_end)
4402 break;
4403
4404 } // Looping through entries in the splitter vector
4405
4406 // Ok, now resize the event list to reflect the fact that it (probably) shrank
4407 events.resize(std::size_t(std::distance(events.begin(), itOut)));
4408}
4409
4414void EventList::initializePartials(std::map<int, EventList *> partials) const {
4415
4416 // collect the state from events which is to be transferred to the partials
4417 bool removeDetIDs{true};
4418 const auto histogramLocal = this->getHistogram();
4419 const auto eventTypeLocal = this->getEventType();
4420
4421 // lambda expression initializing one partial
4422 auto initPartial = [&](EventList *partial) {
4423 partial->clear(removeDetIDs);
4424 partial->copyInfoFrom(*this);
4425 partial->setHistogram(histogramLocal);
4426 partial->switchTo(eventTypeLocal);
4427 };
4428
4429 // iterate over the partials
4430 std::for_each(partials.cbegin(), partials.cend(),
4431 [&](const std::pair<int, EventList *> &pair) { initPartial(pair.second); });
4432}
4433
4443void getEventsFrom(EventList &el, std::vector<TofEvent> *&events) { events = &el.getEvents(); }
4444void getEventsFrom(const EventList &el, std::vector<TofEvent> const *&events) { events = &el.getEvents(); }
4445
4455void getEventsFrom(EventList &el, std::vector<WeightedEvent> *&events) { events = &el.getWeightedEvents(); }
4456void getEventsFrom(const EventList &el, std::vector<WeightedEvent> const *&events) { events = &el.getWeightedEvents(); }
4457
4467void getEventsFrom(EventList &el, std::vector<WeightedEventNoTime> *&events) { events = &el.getWeightedEventsNoTime(); }
4468void getEventsFrom(const EventList &el, std::vector<WeightedEventNoTime> const *&events) {
4469 events = &el.getWeightedEventsNoTime();
4470}
4471
4479template <class T>
4480void EventList::convertUnitsViaTofHelper(typename std::vector<T> &events, Mantid::Kernel::Unit const *fromUnit,
4481 Mantid::Kernel::Unit const *toUnit) {
4482 for (auto &itev : events) {
4483 // Conver to TOF
4484 const double tof = fromUnit->singleToTOF(itev.m_tof);
4485 // And back from TOF to whatever
4486 itev.m_tof = toUnit->singleFromTOF(tof);
4487 }
4488}
4489
4490//--------------------------------------------------------------------------
4499 // Check for initialized
4500 if (!fromUnit || !toUnit)
4501 throw std::runtime_error("EventList::convertUnitsViaTof(): one of the units is NULL!");
4502 if (!fromUnit->isInitialized())
4503 throw std::runtime_error("EventList::convertUnitsViaTof(): fromUnit is not initialized!");
4504 if (!toUnit->isInitialized())
4505 throw std::runtime_error("EventList::convertUnitsViaTof(): toUnit is not initialized!");
4506
4507 switch (eventType) {
4508 case TOF:
4509 convertUnitsViaTofHelper(*this->events, fromUnit, toUnit);
4510 break;
4511 case WEIGHTED:
4512 convertUnitsViaTofHelper(*this->weightedEvents, fromUnit, toUnit);
4513 break;
4514 case WEIGHTED_NOTIME:
4515 convertUnitsViaTofHelper(*this->weightedEventsNoTime, fromUnit, toUnit);
4516 break;
4517 }
4518}
4519
4520//--------------------------------------------------------------------------
4527template <class T>
4528void EventList::convertUnitsQuicklyHelper(typename std::vector<T> &events, const double &factor, const double &power) {
4529 for (auto &event : events) {
4530 // Output unit = factor * (input) ^ power
4531 event.m_tof = factor * std::pow(event.m_tof, power);
4532 }
4533}
4534
4535//--------------------------------------------------------------------------
4541void EventList::convertUnitsQuickly(const double &factor, const double &power) {
4542 switch (eventType) {
4543 case TOF:
4544 convertUnitsQuicklyHelper(*this->events, factor, power);
4545 break;
4546 case WEIGHTED:
4547 convertUnitsQuicklyHelper(*this->weightedEvents, factor, power);
4548 break;
4549 case WEIGHTED_NOTIME:
4550 convertUnitsQuicklyHelper(*this->weightedEventsNoTime, factor, power);
4551 break;
4552 }
4553}
4554
4555HistogramData::Histogram &EventList::mutableHistogramRef() {
4556 if (mru)
4557 mru->deleteIndex(this);
4558 return m_histogram;
4559}
4560
4561void EventList::checkAndSanitizeHistogram(HistogramData::Histogram &histogram) {
4562 if (histogram.xMode() != HistogramData::Histogram::XMode::BinEdges)
4563 throw std::runtime_error("EventList: setting histogram with storage mode "
4564 "other than BinEdges is not possible");
4565 if (histogram.sharedY() || histogram.sharedE())
4566 throw std::runtime_error("EventList: setting histogram data with non-null "
4567 "Y or E data is not possible");
4568 // Avoid flushing of YMode: we only change X but YMode depends on events.
4569 if (histogram.yMode() == HistogramData::Histogram::YMode::Uninitialized)
4570 histogram.setYMode(m_histogram.yMode());
4571 if (histogram.yMode() != m_histogram.yMode())
4572 throw std::runtime_error("EventList: setting histogram data with different "
4573 "YMode is not possible");
4574}
4575
4577 throw std::runtime_error("EventList: setting Points as X data is not "
4578 "possible, only BinEdges are supported");
4579}
4580
4582 throw std::runtime_error("EventList: Cannot set Y or E data, these data are "
4583 "generated automatically based on the events");
4584}
4585
4586} // namespace Mantid::DataObjects
gsl_vector * tmp
std::vector< float > m_tof
sum of all time-of-flight within the bin
const std::vector< double > & rhs
const double m_tofFactor
Definition EventList.cpp:70
const double m_tofShift
Definition EventList.cpp:71
double value
The value of the point.
Definition FitMW.cpp:51
double error
double left
double right
#define fabs(x)
Definition Matrix.cpp:22
int count
counter
Definition Matrix.cpp:37
#define PARALLEL_THREAD_NUMBER
double tolerance
#define MSVC_DIAG_ON(x)
#define MSVC_DIAG_OFF(x)
#define GNU_DIAG_ON(x)
#define GNU_DIAG_OFF(x)
This is a collection of macros for turning compiler warnings off in a controlled manner.
IEventList : Interface to Mantid::DataObjects::EventList class, used to expose to PythonAPI.
Definition IEventList.h:27
A "spectrum" is an object that holds the data for a particular spectrum, in particular:
Definition ISpectrum.h:38
void setHistogram(T &&...data)
Sets the Histogram associated with this spectrum.
Definition ISpectrum.h:110
void addDetectorIDs(const std::set< detid_t > &detIDs)
Add a set of detector IDs to the set of detector IDs.
Definition ISpectrum.cpp:62
const HistogramData::HistogramX & x() const
Definition ISpectrum.h:185
virtual const HistogramData::HistogramE & e() const
Definition ISpectrum.h:187
void setDetectorIDs(const std::set< detid_t > &detIDs)
Set the detector IDs to be the set given.
Definition ISpectrum.cpp:94
void clearDetectorIDs()
Clear the detector IDs set.
void copyInfoFrom(const ISpectrum &other)
Copy spectrum number and detector IDs, but not X vector, from another ISpectrum.
Definition ISpectrum.cpp:24
const HistogramData::HistogramDx & dx() const
Definition ISpectrum.h:188
Kernel::cow_ptr< HistogramData::HistogramX > sharedX() const
Definition ISpectrum.h:199
virtual const HistogramData::HistogramY & y() const
Definition ISpectrum.h:186
const std::set< detid_t > & getDetectorIDs() const
Get a const reference to the detector IDs set.
virtual void copyDataInto(DataObjects::EventList &) const
Override in child classes for polymorphic copying of data.
void setSharedX(const Kernel::cow_ptr< HistogramData::HistogramX > &x) &
Definition ISpectrum.h:203
HistogramData::HistogramX & mutableX() &
Definition ISpectrum.h:189
A class for holding :
Definition EventList.h:58
void sortTimeAtSample(const double &tofFactor, const double &tofShift, bool forceResort=false) const
Sort events by time at sample.
size_t getMemorySize() const override
Memory used by this event list.
static void getWeightErrorsHelper(const std::vector< T > &events, std::vector< double > &weightErrors)
Get the weight error member of all events in a list.
void addPulsetimes(const std::vector< double > &seconds) override
Add an offset to the pulsetime (wall-clock time) of each event in the list.
Mantid::Types::Core::DateAndTime getTimeAtSampleMax(const double &tofFactor, const double &tofOffset) const override
Get the maximum time at sample.
std::vector< Types::Core::DateAndTime > getPulseTimes() const override
Get the pulse times of each event in this EventList.
static std::optional< size_t > findLinearBin(std::span< double const > X, const double tof, const double divisor, const double offset, const bool findExact=true)
Find the bin which this TOF value falls in with linear binning, assumes TOF is in range of X.
void setTofs(const MantidVec &tofs) override
Set a list of TOFs to the current event list.
void maskTof(const double tofMin, const double tofMax) override
Mask out events that have a tof between tofMin and tofMax (inclusively).
EventList(const Mantid::API::EventType event_type=Mantid::API::EventType::TOF)
Constructor (empty)
void multiply(const double value, const double error=0.0) override
Multiply the weights in this event list by a scalar variable with an error; though the error can be 0...
void checkWorksWithPoints() const override
HistogramData::CountStandardDeviations countStandardDeviations() const override
void switchToWeightedEvents()
Switch the EventList to use WeightedEvents instead of TofEvent.
static size_t findExactBin(std::span< double const > X, const double tof, const size_t n_bin)
Find the exact bin which a TOF falls in starting from the provided estimated one.
EventList & operator=(const EventList &)
Copy into this event list from another.
void compressFatEvents(const double tolerance, const Types::Core::DateAndTime &timeStart, const double seconds, EventList *destination)
HistogramData::Histogram & mutableHistogramRef() override
void convertUnitsViaTofHelper(typename std::vector< T > &events, Mantid::Kernel::Unit const *fromUnit, Mantid::Kernel::Unit const *toUnit)
Helper function for the conversion to TOF.
static void histogramForWeightsHelper(const std::vector< T > &events, std::span< double const > X, MantidVec &Y, MantidVec &E)
Generates both the Y and E (error) histograms for an EventList with WeightedEvents.
HistogramData::Counts counts() const override
std::vector< double > getWeights() const override
Return the list of event weight values.
void convertUnitsQuickly(const double &factor, const double &power)
Convert the event's TOF (x) value according to a simple output = a * (input^b) relationship.
void compressEvents(double tolerance, EventList *destination)
Compress the event list by grouping events with the same TOF (within a given tolerance).
double getTofMax() const override
static void filterByTimeROIHelper(std::vector< T > &events, const std::vector< Kernel::TimeInterval > &intervals, EventList *output)
Filter a vector of events into another based on TimeROI.
virtual size_t histogram_size() const
Return the size of the histogram data.
HistogramData::Histogram histogram() const override
Returns the Histogram associated with this spectrum.
~EventList() override
Destructor.
std::vector< Types::Event::TofEvent > & getEvents()
Return the list of TofEvents contained.
HistogramData::FrequencyStandardDeviations frequencyStandardDeviations() const override
void setSortOrder(const EventSortType order) const
Manually set the event list sort order value.
static void processWeightedEvents(const std::vector< T > &events, std::vector< WeightedEventNoTime > &out, const std::shared_ptr< std::vector< double > > histogram_bin_edges, struct FindBin findBin)
HistogramData::Histogram m_histogram
Histogram object holding the histogram data. Currently only X.
Definition EventList.h:347
void generateCountsHistogramPulseTime(const double &xMin, const double &xMax, MantidVec &Y, const double TofMin=std::numeric_limits< double >::lowest(), const double TofMax=std::numeric_limits< double >::max()) const
With respect to PulseTime fill a histogram given equal histogram bins.
double getTofMin() const override
static void compressEventsHelper(const std::vector< T > &events, std::vector< WeightedEventNoTime > &out, double tolerance)
Compress the event list by grouping events with the same TOF.
Mantid::API::EventType eventType
What type of event is in our list.
Definition EventList.h:359
EventSortType getSortType() const
Return the type of sorting used in this event list.
void switchToWeightedEventsNoTime()
Switch the EventList to use WeightedEventNoTime's instead of TofEvent.
static void setTofsHelper(std::vector< T > &events, const std::vector< double > &tofs)
Set a list of TOFs to the current event list.
void addTof(const double offset) override
Add an offset to the TOF of each event in the list.
bool equals(const EventList &rhs, const double tolTof, const double tolWeight, const int64_t tolPulse) const
Kernel::cow_ptr< HistogramData::HistogramE > sharedE() const override
static std::size_t maskConditionHelper(std::vector< T > &events, const std::vector< bool > &mask)
Mask out events by the condition vector.
EventList & operator/=(const double value)
Operator to divide the weights in this EventList by an error-less scalar.
void checkAndSanitizeHistogram(HistogramData::Histogram &histogram) override
void copyDataInto(EventList &sink) const override
Used by copyDataFrom for dynamic dispatch for its source.
const HistogramData::HistogramY & y() const override
HistogramData::Frequencies frequencies() const override
void filterInPlace(const Kernel::TimeROI *timeRoi)
Use a SplittingIntervalVec to filter the event list in place.
std::size_t getNumberEvents() const override
Return the number of events in the list.
Mantid::API::EventType getEventType() const override
Return the type of Event vector contained within.
const MantidVec & readDx() const override
Deprecated, use dx() instead.
void scaleTof(const double factor) override
Convert the units in the TofEvent's m_tof field to some other value, by scaling by a multiplier.
void getPulseTimeMinMax(Mantid::Types::Core::DateAndTime &tMin, Mantid::Types::Core::DateAndTime &tM) const
static void multiplyHelper(std::vector< T > &events, const double value, const double error=0.0)
Helper method for multiplying an event list by a scalar value with/without error.
HistogramData::FrequencyVariances frequencyVariances() const override
std::vector< WeightedEventNoTime > & getWeightedEventsNoTime()
Return the list of WeightedEvent contained.
void setMRU(EventWorkspaceMRU *newMRU)
Sets the MRU list for this event list.
EventWorkspaceMRU * mru
MRU lists of the parent EventWorkspace.
Definition EventList.h:365
void convertUnitsViaTof(Mantid::Kernel::Unit const *fromUnit, Mantid::Kernel::Unit const *toUnit)
Converts the X units in each event by going through TOF.
static void integrateHelper(std::vector< T > &events, const double minX, const double maxX, const bool entireRange, double &sum, double &error)
Integrate the events between a range of X values, or all events.
Mantid::Types::Core::DateAndTime getPulseTimeMax() const override
Kernel::cow_ptr< HistogramData::HistogramX > ptrX() const override
Deprecated, use sharedX() instead. Returns a pointer to the x data.
void generateHistogram(std::span< double const > X, MantidVec &Y, MantidVec &E, bool skipError=false) const override
Generates both the Y and E (error) histograms w.r.t TOF for an EventList with or without WeightedEven...
static void createWeightedEvents(std::vector< WeightedEventNoTime > &out, const std::vector< double > &tof, const std::vector< T > &weight, const std::vector< T > &error)
MantidVec & dataX() override
Deprecated, use mutableX() instead.
void generateCountsHistogram(std::span< double const > X, MantidVec &Y) const
Fill a histogram given specified histogram bounds.
void generateHistogramPulseTime(std::span< double const > X, MantidVec &Y, MantidVec &E, bool skipError=false) const override
Generates both the Y and E (error) histograms w.r.t Pulse Time for an EventList with or without Weigh...
void clearUnused()
Clear any unused event lists (the ones that do not match the currently used type).
MantidVec & dataDx() override
Deprecated, use mutableDx() instead.
bool operator==(const EventList &rhs) const
Equality operator between EventList's.
void generateCountsHistogramTimeAtSample(std::span< double const > X, MantidVec &Y, const double &tofFactor, const double &tofOffset) const
With respect to Time at Sample, fill a histogram given specified histogram bounds.
EventSortType order
Last sorting order.
Definition EventList.h:362
static void divideHistogramHelper(std::vector< T > &events, std::span< double const > X, std::span< double const > Y, std::span< double const > E)
Helper method for dividing an event list by a histogram with error.
void convertTofHelper(std::vector< T > &events, const std::function< double(double)> &func)
void generateHistogramTimeAtSample(std::span< double const > X, MantidVec &Y, MantidVec &E, const double &tofFactor, const double &tofOffset, bool skipError=false) const override
Generates both the Y and E (error) histograms w.r.t Time at sample position.
void clearData() override
Mask the spectrum to this value. Removes all events.
EventList & operator+=(const Types::Event::TofEvent &event)
Append an event to the histogram.
static void minusHelper(std::vector< T1 > &events, const std::vector< T2 > &more_events)
SUBTRACT another EventList from this event list.
void addPulsetimesHelper(std::vector< T > &events, const std::vector< double > &seconds)
Add an offset per event to the pulsetime (wall-clock time) of each event in the list.
static std::optional< size_t > findLogBin(std::span< double const > X, const double tof, const double divisor, const double offset, const bool findExact=true)
Find the bin which this TOF value falls in with log binning, assumes TOF is in range of X.
void addPulsetime(const double seconds) override
Add an offset to the pulsetime (wall-clock time) of each event in the list.
HistogramData::Histogram getHistogram() const
Returns a copy of the Histogram associated with this spectrum.
void reserve(size_t num) override
Reserve a certain number of entries in event list of the specified eventType.
void setX(const Kernel::cow_ptr< HistogramData::HistogramX > &X) override
Deprecated, use setSharedX() instead.
std::vector< double > getWeightErrors() const override
Return the list of event weight error values.
void sortTof() const
Sort events by TOF in one thread.
HistogramData::CountVariances countVariances() const override
MantidVec & dataE() override
Deprecated, use mutableE() instead.
Definition EventList.h:200
static void compressFatEventsHelper(const std::vector< T > &events, std::vector< WeightedEvent > &out, const double tolerance, const Mantid::Types::Core::DateAndTime &timeStart, const double seconds)
MantidVec & dataY() override
Deprecated, use mutableY() instead.
Definition EventList.h:194
Mantid::Types::Core::DateAndTime getTimeAtSampleMin(const double &tofFactor, const double &tofOffset) const override
Get the minimum time at sample.
const HistogramData::HistogramE & e() const override
std::vector< Types::Core::DateAndTime > getPulseTOFTimes() const
Get the Pulse-time + TOF for each event in this EventList.
bool isSortedByTof() const override
Return true if the event list is sorted by TOF.
bool operator!=(const EventList &rhs) const
Inequality comparator.
std::mutex m_sortMutex
Mutex that is locked while sorting an event list.
Definition EventList.h:368
static void filterByPulseTimeHelper(std::vector< T > &events, Types::Core::DateAndTime start, Types::Core::DateAndTime stop, std::vector< T > &output)
std::unique_ptr< std::vector< WeightedEvent > > weightedEvents
List of WeightedEvent's.
Definition EventList.h:353
void initializePartials(std::map< int, EventList * > partials) const
Initialize the detector ID's and event type of the destination event lists when splitting this list.
void reverse()
Reverse the histogram boundaries and the associated events if they are sorted by time-of-flight.
std::vector< Types::Core::DateAndTime > eventTimesCalculator(const UnaryOperation &timesCalc) const
Compute a time (for instance, pulse-time plus TOF) associated to each event in the list.
void switchTo(Mantid::API::EventType newType) override
Switch the EventList to use the given EventType (TOF, WEIGHTED, or WEIGHTED_NOTIME)
void filterInPlaceHelper(Kernel::TimeROI const *timeRoi, typename std::vector< T > &events)
Perform an in-place filtering on a vector of either TofEvent's or WeightedEvent's.
bool empty() const
Much like stl containers, returns true if there is nothing in the event list.
const MantidVec & readX() const override
Deprecated, use x() instead. Returns the x data const.
WeightedEvent getEvent(size_t event_number)
Return the given event in the list.
void integrate(const double minX, const double maxX, const bool entireRange, double &sum, double &error) const
Integrate the events between a range of X values, or all events.
Kernel::cow_ptr< HistogramData::HistogramY > sharedY() const override
void createFromHistogram(const ISpectrum *inSpec, bool GenerateZeros, bool GenerateMultipleEvents, int MaxEventsPerBin)
Create an EventList from a histogram.
void convertUnitsQuicklyHelper(typename std::vector< T > &events, const double &factor, const double &power)
Convert the event's TOF (x) value according to a simple output = a * (input^b) relationship.
MantidVec * makeDataE() const
Calculates and returns a pointer to the E histogrammed data.
void checkIsYAndEWritable() const override
void copyDataFrom(const ISpectrum &source) override
Copy data from another EventList, via ISpectrum reference.
void filterByPulseTime(Types::Core::DateAndTime start, Types::Core::DateAndTime stop, EventList &output) const
Filter this EventList into an output EventList, using keeping only events within the >= start and < e...
void addPulsetimeHelper(std::vector< T > &events, const double seconds)
Add an offset to the pulsetime (wall-clock time) of each event in the list.
std::vector< double > getTofs() const override
Get the times-of-flight of each event in this EventList.
static void multiplyHistogramHelper(std::vector< T > &events, std::span< double const > X, std::span< double const > Y, std::span< double const > E)
Helper method for multiplying an event list by a histogram with error.
EventList & operator*=(const double value)
Operator to multiply the weights in this EventList by an error-less scalar.
void convertTof(std::function< double(double)> func, const int sorting=0) override
Mantid::Types::Core::DateAndTime getPulseTimeMin() const override
std::unique_ptr< std::vector< WeightedEventNoTime > > weightedEventsNoTime
List of WeightedEvent's.
Definition EventList.h:356
static void getWeightsHelper(const std::vector< T > &events, std::vector< double > &weights)
Get the weight member of all events in a list.
void maskCondition(const std::vector< bool > &mask) override
Mask out events by the condition vector.
void sortPulseTime() const
Sort events by Frame.
std::vector< T >::const_iterator findFirstTimeAtSampleEvent(const std::vector< T > &events, const double seek_time, const double &tofFactor, const double &tofOffset) const
Utility function: Returns the iterator into events of the first TofEvent with time at sample > seek_t...
void sortPulseTimeTOFDelta(const Types::Core::DateAndTime &start, const double seconds) const
Sort by the pulse time with a tolerance.
static std::size_t maskTofHelper(std::vector< T > &events, const double tofMin, const double tofMax)
Mask out events that have a tof between tofMin and tofMax (inclusively).
void addEventQuickly(const Types::Event::TofEvent &event)
Append an event to the histogram, without clearing the cache, to make it faster.
Definition EventList.h:106
static std::vector< T >::const_iterator findFirstPulseEvent(const std::vector< T > &events, const double seek_pulsetime)
Utility function: Returns the iterator into events of the first TofEvent with pulsetime() > seek_puls...
std::vector< Types::Core::DateAndTime > getPulseTOFTimesAtSample(const double &factor, const double &shift) const
Get the Pulse-time + time-of-flight of the neutron up to the sample, for each event in this EventList...
void generateErrorsHistogram(std::span< double const > Y, MantidVec &E) const
Generate the Error histogram for the provided counts histogram.
EventList & operator-=(const EventList &more_events)
SUBTRACT another EventList from this event list.
void sort(const EventSortType order) const
Sort events by TOF or Frame.
MantidVec * makeDataY() const
Calculates and returns a pointer to the Y histogrammed data.
std::unique_ptr< std::vector< Types::Event::TofEvent > > events
List of TofEvent (no weights).
Definition EventList.h:350
static void getTofsHelper(const std::vector< T > &events, std::vector< double > &tofs)
Get the m_tof member of all events in a list.
void clear(const bool removeDetIDs=true) override
Clear the list of events and any associated detector ID's.
std::vector< WeightedEvent > & getWeightedEvents()
Return the list of WeightedEvent contained.
void divide(const double value, const double error=0.0) override
Divide the weights in this event list by a scalar with an (optional) error.
This is a container for the MRU (most-recently-used) list of generated histograms.
void insertY(size_t thread_num, YType data, const EventList *index)
Insert a new histogram into the MRU.
EType findE(size_t thread_num, const EventList *index)
Find a Y histogram in the MRU.
void insertE(size_t thread_num, EType data, const EventList *index)
Insert a new histogram into the MRU.
void ensureEnoughBuffersE(size_t thread_num) const
This function makes sure that there are enough data buffers (MRU's) for E for the number of threads r...
YType findY(size_t thread_num, const EventList *index)
Find a Y histogram in the MRU.
void deleteIndex(const EventList *index)
Delete any entries in the MRU at the given index.
void ensureEnoughBuffersY(size_t thread_num) const
This function makes sure that there are enough data buffers (MRU's) for Y for the number of threads r...
1D histogram implementation.
Definition Histogram1D.h:19
Info about a single neutron detection event, including a weight and error value, but excluding the pu...
Definition Events.h:91
Info about a single neutron detection event, including a weight and error value:
Definition Events.h:39
TimeROI : Object that holds information about when the time measurement was active.
Definition TimeROI.h:18
const std::vector< Kernel::TimeInterval > toTimeIntervals() const
This method is to lend itself to helping with transition.
Definition TimeROI.cpp:557
bool useAll() const
TimeROI selects all time to be used.
Definition TimeROI.cpp:693
The base units (abstract) class.
Definition Unit.h:42
virtual double singleToTOF(const double x) const =0
Convert a single X value to TOF.
virtual double singleFromTOF(const double tof) const =0
Convert a single tof value to this unit.
bool isInitialized() const
Definition Unit.h:186
Implements a copy on write data template.
Definition cow_ptr.h:41
EventType
What kind of event list is being stored.
Definition IEventList.h:19
std::size_t numEvents(Nexus::File &file, bool &hasTotalCounts, bool &oldNeXusFileNames, const std::string &prefix)
Get the number of events in the currently opened group.
DLLExport void getEventsFrom(EventList &el, std::vector< Types::Event::TofEvent > *&events)
bool compareEventPulseTime(const TofEvent &e1, const TofEvent &e2)
Compare two events' FRAME id, return true if e1 should be before e2.
static std::vector< T >::const_iterator findFirstEvent(const std::vector< T > &events, T seek_tof)
Utility function: Returns the iterator into events of the first TofEvent with tof() > seek_tof Will r...
EventSortType
How the event list is sorted.
Definition EventList.h:33
bool compareEventPulseTimeTOF(const TofEvent &e1, const TofEvent &e2)
Compare two events' FRAME id, return true if e1 should be before e2.
MANTID_KERNEL_DLL Types::Core::DateAndTime averageSorted(const std::vector< Types::Core::DateAndTime > &times)
averageSorted Assuming that the vector is sorted, find the average time
std::vector< double > MantidVec
typedef for the data storage used in Mantid matrix workspaces
Definition cow_ptr.h:172
int32_t specnum_t
Typedef for a spectrum Number.
Definition IDTypes.h:14
std::optional< size_t > operator()(std::span< double const > X, const double tof, const bool findExact)
FindBin(double step, double xmin)
std::optional< size_t >(* findBin)(std::span< double const >, const double, const double, const double, const bool)
bool operator()(const TofEvent &e1, const TofEvent &e2)
comparePulseTimeTOFDelta(const Types::Core::DateAndTime &start, const double seconds)