back to rolling mixing

Author: JinjooSeo

PWGDQ/Core/MixingHandler.cxx

@@ -194,18 +194,3 @@ int MixingHandler::GetBinFromCategory(VarManager::Variables var, int category) c
   truncatedCategory /= norm;
   return truncatedCategory % nBins[ivar];
 }
-
-//_________________________________________________________________________
-void MixingHandler::SetCategoryBinCenters(int category, float* values) const
-{
-  //
-  // set the mixing variables to the bin centers of this category (used for the leftover mixing)
-  //
-  for (auto const& [var, pos] : fVariables) {
-    int bin = GetBinFromCategory(static_cast<VarManager::Variables>(var), category);
-    if (bin < 0) {
-      continue;
-    }
-    values[var] = 0.5 * (fVariableLimits[pos][bin] + fVariableLimits[pos][bin + 1]);
-  }
-}

PWGDQ/Core/MixingHandler.h

@@ -23,7 +23,6 @@
 
 #include <Rtypes.h>
 
-#include <algorithm>
 #include <array>
 #include <cstdint>
 #include <iostream>
@@ -36,8 +35,6 @@ class MixingHandler : public TNamed
  public:
   // number of track cuts which fit in the 32-bit filtering masks
   static constexpr int NMaxCuts = 32;
-  // smallest pool depth for which a mixed pair can be built
-  static constexpr int16_t MinPoolDepth = 2;
 
   // Struct to define track properties relevant for mixing and few utility functions
   struct MixingTrack {
@@ -48,6 +45,8 @@ class MixingHandler : public TNamed
     // globalIndex is unique only within a dataframe, so the dataframe sequence is part of the track identity
     uint64_t dataFrameSequence = 0;
     uint64_t trackGlobalIndex = 0;
+    // electric charge of the track; 0 means "not set" and disables the charge dependent pair variables
+    int8_t sign = 0;
     bool IsSamePhysicalTrack(const MixingTrack& other) const
     {
       return dataFrameSequence == other.dataFrameSequence && trackGlobalIndex == other.trackGlobalIndex;
@@ -56,7 +55,8 @@ class MixingHandler : public TNamed
     void ClearBit(uint32_t mask) { filteringFlags &= ~mask; }
     void Print() const
     {
-      std::cout << "pt: " << pt << ", eta: " << eta << ", phi: " << phi << ", filteringFlags: " << filteringFlags
+      std::cout << "pt: " << pt << ", eta: " << eta << ", phi: " << phi << ", sign: " << static_cast<int>(sign)
+                << ", filteringFlags: " << filteringFlags
                 << ", dataframe: " << dataFrameSequence << ", track: " << trackGlobalIndex << std::endl;
     }
   };
@@ -84,19 +84,6 @@ class MixingHandler : public TNamed
     }
     // Clear bits in the filtering mask.
     void ClearFilteringMask(uint32_t mask) { filteringMask &= ~mask; }
-    // clear the cut bits from all tracks and the filtering mask; remove tracks with no active bits left
-    void ClearBits(uint32_t mask)
-    {
-      for (auto& track : tracks1) {
-        track.ClearBit(mask);
-      }
-      tracks1.erase(std::remove_if(tracks1.begin(), tracks1.end(), [](auto const& track) { return track.filteringFlags == 0; }), tracks1.end());
-      for (auto& track : tracks2) {
-        track.ClearBit(mask);
-      }
-      tracks2.erase(std::remove_if(tracks2.begin(), tracks2.end(), [](auto const& track) { return track.filteringFlags == 0; }), tracks2.end());
-      ClearFilteringMask(mask);
-    }
     // 1) increment the counters for a given track cut bit mask and if the counters reached the pool depth,
     // 2) clear the corresponding bit in the tracks filtering flags to exclude them from further mixing
     // 3) for each track, if there are no more active bits in the filtering mask, then remove the track from the event
@@ -184,34 +171,16 @@ class MixingHandler : public TNamed
       CleanPool();
       events.push_back(event);
     }
