Skip to content

[PWGLF] Random cone to compute delta pT distributions #16471

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
romainschotter merged 3 commits into from
May 30, 2026
+87 −40
Merged

[PWGLF] Random cone to compute delta pT distributions #16471

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
ea2aebe
Add pT cut in thermal toy model
lorber98 May 30, 2026
c6c2506
Add delta pT distribution in data with random cone + QC plots for rho…
lorber98 May 30, 2026
034deaf
Configurable to plot rho distributions computed with perpendicular co…
lorber98 May 30, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
127 changes: 87 additions & 40 deletions PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -152,6 +152,7 @@
Configurable<bool> calculateFeeddownMatrix{"calculateFeeddownMatrix", true, "Fill feeddown matrix for Lambda if MC"};
Configurable<bool> useV0inJetRec{"useV0inJetRec", true, "Include V0s in jet reconstruction"};
Configurable<bool> doThermalToyModel{"doThermalToyModel", false, "Use the thermal toy model to embed background particles to the jet finder input list"};
Configurable<bool> doPlotRho{"doPlotRho", false, "Plot rho distributions computed with perpendicular cone and area median method."};
Configurable<bool> saveChargedParticleMB{"saveChargedParticleMB", false, "Store charged particle information to build inclusive spectra."};

// Event selection
Expand Down Expand Up @@ -280,6 +281,7 @@
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<int>& particleOnOff = particleOfInterest;
Expand Down Expand Up @@ -330,8 +332,11 @@
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});

// Armenteros-Podolanski plot
// registryQC.add("ArmenterosPreSel_DATA", "ArmenterosPreSel_DATA", HistType::kTH2F, {alphaArmAxis, qtarmAxis});
Expand Down Expand Up @@ -395,10 +400,6 @@
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});
Expand Down Expand Up @@ -500,9 +501,6 @@
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});
Expand Down Expand Up @@ -756,6 +754,54 @@
u2.SetXYZ(u2x, u2y, pz);
}

void computeRandomConeDeltaPt(const std::vector<fastjet::PseudoJet>& fjParticles,
const std::vector<fastjet::PseudoJet>& 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);

Check failure on line 775 in PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number] 

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
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 <typename TrackIts>
bool hasITSHitOnLayer(const TrackIts& track, int layer)
Expand Down Expand Up @@ -1471,7 +1517,7 @@
continue;

if (std::abs(motherPos.eta()) > 0.8)
continue;

Check failure on line 1520 in PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number] 

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.

// Vertex position vector
TVector3 vtxPos(collision.posX(), collision.posY(), collision.posZ());
Expand Down Expand Up @@ -1570,7 +1616,7 @@
continue;

if (std::abs(motherPos.eta()) > 0.8)
continue;

Check failure on line 1619 in PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number] 

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.

// Xi+
if (passedXiSelection(casc, pos, neg, bach, collision) && bach.sign() > 0 && motherBach.pdgCode() == kXiPlusBar) {
Expand Down Expand Up @@ -1608,7 +1654,7 @@
}

if (std::abs(mcParticle.eta()) > 0.8) {
continue;

Check failure on line 1657 in PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number] 

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
}

switch (mcParticle.pdgCode()) {
Expand Down Expand Up @@ -1665,7 +1711,7 @@
rapidityXi = RecoDecay::y(std::array{v0mother.px(), v0mother.py(), v0mother.pz()}, o2::constants::physics::MassXi0);

if (std::fabs(rapidityXi) > 0.5f)
return; // not a valid mother rapidity (PDG selection is later)

Check failure on line 1714 in PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number] 

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.

// __________________________________________
if (passedLambda) {
Expand Down Expand Up @@ -1939,19 +1985,20 @@
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());
// LOG(info) << "[AddV0sForJetReconstructionMCP] Mother pdg code:" << motherPdg;
if (motherPdg == kK0Short || motherPdg == kLambda0) {
isTrackReplaced[i] = true;
// LOG(info) << "[AddV0sForJetReconstructionMCP] V0 daughter particle found in fjParticleObj.";
}
}

Expand Down Expand Up @@ -2035,6 +2082,9 @@
double phi = fRng.Uniform(0.0, TMath::TwoPi());
double eta = fRng.Uniform(configTracks.etaMin, configTracks.etaMax);

if (pt < 0.1f)
continue;

Check failure on line 2086 in PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number] 

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.

double px = pt * std::cos(phi);
double py = pt * std::sin(phi);
double pz = pt * std::sinh(eta);
Expand Down Expand Up @@ -2149,13 +2199,6 @@
std::vector<fastjet::PseudoJet> jets = fastjet::sorted_by_pt(cs.inclusive_jets());
auto [rhoPerp, rhoMPerp] = jetutilities::estimateRhoPerpCone(fjParticles, jets[0], rJet);

// Jet selection
bool isAtLeastOneJetSelected = false;
std::vector<TVector3> selectedJet;
std::vector<TVector3> ue1;
std::vector<TVector3> ue2;
std::vector<double> jetPt;

// Event multiplicity
float multiplicity;
if (centrEstimator == 0) {
Expand All @@ -2164,6 +2207,22 @@
multiplicity = collision.centFT0M();
}

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);

// Jet selection
bool isAtLeastOneJetSelected = false;
std::vector<TVector3> selectedJet;
std::vector<TVector3> ue1;
std::vector<TVector3> ue2;
std::vector<double> jetPt;

// Loop over reconstructed jets
for (const auto& jet : jets) {

Expand All @@ -2178,10 +2237,6 @@
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);
Expand Down Expand Up @@ -2428,7 +2483,7 @@
}

if (multiplicity == -999)
continue; */

Check failure on line 2486 in PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number] 

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.

// Clear containers at the start of the event loop
fjParticles.clear();
Expand Down Expand Up @@ -2535,10 +2590,6 @@
// 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)
Expand Down Expand Up @@ -2847,10 +2898,6 @@
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);
Expand Down Expand Up @@ -2883,7 +2930,7 @@
// --- Generated hadrons in reconstructed jets ----
for (const auto& particle : mcParticlesPerColl) {
if (!particle.isPhysicalPrimary() || std::abs(particle.eta()) > 0.8)
continue;

Check failure on line 2933 in PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number] 

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.

int absPdg = std::abs(particle.pdgCode());
if (absPdg != kK0Short && absPdg != kLambda0)
Expand Down
Loading