Add mother decay radius and remove track tuner

Francesco Mazzaschi · Francesco Mazzaschi · commit ebdb0e17a8e1 · 2025-06-06T17:12:18.000+02:00
diff --git a/PWGLF/DataModel/LFSlimNucleiTables.h b/PWGLF/DataModel/LFSlimNucleiTables.h
@@ -44,12 +44,13 @@ DECLARE_SOA_COLUMN(ITSclsMap, itsClsMap, uint8_t);
 DECLARE_SOA_COLUMN(TPCnCls, tpcNCls, uint8_t);
 DECLARE_SOA_COLUMN(TPCnClsShared, tpcNClsShared, uint8_t);
 DECLARE_SOA_COLUMN(ITSclusterSizes, itsClusterSizes, uint32_t);
+DECLARE_SOA_COLUMN(SurvivedEventSelection, survivedEventSelection, bool);
 DECLARE_SOA_COLUMN(gPt, genPt, float);
 DECLARE_SOA_COLUMN(gEta, genEta, float);
 DECLARE_SOA_COLUMN(gPhi, genPhi, float);
 DECLARE_SOA_COLUMN(PDGcode, pdgCode, int);
 DECLARE_SOA_COLUMN(MotherPDGcode, MotherpdgCode, int);
-DECLARE_SOA_COLUMN(SurvivedEventSelection, survivedEventSelection, bool);
+DECLARE_SOA_COLUMN(MotherDecRad, motherDecRad, float);
 DECLARE_SOA_COLUMN(AbsoDecL, absoDecL, float);
 
 } // namespace NucleiTableNS
