[PWGLF] added configurable to select physical primaries, tracks are now required to have an associated mc particle (#13465)

Co-authored-by: ALICE Action Bot <alibuild@cern.ch>

Author: GiorgioAlbertoLucia

PWGLF/TableProducer/QC/nucleiQC.cxx

@@ -76,6 +76,7 @@ struct nucleiQC {
 
   Configurable<bool> cfgFillTable{"cfgFillTable", true, "Fill output tree"};
   Configurable<bool> cfgDoCheckPdgCode{"cfgDoCheckPdgCode", true, "Should you only select tracks associated to a mc particle with the correct PDG code?"};
+  Configurable<bool> cfgFillOnlyPhysicalPrimaries{"cfgFillOnlyPhysicalPrimaries", true, "Should you only select physical primary particles?"};
   Configurable<LabeledArray<int>> cfgSpeciesToProcess{"cfgSpeciesToProcess", {nuclei::speciesToProcessDefault[0], nuclei::Species::kNspecies, 1, nuclei::names, {"processNucleus"}}, "Nuclei to process"};
   Configurable<LabeledArray<int>> cfgEventSelections{"cfgEventSelections", {nuclei::EvSelDefault[0], 8, 1, nuclei::eventSelectionLabels, nuclei::eventSelectionTitle}, "Event selections"};
   Configurable<int> cfgCentralityEstimator{"cfgCentralityEstimator", 0, "Centrality estimator (FV0A: 0, FT0M: 1, FT0A: 2, FT0C: 3)"};
@@ -136,18 +137,16 @@ struct nucleiQC {
     }
 
     for (int iSpecies = 0; iSpecies < static_cast<int>(nuclei::Species::kNspecies); iSpecies++) {
-      if (cfgSpeciesToProcess->get(iSpecies) == 1) {
+      if (cfgSpeciesToProcess->get(iSpecies) == 1)
         mSpeciesToProcess.emplace_back(iSpecies);
-      }
     }
 
     static_for<0, nuclei::kNspecies - 1>([&](auto iSpecies) {
       constexpr int kSpeciesCt = decltype(iSpecies)::value;
       const int kSpeciesRt = kSpeciesCt;
 
-      if (std::find(mSpeciesToProcess.begin(), mSpeciesToProcess.end(), kSpeciesCt) == mSpeciesToProcess.end()) {
+      if (std::find(mSpeciesToProcess.begin(), mSpeciesToProcess.end(), kSpeciesCt) == mSpeciesToProcess.end())
         return;
-      }
 
       float tpcBetheBlochParams[6];
       for (int iParam = 0; iParam < 6; iParam++) {
@@ -166,9 +165,9 @@ struct nucleiQC {
 
   void initCCDB(const aod::BCsWithTimestamps::iterator& bc)
   {
-    if (mRunNumber == bc.runNumber()) {
+    if (mRunNumber == bc.runNumber())
       return;
-    }
+
     auto timestamp = bc.timestamp();
     mRunNumber = bc.runNumber();
 
@@ -199,9 +198,8 @@ struct nucleiQC {
   bool pidSelection(const Ttrack& track, const Tcollision& collision)
   {
     constexpr int kIndex = iSpecies;
-    if (!nuclei::checkSpeciesValidity(kIndex)) {
+    if (!nuclei::checkSpeciesValidity(kIndex))
       std::runtime_error("species contains invalid nucleus kIndex");
-    }
 
     float centrality = nuclei::getCentrality(collision, cfgCentralityEstimator);
     float nsigmaTPC = mPidManagers[kIndex].getNSigmaTPC(track);
@@ -266,15 +264,14 @@ struct nucleiQC {
                                                                                  : iSpecies == nuclei::Species::kAl ? nuclei::Flags::kHe4
                                                                                                                     : 0;
 
-    if (track.hasTOF()) {
+    if (track.hasTOF())
       candidate.flags |= nuclei::Flags::kHasTOF;
-    }
-    if (track.hasTRD()) {
+
+    if (track.hasTRD())
       candidate.flags |= nuclei::Flags::kHasTRD;
-    }
-    if (!collision.selection_bit(o2::aod::evsel::kNoITSROFrameBorder)) {
+
+    if (!collision.selection_bit(o2::aod::evsel::kNoITSROFrameBorder))
       candidate.flags |= nuclei::Flags::kITSrof;
-    }
   }
 
   template <typename Tparticle>
@@ -304,9 +301,9 @@ struct nucleiQC {
   template <const bool isMc, typename Tcollision, typename Ttrack>
   nuclei::SlimCandidate fillCandidate(const int iSpecies, Tcollision const& collision, Ttrack const& track)
   {
-    if (!nuclei::checkSpeciesValidity(iSpecies)) {
+    if (!nuclei::checkSpeciesValidity(iSpecies))
       std::runtime_error("species contains invalid nucleus index");
-    }
+
     nuclei::SlimCandidate candidate = {.pt = track.pt() * track.sign(),
                                        .eta = track.eta(),
                                        .phi = track.phi(),
@@ -336,7 +333,6 @@ struct nucleiQC {
       }
     }
 
-    mNucleiCandidates.emplace_back(candidate);
     return candidate;
   }
 
@@ -363,9 +359,8 @@ struct nucleiQC {
   void fillHistograms(const nuclei::SlimCandidate& candidate)
   {
     constexpr int kIndex = iSpecies;
-    if (!nuclei::checkSpeciesValidity(kIndex)) {
+    if (!nuclei::checkSpeciesValidity(kIndex))
       std::runtime_error("species contains invalid nucleus kIndex");
-    }
 
     if (isGenerated) {
       const float ptGenerated = (kIndex == nuclei::Species::kPr || kIndex == nuclei::Species::kDe || kIndex == nuclei::Species::kTr) ? candidate.ptGenerated : candidate.ptGenerated / 2.f;
@@ -392,31 +387,38 @@ struct nucleiQC {
     auto bc = collision.template bc_as<aod::BCsWithTimestamps>();
     initCCDB(bc);
 
-    if (!nuclei::eventSelection(collision, mHistograms, cfgEventSelections, cfgCutVertex)) {
+    if (!nuclei::eventSelection(collision, mHistograms, cfgEventSelections, cfgCutVertex))
       return;
-    }
 
     for (const auto& track : tracks) {
 
       static_for<0, nuclei::kNspecies - 1>([&](auto iSpecies) {
         constexpr int kSpeciesCt = decltype(iSpecies)::value;
         const int kSpeciesRt = kSpeciesCt;
 
-        if (std::find(mSpeciesToProcess.begin(), mSpeciesToProcess.end(), kSpeciesRt) == mSpeciesToProcess.end()) {
+        if (std::find(mSpeciesToProcess.begin(), mSpeciesToProcess.end(), kSpeciesRt) == mSpeciesToProcess.end())
+          return;
+
+        if (!track.has_mcParticle())
           return;
-        }
 
-        if (track.has_mcParticle()) {
-          const auto& particle = track.mcParticle();
-          if (cfgDoCheckPdgCode) {
-            if (std::abs(particle.pdgCode()) != nuclei::pdgCodes[kSpeciesRt])
-              return;
-          }
-          if ((particle.y() - cfgRapidityCenterMass) < cfgRapidityMin || (particle.y() - cfgRapidityCenterMass) > cfgRapidityMax) {
+        const auto& particle = track.mcParticle();
+        if (cfgDoCheckPdgCode) {
+          if (std::abs(particle.pdgCode()) != nuclei::pdgCodes[kSpeciesRt])
             return;
-          }
         }
 
+        if ((particle.y() - cfgRapidityCenterMass) < cfgRapidityMin || (particle.y() - cfgRapidityCenterMass) > cfgRapidityMax)
+          return;
+
+        if (cfgFillOnlyPhysicalPrimaries && !particle.isPhysicalPrimary())
+          return;
+
+        nuclei::SlimCandidate candidate;
+        candidate = fillCandidate</*isMc*/ true>(kSpeciesCt, collision, track);
+        if ((candidate.flags >> 10) & 0b1)
+          LOG(info) << "track from material before track selection";
+
         mHistograms.fill(HIST(nuclei::cNames[kSpeciesCt]) + HIST("/hTrackSelections"), nuclei::trackSelection::kNoCuts);
         if (!trackSelection(track))
           return;
@@ -426,35 +428,35 @@ struct nucleiQC {
           return;
         mHistograms.fill(HIST(nuclei::cNames[kSpeciesCt]) + HIST("/hTrackSelections"), nuclei::trackSelection::kPidCuts);
 
-        nuclei::SlimCandidate candidate;
-        if (track.has_mcParticle()) {
-          candidate = fillCandidate</*isMc*/ true>(kSpeciesCt, collision, track);
-          dispatchFillHistograms</*isGenerated*/ true>(kSpeciesRt, candidate);
-        } else {
-          candidate = fillCandidate</*isMc*/ true>(kSpeciesCt, collision, track);
-        }
+        // nuclei::SlimCandidate candidate;
+        // candidate = fillCandidate</*isMc*/ true>(kSpeciesCt, collision, track);
 
+        mNucleiCandidates.emplace_back(candidate);
+        dispatchFillHistograms</*isGenerated*/ true>(kSpeciesRt, candidate);
         dispatchFillHistograms</*isGenerated*/ false>(kSpeciesRt, candidate);
       });
     }
 
     const int mcCollisionId = collision.mcCollisionId();
     auto mcParticlesThisCollision = mcParticles.sliceBy(mMcParticlesPerCollision, mcCollisionId);
     mcParticlesThisCollision.bindExternalIndices(&mcParticles);
+
     for (const auto& particle : mcParticlesThisCollision) {
-      if (std::find(mFilledMcParticleIds.begin(), mFilledMcParticleIds.end(), particle.globalIndex()) != mFilledMcParticleIds.end()) {
+      if (std::find(mFilledMcParticleIds.begin(), mFilledMcParticleIds.end(), particle.globalIndex()) != mFilledMcParticleIds.end())
         continue;
-      }
+
+      if (cfgFillOnlyPhysicalPrimaries && !particle.isPhysicalPrimary())
+        continue;
+
       int iSpecies = nuclei::getSpeciesFromPdg(particle.pdgCode());
-      if (std::find(mSpeciesToProcess.begin(), mSpeciesToProcess.end(), iSpecies) == mSpeciesToProcess.end()) {
+      if (std::find(mSpeciesToProcess.begin(), mSpeciesToProcess.end(), iSpecies) == mSpeciesToProcess.end())
         continue;
-      }
 
       nuclei::SlimCandidate candidate;
       fillNucleusFlagsPdgsMc(particle, candidate);
       fillNucleusGeneratedVariables(particle, candidate);
-      mNucleiCandidates.emplace_back(candidate);
 
+      mNucleiCandidates.emplace_back(candidate);
       dispatchFillHistograms</*isGenerated*/ true>(iSpecies, candidate);
     }