Skip to content

Commit be9f1ea

Browse files
lorber98lorber98
committed
Prepare histograms for unfolding + delta pT distributions for MC
1 parent 978b850 commit be9f1ea

1 file changed

Lines changed: 227 additions & 3 deletions

File tree

‎PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx‎

Lines changed: 227 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -269,6 +269,8 @@ struct StrangenessInJetsIons {
269269

270270
const AxisSpec ptAxis{500, 0.0, 50.0, "#it{p}_{T} (GeV/#it{c})"};
271271
const AxisSpec ptJetAxis{100, 0.0, 100.0, "#it{p}_{T,jet} (GeV/#it{c})"};
272+
const AxisSpec ptJetRecoAxisUnfold{200, 0.0, 200.0, "#it{p}_{T,jet}^{rec} (GeV/#it{c})"};
273+
const AxisSpec ptJetGenAxisUnfold{200, 0.0, 200.0, "#it{p}_{T,jet}^{gen} (GeV/#it{c})"};
272274
const AxisSpec ptChargedAxis{10000, 0.0, 100.0, "#it{p}_{T} (GeV/#it{c})"};
273275
const AxisSpec numJets{21, -0.5, 20.5, "Number of jets per collision"};
274276
const AxisSpec invMassK0sAxis{200, 0.44, 0.56, "m_{#pi#pi} (GeV/#it{c}^{2})"};
@@ -413,6 +415,9 @@ struct StrangenessInJetsIons {
413415
registryMC.add("thermalToy_nBkg_gen", "thermalToy_nBkg_gen", HistType::kTH2F, {multAxis, numBkgParticles});
414416
}
415417

418+
registryMC.add("h2_centrality_deltaPt_RandomCone_gen", "h2_centrality_deltaPt_RandomCone_gen", HistType::kTH2F, {multAxis, deltaPtAxis});
419+
registryMC.add("h2_centrality_rhoPerp_gen", "h2_centrality_rhoPerp_gen", HistType::kTH2F, {multAxis, rhoAxis});
420+
416421
// Histograms for analysis
417422
if (particleOfInterestDict[ParticleOfInterest::kV0Particles]) {
418423
registryMC.add("K0s_generated_jet", "K0s_generated_jet", HistType::kTH2F, {multAxis, ptAxis});
@@ -522,6 +527,10 @@ struct StrangenessInJetsIons {
522527
registryMC.add("charged_embedJets_rec_jet", "charged_embedJets_rec_jet", HistType::kTH2F, {multAxis, ptAxis});
523528
registryMC.add("charged_embedJets_rec_ue", "charged_embedJets_rec_ue", HistType::kTH2F, {multAxis, ptAxis});
524529
}
530+
531+
registryMC.add("h2_centrality_deltaPt_RandomCone_rec", "h2_centrality_deltaPt_RandomCone_rec", HistType::kTH2F, {multAxis, deltaPtAxis});
532+
registryMC.add("h2_centrality_rhoPerp_rec", "h2_centrality_rhoPerp_rec", HistType::kTH2F, {multAxis, rhoAxis});
533+
525534
// Armenteros-Podolanski plot
526535
// registryQC.add("ArmenterosPreSel_REC", "ArmenterosPreSel_REC", HistType::kTH2F, {alphaArmAxis, qtarmAxis});
527536

@@ -669,6 +678,22 @@ struct StrangenessInJetsIons {
669678
registryDataMB.add("ChargedTrack_in_MB", "ChargedTrack_in_MB", HistType::kTH2F, {multAxis, ptChargedAxis});
670679
}
671680
}
681+
682+
if (doprocessPrepareUnfolding) {
683+
// Response matrix detector (only mathced)
684+
// Assi: X = Centrality, Y = pT_gen, Z = pT_reco
685+
registryMC.add("hDetectorResponse", "hDetectorResponse", HistType::kTH3F, {multAxis, ptJetGenAxisUnfold, ptJetRecoAxisUnfold});
686+
687+
// Efficiency components (GEN level)
688+
// Assi: X = Centrality, Y = pT_gen
689+
registryMC.add("hJetPtGenTotal", "hJetPtGenTotal", HistType::kTH2F, {multAxis, ptJetGenAxisUnfold});
690+
registryMC.add("hJetPtGenMatched", "hJetPtGenMatched", HistType::kTH2F, {multAxis, ptJetGenAxisUnfold});
691+
692+
// Purity components (RECO level)
693+
// Assi: X = Centrality, Y = pT_reco
694+
registryMC.add("hJetPtRecoTotal", "hJetPtRecoTotal", HistType::kTH2F, {multAxis, ptJetRecoAxisUnfold});
695+
registryMC.add("hJetPtRecoMatched", "hJetPtRecoMatched", HistType::kTH2F, {multAxis, ptJetRecoAxisUnfold});
696+
}
672697
}
673698

