NucleiTask - Change reco condition on MC

Author: giovannimalfattore

PWGLF/Tasks/Nuspex/lfNucleiBATask.cxx

@@ -77,8 +77,8 @@ struct lfNucleiBATask {
 
   // Efficiency configurator
   std::unordered_set<int> effEvtSet;
-  bool effEvtSetReady = false;
-  Configurable<bool> enableEffEvtSet{"enableEffEvtSet", true, "If true, MCGen uses the event-set built by MCReco; if false, MCGen runs stand-alone."};
+  // bool effEvtSetReady = false;
+  Configurable<bool> requireRecoEventForGen{"requireRecoEventForGen", true, "Require at least one reco collision passing the event selection in processMCGen"};
 
   HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject};
   HistogramRegistry spectraGen{"spectraGen", {}, OutputObjHandlingPolicy::AnalysisObject, false, true};
@@ -112,6 +112,7 @@ struct lfNucleiBATask {
     Configurable<bool> useINELgt0cut{"useINELgt0cut", false, "Apply INEL>0 Event Selection in processData o processMC"};
     Configurable<bool> useINELgt1cut{"useINELgt1cut", false, "Apply INEL>1 Event Selection in processData o processMC"};
     Configurable<bool> enableGenVzCut{"enableGenVzCut", true, "Apply cut in z-Vertex in the processGen function"};
+    Configurable<bool> askGenFT0Coincidence{"askGenFT0Coincidence", false, "Require a generated charged-particle coincidence in FT0A and FT0C"};
   } evselOptions;
 
   // Set the multiplity event limits
