KVEventMixer.h

#ifndef KVEventMixer_H
#define KVEventMixer_H
 
#include <vector>
#include <deque>
 
template<typename ParticleInfoStruct, int NumBins = 1>
class KVEventMixer {
   struct event {
      std::vector<ParticleInfoStruct> particles;
 
      event() = default;
      event(event&&) = default;
      event(const event&) = default;
      event(const ParticleInfoStruct& p)
      {
         add(p);
      }
      void add(const ParticleInfoStruct& p)
      {
         particles.push_back(p);
      }
   };
 
   struct bin_data {
      std::deque<event> partA_events;
      std::deque<event> partB_events;
   };
   bin_data bins[NumBins];
 
   // number of events used for decorrelation (event mixing)
   typename std::deque<event>::size_type decor_events = 10;
 
public:
   KVEventMixer(typename std::deque<event>::size_type number_of_events_to_mix = 10)
      : decor_events{number_of_events_to_mix}
   {}
 
   template<typename TreatCorPairFunc, typename TreatNCorPairFunc, typename partA_iterator, typename partB_iterator>
   void ProcessEvent(int bin_number, partA_iterator iter_A, partB_iterator iter_B, TreatCorPairFunc TreatCorPair, TreatNCorPairFunc TreatNCorPair)
   {
      int n_partA(0), n_partB(0);
 
      for (auto& partA : iter_A) {
         // store partA in new correlated event in list
         ++n_partA;
         if (n_partA == 1) bins[bin_number].partA_events.push_back(ParticleInfoStruct(partA));
         else bins[bin_number].partA_events.back().add(ParticleInfoStruct(partA));
 
         for (auto& partB : iter_B) {
            if (n_partA == 1) { // each partA will loop over all partB, only add partB event once!
               // store partB in new correlated event in list
               ++n_partB;
               if (n_partB == 1) bins[bin_number].partB_events.push_back(ParticleInfoStruct(partB));
               else bins[bin_number].partB_events.back().add(ParticleInfoStruct(partB));
            }
            TreatCorPair(bin_number, partA, partB);
         }
      }
      if (!n_partA) {
         // treat events with no A particles, which may still have B particles
         // in this case the B particles should be added to the mixing buffers
         for (auto& partB : iter_B) {
            // store partB in new correlated event in list
            ++n_partB;
            if (n_partB == 1) bins[bin_number].partB_events.push_back(ParticleInfoStruct(partB));
            else bins[bin_number].partB_events.back().add(ParticleInfoStruct(partB));
         }
      }
 
      // DECORRELATION
      if (n_partA
            &&
            ((bins[bin_number].partB_events.size() == decor_events && !n_partB)
             || bins[bin_number].partB_events.size() == (decor_events + 1) && n_partB)) {
         // uncorrelated spectra using each partA of this event and all partB in each of last decor_events
         for (auto& partA : iter_A) {
            for (int i = 0; i < decor_events; ++i) {
               auto& e = bins[bin_number].partB_events[i];
               for (auto& partB : e.particles) {
                  TreatNCorPair(bin_number, partA, partB);
               }
            }
         }
      }
      if (n_partB           &&
            ((bins[bin_number].partA_events.size() == decor_events && !n_partA)
             || bins[bin_number].partA_events.size() == (decor_events + 1) && n_partA)) {
         // uncorrelated spectra using each partB of this event and all partA in each of last decor_events
         for (auto& partB : iter_B) {
            for (int i = 0; i < decor_events; ++i) {
               auto& e = bins[bin_number].partA_events[i];
               for (auto& partA : e.particles) {
                  TreatNCorPair(bin_number, partB, partA);
               }
            }
         }
      }
 
      // keep only last decor_events events in lists
      if (bins[bin_number].partA_events.size() > decor_events) {
         // remove oldest (first) event
         bins[bin_number].partA_events.pop_front();
      }
      if (bins[bin_number].partB_events.size() > decor_events) {
         // remove oldest (first) event
         bins[bin_number].partB_events.pop_front();
      }
 
   }
};
#endif // KVEventMixer_H