KVEventMixer.h

Go to the documentation of this file.

#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));
         }
      }
 
      // 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();
      }
 
      // DECORRELATION
      if (n_partA && bins[bin_number].partB_events.size() == decor_events) {
         // uncorrelated spectra using each partA of this event and all partB in each of last decor_events
         for (auto& partA : iter_A) {
            for (auto& e : bins[bin_number].partB_events) {
               for (auto& partB : e.particles) {
                  TreatNCorPair(bin_number, partA, partB);
               }
            }
         }
      }
      if (n_partB && bins[bin_number].partA_events.size() == decor_events) {
         // uncorrelated spectra using each partB of this event and all partA in each of last decor_events
         for (auto& partB : iter_B) {
            for (auto& e : bins[bin_number].partA_events) {
               for (auto& partA : e.particles) {
                  TreatNCorPair(bin_number, partB, partA);
               }
            }
         }
      }
 
   }
};
#endif // KVEventMixer_H