From 3b1968405bddb5f047e30c09355c71664e9fc5a0 Mon Sep 17 00:00:00 2001
From: Roberta Ferioli <roberta.ferioli@edu.unito.it>
Date: Fri, 25 Jul 2025 12:07:56 +0200
Subject: [PATCH] PWGLF: Add possibility to process MC as Data, saving info

---
 PWGLF/TableProducer/Nuspex/nucleiSpectra.cxx | 55 ++++++++++++++++++++
 1 file changed, 55 insertions(+)

diff --git a/PWGLF/TableProducer/Nuspex/nucleiSpectra.cxx b/PWGLF/TableProducer/Nuspex/nucleiSpectra.cxx
index ac68cf2c9a0..f1c60d198e1 100644
--- a/PWGLF/TableProducer/Nuspex/nucleiSpectra.cxx
+++ b/PWGLF/TableProducer/Nuspex/nucleiSpectra.cxx
@@ -1052,6 +1052,61 @@ struct nucleiSpectra {
   }
 
   PROCESS_SWITCH(nucleiSpectra, processMatching, "Matching analysis", false);
+
+  void processMCasData(soa::Join<aod::Collisions, aod::EvSels, aod::McCollisionLabels> const& collisions, aod::McCollisions const& mcCollisions, soa::Join<TrackCandidates, aod::McTrackLabels> const& tracks, aod::McParticles const& particlesMC, aod::BCsWithTimestamps const&)
+  {
+    nuclei::candidates.clear();
+    std::vector<bool> goodCollisions(mcCollisions.size(), false);
+    for (auto& collision : collisions) {
+      if (!eventSelection(collision)) {
+        continue;
+      }
+      goodCollisions[collision.mcCollisionId()] = true;
+      const auto& slicedTracks = tracks.sliceBy(tracksPerCollisions, collision.globalIndex());
+      fillDataInfo(collision, slicedTracks);
+    }
+    std::vector<bool> isReconstructed(particlesMC.size(), false);
+    for (size_t i{0}; i < nuclei::candidates.size(); ++i) {
+      auto& c = nuclei::candidates[i];
+      if (c.fillTree) {
+        auto label = tracks.iteratorAt(c.globalIndex);
+        if (label.mcParticleId() < -1 || label.mcParticleId() >= particlesMC.size()) {
+          continue;
+        }
+        auto particle = particlesMC.iteratorAt(label.mcParticleId());
+        int motherPdgCode = 0;
+        float motherDecRadius = -1;
+        isReconstructed[particle.globalIndex()] = true;
+        if (particle.isPhysicalPrimary()) {
+          c.flags |= kIsPhysicalPrimary;
+          if (particle.has_mothers()) {
+            for (auto& motherparticle : particle.mothers_as<aod::McParticles>()) {
+              if (std::find(nuclei::hfMothCodes.begin(), nuclei::hfMothCodes.end(), std::abs(motherparticle.pdgCode())) != nuclei::hfMothCodes.end()) {
+                c.flags |= kIsSecondaryFromWeakDecay;
+                motherPdgCode = motherparticle.pdgCode();
+                motherDecRadius = std::hypot(particle.vx() - motherparticle.vx(), particle.vy() - motherparticle.vy());
+                break;
+              }
+            }
+          }
+        } else if (particle.has_mothers()) {
+          c.flags |= kIsSecondaryFromWeakDecay;
+          for (auto& motherparticle : particle.mothers_as<aod::McParticles>()) {
+            motherPdgCode = motherparticle.pdgCode();
+            motherDecRadius = std::hypot(particle.vx() - motherparticle.vx(), particle.vy() - motherparticle.vy());
+          }
+        } else {
+          c.flags |= kIsSecondaryFromMaterial;
+        }
+
+        isReconstructed[particle.globalIndex()] = true;
+        float absoDecL = computeAbsoDecL(particle);
+
+        nucleiTableMC(c.pt, c.eta, c.phi, c.tpcInnerParam, c.beta, c.zVertex, c.nContrib, c.DCAxy, c.DCAz, c.TPCsignal, c.ITSchi2, c.TPCchi2, c.TOFchi2, c.flags, c.TPCfindableCls, c.TPCcrossedRows, c.ITSclsMap, c.TPCnCls, c.TPCnClsShared, c.clusterSizesITS, goodCollisions[particle.mcCollisionId()], particle.pt(), particle.eta(), particle.phi(), particle.pdgCode(), motherPdgCode, motherDecRadius, absoDecL);
+      }
+    }
+  }
+  PROCESS_SWITCH(nucleiSpectra, processMCasData, "MC as data analysis", false);
 };
 
 WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)