@@ -167,9 +168,9 @@ struct lfNucleiBATask {
   // Set the axis used in this task
   ConfigurableAxis binsPercentile{"binsPercentile", {100, 0, 100}, "Centrality FT0M"};
   ConfigurableAxis binsPt{"binsPt", {VARIABLE_WIDTH, 0.0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.425, 0.45, 0.475, 0.5, 0.5125, 0.525, 0.5375, 0.55, 0.5625, 0.575, 0.5875, 0.6, 0.6125, 0.625, 0.6375, 0.65, 0.6625, 0.675, 0.6875, 0.7, 0.7125, 0.725, 0.7375, 0.75, 0.7625, 0.775, 0.7875, 0.8, 0.8125, 0.825, 0.8375, 0.85, 0.8625, 0.875, 0.8875, 0.9, 0.9125, 0.925, 0.9375, 0.95, 0.9625, 0.975, 0.9875, 1.0, 1.0125, 1.025, 1.0375, 1.05, 1.0625, 1.075, 1.0875, 1.1, 1.1125, 1.125, 1.1375, 1.15, 1.1625, 1.175, 1.1875, 1.2, 1.2125, 1.225, 1.2375, 1.25, 1.2625, 1.275, 1.2875, 1.3, 1.3125, 1.325, 1.3375, 1.35, 1.3625, 1.375, 1.3875, 1.4, 1.4125, 1.425, 1.4375, 1.45, 1.4625, 1.475, 1.4875, 1.5, 1.5125, 1.525, 1.5375, 1.55, 1.5625, 1.575, 1.5875, 1.6, 1.6125, 1.625, 1.6375, 1.65, 1.6625, 1.675, 1.6875, 1.7, 1.7125, 1.725, 1.7375, 1.75, 1.7625, 1.775, 1.7875, 1.8, 1.8125, 1.825, 1.8375, 1.85, 1.8625, 1.875, 1.8875, 1.9, 1.9125, 1.925, 1.9375, 1.95, 1.9625, 1.975, 1.9875, 2.0, 2.0625, 2.125, 2.1875, 2.25, 2.3125, 2.375, 2.4375, 2.5, 2.625, 2.75, 2.875, 3.0, 3.25, 3.5, 3.75, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0}, ""};
-  ConfigurableAxis binsPtHe{"binsPtHe", {VARIABLE_WIDTH, 1.0, 1.25, 1.50, 1.75, 2.0, 2.25, 2.50, 2.75, 3.0, 3.25, 3.50, 3.75, 4.0, 4.50, 5.0, 6.0, 7.0, 8.0}, ""};
-  ConfigurableAxis binsPtZHe{"binsPtZHe", {VARIABLE_WIDTH, 0.5, 0.625, 0.75, 0.875, 1.0, 1.125, 1.25, 1.375, 1.5, 1.625, 1.75, 1.875, 2.0, 2.25, 2.5, 3.0, 3.5, 4.0}, ""};
-  ConfigurableAxis binsPtHeMigration{"binsPtHeMigration", {VARIABLE_WIDTH, 0.10, 0.20, 0.30, 0.40, 0.50, 0.60, 0.70, 0.80, 0.90, 1.00, 1.10, 1.20, 1.30, 1.40, 1.50, 1.60, 1.70, 1.80, 1.90, 2.00, 2.10, 2.20, 2.30, 2.40, 2.50, 2.60, 2.70, 2.80, 2.90, 3.00, 3.25, 3.50, 3.75, 4.00, 4.50, 5.00, 6.00, 7.00, 8.00}, ""};
+  ConfigurableAxis binsPtHe{"binsPtHe", {VARIABLE_WIDTH, 0.0, 0.25, 0.5, 0.75, 1.0, 1.25, 1.50, 1.75, 2.0, 2.25, 2.50, 2.75, 3.0, 3.25, 3.50, 3.75, 4.0, 4.50, 5.0, 6.0, 7.0, 8.0}, ""};
+  ConfigurableAxis binsPtZHe{"binsPtZHe", {VARIABLE_WIDTH, 0.0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1.0, 1.125, 1.25, 1.375, 1.5, 1.625, 1.75, 1.875, 2.0, 2.25, 2.5, 3.0, 3.5, 4.0}, ""};
+  ConfigurableAxis binsPtHeMigration{"binsPtHeMigration", {VARIABLE_WIDTH, 0.0, 0.10, 0.20, 0.30, 0.40, 0.50, 0.60, 0.70, 0.80, 0.90, 1.00, 1.10, 1.20, 1.30, 1.40, 1.50, 1.60, 1.70, 1.80, 1.90, 2.00, 2.10, 2.20, 2.30, 2.40, 2.50, 2.60, 2.70, 2.80, 2.90, 3.00, 3.25, 3.50, 3.75, 4.00, 4.50, 5.00, 6.00, 7.00, 8.00}, ""};
 
   ConfigurableAxis binsdEdx{"binsdEdx", {600, 0.f, 3000.f}, ""};
   ConfigurableAxis binsBeta{"binsBeta", {120, 0.0, 1.2}, ""};
@@ -323,6 +324,38 @@ struct lfNucleiBATask {
     return 0.5f * std::log((e + pz) / (e - pz));
   }
 
+  // Added TVX check on MC
+  bool hasGenTVX(const aod::McParticles& particles, int mcCollisionIndex)
+  {
+    bool hasHitFT0A = false;
+    bool hasHitFT0C = false;
+
+    for (const auto& particle : particles) {
+      if (particle.mcCollisionId() != mcCollisionIndex)
+        continue;
+
+      if (!particle.isPhysicalPrimary())
+        continue;
+
+      auto* pdgInfo = pdgDB->GetParticle(particle.pdgCode());
+      if (!pdgInfo || pdgInfo->Charge() == 0.)
+        continue;
+
+      const float eta = particle.eta();
+
+      if (eta > 3.5f && eta < 4.9f)
+        hasHitFT0A = true;
+
+      if (eta > -3.3f && eta < -2.1f)
+        hasHitFT0C = true;
+
+      if (hasHitFT0A && hasHitFT0C)
+        return true;
+    }
+
+    return false;
+  }
+
   void init(o2::framework::InitContext& context)
   {
     if (initITSPID) {
@@ -333,7 +366,7 @@ struct lfNucleiBATask {
     }
 
     effEvtSet.clear();
-    effEvtSetReady = false;
+    // effEvtSetReady = false;
 
     const AxisSpec pAxis{binsPt, "#it{p} (GeV/#it{c})"};
     const AxisSpec ptAxis{binsPt, "#it{p}_{T} (GeV/#it{c})"};
@@ -2279,6 +2312,21 @@ struct lfNucleiBATask {
     }
     LOG(info) << " MC mode ";
     // MC histograms  -   all, primary, sec. from weak decay, sec. from material
+
+    if (enableDebug) {
+      spectraGen.add("hEffEvtSetVsSG", "effEvtSet vs SmallGroups;in effEvtSet;has accepted reco collision", HistType::kTH2D, {{2, -0.5, 1.5}, {2, -0.5, 1.5}});
+      auto hEffCheck = spectraGen.get<TH2>(HIST("hEffEvtSetVsSG"));
+      hEffCheck->GetXaxis()->SetBinLabel(1, "Not in effEvtSet");
+      hEffCheck->GetXaxis()->SetBinLabel(2, "In effEvtSet");
+      hEffCheck->GetYaxis()->SetBinLabel(1, "No accepted reco");
+      hEffCheck->GetYaxis()->SetBinLabel(2, "Accepted reco");
+
+      spectraGen.add("hEffEvtSetFalsePositive", "False positives in effEvtSet", HistType::kTH1D, {{2, 0., 2.}});
+      auto hFalsePos = spectraGen.get<TH1>(HIST("hEffEvtSetFalsePositive"));
+      hFalsePos->GetXaxis()->SetBinLabel(1, "No associated reco");
+      hFalsePos->GetXaxis()->SetBinLabel(2, "Reco exists, fails cuts");
+    }
+
     if (enableCentrality)
       spectraGen.add("histGenVetxZ", "PosZ generated events", HistType::kTH2F, {{1500, -15.f, 15.f, "Vertex Z (cm)"}, {binsPercentile, "Centrality FT0M"}});
     else
@@ -2415,7 +2463,6 @@ struct lfNucleiBATask {
     }
 
     LOG(info) << "MC Histograms defined";
-    // LOG(info) << "Histograms of lfNucleiBATask:";
     LOG(info) << "Histograms of lfNucleiBATask:";
     histos.print();
     spectraGen.print();
@@ -6382,12 +6429,13 @@ struct lfNucleiBATask {
     }
 
     float mcCentFT0M = -1.f;
+    const bool passesRecoVertex = event.posZ() >= cfgVzCutLow && event.posZ() <= cfgVzCutHigh;
 
-    if (doRecoSep && event.has_mcCollision()) {
+    if (doRecoSep && passesRecoVertex && event.has_mcCollision()) {
       const int mcIdx = event.mcCollisionId();
       if (mcIdx >= 0) {
         effEvtSet.insert(mcIdx);
-        effEvtSetReady = true;
+        // effEvtSetReady = true;
 
         auto mcColIter = mcCollisions.iteratorAt(mcIdx);
         mcCentFT0M = mcColIter.centFT0M();
@@ -6423,12 +6471,13 @@ struct lfNucleiBATask {
     }
 
     float mcCentFT0M = -1.f;
+    const bool passesRecoVertex = event.posZ() >= cfgVzCutLow && event.posZ() <= cfgVzCutHigh;
 
-    if (doRecoSep && event.has_mcCollision()) {
+    if (doRecoSep && passesRecoVertex && event.has_mcCollision()) {
       const int mcIdx = event.mcCollisionId();
       if (mcIdx >= 0) {
         effEvtSet.insert(mcIdx);
-        effEvtSetReady = true;
+        // effEvtSetReady = true;
 
         auto mcColIter = mcCollisions.iteratorAt(mcIdx);
         mcCentFT0M = mcColIter.centFT0M();
@@ -6582,21 +6631,53 @@ struct lfNucleiBATask {
   ////////////
 
   // LOOP OVER GENERATED MC PARTICLES
-  void processMCGen(soa::Join<aod::McCollisions, aod::McCentFT0Ms>::iterator const& mcCollision,
-                    aod::McParticles const& mcParticles)
+  // void processMCGen(soa::Join<aod::McCollisions, aod::McCentFT0Ms>::iterator const& mcCollision,
+  //                   aod::McParticles const& mcParticles)
+  void processMCGen(
+    soa::Join<aod::McCollisions, aod::McCentFT0Ms>::iterator const& mcCollision,
+    const soa::SmallGroups<soa::Join<EventCandidatesMC, o2::aod::PVMults>>& collisions,
+    aod::McParticles const& mcParticles)
   {
-    // Only events that are reconstructed
     const int mcIdx = mcCollision.globalIndex();
-    if (enableEffEvtSet) {
-      if (!effEvtSetReady)
+
+    // Add GenTVX condition
+    if (evselOptions.askGenFT0Coincidence) {
+      if (!hasGenTVX(mcParticles, mcIdx))
         return;
-      if (!effEvtSet.contains(mcIdx))
+    }
+
+    // Require Reco conditions (as in signal loss)
+    if (requireRecoEventForGen) {
+      bool hasAcceptedRecoCollision = false;
+      for (const auto& collision : collisions) {
+        const bool hasTVX = collision.selection_bit(aod::evsel::kIsTriggerTVX);
+        const bool hasNoTFB = collision.selection_bit(aod::evsel::kNoTimeFrameBorder);
+        const bool hasNoItsRofFB = collision.selection_bit(aod::evsel::kNoITSROFrameBorder);
+        if (evselOptions.useSel8 && !collision.sel8())
+          continue;
+        if (evselOptions.useTVXtrigger && !hasTVX)
+          continue;
+        if (evselOptions.removeTFBorder && !hasNoTFB)
+          continue;
+        if (evselOptions.removeITSROFBorder && !hasNoItsRofFB)
+          continue;
+        if (collision.posZ() < cfgVzCutLow || collision.posZ() > cfgVzCutHigh)
+          continue;
+        if (evselOptions.useINELgt0cut && !collision.isInelGt0())
+          continue;
+        if (evselOptions.useINELgt1cut && !collision.isInelGt1())
+          continue;
+        hasAcceptedRecoCollision = true;
+        break;
+      }
+      if (!hasAcceptedRecoCollision)
         return;
     }
 
     if (enableCentrality && (mcCollision.centFT0M() < cfgMultCutLow || mcCollision.centFT0M() > cfgMultCutHigh))
       return;
 
+    // Vz cut at Gen level
     if (evselOptions.enableGenVzCut) {
       if (mcCollision.posZ() < cfgVzCutLow || mcCollision.posZ() > cfgVzCutHigh)
         return;
@@ -6607,9 +6688,6 @@ struct lfNucleiBATask {
     else
       spectraGen.fill(HIST("histGenVetxZ"), mcCollision.posZ());
 
-    // const auto& particlesInCollision = mcParticles.sliceByCached(aod::mcparticle::mcCollisionId, mcIdx, cache);
-    // const auto particlesInCollision = mcParticles.sliceBy(perMCCol, mcIdx);
-
     for (const auto& mcParticleGen : mcParticles) {
       if (mcParticleGen.mcCollisionId() != mcIdx)
         continue;
@@ -6958,17 +7036,25 @@ struct lfNucleiBATask {
     const soa::SmallGroups<soa::Join<EventCandidatesMC, o2::aod::PVMults>>& collisions,
     aod::McParticles const& mcParticles)
   {
+    const int mcIdx = mcCollision.globalIndex();
+    const bool genFT0Coincidence = hasGenTVX(mcParticles, mcIdx);
     bool isINELgt0true = pwglf::isINELgtNmc(mcParticles, 0, pdgDB);
 
+    // GenFT0 coincidence
+    if (evselOptions.askGenFT0Coincidence && !genFT0Coincidence)
+      return;
+
     // EVENT LOSS DENOMINATOR
     // No cuts
     histoGen.fill(HIST("events/hMCGen"), 0.5);
     if (enableCentrality)
       histoGen.fill(HIST("events/hMCGenVsMult"), 0.5, mcCollision.centFT0M());
 
     // Vtz cut
-    if (mcCollision.posZ() < cfgVzCutLow || mcCollision.posZ() > cfgVzCutHigh)
-      return;
+    if (evselOptions.enableGenVzCut) {
+      if (mcCollision.posZ() < cfgVzCutLow || mcCollision.posZ() > cfgVzCutHigh)
+        return;
+    }
 
     histoGen.fill(HIST("events/hMCGen"), 1.5);
     if (enableCentrality)
@@ -7159,6 +7245,30 @@ struct lfNucleiBATask {
       histoGen.fill(HIST("events/hNRecoCollPerMCCollVsMult_INEL"), recoIdxINEL, mcCollision.centFT0M());
       histoGen.fill(HIST("events/hNRecoCollPerMCCollVsMult_INELgt0"), recoIdxINELgt0, mcCollision.centFT0M());
     }
+    // Add checker for EffEvtSet
+    const int mcIdxCheck = mcCollision.globalIndex();
+    const bool inEffEvtSet = effEvtSet.contains(mcIdxCheck);
+
+    if (enableDebug) {
+      double effEvtSetValue = 0.0;
+      double inelValue = 0.0;
+
+      if (inEffEvtSet) {
+        effEvtSetValue = 1.0;
+      }
+      if (atLeastOneINEL) {
+        inelValue = 1.0;
+      }
+
+      spectraGen.fill(HIST("hEffEvtSetVsSG"), effEvtSetValue, inelValue);
+
+      if (inEffEvtSet && !atLeastOneINEL) {
+        if (collisions.size() == 0)
+          spectraGen.fill(HIST("hEffEvtSetFalsePositive"), 0.5);
+        else
+          spectraGen.fill(HIST("hEffEvtSetFalsePositive"), 1.5);
+      }
+    }
 
     if (!atLeastOneINEL)
       return;