[PWGDQ] Write eventfilter to 32 highest bits of tag column in reducedEvents (#3922)

Author: sigurdnese

PWGDQ/Core/HistogramsLibrary.cxx

@@ -152,6 +152,11 @@ void o2::aod::dqhistograms::DefineHistograms(HistogramManager* hm, const char* h
       hm->AddHistogram(histClass, "R2EP_CentFT0C", "", true, 9, 0.0, 90.0, VarManager::kCentFT0C, 100, -1.0, 1.0, VarManager::kR2EP);
       hm->AddHistogram(histClass, "R3EP_CentFT0C", "", true, 9, 0.0, 90.0, VarManager::kCentFT0C, 100, -1.0, 1.0, VarManager::kR3EP);
     }
+    if (subGroupStr.Contains("filter")) {
+      hm->AddHistogram(histClass, "IsDoubleGap", "Is double gap", false, 2, -0.5, 1.5, VarManager::kIsDoubleGap);
+      hm->AddHistogram(histClass, "IsSingleGapA", "Is single gap on side A", false, 2, -0.5, 1.5, VarManager::kIsSingleGapA);
+      hm->AddHistogram(histClass, "IsSingleGapC", "Is single gap on side C", false, 2, -0.5, 1.5, VarManager::kIsSingleGapC);
+    }
   }
 
   if (groupStr.Contains("track")) {

PWGDQ/Core/VarManager.cxx

@@ -660,4 +660,10 @@ void VarManager::SetDefaultVarNames()
   fgVariableUnits[kRapCharmHadron] = " ";
   fgVariableNames[kPhiCharmHadron] = "#varphi (charm hadron)";
   fgVariableUnits[kPhiCharmHadron] = "rad.";
+  fgVariableNames[kIsDoubleGap] = "is double gap event";
+  fgVariableUnits[kIsDoubleGap] = "";
+  fgVariableNames[kIsSingleGapA] = "is single gap event side A";
+  fgVariableUnits[kIsSingleGapA] = "";
+  fgVariableNames[kIsSingleGapC] = "is single gap event side C";
+  fgVariableUnits[kIsSingleGapC] = "";
 }

PWGDQ/Core/VarManager.h

@@ -80,6 +80,7 @@ class VarManager : public TObject
     ReducedEventQvector = BIT(9),
     CollisionCent = BIT(10),
     CollisionMult = BIT(11),
+    EventFilter = BIT(12),
     Track = BIT(0),
     TrackCov = BIT(1),
     TrackExtra = BIT(2),
@@ -193,6 +194,10 @@ class VarManager : public TObject
     kR3SP,
     kR2EP,
     kR3EP,
+    kIsDoubleGap,  // Double rapidity gap
+    kIsSingleGapA, // Rapidity gap on side A
+    kIsSingleGapC, // Rapidity gap on side C
+    kIsSingleGap,  // Rapidity gap on either side
     kNEventWiseVariables,
 
     // Basic track/muon/pair wise variables
@@ -455,6 +460,12 @@ class VarManager : public TObject
     kD0barToKPi
   };
 
+  enum EventFilters {
+    kDoubleGap = 1,
+    kSingleGapA,
+    kSingleGapC
+  };
+
   static TString fgVariableNames[kNVars]; // variable names
   static TString fgVariableUnits[kNVars]; // variable units
   static void SetDefaultVarNames();
@@ -1005,6 +1016,13 @@ void VarManager::FillEvent(T const& event, float* values)
     values[kMCEventImpParam] = event.impactParameter();
   }
 
+  if constexpr ((fillMap & EventFilter) > 0) {
+    values[kIsDoubleGap] = (event.eventFilter() & (uint64_t(1) << kDoubleGap)) > 0;
+    values[kIsSingleGapA] = (event.eventFilter() & (uint64_t(1) << kSingleGapA)) > 0;
+    values[kIsSingleGapC] = (event.eventFilter() & (uint64_t(1) << kSingleGapC)) > 0;
+    values[kIsSingleGap] = values[kIsSingleGapA] || values[kIsSingleGapC];
+  }
+
   FillEventDerived(values);
 }
 

PWGDQ/TableProducer/tableMaker.cxx

