From d5b5a40fa803a0b9fdef059d13c92f6e6b7a98f2 Mon Sep 17 00:00:00 2001
From: rspijkers <rik.spijkers@gmail.com>
Date: Tue, 1 Jul 2025 16:45:51 +0200
Subject: [PATCH] Add gen level histograms with matched reconstructed event

---
 .../Tasks/Strangeness/cascadecorrelations.cxx | 72 +++++++++++++++----
 1 file changed, 57 insertions(+), 15 deletions(-)

diff --git a/PWGLF/Tasks/Strangeness/cascadecorrelations.cxx b/PWGLF/Tasks/Strangeness/cascadecorrelations.cxx
index 8adfd81d430..d1382b57370 100644
--- a/PWGLF/Tasks/Strangeness/cascadecorrelations.cxx
+++ b/PWGLF/Tasks/Strangeness/cascadecorrelations.cxx
@@ -246,13 +246,22 @@ struct CascadeSelector {
       registry.add("gen/hXiPlus", "hXiPlus", HistType::kTH2F, {ptAxis, rapidityAxis});
       registry.add("gen/hOmegaMinus", "hOmegaMinus", HistType::kTH2F, {ptAxis, rapidityAxis});
       registry.add("gen/hOmegaPlus", "hOmegaPlus", HistType::kTH2F, {ptAxis, rapidityAxis});
+
+      registry.add("genwithrec/hXiMinus", "hXiMinus", HistType::kTH2F, {ptAxis, rapidityAxis});
+      registry.add("genwithrec/hXiPlus", "hXiPlus", HistType::kTH2F, {ptAxis, rapidityAxis});
+      registry.add("genwithrec/hOmegaMinus", "hOmegaMinus", HistType::kTH2F, {ptAxis, rapidityAxis});
+      registry.add("genwithrec/hOmegaPlus", "hOmegaPlus", HistType::kTH2F, {ptAxis, rapidityAxis});
+
+      registry.add("genwithrec/hNevents", "hNevents", HistType::kTH1F, {{1, 0, 1, "N generated events with reconstructed event"}});
+      registry.add("gen/hNevents", "hNevents", HistType::kTH1F, {{1, 0, 1, "N generated events"}});
     }
   }
 
-  bool eventSelection(MyCollisions::iterator const& collision)
+  template <typename TCollision>
+  bool eventSelection(TCollision const& collision)
   {
     if (useTrigger) {
-      auto bc = collision.bc_as<aod::BCsWithTimestamps>();
+      auto bc = collision.template bc_as<aod::BCsWithTimestamps>();
       zorro.initCCDB(ccdb.service, bc.runNumber(), bc.timestamp(), triggerList);
       bool eventTrigger = zorro.isSelected(bc.globalBC());
       if (eventTrigger) {
@@ -282,7 +291,8 @@ struct CascadeSelector {
     return true;
   }
 
-  void fillMatchedHistos(LabeledCascades::iterator rec, int flag, MyCollisions::iterator collision)
+  template <typename TCollision>
+  void fillMatchedHistos(LabeledCascades::iterator rec, int flag, TCollision collision)
   {
     if (flag == 0)
       return;
@@ -335,8 +345,8 @@ struct CascadeSelector {
     }
   }
 
-  template <typename TCascade>
-  int processCandidate(TCascade const& casc, MyCollisions::iterator const& collision)
+  template <typename TCascade, typename TCollision>
+  int processCandidate(TCascade const& casc, TCollision const& collision)
   {
     // these are the tracks:
     auto bachTrack = casc.template bachelor_as<FullTracksExtIUWithPID>();
@@ -477,8 +487,11 @@ struct CascadeSelector {
     return 0;
   } // processCandidate
 
-  void processGenMC(aod::McCollision const&, soa::SmallGroups<soa::Join<aod::McCollisionLabels, MyCollisionsMult>> const&, aod::McParticles const& mcParticles)
+  void processGenMC(aod::McCollision const&, soa::SmallGroups<soa::Join<aod::McCollisionLabels, MyCollisions>> const& collisions, aod::McParticles const& mcParticles)
   {
+    // N gen events without any event selection or matched reco event
+    registry.fill(HIST("gen/hNevents"), 0);
+
     for (auto const& mcPart : mcParticles) {
       if (!mcPart.isPhysicalPrimary())
         continue;
@@ -501,15 +514,44 @@ struct CascadeSelector {
       }
     }
 
-    // if (matchedCollisions.size() < 1) {
-    //   return;
-    // } else if (matchedCollisions.size() == 1) {
-    //   for (auto const& collision : matchedCollisions) { // not really a loop, as there is only one collision
-    //   }
-    // } else if (matchedCollisions.size() > 1) {
-    //   registry.fill(HIST("MC/hSplitEvents"), matchedCollisions.size());
-    //   return;
-    // }
+    // Do the same thing, but now making sure there is at least one matched reconstructed event:
+    if (collisions.size() < 1) {
+      return;
+    } else {
+      bool evSel = false; // will be true if at least one rec. collision passes evsel
+      for (auto const& collision : collisions) {
+        // can be more than 1 rec. collisions due to event splitting
+        evSel = eventSelection(collision);
+        if (evSel) // exit loop if we find 1 rec. event that passes evsel
+          break;
+      }
+      if (evSel) {
+        // N gen events with a reconstructed event
+        registry.fill(HIST("genwithrec/hNevents"), 0);
+
+        for (auto const& mcPart : mcParticles) {
+          if (!mcPart.isPhysicalPrimary())
+            continue;
+          if (TMath::Abs(mcPart.eta()) > etaCascades)
+            continue;
+
+          switch (mcPart.pdgCode()) {
+            case 3312:
+              registry.fill(HIST("genwithrec/hXiMinus"), mcPart.pt(), mcPart.y());
+              break;
+            case -3312:
+              registry.fill(HIST("genwithrec/hXiPlus"), mcPart.pt(), mcPart.y());
+              break;
+            case 3334:
+              registry.fill(HIST("genwithrec/hOmegaMinus"), mcPart.pt(), mcPart.y());
+              break;
+            case -3334:
+              registry.fill(HIST("genwithrec/hOmegaPlus"), mcPart.pt(), mcPart.y());
+              break;
+          }
+        }
+      }
+    }
   } // processGen
 
   // wrappers for data/MC processes on reco level