From ea2aebeb3e1e9f114f3cbd05d036a2c00cac8ac8 Mon Sep 17 00:00:00 2001 From: lorber98 <95907752+lorber98@users.noreply.github.com> Date: Sat, 30 May 2026 17:00:24 +0200 Subject: [PATCH 1/3] Add pT cut in thermal toy model --- .../Strangeness/strangenessInJetsIons.cxx | 27 ++++++++++--------- 1 file changed, 15 insertions(+), 12 deletions(-) diff --git a/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx b/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx index e39ceb4020b..ce62123cefd 100644 --- a/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx +++ b/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx @@ -1939,19 +1939,19 @@ struct StrangenessInJetsIons { pj.set_user_index(p.pdgCode()); v0PseudoJets.push_back(pj); // LOG(info) << "[AddV0sForJetReconstructionMCP] Add V0 as input for jet finder."; + } + } - // Remove V0 daughter particles if already in the input list for the jet finder - for (long unsigned int i = 0; i < fjParticleObj.size(); ++i) { - const auto& mcPart = fjParticleObj[i]; - if (!mcPart.has_mothers()) - continue; - auto mother = mcParticles.iteratorAt(mcPart.mothersIds()[0]); - int motherPdg = std::abs(mother.pdgCode()); - if (motherPdg == kK0Short || motherPdg == kLambda0) { - isTrackReplaced[i] = true; - // LOG(info) << "[AddV0sForJetReconstructionMCP] V0 daughter particle found in fjParticleObj."; - } - } + // Remove V0 daughter particles if already in the input list for the jet finder + for (long unsigned int i = 0; i < fjParticleObj.size(); ++i) { + const auto& mcPart = fjParticleObj[i]; + if (!mcPart.has_mothers()) + continue; + auto mother = mcParticles.iteratorAt(mcPart.mothersIds()[0]); + int motherPdg = std::abs(mother.pdgCode()); + if (motherPdg == kK0Short || motherPdg == kLambda0) { + isTrackReplaced[i] = true; + // LOG(info) << "[AddV0sForJetReconstructionMCP] V0 daughter particle found in fjParticleObj."; } } @@ -2035,6 +2035,9 @@ struct StrangenessInJetsIons { double phi = fRng.Uniform(0.0, TMath::TwoPi()); double eta = fRng.Uniform(configTracks.etaMin, configTracks.etaMax); + if (pt < 0.1f) + continue; + double px = pt * std::cos(phi); double py = pt * std::sin(phi); double pz = pt * std::sinh(eta); From c6c2506ada31099c1c44c3b71f9fd5b76f711851 Mon Sep 17 00:00:00 2001 From: lorber98 <95907752+lorber98@users.noreply.github.com> Date: Sat, 30 May 2026 22:30:43 +0200 Subject: [PATCH 2/3] Add delta pT distribution in data with random cone + QC plots for rho estimated with area median and perpendicular cone methods --- .../Strangeness/strangenessInJetsIons.cxx | 84 ++++++++++++++++--- 1 file changed, 73 insertions(+), 11 deletions(-) diff --git a/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx b/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx index ce62123cefd..b4d3cfe55fa 100644 --- a/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx +++ b/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx @@ -280,6 +280,7 @@ struct StrangenessInJetsIons { const AxisSpec qtarmAxis{1000, 0.0f, 0.30f, "q_{T}^{arm}"}; const AxisSpec numBkgParticles{500, 0.0, 500.0, "Number of particles"}; const AxisSpec deltaPtAxis{2400, -60.0, 60.0, "#delta #it{p}_{T} (GeV/#it{c})"}; + const AxisSpec rhoAxis{100, 0.0, 50.0, "#it{#rho} (GeV/#it{c})"}; // Join enum ParticleOfInterest and the configurable vector particlesOfInterest in a map particleOfInterestDict const std::vector& particleOnOff = particleOfInterest; @@ -330,8 +331,13 @@ struct StrangenessInJetsIons { registryData.add("n_jets_vs_mult", "n_jets_vs_mult", HistType::kTH2F, {multAxis, numJets}); // Delta pT distribution - registryData.add("delta_pT_data", "delta_pT_data", HistType::kTH2F, {multAxis, deltaPtAxis}); - registryData.add("delta_pT_RC_data", "delta_pT_RC_data", HistType::kTH2F, {multAxis, deltaPtAxis}); + // registryData.add("delta_pT_data", "delta_pT_data", HistType::kTH2F, {multAxis, deltaPtAxis}); + // registryData.add("delta_pT_RC_data", "delta_pT_RC_data", HistType::kTH2F, {multAxis, deltaPtAxis}); + registryData.add("h2_centrality_deltaPt_RandomCone", "h2_centrality_deltaPt_RandomCone", HistType::kTH2F, {multAxis, deltaPtAxis}); + registryData.add("h2_centrality_rhoPerp", "h2_centrality_rhoPerp", HistType::kTH2F, {multAxis, rhoAxis}); + + registryData.add("rho_perp", "rho_perp", HistType::kTH2F, {multAxis, rhoAxis}); + registryData.add("rho_median", "rho_median", HistType::kTH2F, {multAxis, rhoAxis}); // Armenteros-Podolanski plot // registryQC.add("ArmenterosPreSel_DATA", "ArmenterosPreSel_DATA", HistType::kTH2F, {alphaArmAxis, qtarmAxis}); @@ -756,6 +762,54 @@ struct StrangenessInJetsIons { u2.SetXYZ(u2x, u2y, pz); } + void computeRandomConeDeltaPt(const std::vector& fjParticles, + const std::vector& jets, + float multiplicity, double rhoPerp) + { + // Generate eta and phi for random cone in acceptance region + double randomConeEta = fRng.Uniform(configTracks.etaMin + rJet, configTracks.etaMax - rJet); + double randomConePhi = fRng.Uniform(0.0, TMath::TwoPi()); + + // Exclude leading jet region + if (!jets.empty()) { + float dPhiLeadingJet = getDeltaPhi(jets[0].phi(), randomConePhi); + float dEtaLeadingJet = jets[0].eta() - randomConeEta; + float dRLeadingJet = std::sqrt(dEtaLeadingJet * dEtaLeadingJet + dPhiLeadingJet * dPhiLeadingJet); + + int attempts = 0; + const int maxAttempts = 100; + // If the random cone overlaps with the leading jet (distance < 1.5), then generate the coordinates again + while (dRLeadingJet < 1.5 && attempts < maxAttempts) { + randomConeEta = fRng.Uniform(configTracks.etaMin + rJet, configTracks.etaMax - rJet); + randomConePhi = fRng.Uniform(0.0, TMath::TwoPi()); + + dPhiLeadingJet = getDeltaPhi(jets[0].phi(), randomConePhi); + dEtaLeadingJet = jets[0].eta() - randomConeEta; + dRLeadingJet = std::sqrt(dEtaLeadingJet * dEtaLeadingJet + dPhiLeadingJet * dPhiLeadingJet); + attempts++; + } + } + + // Sum pT of all the particles (charged tracks + V0 if included) falling in the random cone + float randomConePt = 0.0; + for (auto const& part : fjParticles) { + float dPhi = getDeltaPhi(part.phi(), randomConePhi); + float dEta = part.eta() - randomConeEta; + float dR = std::sqrt(dEta * dEta + dPhi * dPhi); + + if (dR < rJet) { + randomConePt += part.pt(); + } + } + + // Compute delta pT: pT_cone - (Area_cone * rho) + double coneArea = TMath::Pi() * rJet * rJet; + double deltaPtRandomCone = randomConePt - (coneArea * rhoPerp); + + registryData.fill(HIST("h2_centrality_deltaPt_RandomCone"), multiplicity, deltaPtRandomCone); + registryData.fill(HIST("h2_centrality_rhoPerp"), multiplicity, rhoPerp); + } + // Find ITS hit template bool hasITSHitOnLayer(const TrackIts& track, int layer) @@ -1949,6 +2003,7 @@ struct StrangenessInJetsIons { continue; auto mother = mcParticles.iteratorAt(mcPart.mothersIds()[0]); int motherPdg = std::abs(mother.pdgCode()); + LOG(info) << "[AddV0sForJetReconstructionMCP] Mother pdg code:" << motherPdg; if (motherPdg == kK0Short || motherPdg == kLambda0) { isTrackReplaced[i] = true; // LOG(info) << "[AddV0sForJetReconstructionMCP] V0 daughter particle found in fjParticleObj."; @@ -2151,13 +2206,7 @@ struct StrangenessInJetsIons { fastjet::ClusterSequenceArea cs(fjParticles, jetDef, areaDef); std::vector jets = fastjet::sorted_by_pt(cs.inclusive_jets()); auto [rhoPerp, rhoMPerp] = jetutilities::estimateRhoPerpCone(fjParticles, jets[0], rJet); - - // Jet selection - bool isAtLeastOneJetSelected = false; - std::vector selectedJet; - std::vector ue1; - std::vector ue2; - std::vector jetPt; + auto [rhoMedian, rhoMMedian] = backgroundSub.estimateRhoAreaMedian(fjParticles, true); // Event multiplicity float multiplicity; @@ -2167,6 +2216,19 @@ struct StrangenessInJetsIons { multiplicity = collision.centFT0M(); } + registryData.fill(HIST("rho_perp"), multiplicity, rhoPerp); + registryData.fill(HIST("rho_median"), multiplicity, rhoMedian); + + // Delta pT distributions with random cone technique + computeRandomConeDeltaPt(fjParticles, jets, multiplicity, rhoPerp); + + // Jet selection + bool isAtLeastOneJetSelected = false; + std::vector selectedJet; + std::vector ue1; + std::vector ue2; + std::vector jetPt; + // Loop over reconstructed jets for (const auto& jet : jets) { @@ -2182,8 +2244,8 @@ struct StrangenessInJetsIons { isAtLeastOneJetSelected = true; // delta pT distribution after jet selection by pT - double deltaPt = jet.pt() - jetMinusBkg.pt(); - registryData.fill(HIST("delta_pT_data"), multiplicity, deltaPt); + // double deltaPt = jet.pt() - jetMinusBkg.pt(); + // registryData.fill(HIST("delta_pT_data"), multiplicity, deltaPt); // Calculation of perpendicular cones TVector3 jetAxis(jet.px(), jet.py(), jet.pz()); From 034deafa682f79cf7587c180e10fba301c2e417f Mon Sep 17 00:00:00 2001 From: lorber98 <95907752+lorber98@users.noreply.github.com> Date: Sat, 30 May 2026 22:41:01 +0200 Subject: [PATCH 3/3] Configurable to plot rho distributions computed with perpendicular cone and area median method --- .../Strangeness/strangenessInJetsIons.cxx | 34 +++++-------------- 1 file changed, 8 insertions(+), 26 deletions(-) diff --git a/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx b/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx index b4d3cfe55fa..4fb4f0cffdb 100644 --- a/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx +++ b/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx @@ -152,6 +152,7 @@ struct StrangenessInJetsIons { Configurable calculateFeeddownMatrix{"calculateFeeddownMatrix", true, "Fill feeddown matrix for Lambda if MC"}; Configurable useV0inJetRec{"useV0inJetRec", true, "Include V0s in jet reconstruction"}; Configurable doThermalToyModel{"doThermalToyModel", false, "Use the thermal toy model to embed background particles to the jet finder input list"}; + Configurable doPlotRho{"doPlotRho", false, "Plot rho distributions computed with perpendicular cone and area median method."}; Configurable saveChargedParticleMB{"saveChargedParticleMB", false, "Store charged particle information to build inclusive spectra."}; // Event selection @@ -331,11 +332,9 @@ struct StrangenessInJetsIons { registryData.add("n_jets_vs_mult", "n_jets_vs_mult", HistType::kTH2F, {multAxis, numJets}); // Delta pT distribution - // registryData.add("delta_pT_data", "delta_pT_data", HistType::kTH2F, {multAxis, deltaPtAxis}); - // registryData.add("delta_pT_RC_data", "delta_pT_RC_data", HistType::kTH2F, {multAxis, deltaPtAxis}); registryData.add("h2_centrality_deltaPt_RandomCone", "h2_centrality_deltaPt_RandomCone", HistType::kTH2F, {multAxis, deltaPtAxis}); registryData.add("h2_centrality_rhoPerp", "h2_centrality_rhoPerp", HistType::kTH2F, {multAxis, rhoAxis}); - + registryData.add("rho_perp", "rho_perp", HistType::kTH2F, {multAxis, rhoAxis}); registryData.add("rho_median", "rho_median", HistType::kTH2F, {multAxis, rhoAxis}); @@ -401,10 +400,6 @@ struct StrangenessInJetsIons { registryMC.add("n_jets_vs_mult_pT_mc_gen", "n_jets_vs_mult_pT_mc_gen", HistType::kTH2F, {multAxis, ptJetAxis}); registryMC.add("n_jets_vs_mult_mc_gen", "n_jets_vs_mult_mc_gen", HistType::kTH2F, {multAxis, numJets}); - // Delta pT distribution - registryMC.add("delta_pT_gen", "delta_pT_gen", HistType::kTH2F, {multAxis, deltaPtAxis}); - // registryMC.add("delta_pT_RC_gen", "delta_pT_RC_gen", HistType::kTH2F, {multAxis, deltaPtAxis}); - if (doThermalToyModel) { registryMC.add("thermalToy_pT_gen", "thermalToy_pT_gen", HistType::kTH2F, {multAxis, ptAxis}); registryMC.add("thermalToy_nBkg_gen", "thermalToy_nBkg_gen", HistType::kTH2F, {multAxis, numBkgParticles}); @@ -506,9 +501,6 @@ struct StrangenessInJetsIons { registryMC.add("n_jets_vs_mult_pT_mc_rec", "n_jets_vs_mult_pT_mc_rec", HistType::kTH2F, {multAxis, ptJetAxis}); registryMC.add("n_jets_vs_mult_mc_rec", "n_jets_vs_mult_mc_rec", HistType::kTH2F, {multAxis, numJets}); - // Delta pT distribution - registryMC.add("delta_pT_rec", "delta_pT_rec", HistType::kTH2F, {multAxis, deltaPtAxis}); - if (doThermalToyModel) { registryMC.add("thermalToy_pT_rec", "thermalToy_pT_rec", HistType::kTH2F, {multAxis, ptAxis}); registryMC.add("thermalToy_nBkg_rec", "thermalToy_nBkg_rec", HistType::kTH2F, {multAxis, numBkgParticles}); @@ -2003,7 +1995,7 @@ struct StrangenessInJetsIons { continue; auto mother = mcParticles.iteratorAt(mcPart.mothersIds()[0]); int motherPdg = std::abs(mother.pdgCode()); - LOG(info) << "[AddV0sForJetReconstructionMCP] Mother pdg code:" << motherPdg; + // LOG(info) << "[AddV0sForJetReconstructionMCP] Mother pdg code:" << motherPdg; if (motherPdg == kK0Short || motherPdg == kLambda0) { isTrackReplaced[i] = true; // LOG(info) << "[AddV0sForJetReconstructionMCP] V0 daughter particle found in fjParticleObj."; @@ -2206,7 +2198,6 @@ struct StrangenessInJetsIons { fastjet::ClusterSequenceArea cs(fjParticles, jetDef, areaDef); std::vector jets = fastjet::sorted_by_pt(cs.inclusive_jets()); auto [rhoPerp, rhoMPerp] = jetutilities::estimateRhoPerpCone(fjParticles, jets[0], rJet); - auto [rhoMedian, rhoMMedian] = backgroundSub.estimateRhoAreaMedian(fjParticles, true); // Event multiplicity float multiplicity; @@ -2216,8 +2207,11 @@ struct StrangenessInJetsIons { multiplicity = collision.centFT0M(); } - registryData.fill(HIST("rho_perp"), multiplicity, rhoPerp); - registryData.fill(HIST("rho_median"), multiplicity, rhoMedian); + if (doPlotRho) { + auto [rhoMedian, rhoMMedian] = backgroundSub.estimateRhoAreaMedian(fjParticles, true); + registryData.fill(HIST("rho_perp"), multiplicity, rhoPerp); + registryData.fill(HIST("rho_median"), multiplicity, rhoMedian); + } // Delta pT distributions with random cone technique computeRandomConeDeltaPt(fjParticles, jets, multiplicity, rhoPerp); @@ -2243,10 +2237,6 @@ struct StrangenessInJetsIons { continue; isAtLeastOneJetSelected = true; - // delta pT distribution after jet selection by pT - // double deltaPt = jet.pt() - jetMinusBkg.pt(); - // registryData.fill(HIST("delta_pT_data"), multiplicity, deltaPt); - // Calculation of perpendicular cones TVector3 jetAxis(jet.px(), jet.py(), jet.pz()); TVector3 ueAxis1(0, 0, 0), ueAxis2(0, 0, 0); @@ -2600,10 +2590,6 @@ struct StrangenessInJetsIons { // Fill jet counter registryMC.fill(HIST("n_jets_vs_mult_pT_mc_gen"), genMultiplicity, jetMinusBkg.pt()); - // delta pT distribution after jet selection by pT - double deltaPt = jet.pt() - jetMinusBkg.pt(); - registryMC.fill(HIST("delta_pT_gen"), genMultiplicity, deltaPt); - // 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); // TODO: replace with rJet (similar results) @@ -2912,10 +2898,6 @@ struct StrangenessInJetsIons { continue; isAtLeastOneJetSelected = true; - // delta pT distribution after jet selection by pT - double deltaPt = jet.pt() - jetMinusBkg.pt(); - registryMC.fill(HIST("delta_pT_rec"), multiplicity, deltaPt); - // Perpendicular cones TVector3 jetAxis(jet.px(), jet.py(), jet.pz()); TVector3 ueAxis1(0, 0, 0), ueAxis2(0, 0, 0);