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"
16
17#ifdef _MSC_VER
18// qualifier applied to function type has no meaning; ignored
19#pragma warning(disable : 4180)
20#endif
21#include "tbb/parallel_sort.h"
22#ifdef _MSC_VER
23#pragma warning(default : 4180)
24#endif
25
26#include <algorithm>
27#include <cfloat>
28#include <cmath>
29#include <functional>
30#include <limits>
31#include <stdexcept>
32
33using std::ostream;
34using std::runtime_error;
35using std::size_t;
36using std::vector;
37
38namespace Mantid::DataObjects {
39using Types::Core::DateAndTime;
40using Types::Event::TofEvent;
41using namespace Mantid::API;
42
43namespace {
44
45constexpr double SEC_TO_NANO{1.e9};
46
47// minimum event vector length to use tbb::parallel_sort
48// this is 4x what parallel_sort uses in the indidividual blocks
49constexpr size_t MIN_VEC_LENGTH_PARALLEL_SORT{2000};
50
58template <typename EventType>
59int64_t calculateCorrectedFullTime(const EventType &event, const double tofFactor, const double tofShift) {
60 return event.pulseTime().totalNanoseconds() +
61 static_cast<int64_t>(tofFactor * (event.tof() * 1.0E3) + (tofShift * 1.0E9));
62}
63
67template <typename EventType> class CompareTimeAtSample {
68private:
69 const double m_tofFactor;
70 const double m_tofShift;
71
72public:
73 CompareTimeAtSample(const double tofFactor, const double tofShift) : m_tofFactor(tofFactor), m_tofShift(tofShift) {}
74
84 bool operator()(const EventType &e1, const EventType &e2) const {
85 const auto tAtSample1 = calculateCorrectedFullTime(e1, m_tofFactor, m_tofShift);
86 const auto tAtSample2 = calculateCorrectedFullTime(e2, m_tofFactor, m_tofShift);
87 return (tAtSample1 < tAtSample2);
88 }
89};
90} // namespace
91//==========================================================================
94//==========================================================================
99bool compareEventPulseTime(const TofEvent &e1, const TofEvent &e2) { return (e1.pulseTime() < e2.pulseTime()); }
100
107bool compareEventPulseTimeTOF(const TofEvent &e1, const TofEvent &e2) {
108
109 if (e1.pulseTime() < e2.pulseTime()) {
110 return true;
111 } else if ((e1.pulseTime() == e2.pulseTime()) && (e1.tof() < e2.tof())) {
112 return true;
113 }
114
115 return false;
116}
117
118// comparator for pulse time with tolerance
120 explicit comparePulseTimeTOFDelta(const Types::Core::DateAndTime &start, const double seconds)
121 : startNano(start.totalNanoseconds()), deltaNano(static_cast<int64_t>(seconds * SEC_TO_NANO)) {}
122
123 bool operator()(const TofEvent &e1, const TofEvent &e2) {
124 // get the pulse times converted into bin number from start time
125 const int64_t e1Pulse = (e1.pulseTime().totalNanoseconds() - startNano) / deltaNano;
126 const int64_t e2Pulse = (e2.pulseTime().totalNanoseconds() - startNano) / deltaNano;
127
128 // compare with the calculated bin information
129 if (e1Pulse < e2Pulse) {
130 return true;
131 } else if ((e1Pulse == e2Pulse) && (e1.tof() < e2.tof())) {
132 return true;
133 }
134
135 return false;
136 }
137
138 int64_t startNano;
139 int64_t deltaNano;
140};
141
142struct FindBin {
143 double divisor;
144 double offset;
145 std::optional<size_t> (*findBin)(const Mantid::MantidVec &, const double, const double, const double, const bool);
146 FindBin(double step, double xmin) {
147 if (step < 0) {
149 divisor = 1. / log1p(abs(step)); // use this to do change of base
150 offset = log(xmin) * divisor;
151 } else {
153 divisor = 1. / step;
154 offset = xmin * divisor;
155 }
156 }
157
158 std::optional<size_t> operator()(const Mantid::MantidVec &X, const double tof, const bool findExact) {
159 return findBin(X, tof, divisor, offset, findExact);
160 }
161};
162
164// EventWorkspace is always histogram data and so is thus EventList
166 : m_histogram(HistogramData::Histogram::XMode::BinEdges, HistogramData::Histogram::YMode::Counts),
167 eventType(event_type), order(UNSORTED), mru(nullptr) {
168 switch (eventType) {
169 case TOF:
170 this->events = std::make_unique<std::vector<Mantid::Types::Event::TofEvent>>();
171 this->weightedEvents = nullptr;
172 this->weightedEventsNoTime = nullptr;
173 break;
174
175 case WEIGHTED:
176 this->events = nullptr;
177 this->weightedEvents = std::make_unique<std::vector<WeightedEvent>>();
178 this->weightedEventsNoTime = nullptr;
179 break;
180
181 case WEIGHTED_NOTIME:
182 this->events = nullptr;
183 this->weightedEventsNoTime = std::make_unique<std::vector<WeightedEventNoTime>>();
184 this->weightedEventsNoTime = nullptr;
185 break;
186 }
187}
188
194 : IEventList(specNo),
195 m_histogram(HistogramData::Histogram::XMode::BinEdges, HistogramData::Histogram::YMode::Counts),
196 weightedEvents(nullptr), weightedEventsNoTime(nullptr), eventType(TOF), order(UNSORTED), mru(mru) {
197 this->events = std::make_unique<std::vector<Mantid::Types::Event::TofEvent>>();
198}
199
202EventList::EventList(const EventList &rhs) : IEventList(rhs), m_histogram(rhs.m_histogram), mru{nullptr} {
203 // Note that operator= also assigns m_histogram, but the above use of the copy
204 // constructor avoid a memory allocation and is thus faster.
205 this->operator=(rhs);
206}
207
210EventList::EventList(const std::vector<Types::Event::TofEvent> &events)
211 : m_histogram(HistogramData::Histogram::XMode::BinEdges, HistogramData::Histogram::YMode::Counts),
212 weightedEvents(nullptr), weightedEventsNoTime(nullptr), eventType(TOF), mru(nullptr) {
213 this->events = std::make_unique<std::vector<Mantid::Types::Event::TofEvent>>(events.cbegin(), events.cend());
214 this->eventType = TOF;
215 this->order = UNSORTED;
216}
217
220EventList::EventList(const std::vector<WeightedEvent> &events)
221 : m_histogram(HistogramData::Histogram::XMode::BinEdges, HistogramData::Histogram::YMode::Counts), events(nullptr),
222 weightedEventsNoTime(nullptr), mru(nullptr) {
223 this->weightedEvents = std::make_unique<std::vector<WeightedEvent>>(events.cbegin(), events.cend());
224 this->eventType = WEIGHTED;
225 this->order = UNSORTED;
226}
227
230EventList::EventList(const std::vector<WeightedEventNoTime> &events)
231 : m_histogram(HistogramData::Histogram::XMode::BinEdges, HistogramData::Histogram::YMode::Counts), events(nullptr),
232 weightedEvents(nullptr), mru(nullptr) {
233 this->weightedEventsNoTime = std::make_unique<std::vector<WeightedEventNoTime>>(events.cbegin(), events.cend());
235 this->order = UNSORTED;
236}
237
240 // clear this out of the MRU (copy of code from EventList::clear()
241 if (mru) {
242 try {
243 mru->deleteIndex(this);
244 } catch (const std::runtime_error &) {
245 // this is an ignorable error
246 }
247 }
248
249 // set all member vectors to nullptr
250 this->events.reset();
251 this->weightedEvents.reset();
252 this->weightedEventsNoTime.reset();
253}
254
256void EventList::copyDataFrom(const ISpectrum &source) { source.copyDataInto(*this); }
257
261 if (events)
262 sink.events = std::make_unique<std::vector<Types::Event::TofEvent>>(events->cbegin(), events->cend());
263 else if (sink.events)
264 sink.events = std::make_unique<std::vector<Types::Event::TofEvent>>();
265 if (weightedEvents)
266 sink.weightedEvents =
267 std::make_unique<std::vector<WeightedEvent>>(weightedEvents->cbegin(), weightedEvents->cend());
268 else if (sink.weightedEvents)
269 sink.weightedEvents = std::make_unique<std::vector<WeightedEvent>>();
271 sink.weightedEventsNoTime = std::make_unique<std::vector<WeightedEventNoTime>>(weightedEventsNoTime->cbegin(),
272 weightedEventsNoTime->cend());
273 else if (sink.weightedEventsNoTime)
274 sink.weightedEventsNoTime = std::make_unique<std::vector<WeightedEventNoTime>>();
275
276 sink.eventType = eventType;
277 sink.order = order;
278}
279
282
283// --------------------------------------------------------------------------
295void EventList::createFromHistogram(const ISpectrum *inSpec, bool GenerateZeros, bool GenerateMultipleEvents,
296 int MaxEventsPerBin) {
297 // Fresh start
298 this->clear(true);
299
300 // Get the input histogram
301 Mantid::HistogramData::HistogramX const &X = inSpec->x();
302 Mantid::HistogramData::HistogramY const &Y = inSpec->y();
303 Mantid::HistogramData::HistogramE const &E = inSpec->e();
304 if (Y.size() + 1 != X.size()) {
305 throw std::runtime_error("Expected a histogram (X vector should be 1 longer than the Y vector)");
306 }
307
308 // Copy detector IDs and spectra
309 this->copyInfoFrom(*inSpec);
310 // We need weights but have no way to set the time. So use weighted, no time
312 if (GenerateZeros)
313 this->weightedEventsNoTime->reserve(Y.size());
314
315 for (size_t i = 0; i < X.size() - 1; i++) {
316 double weight = Y[i];
317 if ((weight != 0.0 || GenerateZeros) && std::isfinite(weight)) {
318 double error = E[i];
319 // Also check that the error is not a bad number
320 if (std::isfinite(error)) {
321 if (GenerateMultipleEvents) {
322 // --------- Multiple events per bin ----------
323 double errorSquared = error * error;
324 // Find how many events to fake
325 double val = weight / E[i];
326 val *= val;
327 // Convert to int with slight rounding up. This is to avoid rounding
328 // errors
329 auto numEvents = int(val + 0.2);
330 if (numEvents < 1)
331 numEvents = 1;
332 if (numEvents > MaxEventsPerBin)
333 numEvents = MaxEventsPerBin;
334 // Scale the weight and error for each
335 weight /= numEvents;
336 errorSquared /= numEvents;
337
338 // Spread the TOF. e.g. 2 events = 0.25, 0.75.
339 double tofStep = (X[i + 1] - X[i]) / (numEvents);
340 for (size_t j = 0; j < size_t(numEvents); j++) {
341 double tof = X[i] + tofStep * (0.5 + double(j));
342 // Create and add the event
343 // TODO: try emplace_back() here.
344 weightedEventsNoTime->emplace_back(tof, weight, errorSquared);
345 }
346 } else {
347 // --------- Single event per bin ----------
348 // TOF = midpoint of the bin
349 double tof = (X[i] + X[i + 1]) / 2.0;
350 // Error squared is carried in the event
351 double errorSquared = E[i];
352 errorSquared *= errorSquared;
353 // Create and add the event
354 weightedEventsNoTime->emplace_back(tof, weight, errorSquared);
355 }
356 } // error is nont NAN or infinite
357 } // weight is non-zero, not NAN, and non-infinite
358 } // (each bin)
359
360 // Set the X binning parameters
361 this->setSharedX(inSpec->sharedX());
362
363 // Manually set that this is sorted by TOF, since it is. This will make it
364 // "threadSafe" in other algos.
365 this->setSortOrder(TOF_SORT);
366}
367
368// --------------------------------------------------------------------------
369// --- Operators
370// -------------------------------------------------------------------
371
377 // Note that we are NOT copying the MRU pointer
378 // the EventWorkspace that possesses the EventList has already configured the mru
379 IEventList::operator=(rhs);
380 m_histogram = rhs.m_histogram;
381 rhs.copyDataInto(*this);
382 return *this;
383}
384
385// --------------------------------------------------------------------------
390EventList &EventList::operator+=(const Types::Event::TofEvent &event) {
391
392 switch (this->eventType) {
393 case TOF:
394 // Simply push the events
395 this->events->emplace_back(event);
396 break;
397
398 case WEIGHTED:
399 this->weightedEvents->emplace_back(event);
400 break;
401
402 case WEIGHTED_NOTIME:
403 this->weightedEventsNoTime->emplace_back(event);
404 break;
405 }
406
407 this->order = UNSORTED;
408 return *this;
409}
410
411// --------------------------------------------------------------------------
418EventList &EventList::operator+=(const std::vector<Types::Event::TofEvent> &more_events) {
419 switch (this->eventType) {
420 case TOF:
421 // Simply push the events
422 this->events->insert(this->events->end(), more_events.cbegin(), more_events.cend());
423 break;
424
425 case WEIGHTED:
426 // Add default weights to all the un-weighted incoming events from the list.
427 // and append to the list
428 this->weightedEvents->reserve(this->weightedEvents->size() + more_events.size());
429 std::copy(more_events.cbegin(), more_events.cend(), std::back_inserter(*this->weightedEvents));
430 break;
431
432 case WEIGHTED_NOTIME:
433 // Add default weights to all the un-weighted incoming events from the list.
434 // and append to the list
435 this->weightedEventsNoTime->reserve(this->weightedEventsNoTime->size() + more_events.size());
436 std::copy(more_events.cbegin(), more_events.cend(), std::back_inserter(*this->weightedEventsNoTime));
437 break;
438 }
439
440 this->order = UNSORTED;
441 return *this;
442}
443
444// --------------------------------------------------------------------------
453 this->switchTo(WEIGHTED);
454 this->weightedEvents->emplace_back(event);
455 this->order = UNSORTED;
456 return *this;
457}
458
459// --------------------------------------------------------------------------
467EventList &EventList::operator+=(const std::vector<WeightedEvent> &more_events) {
468 switch (this->eventType) {
469 case TOF:
470 // Need to switch to weighted
471 this->switchTo(WEIGHTED);
472 // Fall through to the insertion!
473
474 case WEIGHTED:
475 // Append the two lists
476 this->weightedEvents->insert(weightedEvents->end(), more_events.cbegin(), more_events.cend());
477 break;
478
479 case WEIGHTED_NOTIME:
480 // Add default weights to all the un-weighted incoming events from the list.
481 // and append to the list
482 this->weightedEventsNoTime->reserve(this->weightedEventsNoTime->size() + more_events.size());
483 std::copy(more_events.cbegin(), more_events.cend(), std::back_inserter(*this->weightedEventsNoTime));
484 break;
485 }
486
487 this->order = UNSORTED;
488 return *this;
489}
490
491// --------------------------------------------------------------------------
499EventList &EventList::operator+=(const std::vector<WeightedEventNoTime> &more_events) {
500 switch (this->eventType) {
501 case TOF:
502 case WEIGHTED:
503 // Need to switch to weighted with no time
505 // Fall through to the insertion!
506
507 case WEIGHTED_NOTIME:
508 // Simple appending of the two lists
509 this->weightedEventsNoTime->insert(weightedEventsNoTime->end(), more_events.cbegin(), more_events.cend());
510 break;
511 }
512
513 this->order = UNSORTED;
514 return *this;
515}
516
517// --------------------------------------------------------------------------
527 if (!more_events.empty()) {
528 // We'll let the += operator for the given vector of event lists handle it
529 switch (more_events.getEventType()) {
530 case TOF:
531 this->operator+=(*more_events.events);
532 break;
533
534 case WEIGHTED:
535 this->operator+=(*more_events.weightedEvents);
536 break;
537
538 case WEIGHTED_NOTIME:
539 this->operator+=(*more_events.weightedEventsNoTime);
540 break;
541 }
542
543 // No guaranteed order
544 if (this->empty()) {
545 this->order = more_events.order;
546 } else {
547 this->order = UNSORTED;
548 }
549 }
550
551 // Do a union between the detector IDs of both lists
552 addDetectorIDs(more_events.getDetectorIDs());
553
554 return *this;
555}
556
557// --------------------------------------------------------------------------
566template <class T1, class T2> void EventList::minusHelper(std::vector<T1> &events, const std::vector<T2> &more_events) {
567 // Make the end vector big enough in one go (avoids repeated re-allocations).
568 events.reserve(events.size() + more_events.size());
569 /* In the event of subtracting in place, calling the end() vector would make
570 * it point at the wrong place
571 * Using it caused a segault, Ticket #2306.
572 * So we cache the end (this speeds up too).
573 */
574 // We call the constructor for T1. In the case of WeightedEventNoTime, the pulse time will just be ignored.
575 std::transform(more_events.cbegin(), more_events.cend(), std::back_inserter(events),
576 [](const auto &ev) { return T1(ev.tof(), ev.pulseTime(), ev.weight() * (-1.0), ev.errorSquared()); });
577}
578
579// --------------------------------------------------------------------------
588 if (this == &more_events) {
589 // Special case, ticket #3844 part 2.
590 // When doing this = this - this,
591 // simply clear the input event list. Saves memory!
592 this->clearData();
593 return *this;
594 }
595
596 // We'll let the -= operator for the given vector of event lists handle it
597 switch (this->getEventType()) {
598 case TOF:
599 this->switchTo(WEIGHTED);
600 // Fall through
601
602 case WEIGHTED:
603 switch (more_events.getEventType()) {
604 case TOF:
605 minusHelper(*this->weightedEvents, *more_events.events);
606 break;
607 case WEIGHTED:
608 minusHelper(*this->weightedEvents, *more_events.weightedEvents);
609 break;
610 case WEIGHTED_NOTIME:
611 // TODO: Should this throw?
612 minusHelper(*this->weightedEvents, *more_events.weightedEventsNoTime);
613 break;
614 }
615 break;
616
617 case WEIGHTED_NOTIME:
618 switch (more_events.getEventType()) {
619 case TOF:
620 minusHelper(*this->weightedEventsNoTime, *more_events.events);
621 break;
622 case WEIGHTED:
623 minusHelper(*this->weightedEventsNoTime, *more_events.weightedEvents);
624 break;
625 case WEIGHTED_NOTIME:
627 break;
628 }
629 break;
630 }
631
632 // No guaranteed order
633 this->order = UNSORTED;
634
635 // NOTE: What to do about detector ID's?
636 return *this;
637}
638
639namespace {
640/*
641 * Both can be nullptr, or the values can be equal, but do not have one nullptr
642 */
643template <typename T>
644bool vectorPtrEquals(const std::unique_ptr<std::vector<T>> &left, const std::unique_ptr<std::vector<T>> &right) {
645 if (left && right) {
646 return (*left == *right);
647 ;
648 } else if ((left && !right) || (right && !left)) {
649 return false;
650 }
651 return true;
652}
653} // anonymous namespace
654
655// --------------------------------------------------------------------------
661 if (this->getNumberEvents() != rhs.getNumberEvents())
662 return false;
663 if (this->eventType != rhs.eventType)
664 return false;
665 if (this->empty())
666 return true;
667 // Check all event lists; The empty ones will compare equal
668 if (!vectorPtrEquals(events, rhs.events))
669 return false;
670 if (!vectorPtrEquals(weightedEvents, rhs.weightedEvents))
671 return false;
672 if (!vectorPtrEquals(weightedEventsNoTime, rhs.weightedEventsNoTime))
673 return false;
674
675 // nothing wasn't equal, so they are equal
676 return true;
677}
678
683bool EventList::operator!=(const EventList &rhs) const { return (!this->operator==(rhs)); }
684
685bool EventList::equals(const EventList &rhs, const double tolTof, const double tolWeight,
686 const int64_t tolPulse) const {
687 // generic checks
688 if (this->getNumberEvents() != rhs.getNumberEvents())
689 return false;
690 if (this->eventType != rhs.eventType)
691 return false;
692 if (this->empty())
693 return true;
694
695 // loop over the events
696 switch (this->eventType) {
697 case TOF: {
698 auto leftIter = this->events->cbegin();
699 auto leftEnd = this->events->cend();
700 auto rightIter = rhs.events->cbegin();
701 while (leftIter != leftEnd) {
702 if (!leftIter->equals(*rightIter, tolTof, tolPulse))
703 return false;
704 leftIter = std::next(leftIter);
705 rightIter = std::next(rightIter);
706 }
707 break;
708 }
709 case WEIGHTED: {
710 auto leftIter = this->weightedEvents->cbegin();
711 auto leftEnd = this->weightedEvents->cend();
712 auto rightIter = rhs.weightedEvents->cbegin();
713 while (leftIter != leftEnd) {
714 if (!leftIter->equals(*rightIter, tolTof, tolWeight, tolPulse))
715 return false;
716 leftIter = std::next(leftIter);
717 rightIter = std::next(rightIter);
718 }
719 break;
720 }
721 case WEIGHTED_NOTIME: {
722 auto leftIter = this->weightedEventsNoTime->cbegin();
723 auto leftEnd = this->weightedEventsNoTime->cend();
724 auto rightIter = rhs.weightedEventsNoTime->cbegin();
725 while (leftIter != leftEnd) {
726 if (!leftIter->equals(*rightIter, tolTof, tolWeight))
727 return false;
728 leftIter = std::next(leftIter);
729 rightIter = std::next(rightIter);
730 }
731 break;
732 }
733 default:
734 break;
735 }
736
737 // anything that gets this far is equal within tolerances
738 return true;
739}
740
741// -----------------------------------------------------------------------------------------------
746
747// -----------------------------------------------------------------------------------------------
752 switch (newType) {
753 case TOF:
754 if (eventType != TOF)
755 throw std::runtime_error("EventList::switchTo() called on an EventList with weights to go down to TofEvent's. "
756 "This would remove weight information and therefore is not possible.");
757 break;
758
759 case WEIGHTED:
761 break;
762
763 case WEIGHTED_NOTIME:
765 break;
766 }
767 // Make sure to free memory
768 this->clearUnused();
769}
770
771// -----------------------------------------------------------------------------------------------
776 switch (eventType) {
777 case WEIGHTED:
778 // Do nothing; it already is weighted
779 return;
780
781 case WEIGHTED_NOTIME:
782 throw std::runtime_error("EventList::switchToWeightedEvents() called on an EventList with WeightedEventNoTime's. "
783 "It has lost the pulse time information and can't go back to WeightedEvent's.");
784 break;
785
786 case TOF:
787 if (events && !events->empty()) {
788 // Convert and copy all TofEvents to the weightedEvents list.
789 weightedEvents = std::make_unique<std::vector<WeightedEvent>>(events->cbegin(), events->cend());
790 // Get rid of the old events
791 events.reset();
792 } else {
793 weightedEvents = std::make_unique<std::vector<WeightedEvent>>();
794 }
796 break;
797 }
798}
799
800// -----------------------------------------------------------------------------------------------
805 switch (eventType) {
806 case WEIGHTED_NOTIME:
807 // Do nothing if already there
808 return;
809
810 case TOF: {
811 if (events && !events->empty()) {
812 // Convert and copy all TofEvents to the weightedEvents list.
813 this->weightedEventsNoTime = std::make_unique<std::vector<WeightedEventNoTime>>(events->cbegin(), events->cend());
814 // Get rid of the old events
815 events.reset();
816 } else {
817 this->weightedEventsNoTime = std::make_unique<std::vector<WeightedEventNoTime>>();
818 }
819 break;
820 }
821
822 case WEIGHTED: {
823 // Convert and copy all TofEvents to the weightedEvents list.
824 if (weightedEvents && !weightedEvents->empty()) {
826 std::make_unique<std::vector<WeightedEventNoTime>>(weightedEvents->cbegin(), weightedEvents->cend());
827 // Get rid of the old events
828 weightedEvents.reset();
829 } else {
830 this->weightedEventsNoTime = std::make_unique<std::vector<WeightedEventNoTime>>();
831 }
832 break;
833 }
834 }
836}
837
838// ==============================================================================================
839// --- Testing functions (mostly)
840// ---------------------------------------------------------------
841// ==============================================================================================
842
849WeightedEvent EventList::getEvent(size_t event_number) {
850 switch (eventType) {
851 case TOF:
852 return WeightedEvent(events->at(event_number));
853 case WEIGHTED:
854 return weightedEvents->at(event_number);
855 case WEIGHTED_NOTIME: {
856 const auto event = weightedEventsNoTime->at(event_number);
857 return WeightedEvent(event.tof(), 0, event.weight(), event.errorSquared());
858 }
859 }
860 throw std::runtime_error("EventList: invalid event type value was found.");
861}
862
863// ==============================================================================================
864// --- Handling the event list
865// -------------------------------------------------------------------
866// ==============================================================================================
867
875const std::vector<TofEvent> &EventList::getEvents() const {
876 if (eventType != TOF)
877 throw std::runtime_error("EventList::getEvents() called for an EventList that has weights. Use getWeightedEvents() "
878 "or getWeightedEventsNoTime().");
879 if (this->events)
880 return *this->events;
881 else
882 throw std::runtime_error("unweighted event vector is not initialized");
883}
884
891std::vector<TofEvent> &EventList::getEvents() {
892 if (eventType != TOF)
893 throw std::runtime_error("EventList::getEvents() called for an EventList that has weights. Use getWeightedEvents() "
894 "or getWeightedEventsNoTime().");
895 if (this->events)
896 return *this->events;
897 else
898 throw std::runtime_error("unweighted event vector is not initialized");
899}
900
908std::vector<WeightedEvent> &EventList::getWeightedEvents() {
909 if (eventType != WEIGHTED)
910 throw std::runtime_error("EventList::getWeightedEvents() called for an EventList not of type WeightedEvent. Use "
911 "getEvents() or getWeightedEventsNoTime().");
912 if (this->weightedEvents)
913 return *this->weightedEvents;
914 else
915 throw std::runtime_error("weighted event vector is not initialized");
916}
917
925const std::vector<WeightedEvent> &EventList::getWeightedEvents() const {
926 if (eventType != WEIGHTED)
927 throw std::runtime_error("EventList::getWeightedEvents() called for an EventList not of type WeightedEvent. Use "
928 "getEvents() or getWeightedEventsNoTime().");
929 if (this->weightedEvents)
930 return *this->weightedEvents;
931 else
932 throw std::runtime_error("weighted event vector is not initialed");
933}
934
940std::vector<WeightedEventNoTime> &EventList::getWeightedEventsNoTime() {
942 throw std::runtime_error("EventList::getWeightedEventsNoTime() called for an EventList not of type "
943 "WeightedEventNoTime. Use getEvents() or getWeightedEvents().");
944 if (this->weightedEventsNoTime)
945 return *this->weightedEventsNoTime;
946 else
947 throw std::runtime_error("weighted event no time vector is not initialed");
948}
949
955const std::vector<WeightedEventNoTime> &EventList::getWeightedEventsNoTime() const {
957 throw std::runtime_error("EventList::getWeightedEventsNoTime() called for an EventList not of type "
958 "WeightedEventNoTime. Use getEvents() or getWeightedEvents().");
959 if (this->weightedEventsNoTime)
960 return *this->weightedEventsNoTime;
961 else
962 throw std::runtime_error("weighted event no time vector is not initialed");
963}
964
968void EventList::clear(const bool removeDetIDs) {
969 if (mru) {
970 try {
971 mru->deleteIndex(this);
972 } catch (const std::runtime_error &) {
973 // this is an ignorable error
974 }
975 }
976 // clear representations that aren't for the current type
977 this->clearUnused();
978
979 // release unused memory or allocate new vector
980 // rather than creating a new object, reset existing pointer
981 if (!this->empty()) {
982 if (this->events && eventType == TOF) {
983 this->events->clear();
984 std::vector<TofEvent>().swap(*this->events); // STL Trick to release memory
985 }
986 if (this->weightedEvents && eventType == WEIGHTED) {
987 this->weightedEvents->clear();
988 std::vector<WeightedEvent>().swap(*this->weightedEvents); // STL Trick to release memory
989 }
991 this->weightedEventsNoTime->clear();
992 std::vector<WeightedEventNoTime>().swap(*this->weightedEventsNoTime); // STL Trick to release memory
993 }
994 }
995 if (removeDetIDs)
996 this->clearDetectorIDs();
997}
998
1004 if (eventType != TOF && (this->events)) {
1005 this->events.reset();
1006 }
1007 if (eventType != WEIGHTED && (this->weightedEvents)) {
1008 this->weightedEvents.reset();
1009 }
1011 this->weightedEventsNoTime.reset();
1012 }
1013}
1014
1016void EventList::clearData() { this->clear(false); }
1017
1022void EventList::setMRU(EventWorkspaceMRU *newMRU) { mru = newMRU; }
1023
1031void EventList::reserve(size_t num) {
1032 switch (this->eventType) {
1033 case TOF:
1034 this->events->reserve(num);
1035 break;
1036 case WEIGHTED:
1037 this->weightedEvents->reserve(num);
1038 break;
1039 case WEIGHTED_NOTIME:
1040 this->weightedEventsNoTime->reserve(num);
1041 break;
1042 }
1043}
1044
1045// ==============================================================================================
1046// --- Sorting functions -----------------------------------------------------
1047// ==============================================================================================
1048
1049// --------------------------------------------------------------------------
1053void EventList::sort(const EventSortType order) const {
1054 if (order == UNSORTED) {
1055 return; // don't bother doing anything. Why did you ask to unsort?
1056 } else if (order == TOF_SORT) {
1057 this->sortTof();
1058 } else if (order == PULSETIME_SORT) {
1059 this->sortPulseTime();
1060 } else if (order == PULSETIMETOF_SORT) {
1061 this->sortPulseTimeTOF();
1062 } else if (order == PULSETIMETOF_DELTA_SORT) {
1063 throw std::invalid_argument("sorting by pulse time with delta requires "
1064 "extra parameters. Use sortPulseTimeTOFDelta "
1065 "instead.");
1066 } else if (order == TIMEATSAMPLE_SORT) {
1067 throw std::invalid_argument("sorting by time at sample requires extra "
1068 "parameters. Use sortTimeAtSample instead.");
1069 } else {
1070 throw runtime_error("Invalid sort type in EventList::sort(EventSortType)");
1071 }
1072}
1073
1074// --------------------------------------------------------------------------
1079void EventList::setSortOrder(const EventSortType order) const { this->order = order; }
1080
1081namespace {
1082// these are abstractions
1083template <class RandomIt> void switchable_sort(RandomIt first, RandomIt last) {
1084 const auto vec_size = static_cast<size_t>(std::distance(first, last));
1085 if (vec_size < 2)
1086 return;
1087 else if (vec_size < MIN_VEC_LENGTH_PARALLEL_SORT)
1088 std::sort(first, last);
1089 else
1090 tbb::parallel_sort(first, last);
1091}
1092
1093template <class RandomIt, class Compare> void switchable_sort(RandomIt first, RandomIt last, Compare comp) {
1094 const auto vec_size = static_cast<size_t>(std::distance(first, last));
1095 if (vec_size < 2)
1096 return;
1097 else if (vec_size < MIN_VEC_LENGTH_PARALLEL_SORT)
1098 std::sort(first, last, std::move(comp));
1099 else
1100 tbb::parallel_sort(first, last, comp);
1101}
1102} // anonymous namespace
1103
1104// --------------------------------------------------------------------------
1107 // nothing to do
1108 if (this->order == TOF_SORT)
1109 return;
1110
1111 // Avoid sorting from multiple threads
1112 std::lock_guard<std::mutex> _lock(m_sortMutex);
1113 // If the list was sorted while waiting for the lock, return.
1114 if (this->order == TOF_SORT) // cppcheck-suppress identicalConditionAfterEarlyExit
1115 return;
1116
1117 switch (eventType) {
1118 case TOF:
1119 switchable_sort(events->begin(), events->end());
1120 break;
1121 case WEIGHTED:
1122 switchable_sort(weightedEvents->begin(), weightedEvents->end());
1123 break;
1124 case WEIGHTED_NOTIME:
1125 switchable_sort(weightedEventsNoTime->begin(), weightedEventsNoTime->end());
1126 break;
1127 }
1128 // Save the order to avoid unnecessary re-sorting.
1129 this->order = TOF_SORT;
1130}
1131
1132// --------------------------------------------------------------------------
1141void EventList::sortTimeAtSample(const double &tofFactor, const double &tofShift, bool forceResort) const {
1142 // Check pre-cached sort flag.
1143 if (this->order == TIMEATSAMPLE_SORT && !forceResort)
1144 return;
1145
1146 // Avoid sorting from multiple threads
1147 std::lock_guard<std::mutex> _lock(m_sortMutex);
1148 // If the list was sorted while waiting for the lock, return.
1149 if (this->order == TIMEATSAMPLE_SORT && !forceResort)
1150 return;
1151
1152 // Perform sort.
1153 switch (eventType) {
1154 case TOF: {
1155 CompareTimeAtSample<TofEvent> comparitor(tofFactor, tofShift);
1156 switchable_sort(events->begin(), events->end(), comparitor);
1157 } break;
1158 case WEIGHTED: {
1159 CompareTimeAtSample<WeightedEvent> comparitor(tofFactor, tofShift);
1160 switchable_sort(weightedEvents->begin(), weightedEvents->end(), comparitor);
1161 } break;
1162 case WEIGHTED_NOTIME: {
1163 CompareTimeAtSample<WeightedEventNoTime> comparitor(tofFactor, tofShift);
1164 switchable_sort(weightedEventsNoTime->begin(), weightedEventsNoTime->end(), comparitor);
1165 } break;
1166 }
1167 // Save the order to avoid unnecessary re-sorting.
1168 this->order = TIMEATSAMPLE_SORT;
1169}
1170
1171// --------------------------------------------------------------------------
1174 if (this->order == PULSETIME_SORT || this->order == PULSETIMETOF_SORT)
1175 return; // nothing to do
1176
1177 // Avoid sorting from multiple threads
1178 std::lock_guard<std::mutex> _lock(m_sortMutex);
1179 // If the list was sorted while waiting for the lock, return.
1180 if (this->order == PULSETIME_SORT)
1181 return;
1182
1183 // Perform sort.
1184 switch (eventType) {
1185 case TOF:
1186 switchable_sort(events->begin(), events->end(), compareEventPulseTime);
1187 break;
1188 case WEIGHTED:
1189 switchable_sort(weightedEvents->begin(), weightedEvents->end(), compareEventPulseTime);
1190 break;
1191 case WEIGHTED_NOTIME:
1192 // Do nothing; there is no time to sort
1193 break;
1194 }
1195 // Save the order to avoid unnecessary re-sorting.
1196 this->order = PULSETIME_SORT;
1197}
1198
1199/*
1200 * Sort events by pulse time + TOF
1201 * (the absolute time)
1202 */
1204 if (this->order == PULSETIMETOF_SORT)
1205 return; // already ordered
1206
1207 // Avoid sorting from multiple threads
1208 std::lock_guard<std::mutex> _lock(m_sortMutex);
1209 // If the list was sorted while waiting for the lock, return.
1210 if (this->order == PULSETIMETOF_SORT) // cppcheck-suppress identicalConditionAfterEarlyExit
1211 return;
1212
1213 switch (eventType) {
1214 case TOF:
1215 switchable_sort(events->begin(), events->end(), compareEventPulseTimeTOF);
1216 break;
1217 case WEIGHTED:
1218 switchable_sort(weightedEvents->begin(), weightedEvents->end(), compareEventPulseTimeTOF);
1219 break;
1220 case WEIGHTED_NOTIME:
1221 // Do nothing; there is no time to sort
1222 break;
1223 }
1224
1225 // Save
1226 this->order = PULSETIMETOF_SORT;
1227}
1228
1236void EventList::sortPulseTimeTOFDelta(const Types::Core::DateAndTime &start, const double seconds) const {
1237 // Avoid sorting from multiple threads
1238 std::lock_guard<std::mutex> _lock(m_sortMutex);
1239
1240 std::function<bool(const TofEvent &, const TofEvent &)> comparator = comparePulseTimeTOFDelta(start, seconds);
1241
1242 switch (eventType) {
1243 case TOF:
1244 switchable_sort(events->begin(), events->end(), std::move(comparator));
1245 break;
1246 case WEIGHTED:
1247 switchable_sort(weightedEvents->begin(), weightedEvents->end(), std::move(comparator));
1248 break;
1249 case WEIGHTED_NOTIME:
1250 // Do nothing; there is no time to sort
1251 break;
1252 }
1253
1254 this->order = UNSORTED; // so the function always re-runs
1255}
1256
1257// --------------------------------------------------------------------------
1259bool EventList::isSortedByTof() const { return (this->order == TOF_SORT); }
1260
1261// --------------------------------------------------------------------------
1264
1265// --------------------------------------------------------------------------
1272 // reverse the histogram bin parameters
1273 MantidVec &x = dataX();
1274 std::reverse(x.begin(), x.end());
1275
1276 // flip the events if they are tof sorted
1277 if (this->isSortedByTof()) {
1278 switch (eventType) {
1279 case TOF:
1280 std::reverse(this->events->begin(), this->events->end());
1281 break;
1282 case WEIGHTED:
1283 std::reverse(this->weightedEvents->begin(), this->weightedEvents->end());
1284 break;
1285 case WEIGHTED_NOTIME:
1286 std::reverse(this->weightedEventsNoTime->begin(), this->weightedEventsNoTime->end());
1287 break;
1288 }
1289 // And we are still sorted! :)
1290 }
1291 // Otherwise, do nothing. If it was sorted by pulse time, then it still is
1292}
1293
1294// --------------------------------------------------------------------------
1303 switch (eventType) {
1304 case TOF:
1305 return (this->events) ? this->events->size() : 0;
1306 case WEIGHTED:
1307 return (this->weightedEvents) ? this->weightedEvents->size() : 0;
1308 case WEIGHTED_NOTIME:
1309 return (this->weightedEventsNoTime) ? this->weightedEventsNoTime->size() : 0;
1310 }
1311 throw std::runtime_error("EventList: invalid event type value was found.");
1312}
1313
1317bool EventList::empty() const {
1318 switch (eventType) {
1319 case TOF:
1320 if (this->events)
1321 return this->events->empty();
1322 else
1323 throw std::runtime_error("TOF events is nullptr");
1324 case WEIGHTED:
1325 if (this->weightedEvents)
1326 return this->weightedEvents->empty();
1327 else
1328 throw std::runtime_error("WEIGHTED events is nullptr");
1329 case WEIGHTED_NOTIME:
1330 if (this->weightedEventsNoTime)
1331 return this->weightedEventsNoTime->empty();
1332 else
1333 throw std::runtime_error("WEIGHTED_NOTIME events is nullptr");
1334 }
1335 throw std::runtime_error("EventList: invalid event type value was found.");
1336}
1337
1338// --------------------------------------------------------------------------
1346 switch (eventType) {
1347 case TOF:
1348 return this->events->capacity() * sizeof(TofEvent) + sizeof(EventList);
1349 case WEIGHTED:
1350 return this->weightedEvents->capacity() * sizeof(WeightedEvent) + sizeof(EventList);
1351 case WEIGHTED_NOTIME:
1352 return this->weightedEventsNoTime->capacity() * sizeof(WeightedEventNoTime) + sizeof(EventList);
1353 }
1354 throw std::runtime_error("EventList: invalid event type value was found.");
1355}
1356
1357// --------------------------------------------------------------------------
1361 size_t x_size = x().size();
1362 if (x_size > 1)
1363 return x_size - 1;
1364 else
1365 return 0;
1366}
1367
1368// ==============================================================================================
1369// --- Setting the Histogram X axis, without recalculating the histogram
1370// -----------------------
1371// ==============================================================================================
1372
1378 m_histogram.setSharedX(X);
1379 if (mru)
1380 mru->deleteIndex(this);
1381}
1382
1387 if (mru)
1388 mru->deleteIndex(this);
1389 return m_histogram.dataX();
1390}
1391
1394const MantidVec &EventList::dataX() const { return m_histogram.dataX(); }
1395
1397const MantidVec &EventList::readX() const { return m_histogram.x().rawData(); }
1398
1401
1405const MantidVec &EventList::dataDx() const { return m_histogram.dataDx(); }
1407const MantidVec &EventList::readDx() const { return m_histogram.readDx(); }
1408
1409// ==============================================================================================
1410// --- Return Data Vectors --------------------------------------------------
1411// ==============================================================================================
1412
1419 auto Y = new MantidVec();
1420 MantidVec E;
1421 // Generate the Y histogram while skipping the E if possible.
1422 generateHistogram(readX(), *Y, E, true);
1423 return Y;
1424}
1425
1432 MantidVec Y;
1433 auto E = new MantidVec();
1434 generateHistogram(readX(), Y, *E);
1435 // Y is unused.
1436 return E;
1437}
1438
1440HistogramData::Histogram EventList::getHistogram() const { return m_histogram; }
1441
1442HistogramData::Histogram EventList::histogram() const {
1443 HistogramData::Histogram ret(m_histogram);
1444 ret.setSharedY(sharedY());
1445 ret.setSharedE(sharedE());
1446 return ret;
1447}
1448
1449HistogramData::Counts EventList::counts() const { return histogram().counts(); }
1450
1451HistogramData::CountVariances EventList::countVariances() const { return histogram().countVariances(); }
1452
1453HistogramData::CountStandardDeviations EventList::countStandardDeviations() const {
1454 return histogram().countStandardDeviations();
1455}
1456
1457HistogramData::Frequencies EventList::frequencies() const { return histogram().frequencies(); }
1458
1459HistogramData::FrequencyVariances EventList::frequencyVariances() const { return histogram().frequencyVariances(); }
1460
1461HistogramData::FrequencyStandardDeviations EventList::frequencyStandardDeviations() const {
1462 return histogram().frequencyStandardDeviations();
1463}
1464
1465const HistogramData::HistogramY &EventList::y() const {
1466 if (!mru)
1467 throw std::runtime_error("'EventList::y()' called with no MRU set. This is not allowed.");
1468
1469 return *sharedY();
1470}
1471const HistogramData::HistogramE &EventList::e() const {
1472 if (!mru)
1473 throw std::runtime_error("'EventList::e()' called with no MRU set. This is not allowed.");
1474
1475 return *sharedE();
1476}
1478 // This is the thread number from which this function was called.
1479 const int thread = PARALLEL_THREAD_NUMBER;
1480
1482
1483 // Is the data in the mrulist?
1484 if (mru) {
1485 mru->ensureEnoughBuffersY(static_cast<size_t>(thread));
1486 yData = mru->findY(static_cast<size_t>(thread), this);
1487 }
1488
1489 if (!yData) {
1490 MantidVec Y;
1491 MantidVec E;
1492 this->generateHistogram(readX(), Y, E);
1493
1494 // Create the MRU object
1495 yData = Kernel::make_cow<HistogramData::HistogramY>(std::move(Y));
1496
1497 // Lets save it in the MRU
1498 if (mru) {
1499 mru->insertY(thread, yData, this);
1500 auto eData = Kernel::make_cow<HistogramData::HistogramE>(std::move(E));
1501 mru->ensureEnoughBuffersE(thread);
1502 mru->insertE(thread, eData, this);
1503 }
1504 }
1505 return yData;
1506}
1508 // This is the thread number from which this function was called.
1509 const auto thread = static_cast<size_t>(PARALLEL_THREAD_NUMBER);
1510
1512
1513 // Is the data in the mrulist?
1514 if (mru) {
1515 mru->ensureEnoughBuffersE(thread);
1516 eData = mru->findE(thread, this);
1517 }
1518
1519 if (!eData) {
1520 // Now use that to get E -- Y values are generated from another function
1521 MantidVec Y_ignored;
1522 MantidVec E;
1523 this->generateHistogram(readX(), Y_ignored, E);
1524 eData = Kernel::make_cow<HistogramData::HistogramE>(std::move(E));
1525
1526 // Lets save it in the MRU
1527 if (mru)
1528 mru->insertE(thread, eData, this);
1529 }
1530 return eData;
1531}
1538 if (!mru)
1539 throw std::runtime_error("'EventList::dataY()' called with no MRU set. This is not allowed.");
1540
1541 // WARNING: The Y data of sharedY() is stored in MRU, returning reference fine
1542 // as long as it stays there.
1543 return sharedY()->rawData();
1544}
1545
1552 if (!mru)
1553 throw std::runtime_error("'EventList::dataE()' called with no MRU set. This is not allowed.");
1554
1555 // WARNING: The E data of sharedE() is stored in MRU, returning reference fine
1556 // as long as it stays there.
1557 return sharedE()->rawData();
1558}
1559
1560namespace {
1561inline double calcNorm(const double errorSquared) {
1562 if (errorSquared == 0.)
1563 return 0;
1564 else if (errorSquared == 1.)
1565 return 1.;
1566 else
1567 return 1. / std::sqrt(errorSquared);
1568}
1569} // namespace
1570
1571// --------------------------------------------------------------------------
1580template <class T>
1581inline void EventList::compressEventsHelper(const std::vector<T> &events, std::vector<WeightedEventNoTime> &out,
1582 double tolerance) {
1583 // Clear the output. We can't know ahead of time how much space to reserve :(
1584 out.clear();
1585 // We will make a starting guess of 1/20th of the number of input events.
1586 out.reserve(events.size() / 20);
1587
1588 // The last TOF to which we are comparing.
1589 double lastTof = events.front().m_tof;
1590 // For getting an accurate average TOF
1591 double totalTof = 0;
1592 int num = 0;
1593 // Carrying weight, error, and normalization
1594 double weight = 0;
1595 double errorSquared = 0;
1596 double normalization = 0.;
1597
1598 double bin_end = lastTof;
1599 std::function<bool(const double, const double)> compareTof;
1600 std::function<double(const double, double)> next_bin;
1601
1602 if (tolerance < 0) { // log
1603 if (lastTof < 0)
1604 throw std::runtime_error("compressEvents with log binning doesn't work with negative TOF");
1605
1606 if (lastTof == 0)
1607 bin_end = fabs(tolerance);
1608
1609 // for log we do "less than" so that is matches the log binning of the Rebin algorithm
1610 compareTof = [](const double lhs, const double rhs) { return lhs < rhs; };
1611 next_bin = [tolerance](const double lastTof, double bin_end) {
1612 // advance the bin_end until we find the one that this next event falls into
1613 while (lastTof >= bin_end)
1614 bin_end = bin_end * (1 - tolerance);
1615 return bin_end;
1616 };
1617 } else { // linear
1618 // for linear we do "less than or equals" because that is how it was originally implemented
1619 compareTof = [](const double lhs, const double rhs) { return lhs <= rhs; };
1620 next_bin = [tolerance](const double lastTof, double) { return lastTof + tolerance; };
1621 }
1622
1623 // get first bin_end
1624 bin_end = next_bin(lastTof, bin_end);
1625
1626 for (auto it = events.cbegin(); it != events.cend(); it++) {
1627 if (compareTof(it->m_tof, bin_end)) {
1628 // Carry the error and weight
1629 weight += it->weight();
1630 errorSquared += it->errorSquared();
1631 // Track the average tof
1632 num++;
1633 const double norm = calcNorm(it->errorSquared());
1634 normalization += norm;
1635 totalTof += it->m_tof * norm;
1636 } else {
1637 // We exceeded the tolerance
1638 // Create a new event with the average TOF and summed weights and
1639 // squared errors.
1640 if (num == 1) {
1641 // last time-of-flight is the only one contributing
1642 out.emplace_back(lastTof, weight, errorSquared);
1643 } else if (num > 1) {
1644 out.emplace_back(totalTof / normalization, weight, errorSquared);
1645 }
1646 // Start a new combined object
1647 num = 1;
1648 const double norm = calcNorm(it->errorSquared());
1649 normalization = norm;
1650 totalTof = it->m_tof * norm;
1651 weight = it->weight();
1652 errorSquared = it->errorSquared();
1653 lastTof = it->m_tof;
1654
1655 bin_end = next_bin(lastTof, bin_end);
1656 }
1657 }
1658
1659 // Put the last event in there too with the average TOF and summed weights and
1660 // squared errors.
1661 if (num == 1) {
1662 // last time-of-flight is the only one contributing
1663 out.emplace_back(lastTof, weight, errorSquared);
1664 } else if (num > 1) {
1665 out.emplace_back(totalTof / normalization, weight, errorSquared);
1666 }
1667
1668 // If you have over-allocated by more than 5%, reduce the size.
1669 size_t excess_limit = out.size() / 20;
1670 if ((out.capacity() - out.size()) > excess_limit) {
1671 out.shrink_to_fit();
1672 }
1673}
1674
1675template <class T>
1676inline void EventList::compressFatEventsHelper(const std::vector<T> &events, std::vector<WeightedEvent> &out,
1677 const double tolerance, const Types::Core::DateAndTime &timeStart,
1678 const double seconds) {
1679 // Clear the output. We can't know ahead of time how much space to reserve :(
1680 out.clear();
1681 // We will make a starting guess of 1/20th of the number of input events.
1682 out.reserve(events.size() / 20);
1683
1684 // The last TOF to which we are comparing.
1685 double lastTof = events.front().m_tof;
1686 // For getting an accurate average TOF
1687 double totalTof = 0;
1688
1689 // pulsetime bin information - stored as int nanoseconds because it
1690 // is the implementation type for DateAndTime object
1691 const int64_t pulsetimeStart = timeStart.totalNanoseconds();
1692 const auto pulsetimeDelta = static_cast<int64_t>(seconds * SEC_TO_NANO);
1693
1694 // pulsetime information
1695 std::vector<DateAndTime> pulsetimes; // all the times for new event
1696 std::vector<double> pulsetimeWeights;
1697
1698 // Carrying weight and error
1699 double weight = 0.;
1700 double errorSquared = 0.;
1701 double tofNormalization = 0.;
1702
1703 // Move up to first event that has a large enough pulsetime. This is just in case someone starts from after the
1704 // starttime of the run. It is expected that users will normally use the default which means this will only check the
1705 // first event.
1706 auto it = events.cbegin();
1707 for (; it != events.cend(); ++it) {
1708 if (it->m_pulsetime >= timeStart)
1709 break;
1710 }
1711
1712 if (it == events.cend())
1713 throw std::runtime_error("failed to find first pulse time in the events");
1714
1715 // bin if the pulses are histogrammed
1716 int64_t lastPulseBin = (it->m_pulsetime.totalNanoseconds() - pulsetimeStart) / pulsetimeDelta;
1717
1718 double bin_end = lastTof;
1719 double tof_min{0};
1720 std::function<bool(const double, const double)> compareTof;
1721 std::function<double(const double, double)> next_bin;
1722
1723 if (tolerance < 0) { // log
1724 // for log we do "less than" so that is matches the log binning of the Rebin algorithm
1725 compareTof = [](const double lhs, const double rhs) { return lhs < rhs; };
1726 next_bin = [tolerance](const double lastTof, double bin_end) {
1727 // advance the bin_end until we find the one that this next event falls into
1728 while (lastTof >= bin_end)
1729 bin_end = bin_end * (1 - tolerance);
1730 return bin_end;
1731 };
1732
1733 // get minimum Tof so that binning is consistent across all pulses
1734 const auto event_min = std::min_element(
1735 events.cbegin(), events.cend(), [](const auto &left, const auto &right) { return left.tof() < right.tof(); });
1736 bin_end = tof_min = event_min->tof();
1737
1738 if (tof_min < 0)
1739 throw std::runtime_error("compressEvents with log binning doesn't work with negative TOF");
1740
1741 // can't start at 0 as this will create an infinite loop
1742 if (tof_min == 0)
1743 bin_end = tof_min = fabs(tolerance);
1744
1745 } else { // linear
1746 // for linear we do "less than or equals" because that is how it was originally implemented
1747 compareTof = [](const double lhs, const double rhs) { return lhs <= rhs; };
1748 next_bin = [tolerance](const double lastTof, double) { return lastTof + tolerance; };
1749 }
1750
1751 // get first bin_end
1752 bin_end = next_bin(lastTof, bin_end);
1753
1754 // loop through events and accumulate weight
1755 for (; it != events.cend(); ++it) {
1756 const int64_t eventPulseBin = (it->m_pulsetime.totalNanoseconds() - pulsetimeStart) / pulsetimeDelta;
1757 if ((eventPulseBin <= lastPulseBin) && compareTof(it->m_tof, bin_end)) {
1758 // Carry the error and weight
1759 weight += it->weight();
1760 errorSquared += it->errorSquared();
1761 double norm = calcNorm(it->errorSquared());
1762 tofNormalization += norm;
1763 // Track the average tof
1764 totalTof += it->m_tof * norm;
1765 // Accumulate the pulse times
1766 pulsetimes.emplace_back(it->m_pulsetime);
1767 pulsetimeWeights.emplace_back(norm);
1768 } else {
1769 // We exceeded the tolerance
1770 if (!pulsetimes.empty()) {
1771 // Create a new event with the average TOF and summed weights and
1772 // squared errors. 1 event used doesn't need to average
1773 if (pulsetimes.size() == 1) {
1774 out.emplace_back(lastTof, pulsetimes.front(), weight, errorSquared);
1775 } else {
1776 out.emplace_back(totalTof / tofNormalization,
1777 Kernel::DateAndTimeHelpers::averageSorted(pulsetimes, pulsetimeWeights), weight,
1778 errorSquared);
1779 }
1780 }
1781 if (tolerance < 0 && eventPulseBin != lastPulseBin)
1782 // reset the bin_end for the new pulse bin
1783 bin_end = tof_min;
1784
1785 // Start a new combined object
1786 double norm = calcNorm(it->errorSquared());
1787 totalTof = it->m_tof * norm;
1788 weight = it->weight();
1789 errorSquared = it->errorSquared();
1790 tofNormalization = norm;
1791 lastTof = it->m_tof;
1792 lastPulseBin = eventPulseBin;
1793 pulsetimes.clear();
1794 pulsetimes.emplace_back(it->m_pulsetime);
1795 pulsetimeWeights.clear();
1796 pulsetimeWeights.emplace_back(norm);
1797
1798 bin_end = next_bin(lastTof, bin_end);
1799 }
1800 }
1801
1802 // Put the last event in there too.
1803 if (!pulsetimes.empty()) {
1804 // Create a new event with the average TOF and summed weights and
1805 // squared errors. 1 event used doesn't need to average
1806 if (pulsetimes.size() == 1) {
1807 out.emplace_back(lastTof, pulsetimes.front(), weight, errorSquared);
1808 } else {
1809 out.emplace_back(totalTof / tofNormalization,
1810 Kernel::DateAndTimeHelpers::averageSorted(pulsetimes, pulsetimeWeights), weight, errorSquared);
1811 }
1812 }
1813
1814 // If you have over-allocated by more than 5%, reduce the size.
1815 size_t excess_limit = out.size() / 20;
1816 if ((out.capacity() - out.size()) > excess_limit) {
1817 out.shrink_to_fit();
1818 }
1819}
1820
1821// --------------------------------------------------------------------------
1832 if (this->empty()) {
1833 // allocate memory in correct vector
1835 destination->weightedEventsNoTime = std::make_unique<std::vector<WeightedEventNoTime>>();
1836 } else {
1837 this->sortTof();
1838 switch (eventType) {
1839 case TOF:
1840 // if (parallel)
1841 // compressEventsParallelHelper(this->events,
1842 // destination->weightedEventsNoTime, tolerance);
1843 // else
1844 destination->weightedEventsNoTime = std::make_unique<std::vector<WeightedEventNoTime>>();
1846 break;
1847
1848 case WEIGHTED:
1849 // if (parallel)
1850 // compressEventsParallelHelper(this->weightedEvents,
1851 // destination->weightedEventsNoTime, tolerance);
1852 // else
1853 destination->weightedEventsNoTime = std::make_unique<std::vector<WeightedEventNoTime>>();
1855
1856 break;
1857
1858 case WEIGHTED_NOTIME:
1859 if (destination == this) {
1860 // Put results in a temp output
1861 auto out = std::make_unique<std::vector<WeightedEventNoTime>>();
1862 // if (parallel)
1863 // compressEventsParallelHelper(this->weightedEventsNoTime,
1864 // out,
1865 // tolerance);
1866 // else
1868 // Put it back
1869 this->weightedEventsNoTime.swap(out);
1870 } else {
1871 // if (parallel)
1872 // compressEventsParallelHelper(this->weightedEventsNoTime,
1873 // destination->weightedEventsNoTime, tolerance);
1874 // else
1875 destination->weightedEventsNoTime = std::make_unique<std::vector<WeightedEventNoTime>>();
1877 }
1878 break;
1879 }
1880 }
1881 // In all cases, you end up WEIGHTED_NOTIME.
1882 destination->eventType = WEIGHTED_NOTIME;
1883 // The sort is still valid!
1884 destination->order = TOF_SORT;
1885 // Empty out storage for vectors that are now unused.
1886 destination->clearUnused();
1887}
1888
1889template <class T>
1890inline void EventList::createWeightedEvents(std::vector<WeightedEventNoTime> &out, const std::vector<double> &tof,
1891 const std::vector<T> &weight, const std::vector<T> &error) {
1892 out.clear();
1893 for (size_t i = 0; i < weight.size(); ++i) {
1894 const auto errors = static_cast<float>(error[i]);
1895 if (errors > 0)
1896 out.emplace_back(tof[i], static_cast<float>(weight[i]), errors);
1897 }
1898}
1899
1900template <class T>
1901inline void EventList::processWeightedEvents(const std::vector<T> &events, std::vector<WeightedEventNoTime> &out,
1902 const std::shared_ptr<std::vector<double>> histogram_bin_edges,
1903 struct FindBin findBin) {
1904 const auto NUM_BINS = histogram_bin_edges->size() - 1;
1905 std::vector<double> tof(NUM_BINS, 0.);
1906 std::vector<double> normalization(NUM_BINS, 0.);
1907 std::vector<float> weight(NUM_BINS, 0.);
1908 std::vector<float> error(NUM_BINS, 0.);
1909 for (const auto &ev : events) {
1910 const auto &bin_optional = findBin(*histogram_bin_edges.get(), ev.m_tof, false);
1911 if (bin_optional) {
1912 const auto bin = bin_optional.value();
1913 const double norm = calcNorm(ev.m_errorSquared);
1914 tof[bin] += ev.m_tof * norm;
1915 normalization[bin] += norm;
1916 weight[bin] += ev.m_weight;
1917 error[bin] += ev.m_errorSquared;
1918 }
1919 }
1920
1921 // normalize TOFs
1922 std::transform(tof.begin(), tof.end(), normalization.begin(), tof.begin(), std::divides<double>());
1923
1924 createWeightedEvents(out, tof, weight, error);
1925}
1926
1928 const std::shared_ptr<std::vector<double>> histogram_bin_edges) {
1929 if (this->empty()) {
1930 // allocate memory in correct vector
1932 destination->weightedEventsNoTime = std::make_unique<std::vector<WeightedEventNoTime>>();
1933 } else {
1934 const auto NUM_BINS = histogram_bin_edges->size() - 1;
1935 const auto xmin = static_cast<double>(histogram_bin_edges->front());
1936
1937 auto findBin = FindBin(tolerance, xmin);
1938
1939 switch (eventType) {
1940 case TOF: {
1941 std::vector<double> tof(NUM_BINS, 0);
1942 std::vector<uint32_t> count(NUM_BINS, 0);
1943 for (const auto &ev : *this->events) {
1944 const auto &bin_optional = findBin(*histogram_bin_edges.get(), ev.m_tof, false);
1945 if (bin_optional) {
1946 const auto bin = bin_optional.value();
1947 count[bin]++;
1948 tof[bin] += ev.m_tof;
1949 }
1950 }
1951
1952 // average TOFs
1953 std::transform(tof.begin(), tof.end(), count.begin(), tof.begin(), std::divides<double>());
1954
1955 destination->weightedEventsNoTime = std::make_unique<std::vector<WeightedEventNoTime>>();
1957 break;
1958 }
1959
1960 case WEIGHTED: {
1961 destination->weightedEventsNoTime = std::make_unique<std::vector<WeightedEventNoTime>>();
1962 processWeightedEvents(*this->weightedEvents, *destination->weightedEventsNoTime, histogram_bin_edges, findBin);
1963 break;
1964 }
1965 case WEIGHTED_NOTIME:
1966 if (destination == this) {
1967 // Put results in a temp output
1968 auto out = std::make_unique<std::vector<WeightedEventNoTime>>();
1969 processWeightedEvents(*this->weightedEventsNoTime, *out, histogram_bin_edges, findBin);
1970 // Put it back
1971 this->weightedEventsNoTime.swap(out);
1972 } else {
1973 destination->weightedEventsNoTime = std::make_unique<std::vector<WeightedEventNoTime>>();
1974 processWeightedEvents(*this->weightedEventsNoTime, *destination->weightedEventsNoTime, histogram_bin_edges,
1975 findBin);
1976 }
1977 break;
1978 }
1979 }
1980
1981 // In all cases, you end up WEIGHTED_NOTIME.
1982 destination->eventType = WEIGHTED_NOTIME;
1983 // The result will be sorted
1984 destination->order = TOF_SORT;
1985 // Empty out storage for vectors that are now unused.
1986 destination->clearUnused();
1987}
1988
1989void EventList::compressFatEvents(const double tolerance, const Mantid::Types::Core::DateAndTime &timeStart,
1990 const double seconds, EventList *destination) {
1991 if (this->empty()) {
1992 // allocate memory in correct vector
1993 if (eventType != WEIGHTED)
1994 destination->weightedEvents = std::make_unique<std::vector<WeightedEvent>>();
1995 } else {
1996 switch (eventType) {
1997 case WEIGHTED_NOTIME:
1998 throw std::invalid_argument("Cannot compress events that do not have pulsetime");
1999 case TOF:
2000 this->sortPulseTimeTOFDelta(timeStart, seconds);
2001 destination->weightedEvents = std::make_unique<std::vector<WeightedEvent>>();
2002 compressFatEventsHelper(*this->events, *destination->weightedEvents, tolerance, timeStart, seconds);
2003 break;
2004 case WEIGHTED:
2005 this->sortPulseTimeTOFDelta(timeStart, seconds);
2006 if (destination == this) {
2007 // Put results in a temp output
2008 auto out = std::make_unique<std::vector<WeightedEvent>>();
2009 compressFatEventsHelper(*this->weightedEvents, *out, tolerance, timeStart, seconds);
2010 // Put it back
2011 this->weightedEvents.swap(out);
2012 } else {
2013 destination->weightedEvents = std::make_unique<std::vector<WeightedEvent>>();
2014 compressFatEventsHelper(*this->weightedEvents, *destination->weightedEvents, tolerance, timeStart, seconds);
2015 }
2016 break;
2017 }
2018 }
2019 // In all cases, you end up WEIGHTED_NOTIME.
2020 destination->eventType = WEIGHTED;
2021 // The sort order is pulsetimetof as we've compressed out the tolerance
2022 destination->order = PULSETIMETOF_SORT;
2023 // Empty out storage for vectors that are now unused.
2024 destination->clearUnused();
2025}
2026
2027// --------------------------------------------------------------------------
2037template <class T>
2038typename std::vector<T>::const_iterator static findFirstEvent(const std::vector<T> &events, T seek_tof) {
2039 return std::find_if_not(events.cbegin(), events.cend(), [seek_tof](const T &x) { return x < seek_tof; });
2040}
2041
2042// --------------------------------------------------------------------------
2052template <class T>
2053typename std::vector<T>::const_iterator EventList::findFirstPulseEvent(const std::vector<T> &events,
2054 const double seek_pulsetime) {
2055 auto itev = events.cbegin();
2056 auto itev_end = events.cend(); // cache for speed
2057
2058 // if tof < X[0], that means that you need to skip some events
2059 while ((itev != itev_end) && (static_cast<double>(itev->pulseTime().totalNanoseconds()) < seek_pulsetime))
2060 itev++;
2061 // Better fix would be to use a binary search instead of the linear one used
2062 // here.
2063 return itev;
2064}
2065
2066// --------------------------------------------------------------------------
2079template <class T>
2080typename std::vector<T>::const_iterator
2081EventList::findFirstTimeAtSampleEvent(const std::vector<T> &events, const double seek_time, const double &tofFactor,
2082 const double &tofOffset) const {
2083 auto itev = events.cbegin();
2084 auto itev_end = events.cend(); // cache for speed
2085
2086 // if tof < X[0], that means that you need to skip some events
2087 while ((itev != itev_end) &&
2088 (static_cast<double>(calculateCorrectedFullTime(*itev, tofFactor, tofOffset)) < seek_time))
2089 itev++;
2090 // Better fix would be to use a binary search instead of the linear one used
2091 // here.
2092 return itev;
2093}
2094
2095// --------------------------------------------------------------------------
2105template <class T> typename std::vector<T>::iterator static findFirstEvent(std::vector<T> &events, T seek_tof) {
2106 return std::find_if_not(events.begin(), events.end(), [seek_tof](const T &x) { return x < seek_tof; });
2107}
2108
2109// --------------------------------------------------------------------------
2119template <class T>
2120void EventList::histogramForWeightsHelper(const std::vector<T> &events, const MantidVec &X, MantidVec &Y,
2121 MantidVec &E) {
2122 // For slight speed=up.
2123 size_t x_size = X.size();
2124
2125 if (x_size <= 1) {
2126 // X was not set. Return an empty array.
2127 Y.resize(0, 0);
2128 return;
2129 }
2130
2131 // If the sizes are the same, then the "resize" command will NOT clear the original values.
2132 bool mustFill = (Y.size() == x_size - 1);
2133 // Clear the Y data, assign all to 0.
2134 Y.resize(x_size - 1, 0.0);
2135 // Clear the Error data, assign all to 0.
2136 // Note: Errors will be squared until the last step.
2137 E.resize(x_size - 1, 0.0);
2138
2139 if (mustFill) {
2140 // We must make sure the starting point is 0.0
2141 std::fill(Y.begin(), Y.end(), 0.0);
2142 std::fill(E.begin(), E.end(), 0.0);
2143 }
2144
2145 //---------------------- Histogram without weights
2146 //---------------------------------
2147
2148 // Do we even have any events to do?
2149 if (!events.empty()) {
2150 // Iterate through all events (sorted by tof)
2151 auto itev = findFirstEvent(events, T(X[0]));
2152 auto itev_end = events.cend();
2153 // The above can still take you to end() if no events above X[0], so check
2154 // again.
2155 if (itev == itev_end)
2156 return;
2157
2158 // Find the first bin
2159 size_t bin = 0;
2160 // The tof is greater the first bin boundary, so we need to find the first
2161 // bin
2162 double tof = itev->tof();
2163 while (bin < x_size - 1) {
2164 // Within range?
2165 if ((tof >= X[bin]) && (tof < X[bin + 1])) {
2166 // Add up the weight (convert to double before adding, to preserve
2167 // precision)
2168 Y[bin] += double(itev->m_weight);
2169 E[bin] += double(itev->m_errorSquared); // square of error
2170 break;
2171 }
2172 ++bin;
2173 }
2174 // Go to the next event, we've already binned this first one.
2175 ++itev;
2176
2177 // Keep going through all the events
2178 while ((itev != itev_end) && (bin < x_size - 1)) {
2179 tof = itev->tof();
2180 while (bin < x_size - 1) {
2181 // Within range? Since both events and X are sorted, they are going to
2182 // have
2183 // tof >= X[bin] because the previous event was.
2184 if (tof < X[bin + 1]) {
2185 // Add up the weight (convert to double before adding, to preserve
2186 // precision)
2187 Y[bin] += double(itev->m_weight);
2188 E[bin] += double(itev->m_errorSquared); // square of error
2189 break;
2190 }
2191 ++bin;
2192 }
2193 ++itev;
2194 }
2195 } // end if (there are any events to histogram)
2196
2197 // Now do the sqrt of all errors
2198 std::transform(E.cbegin(), E.cend(), E.begin(), static_cast<double (*)(double)>(sqrt));
2199}
2200
2201// --------------------------------------------------------------------------
2215template <class T>
2216void EventList::histogramForWeightsHelper(const std::vector<T> &events, const double step, const MantidVec &X,
2217 MantidVec &Y, MantidVec &E) {
2218 size_t x_size = X.size();
2219
2220 if (x_size <= 1) {
2221 // X was not set. Return an empty array.
2222 Y.resize(0, 0);
2223 return;
2224 }
2225
2226 // If the sizes are the same, then the "resize" command will NOT clear the original values.
2227 bool mustFill = (Y.size() == x_size - 1);
2228 Y.resize(x_size - 1, 0.0);
2229 E.resize(x_size - 1, 0.0);
2230 if (mustFill) {
2231 // We must make sure the starting point is 0.0
2232 std::fill(Y.begin(), Y.end(), 0.0);
2233 std::fill(E.begin(), E.end(), 0.0);
2234 }
2235
2236 if (events.empty())
2237 return;
2238
2239 const auto xmin = X.front();
2240 const auto xmax = X.back();
2241
2242 auto findBin = FindBin(step, xmin);
2243
2244 for (const T &ev : events) {
2245 const double tof = ev.tof();
2246 if (tof < xmin || tof >= xmax)
2247 continue;
2248
2249 std::optional<size_t> n_bin = findBin(X, tof, true);
2250
2251 if (n_bin) {
2252 Y[n_bin.value()] += ev.weight();
2253 E[n_bin.value()] += ev.errorSquared();
2254 }
2255 }
2256
2257 // Now do the sqrt of all errors
2258 std::transform(E.cbegin(), E.cend(), E.begin(), static_cast<double (*)(double)>(sqrt));
2259}
2260
2261// --------------------------------------------------------------------------
2271void EventList::generateHistogramPulseTime(const MantidVec &X, MantidVec &Y, MantidVec &E, bool skipError) const {
2272 // All types of weights need to be sorted by Pulse Time
2273 this->sortPulseTime();
2274
2275 switch (eventType) {
2276 case TOF:
2277 // Make the single ones
2279 if (!skipError)
2280 this->generateErrorsHistogram(Y, E);
2281 break;
2282
2283 case WEIGHTED:
2284 throw std::runtime_error("Cannot histogram by pulse time on Weighted "
2285 "Events currently"); // This could be supported.
2286
2287 case WEIGHTED_NOTIME:
2288 throw std::runtime_error("Cannot histogram by pulse time on Weighted Events NoTime");
2289 }
2290}
2291
2304 const double &tofOffset, bool skipError) const {
2305 // All types of weights need to be sorted by time at sample
2306 this->sortTimeAtSample(tofFactor, tofOffset);
2307
2308 switch (eventType) {
2309 case TOF:
2310 // Make the single ones
2311 this->generateCountsHistogramTimeAtSample(X, Y, tofFactor, tofOffset);
2312 if (!skipError)
2313 this->generateErrorsHistogram(Y, E);
2314 break;
2315
2316 case WEIGHTED:
2317 throw std::runtime_error("Cannot histogram by time at sample on Weighted "
2318 "Events currently"); // This could be supported.
2319
2320 case WEIGHTED_NOTIME:
2321 throw std::runtime_error("Cannot histogram by time at sample on Weighted Events NoTime");
2322 }
2323}
2324
2325// --------------------------------------------------------------------------
2335void EventList::generateHistogram(const MantidVec &X, MantidVec &Y, MantidVec &E, bool skipError) const {
2336 // All types of weights need to be sorted by TOF
2337
2338 this->sortTof();
2339
2340 switch (eventType) {
2341 case TOF:
2342 // Make the single ones
2343 this->generateCountsHistogram(X, Y);
2344 if (!skipError)
2345 this->generateErrorsHistogram(Y, E);
2346 break;
2347
2348 case WEIGHTED:
2350 break;
2351
2352 case WEIGHTED_NOTIME:
2354 break;
2355 }
2356}
2357
2358// --------------------------------------------------------------------------
2373void EventList::generateHistogram(const double step, const MantidVec &X, MantidVec &Y, MantidVec &E,
2374 bool skipError) const {
2375 // if events are already sorted, use faster sorted histogram method
2376 if (isSortedByTof() || empty())
2377 return generateHistogram(X, Y, E, skipError);
2378
2379 switch (eventType) {
2380 case TOF:
2381 this->generateCountsHistogram(step, X, Y);
2382 if (!skipError)
2383 this->generateErrorsHistogram(Y, E);
2384 break;
2385
2386 case WEIGHTED:
2387 histogramForWeightsHelper(*this->weightedEvents, step, X, Y, E);
2388 break;
2389
2390 case WEIGHTED_NOTIME:
2392 break;
2393 }
2394}
2395
2396// --------------------------------------------------------------------------
2404 // For slight speed=up.
2405 size_t x_size = X.size();
2406
2407 if (x_size <= 1) {
2408 // X was not set. Return an empty array.
2409 Y.resize(0, 0);
2410 return;
2411 }
2412
2413 // Sort the events by pulsetime
2414 this->sortPulseTime();
2415 // Clear the Y data, assign all to 0.
2416 Y.resize(x_size - 1, 0);
2417
2418 //---------------------- Histogram without weights
2419 //---------------------------------
2420
2421 if (!this->events->empty()) {
2422 // Iterate through all events (sorted by pulse time)
2423 auto itev = findFirstPulseEvent(*this->events, X[0]);
2424 auto itev_end = events->cend(); // cache for speed
2425 // The above can still take you to end() if no events above X[0], so check
2426 // again.
2427 if (itev == itev_end)
2428 return;
2429
2430 // Find the first bin
2431 size_t bin = 0;
2432
2433 // The tof is greater the first bin boundary, so we need to find the first
2434 // bin
2435 double pulsetime = static_cast<double>(itev->pulseTime().totalNanoseconds());
2436 while (bin < x_size - 1) {
2437 // Within range?
2438 if ((pulsetime >= X[bin]) && (pulsetime < X[bin + 1])) {
2439 Y[bin]++;
2440 break;
2441 }
2442 ++bin;
2443 }
2444 // Go to the next event, we've already binned this first one.
2445 ++itev;
2446
2447 // Keep going through all the events
2448 while ((itev != itev_end) && (bin < x_size - 1)) {
2449 pulsetime = static_cast<double>(itev->pulseTime().totalNanoseconds());
2450 while (bin < x_size - 1) {
2451 // Within range?
2452 if ((pulsetime >= X[bin]) && (pulsetime < X[bin + 1])) {
2453 Y[bin]++;
2454 break;
2455 }
2456 ++bin;
2457 }
2458 ++itev;
2459 }
2460 } // end if (there are any events to histogram)
2461}
2462
2477void EventList::generateCountsHistogramPulseTime(const double &xMin, const double &xMax, MantidVec &Y,
2478 const double TOF_min, const double TOF_max) const {
2479
2480 if (this->events->empty())
2481 return;
2482
2483 size_t nBins = Y.size();
2484
2485 if (nBins == 0)
2486 return;
2487
2488 double step = (xMax - xMin) / static_cast<double>(nBins);
2489
2490 for (const TofEvent &ev : *this->events) {
2491 double pulsetime = static_cast<double>(ev.pulseTime().totalNanoseconds());
2492 if (pulsetime < xMin || pulsetime >= xMax)
2493 continue;
2494 if (ev.tof() < TOF_min || ev.tof() >= TOF_max)
2495 continue;
2496
2497 auto n_bin = static_cast<size_t>((pulsetime - xMin) / step);
2498 Y[n_bin]++;
2499 }
2500}
2501
2502// --------------------------------------------------------------------------
2512 const double &tofOffset) const {
2513 // For slight speed=up.
2514 const size_t x_size = X.size();
2515
2516 if (x_size <= 1) {
2517 // X was not set. Return an empty array.
2518 Y.resize(0, 0);
2519 return;
2520 }
2521
2522 // Sort the events by pulsetime
2523 this->sortTimeAtSample(tofFactor, tofOffset);
2524 // Clear the Y data, assign all to 0.
2525 Y.resize(x_size - 1, 0);
2526
2527 //---------------------- Histogram without weights
2528 //---------------------------------
2529
2530 if (!this->events->empty()) {
2531 // Iterate through all events (sorted by pulse time)
2532 auto itev = findFirstTimeAtSampleEvent(*this->events, X[0], tofFactor, tofOffset);
2533 std::vector<TofEvent>::const_iterator itev_end = events->end(); // cache for speed
2534 // The above can still take you to end() if no events above X[0], so check
2535 // again.
2536 if (itev == itev_end)
2537 return;
2538
2539 // Find the first bin
2540 size_t bin = 0;
2541
2542 auto tAtSample = static_cast<double>(calculateCorrectedFullTime(*itev, tofFactor, tofOffset));
2543 while (bin < x_size - 1) {
2544 // Within range?
2545 if ((tAtSample >= X[bin]) && (tAtSample < X[bin + 1])) {
2546 Y[bin]++;
2547 break;
2548 }
2549 ++bin;
2550 }
2551 // Go to the next event, we've already binned this first one.
2552 ++itev;
2553
2554 // Keep going through all the events
2555 while ((itev != itev_end) && (bin < x_size - 1)) {
2556 tAtSample = static_cast<double>(calculateCorrectedFullTime(*itev, tofFactor, tofOffset));
2557 while (bin < x_size - 1) {
2558 // Within range?
2559 if ((tAtSample >= X[bin]) && (tAtSample < X[bin + 1])) {
2560 Y[bin]++;
2561 break;
2562 }
2563 ++bin;
2564 }
2565 ++itev;
2566 }
2567 } // end if (there are any events to histogram)
2568}
2569
2570// --------------------------------------------------------------------------
2577 // For slight speed=up.
2578 size_t x_size = X.size();
2579
2580 if (x_size <= 1) {
2581 // X was not set. Return an empty array.
2582 Y.resize(0, 0);
2583 return;
2584 }
2585
2586 // Sort the events by tof
2587 this->sortTof();
2588 // Clear the Y data, assign all to 0.
2589 Y.resize(x_size - 1, 0);
2590
2591 //---------------------- Histogram without weights
2592 //---------------------------------
2593
2594 // Do we even have any events to do?
2595 if (!this->events->empty()) {
2596 // Iterate through all events (sorted by tof) placing them in the correct
2597 // bin.
2598 auto itev = findFirstEvent(*this->events, TofEvent(X[0]));
2599 const auto itend = this->events->end();
2600 // Go through all the events,
2601 for (auto itx = X.cbegin(); itev != itend; ++itev) {
2602 const double tof = itev->tof();
2603 itx = std::find_if(itx, X.cend(), [tof](const double x) { return tof < x; });
2604 if (itx == X.cend()) {
2605 break;
2606 }
2607 const auto bin = static_cast<size_t>(std::max(std::distance(X.cbegin(), itx) - 1, std::ptrdiff_t{0}));
2608 ++Y[bin];
2609 }
2610 } // end if (there are any events to histogram)
2611}
2612
2623std::optional<size_t> EventList::findLinearBin(const MantidVec &X, const double tof, const double divisor,
2624 const double offset, const bool findExact) {
2625 const auto bin = static_cast<size_t>(tof * divisor - offset);
2626 if (bin >= X.size())
2627 return std::nullopt;
2628 else if (findExact)
2629 return findExactBin(X, tof, bin);
2630 else
2631 return bin;
2632}
2633
2652std::optional<size_t> EventList::findLogBin(const MantidVec &X, const double tof, const double divisor,
2653 const double offset, const bool findExact) {
2654 const auto bin = static_cast<size_t>(log(tof) * divisor - offset);
2655 if (bin >= X.size())
2656 return std::nullopt;
2657 else if (findExact)
2658 return findExactBin(X, tof, bin);
2659 else
2660 return bin;
2661}
2662
2670size_t EventList::findExactBin(const MantidVec &X, const double tof, const size_t n_bin) {
2671 // is tof slower than suggested bin
2672 auto tof_of_bin = X.cbegin() + n_bin; // boundary suggested
2673 if (tof < *tof_of_bin)
2674 return std::move(n_bin - 1);
2675
2676 // is tof higher than suggested bin
2677 ++tof_of_bin; // move to next boundary
2678 if (tof >= *tof_of_bin)
2679 return std::move(n_bin + 1);
2680
2681 // tof is in the bin
2682 return std::move(n_bin);
2683}
2684
2685// --------------------------------------------------------------------------
2696void EventList::generateCountsHistogram(const double step, const MantidVec &X, MantidVec &Y) const {
2697 // For slight speed=up.
2698 size_t x_size = X.size();
2699
2700 if (x_size <= 1) {
2701 // X was not set. Return an empty array.
2702 Y.resize(0, 0);
2703 return;
2704 }
2705
2706 // If the sizes are the same, then the "resize" command will NOT clear the original values.
2707 bool mustFill = (Y.size() == x_size - 1);
2708 // Clear the Y data, assign all to 0.
2709 Y.resize(x_size - 1, 0);
2710 if (mustFill) // starting point is no counts
2711 std::fill(Y.begin(), Y.end(), 0.0);
2712
2713 // Do we even have any events to do?
2714 if (this->events->empty())
2715 return;
2716
2717 const auto xmin = X.front();
2718 const auto xmax = X.back();
2719
2720 auto findBin = FindBin(step, xmin);
2721
2722 for (const TofEvent &ev : *this->events) {
2723 const double tof = ev.tof();
2724 if (tof < xmin || tof >= xmax)
2725 continue;
2726
2727 const std::optional<size_t> n_bin = findBin(X, tof, true);
2728
2729 if (n_bin)
2730 Y[n_bin.value()]++;
2731 }
2732}
2733
2734// --------------------------------------------------------------------------
2743 // Fill the vector for the errors, containing sqrt(count)
2744 E.resize(Y.size(), 0);
2745
2746 // windows can get confused about std::sqrt
2747 std::transform(Y.cbegin(), Y.cend(), E.begin(), static_cast<double (*)(double)>(sqrt));
2748
2749} //----------------------------------------------------------------------------------
2750
2761template <class T>
2762void EventList::integrateHelper(std::vector<T> &events, const double minX, const double maxX, const bool entireRange,
2763 double &sum, double &error) {
2764 sum = 0;
2765 error = 0;
2766 // Nothing in the list?
2767 if (events.empty())
2768 return;
2769
2770 // Iterators for limits - whole range by default
2771 auto lowit = events.cbegin();
2772 auto highit = events.cend();
2773
2774 // But maybe we don't want the entire range?
2775 if (!entireRange) {
2776 // If a silly range was given, return 0.
2777 if (maxX < minX)
2778 return;
2779
2780 // If the first element is lower that the xmin then search for new lowit
2781 if (lowit->tof() < minX)
2782 lowit = std::lower_bound(events.cbegin(), events.cend(), minX);
2783 // If the last element is higher that the xmax then search for new lowit
2784 if ((highit - 1)->tof() > maxX) {
2785 highit = std::upper_bound(lowit, events.cend(), T(maxX));
2786 }
2787 }
2788
2789 // Sum up all the weights
2790 for (auto it = lowit; it != highit; ++it) {
2791 sum += it->weight();
2792 error += it->errorSquared();
2793 }
2794 error = std::sqrt(error);
2795}
2796
2797// --------------------------------------------------------------------------
2806double EventList::integrate(const double minX, const double maxX, const bool entireRange) const {
2807 double sum(0), error(0);
2808 integrate(minX, maxX, entireRange, sum, error);
2809 return sum;
2810}
2811
2821void EventList::integrate(const double minX, const double maxX, const bool entireRange, double &sum,
2822 double &error) const {
2823 sum = 0;
2824 error = 0;
2825 if (!entireRange) {
2826 // The event list must be sorted by TOF!
2827 this->sortTof();
2828 }
2829
2830 // Convert the list
2831 switch (eventType) {
2832 case TOF:
2833 integrateHelper(*this->events, minX, maxX, entireRange, sum, error);
2834 break;
2835 case WEIGHTED:
2836 integrateHelper(*this->weightedEvents, minX, maxX, entireRange, sum, error);
2837 break;
2838 case WEIGHTED_NOTIME:
2839 integrateHelper(*this->weightedEventsNoTime, minX, maxX, entireRange, sum, error);
2840 break;
2841 default:
2842 throw std::runtime_error("EventList: invalid event type value was found.");
2843 }
2844}
2845
2846// ==============================================================================================
2847// ----------- Conversion Functions (changing tof values)
2848// ---------------------------------------
2849// ==============================================================================================
2850
2857void EventList::convertTof(std::function<double(double)> func, const int sorting) {
2858 // fix the histogram parameter
2859 MantidVec &x = dataX();
2860 transform(x.cbegin(), x.cend(), x.begin(), func);
2861
2862 // do nothing if sorting > 0
2863 if (sorting == 0) {
2864 this->setSortOrder(UNSORTED);
2865 } else if ((sorting < 0) && (this->getSortType() == TOF_SORT)) {
2866 this->reverse();
2867 }
2868
2869 if (this->getNumberEvents() == 0)
2870 return;
2871
2872 // Convert the list
2873 switch (eventType) {
2874 case TOF:
2875 this->convertTofHelper(*this->events, func);
2876 break;
2877 case WEIGHTED:
2878 this->convertTofHelper(*this->weightedEvents, func);
2879 break;
2880 case WEIGHTED_NOTIME:
2881 this->convertTofHelper(*this->weightedEventsNoTime, func);
2882 break;
2883 }
2884}
2885
2890template <class T> void EventList::convertTofHelper(std::vector<T> &events, const std::function<double(double)> &func) {
2891 // iterate through all events
2892 for (auto &ev : events)
2893 ev.m_tof = func(ev.m_tof);
2894}
2895
2896// --------------------------------------------------------------------------
2902void EventList::convertTof(const double factor, const double offset) {
2903 // fix the histogram parameter
2904 auto &x = mutableX();
2905 x *= factor;
2906 x += offset;
2907
2908 if ((factor < 0.) && (this->getSortType() == TOF_SORT))
2909 this->reverse();
2910
2911 if (this->getNumberEvents() == 0)
2912 return;
2913
2914 // Convert the list
2915 switch (eventType) {
2916 case TOF:
2917 this->convertTofHelper(*this->events, factor, offset);
2918 break;
2919 case WEIGHTED:
2920 this->convertTofHelper(*this->weightedEvents, factor, offset);
2921 break;
2922 case WEIGHTED_NOTIME:
2923 this->convertTofHelper(*this->weightedEventsNoTime, factor, offset);
2924 break;
2925 }
2926}
2927
2928// --------------------------------------------------------------------------
2937template <class T> void EventList::convertTofHelper(std::vector<T> &events, const double factor, const double offset) {
2938 // iterate through all events
2939 for (auto &event : events) {
2940 event.m_tof = event.m_tof * factor + offset;
2941 }
2942}
2943
2944// --------------------------------------------------------------------------
2951void EventList::scaleTof(const double factor) { this->convertTof(factor, 0.0); }
2952
2953// --------------------------------------------------------------------------
2958void EventList::addTof(const double offset) { this->convertTof(1.0, offset); }
2959
2960// --------------------------------------------------------------------------
2966template <class T> void EventList::addPulsetimeHelper(std::vector<T> &events, const double seconds) {
2967 // iterate through all events
2968 for (auto &event : events) {
2969 event.m_pulsetime += seconds;
2970 }
2971}
2972
2979template <class T> void EventList::addPulsetimesHelper(std::vector<T> &events, const std::vector<double> &seconds) {
2980 auto eventIterEnd{events.end()};
2981 auto secondsIter{seconds.cbegin()};
2982 for (auto eventIter = events.begin(); eventIter < eventIterEnd; ++eventIter, ++secondsIter) {
2983 eventIter->m_pulsetime += *secondsIter;
2984 }
2985}
2986
2987// --------------------------------------------------------------------------
2992void EventList::addPulsetime(const double seconds) {
2993 if (this->getNumberEvents() == 0)
2994 return;
2995
2996 // Convert the list
2997 switch (eventType) {
2998 case TOF:
2999 this->addPulsetimeHelper(*this->events, seconds);
3000 break;
3001 case WEIGHTED:
3002 this->addPulsetimeHelper(*this->weightedEvents, seconds);
3003 break;
3004 case WEIGHTED_NOTIME:
3005 throw std::runtime_error("EventList::addPulsetime() called on an event "
3006 "list with no pulse times. You must call this "
3007 "algorithm BEFORE CompressEvents.");
3008 break;
3009 }
3010}
3011
3012// --------------------------------------------------------------------------
3017void EventList::addPulsetimes(const std::vector<double> &seconds) {
3018 if (this->getNumberEvents() == 0)
3019 return;
3020 if (this->getNumberEvents() != seconds.size()) {
3021 throw std::runtime_error("");
3022 }
3023
3024 // Convert the list
3025 switch (eventType) {
3026 case TOF:
3027 this->addPulsetimesHelper(*this->events, seconds);
3028 break;
3029 case WEIGHTED:
3030 this->addPulsetimesHelper(*this->weightedEvents, seconds);
3031 break;
3032 case WEIGHTED_NOTIME:
3033 throw std::runtime_error("EventList::addPulsetime() called on an event "
3034 "list with no pulse times. You must call this "
3035 "algorithm BEFORE CompressEvents.");
3036 break;
3037 }
3038}
3039
3040// --------------------------------------------------------------------------
3048template <class T>
3049std::size_t EventList::maskTofHelper(std::vector<T> &events, const double tofMin, const double tofMax) {
3050 // quick checks to make sure that the masking range is even in the data
3051 if (tofMin > events.crbegin()->tof())
3052 return 0;
3053 if (tofMax < events.cbegin()->tof())
3054 return 0;
3055
3056 // Find the index of the first tofMin
3057 auto it_first = std::lower_bound(events.begin(), events.end(), tofMin);
3058 if ((it_first != events.end()) && (it_first->tof() < tofMax)) {
3059 // Something was found
3060 // Look for the first one > tofMax
3061 auto it_last = std::upper_bound(it_first, events.end(), T(tofMax));
3062
3063 if (it_first >= it_last) {
3064 throw std::runtime_error("Event filter is all messed up"); // TODO
3065 }
3066
3067 size_t tmp = std::size_t(std::distance(it_first, it_last));
3068 // it_last will either be at the end (if not found) or before it.
3069 // Erase this range from the vector
3070 events.erase(it_first, it_last);
3071
3072 // Done! Sorting is still valid, no need to redo.
3073 return tmp; //(it_last - it_first); the iterators get invalid after erase
3074 }
3075 return 0; // didn't remove any events
3076}
3077
3078// --------------------------------------------------------------------------
3085void EventList::maskTof(const double tofMin, const double tofMax) {
3086 if (tofMax <= tofMin)
3087 throw std::runtime_error("EventList::maskTof: tofMax must be > tofMin");
3088
3089 // don't do anything with an emply list
3090 if (this->getNumberEvents() == 0)
3091 return;
3092
3093 // Start by sorting by tof
3094 this->sortTof();
3095
3096 // Convert the list
3097 size_t numOrig = 0;
3098 size_t numDel = 0;
3099 switch (eventType) {
3100 case TOF:
3101 numOrig = this->events->size();
3102 numDel = this->maskTofHelper(*this->events, tofMin, tofMax);
3103 break;
3104 case WEIGHTED:
3105 numOrig = this->weightedEvents->size();
3106 numDel = this->maskTofHelper(*this->weightedEvents, tofMin, tofMax);
3107 break;
3108 case WEIGHTED_NOTIME:
3109 numOrig = this->weightedEventsNoTime->size();
3110 numDel = this->maskTofHelper(*this->weightedEventsNoTime, tofMin, tofMax);
3111 break;
3112 }
3113
3114 if (numDel >= numOrig)
3115 this->clear(false);
3116}
3117
3118// --------------------------------------------------------------------------
3125template <class T> std::size_t EventList::maskConditionHelper(std::vector<T> &events, const std::vector<bool> &mask) {
3126
3127 // runs through the two synchronized vectors and delete elements
3128 // for condition false
3129 auto itm = std::find(mask.begin(), mask.end(), false);
3130 auto first = events.begin() + (itm - mask.begin());
3131
3132 if (itm != mask.end()) {
3133 for (auto ite = first; ++ite != events.end() && ++itm != mask.end();) {
3134 if (*itm != false) {
3135 *first++ = std::move(*ite);
3136 }
3137 }
3138 }
3139
3140 const auto n = static_cast<size_t>(events.end() - first);
3141 if (n != 0)
3142 events.erase(first, events.end());
3143
3144 return n;
3145}
3146
3147// --------------------------------------------------------------------------
3153void EventList::maskCondition(const std::vector<bool> &mask) {
3154
3155 // mask size must match the number of events
3156 if (this->getNumberEvents() != mask.size())
3157 throw std::runtime_error("EventList::maskTof: tofMax must be > tofMin");
3158
3159 // don't do anything with an emply list
3160 if (this->getNumberEvents() == 0)
3161 return;
3162
3163 // Convert the list
3164 size_t numOrig = 0;
3165 size_t numDel = 0;
3166 switch (eventType) {
3167 case TOF:
3168 numOrig = this->events->size();
3169 numDel = this->maskConditionHelper(*this->events, mask);
3170 break;
3171 case WEIGHTED:
3172 numOrig = this->weightedEvents->size();
3173 numDel = this->maskConditionHelper(*this->weightedEvents, mask);
3174 break;
3175 case WEIGHTED_NOTIME:
3176 numOrig = this->weightedEventsNoTime->size();
3177 numDel = this->maskConditionHelper(*this->weightedEventsNoTime, mask);
3178 break;
3179 }
3180
3181 if (numDel >= numOrig)
3182 this->clear(false);
3183}
3184
3185// --------------------------------------------------------------------------
3191template <class T> void EventList::getTofsHelper(const std::vector<T> &events, std::vector<double> &tofs) {
3192 tofs.clear();
3193 for (auto itev = events.cbegin(); itev != events.cend(); ++itev)
3194 tofs.emplace_back(itev->m_tof);
3195}
3196
3200void EventList::getTofs(std::vector<double> &tofs) const {
3201 // Set the capacity of the vector to avoid multiple resizes
3202 tofs.reserve(this->getNumberEvents());
3203
3204 // Convert the list
3205 switch (eventType) {
3206 case TOF:
3207 this->getTofsHelper(*this->events, tofs);
3208 break;
3209 case WEIGHTED:
3210 this->getTofsHelper(*this->weightedEvents, tofs);
3211 break;
3212 case WEIGHTED_NOTIME:
3213 this->getTofsHelper(*this->weightedEventsNoTime, tofs);
3214 break;
3215 }
3216}
3217
3222std::vector<double> EventList::getTofs() const {
3223 std::vector<double> tofs;
3224 this->getTofs(tofs);
3225 return tofs;
3226}
3227
3228// --------------------------------------------------------------------------
3234template <class T> void EventList::getWeightsHelper(const std::vector<T> &events, std::vector<double> &weights) {
3235 weights.clear();
3236 weights.reserve(events.size());
3237 std::transform(events.cbegin(), events.cend(), std::back_inserter(weights),
3238 [](const auto &event) { return event.weight(); });
3239}
3240
3244void EventList::getWeights(std::vector<double> &weights) const {
3245 // Set the capacity of the vector to avoid multiple resizes
3246 weights.reserve(this->getNumberEvents());
3247
3248 // Convert the list
3249 switch (eventType) {
3250 case WEIGHTED:
3251 this->getWeightsHelper(*this->weightedEvents, weights);
3252 break;
3253 case WEIGHTED_NOTIME:
3254 this->getWeightsHelper(*this->weightedEventsNoTime, weights);
3255 break;
3256 default:
3257 // not a weighted event type, return 1.0 for all.
3258 weights.assign(this->getNumberEvents(), 1.0);
3259 break;
3260 }
3261}
3262
3267std::vector<double> EventList::getWeights() const {
3268 std::vector<double> weights;
3269 this->getWeights(weights);
3270 return weights;
3271}
3272
3273// --------------------------------------------------------------------------
3279template <class T>
3280void EventList::getWeightErrorsHelper(const std::vector<T> &events, std::vector<double> &weightErrors) {
3281 weightErrors.clear();
3282 weightErrors.reserve(events.size());
3283 std::transform(events.cbegin(), events.cend(), std::back_inserter(weightErrors),
3284 [](const auto &event) { return event.error(); });
3285}
3286
3290void EventList::getWeightErrors(std::vector<double> &weightErrors) const {
3291 // Set the capacity of the vector to avoid multiple resizes
3292 weightErrors.reserve(this->getNumberEvents());
3293
3294 // Convert the list
3295 switch (eventType) {
3296 case WEIGHTED:
3297 this->getWeightErrorsHelper(*this->weightedEvents, weightErrors);
3298 break;
3299 case WEIGHTED_NOTIME:
3300 this->getWeightErrorsHelper(*this->weightedEventsNoTime, weightErrors);
3301 break;
3302 default:
3303 // not a weighted event type, return 1.0 for all.
3304 weightErrors.assign(this->getNumberEvents(), 1.0);
3305 break;
3306 }
3307}
3308
3313std::vector<double> EventList::getWeightErrors() const {
3314 std::vector<double> weightErrors;
3315 this->getWeightErrors(weightErrors);
3316 return weightErrors;
3317}
3318
3324template <typename UnaryOperation>
3325std::vector<DateAndTime> EventList::eventTimesCalculator(const UnaryOperation &timesCalc) const {
3326 std::vector<DateAndTime> times;
3327 switch (eventType) {
3328 case TOF:
3329 times.reserve(events->size());
3330 std::transform(events->cbegin(), events->cend(), std::back_inserter(times), timesCalc);
3331 break;
3332 case WEIGHTED:
3333 times.reserve(weightedEvents->size());
3334 std::transform(weightedEvents->cbegin(), weightedEvents->cend(), std::back_inserter(times), timesCalc);
3335 break;
3336 case WEIGHTED_NOTIME:
3337 times.reserve(weightedEventsNoTime->size());
3338 std::transform(weightedEventsNoTime->cbegin(), weightedEventsNoTime->cend(), std::back_inserter(times), timesCalc);
3339 break;
3340 }
3341 return times;
3342}
3343
3348std::vector<Mantid::Types::Core::DateAndTime> EventList::getPulseTimes() const {
3349 auto timeCalc = [](const auto &event) { return event.pulseTime(); };
3350 return eventTimesCalculator(timeCalc);
3351}
3352
3354std::vector<DateAndTime> EventList::getPulseTOFTimes() const {
3355 auto timeCalc = [](const auto &event) { return event.pulseTOFTime(); };
3356 return eventTimesCalculator(timeCalc);
3357}
3358
3363std::vector<DateAndTime> EventList::getPulseTOFTimesAtSample(const double &factor, const double &shift) const {
3364 auto timeCalc = [factor, shift](const auto &event) { return event.pulseTOFTimeAtSample(factor, shift); };
3365 return eventTimesCalculator(timeCalc);
3366}
3367
3368// --------------------------------------------------------------------------
3369
3370namespace { // anonymous namespace
3371template <class T> double getTofMinimumHelper(const std::vector<T> &events) {
3372 const auto result = std::min_element(events.cbegin(), events.cend(),
3373 [](const auto &left, const auto &right) { return left.tof() < right.tof(); });
3374 return result->tof();
3375}
3376
3377template <class T> double getTofMaximumHelper(const std::vector<T> &events) {
3378 const auto result = std::max_element(events.cbegin(), events.cend(),
3379 [](const auto &left, const auto &right) { return left.tof() < right.tof(); });
3380 return result->tof();
3381}
3382} // anonymous namespace
3383
3387double EventList::getTofMin() const {
3388 // set up as the maximum available double
3389 double tMin = std::numeric_limits<double>::max();
3390
3391 // no events is a soft error
3392 if (this->empty())
3393 return tMin;
3394
3395 // when events are ordered by tof just need the first value
3396 if (this->order == TOF_SORT) {
3397 switch (eventType) {
3398 case TOF:
3399 return this->events->front().tof();
3400 case WEIGHTED:
3401 return this->weightedEvents->front().tof();
3402 case WEIGHTED_NOTIME:
3403 return this->weightedEventsNoTime->front().tof();
3404 }
3405 }
3406
3407 // now we are stuck with a linear search
3408 switch (eventType) {
3409 case TOF: {
3410 tMin = getTofMinimumHelper(*this->events);
3411 break;
3412 }
3413 case WEIGHTED: {
3414 tMin = getTofMinimumHelper(*this->weightedEvents);
3415 break;
3416 }
3417 case WEIGHTED_NOTIME: {
3418 tMin = getTofMinimumHelper(*this->weightedEventsNoTime);
3419 break;
3420 }
3421 }
3422
3423 return tMin;
3424}
3425
3429double EventList::getTofMax() const {
3430 // set up as the minimum available double
3431 double tMax = std::numeric_limits<double>::lowest();
3432
3433 // no events is a soft error
3434 if (this->empty())
3435 return tMax;
3436
3437 // when events are ordered by tof just need the first value
3438 if (this->order == TOF_SORT) {
3439 switch (eventType) {
3440 case TOF:
3441 return this->events->back().tof();
3442 case WEIGHTED:
3443 return this->weightedEvents->back().tof();
3444 case WEIGHTED_NOTIME:
3445 return this->weightedEventsNoTime->back().tof();
3446 }
3447 }
3448
3449 // now we are stuck with a linear search
3450 switch (eventType) {
3451 case TOF: {
3452 tMax = getTofMaximumHelper(*this->events);
3453 break;
3454 }
3455 case WEIGHTED: {
3456 tMax = getTofMaximumHelper(*this->weightedEvents);
3457 break;
3458 }
3459 case WEIGHTED_NOTIME: {
3460 tMax = getTofMaximumHelper(*this->weightedEventsNoTime);
3461 break;
3462 }
3463 }
3464
3465 return tMax;
3466}
3467
3468// --------------------------------------------------------------------------
3469namespace { // anonymous namespace
3470template <class T> DateAndTime getPulseMinimumHelper(const std::vector<T> &events) {
3471 const auto result = std::min_element(events.cbegin(), events.cend(), [](const auto &left, const auto &right) {
3472 return left.pulseTime() < right.pulseTime();
3473 });
3474 return result->pulseTime();
3475}
3476
3477template <class T> DateAndTime getPulseMaximumHelper(const std::vector<T> &events) {
3478 const auto result = std::max_element(events.cbegin(), events.cend(), [](const auto &left, const auto &right) {
3479 return left.pulseTime() < right.pulseTime();
3480 });
3481 return result->pulseTime();
3482}
3483} // anonymous namespace
3484
3488DateAndTime EventList::getPulseTimeMin() const {
3489 // no events is a soft error
3490 if (this->empty())
3491 return DateAndTime::maximum();
3492
3493 // when events are ordered by pulse time just need the first value
3494 if (this->order == PULSETIME_SORT) {
3495 switch (eventType) {
3496 case TOF:
3497 return this->events->front().pulseTime();
3498 case WEIGHTED:
3499 return this->weightedEvents->front().pulseTime();
3500 case WEIGHTED_NOTIME:
3501 return this->weightedEventsNoTime->front().pulseTime();
3502 }
3503 }
3504
3505 // now we are stuck with a linear search
3506 switch (eventType) {
3507 case TOF:
3508 return getPulseMinimumHelper(*this->events);
3509 case WEIGHTED:
3510 return getPulseMinimumHelper(*this->weightedEvents);
3511 case WEIGHTED_NOTIME:
3512 return getPulseMinimumHelper(*this->weightedEventsNoTime);
3513 }
3514
3515 return DateAndTime::maximum();
3516}
3517
3521DateAndTime EventList::getPulseTimeMax() const {
3522 // no events is a soft error
3523 if (this->empty())
3524 return DateAndTime::minimum();
3525
3526 // when events are ordered by pulse time just need the first value
3527 if (this->order == PULSETIME_SORT) {
3528 switch (eventType) {
3529 case TOF:
3530 return this->events->back().pulseTime();
3531 case WEIGHTED:
3532 return this->weightedEvents->back().pulseTime();
3533 case WEIGHTED_NOTIME:
3534 return this->weightedEventsNoTime->back().pulseTime();
3535 }
3536 }
3537
3538 // now we are stuck with a linear search
3539 switch (eventType) {
3540 case TOF:
3541 return getPulseMaximumHelper(*this->events);
3542 case WEIGHTED:
3543 return getPulseMaximumHelper(*this->weightedEvents);
3544 case WEIGHTED_NOTIME:
3545 return getPulseMaximumHelper(*this->weightedEventsNoTime);
3546 }
3547
3548 return DateAndTime::minimum();
3549}
3550
3551void EventList::getPulseTimeMinMax(Mantid::Types::Core::DateAndTime &tMin,
3552 Mantid::Types::Core::DateAndTime &tMax) const {
3553 // set up as the minimum available date time.
3554 tMax = DateAndTime::minimum();
3555 tMin = DateAndTime::maximum();
3556
3557 // no events is a soft error
3558 if (this->empty())
3559 return;
3560
3561 // when events are ordered by pulse time just need the first/last values
3562 if (this->order == PULSETIME_SORT) {
3563 switch (eventType) {
3564 case TOF:
3565 tMin = this->events->front().pulseTime();
3566 tMax = this->events->back().pulseTime();
3567 return;
3568 case WEIGHTED:
3569 tMin = this->weightedEvents->front().pulseTime();
3570 tMax = this->weightedEvents->back().pulseTime();
3571 return;
3572 case WEIGHTED_NOTIME:
3573 tMin = this->weightedEventsNoTime->front().pulseTime();
3574 tMax = this->weightedEventsNoTime->back().pulseTime();
3575 return;
3576 }
3577 }
3578
3579 // now we are stuck with a linear search
3580 // could this be done more efficiently than using ->at?
3581 size_t numEvents = this->getNumberEvents();
3582 DateAndTime temp = tMax; // start with the smallest possible value
3583 for (size_t i = 0; i < numEvents; i++) {
3584 switch (eventType) {
3585 case TOF:
3586 temp = this->events->at(i).pulseTime();
3587 break;
3588 case WEIGHTED:
3589 temp = this->weightedEvents->at(i).pulseTime();
3590 break;
3591 case WEIGHTED_NOTIME:
3592 temp = this->weightedEventsNoTime->at(i).pulseTime();
3593 break;
3594 }
3595 if (temp > tMax)
3596 tMax = temp;
3597 if (temp < tMin)
3598 tMin = temp;
3599 }
3600}
3601
3602DateAndTime EventList::getTimeAtSampleMax(const double &tofFactor, const double &tofOffset) const {
3603 // set up as the minimum available date time.
3604 DateAndTime tMax = DateAndTime::minimum();
3605
3606 // no events is a soft error
3607 if (this->empty())
3608 return tMax;
3609
3610 // when events are ordered by time at sample just need the first value
3611 if (this->order == TIMEATSAMPLE_SORT) {
3612 switch (eventType) {
3613 case TOF:
3614 return calculateCorrectedFullTime(this->events->back(), tofFactor, tofOffset);
3615 case WEIGHTED:
3616 return calculateCorrectedFullTime(this->weightedEvents->back(), tofFactor, tofOffset);
3617 case WEIGHTED_NOTIME:
3618 return calculateCorrectedFullTime(this->weightedEventsNoTime->back(), tofFactor, tofOffset);
3619 }
3620 }
3621
3622 // now we are stuck with a linear search
3623 size_t numEvents = this->getNumberEvents();
3624 DateAndTime temp = tMax; // start with the smallest possible value
3625 for (size_t i = 0; i < numEvents; i++) {
3626 switch (eventType) {
3627 case TOF:
3628 temp = calculateCorrectedFullTime(this->events->at(i), tofFactor, tofOffset);
3629 break;
3630 case WEIGHTED:
3631 temp = calculateCorrectedFullTime(this->weightedEvents->at(i), tofFactor, tofOffset);
3632 break;
3633 case WEIGHTED_NOTIME:
3634 temp = calculateCorrectedFullTime(this->weightedEventsNoTime->at(i), tofFactor, tofOffset);
3635 break;
3636 }
3637 if (temp > tMax)
3638 tMax = temp;
3639 }
3640 return tMax;
3641}
3642
3643DateAndTime EventList::getTimeAtSampleMin(const double &tofFactor, const double &tofOffset) const {
3644 // set up as the minimum available date time.
3645 DateAndTime tMin = DateAndTime::maximum();
3646
3647 // no events is a soft error
3648 if (this->empty())
3649 return tMin;
3650
3651 // when events are ordered by time at sample just need the first value
3652 if (this->order == TIMEATSAMPLE_SORT) {
3653 switch (eventType) {
3654 case TOF:
3655 return calculateCorrectedFullTime(this->events->front(), tofFactor, tofOffset);
3656 case WEIGHTED:
3657 return calculateCorrectedFullTime(this->weightedEvents->front(), tofFactor, tofOffset);
3658 case WEIGHTED_NOTIME:
3659 return calculateCorrectedFullTime(this->weightedEventsNoTime->front(), tofFactor, tofOffset);
3660 }
3661 }
3662
3663 // now we are stuck with a linear search
3664 size_t numEvents = this->getNumberEvents();
3665 DateAndTime temp = tMin; // start with the smallest possible value
3666 for (size_t i = 0; i < numEvents; i++) {
3667 switch (eventType) {
3668 case TOF:
3669 temp = calculateCorrectedFullTime(this->events->at(i), tofFactor, tofOffset);
3670 break;
3671 case WEIGHTED:
3672 temp = calculateCorrectedFullTime(this->weightedEvents->at(i), tofFactor, tofOffset);
3673 break;
3674 case WEIGHTED_NOTIME:
3675 temp = calculateCorrectedFullTime(this->weightedEventsNoTime->at(i), tofFactor, tofOffset);
3676 break;
3677 }
3678 if (temp < tMin)
3679 tMin = temp;
3680 }
3681 return tMin;
3682}
3683
3684// --------------------------------------------------------------------------
3690template <class T> void EventList::setTofsHelper(std::vector<T> &events, const std::vector<double> &tofs) {
3691 if (tofs.empty())
3692 return;
3693
3694 size_t x_size = tofs.size();
3695 if (events.size() != x_size)
3696 return; // should this throw an exception?
3697
3698 for (size_t i = 0; i < x_size; ++i)
3699 events[i].m_tof = tofs[i];
3700}
3701
3702// --------------------------------------------------------------------------
3709 this->order = UNSORTED;
3710
3711 // Convert the list
3712 switch (eventType) {
3713 case TOF:
3714 this->setTofsHelper(*this->events, tofs);
3715 break;
3716 case WEIGHTED:
3717 this->setTofsHelper(*this->weightedEvents, tofs);
3718 break;
3719 case WEIGHTED_NOTIME:
3720 this->setTofsHelper(*this->weightedEventsNoTime, tofs);
3721 break;
3722 }
3723}
3724
3725// ==============================================================================================
3726// ----------- MULTIPLY AND DIVIDE ---------------------------------------
3727// ==============================================================================================
3728
3729//------------------------------------------------------------------------------------------------
3737template <class T> void EventList::multiplyHelper(std::vector<T> &events, const double value, const double error) {
3738 // Square of the value
3739 const double valueSquared = value * value;
3740
3741 auto itev_end = events.end();
3742
3743 if (error == 0) {
3744 // Error-less calculation
3745 for (auto itev = events.begin(); itev != itev_end; itev++) {
3746 itev->m_errorSquared = static_cast<float>(itev->m_errorSquared * valueSquared);
3747 itev->m_weight *= static_cast<float>(value);
3748 }
3749 } else {
3750 // Carry the scalar error
3751 const double errorSquared = error * error; // Square of the value's error
3752 for (auto itev = events.begin(); itev != itev_end; itev++) {
3753 itev->m_errorSquared =
3754 static_cast<float>(itev->m_errorSquared * valueSquared + errorSquared * itev->m_weight * itev->m_weight);
3755 itev->m_weight *= static_cast<float>(value);
3756 }
3757 }
3758}
3759
3760//------------------------------------------------------------------------------------------------
3773 this->multiply(value);
3774 return *this;
3775}
3776
3777//------------------------------------------------------------------------------------------------
3806void EventList::multiply(const double value, const double error) {
3807 // Do nothing if multiplying by exactly one and there is no error
3808 if ((value == 1.0) && (error == 0.0))
3809 return;
3810
3811 switch (eventType) {
3812 case TOF:
3813 // Switch to weights if needed.
3814 this->switchTo(WEIGHTED);
3815 // Fall through
3816
3817 case WEIGHTED:
3818 multiplyHelper(*this->weightedEvents, value, error);
3819 break;
3820
3821 case WEIGHTED_NOTIME:
3823 break;
3824 }
3825}
3826
3827//------------------------------------------------------------------------------------------------
3836template <class T>
3837void EventList::multiplyHistogramHelper(std::vector<T> &events, const MantidVec &X, const MantidVec &Y,
3838 const MantidVec &E) {
3839 // Validate inputs
3840 if ((X.size() < 2) || (Y.size() != E.size()) || (X.size() != 1 + Y.size())) {
3841 std::stringstream msg;
3842 msg << "EventList::multiply() was given invalid size or "
3843 "inconsistent histogram arrays: X["
3844 << X.size() << "] "
3845 << "Y[" << Y.size() << " E[" << E.size() << "]";
3846 throw std::invalid_argument(msg.str());
3847 }
3848
3849 size_t x_size = X.size();
3850
3851 // Iterate through all events (sorted by tof)
3852 auto itev = findFirstEvent(events, T(X[0]));
3853 auto itev_end = events.end();
3854 // The above can still take you to end() if no events above X[0], so check
3855 // again.
3856 if (itev == itev_end)
3857 return;
3858
3859 // Find the first bin
3860 size_t bin = 0;
3861
3862 // Multiplier values
3863 double value;
3864 double error;
3865 double valueSquared;
3866 double errorSquared;
3867
3868 // If the tof is greater the first bin boundary, so we need to find the first
3869 // bin
3870 double tof = itev->tof();
3871 while (bin < x_size - 1) {
3872 // Within range?
3873 if ((tof >= X[bin]) && (tof < X[bin + 1]))
3874 break; // Stop increasing bin
3875 ++bin;
3876 }
3877
3878 // New bin! Find what you are multiplying!
3879 value = Y[bin];
3880 error = E[bin];
3881 valueSquared = value * value;
3882 errorSquared = error * error;
3883
3884 // Keep going through all the events
3885 while ((itev != itev_end) && (bin < x_size - 1)) {
3886 tof = itev->tof();
3887 while (bin < x_size - 1) {
3888 // Event is Within range?
3889 if ((tof >= X[bin]) && (tof < X[bin + 1])) {
3890 // Process this event. Multiply and calculate error.
3891 itev->m_errorSquared =
3892 static_cast<float>(itev->m_errorSquared * valueSquared + errorSquared * itev->m_weight * itev->m_weight);
3893 itev->m_weight *= static_cast<float>(value);
3894 break; // out of the bin-searching-while-loop
3895 }
3896 ++bin;
3897 if (bin >= x_size - 1)
3898 break;
3899
3900 // New bin! Find what you are multiplying!
3901 value = Y[bin];
3902 error = E[bin];
3903 valueSquared = value * value;
3904 errorSquared = error * error;
3905 }
3906 ++itev;
3907 }
3908}
3909
3910//------------------------------------------------------------------------------------------------
3931void EventList::multiply(const MantidVec &X, const MantidVec &Y, const MantidVec &E) {
3932 switch (eventType) {
3933 case TOF:
3934 // Switch to weights if needed.
3935 this->switchTo(WEIGHTED);
3936 // Fall through
3937
3938 case WEIGHTED:
3939 // Sorting by tof is necessary for the algorithm
3940 this->sortTof();
3942 break;
3943
3944 case WEIGHTED_NOTIME:
3945 // Sorting by tof is necessary for the algorithm
3946 this->sortTof();
3948 break;
3949 }
3950}
3951
3952//------------------------------------------------------------------------------------------------
3961template <class T>
3962void EventList::divideHistogramHelper(std::vector<T> &events, const MantidVec &X, const MantidVec &Y,
3963 const MantidVec &E) {
3964 // Validate inputs
3965 if ((X.size() < 2) || (Y.size() != E.size()) || (X.size() != 1 + Y.size())) {
3966 std::stringstream msg;
3967 msg << "EventList::divide() was given invalid size or "
3968 "inconsistent histogram arrays: X["
3969 << X.size() << "] "
3970 << "Y[" << Y.size() << " E[" << E.size() << "]";
3971 throw std::invalid_argument(msg.str());
3972 }
3973
3974 size_t x_size = X.size();
3975
3976 // Iterate through all events (sorted by tof)
3977 auto itev = findFirstEvent(events, T(X[0]));
3978 auto itev_end = events.end();
3979 // The above can still take you to end() if no events above X[0], so check
3980 // again.
3981 if (itev == itev_end)
3982 return;
3983
3984 // Find the first bin
3985 size_t bin = 0;
3986
3987 // Multiplier values
3988 double value;
3989 double error;
3990 double valError_over_value_squared;
3991
3992 // If the tof is greater the first bin boundary, so we need to find the first
3993 // bin
3994 double tof = itev->tof();
3995 while (bin < x_size - 1) {
3996 // Within range?
3997 if ((tof >= X[bin]) && (tof < X[bin + 1]))
3998 break; // Stop increasing bin
3999 ++bin;
4000 }
4001
4002 // New bin! Find what you are multiplying!
4003 value = Y[bin];
4004 error = E[bin];
4005
4006 // --- Division case ---
4007 if (value == 0) {
4008 value = std::numeric_limits<float>::quiet_NaN(); // Avoid divide by zero
4009 valError_over_value_squared = 0;
4010 } else
4011 valError_over_value_squared = error * error / (value * value);
4012
4013 // Keep going through all the events
4014 while ((itev != events.end()) && (bin < x_size - 1)) {
4015 tof = itev->tof();
4016 while (bin < x_size - 1) {
4017 // Event is Within range?
4018 if ((tof >= X[bin]) && (tof < X[bin + 1])) {
4019 // Process this event. Divide and calculate error.
4020 double newWeight = itev->m_weight / value;
4021 itev->m_errorSquared = static_cast<float>(
4022 newWeight * newWeight *
4023 ((itev->m_errorSquared / (itev->m_weight * itev->m_weight)) + valError_over_value_squared));
4024 itev->m_weight = static_cast<float>(newWeight);
4025 break; // out of the bin-searching-while-loop
4026 }
4027 ++bin;
4028 if (bin >= x_size - 1)
4029 break;
4030
4031 // New bin! Find what you are multiplying!
4032 value = Y[bin];
4033 error = E[bin];
4034
4035 // --- Division case ---
4036 if (value == 0) {
4037 value = std::numeric_limits<float>::quiet_NaN(); // Avoid divide by zero
4038 valError_over_value_squared = 0;
4039 } else
4040 valError_over_value_squared = error * error / (value * value);
4041 }
4042 ++itev;
4043 }
4044}
4045
4046//------------------------------------------------------------------------------------------------
4068void EventList::divide(const MantidVec &X, const MantidVec &Y, const MantidVec &E) {
4069 switch (eventType) {
4070 case TOF:
4071 // Switch to weights if needed.
4072 this->switchTo(WEIGHTED);
4073 // Fall through
4074
4075 case WEIGHTED:
4076 // Sorting by tof is necessary for the algorithm
4077 this->sortTof();
4078 divideHistogramHelper(*this->weightedEvents, X, Y, E);
4079 break;
4080
4081 case WEIGHTED_NOTIME:
4082 // Sorting by tof is necessary for the algorithm
4083 this->sortTof();
4085 break;
4086 }
4087}
4088
4089//------------------------------------------------------------------------------------------------
4099 if (value == 0.0)
4100 throw std::invalid_argument("EventList::divide() called with value of 0.0. Cannot divide by zero.");
4101 this->multiply(1.0 / value, 0.0);
4102 return *this;
4103}
4104
4105//------------------------------------------------------------------------------------------------
4115void EventList::divide(const double value, const double error) {
4116 if (value == 0.0)
4117 throw std::invalid_argument("EventList::divide() called with value of 0.0. Cannot divide by zero.");
4118 // Do nothing if dividing by exactly 1.0, no error
4119 else if (value == 1.0 && error == 0.0)
4120 return;
4121
4122 // We'll multiply by 1/value
4123 double invValue = 1.0 / value;
4124 // Relative error remains the same
4125 double invError = (error / value) * invValue;
4126
4127 this->multiply(invValue, invError);
4128}
4129
4130// ==============================================================================================
4131// ----------- SPLITTING AND FILTERING ---------------------------------------
4132// ==============================================================================================
4133//------------------------------------------------------------------------------------------------
4143void EventList::filterByPulseTime(Types::Core::DateAndTime start, Types::Core::DateAndTime stop,
4144 EventList &output) const {
4145 if (this == &output) {
4146 throw std::invalid_argument("In-place filtering is not allowed");
4147 }
4148
4149 // Start by sorting the event list by pulse time.
4150 this->sortPulseTime();
4151 // Clear the output
4152 output.clear();
4153 // Has to match the given type
4154 output.switchTo(eventType);
4155 output.setDetectorIDs(this->getDetectorIDs());
4156 output.setHistogram(m_histogram);
4157 output.setSortOrder(this->order);
4158
4159 // Iterate through all events (sorted by pulse time)
4160 switch (eventType) {
4161 case TOF:
4162 filterByPulseTimeHelper(*this->events, start, stop, *output.events);
4163 break;
4164 case WEIGHTED:
4165 filterByPulseTimeHelper(*this->weightedEvents, start, stop, *output.weightedEvents);
4166 break;
4167 case WEIGHTED_NOTIME:
4168 throw std::runtime_error("EventList::filterByPulseTime() called on an "
4169 "EventList that no longer has time information.");
4170 break;
4171 }
4172}
4173
4186void EventList::filterByPulseTime(Kernel::TimeROI const *timeRoi, EventList *output) const {
4187
4188 this->sortPulseTime();
4189 // Clear the output
4190
4191 output->clear();
4192 output->setDetectorIDs(this->getDetectorIDs());
4193 output->setHistogram(m_histogram);
4194 // Has to match the given type
4195 output->switchTo(eventType);
4196
4197 if ((timeRoi == nullptr) || (timeRoi->useAll())) {
4198 throw std::invalid_argument("TimeROI can not use all time");
4199 }
4200 const auto &intervals = timeRoi->toTimeIntervals();
4201 if (intervals.empty())
4202 return;
4203
4204 switch (eventType) {
4205 case TOF:
4206 filterByTimeROIHelper(*this->events, intervals, output);
4207 break;
4208 case WEIGHTED:
4209 filterByTimeROIHelper(*this->weightedEvents, intervals, output);
4210 break;
4211 case WEIGHTED_NOTIME:
4212 throw std::runtime_error("EventList::filterByPulseTime() called on an "
4213 "EventList that no longer has time information.");
4214 break;
4215 }
4216}
4217
4225template <class T>
4226void EventList::filterByPulseTimeHelper(std::vector<T> &events, DateAndTime start, DateAndTime stop,
4227 std::vector<T> &output) {
4228 std::copy_if(events.begin(), events.end(), std::back_inserter(output),
4229 [start, stop](const T &t) { return (t.m_pulsetime >= start) && (t.m_pulsetime < stop); });
4230}
4231
4238template <class T>
4239void EventList::filterByTimeROIHelper(std::vector<T> &events, const std::vector<Kernel::TimeInterval> &intervals,
4240 EventList *output) {
4241 // Iterate through the splitter at the same time
4242 auto itspl = intervals.cbegin();
4243 auto itspl_end = intervals.cend();
4244 // Iterate through all events (sorted by tof)
4245 auto itev = events.cbegin();
4246 auto itev_end = events.cend();
4247
4248 // This is the time of the first section. Anything before is thrown out.
4249 while (itspl != itspl_end) {
4250 // Get the splitting interval times and destination
4251 DateAndTime start = itspl->start();
4252 DateAndTime stop = itspl->stop();
4253 // Skip the events before the start of the time
4254 while ((itev != itev_end) && (itev->m_pulsetime < start))
4255 itev++;
4256
4257 // Go through all the events that are in the interval (if any)
4258 while ((itev != itev_end) && (itev->m_pulsetime < stop)) {
4259 // Copy the event into another
4260 const T eventCopy(*itev);
4261 output->addEventQuickly(eventCopy);
4262 ++itev;
4263 }
4264
4265 // Go to the next interval
4266 ++itspl;
4267 // But if we reached the end, then we are done.
4268 if (itspl == itspl_end)
4269 break;
4270
4271 // No need to keep looping through the filter if we are out of events
4272 if (itev == itev_end)
4273 break;
4274 }
4275 // Done!
4276}
4277
4283 if (timeRoi == nullptr) {
4284 throw std::runtime_error("TimeROI can not be a nullptr\n");
4285 }
4286 if (timeRoi->useAll()) {
4287 throw std::invalid_argument("TimeROI can not be empty\n");
4288 }
4289 // Start by sorting the event list by pulse time.
4290 this->sortPulseTime();
4291
4292 // Iterate through all events (sorted by pulse time)
4293 switch (eventType) {
4294 case TOF:
4295 filterInPlaceHelper(timeRoi, *this->events);
4296 break;
4297 case WEIGHTED:
4298 filterInPlaceHelper(timeRoi, *this->weightedEvents);
4299 break;
4300 case WEIGHTED_NOTIME:
4301 throw std::runtime_error("EventList::filterInPlace() called on an "
4302 "EventList that no longer has time information.");
4303 break;
4304 }
4305}
4306
4333template <class T>
4334void EventList::filterInPlaceHelper(Kernel::TimeROI const *timeRoi, typename std::vector<T> &events) {
4335
4336 const auto splitter = timeRoi->toTimeIntervals();
4337 // Iterate through the splitter at the same time
4338 auto itspl = splitter.cbegin();
4339 auto itspl_end = splitter.cend();
4340 DateAndTime start, stop;
4341
4342 // Iterate for the input
4343 auto itev = events.begin();
4344 auto itev_end = events.end();
4345
4346 // Iterator for the outputted list; will follow the input except when events
4347 // are dropped.
4348 auto itOut = events.begin();
4349
4350 // This is the time of the first section. Anything before is thrown out.
4351 while (itspl != itspl_end) {
4352 // Get the splitting interval times and destination
4353 start = itspl->start();
4354 stop = itspl->stop();
4355 // Skip the events before the start of the time
4356 while ((itev != itev_end) && (itev->m_pulsetime < start))
4357 itev++;
4358
4359 // Are we aligned in the input vs output?
4360 bool copyingInPlace = (itOut == itev);
4361 if (copyingInPlace) {
4362 while ((itev != itev_end) && (itev->m_pulsetime < stop))
4363 ++itev;
4364 // Make sure the iterators still match
4365 itOut = itev;
4366 } else {
4367 // Go through all the events that are in the interval (if any)
4368 while ((itev != itev_end) && (itev->m_pulsetime < stop)) {
4369 *itOut = *itev;
4370 ++itOut;
4371 ++itev;
4372 }
4373 }
4374
4375 // Go to the next interval
4376 ++itspl;
4377 // But if we reached the end, then we are done.
4378 if (itspl == itspl_end)
4379 break;
4380
4381 // No need to keep looping through the filter if we are out of events
4382 if (itev == itev_end)
4383 break;
4384
4385 } // Looping through entries in the splitter vector
4386
4387 // Ok, now resize the event list to reflect the fact that it (probably) shrank
4388 events.resize(std::size_t(std::distance(events.begin(), itOut)));
4389}
4390
4395void EventList::initializePartials(std::map<int, EventList *> partials) const {
4396
4397 // collect the state from events which is to be transferred to the partials
4398 bool removeDetIDs{true};
4399 const auto histogramLocal = this->getHistogram();
4400 const auto eventTypeLocal = this->getEventType();
4401
4402 // lambda expression initializing one partial
4403 auto initPartial = [&](EventList *partial) {
4404 partial->clear(removeDetIDs);
4405 partial->copyInfoFrom(*this);
4406 partial->setHistogram(histogramLocal);
4407 partial->switchTo(eventTypeLocal);
4408 };
4409
4410 // iterate over the partials
4411 std::for_each(partials.cbegin(), partials.cend(),
4412 [&](const std::pair<int, EventList *> &pair) { initPartial(pair.second); });
4413}
4414
4424void getEventsFrom(EventList &el, std::vector<TofEvent> *&events) { events = &el.getEvents(); }
4425void getEventsFrom(const EventList &el, std::vector<TofEvent> const *&events) { events = &el.getEvents(); }
4426
4436void getEventsFrom(EventList &el, std::vector<WeightedEvent> *&events) { events = &el.getWeightedEvents(); }
4437void getEventsFrom(const EventList &el, std::vector<WeightedEvent> const *&events) { events = &el.getWeightedEvents(); }
4438
4448void getEventsFrom(EventList &el, std::vector<WeightedEventNoTime> *&events) { events = &el.getWeightedEventsNoTime(); }
4449void getEventsFrom(const EventList &el, std::vector<WeightedEventNoTime> const *&events) {
4450 events = &el.getWeightedEventsNoTime();
4451}
4452
4460template <class T>
4461void EventList::convertUnitsViaTofHelper(typename std::vector<T> &events, Mantid::Kernel::Unit const *fromUnit,
4462 Mantid::Kernel::Unit const *toUnit) {
4463 for (auto &itev : events) {
4464 // Conver to TOF
4465 const double tof = fromUnit->singleToTOF(itev.m_tof);
4466 // And back from TOF to whatever
4467 itev.m_tof = toUnit->singleFromTOF(tof);
4468 }
4469}
4470
4471//--------------------------------------------------------------------------
4480 // Check for initialized
4481 if (!fromUnit || !toUnit)
4482 throw std::runtime_error("EventList::convertUnitsViaTof(): one of the units is NULL!");
4483 if (!fromUnit->isInitialized())
4484 throw std::runtime_error("EventList::convertUnitsViaTof(): fromUnit is not initialized!");
4485 if (!toUnit->isInitialized())
4486 throw std::runtime_error("EventList::convertUnitsViaTof(): toUnit is not initialized!");
4487
4488 switch (eventType) {
4489 case TOF:
4490 convertUnitsViaTofHelper(*this->events, fromUnit, toUnit);
4491 break;
4492 case WEIGHTED:
4493 convertUnitsViaTofHelper(*this->weightedEvents, fromUnit, toUnit);
4494 break;
4495 case WEIGHTED_NOTIME:
4496 convertUnitsViaTofHelper(*this->weightedEventsNoTime, fromUnit, toUnit);
4497 break;
4498 }
4499}
4500
4501//--------------------------------------------------------------------------
4508template <class T>
4509void EventList::convertUnitsQuicklyHelper(typename std::vector<T> &events, const double &factor, const double &power) {
4510 for (auto &event : events) {
4511 // Output unit = factor * (input) ^ power
4512 event.m_tof = factor * std::pow(event.m_tof, power);
4513 }
4514}
4515
4516//--------------------------------------------------------------------------
4522void EventList::convertUnitsQuickly(const double &factor, const double &power) {
4523 switch (eventType) {
4524 case TOF:
4525 convertUnitsQuicklyHelper(*this->events, factor, power);
4526 break;
4527 case WEIGHTED:
4528 convertUnitsQuicklyHelper(*this->weightedEvents, factor, power);
4529 break;
4530 case WEIGHTED_NOTIME:
4531 convertUnitsQuicklyHelper(*this->weightedEventsNoTime, factor, power);
4532 break;
4533 }
4534}
4535
4536HistogramData::Histogram &EventList::mutableHistogramRef() {
4537 if (mru)
4538 mru->deleteIndex(this);
4539 return m_histogram;
4540}
4541
4542void EventList::checkAndSanitizeHistogram(HistogramData::Histogram &histogram) {
4543 if (histogram.xMode() != HistogramData::Histogram::XMode::BinEdges)
4544 throw std::runtime_error("EventList: setting histogram with storage mode "
4545 "other than BinEdges is not possible");
4546 if (histogram.sharedY() || histogram.sharedE())
4547 throw std::runtime_error("EventList: setting histogram data with non-null "
4548 "Y or E data is not possible");
4549 // Avoid flushing of YMode: we only change X but YMode depends on events.
4550 if (histogram.yMode() == HistogramData::Histogram::YMode::Uninitialized)
4551 histogram.setYMode(m_histogram.yMode());
4552 if (histogram.yMode() != m_histogram.yMode())
4553 throw std::runtime_error("EventList: setting histogram data with different "
4554 "YMode is not possible");
4555}
4556
4558 throw std::runtime_error("EventList: setting Points as X data is not "
4559 "possible, only BinEdges are supported");
4560}
4561
4563 throw std::runtime_error("EventList: Cannot set Y or E data, these data are "
4564 "generated automatically based on the events");
4565}
4566
4567} // 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:69
const double m_tofShift
Definition EventList.cpp:70
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
IEventList : Interface to Mantid::DataObjects::EventList class, used to expose to PythonAPI.
Definition IEventList.h:26
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:96
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:171
virtual const HistogramData::HistogramE & e() const
Definition ISpectrum.h:173
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
Kernel::cow_ptr< HistogramData::HistogramX > sharedX() const
Definition ISpectrum.h:185
virtual const HistogramData::HistogramY & y() const
Definition ISpectrum.h:172
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:189
HistogramData::HistogramX & mutableX() &
Definition ISpectrum.h:175
A class for holding :
Definition EventList.h:57
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 size_t findExactBin(const MantidVec &X, const double tof, const size_t n_bin)
Find the exact bin which a TOF falls in starting from the provided estimated one.
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.
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 std::optional< size_t > findLinearBin(const MantidVec &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.
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.
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.
void generateHistogram(const MantidVec &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...
~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:330
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.
void generateErrorsHistogram(const MantidVec &Y, MantidVec &E) const
Generate the Error histogram for the provided counts histogram.
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:342
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.
void generateCountsHistogram(const MantidVec &X, MantidVec &Y) const
Fill a histogram given specified histogram bounds.
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:348
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.
static void multiplyHistogramHelper(std::vector< T > &events, const MantidVec &X, const MantidVec &Y, const MantidVec &E)
Helper method for multiplying an event list by a histogram with error.
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.
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 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.
EventSortType order
Last sorting order.
Definition EventList.h:345
void convertTofHelper(std::vector< T > &events, const std::function< double(double)> &func)
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.
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:187
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:184
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:351
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:336
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 generateCountsHistogramTimeAtSample(const MantidVec &X, MantidVec &Y, const double &tofFactor, const double &tofOffset) const
With respect to Time at Sample, fill a histogram given specified histogram bounds.
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)
static void divideHistogramHelper(std::vector< T > &events, const MantidVec &X, const MantidVec &Y, const MantidVec &E)
Helper method for dividing an event list by a histogram with error.
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.
static void histogramForWeightsHelper(const std::vector< T > &events, const MantidVec &X, MantidVec &Y, MantidVec &E)
Generates both the Y and E (error) histograms for an EventList with WeightedEvents.
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.
static std::optional< size_t > findLogBin(const MantidVec &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 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 generateHistogramTimeAtSample(const MantidVec &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 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 generateHistogramPulseTime(const MantidVec &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 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.
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:339
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:105
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...
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:333
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:18
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:41
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:141
Implements a copy on write data template.
Definition cow_ptr.h:41
EventType
What kind of event list is being stored.
Definition IEventList.h:18
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.
Definition EventList.cpp:99
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:32
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
FindBin(double step, double xmin)
std::optional< size_t > operator()(const Mantid::MantidVec &X, const double tof, const bool findExact)
std::optional< size_t >(* findBin)(const Mantid::MantidVec &, 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)