Skip to content

Commit 4133d27

Browse files
lorber98lorber98
authored
[PWGLF] Random cone to compute delta pT distributions (#16471)
1 parent d98d20e commit 4133d27

1 file changed

Lines changed: 87 additions & 40 deletions

File tree

‎PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx‎

Lines changed: 87 additions & 40 deletions
Original file line numberDiff line numberDiff line change
@@ -153,6 +153,7 @@ struct StrangenessInJetsIons {
153153
Configurable<bool> calculateFeeddownMatrix{"calculateFeeddownMatrix", true, "Fill feeddown matrix for Lambda if MC"};
154154
Configurable<bool> useV0inJetRec{"useV0inJetRec", true, "Include V0s in jet reconstruction"};
155155
Configurable<bool> doThermalToyModel{"doThermalToyModel", false, "Use the thermal toy model to embed background particles to the jet finder input list"};
156+
Configurable<bool> doPlotRho{"doPlotRho", false, "Plot rho distributions computed with perpendicular cone and area median method."};
156157
Configurable<bool> saveChargedParticleMB{"saveChargedParticleMB", false, "Store charged particle information to build inclusive spectra."};
157158

158159
// Event selection
@@ -281,6 +282,7 @@ struct StrangenessInJetsIons {
281282
const AxisSpec qtarmAxis{1000, 0.0f, 0.30f, "q_{T}^{arm}"};
282283
const AxisSpec numBkgParticles{500, 0.0, 500.0, "Number of particles"};
283284
const AxisSpec deltaPtAxis{2400, -60.0, 60.0, "#delta #it{p}_{T} (GeV/#it{c})"};
285+
const AxisSpec rhoAxis{100, 0.0, 50.0, "#it{#rho} (GeV/#it{c})"};
284286

285287
// Join enum ParticleOfInterest and the configurable vector particlesOfInterest in a map particleOfInterestDict
286288
const std::vector<int>& particleOnOff = particleOfInterest;
@@ -331,8 +333,11 @@ struct StrangenessInJetsIons {
331333
registryData.add("n_jets_vs_mult", "n_jets_vs_mult", HistType::kTH2F, {multAxis, numJets});
332334

333335
// Delta pT distribution
334-
registryData.add("delta_pT_data", "delta_pT_data", HistType::kTH2F, {multAxis, deltaPtAxis});
335-
registryData.add("delta_pT_RC_data", "delta_pT_RC_data", HistType::kTH2F, {multAxis, deltaPtAxis});
336+
registryData.add("h2_centrality_deltaPt_RandomCone", "h2_centrality_deltaPt_RandomCone", HistType::kTH2F, {multAxis, deltaPtAxis});
337+
registryData.add("h2_centrality_rhoPerp", "h2_centrality_rhoPerp", HistType::kTH2F, {multAxis, rhoAxis});
338+
339+
registryData.add("rho_perp", "rho_perp", HistType::kTH2F, {multAxis, rhoAxis});
340+
registryData.add("rho_median", "rho_median", HistType::kTH2F, {multAxis, rhoAxis});
336341

337342
// Armenteros-Podolanski plot
338343
// registryQC.add("ArmenterosPreSel_DATA", "ArmenterosPreSel_DATA", HistType::kTH2F, {alphaArmAxis, qtarmAxis});
@@ -396,10 +401,6 @@ struct StrangenessInJetsIons {
396401
registryMC.add("n_jets_vs_mult_pT_mc_gen", "n_jets_vs_mult_pT_mc_gen", HistType::kTH2F, {multAxis, ptJetAxis});
397402
registryMC.add("n_jets_vs_mult_mc_gen", "n_jets_vs_mult_mc_gen", HistType::kTH2F, {multAxis, numJets});
398403

399-
// Delta pT distribution
400-
registryMC.add("delta_pT_gen", "delta_pT_gen", HistType::kTH2F, {multAxis, deltaPtAxis});
401-
// registryMC.add("delta_pT_RC_gen", "delta_pT_RC_gen", HistType::kTH2F, {multAxis, deltaPtAxis});
402-
403404
if (doThermalToyModel) {
404405
registryMC.add("thermalToy_pT_gen", "thermalToy_pT_gen", HistType::kTH2F, {multAxis, ptAxis});
405406
registryMC.add("thermalToy_nBkg_gen", "thermalToy_nBkg_gen", HistType::kTH2F, {multAxis, numBkgParticles});
@@ -501,9 +502,6 @@ struct StrangenessInJetsIons {
501502
registryMC.add("n_jets_vs_mult_pT_mc_rec", "n_jets_vs_mult_pT_mc_rec", HistType::kTH2F, {multAxis, ptJetAxis});
502503
registryMC.add("n_jets_vs_mult_mc_rec", "n_jets_vs_mult_mc_rec", HistType::kTH2F, {multAxis, numJets});
503504

504-
// Delta pT distribution
505-
registryMC.add("delta_pT_rec", "delta_pT_rec", HistType::kTH2F, {multAxis, deltaPtAxis});
506-
507505
if (doThermalToyModel) {
508506
registryMC.add("thermalToy_pT_rec", "thermalToy_pT_rec", HistType::kTH2F, {multAxis, ptAxis});
509507
registryMC.add("thermalToy_nBkg_rec", "thermalToy_nBkg_rec", HistType::kTH2F, {multAxis, numBkgParticles});
@@ -757,6 +755,54 @@ struct StrangenessInJetsIons {
757755
u2.SetXYZ(u2x, u2y, pz);
758756
}
759757

758+
void computeRandomConeDeltaPt(const std::vector<fastjet::PseudoJet>& fjParticles,
759+
const std::vector<fastjet::PseudoJet>& jets,
760+
float multiplicity, double rhoPerp)
761+
{
762+
// Generate eta and phi for random cone in acceptance region
763+
double randomConeEta = fRng.Uniform(configTracks.etaMin + rJet, configTracks.etaMax - rJet);
764+
double randomConePhi = fRng.Uniform(0.0, TMath::TwoPi());
765+
766+
// Exclude leading jet region
767+
if (!jets.empty()) {
768+
float dPhiLeadingJet = getDeltaPhi(jets[0].phi(), randomConePhi);
769+
float dEtaLeadingJet = jets[0].eta() - randomConeEta;
770+
float dRLeadingJet = std::sqrt(dEtaLeadingJet * dEtaLeadingJet + dPhiLeadingJet * dPhiLeadingJet);
771+
772+
int attempts = 0;
773+
const int maxAttempts = 100;
774+
// If the random cone overlaps with the leading jet (distance < 1.5), then generate the coordinates again
775+
while (dRLeadingJet < 1.5 && attempts < maxAttempts) {
776+
randomConeEta = fRng.Uniform(configTracks.etaMin + rJet, configTracks.etaMax - rJet);
777+
randomConePhi = fRng.Uniform(0.0, TMath::TwoPi());
778+
779+
dPhiLeadingJet = getDeltaPhi(jets[0].phi(), randomConePhi);
780+
dEtaLeadingJet = jets[0].eta() - randomConeEta;
781+
dRLeadingJet = std::sqrt(dEtaLeadingJet * dEtaLeadingJet + dPhiLeadingJet * dPhiLeadingJet);
782+
attempts++;
783+
}
784+
}
785+
786+
// Sum pT of all the particles (charged tracks + V0 if included) falling in the random cone
787+
float randomConePt = 0.0;
788+
for (auto const& part : fjParticles) {
789+
float dPhi = getDeltaPhi(part.phi(), randomConePhi);
790+
float dEta = part.eta() - randomConeEta;
791+
float dR = std::sqrt(dEta * dEta + dPhi * dPhi);
792+
793+
if (dR < rJet) {
794+
randomConePt += part.pt();
795+
}
796+
}
797+
798+
// Compute delta pT: pT_cone - (Area_cone * rho)
799+
double coneArea = TMath::Pi() * rJet * rJet;
800+
double deltaPtRandomCone = randomConePt - (coneArea * rhoPerp);
801+
802+
registryData.fill(HIST("h2_centrality_deltaPt_RandomCone"), multiplicity, deltaPtRandomCone);
803+
registryData.fill(HIST("h2_centrality_rhoPerp"), multiplicity, rhoPerp);
804+
}
805+
760806
// Find ITS hit
761807
template <typename TrackIts>
762808
bool hasITSHitOnLayer(const TrackIts& track, int layer)
@@ -1940,19 +1986,20 @@ struct StrangenessInJetsIons {
19401986
pj.set_user_index(p.pdgCode());
19411987
v0PseudoJets.push_back(pj);
19421988
// LOG(info) << "[AddV0sForJetReconstructionMCP] Add V0 as input for jet finder.";
1989+
}
1990+
}
19431991

1944-
// Remove V0 daughter particles if already in the input list for the jet finder
1945-
for (long unsigned int i = 0; i < fjParticleObj.size(); ++i) {
1946-
const auto& mcPart = fjParticleObj[i];
1947-
if (!mcPart.has_mothers())
1948-
continue;
1949-
auto mother = mcParticles.iteratorAt(mcPart.mothersIds()[0]);
1950-
int motherPdg = std::abs(mother.pdgCode());
1951-
if (motherPdg == kK0Short || motherPdg == kLambda0) {
1952-
isTrackReplaced[i] = true;
1953-
// LOG(info) << "[AddV0sForJetReconstructionMCP] V0 daughter particle found in fjParticleObj.";
1954-
}
1955-
}
1992+
// Remove V0 daughter particles if already in the input list for the jet finder
1993+
for (long unsigned int i = 0; i < fjParticleObj.size(); ++i) {
1994+
const auto& mcPart = fjParticleObj[i];
1995+
if (!mcPart.has_mothers())
1996+
continue;
1997+
auto mother = mcParticles.iteratorAt(mcPart.mothersIds()[0]);
1998+
int motherPdg = std::abs(mother.pdgCode());
1999+
// LOG(info) << "[AddV0sForJetReconstructionMCP] Mother pdg code:" << motherPdg;
2000+
if (motherPdg == kK0Short || motherPdg == kLambda0) {
2001+
isTrackReplaced[i] = true;
2002+
// LOG(info) << "[AddV0sForJetReconstructionMCP] V0 daughter particle found in fjParticleObj.";
19562003
}
19572004
}
19582005

@@ -2036,6 +2083,9 @@ struct StrangenessInJetsIons {
20362083
double phi = fRng.Uniform(0.0, TMath::TwoPi());
20372084
double eta = fRng.Uniform(configTracks.etaMin, configTracks.etaMax);
20382085

2086+
if (pt < 0.1f)
2087+
continue;
2088+
20392089
double px = pt * std::cos(phi);
20402090
double py = pt * std::sin(phi);
20412091
double pz = pt * std::sinh(eta);
@@ -2150,13 +2200,6 @@ struct StrangenessInJetsIons {
21502200
std::vector<fastjet::PseudoJet> jets = fastjet::sorted_by_pt(cs.inclusive_jets());
21512201
auto [rhoPerp, rhoMPerp] = jetutilities::estimateRhoPerpCone(fjParticles, jets[0], rJet);
21522202

2153-
// Jet selection
2154-
bool isAtLeastOneJetSelected = false;
2155-
std::vector<TVector3> selectedJet;
2156-
std::vector<TVector3> ue1;
2157-
std::vector<TVector3> ue2;
2158-
std::vector<double> jetPt;
2159-
21602203
// Event multiplicity
21612204
float multiplicity;
21622205
if (centrEstimator == 0) {
@@ -2165,6 +2208,22 @@ struct StrangenessInJetsIons {
21652208
multiplicity = collision.centFT0M();
21662209
}
21672210

2211+
if (doPlotRho) {
2212+
auto [rhoMedian, rhoMMedian] = backgroundSub.estimateRhoAreaMedian(fjParticles, true);
2213+
registryData.fill(HIST("rho_perp"), multiplicity, rhoPerp);
2214+
registryData.fill(HIST("rho_median"), multiplicity, rhoMedian);
2215+
}
2216+
2217+
// Delta pT distributions with random cone technique
2218+
computeRandomConeDeltaPt(fjParticles, jets, multiplicity, rhoPerp);
2219+
2220+
// Jet selection
2221+
bool isAtLeastOneJetSelected = false;
2222+
std::vector<TVector3> selectedJet;
2223+
std::vector<TVector3> ue1;
2224+
std::vector<TVector3> ue2;
2225+
std::vector<double> jetPt;
2226+
21682227
// Loop over reconstructed jets
21692228
for (const auto& jet : jets) {
21702229

@@ -2179,10 +2238,6 @@ struct StrangenessInJetsIons {
21792238
continue;
21802239
isAtLeastOneJetSelected = true;
21812240

2182-
// delta pT distribution after jet selection by pT
2183-
double deltaPt = jet.pt() - jetMinusBkg.pt();
2184-
registryData.fill(HIST("delta_pT_data"), multiplicity, deltaPt);
2185-
21862241
// Calculation of perpendicular cones
21872242
TVector3 jetAxis(jet.px(), jet.py(), jet.pz());
21882243
TVector3 ueAxis1(0, 0, 0), ueAxis2(0, 0, 0);
@@ -2536,10 +2591,6 @@ struct StrangenessInJetsIons {
25362591
// Fill jet counter
25372592
registryMC.fill(HIST("n_jets_vs_mult_pT_mc_gen"), genMultiplicity, jetMinusBkg.pt());
25382593

2539-
// delta pT distribution after jet selection by pT
2540-
double deltaPt = jet.pt() - jetMinusBkg.pt();
2541-
registryMC.fill(HIST("delta_pT_gen"), genMultiplicity, deltaPt);
2542-
25432594
// Set up two perpendicular cone axes for underlying event estimation
25442595
TVector3 jetAxis(jet.px(), jet.py(), jet.pz());
25452596
double coneRadius = std::sqrt(jet.area() / PI); // TODO: replace with rJet (similar results)
@@ -2848,10 +2899,6 @@ struct StrangenessInJetsIons {
28482899
continue;
28492900
isAtLeastOneJetSelected = true;
28502901

2851-
// delta pT distribution after jet selection by pT
2852-
double deltaPt = jet.pt() - jetMinusBkg.pt();
2853-
registryMC.fill(HIST("delta_pT_rec"), multiplicity, deltaPt);
2854-
28552902
// Perpendicular cones
28562903
TVector3 jetAxis(jet.px(), jet.py(), jet.pz());
28572904
TVector3 ueAxis1(0, 0, 0), ueAxis2(0, 0, 0);

0 commit comments

Comments
 (0)