@@ -108,6 +108,7 @@ DECLARE_SOA_TABLE(AmbiguousTracksFwd, "AOD", "AMBIGUOUSFWDTR", //! Table for Fwd
 
 constexpr static uint32_t gkEventFillMap = VarManager::ObjTypes::BC | VarManager::ObjTypes::Collision;
 constexpr static uint32_t gkEventFillMapWithMult = VarManager::ObjTypes::BC | VarManager::ObjTypes::Collision | VarManager::ObjTypes::CollisionMult;
+constexpr static uint32_t gkEventFillMapWithMultsAndEventFilter = VarManager::ObjTypes::BC | VarManager::ObjTypes::Collision | VarManager::ObjTypes::CollisionMult | VarManager::ObjTypes::EventFilter;
 constexpr static uint32_t gkEventFillMapWithCent = VarManager::ObjTypes::BC | VarManager::ObjTypes::Collision | VarManager::ObjTypes::CollisionCent;
 constexpr static uint32_t gkEventFillMapWithCentAndMults = VarManager::ObjTypes::BC | VarManager::ObjTypes::Collision | VarManager::ObjTypes::CollisionCent | VarManager::CollisionMult;
 // constexpr static uint32_t gkEventFillMapWithCentRun2 = VarManager::ObjTypes::BC | VarManager::ObjTypes::Collision | VarManager::ObjTypes::CollisionCentRun2; // Unused variable
@@ -368,7 +369,10 @@ struct TableMaker {
     if (collision.sel7()) {
       tag |= (uint64_t(1) << evsel::kNsel); //! SEL7 stored at position kNsel in the tag bit map
     }
-    // TODO: Add the event level decisions from the filtering task into the tag
+    // Put the 32 first bits of the event filter in the last 32 bits of the tag
+    if constexpr ((TEventFillMap & VarManager::ObjTypes::EventFilter) > 0) {
+      tag |= (collision.eventFilter() << 32);
+    }
 
     VarManager::ResetValues(0, VarManager::kNEventWiseVariables);
     // TODO: These variables cannot be filled in the VarManager for the moment as long as BCsWithTimestamps are used.
@@ -1208,7 +1212,7 @@ struct TableMaker {
     }
     (reinterpret_cast<TH2I*>(fStatsList->At(0)))->Fill(1.0, static_cast<float>(kNaliases));
     if (collision.eventFilter()) {
-      fullSkimming<gkEventFillMapWithMult, gkTrackFillMap, 0u>(collision, bcs, tracksBarrel, nullptr, nullptr, nullptr);
+      fullSkimming<gkEventFillMapWithMultsAndEventFilter, gkTrackFillMap, 0u>(collision, bcs, tracksBarrel, nullptr, nullptr, nullptr);
     }
   }
 

PWGDQ/Tasks/filterPbPb.cxx

@@ -48,7 +48,7 @@ void DefineHistograms(HistogramManager* histMan, TString histClasses);
 struct DQFilterPbPbTask {
   Produces<aod::DQEventFilter> eventFilter;
   OutputObj<TH1I> fStats{"Statistics"};
-  OutputObj<TH1F> fIsEventDGOutcome{TH1F("Filter outcome", "Filter outcome", 14, -1.5, 6.5)};
+  OutputObj<TH1F> fIsEventDGOutcome{TH1F("Filter outcome", "Filter outcome", 7, -0.5, 6.5)};
 
   Configurable<std::string> fConfigBarrelSelections{"cfgBarrelSels", "jpsiPID1:2:5", "<track-cut>:<nmin>:<nmax>,[<track-cut>:<nmin>:<nmax>],..."};
   Configurable<int> fConfigNDtColl{"cfgNDtColl", 4, "Number of standard deviations to consider in BC range"};
@@ -64,7 +64,6 @@ struct DQFilterPbPbTask {
   Configurable<bool> fConfigUseFV0{"cfgUseFV0", true, "Whether to use FV0 for veto"};
   Configurable<bool> fConfigUseFT0{"cfgUseFT0", true, "Whether to use FT0 for veto"};
   Configurable<bool> fConfigUseFDD{"cfgUseFDD", true, "Whether to use FDD for veto"};
-  Configurable<std::string> fConfigFITSides{"cfgFITSides", "both", "both, either, A, C, neither"};
   Configurable<bool> fConfigVetoForward{"cfgVetoForward", true, "Whether to veto on forward tracks"};
   Configurable<bool> fConfigVetoBarrel{"cfgVetoBarrel", false, "Whether to veto on barrel tracks"};
 
@@ -74,15 +73,12 @@ struct DQFilterPbPbTask {
   std::vector<int> fBarrelNminTracks; // minimal number of tracks in barrel
   std::vector<int> fBarrelNmaxTracks; // maximal number of tracks in barrel
 
-  int FITVetoSides = -1; // Integer to encode which side(s) of the FIT to use for veto
-  std::vector<std::string> FITVetoSidesOptions = {"both", "either", "A", "C", "neither"};
-
   // Helper function for selecting DG events
   template <typename TEvent, typename TBCs, typename TTracks, typename TMuons>
-  int isEventDG(TEvent const& collision, TBCs const& bcs, TTracks const& tracks, TMuons const& muons,
-                aod::FT0s& ft0s, aod::FV0As& fv0as, aod::FDDs& fdds,
-                std::vector<float> FITAmpLimits, int nDtColl, int minNBCs, int minNPVCs, int maxNPVCs, float maxFITTime,
-                bool useFV0, bool useFT0, bool useFDD, int FITSide, bool doVetoFwd, bool doVetoBarrel)
+  uint64_t isEventDG(TEvent const& collision, TBCs const& bcs, TTracks const& tracks, TMuons const& muons,
+                     aod::FT0s& ft0s, aod::FV0As& fv0as, aod::FDDs& fdds,
+                     std::vector<float> FITAmpLimits, int nDtColl, int minNBCs, int minNPVCs, int maxNPVCs, float maxFITTime,
+                     bool useFV0, bool useFT0, bool useFDD, bool doVetoFwd, bool doVetoBarrel)
   {
     fIsEventDGOutcome->Fill(0., 1.);
     // Find BC associated with collision
@@ -210,51 +206,54 @@ struct DQFilterPbPbTask {
     }
 
     // Compute FIT decision
-    bool FITDecision = 0;
-    if (FITSide == 0) {
-      FITDecision = isSideAClean && isSideCClean;
-    } else if (FITSide == 1) {
-      FITDecision = isSideAClean ^ isSideCClean;
-    } else if (FITSide == 2) {
-      FITDecision = isSideAClean;
-    } else if (FITSide == 3) {
-      FITDecision = isSideCClean;
-    } else if (FITSide == 4) {
-      FITDecision = 1;
+    uint64_t FITDecision = 0;
+    if (isSideAClean && isSideCClean) {
+      FITDecision |= (uint64_t(1) << VarManager::kDoubleGap);
+    } else if (isSideAClean && !isSideCClean) {
+      FITDecision |= (uint64_t(1) << VarManager::kSingleGapA);
+    } else if (!isSideAClean && isSideCClean) {
+      FITDecision |= (uint64_t(1) << VarManager::kSingleGapC);
     }
     if (!FITDecision) {
-      return 1;
+      fIsEventDGOutcome->Fill(2, 1);
+      return 0;
     }
 
     // Veto on activiy in either forward or barrel region
     if (doVetoFwd) {
       for (auto& muon : muons) {
         // Only care about muons with good timing (MID)
         if (muon.trackType() == 0 || muon.trackType() == 3) {
-          return 2;
+          fIsEventDGOutcome->Fill(3, 1);
+          return 0;
         }
       }
     }
     if (doVetoBarrel) {
       if (tracks.size() > 0) {
-        return 3;
+        fIsEventDGOutcome->Fill(4, 1);
+        return 0;
       }
     }
 
     // No global tracks which are not vtx tracks
     for (auto& track : tracks) {
       if (track.isGlobalTrack() && !track.isPVContributor()) {
-        return 4;
+        fIsEventDGOutcome->Fill(5, 1);
+        return 0;
       }
     }
 
     // Number of primary vertex contributors
     if (collision.numContrib() < minNPVCs || collision.numContrib() > maxNPVCs) {
-      return 5;
+      fIsEventDGOutcome->Fill(6, 1);
+      return 0;
     }
 
     // If we made it here, the event passed
-    return 0;
+    fIsEventDGOutcome->Fill(1, 1);
+    // Return filter bitmap corresponding to FIT decision
+    return FITDecision;
   }
 
   void DefineCuts()
@@ -290,15 +289,6 @@ struct DQFilterPbPbTask {
   void init(o2::framework::InitContext&)
   {
     DefineCuts();
-
-    for (int i = 0; i < 5; i++) {
-      if (fConfigFITSides.value == FITVetoSidesOptions[i]) {
-        FITVetoSides = i;
-      }
-    }
-    if (FITVetoSides == -1) {
-      LOGF(fatal, "Invalid choice of FIT side(s) for veto: %s", fConfigFITSides.value);
-    }
   }
 
   template <uint32_t TEventFillMap>
@@ -308,18 +298,13 @@ struct DQFilterPbPbTask {
     fStats->Fill(-2.0);
 
     std::vector<float> FITAmpLimits = {fConfigFV0AmpLimit, fConfigFT0AAmpLimit, fConfigFT0CAmpLimit, fConfigFDDAAmpLimit, fConfigFDDCAmpLimit};
-    int filterOutcome = isEventDG(collision, bcs, tracks, muons, ft0s, fv0as, fdds,
-                                  FITAmpLimits, fConfigNDtColl, fConfigMinNBCs, fConfigMinNPVCs, fConfigMaxNPVCs, fConfigMaxFITTime,
-                                  fConfigUseFV0, fConfigUseFT0, fConfigUseFDD, FITVetoSides, fConfigVetoForward, fConfigVetoBarrel);
-    fIsEventDGOutcome->Fill(filterOutcome + 1, 1);
-
-    // Don't need info on filter outcome anymore; reduce to a boolean
-    bool isDG = !filterOutcome;
-    fStats->Fill(-1.0, isDG);
-
-    // Compute event filter
-    uint64_t filter = 0;
-    filter |= isDG;
+    uint64_t filter = isEventDG(collision, bcs, tracks, muons, ft0s, fv0as, fdds,
+                                FITAmpLimits, fConfigNDtColl, fConfigMinNBCs, fConfigMinNPVCs, fConfigMaxNPVCs, fConfigMaxFITTime,
+                                fConfigUseFV0, fConfigUseFT0, fConfigUseFDD, fConfigVetoForward, fConfigVetoBarrel);
+
+    bool isSelected = filter;
+    fStats->Fill(-1.0, isSelected);
+
     eventFilter(filter);
   }