@@ -151,12 +152,13 @@ DECLARE_SOA_TABLE(NucleiTableMC, "AOD", "NUCLEITABLEMC",
                   NucleiTableNS::TPCnCls,
                   NucleiTableNS::TPCnClsShared,
                   NucleiTableNS::ITSclusterSizes,
+                  NucleiTableNS::SurvivedEventSelection,
                   NucleiTableNS::gPt,
                   NucleiTableNS::gEta,
                   NucleiTableNS::gPhi,
                   NucleiTableNS::PDGcode,
                   NucleiTableNS::MotherPDGcode,
-                  NucleiTableNS::SurvivedEventSelection,
+                  NucleiTableNS::MotherDecRad,
                   NucleiTableNS::AbsoDecL);
 
 DECLARE_SOA_TABLE(NucleiPairTable, "AOD", "NUCLEIPAIRTABLE",
diff --git a/PWGLF/TableProducer/Nuspex/nucleiSpectra.cxx b/PWGLF/TableProducer/Nuspex/nucleiSpectra.cxx
@@ -41,7 +41,6 @@
 #include "Common/TableProducer/PID/pidTOFBase.h"
 #include "Common/Core/EventPlaneHelper.h"
 #include "Common/DataModel/Qvectors.h"
-#include "Common/Tools/TrackTuner.h"
 #include "Common/Core/RecoDecay.h"
 
 #include "DataFormatsParameters/GRPMagField.h"
@@ -273,7 +272,6 @@ struct nucleiSpectra {
   Service<o2::ccdb::BasicCCDBManager> ccdb;
   Zorro zorro;
   OutputObj<ZorroSummary> zorroSummary{"zorroSummary"};
-  TrackTuner trackTunerObj;
 
   Configurable<bool> cfgCompensatePIDinTracking{"cfgCompensatePIDinTracking", false, "If true, divide tpcInnerParam by the electric charge"};
 
@@ -282,8 +280,8 @@ struct nucleiSpectra {
   Configurable<float> cfgCutVertex{"cfgCutVertex", 10.0f, "Accepted z-vertex range"};
   Configurable<float> cfgCutEta{"cfgCutEta", 0.8f, "Eta range for tracks"};
   Configurable<float> cfgCutTpcMom{"cfgCutTpcMom", 0.2f, "Minimum TPC momentum for tracks"};
-  Configurable<float> cfgCutRapidityMin{"cfgCutRapidityMin", -0.5, "Minimum rapidity for tracks"};
-  Configurable<float> cfgCutRapidityMax{"cfgCutRapidityMax", 0.5, "Maximum rapidity for tracks"};
+  Configurable<float> cfgCutRapidityMin{"cfgCutRapidityMin", -1., "Minimum rapidity for tracks"};
+  Configurable<float> cfgCutRapidityMax{"cfgCutRapidityMax", 1., "Maximum rapidity for tracks"};
   Configurable<bool> cfgCutOnReconstructedRapidity{"cfgCutOnReconstructedRapidity", false, "Cut on reconstructed rapidity"};
   Configurable<float> cfgCutNclusITS{"cfgCutNclusITS", 5, "Minimum number of ITS clusters"};
   Configurable<float> cfgCutNclusTPC{"cfgCutNclusTPC", 70, "Minimum number of TPC clusters"};
@@ -327,9 +325,7 @@ struct nucleiSpectra {
 
   Configurable<bool> cfgSkimmedProcessing{"cfgSkimmedProcessing", false, "Skimmed dataset processing"};
 
-  // configurables for track tuner
-  Configurable<bool> cfgUseTrackTuner{"cfgUseTrackTuner", false, "Apply track tuner corrections to MC tracks"};
-  Configurable<std::string> cfgTrackTunerParams{"cfgTrackTunerParams", "debugInfo=0|updateTrackDCAs=1|updateTrackCovMat=1|updateCurvature=0|updateCurvatureIU=0|updatePulls=1|isInputFileFromCCDB=1|pathInputFile=Users/m/mfaggin/test/inputsTrackTuner/pp2023/smoothHighPtMC|nameInputFile=trackTuner_DataLHC23fPass1_McLHC23k4b_run535085.root|pathFileQoverPt=Users/h/hsharma/qOverPtGraphs|nameFileQoverPt=D0sigma_Data_removal_itstps_MC_LHC22b1b.root|usePvRefitCorrections=0|qOverPtMC=-1.|qOverPtData=-1.", "TrackTuner parameter initialization (format: <name>=<value>|<name>=<value>)"};
+ 
   // running variables for track tuner
   o2::dataformats::DCA mDcaInfoCov;
   o2::track::TrackParametrizationWithError<float> mTrackParCov;
@@ -516,7 +512,7 @@ struct nucleiSpectra {
     spectra.add("hTpcSignalData", "Specific energy loss", HistType::kTH2F, {{600, -6., 6., "#it{p} (GeV/#it{c})"}, {1400, 0, 1400, "d#it{E} / d#it{X} (a. u.)"}});
     spectra.add("hTpcSignalDataSelected", "Specific energy loss for selected particles", HistType::kTH2F, {{600, -6., 6., "#it{p} (GeV/#it{c})"}, {1400, 0, 1400, "d#it{E} / d#it{X} (a. u.)"}});
     spectra.add("hTofSignalData", "TOF beta", HistType::kTH2F, {{500, 0., 5., "#it{p} (GeV/#it{c})"}, {750, 0, 1.5, "TOF #beta"}});
-    spectra.add("hMCDCAxyVsDecRadius", "DCA xy vs DecayRadius", HistType::kTH2F, {{20, 0., 0.2, "Decay radius"}, {150, -0.05, 0.05, "DCA_{xy} (cm)"}});
+
     for (int iC{0}; iC < 2; ++iC) {
       nuclei::hGloTOFtracks[iC] = spectra.add<TH2>(fmt::format("hTPCTOFtracks{}", nuclei::matter[iC]).data(), fmt::format("Global vs TOF matched {} tracks in a collision", nuclei::chargeLabelNames[iC]).data(), HistType::kTH2D, {{300, -0.5, 300.5, "Number of global tracks"}, {300, -0.5, 300.5, "Number of TOF matched tracks"}});
 
@@ -562,12 +558,6 @@ struct nucleiSpectra {
     }
 
     nuclei::lut = o2::base::MatLayerCylSet::rectifyPtrFromFile(ccdb->get<o2::base::MatLayerCylSet>("GLO/Param/MatLUT"));
-    // TrackTuner initialization
-    if (cfgUseTrackTuner) {
-      std::string outputStringParams = trackTunerObj.configParams(cfgTrackTunerParams);
-      spectra.add("hTrackTunedTracks", outputStringParams.c_str(), HistType::kTH1F, {{1, 0.5, 1.5, ""}});
-      trackTunerObj.getDcaGraphs();
-    }
   }
 
   template <typename Tcoll>
@@ -665,20 +655,6 @@ struct nucleiSpectra {
       mDcaInfoCov.set(999, 999, 999, 999, 999);
       setTrackParCov(track, mTrackParCov);
       mTrackParCov.setPID(track.pidForTracking());
-      if constexpr (
-        requires {
-          track.has_mcParticle();
-        }) {
-        if (cfgUseTrackTuner) {
-          bool hasMcParticle = track.has_mcParticle();
-          if (hasMcParticle) {
-            spectra.get<TH1>(HIST("hTrackTunedTracks"))->Fill(1); // all tracks
-            auto mcParticle = track.mcParticle();
-            trackTunerObj.tuneTrackParams(mcParticle, mTrackParCov, matCorr, &mDcaInfoCov, spectra.get<TH1>(HIST("hTrackTunedTracks")));
-          }
-        }
-      }
-
       std::array<float, 2> dcaInfo;
       o2::base::Propagator::Instance()->propagateToDCA(collVtx, mTrackParCov, mBz, 2.f, static_cast<o2::base::Propagator::MatCorrType>(cfgMaterialCorrection.value), &dcaInfo);
 
@@ -957,24 +933,6 @@ struct nucleiSpectra {
           if (!c.correctPV) {
             c.flags |= kIsAmbiguous;
           }
-          if (!particle.isPhysicalPrimary()) {
-            c.isSecondary = true;
-            if (particle.getProcess() == 4) {
-              c.fromWeakDecay = true;
-            }
-          } else {
-            // if the particle has a hf mother it is flagged as secondary
-            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.isSecondary = true;
-                  c.fromWeakDecay = true;
-                  break;
-                }
-              }
-            }
-          }
-
           if (c.fillDCAHist && cfgDCAHists->get(iS, c.pt < 0)) {
             nuclei::hDCAHists[c.pt < 0][iS]->Fill(std::abs(c.pt), c.DCAxy, c.DCAz, c.nSigmaTPC[iS], c.tofMasses[iS], c.ITSnCls, c.TPCnCls, c.correctPV, c.isSecondary, c.fromWeakDecay);
           }
@@ -983,31 +941,38 @@ struct nucleiSpectra {
       if (!storeIt) {
         continue;
       }
-      int MotherpdgCode = 0;
+      if (particle.y() < cfgCutRapidityMin || particle.y() > cfgCutRapidityMax) {
+        continue;
+      }
+
+      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();
-              float decRadius = std::hypot(particle.vx() - motherparticle.vx(), particle.vy() - motherparticle.vy());
-              spectra.fill(HIST("hMCDCAxyVsDecRadius"), decRadius, c.DCAxy);
+              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();
+          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.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, particle.pt(), particle.eta(), particle.phi(), particle.pdgCode(), MotherpdgCode, goodCollisions[particle.mcCollisionId()], absoDecL);
+      nucleiTableMC(c.pt, c.eta, c.phi, c.tpcInnerParam, c.beta, c.zVertex, 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);
     }
 
     int index{0};
@@ -1023,6 +988,7 @@ struct nucleiSpectra {
 
         uint16_t flags = 0;
         int motherPdgCode = 0;
+        float motherDecRadius = -1;
         if (particle.isPhysicalPrimary()) {
           flags |= kIsPhysicalPrimary;
           nuclei::hGenNuclei[iS][particle.pdgCode() < 0]->Fill(1., particle.pt());
@@ -1031,6 +997,7 @@ struct nucleiSpectra {
               if (std::find(nuclei::hfMothCodes.begin(), nuclei::hfMothCodes.end(), std::abs(motherparticle.pdgCode())) != nuclei::hfMothCodes.end()) {
                 flags |= kIsSecondaryFromWeakDecay;
                 motherPdgCode = motherparticle.pdgCode();
+                motherDecRadius = std::hypot(particle.vx() - motherparticle.vx(), particle.vy() - motherparticle.vy());
                 break;
               }
             }
@@ -1039,6 +1006,7 @@ struct nucleiSpectra {
           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 {
           flags |= kIsSecondaryFromMaterial;
@@ -1047,7 +1015,7 @@ struct nucleiSpectra {
         if (!isReconstructed[index] && (cfgTreeConfig->get(iS, 0u) || cfgTreeConfig->get(iS, 1u))) {
           float absDecL = computeAbsoDecL(particle);
 
-          nucleiTableMC(999., 999., 999., 0., 0., 999., 999., 999., -1, -1, -1, -1, flags, 0, 0, 0, 0, 0, 0, particle.pt(), particle.eta(), particle.phi(), particle.pdgCode(), motherPdgCode, goodCollisions[particle.mcCollisionId()], absDecL);
+          nucleiTableMC(999., 999., 999., 0., 0., 999., 999., 999., -1, -1, -1, -1, flags, 0, 0, 0, 0, 0, 0, goodCollisions[particle.mcCollisionId()], particle.pt(), particle.eta(), particle.phi(), particle.pdgCode(), motherPdgCode, motherDecRadius, absDecL);
         }
         break;
       }
