diff --git a/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx b/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx index e39ceb4020b..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 @@ -280,6 +281,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 +332,11 @@ 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}); // Armenteros-Podolanski plot // registryQC.add("ArmenterosPreSel_DATA", "ArmenterosPreSel_DATA", HistType::kTH2F, {alphaArmAxis, qtarmAxis}); @@ -395,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}); @@ -500,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}); @@ -756,6 +754,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) @@ -1939,19 +1985,20 @@ 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()); + // LOG(info) << "[AddV0sForJetReconstructionMCP] Mother pdg code:" << motherPdg; + if (motherPdg == kK0Short || motherPdg == kLambda0) { + isTrackReplaced[i] = true; + // LOG(info) << "[AddV0sForJetReconstructionMCP] V0 daughter particle found in fjParticleObj."; } } @@ -2035,6 +2082,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); @@ -2149,13 +2199,6 @@ struct StrangenessInJetsIons { 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; - // Event multiplicity float multiplicity; if (centrEstimator == 0) { @@ -2164,6 +2207,22 @@ struct StrangenessInJetsIons { 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 selectedJet; + std::vector ue1; + std::vector ue2; + std::vector jetPt; + // Loop over reconstructed jets for (const auto& jet : jets) { @@ -2178,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); @@ -2535,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) @@ -2847,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);