[PWGLF] Select primary and from HF for jet clustering,improved perp cone calculation

Author: alcaliva

PWGLF/Tasks/Nuspex/antinucleiInJets.cxx

@@ -12,7 +12,7 @@
 /// \file antinucleiInJets.cxx
 ///
 /// \brief task for analysis of antinuclei in jets using Fastjet
-/// \author Alberto Caliva (alberto.caliva@cern.ch), Chiara Pinto (chiara.pinto@cern.ch)
+/// \author Alberto Caliva (alberto.caliva@cern.ch), Chiara Pinto (chiara.pinto@cern.ch), Francesca Casillo (francesca.casillo@cern.ch)
 /// \since February 13, 2025
 
 #include "PWGJE/Core/JetBkgSubUtils.h"
@@ -111,14 +111,14 @@ struct AntinucleiInJets {
   Configurable<double> ptLeadingMin{"ptLeadingMin", 5.0, "pt Leading Min"};
   Configurable<double> rJet{"rJet", 0.3, "Jet resolution parameter R"};
   Configurable<double> zVtx{"zVtx", 10.0, "Maximum zVertex"};
-  Configurable<bool> applyAreaCut{"applyAreaCut", false, "apply area cut"};
-  Configurable<double> maxNormalizedJetArea{"maxNormalizedJetArea", 0.6, "area cut"};
+  Configurable<bool> applyAreaCut{"applyAreaCut", true, "apply area cut"};
+  Configurable<double> maxNormalizedJetArea{"maxNormalizedJetArea", 1.0, "area cut"};
   Configurable<double> deltaEtaEdge{"deltaEtaEdge", 0.05, "eta gap from the edge"};
   Configurable<int> nSyst{"nSyst", 50, "number of systematic variations"};
 
   // Track quality, kinematic, and PID selection parameters
   Configurable<bool> requirePvContributor{"requirePvContributor", false, "require that the track is a PV contributor"};
-  Configurable<bool> applyItsPid{"applyItsPid", true, "apply ITS PID"};
+  Configurable<bool> applyItsPid{"applyItsPid", false, "apply ITS PID"};
   Configurable<int> minItsNclusters{"minItsNclusters", 5, "minimum number of ITS clusters"};
   Configurable<int> minTpcNcrossedRows{"minTpcNcrossedRows", 100, "minimum number of TPC crossed pad rows"};
   Configurable<double> minChiSquareTpc{"minChiSquareTpc", 0.0, "minimum TPC chi^2/Ncls"};
@@ -138,18 +138,21 @@ struct AntinucleiInJets {
   Configurable<double> ptMaxItsPidHel{"ptMaxItsPidHel", 1.0, "maximum pt for ITS PID for helium"};
   Configurable<double> nSigmaItsMin{"nSigmaItsMin", -3.0, "nSigmaITS min"};
   Configurable<double> nSigmaItsMax{"nSigmaItsMax", +3.0, "nSigmaITS max"};
-  Configurable<bool> setMCDefaultItsParams{"setMCDefaultItsParams", false, "set MC default parameters"};
+  Configurable<bool> setMCDefaultItsParams{"setMCDefaultItsParams", true, "set MC default parameters"};
   Configurable<bool> cfgCompensatePIDinTracking{"cfgCompensatePIDinTracking", false, "If true, divide tpcInnerParam by the electric charge"};
   Configurable<std::array<double, 5>> cfgBetheBlochParams{"cfgBetheBlochParams", {0.6539, 1.591, 0.8225, 2.363, 0.09}, "TPC Bethe-Bloch parameterisation for He3"};
 
   // Configuration parameters for CCDB access and reweighting input files
+  Configurable<bool> applyReweighting{"applyReweighting", true, "enable reweighting for efficiency"};
   Configurable<std::string> urlToCcdb{"urlToCcdb", "http://alice-ccdb.cern.ch", "url of the personal ccdb"};
-  Configurable<std::string> pathToFile{"pathToFile", "Users/a/alcaliva/PrimaryFractionAntip", "path to file with reweighting"};
+  Configurable<std::string> pathToFile{"pathToFile", "Users/a/alcaliva/reweightingHistogramsAntipInJet", "path to file"};
   Configurable<std::string> weightsProton{"weightsProton", "", "weightsProton"};
   Configurable<std::string> weightsLambda{"weightsLambda", "", "weightsLambda"};
   Configurable<std::string> weightsSigma{"weightsSigma", "", "weightsSigma"};
   Configurable<std::string> weightsXi{"weightsXi", "", "weightsXi"};
   Configurable<std::string> weightsOmega{"weightsOmega", "", "weightsOmega"};
+  Configurable<std::string> weightsJet{"weightsJet", "", "weightsJet"};
+  Configurable<std::string> weightsUe{"weightsUe", "", "weightsUe"};
 
   // Number of events
   Configurable<int> shrinkInterval{"shrinkInterval", 1000, "variable that controls how often shrinking happens"};
@@ -160,6 +163,8 @@ struct AntinucleiInJets {
   TH1F* primaryAntiSigma;
   TH1F* primaryAntiXi;
   TH1F* primaryAntiOmega;
+  TH1F* antiprotonsInsideJets;
+  TH1F* antiprotonsPerpCone;
 
   // CCDB manager service for accessing condition data
   Service<o2::ccdb::BasicCCDBManager> ccdb;
@@ -205,7 +210,7 @@ struct AntinucleiInJets {
       ccdb->setLocalObjectValidityChecking();
       ccdb->setCreatedNotAfter(std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count());
       ccdb->setFatalWhenNull(false);
-      getReweightingHistograms(ccdb, TString(pathToFile), TString(weightsProton), TString(weightsLambda), TString(weightsSigma), TString(weightsXi), TString(weightsOmega));
+      getReweightingHistograms(ccdb, TString(pathToFile), TString(weightsProton), TString(weightsLambda), TString(weightsSigma), TString(weightsXi), TString(weightsOmega), TString(weightsJet), TString(weightsUe));
     }
 
     // Binning
@@ -295,6 +300,9 @@ struct AntinucleiInJets {
       registryMC.add("genEvents", "number of generated events in mc", HistType::kTH1F, {{10, 0, 10, "counter"}});
       registryMC.add("genJets", "number of generated jets", HistType::kTH1F, {{10, 0, 10, "counter"}});
 
+      // Size of particle array
+      registryMC.add("sizeParticleArray", "number of particles", HistType::kTH1F, {{1000, 0, 1000, "#it{N}_{part}"}});
+
       // Generated spectra of antiprotons
       registryMC.add("antiproton_gen_jet", "antiproton_gen_jet", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}});
       registryMC.add("antiproton_gen_ue", "antiproton_gen_ue", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}});
