[PWGDQ] Adding protections to the prefiltering task; fix problem with exclusive MCSignals (AliceO2Group#8831)

Co-authored-by: Ionut Cristian Arsene <[email protected]>

Author: iarsene

PWGDQ/Core/MCSignal.cxx

@@ -27,7 +27,8 @@ MCSignal::MCSignal() : TNamed("", ""),
                        fExcludeCommonAncestor(false),
                        fTempAncestorLabel(-1),
                        fDecayChannelIsExclusive(false),
-                       fDecayChannelIsNotExclusive(false)
+                       fDecayChannelIsNotExclusive(false),
+                       fNAncestorDirectProngs(0)
 {
 }
 
@@ -39,7 +40,8 @@ MCSignal::MCSignal(int nProngs, const char* name /*= ""*/, const char* title /*=
                                                                                          fExcludeCommonAncestor(false),
                                                                                          fTempAncestorLabel(-1),
                                                                                          fDecayChannelIsExclusive(false),
-                                                                                         fDecayChannelIsNotExclusive(false)
+                                                                                         fDecayChannelIsNotExclusive(false),
+                                                                                         fNAncestorDirectProngs(0)
 {
   fProngs.reserve(nProngs);
 }
@@ -52,7 +54,8 @@ MCSignal::MCSignal(const char* name, const char* title, std::vector<MCProng> pro
                                                                                                                                                         fExcludeCommonAncestor(excludeCommonAncestor),
                                                                                                                                                         fTempAncestorLabel(-1),
                                                                                                                                                         fDecayChannelIsExclusive(false),
-                                                                                                                                                        fDecayChannelIsNotExclusive(false)
+                                                                                                                                                        fDecayChannelIsNotExclusive(false),
+                                                                                                                                                        fNAncestorDirectProngs(0)
 {
 }
 
@@ -82,6 +85,7 @@ void MCSignal::PrintConfig()
   cout << "Exclude common ancestor combinations: " << fExcludeCommonAncestor << endl;
   cout << "Decay channel is exclusive: " << fDecayChannelIsExclusive << endl;
   cout << "Decay channel is not exclusive: " << fDecayChannelIsNotExclusive << endl;
+  cout << "Decay channel direct prongs for the common ancestor: " << fNAncestorDirectProngs << endl;
   cout << "Printing " << fNProngs << "/" << fProngs.size() << " prongs:" << endl;
   int i = 0;
   for (auto& pr : fProngs) {

PWGDQ/Core/MCSignal.h

@@ -73,13 +73,15 @@ class MCSignal : public TNamed
 
   void SetProngs(std::vector<MCProng> prongs, std::vector<int8_t> commonAncestors);
   void AddProng(MCProng prong, int8_t commonAncestor = -1);
-  void SetDecayChannelIsExclusive(bool option = true)
+  void SetDecayChannelIsExclusive(int nProngs, bool option = true)
   {
     fDecayChannelIsExclusive = option;
+    fNAncestorDirectProngs = nProngs;
   }
-  void SetDecayChannelIsNotExclusive(bool option = true)
+  void SetDecayChannelIsNotExclusive(int nProngs, bool option = true)
   {
     fDecayChannelIsNotExclusive = option;
+    fNAncestorDirectProngs = nProngs;
   }
 
   int GetNProngs() const
@@ -98,6 +100,10 @@ class MCSignal : public TNamed
   {
     return fDecayChannelIsNotExclusive;
   }
+  int GetNAncestorDirectProngs() const
+  {
+    return fNAncestorDirectProngs;
+  }
 
   template <typename... T>
   bool CheckSignal(bool checkSources, const T&... args)
@@ -117,8 +123,9 @@ class MCSignal : public TNamed
   unsigned int fNProngs;                   // number of prongs
   std::vector<int8_t> fCommonAncestorIdxs; // index of the most recent ancestor, relative to each prong's history
   bool fExcludeCommonAncestor;             // explicitly request that there is no common ancestor
-  bool fDecayChannelIsExclusive;           // if true, then the indicated mother particle has a number of daughters which is equal to the number of prongs defined in this MC signal
-  bool fDecayChannelIsNotExclusive;        // if true, then the indicated mother particle has a number of daughters which is larger than the number of prongs defined in this MC signal
+  bool fDecayChannelIsExclusive;           // if true, then the indicated mother particle has a number of daughters which is equal to the number of direct prongs defined in this MC signal
+  bool fDecayChannelIsNotExclusive;        // if true, then the indicated mother particle has a number of daughters which is larger than the number of direct prongs defined in this MC signal
+  int fNAncestorDirectProngs;              // number of direct prongs belonging to the common ancestor specified by this signal
   int fTempAncestorLabel;
 
   template <typename T>
@@ -162,10 +169,10 @@ bool MCSignal::CheckProng(int i, bool checkSources, const T& track)
         //  If these numbers are equal, it means this decay MCSignal match is exclusive (there are no additional prongs for this mother besides the
         //     prongs defined here).
         if (currentMCParticle.has_daughters()) {
-          if (fDecayChannelIsExclusive && currentMCParticle.daughtersIds()[1] - currentMCParticle.daughtersIds()[0] + 1 != fNProngs) {
+          if (fDecayChannelIsExclusive && currentMCParticle.daughtersIds()[1] - currentMCParticle.daughtersIds()[0] + 1 != fNAncestorDirectProngs) {
             return false;
           }
-          if (fDecayChannelIsNotExclusive && currentMCParticle.daughtersIds()[1] - currentMCParticle.daughtersIds()[0] + 1 == fNProngs) {
+          if (fDecayChannelIsNotExclusive && currentMCParticle.daughtersIds()[1] - currentMCParticle.daughtersIds()[0] + 1 == fNAncestorDirectProngs) {
             return false;
           }
         }

PWGDQ/Core/MCSignalLibrary.cxx

@@ -712,14 +712,14 @@ MCSignal* o2::aod::dqmcsignals::GetMCSignal(const char* name)
     MCProng prong(2, {11, 443}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false});
     prong.SetSourceBit(0, MCProng::kPhysicalPrimary);
     signal = new MCSignal(name, "ee pairs from j/psi decays", {prong, prong}, {1, 1}); // signal at pair level
-    signal->SetDecayChannelIsExclusive(true);
+    signal->SetDecayChannelIsExclusive(2, true);
     return signal;
   }
   if (!nameStr.compare("eeFromJpsiNotExclusive")) {
     MCProng prong(2, {11, 443}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false});
     prong.SetSourceBit(0, MCProng::kPhysicalPrimary);
     signal = new MCSignal(name, "ee pairs from j/psi decays", {prong, prong}, {1, 1}); // signal at pair level
-    signal->SetDecayChannelIsNotExclusive(true);
+    signal->SetDecayChannelIsNotExclusive(2, true);
     return signal;
   }
   if (!nameStr.compare("eePrimaryFromPromptJPsi")) {
@@ -1145,15 +1145,15 @@ MCSignal* o2::aod::dqmcsignals::GetMCSignal(const char* name)
     MCProng pronge(3, {11, 443, 521}, {true, true, true}, {false, false, false}, {0, 0, 0}, {0, 0, 0}, {false, false, false});
     MCProng prongKaon(2, {321, 521}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false});
     signal = new MCSignal(name, "Kaon and electron pair from B+", {pronge, pronge, prongKaon}, {2, 2, 1});
-    signal->SetDecayChannelIsExclusive(true);
+    signal->SetDecayChannelIsExclusive(2, true);
     return signal;
   }
 
   if (!nameStr.compare("eeKaonFromBplusNotExclusive")) {
     MCProng pronge(3, {11, 443, 521}, {true, true, true}, {false, false, false}, {0, 0, 0}, {0, 0, 0}, {false, false, false});
     MCProng prongKaon(2, {321, 521}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false});
     signal = new MCSignal(name, "Kaon and electron pair from B+", {pronge, pronge, prongKaon}, {2, 2, 1});
-    signal->SetDecayChannelIsNotExclusive(true);
+    signal->SetDecayChannelIsNotExclusive(2, true);
     return signal;
   }
 

PWGDQ/Tasks/dqEfficiency_withAssoc.cxx

@@ -919,7 +919,7 @@ struct AnalysisPrefilterSelection {
   Configurable<std::string> fConfigPrefilterPairCut{"cfgPrefilterPairCut", "", "Prefilter pair cut"};
   Configurable<std::string> fConfigTrackCuts{"cfgTrackCuts", "", "Track cuts for which to run the prefilter"};
   // Track related options
-  Configurable<bool> fPropTrack{"cfgPropTrack", true, "Propgate tracks to associated collision to recalculate DCA and momentum vector"};
+  Configurable<bool> fPropTrack{"cfgPropTrack", false, "Propgate tracks to associated collision to recalculate DCA and momentum vector"};
 
   std::map<uint32_t, uint32_t> fPrefilterMap;
   AnalysisCompositeCut* fPairCut;
@@ -934,23 +934,39 @@ struct AnalysisPrefilterSelection {
       return;
     }
 
+    bool runPrefilter = true;
     // get the list of track cuts to be prefiltered
     TString trackCutsStr = fConfigTrackCuts.value;
     TObjArray* objArrayTrackCuts = nullptr;
     if (!trackCutsStr.IsNull()) {
       objArrayTrackCuts = trackCutsStr.Tokenize(",");
+      if (objArrayTrackCuts == nullptr) {
+        runPrefilter = false;
+      }
+    } else {
+      LOG(warn) << " No track cuts to prefilter! Prefilter will not be run";
+      runPrefilter = false;
+    }
+    // get the cut to be used as loose selection
+    TString prefilterTrackCutStr = fConfigPrefilterTrackCut.value;
+    if (prefilterTrackCutStr.IsNull()) {
+      LOG(warn) << " No prefilter loose selection specified! Prefilter will not be run";
+      runPrefilter = false;
     }
 
-    // get the list of cuts that were computed in the barrel track-selection task and create a bit mask
-    //  to mark just the ones we want to apply a prefilter on
     fPrefilterMask = 0;
     fPrefilterCutBit = -1;
-    string trackCuts;
-    getTaskOptionValue<string>(context, "analysis-track-selection", "cfgTrackCuts", trackCuts, false);
-    TString allTrackCutsStr = trackCuts;
-    TString prefilterTrackCutStr = fConfigPrefilterTrackCut.value;
-    if (!trackCutsStr.IsNull()) {
+    if (runPrefilter) {
+      // get the list of cuts that were computed in the barrel track-selection task and create a bit mask
+      //  to mark just the ones we want to apply a prefilter on
+      string trackCuts;
+      getTaskOptionValue<string>(context, "analysis-track-selection", "cfgTrackCuts", trackCuts, false);
+      TString allTrackCutsStr = trackCuts;
+
       std::unique_ptr<TObjArray> objArray(allTrackCutsStr.Tokenize(","));
+      if (objArray == nullptr) {
+        LOG(fatal) << " Not getting any track cuts from the barrel-track-selection ";
+      }
       if (objArray->FindObject(prefilterTrackCutStr.Data()) == nullptr) {
         LOG(fatal) << " Prefilter track cut not among the cuts calculated by the track-selection task! ";
       }
@@ -963,18 +979,15 @@ struct AnalysisPrefilterSelection {
           fPrefilterCutBit = icut;
         }
       }
+      // setup the prefilter pair cut
+      fPairCut = new AnalysisCompositeCut(true);
+      TString pairCutStr = fConfigPrefilterPairCut.value;
+      if (!pairCutStr.IsNull()) {
+        fPairCut = dqcuts::GetCompositeCut(pairCutStr.Data());
+      }
     }
-    // NOTE: If no prefilter loose cut is specified to be "prefiltered" or no track cuts are asked to be prefiltered,
-    //        then make a warning. In the processing, the produced table will just lead to a neutral behaviour (no prefilter applied)
     if (fPrefilterMask == static_cast<uint32_t>(0) || fPrefilterCutBit < 0) {
-      LOG(warn) << "No loose cut or track cuts for prefiltering are specified. This task will do nothing.";
-    }
-
-    // setup the prefilter pair cut
-    fPairCut = new AnalysisCompositeCut(true);
-    TString pairCutStr = fConfigPrefilterPairCut.value;
-    if (!pairCutStr.IsNull()) {
-      fPairCut = dqcuts::GetCompositeCut(pairCutStr.Data());
+      LOG(warn) << "No specified loose cut or track cuts for prefiltering. This task will do nothing.";
     }
 
     VarManager::SetUseVars(AnalysisCut::fgUsedVars); // provide the list of required variables so that VarManager knows what to fill
@@ -1039,18 +1052,26 @@ struct AnalysisPrefilterSelection {
         runPrefilter<gkTrackFillMap>(event, groupedAssocs, tracks);
       }
     }
+
     uint32_t mymap = -1;
-    for (auto& assoc : assocs) {
-      auto track = assoc.template reducedtrack_as<MyBarrelTracks>();
-      mymap = -1;
-      // If cuts were not configured, then produce a map with all 1's
-      if (fPrefilterCutBit < 0 || fPrefilterMask == 0) {
-        prefilter(mymap);
-      } else if (fPrefilterMap.find(track.globalIndex()) != fPrefilterMap.end()) {
-        // NOTE: publish the bitwise negated bits (~), so there will be zeroes for cuts that failed the prefiltering and 1 everywhere else
-        mymap = ~fPrefilterMap[track.globalIndex()];
+    // If cuts were not configured, then produce a map with all 1's and publish it for all associations
+    if (fPrefilterCutBit < 0 || fPrefilterMask == 0) {
+      for (int i = 0; i < assocs.size(); ++i) {
         prefilter(mymap);
       }
+    } else {
+      for (auto& assoc : assocs) {
+        // TODO: just use the index from the assoc (no need to cast the whole track)
+        auto track = assoc.template reducedtrack_as<MyBarrelTracks>();
+        mymap = -1;
+        if (fPrefilterMap.find(track.globalIndex()) != fPrefilterMap.end()) {
+          // NOTE: publish the bitwise negated bits (~), so there will be zeroes for cuts that failed the prefiltering and 1 everywhere else
+          mymap = ~fPrefilterMap[track.globalIndex()];
+          prefilter(mymap);
+        } else {
+          prefilter(mymap); // track did not pass the prefilter selections, so publish just 1's
+        }
+      }
     }
   }
 

PWGDQ/Tasks/tableReader_withAssoc.cxx

@@ -856,32 +856,47 @@ struct AnalysisPrefilterSelection {
   uint32_t fPrefilterMask;
   int fPrefilterCutBit;
 
+  Preslice<aod::ReducedTracksAssoc> trackAssocsPerCollision = aod::reducedtrack_association::reducedeventId;
+
   void init(o2::framework::InitContext& context)
   {
     if (context.mOptions.get<bool>("processDummy")) {
       return;
     }
 
+    bool runPrefilter = true;
     // get the list of track cuts to be prefiltered
     TString trackCutsStr = fConfigTrackCuts.value;
     TObjArray* objArrayTrackCuts = nullptr;
     if (!trackCutsStr.IsNull()) {
       objArrayTrackCuts = trackCutsStr.Tokenize(",");
+      if (objArrayTrackCuts == nullptr) {
+        runPrefilter = false;
+      }
+    } else {
+      LOG(warn) << " No track cuts to prefilter! Prefilter will not be run";
+      runPrefilter = false;
     }
-    if (objArrayTrackCuts->GetEntries() == 0) {
-      LOG(fatal) << " No track cuts to prefilter!";
+    // get the cut to be used as loose selection
+    TString prefilterTrackCutStr = fConfigPrefilterTrackCut.value;
+    if (prefilterTrackCutStr.IsNull()) {
+      LOG(warn) << " No prefilter loose selection specified! Prefilter will not be run";
+      runPrefilter = false;
     }
 
-    // get the list of cuts that were computed in the barrel track-selection task and create a bit mask
-    //  to mark just the ones we want to apply a prefilter on
     fPrefilterMask = static_cast<uint32_t>(0);
     fPrefilterCutBit = -1;
-    string trackCuts;
-    getTaskOptionValue<string>(context, "analysis-track-selection", "cfgTrackCuts", trackCuts, false);
-    TString allTrackCutsStr = trackCuts;
-    TString prefilterTrackCutStr = fConfigPrefilterTrackCut.value;
-    if (!trackCutsStr.IsNull()) {
+    if (runPrefilter) {
+      // get the list of cuts that were computed in the barrel track-selection task and create a bit mask
+      //  to mark just the ones we want to apply a prefilter on
+      string trackCuts;
+      getTaskOptionValue<string>(context, "analysis-track-selection", "cfgTrackCuts", trackCuts, false);
+      TString allTrackCutsStr = trackCuts;
+
       std::unique_ptr<TObjArray> objArray(allTrackCutsStr.Tokenize(","));
+      if (objArray == nullptr) {
+        LOG(fatal) << " Not getting any track cuts from the barrel-track-selection ";
+      }
       if (objArray->FindObject(prefilterTrackCutStr.Data()) == nullptr) {
         LOG(fatal) << " Prefilter track cut not among the cuts calculated by the track-selection task! ";
       }
@@ -894,18 +909,17 @@ struct AnalysisPrefilterSelection {
           fPrefilterCutBit = icut;
         }
       }
+      // setup the prefilter pair cut
+      fPairCut = new AnalysisCompositeCut(true);
+      TString pairCutStr = fConfigPrefilterPairCut.value;
+      if (!pairCutStr.IsNull()) {
+        fPairCut = dqcuts::GetCompositeCut(pairCutStr.Data());
+      }
     }
     if (fPrefilterMask == static_cast<uint32_t>(0) || fPrefilterCutBit < 0) {
       LOG(warn) << "No specified loose cut or track cuts for prefiltering. This task will do nothing.";
     }
 
-    // setup the prefilter pair cut
-    fPairCut = new AnalysisCompositeCut(true);
-    TString pairCutStr = fConfigPrefilterPairCut.value;
-    if (!pairCutStr.IsNull()) {
-      fPairCut = dqcuts::GetCompositeCut(pairCutStr.Data());
-    }
-
     VarManager::SetUseVars(AnalysisCut::fgUsedVars); // provide the list of required variables so that VarManager knows what to fill
     VarManager::SetDefaultVarNames();
 
@@ -955,11 +969,8 @@ struct AnalysisPrefilterSelection {
     } // end loop over combinations
   }
 
-  Preslice<aod::ReducedTracksAssoc> trackAssocsPerCollision = aod::reducedtrack_association::reducedeventId;
-
   void processBarrelSkimmed(MyEvents const& events, soa::Join<aod::ReducedTracksAssoc, aod::BarrelTrackCuts> const& assocs, MyBarrelTracks const& tracks)
   {
-
     fPrefilterMap.clear();
 
     for (auto& event : events) {
@@ -968,18 +979,26 @@ struct AnalysisPrefilterSelection {
         runPrefilter<gkTrackFillMap>(groupedAssocs, tracks);
       }
     }
+
     uint32_t mymap = -1;
-    for (auto& assoc : assocs) {
-      auto track = assoc.template reducedtrack_as<MyBarrelTracks>();
-      mymap = -1;
-      // If cuts were not configured, then produce a map with all 1's
-      if (fPrefilterCutBit < 0 || fPrefilterMask == 0) {
-        prefilter(mymap);
-      } else if (fPrefilterMap.find(track.globalIndex()) != fPrefilterMap.end()) {
-        // NOTE: publish the bitwise negated bits (~), so there will be zeroes for cuts that failed the prefiltering and 1 everywhere else
-        mymap = ~fPrefilterMap[track.globalIndex()];
+    // If cuts were not configured, then produce a map with all 1's and publish it for all associations
+    if (fPrefilterCutBit < 0 || fPrefilterMask == 0) {
+      for (int i = 0; i < assocs.size(); ++i) {
         prefilter(mymap);
       }
+    } else {
+      for (auto& assoc : assocs) {
+        // TODO: just use the index from the assoc (no need to cast the whole track)
+        auto track = assoc.template reducedtrack_as<MyBarrelTracks>();
+        mymap = -1;
+        if (fPrefilterMap.find(track.globalIndex()) != fPrefilterMap.end()) {
+          // NOTE: publish the bitwise negated bits (~), so there will be zeroes for cuts that failed the prefiltering and 1 everywhere else
+          mymap = ~fPrefilterMap[track.globalIndex()];
+          prefilter(mymap);
+        } else {
+          prefilter(mymap); // track did not pass the prefilter selections, so publish just 1's
+        }
+      }
     }
   }