-    // fixed-block mixing: AddEvent() until GetMixingMask() reports full cuts, mix, then ClearBits()
-    void AddEvent(const MixingEvent& event) { events.push_back(event); }
-    // bit mask of the cuts for which at least poolDepth events are in the pool
-    uint32_t GetMixingMask(int16_t poolDepth) const
+    // Same, but the stored events are aged only for the cuts in agingMask. Passing the filtering mask of the
+    // incoming event ages an event only for the cuts for which a mixed pair was actually produced, so that the
+    // pool depth is a number of mixed partners and not a number of arrivals.
+    void UpdatePool(const MixingEvent& event, int16_t poolDepth, uint32_t agingMask)
     {
-      std::array<int16_t, NMaxCuts> counts = {0};
-      for (auto const& event : events) {
-        for (int icut = 0; icut < NMaxCuts; ++icut) {
-          if (event.filteringMask & (static_cast<uint32_t>(1) << icut)) {
-            counts[icut]++;
-          }
-        }
-      }
-      uint32_t fullMask = 0;
-      for (int icut = 0; icut < NMaxCuts; ++icut) {
-        if (counts[icut] >= poolDepth) {
-          fullMask |= static_cast<uint32_t>(1) << icut;
-        }
-      }
-      return fullMask;
-    }
-    // clear the given cut bits from all events in the pool and remove the events with no tracks left
-    void ClearBits(uint32_t mask)
-    {
-      for (auto& event : events) {
-        event.ClearBits(mask);
+      for (auto& poolEvent : events) { // o2-linter: disable=const-ref-in-for-loop (the events are modified)
+        poolEvent.IncrementCounters(agingMask, poolDepth);
       }
       CleanPool();
+      events.push_back(event);
     }
     // getter for the events in the pool
     const std::vector<MixingEvent>& GetEvents() const { return events; }
@@ -240,14 +209,11 @@ class MixingHandler : public TNamed
   // int GetMixingVariable(VarManager::Variables var); // returns the position in the internal varible list of the handler. Useful for checks, mostly
   // std::vector<float> GetMixingVariableLimits(VarManager::Variables var);
   MixingPool& GetPool(int category) { return fPools[category]; }
-  std::map<int, MixingPool>& GetPools() { return fPools; }
   int16_t GetPoolDepth() const { return fPoolDepth; }
 
   void Init();
   int FindEventCategory(float* values);
   int GetBinFromCategory(VarManager::Variables var, int category) const;
-  // set the mixing variables to the bin centers of the given category
-  void SetCategoryBinCenters(int category, float* values) const;
 
  private:
   MixingHandler(const MixingHandler& handler);

PWGDQ/Tasks/tableReader_withAssoc.h

@@ -73,7 +73,6 @@
 #include <cstdio>
 #include <cstdlib>
 #include <cstring>
-#include <functional>
 #include <iterator>
 #include <map>
 #include <memory>
@@ -82,26 +81,6 @@
 #include <variant>
 #include <vector>
 
-// OutputObj with a callback run on the first dereference, i.e. right before the end-of-stream snapshot
-// (tasks cannot register their own EndOfStream callback, CallbackService::set replaces the framework one)
-template <typename T>
-struct FinalizingOutputObj : public o2::framework::OutputObj<T> {
-  using o2::framework::OutputObj<T>::OutputObj;
-  std::function<void()> finalizeBeforeSnapshot;
-  bool finalized = false;
-
-  T& operator*()
-  {
-    if (!finalized) {
-      finalized = true;
-      if (finalizeBeforeSnapshot) {
-        finalizeBeforeSnapshot();
-      }
-    }
-    return o2::framework::OutputObj<T>::operator*();
-  }
-};
-
 // Some definitions
 namespace o2::aod
 {
@@ -1376,7 +1355,7 @@ struct AnalysisSameEventPairing {
   TH1D* ResoFlowEP = nullptr;
   int fCurrentRun = -1; // needed to detect if the run changed and trigger update of calibrations etc.
 
-  FinalizingOutputObj<THashList> fOutputList{"output"};
+  o2::framework::OutputObj<THashList> fOutputList{"output"};
 
   struct : o2::framework::ConfigurableGroup {
     o2::framework::Configurable<std::string> track{"cfgTrackCuts", "jpsiO2MCdebugCuts2", "Comma separated list of barrel track cuts"};
@@ -1741,9 +1720,6 @@ struct AnalysisSameEventPairing {
     }
 
     if (fConfigRunMixingAcrossTFs) {
-      if (fConfigMixingDepth.value < MixingHandler::MinPoolDepth) {
-        LOGF(fatal, "cfgMixingDepth must be at least %d", MixingHandler::MinPoolDepth);
-      }
       if (fNCutsBarrel > MixingHandler::NMaxCuts) {
         LOGF(fatal, "Across-TF mixing supports at most %d barrel track-cut bits, got %d", MixingHandler::NMaxCuts, fNCutsBarrel);
       }
@@ -1849,12 +1825,6 @@ struct AnalysisSameEventPairing {
       o2::aod::dqhistograms::AddHistogramsFromJSON(fHistMan, fConfigAddJSONHistograms.value.c_str()); // ad-hoc histograms via JSON
       VarManager::SetUseVars(fHistMan->GetUsedVars());                                                // provide the list of required variables so that VarManager knows what to fill
       fOutputList.setObject(fHistMan->GetMainHistogramList());
-      // mix the events left in the pools at the end of the stream
-      fOutputList.finalizeBeforeSnapshot = [this]() {
-        if (fConfigRunMixingAcrossTFs) {
-          runLeftoverMixing();
-        }
-      };
     }
   }
 
@@ -1932,97 +1902,6 @@ struct AnalysisSameEventPairing {
     }
   }
 
-  // Mix all the events in the pool for the cuts in mixingMask (event-wise variables are those of the current event)
-  void runEventMixing(MixingHandler::MixingPool& pool, uint32_t mixingMask)
-  {
-    auto const& events = pool.GetEvents();
-    // each pair of events is mixed once
-    for (size_t iev1 = 0; iev1 < events.size(); iev1++) {
-      for (size_t iev2 = iev1 + 1; iev2 < events.size(); iev2++) {
-        auto const& mixingEvent = events[iev1];
-        auto const& poolEvent = events[iev2];
-        if (!(mixingEvent.filteringMask & poolEvent.filteringMask & mixingMask)) {
-          continue;
-        }
-        for (auto const& t1 : mixingEvent.tracks1) {
-          // run +- pairing
-          for (auto const& t2 : poolEvent.tracks2) {
-            // check the two-track filter for the mixed pair
-            uint32_t mixedTwoTrackFilter = t1.filteringFlags & t2.filteringFlags & mixingMask;
-            if (!mixedTwoTrackFilter) {
-              continue;
-            }
-            VarManager::FillPairMEAcrossTFs(t1, t2);
-            for (int icut = 0; icut < fNCutsBarrel; icut++) {
-              if (mixedTwoTrackFilter & (static_cast<uint32_t>(1) << icut)) {
-                fHistMan->FillHistClass(Form("PairsBarrelMEPM_%s", fTrackCuts[icut].Data()), dqtablereader_helpers::varValues());
-              }
-            }
-          }
-          // run ++ pairing
-          for (auto const& t2 : poolEvent.tracks1) {
-            // check the two-track filter for the mixed pair and skip the same track associated to both collisions
-            uint32_t mixedTwoTrackFilter = t1.filteringFlags & t2.filteringFlags & mixingMask;
-            if (!mixedTwoTrackFilter || t1.IsSamePhysicalTrack(t2)) {
-              continue;
-            }
-            VarManager::FillPairMEAcrossTFs(t1, t2);
-            for (int icut = 0; icut < fNCutsBarrel; icut++) {
-              if (mixedTwoTrackFilter & (static_cast<uint32_t>(1) << icut)) {
-                fHistMan->FillHistClass(Form("PairsBarrelMEPP_%s", fTrackCuts[icut].Data()), dqtablereader_helpers::varValues());
-              }
-            }
-          }
-        }
-        for (auto const& t1 : mixingEvent.tracks2) {
-          // run -+ pairing
-          for (auto const& t2 : poolEvent.tracks1) {
-            // check the two-track filter for the mixed pair
-            uint32_t mixedTwoTrackFilter = t1.filteringFlags & t2.filteringFlags & mixingMask;
-            if (!mixedTwoTrackFilter) {
-              continue;
-            }
-            VarManager::FillPairMEAcrossTFs(t1, t2);
-            for (int icut = 0; icut < fNCutsBarrel; icut++) {
-              if (mixedTwoTrackFilter & (static_cast<uint32_t>(1) << icut)) {
-                fHistMan->FillHistClass(Form("PairsBarrelMEPM_%s", fTrackCuts[icut].Data()), dqtablereader_helpers::varValues());
-              }
-            }
-          }
-          // run -- pairing
-          for (auto const& t2 : poolEvent.tracks2) {
-            // check the two-track filter for the mixed pair and skip the same track associated to both collisions
-            uint32_t mixedTwoTrackFilter = t1.filteringFlags & t2.filteringFlags & mixingMask;
-            if (!mixedTwoTrackFilter || t1.IsSamePhysicalTrack(t2)) {
-              continue;
-            }
-            VarManager::FillPairMEAcrossTFs(t1, t2);
-            for (int icut = 0; icut < fNCutsBarrel; icut++) {
-              if (mixedTwoTrackFilter & (static_cast<uint32_t>(1) << icut)) {
-                fHistMan->FillHistClass(Form("PairsBarrelMEMM_%s", fTrackCuts[icut].Data()), dqtablereader_helpers::varValues());
-              }
-            }
-          }
-        }
-      }
-    }
-  }
-
-  // Mix the events left in the pools (end of stream or run change), with the mixing variables set to the category bin centers
-  void runLeftoverMixing()
-  {
-    for (auto& [category, pool] : fMixingHandler.GetPools()) { // o2-linter: disable=const-ref-in-for-loop (the pools are modified)
-      const uint32_t mixingMask = pool.GetMixingMask(MixingHandler::MinPoolDepth);
-      if (!mixingMask) {
-        continue;
-      }
-      VarManager::ResetValues(0, VarManager::kNEventWiseVariables);
-      fMixingHandler.SetCategoryBinCenters(category, dqtablereader_helpers::varValues());
-      runEventMixing(pool, mixingMask);
-      pool.ClearBits(mixingMask);
-    }
-  }
-
   // Template function to run same event pairing (barrel-barrel, muon-muon, barrel-muon)
   template <bool TTwoProngFitter, int TPairType, uint32_t TEventFillMap, uint32_t TTrackFillMap, typename TEvents, typename TTrackAssocs, typename TTracks>
   void runSameEventPairing(TEvents const& events, o2::framework::Preslice<TTrackAssocs>& preslice, TTrackAssocs const& assocs, TTracks const& /*tracks*/)
@@ -2031,9 +1910,6 @@ struct AnalysisSameEventPairing {
       if (fCurrentRun != events.begin().runNumber()) {
         if (fConfigRunMixingAcrossTFs) {
           // do not mix events from different runs
-          if (fCurrentRun >= 0) {
-            runLeftoverMixing();
-          }
           fMixingHandler.ClearPools();
         }
         initParamsFromCCDB(events.begin().timestamp(), events.begin().runNumber(), TTwoProngFitter);
@@ -2636,7 +2512,7 @@ struct AnalysisSameEventPairing {
               continue;
             }
             auto t1 = assoc.template reducedtrack_as<TTracks>();
-            MixingHandler::MixingTrack mixingTrack(t1.pt(), t1.eta(), t1.phi(), trackFilterForMixing, currentMixingDataFrameSequence, static_cast<uint64_t>(assoc.reducedtrackId()));
+            MixingHandler::MixingTrack mixingTrack(t1.pt(), t1.eta(), t1.phi(), trackFilterForMixing, currentMixingDataFrameSequence, static_cast<uint64_t>(assoc.reducedtrackId()), static_cast<int8_t>(t1.sign()));
             if (t1.sign() > 0) {
               mixingEvent.AddTrack1(mixingTrack);
             } else {
@@ -2647,15 +2523,72 @@ struct AnalysisSameEventPairing {
         if (mixingEvent.tracks1.empty() && mixingEvent.tracks2.empty()) {
           continue;
         }
-        // 2) add the event to the pool corresponding to this event
+        // 2) run the mixing with the events in the pool corresponding to this event
         auto& pool = fMixingHandler.GetPool(mixingCategory);
-        pool.AddEvent(mixingEvent);
-        // 3) mix all the events in the pool for the cuts which reached the pool depth
-        uint32_t mixingMask = pool.GetMixingMask(fMixingHandler.GetPoolDepth());
-        if (mixingMask) {
-          runEventMixing(pool, mixingMask);
-          pool.ClearBits(mixingMask);
+        for (auto const& poolEvent : pool.GetEvents()) {
+          for (auto const& t1 : mixingEvent.tracks1) {
+            // run +- pairing
+            for (auto const& t2 : poolEvent.tracks2) {
+              // check the two-track filter for the mixed pair
+              uint32_t mixedTwoTrackFilter = t1.filteringFlags & t2.filteringFlags;
+              if (!mixedTwoTrackFilter) {
+                continue;
+              }
+              VarManager::FillPairMEAcrossTFs(t1, t2);
+              for (int icut = 0; icut < ncuts; icut++) {
+                if (mixedTwoTrackFilter & (static_cast<uint32_t>(1) << icut)) {
+                  fHistMan->FillHistClass(Form("PairsBarrelMEPM_%s", fTrackCuts[icut].Data()), dqtablereader_helpers::varValues());
+                }
+              }
+            }
+            // run ++ pairing
+            for (auto const& t2 : poolEvent.tracks1) {
+              // check the two-track filter for the mixed pair and skip the same track associated to both collisions
+              uint32_t mixedTwoTrackFilter = t1.filteringFlags & t2.filteringFlags;
+              if (!mixedTwoTrackFilter || t1.IsSamePhysicalTrack(t2)) {
+                continue;
+              }
+              VarManager::FillPairMEAcrossTFs(t1, t2);
+              for (int icut = 0; icut < ncuts; icut++) {
+                if (mixedTwoTrackFilter & (static_cast<uint32_t>(1) << icut)) {
+                  fHistMan->FillHistClass(Form("PairsBarrelMEPP_%s", fTrackCuts[icut].Data()), dqtablereader_helpers::varValues());
+                }
+              }
+            }
+          }
+          for (auto const& t1 : mixingEvent.tracks2) {
+            // run -+ pairing
+            for (auto const& t2 : poolEvent.tracks1) {
+              // check the two-track filter for the mixed pair
+              uint32_t mixedTwoTrackFilter = t1.filteringFlags & t2.filteringFlags;
+              if (!mixedTwoTrackFilter) {
+                continue;
+              }
+              VarManager::FillPairMEAcrossTFs(t1, t2);
+              for (int icut = 0; icut < ncuts; icut++) {
+                if (mixedTwoTrackFilter & (static_cast<uint32_t>(1) << icut)) {
+                  fHistMan->FillHistClass(Form("PairsBarrelMEPM_%s", fTrackCuts[icut].Data()), dqtablereader_helpers::varValues());
+                }
+              }
+            }
+            // run -- pairing
+            for (auto const& t2 : poolEvent.tracks2) {
+              // check the two-track filter for the mixed pair and skip the same track associated to both collisions
+              uint32_t mixedTwoTrackFilter = t1.filteringFlags & t2.filteringFlags;
+              if (!mixedTwoTrackFilter || t1.IsSamePhysicalTrack(t2)) {
+                continue;
+              }
+              VarManager::FillPairMEAcrossTFs(t1, t2);
+              for (int icut = 0; icut < ncuts; icut++) {
+                if (mixedTwoTrackFilter & (static_cast<uint32_t>(1) << icut)) {
+                  fHistMan->FillHistClass(Form("PairsBarrelMEMM_%s", fTrackCuts[icut].Data()), dqtablereader_helpers::varValues());
+                }
+              }
+            }
+          }
         }
+        // 3) add the current event to the pool
+        pool.UpdatePool(mixingEvent, fMixingHandler.GetPoolDepth(), mixingEvent.filteringMask);
         // pool.Print();
       }
     } // end loop over events