[PWGLF] Changed delta phi computation + phi in McReco for efficiency (#14565)

Author: scannito

PWGLF/DataModel/LFPhiStrangeCorrelationTables.h

@@ -37,23 +37,23 @@ DECLARE_SOA_DYNAMIC_COLUMN(InMassRegion, inMassRegion,
                            });
 } // namespace lf_selection_phi_candidate
 
-DECLARE_SOA_TABLE(PhimesonCandidatesData, "AOD", "PHICANDIDATESDATA",
+DECLARE_SOA_TABLE(PhimesonCandidatesData, "AOD", "PHICANDDATA",
                   lf_selection_phi_candidate::CollisionId,
                   lf_selection_phi_candidate::M,
                   lf_selection_phi_candidate::Pt,
                   lf_selection_phi_candidate::Y,
                   lf_selection_phi_candidate::Phi,
                   lf_selection_phi_candidate::InMassRegion<lf_selection_phi_candidate::M>);
 
-DECLARE_SOA_TABLE(PhimesonCandidatesMcReco, "AOD", "PHICANDIDATESMCRECO",
+DECLARE_SOA_TABLE(PhimesonCandidatesMcReco, "AOD", "PHICANDMCRECO",
                   lf_selection_phi_candidate::CollisionId,
                   lf_selection_phi_candidate::M,
                   lf_selection_phi_candidate::Pt,
                   lf_selection_phi_candidate::Y,
                   lf_selection_phi_candidate::Phi,
                   lf_selection_phi_candidate::InMassRegion<lf_selection_phi_candidate::M>);
 
-DECLARE_SOA_TABLE(PhimesonCandidatesMcGen, "AOD", "PHICANDIDATESMCGEN",
+DECLARE_SOA_TABLE(PhimesonCandidatesMcGen, "AOD", "PHICANDMCGEN",
                   lf_selection_phi_candidate::CollisionId,
                   lf_selection_phi_candidate::M,
                   lf_selection_phi_candidate::Pt,

PWGLF/TableProducer/Strangeness/phiStrangeCorrelator.cxx

@@ -62,7 +62,7 @@ using namespace o2::framework::expressions;
 struct PhiMesonCandProducer {
   // Produce the table with the phi candidates information
   Produces<aod::PhimesonCandidatesData> phimesonCandidatesData;
-  // Produces<aod::PhimesonCandidatesMcReco> phimesonCandidatesMcReco;
+  Produces<aod::PhimesonCandidatesMcReco> phimesonCandidatesMcReco;
   Produces<aod::PhimesonCandidatesMcGen> phimesonCandidatesMcGen;
 
   HistogramRegistry histos{"phiCandidates", {}, OutputObjHandlingPolicy::AnalysisObject, true, true};
@@ -196,8 +196,9 @@ struct PhiMesonCandProducer {
     }
   }
 
-  PROCESS_SWITCH(PhiMesonCandProducer, processData, "Process function to select Phi meson candidates in Data", true);
+  PROCESS_SWITCH(PhiMesonCandProducer, processData, "Process function to select Phi meson candidates in Data or in McReco (w/o McTruth) analysis", true);
 
+  /*
   void processMCRecoDataLike(SimCollisions::iterator const& collision, FullMCTracks const&)
   {
     auto posThisColl = posMCTracks->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache);
@@ -226,19 +227,30 @@ struct PhiMesonCandProducer {
   }
 
   PROCESS_SWITCH(PhiMesonCandProducer, processMCRecoDataLike, "Process function to select Phi meson candidates in MCReco w/o MC truth", false);
+  */
 
-  /*void processMCReco(SimCollisions::iterator const& collision, FullMCTracks const&)
+  void processMCReco(SimCollisions::iterator const& collision, FullMCTracks const&, aod::McParticles const& mcParticles)
   {
     auto posThisColl = posMCTracks->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache);
     auto negThisColl = negMCTracks->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache);
 
     for (const auto& track1 : posThisColl) {
       if (!selectionTrackResonance(track1) || !selectionPIDKaonpTdependent(track1))
         continue;
+      if (!track1.has_mcParticle())
+        continue;
+      const auto track1McParticle = mcParticles.rawIteratorAt(track1.mcParticleId());
+      if (track1McParticle.pdgCode() != PDG_t::kKPlus || !track1McParticle.isPhysicalPrimary())
+        continue;
 
       for (const auto& track2 : negThisColl) {
         if (!selectionTrackResonance(track2) || !selectionPIDKaonpTdependent(track2))
           continue;
+        if (!track2.has_mcParticle())
+          continue;
+        const auto track2McParticle = mcParticles.rawIteratorAt(track2.mcParticleId());
+        if (track2McParticle.pdgCode() != PDG_t::kKMinus || !track2McParticle.isPhysicalPrimary())
+          continue;
 
         ROOT::Math::PxPyPzMVector recPhi = recMother(track1, track2, massKa, massKa);
 
@@ -249,12 +261,34 @@ struct PhiMesonCandProducer {
         if (std::abs(recPhi.Rapidity()) > phiConfigs.cfgYAcceptance)
           continue;
 
-        phimesonCandidatesMcReco(collision.globalIndex(), recPhi.M(), recPhi.Pt(), recPhi.Rapidity(), recPhi.Phi());
+        const auto track1mcPartMotherIndexes = track1McParticle.mothersIds();
+        const auto track2mcPartMotherIndexes = track2McParticle.mothersIds();
+
+        auto genPhiMaybe = [&]() -> std::optional<aod::McParticles::iterator> {
+          for (const auto& mother1Index : track1mcPartMotherIndexes) {
+            for (const auto& mother2Index : track2mcPartMotherIndexes) {
+              if (mother1Index != mother2Index)
+                continue;
+
+              const auto motherMcParticle = mcParticles.rawIteratorAt(mother1Index);
+              if (std::abs(motherMcParticle.pdgCode()) == o2::constants::physics::Pdg::kPhi)
+                return motherMcParticle;
+            }
+          }
+
+          return std::nullopt;
+        }();
+
+        if (!genPhiMaybe)
+          continue;
+        const auto genPhi = *genPhiMaybe;
+
+        phimesonCandidatesMcReco(collision.globalIndex(), recPhi.M(), genPhi.pt(), genPhi.y(), genPhi.phi());
       }
     }
   }
 
-  PROCESS_SWITCH(PhiMesonCandProducer, processMCReco, "Process function to select Phi meson candidates in MCReco w MC truth", false);*/
+  PROCESS_SWITCH(PhiMesonCandProducer, processMCReco, "Process function to select Phi meson candidates in MCReco w MC truth", false);
 
   void processMCGen(aod::McCollisions::iterator const& mcCollision, aod::McParticles const& mcParticles)
   {

PWGLF/Tasks/Strangeness/phiStrangeCorrelation.cxx

@@ -18,6 +18,7 @@
 #include "PWGLF/DataModel/mcCentrality.h"
 #include "PWGLF/Utils/inelGt.h"
 
+#include "Common/Core/RecoDecay.h"
 #include "Common/Core/TableHelper.h"
 #include "Common/Core/TrackSelection.h"
 #include "Common/Core/TrackSelectionDefaults.h"
@@ -282,7 +283,8 @@ struct PhiStrangenessCorrelation {
     Preslice<SimCollisions> collPerMCCollision = aod::mccollisionlabel::mcCollisionId;
     Preslice<FullMCV0s> v0PerCollision = aod::v0::collisionId;
     Preslice<FullMCTracks> trackPerCollision = aod::track::collisionId;
-    PresliceUnsorted<aod::PhimesonCandidatesData> phiCandPerCollision = aod::lf_selection_phi_candidate::collisionId;
+    // Preslice<aod::PhimesonCandidatesData> phiCandDataPerCollision = aod::lf_selection_phi_candidate::collisionId;
+    PresliceUnsorted<aod::PhimesonCandidatesMcReco> phiCandPerCollision = aod::lf_selection_phi_candidate::collisionId;
 
     // Preslice<aod::McParticles> mcPartPerMCCollision = aod::mcparticle::mcCollisionId;
   } preslices;
@@ -316,8 +318,11 @@ struct PhiStrangenessCorrelation {
 
     histos.add("phi/h3PhiData", "Invariant mass of Phi in Data", kTH3F, {binnedmultAxis, binnedpTPhiAxis, massPhiAxis});
     for (const auto& label : phiMassRegionLabels) {
-      histos.add(fmt::format("phiK0S/h5PhiK0SData2PartCorr{}", label).c_str(), "Deltay vs deltaphi for Phi and K0Short in Data", kTHnSparseF, {binnedmultAxis, binnedpTPhiAxis, binnedpTK0SAxis, deltayAxis, deltaphiAxis});
-      histos.add(fmt::format("phiPi/h5PhiPiData2PartCorr{}", label).c_str(), "Deltay vs deltaphi for Phi and Pion in Data", kTHnSparseF, {binnedmultAxis, binnedpTPhiAxis, binnedpTPiAxis, deltayAxis, deltaphiAxis});
+      histos.add(fmt::format("phiK0S/h5PhiK0SData{}", label).c_str(), "Deltay vs deltaphi for Phi and K0Short in Data", kTHnSparseF, {binnedmultAxis, binnedpTPhiAxis, binnedpTK0SAxis, deltayAxis, deltaphiAxis});
+      histos.add(fmt::format("phiPi/h5PhiPiData{}", label).c_str(), "Deltay vs deltaphi for Phi and Pion in Data", kTHnSparseF, {binnedmultAxis, binnedpTPhiAxis, binnedpTPiAxis, deltayAxis, deltaphiAxis});
+
+      histos.add(fmt::format("phiK0S/h5PhiK0SDataME{}", label).c_str(), "Deltay vs deltaphi for Phi and K0Short in Data ME", kTHnSparseF, {binnedmultAxis, binnedpTPhiAxis, binnedpTK0SAxis, deltayAxis, deltaphiAxis});
+      histos.add(fmt::format("phiPi/h5PhiPiDataME{}", label).c_str(), "Deltay vs deltaphi for Phi and Pion in Data ME", kTHnSparseF, {binnedmultAxis, binnedpTPhiAxis, binnedpTPiAxis, deltayAxis, deltaphiAxis});
     }
 
     // histos.add("phiK0S/h5PhiK0SDataNewProc", "2D Invariant mass of Phi and K0Short in Data", kTHnSparseF, {deltayAxis, binnedmultAxis, binnedpTK0SAxis, massK0SAxis, massPhiAxis});
@@ -380,6 +385,11 @@ struct PhiStrangenessCorrelation {
     return totalEfficiency <= 0.0f ? 1.0f : 1.0f / totalEfficiency;
   }
 
+  float getDeltaPhi(float phiTrigger, float phiAssociated)
+  {
+    return RecoDecay::constrainAngle(phiTrigger - phiAssociated, -o2::constants::math::PIHalf);
+  }
+
   // Single track selection for strangeness sector
   template <typename T>
   bool selectionTrackStrangeness(const T& track)
@@ -522,7 +532,7 @@ struct PhiStrangenessCorrelation {
     return true;
   }
 
-  void processPhiK0SPionDeltayDeltaphiData2D(SelCollisions::iterator const& collision, aod::PhimesonCandidatesData const& phiCandidates, FullTracks const& fullTracks, FullV0s const& V0s, V0DauTracks const&)
+  void processPhiK0SPionData(SelCollisions::iterator const& collision, aod::PhimesonCandidatesData const& phiCandidates, FullTracks const& fullTracks, FullV0s const& V0s, V0DauTracks const&)
   {
     float multiplicity = collision.centFT0M();
 
@@ -555,7 +565,7 @@ struct PhiStrangenessCorrelation {
           /*float weightPhiK0S = computeWeight(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y()),
                                              BoundEfficiencyMap(effMapK0S, multiplicity, v0.pt(), v0.yK0Short()));*/
 
-          histos.fill(HIST("phiK0S/h5PhiK0SData2PartCorr") + HIST(phiMassRegionLabels[i]), multiplicity, phiCand.pt(), v0.pt(), phiCand.y() - v0.yK0Short(), phiCand.phi() - v0.phi(), weightPhiK0S);
+          histos.fill(HIST("phiK0S/h5PhiK0SData") + HIST(phiMassRegionLabels[i]), multiplicity, phiCand.pt(), v0.pt(), phiCand.y() - v0.yK0Short(), getDeltaPhi(phiCand.phi(), v0.phi()), weightPhiK0S);
         }
 
         // Loop over all primary pion candidates
@@ -573,15 +583,62 @@ struct PhiStrangenessCorrelation {
           float weightPhiPion = computeWeight(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y()),
                                               BoundEfficiencyMap(effMapPion, multiplicity, track.pt(), track.rapidity(massPi)));*/
 
-          histos.fill(HIST("phiPi/h5PhiPiData2PartCorr") + HIST(phiMassRegionLabels[i]), multiplicity, phiCand.pt(), track.pt(), phiCand.y() - track.rapidity(massPi), phiCand.phi() - track.phi(), weightPhiPion);
+          histos.fill(HIST("phiPi/h5PhiPiData") + HIST(phiMassRegionLabels[i]), multiplicity, phiCand.pt(), track.pt(), phiCand.y() - track.rapidity(massPi), getDeltaPhi(phiCand.phi(), track.phi()), weightPhiPion);
+        }
+      });
+    }
+  }
+
+  PROCESS_SWITCH(PhiStrangenessCorrelation, processPhiK0SPionData, "Process function for Phi-K0S and Phi-Pion Deltay and Deltaphi 2D Correlations in Data", true);
+
+  /*
+  void processPhiK0SPionDataME(SelCollisions::iterator const& collision, aod::PhimesonCandidatesData const& phiCandidates, FullTracks const& fullTracks, FullV0s const& V0s, V0DauTracks const&)
+  {
+    Pair<SelCollisions, FullTracks, FullV0s, BinningTypeVertexCent> pairPhiK0S{binningOnVertexAndCent, cfgNoMixedEvents, -1, collisions, tracksV0sTuple, &cache};
+    Triple<aod::Collisions, aod::Tracks, aod::V0s, aod::Tracks, BinningType> triple{binningOnPositions, 5, -1, &cache};
+    float multiplicity = collision.centFT0M();
+
+    const std::array<std::pair<float, float>, 2> phiMassRegions = {phiConfigs.rangeMPhiSignal, phiConfigs.rangeMPhiSideband};
+
+    // Loop over all positive tracks
+    for (const auto& phiCand : phiCandidates) {
+      static_for<0, phiMassRegionLabels.size() - 1>([&](auto i_idx) {
+        constexpr unsigned int i = i_idx.value;
+
+        const auto& [minMass, maxMass] = phiMassRegions[i];
+        if (!phiCand.inMassRegion(minMass, maxMass))
+          return;
+
+        // V0 already reconstructed by the builder
+        for (const auto& v0 : V0s) {
+          // Cut on V0 dynamic columns
+          if (!selectionV0<false>(v0, collision))
+            continue;
+
+          float weightPhiK0S = computeWeight(BoundEfficiencyMap(effMaps[Phi], multiplicity, phiCand.pt(), phiCand.y()),
+                                             BoundEfficiencyMap(effMaps[K0S], multiplicity, v0.pt(), v0.yK0Short()));
+
+          histos.fill(HIST("phiK0S/h5PhiK0SDataME") + HIST(phiMassRegionLabels[i]), multiplicity, phiCand.pt(), v0.pt(), phiCand.y() - v0.yK0Short(), getDeltaPhi(phiCand.phi(), v0.phi()), weightPhiK0S);
+        }
+
+        // Loop over all primary pion candidates
+        for (const auto& track : fullTracks) {
+          if (!selectionPion(track))
+            continue;
+
+          float weightPhiPion = computeWeight(BoundEfficiencyMap(effMaps[Phi], multiplicity, phiCand.pt(), phiCand.y()),
+                                              BoundEfficiencyMap(effMaps[Pion], multiplicity, track.pt(), track.rapidity(massPi)));
+
+          histos.fill(HIST("phiPi/h5PhiPiDataME") + HIST(phiMassRegionLabels[i]), multiplicity, phiCand.pt(), track.pt(), phiCand.y() - track.rapidity(massPi), getDeltaPhi(phiCand.phi(), track.phi()), weightPhiPion);
         }
       });
     }
   }
 
-  PROCESS_SWITCH(PhiStrangenessCorrelation, processPhiK0SPionDeltayDeltaphiData2D, "Process function for Phi-K0S and Phi-Pion Deltay and Deltaphi 2D Correlations in Data", true);
+  PROCESS_SWITCH(PhiStrangenessCorrelation, processPhiK0SPionDataME, "Process function for Phi-K0S and Phi-Pion Deltay and Deltaphi 2D Correlations in Data ME", true);
+  */
 
-  void processParticleEfficiency(MCCollisions::iterator const& mcCollision, SimCollisions const& collisions, FullMCTracks const& fullMCTracks, FullMCV0s const& V0s, V0DauMCTracks const&, aod::McParticles const& mcParticles, aod::PhimesonCandidatesData const& phiCandidates)
+  void processParticleEfficiency(MCCollisions::iterator const& mcCollision, SimCollisions const& collisions, FullMCTracks const& fullMCTracks, FullMCV0s const& V0s, V0DauMCTracks const&, aod::McParticles const& mcParticles, aod::PhimesonCandidatesMcReco const& phiCandidatesMcReco)
   {
     uint16_t numberAssocColls{0};
     std::vector<float> zVtxs;
@@ -595,7 +652,7 @@ struct PhiStrangenessCorrelation {
       zVtxs.push_back(collision.posZ());
 
       if (selectionType == 0) {
-        const auto phiCandidatesThisColl = phiCandidates.sliceBy(preslices.phiCandPerCollision, collision.globalIndex());
+        const auto phiCandidatesThisColl = phiCandidatesMcReco.sliceBy(preslices.phiCandPerCollision, collision.globalIndex());
         for (const auto& phiCand : phiCandidatesThisColl) {
           histos.fill(HIST("phi/h4PhiMCReco"), collision.posZ(), mcCollision.centFT0M(), phiCand.pt(), phiCand.y());
         }