674699
// Delta phi calculation
@@ -799,9 +824,11 @@ struct StrangenessInJetsIons {
799824
return std::sqrt(deltaEta * deltaEta + deltaPhi * deltaPhi);
800825
}
801826

827+
// selProcess = 0 (data), = 1 (MC GEN), = 2 (MC REC)
802828
void computeRandomConeDeltaPt(const std::vector<fastjet::PseudoJet>& fjParticles,
803829
const std::vector<fastjet::PseudoJet>& jets,
804-
float multiplicity, double rhoPerp)
830+
float multiplicity, double rhoPerp,
831+
int selProcess = 0)
805832
{
806833
// Generate eta and phi for random cone in acceptance region
807834
double randomConeEta = fRng.Uniform(configTracks.etaMin + rJet, configTracks.etaMax - rJet);
@@ -843,8 +870,17 @@ struct StrangenessInJetsIons {
843870
double coneArea = TMath::Pi() * rJet * rJet;
844871
double deltaPtRandomCone = randomConePt - (coneArea * rhoPerp);
845872

846-
registryData.fill(HIST("h2_centrality_deltaPt_RandomCone"), multiplicity, deltaPtRandomCone);
847-
registryData.fill(HIST("h2_centrality_rhoPerp"), multiplicity, rhoPerp);
873+
if (selProcess == 0) {
874+
registryData.fill(HIST("h2_centrality_deltaPt_RandomCone"), multiplicity, deltaPtRandomCone);
875+
registryData.fill(HIST("h2_centrality_rhoPerp"), multiplicity, rhoPerp);
876+
} else if (selProcess == 1) {
877+
// Ricordati di definire questi istogrammi nel tuo book/registry MC!
878+
registryMC.fill(HIST("h2_centrality_deltaPt_RandomCone_gen"), multiplicity, deltaPtRandomCone);
879+
registryMC.fill(HIST("h2_centrality_rhoPerp_gen"), multiplicity, rhoPerp);
880+
} else if (selProcess == 2) {
881+
registryMC.fill(HIST("h2_centrality_deltaPt_RandomCone_rec"), multiplicity, deltaPtRandomCone);
882+
registryMC.fill(HIST("h2_centrality_rhoPerp_rec"), multiplicity, rhoPerp);
883+
}
848884
}
849885

850886
// Find ITS hit
@@ -2149,6 +2185,48 @@ struct StrangenessInJetsIons {
21492185
}
21502186
}
21512187

