[PWGLF] check the production vertex of nuclei from material (#16793)

Author: GiorgioAlbertoLucia

PWGLF/DataModel/LFSlimNucleiTables.h

@@ -56,8 +56,11 @@ DECLARE_SOA_COLUMN(AbsoDecL, absoDecL, float);
 DECLARE_SOA_COLUMN(McProcess, mcProcess, uint64_t);
 DECLARE_SOA_COLUMN(gEventMask, genEventMask, uint8_t);
 
-DECLARE_SOA_COLUMN(NsigmaTpc, nsigmaTpc, uint8_t);
-DECLARE_SOA_COLUMN(NsigmaTof, nsigmaTof, uint8_t);
+DECLARE_SOA_COLUMN(NsigmaTpc, nsigmaTpc, float);
+DECLARE_SOA_COLUMN(NsigmaTof, nsigmaTof, float);
+
+DECLARE_SOA_COLUMN(Vx, vx, float);
+DECLARE_SOA_COLUMN(Vy, vy, float);
 
 } // namespace NucleiTableNS
 
@@ -245,6 +248,10 @@ DECLARE_SOA_TABLE(NucleiTableMCExtension, "AOD", "NUCTABLEMCSEL",
 DECLARE_SOA_TABLE(NucleiTableExt, "AOD", "NUCLEITABLEEXT",
                   NucleiTableNS::NsigmaTpc,
                   NucleiTableNS::NsigmaTof);
+// Extended table for studies on nuclei from material
+DECLARE_SOA_TABLE(NucleiTableMat, "AOD", "NUCLEITABLEMAT",
+                  NucleiTableNS::Vx,
+                  NucleiTableNS::Vy);
 
 } // namespace o2::aod
 

PWGLF/TableProducer/QC/nucleiQC.cxx

