Skip to content

Commit 3ac5680

Browse files
FrancescaCasilloFrancesca Casilloalibuild
authored
[PWGLF] Add coalescence for corr analysis (#14639)
Co-authored-by: Francesca Casillo <francesca@MacBook-Air-di-Francesca-3.local> Co-authored-by: ALICE Action Bot <alibuild@cern.ch>
1 parent f1a1176 commit 3ac5680

File tree

1 file changed

+238
-3
lines changed

1 file changed

+238
-3
lines changed

‎PWGLF/Tasks/Nuspex/antinucleiInJets.cxx‎

Lines changed: 238 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -20,6 +20,7 @@
2020
#include "PWGJE/Core/JetUtilities.h"
2121
#include "PWGJE/DataModel/Jet.h"
2222
#include "PWGJE/DataModel/JetReducedData.h"
23+
#include "PWGLF/DataModel/mcCentrality.h"
2324

2425
#include "Common/Core/TrackSelection.h"
2526
#include "Common/Core/Zorro.h"
@@ -34,7 +35,6 @@
3435

3536
#include "CCDB/BasicCCDBManager.h"
3637
#include "CCDB/CcdbApi.h"
37-
#include "MathUtils/BetheBlochAleph.h"
3838
#include "Framework/ASoA.h"
3939
#include "Framework/ASoAHelpers.h"
4040
#include "Framework/AnalysisDataModel.h"
@@ -44,6 +44,7 @@
4444
#include "Framework/Logger.h"
4545
#include "Framework/RunningWorkflowInfo.h"
4646
#include "Framework/runDataProcessing.h"
47+
#include "MathUtils/BetheBlochAleph.h"
4748
#include "ReconstructionDataFormats/DCA.h"
4849
#include "ReconstructionDataFormats/PID.h"
4950
#include "ReconstructionDataFormats/Track.h"
@@ -89,10 +90,9 @@ using std::array;
8990
// Define convenient aliases for commonly used table joins
9091
using SelectedCollisions = soa::Join<aod::Collisions, aod::EvSels, aod::CentFT0Ms>;
9192
using RecCollisionsMc = soa::Join<aod::Collisions, aod::EvSels, aod::CentFT0Ms, aod::McCollisionLabels>;
92-
using GenCollisionsMc = aod::McCollisions;
93+
using GenCollisionsMc = soa::Join<aod::McCollisions, aod::McCentFT0Ms>;
9394
using AntiNucleiTracks = soa::Join<aod::Tracks, aod::TracksExtra, aod::TrackSelection, aod::TrackSelectionExtension, aod::TracksDCA, aod::pidTPCFullPr, aod::pidTPCFullDe, aod::pidTPCFullHe, aod::pidTOFFullPr, aod::pidTOFFullDe, aod::pidTOFFullHe>;
9495
using AntiNucleiTracksMc = soa::Join<AntiNucleiTracks, aod::McTrackLabels>;
95-
9696
using LorentzVector = ROOT::Math::PxPyPzEVector;
9797

9898
// Lightweight particle container for fast kinematic access
@@ -489,6 +489,38 @@ struct AntinucleiInJets {
489489
registryMC.add("antiproton_coal_ue", "antiproton_coal_ue", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}});
490490
}
491491

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,6 +3629,209 @@ 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+
} else if (absPdg == PDG_t::kNeutron) { // Store Neutrons
3695+
neutronCandidates.push_back({particle.px(), particle.py(), particle.pz(), pdg, mcId, false});
3696+
}
3697+
}
3698+
}
3699+
3700+
// Reject empty events
3701+
if (protonCandidates.empty() && neutronCandidates.empty())
3702+
continue;
3703+
3704+
registryMC.fill(HIST("genEventsCoalescenceCorr"), 2.5);
3705+
3706+
// Build deuterons
3707+
for (size_t iP = 0; iP < protonCandidates.size(); ++iP) {
3708+
if (protonCandidates[iP].used)
3709+
continue;
3710+
3711+
for (size_t iN = 0; iN < neutronCandidates.size(); ++iN) {
3712+
if (neutronCandidates[iN].used)
3713+
continue;
3714+
3715+
// Physics consistency check
3716+
if (protonCandidates[iP].pdgCode * neutronCandidates[iN].pdgCode < 0)
3717+
continue;
3718+
3719+
if (passDeuteronCoalescence(protonCandidates[iP], neutronCandidates[iN], coalescenceMomentum, mRand)) {
3720+
neutronCandidates[iN].used = true;
3721+
protonCandidates[iP].used = true;
3722+
3723+
int sign = (protonCandidates[iP].pdgCode > 0) ? +1 : -1;
3724+
int deuteronPdg = sign * o2::constants::physics::Pdg::kDeuteron;
3725+
3726+
double pxDeut = protonCandidates[iP].px + neutronCandidates[iN].px;
3727+
double pyDeut = protonCandidates[iP].py + neutronCandidates[iN].py;
3728+
double pzDeut = protonCandidates[iP].pz + neutronCandidates[iN].pz;
3729+
double energyDeut = std::sqrt(pxDeut * pxDeut + pyDeut * pyDeut + pzDeut * pzDeut + massDeut * massDeut);
3730+
LorentzVector pd(pxDeut, pyDeut, pzDeut, energyDeut);
3731+
if (pd.Eta() >= minEta && pd.Eta() <= maxEta && (0.5 * pd.Pt()) >= MinPtPerNucleon) {
3732+
// Store Deuteron
3733+
finalDeuterons.push_back({pxDeut, pyDeut, pzDeut, deuteronPdg, protonCandidates[iP].mcIndex, false});
3734+
}
3735+
3736+
break;
3737+
}
3738+
}
3739+
}
3740+
3741+
// Add unused protons to final vectors
3742+
for (const auto& proton : protonCandidates) {
3743+
if (!proton.used) {
3744+
finalProtons.push_back(proton);
3745+
}
3746+
}
3747+
3748+
// Correlation Analysis
3749+
std::vector<int> nAntiprotonFullEvent(nBins, 0);
3750+
std::vector<int> nAntideuteronFullEvent(nBins, 0);
3751+
int nTotProtonFullEvent(0);
3752+
int nTotDeuteronFullEvent(0);
3753+
int nTotAntiprotonFullEvent(0);
3754+
int nTotAntideuteronFullEvent(0);
3755+
3756+
// Loop over final protons
3757+
for (const auto& part : finalProtons) {
3758+
double pt = std::hypot(part.px, part.py);
3759+
3760+
if (part.eta() < minEta || part.eta() > maxEta)
3761+
continue;
3762+
3763+
// Standard histograms for antiprotons
3764+
if (part.pdgCode == PDG_t::kProtonBar) {
3765+
registryMC.fill(HIST("antiproton_fullEvent_CoalescenceCorr"), pt);
3766+
}
3767+
3768+
// Kinematic selection and Multiplicity counting
3769+
if (pt < ptOverAbins[0] || pt >= ptOverAbins[nBins])
3770+
continue;
3771+
3772+
if (part.pdgCode > 0) {
3773+
nTotProtonFullEvent++;
3774+
} else {
3775+
nTotAntiprotonFullEvent++;
3776+
int ibin = findBin(ptOverAbins, pt);
3777+
if (ibin >= 0 && ibin < nBins)
3778+
nAntiprotonFullEvent[ibin]++;
3779+
}
3780+
}
3781+
3782+
// Loop over final deuterons
3783+
for (const auto& part : finalDeuterons) {
3784+
double pt = std::hypot(part.px, part.py);
3785+
double ptPerNucleon = 0.5 * pt;
3786+
3787+
// Apply detector acceptance cuts (to match real data)
3788+
if (part.eta() < minEta || part.eta() > maxEta)
3789+
continue;
3790+
3791+
// Standard histograms for antideuterons
3792+
if (part.pdgCode == -o2::constants::physics::Pdg::kDeuteron) {
3793+
registryMC.fill(HIST("antideuteron_fullEvent_CoalescenceCorr"), pt);
3794+
}
3795+
3796+
// Kinematic selection and Multiplicity counting
3797+
if (ptPerNucleon < ptOverAbins[0] || ptPerNucleon >= ptOverAbins[nBins])
3798+
continue;
3799+
3800+
if (part.pdgCode > 0) {
3801+
nTotDeuteronFullEvent++;
3802+
} else {
3803+
nTotAntideuteronFullEvent++;
3804+
int ibin = findBin(ptOverAbins, ptPerNucleon);
3805+
if (ibin >= 0 && ibin < nBins)
3806+
nAntideuteronFullEvent[ibin]++;
3807+
}
3808+
}
3809+
3810+
// Fill correlation histograms
3811+
int netProtonFullEvent = nTotProtonFullEvent - nTotAntiprotonFullEvent;
3812+
int netDeuteronFullEvent = nTotDeuteronFullEvent - nTotAntideuteronFullEvent;
3813+
3814+
registryCorr.fill(HIST("rho_fullEvent_CoalescenceCorr"), nTotAntideuteronFullEvent, nTotAntiprotonFullEvent, multiplicity);
3815+
registryCorr.fill(HIST("rho_netP_netD_fullEvent_CoalescenceCorr"), netDeuteronFullEvent, netProtonFullEvent);
3816+
3817+
// Fill efficiency histograms
3818+
for (int i = 0; i < nBins; i++) {
3819+
double ptAcenteri = 0.5 * (ptOverAbins[i] + ptOverAbins[i + 1]);
3820+
3821+
registryCorr.fill(HIST("q1d_fullEvent_CoalescenceCorr"), nAntideuteronFullEvent[i], ptAcenteri, multiplicity);
3822+
registryCorr.fill(HIST("q1p_fullEvent_CoalescenceCorr"), nAntiprotonFullEvent[i], ptAcenteri, multiplicity);
3823+
3824+
for (int j = 0; j < nBins; j++) {
3825+
double ptAcenterj = 0.5 * (ptOverAbins[j] + ptOverAbins[j + 1]);
3826+
3827+
registryCorr.fill(HIST("q1d_square_fullEvent_CoalescenceCorr"), ptAcenteri, ptAcenterj, static_cast<double>(nAntideuteronFullEvent[i] * nAntideuteronFullEvent[j]), multiplicity);
3828+
registryCorr.fill(HIST("q1p_square_fullEvent_CoalescenceCorr"), ptAcenteri, ptAcenterj, static_cast<double>(nAntiprotonFullEvent[i] * nAntiprotonFullEvent[j]), multiplicity);
3829+
registryCorr.fill(HIST("q1d_q1p_fullEvent_CoalescenceCorr"), ptAcenteri, ptAcenterj, static_cast<double>(nAntideuteronFullEvent[i] * nAntiprotonFullEvent[j]), multiplicity);
3830+
}
3831+
}
3832+
}
3833+
}
3834+
PROCESS_SWITCH(AntinucleiInJets, processCoalescenceCorr, "process coalescence correlation", false);
36003835
};
36013836

36023837
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)

0 commit comments

Comments
 (0)