From 1c73631ecf0d8f708b66e6ad37954cade4681370 Mon Sep 17 00:00:00 2001 From: Alberto Caliva Date: Sun, 19 Oct 2025 10:24:08 +0200 Subject: [PATCH 1/2] [PWGLF] Select primary and from HF for jet clustering,improved perp cone calculation --- PWGLF/Tasks/Nuspex/antinucleiInJets.cxx | 217 +++++++++++++++++++----- 1 file changed, 177 insertions(+), 40 deletions(-) diff --git a/PWGLF/Tasks/Nuspex/antinucleiInJets.cxx b/PWGLF/Tasks/Nuspex/antinucleiInJets.cxx index 1258c5cf330..d22d88b3de8 100644 --- a/PWGLF/Tasks/Nuspex/antinucleiInJets.cxx +++ b/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 ptLeadingMin{"ptLeadingMin", 5.0, "pt Leading Min"}; Configurable rJet{"rJet", 0.3, "Jet resolution parameter R"}; Configurable zVtx{"zVtx", 10.0, "Maximum zVertex"}; - Configurable applyAreaCut{"applyAreaCut", false, "apply area cut"}; - Configurable maxNormalizedJetArea{"maxNormalizedJetArea", 0.6, "area cut"}; + Configurable applyAreaCut{"applyAreaCut", true, "apply area cut"}; + Configurable maxNormalizedJetArea{"maxNormalizedJetArea", 1.0, "area cut"}; Configurable deltaEtaEdge{"deltaEtaEdge", 0.05, "eta gap from the edge"}; Configurable nSyst{"nSyst", 50, "number of systematic variations"}; // Track quality, kinematic, and PID selection parameters Configurable requirePvContributor{"requirePvContributor", false, "require that the track is a PV contributor"}; - Configurable applyItsPid{"applyItsPid", true, "apply ITS PID"}; + Configurable applyItsPid{"applyItsPid", false, "apply ITS PID"}; Configurable minItsNclusters{"minItsNclusters", 5, "minimum number of ITS clusters"}; Configurable minTpcNcrossedRows{"minTpcNcrossedRows", 100, "minimum number of TPC crossed pad rows"}; Configurable minChiSquareTpc{"minChiSquareTpc", 0.0, "minimum TPC chi^2/Ncls"}; @@ -138,18 +138,21 @@ struct AntinucleiInJets { Configurable ptMaxItsPidHel{"ptMaxItsPidHel", 1.0, "maximum pt for ITS PID for helium"}; Configurable nSigmaItsMin{"nSigmaItsMin", -3.0, "nSigmaITS min"}; Configurable nSigmaItsMax{"nSigmaItsMax", +3.0, "nSigmaITS max"}; - Configurable setMCDefaultItsParams{"setMCDefaultItsParams", false, "set MC default parameters"}; + Configurable setMCDefaultItsParams{"setMCDefaultItsParams", true, "set MC default parameters"}; Configurable cfgCompensatePIDinTracking{"cfgCompensatePIDinTracking", false, "If true, divide tpcInnerParam by the electric charge"}; Configurable> 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 applyReweighting{"applyReweighting", true, "enable reweighting for efficiency"}; Configurable urlToCcdb{"urlToCcdb", "http://alice-ccdb.cern.ch", "url of the personal ccdb"}; - Configurable pathToFile{"pathToFile", "Users/a/alcaliva/PrimaryFractionAntip", "path to file with reweighting"}; + Configurable pathToFile{"pathToFile", "Users/a/alcaliva/reweightingHistogramsAntipInJet", "path to file"}; Configurable weightsProton{"weightsProton", "", "weightsProton"}; Configurable weightsLambda{"weightsLambda", "", "weightsLambda"}; Configurable weightsSigma{"weightsSigma", "", "weightsSigma"}; Configurable weightsXi{"weightsXi", "", "weightsXi"}; Configurable weightsOmega{"weightsOmega", "", "weightsOmega"}; + Configurable weightsJet{"weightsJet", "", "weightsJet"}; + Configurable weightsUe{"weightsUe", "", "weightsUe"}; // Number of events Configurable 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 ccdb; @@ -205,7 +210,7 @@ struct AntinucleiInJets { ccdb->setLocalObjectValidityChecking(); ccdb->setCreatedNotAfter(std::chrono::duration_cast(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,7 +493,7 @@ struct AntinucleiInJets { } } - void getReweightingHistograms(o2::framework::Service const& ccdbObj, TString filepath, TString antip, TString antilambda, TString antisigma, TString antixi, TString antiomega) + void getReweightingHistograms(o2::framework::Service const& ccdbObj, TString filepath, TString antip, TString antilambda, TString antisigma, TString antixi, TString antiomega, TString jet, TString ue) { TList* list = ccdbObj->get(filepath.Data()); if (!list) { @@ -493,6 +501,7 @@ struct AntinucleiInJets { return; } + // Get reweighting histograms for primary fraction primaryAntiprotons = static_cast(list->FindObject(antip)); primaryAntiLambda = static_cast(list->FindObject(antilambda)); primaryAntiSigma = static_cast(list->FindObject(antisigma)); @@ -500,7 +509,14 @@ struct AntinucleiInJets { primaryAntiOmega = static_cast(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(list->FindObject(jet)); + antiprotonsPerpCone = static_cast(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 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(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 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 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 jetConstituents = jet.constituents(); From 7cb97c94e1ba75fcec56ccf1e836d1834cda8770 Mon Sep 17 00:00:00 2001 From: Alberto Caliva Date: Sun, 19 Oct 2025 10:39:52 +0200 Subject: [PATCH 2/2] fixed o2linter and clang format --- PWGLF/Tasks/Nuspex/antinucleiInJets.cxx | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/PWGLF/Tasks/Nuspex/antinucleiInJets.cxx b/PWGLF/Tasks/Nuspex/antinucleiInJets.cxx index d22d88b3de8..3103595c4de 100644 --- a/PWGLF/Tasks/Nuspex/antinucleiInJets.cxx +++ b/PWGLF/Tasks/Nuspex/antinucleiInJets.cxx @@ -554,15 +554,21 @@ struct AntinucleiInJets { { // Keep only pi, K, p, e, mu int pdg = std::abs(particle.pdgCode()); - if (!(pdg == 211 || pdg == 321 || pdg == 2212 || pdg == 11 || pdg == 13)) + if (!(pdg == PDG_t::kPiPlus || pdg == PDG_t::kKPlus || pdg == PDG_t::kProton || pdg == PDG_t::kElectron || pdg == PDG_t::kMuonMinus)) return false; + // Constants for identifying heavy-flavor (charm and bottom) content from PDG codes + static constexpr int kCharmQuark = 4; + static constexpr int kBottomQuark = 5; + static constexpr int hundreds = 100; + static constexpr int thousands = 1000; + // 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); + fromHF = (motherPdg / hundreds == kCharmQuark || motherPdg / hundreds == kBottomQuark || motherPdg / thousands == kCharmQuark || motherPdg / thousands == kBottomQuark); } // Select only physical primary particles or from heavy-flavor @@ -1760,7 +1766,7 @@ struct AntinucleiInJets { for (const auto& particle : mcParticles) { // Select physical primary particles or HF decay products - if (!isPhysicalPrimaryOrFromHF(particle,mcParticles)) + if (!isPhysicalPrimaryOrFromHF(particle, mcParticles)) continue; // Select particles within acceptance