@@ -485,22 +493,30 @@ struct AntinucleiInJets {
     }
   }
 
-  void getReweightingHistograms(o2::framework::Service<o2::ccdb::BasicCCDBManager> const& ccdbObj, TString filepath, TString antip, TString antilambda, TString antisigma, TString antixi, TString antiomega)
+  void getReweightingHistograms(o2::framework::Service<o2::ccdb::BasicCCDBManager> const& ccdbObj, TString filepath, TString antip, TString antilambda, TString antisigma, TString antixi, TString antiomega, TString jet, TString ue)
   {
     TList* list = ccdbObj->get<TList>(filepath.Data());
     if (!list) {
       LOGP(error, "Could not retrieve the list from CCDB");
       return;
     }
 
+    // Get reweighting histograms for primary fraction
     primaryAntiprotons = static_cast<TH1F*>(list->FindObject(antip));
     primaryAntiLambda = static_cast<TH1F*>(list->FindObject(antilambda));
     primaryAntiSigma = static_cast<TH1F*>(list->FindObject(antisigma));
     primaryAntiXi = static_cast<TH1F*>(list->FindObject(antixi));
     primaryAntiOmega = static_cast<TH1F*>(list->FindObject(antiomega));
 
     if (!primaryAntiprotons || !primaryAntiSigma || !primaryAntiLambda || !primaryAntiXi || !primaryAntiOmega) {
-      LOGP(error, "Missing one or more reweighting histograms in CCDB list");
+      LOGP(error, "Missing one or more reweighting histograms for primary fraction in CCDB list");
+    }
+
+    // Get reweighting histograms for efficiency
+    antiprotonsInsideJets = static_cast<TH1F*>(list->FindObject(jet));
+    antiprotonsPerpCone = static_cast<TH1F*>(list->FindObject(ue));
+    if (!antiprotonsInsideJets || !antiprotonsPerpCone) {
+      LOGP(error, "Missing one or more reweighting histograms for efficiency in CCDB list");
     }
 
     LOGP(info, "Successfully loaded reweighting histograms from CCDB path");
@@ -533,6 +549,27 @@ struct AntinucleiInJets {
     return current;
   }
 
+  // Check if particle is a physical primary or a decay product of a heavy-flavor hadron
+  bool isPhysicalPrimaryOrFromHF(aod::McParticle const& particle, aod::McParticles const& mcParticles)
+  {
+    // Keep only pi, K, p, e, mu
+    int pdg = std::abs(particle.pdgCode());
+    if (!(pdg == 211 || pdg == 321 || pdg == 2212 || pdg == 11 || pdg == 13))
+      return false;
+
+    // Check if particle is from heavy-flavor decay
+    bool fromHF = false;
+    if (particle.has_mothers()) {
+      auto mother = mcParticles.iteratorAt(particle.mothersIds()[0]);
+      int motherPdg = std::abs(mother.pdgCode());
+      fromHF = (motherPdg / 100 == 4 || motherPdg / 100 == 5 || motherPdg / 1000 == 4 || motherPdg / 1000 == 5);
+    }
+
+    // Select only physical primary particles or from heavy-flavor
+    return (particle.isPhysicalPrimary() || fromHF);
+  }
+
+  /*
   // Compute two unit vectors perpendicular to p
   void getPerpendicularAxis(const TVector3& p, TVector3& u, double sign)
   {
@@ -584,6 +621,70 @@ struct AntinucleiInJets {
     double uy = (-pz2 - px * ux) / py;
     u.SetXYZ(ux, uy, pz);
   }
+  */
+
+  // Compute two transverse directions orthogonal to vector p
+  void getPerpendicularDirections(const TVector3& p, TVector3& u1, TVector3& u2)
+  {
+    // Get momentum components
+    double px = p.X();
+    double py = p.Y();
+    double pz = p.Z();
+
+    // Precompute squared terms
+    double px2 = px * px;
+    double py2 = py * py;
+    double pz2 = pz * pz;
+    double pz4 = pz2 * pz2;
+
+    // Case 1: vector along z-axis -> undefined perpendiculars
+    if (px == 0 && py == 0) {
+      u1.SetXYZ(0, 0, 0);
+      u2.SetXYZ(0, 0, 0);
+      return;
+    }
+
+    // Case 2: px = 0 -> avoid division by zero
+    if (px == 0 && py != 0) {
+      double ux = std::sqrt(py2 - pz4 / py2);
+      double uy = -pz2 / py;
+      u1.SetXYZ(ux, uy, pz);
+      u2.SetXYZ(-ux, uy, pz);
+      return;
+    }
+
+    // Case 3: py = 0 -> avoid division by zero
+    if (py == 0 && px != 0) {
+      double ux = -pz2 / px;
+      double uy = std::sqrt(px2 - pz4 / px2);
+      u1.SetXYZ(ux, uy, pz);
+      u2.SetXYZ(ux, -uy, pz);
+      return;
+    }
+
+    // General case: solve quadratic for perpendicular vectors
+    double a = px2 + py2;
+    double b = 2.0 * px * pz2;
+    double c = pz4 - py2 * py2 - px2 * py2;
+    double delta = b * b - 4.0 * a * c;
+
+    // Invalid or degenerate solutions
+    if (delta < 0 || a == 0) {
+      u1.SetXYZ(0, 0, 0);
+      u2.SetXYZ(0, 0, 0);
+      return;
+    }
+
+    // Solution 1
+    double u1x = (-b + std::sqrt(delta)) / (2.0 * a);
+    double u1y = (-pz2 - px * u1x) / py;
+    u1.SetXYZ(u1x, u1y, pz);
+
+    // Solution 2
+    double u2x = (-b - std::sqrt(delta)) / (2.0 * a);
+    double u2y = (-pz2 - px * u2x) / py;
+    u2.SetXYZ(u2x, u2y, pz);
+  }
 
   // Compute delta phi
   double getDeltaPhi(double a1, double a2)
@@ -917,10 +1018,11 @@ struct AntinucleiInJets {
       // Perpendicular cones
       double coneRadius = std::sqrt(jet.area() / PI);
       TVector3 jetAxis(jet.px(), jet.py(), jet.pz());
-      TVector3 ueAxis1(0, 0, 0);
-      TVector3 ueAxis2(0, 0, 0);
-      getPerpendicularAxis(jetAxis, ueAxis1, +1);
-      getPerpendicularAxis(jetAxis, ueAxis2, -1);
+      TVector3 ueAxis1(0, 0, 0), ueAxis2(0, 0, 0);
+      getPerpendicularDirections(jetAxis, ueAxis1, ueAxis2);
+      if (ueAxis1.Mag() == 0 || ueAxis2.Mag() == 0) {
+        continue;
+      }
 
       // Fill histogram with jet effective area / piR^2
       registryData.fill(HIST("jetEffectiveAreaOverPiR2"), jet.area() / (PI * rJet * rJet));
@@ -1274,10 +1376,11 @@ struct AntinucleiInJets {
       std::vector<fastjet::PseudoJet> jetConstituents = jet.constituents();
       TVector3 jetAxis(jet.px(), jet.py(), jet.pz());
       double coneRadius = std::sqrt(jet.area() / PI);
-      TVector3 ueAxis1(0, 0, 0);
-      TVector3 ueAxis2(0, 0, 0);
-      getPerpendicularAxis(jetAxis, ueAxis1, +1);
-      getPerpendicularAxis(jetAxis, ueAxis2, -1);
+      TVector3 ueAxis1(0, 0, 0), ueAxis2(0, 0, 0);
+      getPerpendicularDirections(jetAxis, ueAxis1, ueAxis2);
+      if (ueAxis1.Mag() == 0 || ueAxis2.Mag() == 0) {
+        continue;
+      }
 
       registryQC.fill(HIST("NchJetCone"), static_cast<int>(jetConstituents.size()));
 
@@ -1656,9 +1759,11 @@ struct AntinucleiInJets {
       // Loop over MC particles
       for (const auto& particle : mcParticles) {
 
-        // Select physical primaries within acceptance
-        if (!particle.isPhysicalPrimary())
+        // Select physical primary particles or HF decay products
+        if (!isPhysicalPrimaryOrFromHF(particle,mcParticles))
           continue;
+
+        // Select particles within acceptance
         static constexpr double MinPtParticle = 0.1;
         if (particle.eta() < minEta || particle.eta() > maxEta || particle.pt() < MinPtParticle)
           continue;
@@ -1681,6 +1786,9 @@ struct AntinucleiInJets {
         continue;
       registryMC.fill(HIST("genEvents"), 2.5);
 
+      // Size of particle array
+      registryMC.fill(HIST("sizeParticleArray"), fjParticles.size());
+
       // Cluster MC particles into jets using anti-kt algorithm
       fastjet::ClusterSequenceArea cs(fjParticles, jetDef, areaDef);
       std::vector<fastjet::PseudoJet> jets = fastjet::sorted_by_pt(cs.inclusive_jets());
@@ -1721,19 +1829,28 @@ struct AntinucleiInJets {
           // Fill normalization histogram
           registryMC.fill(HIST("antiproton_deltay_deltaphi_jet"), particle.eta() - jet.eta(), getDeltaPhi(particle.phi(), jet.phi()));
 
+          // Calculate weight
+          double weightJet(1.0);
+          if (applyReweighting && particle.pt() < antiprotonsInsideJets->GetXaxis()->GetXmax()) {
+            int ipt = antiprotonsInsideJets->FindBin(particle.pt());
+            weightJet = antiprotonsInsideJets->GetBinContent(ipt);
+          }
+
           // Fill histogram for generated antiprotons
-          registryMC.fill(HIST("antiproton_gen_jet"), particle.pt());
+          registryMC.fill(HIST("antiproton_gen_jet"), particle.pt(), weightJet);
 
           // Fill 2d (pt,eta) distribution of antiprotons
-          registryMC.fill(HIST("antiproton_eta_pt_jet"), particle.pt(), particle.eta());
+          registryMC.fill(HIST("antiproton_eta_pt_jet"), particle.pt(), particle.eta(), weightJet);
         }
 
         // Set up two perpendicular cone axes for underlying event estimation
         TVector3 jetAxis(jet.px(), jet.py(), jet.pz());
         double coneRadius = std::sqrt(jet.area() / PI);
         TVector3 ueAxis1(0, 0, 0), ueAxis2(0, 0, 0);
-        getPerpendicularAxis(jetAxis, ueAxis1, +1);
-        getPerpendicularAxis(jetAxis, ueAxis2, -1);
+        getPerpendicularDirections(jetAxis, ueAxis1, ueAxis2);
+        if (ueAxis1.Mag() == 0 || ueAxis2.Mag() == 0) {
+          continue;
+        }
 
         // Loop over MC particles to analyze underlying event region
         for (const auto& protonVec : protonMomentum) {
@@ -1757,18 +1874,21 @@ struct AntinucleiInJets {
             continue;
 
           // Fill normalization histogram
-          if (deltaRUe1 < maxConeRadius) {
-            registryMC.fill(HIST("antiproton_deltay_deltaphi_ue"), protonVec.Eta() - ueAxis1.Eta(), getDeltaPhi(protonVec.Phi(), ueAxis1.Phi()));
-          }
-          if (deltaRUe2 < maxConeRadius) {
-            registryMC.fill(HIST("antiproton_deltay_deltaphi_ue"), protonVec.Eta() - ueAxis2.Eta(), getDeltaPhi(protonVec.Phi(), ueAxis2.Phi()));
+          registryMC.fill(HIST("antiproton_deltay_deltaphi_ue"), protonVec.Eta() - ueAxis1.Eta(), getDeltaPhi(protonVec.Phi(), ueAxis1.Phi()));
+          registryMC.fill(HIST("antiproton_deltay_deltaphi_ue"), protonVec.Eta() - ueAxis2.Eta(), getDeltaPhi(protonVec.Phi(), ueAxis2.Phi()));
+
+          // Calculate weight
+          double weightUe(1.0);
+          if (applyReweighting && protonVec.Pt() < antiprotonsPerpCone->GetXaxis()->GetXmax()) {
+            int ipt = antiprotonsPerpCone->FindBin(protonVec.Pt());
+            weightUe = antiprotonsPerpCone->GetBinContent(ipt);
           }
 
           // Fill histogram for antiprotons in the UE
-          registryMC.fill(HIST("antiproton_gen_ue"), protonVec.Pt());
+          registryMC.fill(HIST("antiproton_gen_ue"), protonVec.Pt(), weightUe);
 
           // Fill 2d (pt,eta) distribution of antiprotons
-          registryMC.fill(HIST("antiproton_eta_pt_ue"), protonVec.Pt(), protonVec.Eta());
+          registryMC.fill(HIST("antiproton_eta_pt_ue"), protonVec.Pt(), protonVec.Eta(), weightUe);
         }
       }
       if (isAtLeastOneJetSelected) {
@@ -1910,8 +2030,10 @@ struct AntinucleiInJets {
         double coneRadius = std::sqrt(jet.area() / PI);
         TVector3 jetAxis(jet.px(), jet.py(), jet.pz());
         TVector3 ueAxis1(0, 0, 0), ueAxis2(0, 0, 0);
-        getPerpendicularAxis(jetAxis, ueAxis1, +1);
-        getPerpendicularAxis(jetAxis, ueAxis2, -1);
+        getPerpendicularDirections(jetAxis, ueAxis1, ueAxis2);
+        if (ueAxis1.Mag() == 0 || ueAxis2.Mag() == 0) {
+          continue;
+        }
 
         // Get jet constituents
         std::vector<fastjet::PseudoJet> jetConstituents = jet.constituents();
@@ -1981,11 +2103,18 @@ struct AntinucleiInJets {
           // Fill antiproton spectrum for physical primaries
           registryMC.fill(HIST("antiproton_prim_jet"), pt);
 
+          // Calculate weight
+          double weightJet(1.0);
+          if (applyReweighting && mcparticle.pt() < antiprotonsInsideJets->GetXaxis()->GetXmax()) {
+            int ipt = antiprotonsInsideJets->FindBin(mcparticle.pt());
+            weightJet = antiprotonsInsideJets->GetBinContent(ipt);
+          }
+
           // Fill histograms (TPC and TOF) only for selected candidates
           if (passedItsPidProt && nsigmaTPCPr > minNsigmaTpc && nsigmaTPCPr < maxNsigmaTpc) {
-            registryMC.fill(HIST("antiproton_rec_tpc_jet"), pt);
+            registryMC.fill(HIST("antiproton_rec_tpc_jet"), pt, weightJet);
             if (track.hasTOF() && nsigmaTOFPr > minNsigmaTof && nsigmaTOFPr < maxNsigmaTof) {
-              registryMC.fill(HIST("antiproton_rec_tof_jet"), pt);
+              registryMC.fill(HIST("antiproton_rec_tof_jet"), pt, weightJet);
             }
           }
         }
@@ -2056,11 +2185,18 @@ struct AntinucleiInJets {
           // Fill antiproton spectrum for physical primaries
           registryMC.fill(HIST("antiproton_prim_ue"), pt);
 
+          // Calculate weight
+          double weightUe(1.0);
+          if (applyReweighting && mcparticle.pt() < antiprotonsPerpCone->GetXaxis()->GetXmax()) {
+            int ipt = antiprotonsPerpCone->FindBin(mcparticle.pt());
+            weightUe = antiprotonsPerpCone->GetBinContent(ipt);
+          }
+
           // Fill histograms (TPC and TOF) only for selected candidates
           if (passedItsPidProt && nsigmaTPCPr > minNsigmaTpc && nsigmaTPCPr < maxNsigmaTpc) {
-            registryMC.fill(HIST("antiproton_rec_tpc_ue"), pt);
+            registryMC.fill(HIST("antiproton_rec_tpc_ue"), pt, weightUe);
             if (track.hasTOF() && nsigmaTOFPr > minNsigmaTof && nsigmaTOFPr < maxNsigmaTof) {
-              registryMC.fill(HIST("antiproton_rec_tof_ue"), pt);
+              registryMC.fill(HIST("antiproton_rec_tof_ue"), pt, weightUe);
             }
           }
         }
@@ -2624,10 +2760,11 @@ struct AntinucleiInJets {
       // Perpendicular cones
       double coneRadius = std::sqrt(jet.area() / PI);
       TVector3 jetAxis(jet.px(), jet.py(), jet.pz());
-      TVector3 ueAxis1(0, 0, 0);
-      TVector3 ueAxis2(0, 0, 0);
-      getPerpendicularAxis(jetAxis, ueAxis1, +1);
-      getPerpendicularAxis(jetAxis, ueAxis2, -1);
+      TVector3 ueAxis1(0, 0, 0), ueAxis2(0, 0, 0);
+      getPerpendicularDirections(jetAxis, ueAxis1, ueAxis2);
+      if (ueAxis1.Mag() == 0 || ueAxis2.Mag() == 0) {
+        continue;
+      }
 
       // Get jet constituents
       std::vector<fastjet::PseudoJet> jetConstituents = jet.constituents();