@@ -139,6 +139,27 @@ struct comparePulseTimeTOFDelta {
139139 int64_t deltaNano;
140140};
141141
142+ struct FindBin {
143+ double divisor;
144+ double offset;
145+ boost::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 ) {
148+ findBin = Mantid::DataObjects::EventList::findLogBin;
149+ divisor = 1 . / log1p (abs (step)); // use this to do change of base
150+ offset = log (xmin) * divisor;
151+ } else {
152+ findBin = Mantid::DataObjects::EventList::findLinearBin;
153+ divisor = 1 . / step;
154+ offset = xmin * divisor;
155+ }
156+ }
157+
158+ boost::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+
142163// / Constructor (empty)
143164// EventWorkspace is always histogram data and so is thus EventList
144165EventList::EventList (const EventType event_type)
@@ -1753,6 +1774,90 @@ void EventList::compressEvents(double tolerance, EventList *destination) {
17531774 destination->clearUnused ();
17541775}
17551776
1777+ template <class T >
1778+ inline void EventList::createWeightedEvents (std::vector<WeightedEventNoTime> &out, const std::vector<double > &tof,
1779+ const std::vector<T> &weight, const std::vector<T> &error) {
1780+ out.clear ();
1781+ for (size_t i = 0 ; i < weight.size (); ++i) {
1782+ const auto errors = static_cast <float >(error[i]);
1783+ if (errors > 0 )
1784+ out.emplace_back (tof[i], static_cast <float >(weight[i]), errors);
1785+ }
1786+ }
1787+
1788+ template <class T >
1789+ inline void EventList::processWeightedEvents (const std::vector<T> &events, std::vector<WeightedEventNoTime> &out,
1790+ const std::shared_ptr<std::vector<double >> histogram_bin_edges,
1791+ struct FindBin findBin) {
1792+ const auto NUM_BINS = histogram_bin_edges->size () - 1 ;
1793+ std::vector<double > tof (NUM_BINS, 0 .);
1794+ std::vector<double > normalization (NUM_BINS, 0 .);
1795+ std::vector<float > weight (NUM_BINS, 0 .);
1796+ std::vector<float > error (NUM_BINS, 0 .);
1797+ for (const auto &ev : events) {
1798+ const auto &bin_optional = findBin (*histogram_bin_edges.get (), ev.m_tof , false );
1799+ if (bin_optional) {
1800+ const auto bin = bin_optional.get ();
1801+ const double norm = calcNorm (ev.m_errorSquared );
1802+ tof[bin] += ev.m_tof * norm;
1803+ normalization[bin] += norm;
1804+ weight[bin] += ev.m_weight ;
1805+ error[bin] += ev.m_errorSquared ;
1806+ }
1807+ }
1808+
1809+ // normalize TOFs
1810+ std::transform (tof.begin (), tof.end (), normalization.begin (), tof.begin (), std::divides<double >());
1811+
1812+ createWeightedEvents (out, tof, weight, error);
1813+ }
1814+
1815+ void EventList::compressEvents (double tolerance, EventList *destination,
1816+ std::shared_ptr<std::vector<double >> histogram_bin_edges) {
1817+ if (!this ->empty ()) {
1818+ const auto NUM_BINS = histogram_bin_edges->size () - 1 ;
1819+ const auto xmin = static_cast <double >(histogram_bin_edges->front ());
1820+
1821+ auto findBin = FindBin (tolerance, xmin);
1822+
1823+ switch (eventType) {
1824+ case TOF: {
1825+ std::vector<double > tof (NUM_BINS, 0 );
1826+ std::vector<uint32_t > count (NUM_BINS, 0 );
1827+ for (const auto &ev : this ->events ) {
1828+ const auto &bin_optional = findBin (*histogram_bin_edges.get (), ev.m_tof , false );
1829+ if (bin_optional) {
1830+ const auto bin = bin_optional.get ();
1831+ count[bin]++;
1832+ tof[bin] += ev.m_tof ;
1833+ }
1834+ }
1835+
1836+ // average TOFs
1837+ std::transform (tof.begin (), tof.end (), count.begin (), tof.begin (), std::divides<double >());
1838+
1839+ createWeightedEvents (destination->weightedEventsNoTime , tof, count, count);
1840+ } break ;
1841+
1842+ case WEIGHTED:
1843+ processWeightedEvents (this ->weightedEvents , destination->weightedEventsNoTime , histogram_bin_edges, findBin);
1844+ break ;
1845+
1846+ case WEIGHTED_NOTIME:
1847+ processWeightedEvents (this ->weightedEventsNoTime , destination->weightedEventsNoTime , histogram_bin_edges,
1848+ findBin);
1849+ break ;
1850+ }
1851+ }
1852+
1853+ // In all cases, you end up WEIGHTED_NOTIME.
1854+ destination->eventType = WEIGHTED_NOTIME;
1855+ // The result will be sorted
1856+ destination->order = TOF_SORT;
1857+ // Empty out storage for vectors that are now unused.
1858+ destination->clearUnused ();
1859+ }
1860+
17561861void EventList::compressFatEvents (const double tolerance, const Mantid::Types::Core::DateAndTime &timeStart,
17571862 const double seconds, EventList *destination) {
17581863
@@ -2002,27 +2107,14 @@ void EventList::histogramForWeightsHelper(const std::vector<T> &events, const do
20022107 const auto xmin = X.front ();
20032108 const auto xmax = X.back ();
20042109
2005- double divisor, offset;
2006- boost::optional<size_t > (*findBin)(const MantidVec &, const double , const double , const double ,
2007- const bool ); // function pointer
2008-
2009- // store 1/divisor because multiplication is less flops than division
2010- if (step < 0 ) {
2011- findBin = findLogBin;
2012- divisor = 1 . / log1p (abs (step)); // use this to do change of base
2013- offset = log (xmin) * divisor;
2014- } else {
2015- findBin = findLinearBin;
2016- divisor = 1 . / step;
2017- offset = xmin * divisor;
2018- }
2110+ auto findBin = FindBin (step, xmin);
20192111
20202112 for (const T &ev : events) {
20212113 const double tof = ev.tof ();
20222114 if (tof < xmin || tof >= xmax)
20232115 continue ;
20242116
2025- boost::optional<size_t > n_bin = findBin (X, tof, divisor, offset, true );
2117+ boost::optional<size_t > n_bin = findBin (X, tof, true );
20262118
20272119 if (n_bin) {
20282120 Y[n_bin.get ()] += ev.weight ();
@@ -2486,27 +2578,14 @@ void EventList::generateCountsHistogram(const double step, const MantidVec &X, M
24862578 const auto xmin = X.front ();
24872579 const auto xmax = X.back ();
24882580
2489- double divisor, offset;
2490- boost::optional<size_t > (*findBin)(const MantidVec &, const double , const double , const double ,
2491- const bool ); // function pointer
2492-
2493- // store 1/divisor because multiplication is less flops than division
2494- if (step < 0 ) {
2495- findBin = findLogBin;
2496- divisor = 1 . / log1p (abs (step)); // use this to do change of base
2497- offset = log (xmin) * divisor;
2498- } else {
2499- findBin = findLinearBin;
2500- divisor = 1 . / step;
2501- offset = xmin * divisor;
2502- }
2581+ auto findBin = FindBin (step, xmin);
25032582
25042583 for (const TofEvent &ev : this ->events ) {
25052584 const double tof = ev.tof ();
25062585 if (tof < xmin || tof >= xmax)
25072586 continue ;
25082587
2509- const boost::optional<size_t > n_bin = findBin (X, tof, divisor, offset, true );
2588+ const boost::optional<size_t > n_bin = findBin (X, tof, true );
25102589
25112590 if (n_bin)
25122591 Y[n_bin.get ()]++;
0 commit comments