@@ -162,6 +162,7 @@ struct nucleiQC {
   std::array<bool, nuclei::Species::kNspecies> mFillSpecies{false};
   Produces<aod::NucleiTableRed> mNucleiTableRed;
   Produces<aod::NucleiTableExt> mNucleiTableExt;
+  Produces<aod::NucleiTableMat> mNucleiTableMat;
 
   std::vector<nuclei::SlimCandidate> mNucleiCandidates;
   std::vector<int> mFilledMcParticleIds;
@@ -209,6 +210,7 @@ struct nucleiQC {
 
       nuclei::createHistogramRegistryNucleus<kSpeciesCt>(mHistograms);
       mHistograms.add(fmt::format("{}/hTrackQuality", nuclei::cNames[kSpeciesRt]).c_str(), (fmt::format("{} track quality;", nuclei::cNames[kSpeciesRt]) + std::string("#it{p}_{T} / #it{Z} (GeV/#it{c}); Selection step; Counts")).c_str(), o2::framework::HistType::kTH2D, {{400, -10.0f, 10.0f}, {trackQuality::kNtrackQuality, -0.5f, static_cast<float>(trackQuality::kNtrackQuality) - 0.5f}});
+      mHistograms.add(fmt::format("{}/h2Productionvertex", nuclei::cNames[kSpeciesRt]).c_str(), (fmt::format("{} production vertex;", nuclei::cNames[kSpeciesRt]) + std::string("#it{x} (cm); #it{y} (cm); Counts")).c_str(), o2::framework::HistType::kTH2D, {{400, -100.0f, 100.0f}, {400, -100.0f, 100.0f}});
       for (size_t iSel = 0; iSel < trackQuality::kNtrackQuality; iSel++) {
         mHistograms.get<TH2>(HIST(nuclei::cNames[kSpeciesRt]) + HIST("/hTrackQuality"))->GetYaxis()->SetBinLabel(iSel + 1, trackQualityLabels[iSel].c_str());
       }
@@ -352,7 +354,7 @@ struct nucleiQC {
     return true;
   }
 
-  template <typename Tparticle>
+  template <int iSpecies, typename Tparticle>
   void fillNucleusFlagsPdgsMc(const Tparticle& particle, nuclei::SlimCandidate& candidate)
   {
     candidate.pdgCode = particle.pdgCode();
@@ -383,6 +385,9 @@ struct nucleiQC {
       candidate.flags |= nuclei::QcFlags::kQcIsSecondaryFromWeakDecay;
     } else {
       candidate.flags |= nuclei::QcFlags::kQcIsSecondaryFromMaterial;
+      mHistograms.fill(HIST(nuclei::cNames[iSpecies]) + HIST("/h2Productionvertex"), particle.vx(), particle.vy());
+      candidate.vx = particle.vx();
+      candidate.vy = particle.vy();
     }
   }
 
@@ -443,15 +448,22 @@ struct nucleiQC {
                                        .centrality = nuclei::getCentrality(collision, cfgCentralityEstimator, mHistFailCentrality),
                                        .mcProcess = TMCProcess::kPNoProcess,
                                        .nsigmaTpc = mPidManagers[iSpecies].getNSigmaTPC(track),
-                                       .nsigmaTof = mPidManagers[iSpecies].getNSigmaTOF(track)};
+                                       .nsigmaTof = mPidManagers[iSpecies].getNSigmaTOF(track),
+                                       .vx = -999.f,
+                                       .vy = -999.f};
 
     fillNucleusFlagsPdgs(collision, track, candidate);
 
     if constexpr (isMc) {
       if (track.has_mcParticle()) {
 
         const auto& particle = track.mcParticle();
-        fillNucleusFlagsPdgsMc(particle, candidate);
+        static_for<0, nuclei::kNspecies - 1>([&](auto iSpeciesCtV) {
+          constexpr int kSpeciesCt = decltype(iSpeciesCtV)::value;
+          if (std::abs(particle.pdgCode()) == nuclei::pdgCodes[kSpeciesCt]) {
+            fillNucleusFlagsPdgsMc<kSpeciesCt>(particle, candidate);
+          }
+        });
         fillNucleusGeneratedVariables(particle, candidate);
       }
     }
@@ -527,6 +539,9 @@ struct nucleiQC {
     mNucleiTableExt(
       candidate.nsigmaTpc,
       candidate.nsigmaTof);
+    mNucleiTableMat(
+      candidate.vx,
+      candidate.vy);
   }
 
   void processMc(const Collisions& collisions, const TrackCandidatesMC& tracks, const aod::BCsWithTimestamps&, const aod::McParticles& mcParticles, const aod::McCollisions& /*mcCollisions*/)
@@ -647,7 +662,12 @@ struct nucleiQC {
       const auto& centralityIt = mcCollisionIdToCentrality.find(particle.mcCollisionId());
       candidate.centrality = centralityIt != mcCollisionIdToCentrality.end() ? centralityIt->second : -1.f;
       fillCollisionFlag(particle, candidate, reconstructedCollisions);
-      fillNucleusFlagsPdgsMc(particle, candidate);
+      static_for<0, nuclei::kNspecies - 1>([&](auto iSpeciesCtV) {
+        constexpr int kSpeciesCt = decltype(iSpeciesCtV)::value;
+        if (std::abs(particle.pdgCode()) == nuclei::pdgCodes[kSpeciesCt]) {
+          fillNucleusFlagsPdgsMc<kSpeciesCt>(particle, candidate);
+        }
+      });
       fillNucleusGeneratedVariables(particle, candidate);
 
       writeCandidate(candidate);

PWGLF/Utils/nucleiUtils.h

@@ -122,6 +122,8 @@ struct SlimCandidate {
   uint64_t mcProcess = TMCProcess::kPNoProcess;
   float nsigmaTpc = -999.f;
   float nsigmaTof = -999.f;
+  float vx = -999.f; // production vertex x coordinate
+  float vy = -999.f;
 };
 
 enum Species {