2188+
void ProcessJetMatchingFactorized(const std::vector<fastjet::PseudoJet>& jetsGen,
2189+
const std::vector<fastjet::PseudoJet>& jetsReco,
2190+
float centrality)
2191+
{
2192+
const double maxDeltaR = 0.24;
2193+
std::vector<bool> isRecoJetMatched(jetsReco.size(), false);
2194+
2195+
for (const auto& jetGen : jetsGen) {
2196+
registryMC.fill(HIST("hJetPtGenTotal"), centrality, jetGen.pt());
2197+
}
2198+
2199+
for (const auto& jetReco : jetsReco) {
2200+
registryMC.fill(HIST("hJetPtRecoTotal"), centrality, jetReco.pt());
2201+
}
2202+
2203+
// --- Geometrical matching reco <----> gen
2204+
for (const auto& jetGen : jetsGen) {
2205+
int bestRecoIdx = -1;
2206+
double minDeltaR = maxDeltaR;
2207+
2208+
// Search closest jet RECO in (eta,phi) space
2209+
for (long unsigned int iReco = 0; iReco < jetsReco.size(); ++iReco) {
2210+
if (isRecoJetMatched[iReco])
2211+
continue;
2212+
2213+
double dR = jetGen.delta_R(jetsReco[iReco]);
2214+
if (dR < minDeltaR) {
2215+
minDeltaR = dR;
2216+
bestRecoIdx = iReco;
2217+
}
2218+
}
2219+
2220+
if (bestRecoIdx != -1) {
2221+
isRecoJetMatched[bestRecoIdx] = true; // RECO matched
2222+
2223+
registryMC.fill(HIST("hDetectorResponse"), centrality, jetGen.pt(), jetsReco[bestRecoIdx].pt());
2224+
registryMC.fill(HIST("hJetPtGenMatched"), centrality, jetGen.pt());
2225+
registryMC.fill(HIST("hJetPtRecoMatched"), centrality, jetsReco[bestRecoIdx].pt());
2226+
}
2227+
}
2228+
}
2229+
21522230
// Process data
21532231
void processData(SelCollisions::iterator const& collision, aod::V0Datas const& fullV0s,
21542232
aod::CascDataExt const& Cascades, DaughterTracks const& tracks,
@@ -2636,6 +2714,9 @@ struct StrangenessInJetsIons {
26362714
// Estimate background energy density (rho) in perpendicular cone
26372715
auto [rhoPerp, rhoMPerp] = jetutilities::estimateRhoPerpCone(fjParticles, jets[0], rJet);
26382716

2717+
// Delta pT distributions with random cone technique
2718+
computeRandomConeDeltaPt(fjParticles, jets, genMultiplicity, rhoPerp, 1);
2719+
26392720
// Loop over clustered jets
26402721
int countSelJet = 0; // number of selected jets
26412722
for (const auto& jet : jets) {
@@ -2989,6 +3070,9 @@ struct StrangenessInJetsIons {
29893070
std::vector<fastjet::PseudoJet> jets = fastjet::sorted_by_pt(cs.inclusive_jets());
29903071
auto [rhoPerp, rhoMPerp] = jetutilities::estimateRhoPerpCone(fjParticles, jets[0], rJet);
29913072

3073+
// Delta pT distributions with random cone technique
3074+
computeRandomConeDeltaPt(fjParticles, jets, multiplicity, rhoPerp, 2);
3075+
29923076
// Jet selection
29933077
bool isAtLeastOneJetSelected = false;
29943078
std::vector<double> jetPt;
@@ -3868,6 +3952,146 @@ struct StrangenessInJetsIons {
38683952
}
38693953
}
38703954
PROCESS_SWITCH(StrangenessInJetsIons, processDataMB, "Process data in minimum bias events", false);
3955+
3956+
void processPrepareUnfolding(soa::Join<aod::McCollisions, aod::McCentFT0Ms, aod::McCentFT0Cs> const& genCollisions,
3957+
SimCollisions const& recoCollisions,
3958+
DaughterTracksMC const& mcTracks, soa::Join<aod::V0Datas, aod::McV0Labels> const& fullV0s,
3959+
aod::McParticles const& mcParticles)
3960+
{
3961+
fastjet::JetDefinition jetDef(fastjet::antikt_algorithm, rJet);
3962+
fastjet::AreaDefinition areaDef(fastjet::active_area, fastjet::GhostedAreaSpec(1.0));
3963+
3964+
// Loop over reco collisions
3965+
for (const auto& recoColl : recoCollisions) {
3966+
if (!recoColl.has_mcCollision()) {
3967+
continue;
3968+
}
3969+
3970+
// Event selections
3971+
if (!recoColl.sel8())
3972+
continue;
3973+
if (std::fabs(recoColl.posZ()) > zVtx)
3974+
continue;
3975+
if (requireNoSameBunchPileup && !recoColl.selection_bit(o2::aod::evsel::kNoSameBunchPileup))
3976+
continue;
3977+
if (requireGoodZvtxFT0vsPV && !recoColl.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV))
3978+
continue;
3979+
3980+
// Centrality
3981+
const auto& mcCollision = recoColl.mcCollision_as<soa::Join<aod::McCollisions, aod::McCentFT0Ms, aod::McCentFT0Cs>>();
3982+
float centrality;
3983+
if (centrEstimator == 0) {
3984+
centrality = mcCollision.centFT0C();
3985+
} else {
3986+
centrality = mcCollision.centFT0M();
3987+
}
3988+
3989+
std::vector<fastjet::PseudoJet> jetsGenThisEvent;
3990+
std::vector<fastjet::PseudoJet> jetsRecoThisEvent;
3991+
3992+
// ---- Extract reconstructed jets ----
3993+
// Number of V0 and cascades per collision
3994+
auto v0sPerColl = fullV0s.sliceBy(perCollisionV0, recoColl.globalIndex());
3995+
auto tracksPerColl = mcTracks.sliceBy(perCollisionTrk, recoColl.globalIndex());
3996+
const auto& mcParticlesPerColl = mcParticles.sliceBy(perMCCollision, mcCollision.globalIndex());
3997+
3998+
// Vertex position vector
3999+
TVector3 vtxPos(recoColl.posX(), recoColl.posY(), recoColl.posZ());
4000+
std::vector<fastjet::PseudoJet> fjParticlesReco;
4001+
std::vector<std::decay_t<decltype(*tracksPerColl.begin())>> fjTracks;
4002+
for (auto const& track : tracksPerColl) {
4003+
if (!passedTrackSelectionForJetReconstruction(track))
4004+
continue;
4005+
4006+
fastjet::PseudoJet fourMomentum(track.px(), track.py(), track.pz(), track.energy(o2::constants::physics::MassPionCharged));
4007+
fjParticlesReco.emplace_back(fourMomentum);
4008+
fjTracks.push_back(track);
4009+
}
4010+
4011+
// Include V0s as tracks for jet reconstruction
4012+
if (useV0inJetRec && particleOfInterestDict[ParticleOfInterest::kV0Particles]) {
4013+
AddV0sForJetReconstructionMCD(fjParticlesReco, fjTracks, v0sPerColl, mcParticles, vtxPos);
4014+
}
4015+
fjTracks.clear();
4016+
4017+
if (fjParticlesReco.size() >= 1) { // Reject empty events
4018+
fastjet::ClusterSequenceArea csReco(fjParticlesReco, jetDef, areaDef);
4019+
std::vector<fastjet::PseudoJet> recoJets = fastjet::sorted_by_pt(csReco.inclusive_jets());
4020+
if (recoJets.empty())
4021+
continue;
4022+
auto [rhoPerp, rhoMPerp] = jetutilities::estimateRhoPerpCone(fjParticlesReco, recoJets[0], rJet);
4023+
4024+
for (const auto& jet : recoJets) {
4025+
// jet must be fully contained in the acceptance
4026+
if ((std::fabs(jet.eta()) + rJet) > (configTracks.etaMax - deltaEtaEdge))
4027+
continue;
4028+
4029+
// jet pt must be larger than threshold
4030+
auto jetForSub = jet;
4031+
fastjet::PseudoJet jetMinusBkg = backgroundSub.doRhoAreaSub(jetForSub, rhoPerp, rhoMPerp);
4032+
if (jetMinusBkg.pt() < minJetPt)
4033+
continue;
4034+
4035+
jetsRecoThisEvent.push_back(jetMinusBkg);
4036+
}
4037+
}
4038+
// ----------------------------------------
4039+
4040+
// ---- Extract generated jets ----
4041+
// LOG(info) << "recoCollisions.globalIndex()" << recoColl.globalIndex();
4042+
// LOG(info) << "recoColl.mcCollisionId()" << recoColl.mcCollisionId();
4043+
// LOG(info) << "mcCollision.globalIndex()" << mcCollision.globalIndex();
4044+
std::vector<fastjet::PseudoJet> fjParticlesGen;
4045+
std::vector<std::decay_t<decltype(*mcParticlesPerColl.begin())>> fjParticleObj;
4046+
4047+
for (const auto& particle : mcParticlesPerColl) {
4048+
double minPtParticle = 0.1f;
4049+
if (particle.eta() < configTracks.etaMin || particle.eta() > configTracks.etaMax || particle.pt() < minPtParticle)
4050+
continue;
4051+
4052+
// Select physical primary particles or HF decay products
4053+
if (!isPhysicalPrimaryOrFromHF(particle, mcParticles))
4054+
continue;
4055+
4056+
// Build 4-momentum assuming charged pion mass
4057+
static constexpr float kMassPionChargedSquared = o2::constants::physics::MassPionCharged * o2::constants::physics::MassPionCharged;
4058+
const double energy = std::sqrt(particle.p() * particle.p() + kMassPionChargedSquared);
4059+
fastjet::PseudoJet fourMomentum(particle.px(), particle.py(), particle.pz(), energy);
4060+
fourMomentum.set_user_index(particle.pdgCode());
4061+
fjParticlesGen.emplace_back(fourMomentum);
4062+
fjParticleObj.push_back(particle);
4063+
}
4064+
4065+
if (useV0inJetRec && particleOfInterestDict[ParticleOfInterest::kV0Particles]) {
4066+
AddV0sForJetReconstructionMCP(fjParticlesGen, fjParticleObj, mcParticlesPerColl, mcParticles);
4067+
}
4068+
4069+
if (fjParticlesGen.size() >= 1) { // Skip events with no particles
4070+
fastjet::ClusterSequenceArea csGen(fjParticlesGen, jetDef, areaDef);
4071+
std::vector<fastjet::PseudoJet> genJets = fastjet::sorted_by_pt(csGen.inclusive_jets());
4072+
if (genJets.empty())
4073+
continue;
4074+
auto [rhoPerp, rhoMPerp] = jetutilities::estimateRhoPerpCone(fjParticlesGen, genJets[0], rJet);
4075+
4076+
for (const auto& jet : genJets) {
4077+
if ((std::fabs(jet.eta()) + rJet) > (configTracks.etaMax - deltaEtaEdge))
4078+
continue;
4079+
4080+
auto jetForSub = jet;
4081+
fastjet::PseudoJet jetMinusBkg = backgroundSub.doRhoAreaSub(jetForSub, rhoPerp, rhoMPerp);
4082+
if (jetMinusBkg.pt() < minJetPt)
4083+
continue;
4084+
4085+
jetsGenThisEvent.push_back(jetMinusBkg);
4086+
}
4087+
}
4088+
// ----------------------------------------
4089+
4090+
// Fill histograms for unfolding
4091+
ProcessJetMatchingFactorized(jetsGenThisEvent, jetsRecoThisEvent, centrality);
4092+
}
4093+
}
4094+
PROCESS_SWITCH(StrangenessInJetsIons, processPrepareUnfolding, "Process to prepare hsitograms for the unfolding procedure", false);
38714095
};
38724096

38734097
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)

0 commit comments

Comments
 (0)