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