Skip to content

Commit bac751c

Browse files
Francesca CasilloFrancesca Casillo
Francesca Casillo
authored and
Francesca Casillo
committed
[PWGLF] Add coalescence for corr analysis
1 parent 8f45b4c commit bac751c

File tree

1 file changed

+232
-1
lines changed

1 file changed

+232
-1
lines changed

‎PWGLF/Tasks/Nuspex/antinucleiInJets.cxx‎

Lines changed: 232 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -488,7 +488,39 @@ struct AntinucleiInJets {
488488
registryMC.add("antiproton_coal_jet", "antiproton_coal_jet", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}});
489489
registryMC.add("antiproton_coal_ue", "antiproton_coal_ue", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}});
490490
}
491-
491+
492+
// Coalescence and Correlation analysis
493+
if (doprocessCoalescenceCorr) {
494+
495+
// Axes definitions for multidimensional histogram binning
496+
const AxisSpec multiplicityAxis{100, 0.0, 100.0, "multiplicity percentile"};
497+
const AxisSpec ptPerNucleonAxis{5, 0.4, 0.9, "{p}_{T}/A (GeV/#it{c})"};
498+
const AxisSpec nAntideuteronsAxis{10, 0.0, 10.0, "N_{#bar{d}}"};
499+
const AxisSpec nAntiprotonsAxis{10, 0.0, 10.0, "N_{#bar{p}}"};
500+
const AxisSpec nBarD2Axis{100, 0.0, 100.0, "N_{#bar{d}}^{i} #times N_{#bar{d}}^{j}"};
501+
const AxisSpec nBarP2Axis{100, 0.0, 100.0, "N_{#bar{p}}^{i} #times N_{#bar{p}}^{j}"};
502+
const AxisSpec nBarDnBarPAxis{100, 0.0, 100.0, "N_{#bar{d}}^{i} #times N_{#bar{p}}^{j}"};
503+
504+
registryMC.add("genEventsCoalescenceCorr", "genEventsCoalescenceCorr", HistType::kTH1F, {{20, 0, 20, "counter"}});
505+
registryMC.add("antideuteron_fullEvent_CoalescenceCorr", "antideuteron_fullEvent_CoalescenceCorr", HistType::kTH1F, {{nbins, 2 * min, 2 * max, "#it{p}_{T} (GeV/#it{c})"}});
506+
registryMC.add("antiproton_fullEvent_CoalescenceCorr", "antiproton_fullEvent_CoalescenceCorr", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}});
507+
508+
//Counter histograms
509+
registryCorr.add("eventCounter_CoalescenceCorr", "number of events in Coalescence simulation", HistType::kTH1F, {{20, 0, 20, "counter"}});
510+
registryCorr.add("eventCounter_centrality_fullEvent_CoalescenceCorr", "Number of events per centrality (Full Event) in Coalescence simulation", HistType::kTH1F, {multiplicityAxis});
511+
512+
// Correlation histograms
513+
registryCorr.add("rho_fullEvent_CoalescenceCorr", "rho_fullEvent_CoalescenceCorr", HistType::kTH3F, {nAntideuteronsAxis, nAntiprotonsAxis, multiplicityAxis});
514+
registryCorr.add("rho_netP_netD_fullEvent_CoalescenceCorr", "rho_netP_netD_fullEvent_CoalescenceCorr", HistType::kTH2F, {nAntideuteronsAxis, nAntiprotonsAxis});
515+
516+
// Efficiency histograms full event
517+
registryCorr.add("q1d_fullEvent_CoalescenceCorr", "q1d_fullEvent_CoalescenceCorr", HistType::kTH3F, {nAntideuteronsAxis, ptPerNucleonAxis, multiplicityAxis});
518+
registryCorr.add("q1p_fullEvent_CoalescenceCorr", "q1p_fullEvent_CoalescenceCorr", HistType::kTH3F, {nAntiprotonsAxis, ptPerNucleonAxis, multiplicityAxis});
519+
registryCorr.add("q1d_square_fullEvent_CoalescenceCorr", "q1d_square_fullEvent_CoalescenceCorr", HistType::kTHnSparseD, {ptPerNucleonAxis, ptPerNucleonAxis, nBarD2Axis, multiplicityAxis});
520+
registryCorr.add("q1p_square_fullEvent_CoalescenceCorr", "q1p_square_fullEvent_CoalescenceCorr", HistType::kTHnSparseD, {ptPerNucleonAxis, ptPerNucleonAxis, nBarP2Axis, multiplicityAxis});
521+
registryCorr.add("q1d_q1p_fullEvent_CoalescenceCorr", "q1d_q1p_fullEvent_CoalescenceCorr", HistType::kTHnSparseD, {ptPerNucleonAxis, ptPerNucleonAxis, nBarDnBarPAxis, multiplicityAxis});
522+
}
523+
492524
// Systematic uncertainties (Data)
493525
if (doprocessSystData) {
494526
registryData.add("number_of_events_data_syst", "event counter", HistType::kTH1F, {{20, 0, 20, "counter"}});
@@ -3597,9 +3629,208 @@ struct AntinucleiInJets {
35973629
}
35983630
}
35993631
PROCESS_SWITCH(AntinucleiInJets, processCoalescence, "process coalescence", false);
3632+
3633+
// process Coalescence and Correlation Analysis
3634+
void processCoalescenceCorr(GenCollisionsMc const& collisions, aod::McParticles const& mcParticles)
3635+
{
3636+
// Deuteron Mass and minimum pt
3637+
double massDeut = o2::constants::physics::MassDeuteron;
3638+
static constexpr double MinPtPerNucleon = 0.1; // Cut on pt/A
3639+
3640+
// Containers for candidates (before coalescence)
3641+
std::vector<ReducedParticle> protonCandidates;
3642+
std::vector<ReducedParticle> neutronCandidates;
3643+
3644+
// Final containers for analysis (after coalescence)
3645+
std::vector<ReducedParticle> finalProtons;
3646+
std::vector<ReducedParticle> finalDeuterons;
3647+
3648+
// pt/A bins
3649+
std::vector<double> ptOverAbins = {0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0};
3650+
const int nBins = ptOverAbins.size() - 1;
3651+
3652+
// Loop over all simulated collisions
3653+
for (const auto& collision : collisions) {
3654+
3655+
// Clear containers
3656+
protonCandidates.clear();
3657+
neutronCandidates.clear();
3658+
finalProtons.clear();
3659+
finalDeuterons.clear();
3660+
3661+
// Event counter: before event selection
3662+
registryCorr.fill(HIST("eventCounter_CoalescenceCorr"), 0.5);
3663+
registryMC.fill(HIST("genEventsCoalescenceCorr"), 0.5);
3664+
3665+
// Apply event selection: require vertex position to be within the allowed z range
3666+
if (std::fabs(collision.posZ()) > zVtx)
3667+
continue;
3668+
3669+
// Event counter: after event selection
3670+
registryCorr.fill(HIST("eventCounter_CoalescenceCorr"), 1.5);
3671+
registryMC.fill(HIST("genEventsCoalescenceCorr"), 1.5);
3672+
3673+
// Multiplicity percentile
3674+
const float multiplicity = collision.centFT0M();
3675+
3676+
// Fill event counter vs centrality
3677+
registryCorr.fill(HIST("eventCounter_centrality_fullEvent_CoalescenceCorr"), multiplicity);
3678+
3679+
// Get particles in this MC collision
3680+
const auto mcParticlesThisMcColl = mcParticles.sliceBy(mcParticlesPerMcCollision, collision.globalIndex());
3681+
3682+
// Loop over MC particles
3683+
for (const auto& particle : mcParticlesThisMcColl) {
3684+
3685+
// Monte Carlo index
3686+
int mcId = particle.globalIndex();
3687+
int pdg = particle.pdgCode();
3688+
int absPdg = std::abs(pdg);
3689+
3690+
// Store Protons
3691+
if (particle.isPhysicalPrimary()) {
3692+
if (absPdg == PDG_t::kProton) {
3693+
protonCandidates.push_back({particle.px(), particle.py(), particle.pz(), pdg, mcId, false});
3694+
}
3695+
else if (absPdg == PDG_t::kNeutron) { // Store Neutrons
3696+
neutronCandidates.push_back({particle.px(), particle.py(), particle.pz(), pdg, mcId, false});
3697+
}
3698+
}
3699+
}
3700+
3701+
// Reject empty events
3702+
if (protonCandidates.empty() && neutronCandidates.empty())
3703+
continue;
3704+
3705+
registryMC.fill(HIST("genEventsCoalescenceCorr"), 2.5);
3706+
3707+
// Build deuterons
3708+
for (auto& proton : protonCandidates) {
3709+
if (proton.used) continue;
3710+
3711+
for (auto& neutron : neutronCandidates) {
3712+
if (neutron.used) continue;
3713+
3714+
// Physics consistency check
3715+
if (proton.pdgCode * neutron.pdgCode < 0) continue;
3716+
3717+
if (passDeuteronCoalescence(proton, neutron, coalescenceMomentum, mRand)) {
3718+
3719+
neutron.used = true;
3720+
proton.used = true;
3721+
3722+
int sign = (proton.pdgCode > 0) ? +1 : -1;
3723+
int deuteronPdg = sign * o2::constants::physics::Pdg::kDeuteron;
3724+
3725+
double pxDeut = proton.px + neutron.px;
3726+
double pyDeut = proton.py + neutron.py;
3727+
double pzDeut = proton.pz + neutron.pz;
3728+
double energyDeut = std::sqrt(pxDeut * pxDeut + pyDeut * pyDeut + pzDeut * pzDeut + massDeut * massDeut);
3729+
LorentzVector pd(pxDeut, pyDeut, pzDeut, energyDeut);
3730+
if (pd.Eta() >= minEta && pd.Eta() <= maxEta && (0.5 * pd.Pt()) >= MinPtPerNucleon) {
3731+
// Store Deuteron
3732+
finalDeuterons.push_back({pxDeut, pyDeut, pzDeut, deuteronPdg, proton.mcIndex, false});
3733+
}
3734+
3735+
break;
3736+
}
3737+
}
3738+
}
3739+
3740+
// Add unused protons to final vectors
3741+
for (const auto& proton : protonCandidates) {
3742+
if (!proton.used) {
3743+
finalProtons.push_back(proton);
3744+
}
3745+
}
3746+
3747+
// Correlation Analysis
3748+
std::vector<int> nAntiprotonFullEvent(nBins, 0);
3749+
std::vector<int> nAntideuteronFullEvent(nBins, 0);
3750+
int nTotProtonFullEvent(0);
3751+
int nTotDeuteronFullEvent(0);
3752+
int nTotAntiprotonFullEvent(0);
3753+
int nTotAntideuteronFullEvent(0);
3754+
3755+
// Loop over final protons
3756+
for (const auto& part : finalProtons) {
3757+
double pt = std::hypot(part.px, part.py);
3758+
3759+
if (part.eta() < minEta || part.eta() > maxEta)
3760+
continue;
3761+
3762+
// Standard histograms for antiprotons
3763+
if (part.pdgCode == PDG_t::kProtonBar) {
3764+
registryMC.fill(HIST("antiproton_fullEvent_CoalescenceCorr"), pt);
3765+
}
3766+
3767+
// Kinematic selection and Multiplicity counting
3768+
if (pt < ptOverAbins[0] || pt >= ptOverAbins[nBins]) continue;
3769+
3770+
if (part.pdgCode > 0) {
3771+
nTotProtonFullEvent++;
3772+
} else {
3773+
nTotAntiprotonFullEvent++;
3774+
int ibin = findBin(ptOverAbins, pt);
3775+
if (ibin >= 0 && ibin < nBins) nAntiprotonFullEvent[ibin]++;
3776+
}
3777+
}
3778+
3779+
// Loop over final deuterons
3780+
for (const auto& part : finalDeuterons) {
3781+
double pt = std::hypot(part.px, part.py);
3782+
double ptPerNucleon = 0.5 * pt;
3783+
3784+
// Apply detector acceptance cuts (to match real data)
3785+
if (part.eta() < minEta || part.eta() > maxEta)
3786+
continue;
3787+
3788+
// Standard histograms for antideuterons
3789+
if (part.pdgCode == -o2::constants::physics::Pdg::kDeuteron) {
3790+
registryMC.fill(HIST("antideuteron_fullEvent_CoalescenceCorr"), pt);
3791+
}
3792+
3793+
// Kinematic selection and Multiplicity counting
3794+
if (ptPerNucleon < ptOverAbins[0] || ptPerNucleon >= ptOverAbins[nBins]) continue;
3795+
3796+
if (part.pdgCode > 0) {
3797+
nTotDeuteronFullEvent++;
3798+
} else {
3799+
nTotAntideuteronFullEvent++;
3800+
int ibin = findBin(ptOverAbins, ptPerNucleon);
3801+
if (ibin >= 0 && ibin < nBins) nAntideuteronFullEvent[ibin]++;
3802+
}
3803+
}
3804+
3805+
// Fill correlation histograms
3806+
int netProtonFullEvent = nTotProtonFullEvent - nTotAntiprotonFullEvent;
3807+
int netDeuteronFullEvent = nTotDeuteronFullEvent - nTotAntideuteronFullEvent;
3808+
3809+
registryCorr.fill(HIST("rho_fullEvent_CoalescenceCorr"), nTotAntideuteronFullEvent, nTotAntiprotonFullEvent, multiplicity);
3810+
registryCorr.fill(HIST("rho_netP_netD_fullEvent_CoalescenceCorr"), netDeuteronFullEvent, netProtonFullEvent);
3811+
3812+
// Fill efficiency histograms
3813+
for (int i = 0; i < nBins; i++) {
3814+
double ptAcenteri = 0.5 * (ptOverAbins[i] + ptOverAbins[i + 1]);
3815+
3816+
registryCorr.fill(HIST("q1d_fullEvent_CoalescenceCorr"), nAntideuteronFullEvent[i], ptAcenteri, multiplicity);
3817+
registryCorr.fill(HIST("q1p_fullEvent_CoalescenceCorr"), nAntiprotonFullEvent[i], ptAcenteri, multiplicity);
3818+
3819+
for (int j = 0; j < nBins; j++) {
3820+
double ptAcenterj = 0.5 * (ptOverAbins[j] + ptOverAbins[j + 1]);
3821+
3822+
registryCorr.fill(HIST("q1d_square_fullEvent_CoalescenceCorr"), ptAcenteri, ptAcenterj, (double)(nAntideuteronFullEvent[i] * nAntideuteronFullEvent[j]), multiplicity);
3823+
registryCorr.fill(HIST("q1p_square_fullEvent_CoalescenceCorr"), ptAcenteri, ptAcenterj, (double)(nAntiprotonFullEvent[i] * nAntiprotonFullEvent[j]), multiplicity);
3824+
registryCorr.fill(HIST("q1d_q1p_fullEvent_CoalescenceCorr"), ptAcenteri, ptAcenterj, (double)(nAntideuteronFullEvent[i] * nAntiprotonFullEvent[j]), multiplicity);
3825+
}
3826+
}
3827+
}
3828+
}
3829+
PROCESS_SWITCH(AntinucleiInJets, processCoalescenceCorr, "process coalescence correlation", false);
36003830
};
36013831

36023832
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
36033833
{
36043834
return WorkflowSpec{adaptAnalysisTask<AntinucleiInJets>(cfgc)};
36053835
}
3836+

0 commit comments

Comments
 (0)