Skip to content

[PWGLF] added DCA templates and function to compute multiplicity in different event classes #11768

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
alcaliva merged 2 commits into from
Jun 25, 2025
Merged

[PWGLF] added DCA templates and function to compute multiplicity in different event classes #11768

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
ffa4435
[PWGLF] added DCA templates and function to compute multiplicity in d…
alcaliva Jun 25, 2025
e86773a
fixed clang format: 1 try
alcaliva Jun 25, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
120 changes: 109 additions & 11 deletions PWGLF/Tasks/Nuspex/antinucleiInJets.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -77,6 +77,7 @@ using std::array;

using SelectedCollisions = soa::Join<aod::Collisions, aod::EvSels>;
using SimCollisions = soa::Join<aod::Collisions, aod::EvSels, aod::McCollisionLabels>;
using GenCollisions = aod::McCollisions;

using FullNucleiTracks = soa::Join<aod::Tracks, aod::TracksExtra, aod::TrackSelection, aod::TrackSelectionExtension, aod::TracksDCA, aod::pidTPCFullPr, aod::pidTPCFullDe, aod::pidTPCFullHe, aod::pidTOFFullPr, aod::pidTOFFullDe, aod::pidTOFFullHe>;

Expand All @@ -88,9 +89,11 @@ struct AntinucleiInJets {
HistogramRegistry registryData{"registryData", {}, OutputObjHandlingPolicy::AnalysisObject, true, true};
HistogramRegistry registryMC{"registryMC", {}, OutputObjHandlingPolicy::AnalysisObject, true, true};
HistogramRegistry registryQC{"registryQC", {}, OutputObjHandlingPolicy::AnalysisObject, true, true};
HistogramRegistry registryMult{"registryMult", {}, OutputObjHandlingPolicy::AnalysisObject, true, true};

// global parameters
Configurable<double> minJetPt{"minJetPt", 10.0, "Minimum pt of the jet"};
Configurable<double> ptLeadingMin{"ptLeadingMin", 5.0, "pt Leading Min"};
Configurable<double> rJet{"rJet", 0.3, "Jet resolution parameter R"};
Configurable<double> zVtx{"zVtx", 10.0, "Maximum zVertex"};
Configurable<double> deltaEtaEdge{"deltaEtaEdge", 0.05, "eta gap from the edge"};
Expand Down Expand Up @@ -179,6 +182,11 @@ struct AntinucleiInJets {
registryQC.add("jetPtDifference", "jetPtDifference", HistType::kTH1F, {{200, -1, 1, "#Deltap_{T}^{jet}"}});
}

if (doprocessMultEvents) {
registryMult.add("multiplicityEvtsPtLeading", "multiplicityEvtsPtLeading", HistType::kTH1F, {{1000, 0, 1000, "#it{N}_{ch}"}});
registryMult.add("multiplicityEvtsWithJet", "multiplicityEvtsWithJet", HistType::kTH1F, {{1000, 0, 1000, "#it{N}_{ch}"}});
}

// data histograms
if (doprocessData) {

Expand All @@ -191,8 +199,8 @@ struct AntinucleiInJets {
registryData.add("antiproton_jet_tof", "antiproton_jet_tof", HistType::kTH2F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}, {400, -20.0, 20.0, "n#sigma_{TOF}"}});
registryData.add("antiproton_ue_tpc", "antiproton_ue_tpc", HistType::kTH2F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}, {400, -20.0, 20.0, "n#sigma_{TPC}"}});
registryData.add("antiproton_ue_tof", "antiproton_ue_tof", HistType::kTH2F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}, {400, -20.0, 20.0, "n#sigma_{TOF}"}});
registryData.add("antiproton_dca_jet", "antiproton_dca_jet", HistType::kTH2F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}, {200, -0.5, 0.5, "DCA_{xy} (cm)"}});
registryData.add("antiproton_dca_ue", "antiproton_dca_ue", HistType::kTH2F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}, {200, -0.5, 0.5, "DCA_{xy} (cm)"}});
registryData.add("antiproton_dca_jet", "antiproton_dca_jet", HistType::kTH2F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}, {200, -1, 1, "DCA_{xy} (cm)"}});
registryData.add("antiproton_dca_ue", "antiproton_dca_ue", HistType::kTH2F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}, {200, -1, 1, "DCA_{xy} (cm)"}});

// antideuterons
registryData.add("antideuteron_jet_tpc", "antideuteron_jet_tpc", HistType::kTH2F, {{nbins, min * 2, max * 2, "#it{p}_{T} (GeV/#it{c})"}, {400, -20.0, 20.0, "n#sigma_{TPC}"}});
Expand Down Expand Up @@ -271,6 +279,12 @@ struct AntinucleiInJets {
registryMC.add("antiproton_prim_ue", "antiproton_prim_ue", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}});
registryMC.add("antiproton_incl_ue", "antiproton_incl_ue", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}});

// DCA Templates
registryData.add("antiproton_prim_dca_jet", "antiproton_prim_dca_jet", HistType::kTH2F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}, {200, -1, 1, "DCA_{xy} (cm)"}});
registryData.add("antiproton_prim_dca_ue", "antiproton_prim_dca_ue", HistType::kTH2F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}, {200, -1, 1, "DCA_{xy} (cm)"}});
registryData.add("antiproton_all_dca_jet", "antiproton_all_dca_jet", HistType::kTH2F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}, {200, -1, 1, "DCA_{xy} (cm)"}});
registryData.add("antiproton_all_dca_ue", "antiproton_all_dca_ue", HistType::kTH2F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}, {200, -1, 1, "DCA_{xy} (cm)"}});

// antiproton reweighting
registryMC.add("antiproton_eta_pt_pythia", "antiproton_eta_pt_pythia", HistType::kTH2F, {{200, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}, {20, -1.0, 1.0, "#it{#eta}"}});
registryMC.add("antiproton_eta_pt_jet", "antiproton_eta_pt_jet", HistType::kTH2F, {{200, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}, {20, -1.0, 1.0, "#it{#eta}"}});
Expand Down Expand Up @@ -597,7 +611,8 @@ struct AntinucleiInJets {
// jet pt must be larger than threshold
auto jetForSub = jet;
fastjet::PseudoJet jetMinusBkg = backgroundSub.doRhoAreaSub(jetForSub, rhoPerp, rhoMPerp);
if (getCorrectedPt(jetMinusBkg.pt(), responseMatrix) < minJetPt)
// if (getCorrectedPt(jetMinusBkg.pt(), responseMatrix) < minJetPt)
if (jetMinusBkg.pt() < minJetPt)
continue;
isAtLeastOneJetSelected = true;

Expand Down Expand Up @@ -768,6 +783,68 @@ struct AntinucleiInJets {
}
PROCESS_SWITCH(AntinucleiInJets, processData, "Process Data", true);

void processMultEvents(SelectedCollisions::iterator const& collision, FullNucleiTracks const& tracks)
{
// event selection
if (!collision.sel8() || std::fabs(collision.posZ()) > zVtx)
return;

// Leading Track
double ptMax(0.0);

// loop over reconstructed tracks
int id(-1);
std::vector<fastjet::PseudoJet> fjParticles;
for (auto const& track : tracks) {
id++;
if (!passedTrackSelectionForJetReconstruction(track))
continue;
if (track.pt() > ptMax) {
ptMax = track.pt();
}

// 4-momentum representation of a particle
fastjet::PseudoJet fourMomentum(track.px(), track.py(), track.pz(), track.energy(MassPionCharged));
fourMomentum.set_user_index(id);
fjParticles.emplace_back(fourMomentum);
}
// reject empty events
if (fjParticles.empty()) {
return;
}

if (ptMax > ptLeadingMin) {
registryMult.fill(HIST("multiplicityEvtsPtLeading"), fjParticles.size());
}

// cluster particles using the anti-kt algorithm
fastjet::JetDefinition jetDef(fastjet::antikt_algorithm, rJet);
fastjet::AreaDefinition areaDef(fastjet::active_area, fastjet::GhostedAreaSpec(1.0));
fastjet::ClusterSequenceArea cs(fjParticles, jetDef, areaDef);
std::vector<fastjet::PseudoJet> jets = fastjet::sorted_by_pt(cs.inclusive_jets());
auto [rhoPerp, rhoMPerp] = backgroundSub.estimateRhoPerpCone(fjParticles, jets);

// loop over reconstructed jets
bool isAtLeastOneJetSelected = false;
for (const auto& jet : jets) {

// jet must be fully contained in the acceptance
if ((std::fabs(jet.eta()) + rJet) > (maxEta - deltaEtaEdge))
continue;

// jet pt must be larger than threshold
auto jetForSub = jet;
fastjet::PseudoJet jetMinusBkg = backgroundSub.doRhoAreaSub(jetForSub, rhoPerp, rhoMPerp);
if (jetMinusBkg.pt() < minJetPt)
continue;
isAtLeastOneJetSelected = true;
}
if (isAtLeastOneJetSelected) {
registryMult.fill(HIST("multiplicityEvtsWithJet"), fjParticles.size());
}
}
PROCESS_SWITCH(AntinucleiInJets, processMultEvents, "Process Mult Events", true);

// Process QC
void processQC(SelectedCollisions::iterator const& collision, FullNucleiTracks const& tracks)
{
Expand Down Expand Up @@ -1087,13 +1164,13 @@ struct AntinucleiInJets {
}
PROCESS_SWITCH(AntinucleiInJets, processEfficiency, "process efficiency", false);

void processJetsMCgen(SimCollisions const& collisions, aod::McParticles const& mcParticles)
void processJetsMCgen(GenCollisions const& collisions, aod::McParticles const& mcParticles)
{
// Loop over all simulated collision events
for (const auto& collision : collisions) {

// Apply event selection: require sel8 and vertex position within the allowed z range
if (!collision.sel8() || std::fabs(collision.posZ()) > zVtx)
// Apply event selection: require vertex position within the allowed z range
if (std::fabs(collision.posZ()) > zVtx)
continue;

// Loop over all MC particles and select physical primaries within acceptance
Expand Down Expand Up @@ -1256,7 +1333,8 @@ struct AntinucleiInJets {
fastjet::PseudoJet jetMinusBkg = backgroundSub.doRhoAreaSub(jetForSub, rhoPerp, rhoMPerp);

// Apply jet pT threshold
if (getCorrectedPt(jetMinusBkg.pt(), responseMatrix) < minJetPt)
// if (getCorrectedPt(jetMinusBkg.pt(), responseMatrix) < minJetPt)
if (jetMinusBkg.pt() < minJetPt)
continue;

// Set up two perpendicular cone axes for underlying event estimation
Expand All @@ -1271,11 +1349,20 @@ struct AntinucleiInJets {
auto const& track = mcTracks.iteratorAt(particle.user_index());
if (!track.has_mcParticle())
continue;
const auto mcparticle = track.mcParticle();

// Fill DCA templates
if (mcparticle.pdgCode() == kProtonBar && passedTrackSelection(track) && std::fabs(track.dcaZ()) < maxDcaz) {
if (mcparticle.isPhysicalPrimary()) {
registryMC.fill(HIST("antiproton_prim_dca_jet"), track.pt(), track.dcaXY());
} else {
registryMC.fill(HIST("antiproton_all_dca_jet"), track.pt(), track.dcaXY());
}
}

// Apply standard track quality and PID selection
if (!passedTrackSelection(track) || std::fabs(track.dcaXY()) > maxDcaxy || std::fabs(track.dcaZ()) > maxDcaz)
continue;
const auto mcparticle = track.mcParticle();
if (track.sign() > 0 || mcparticle.pdgCode() != kProtonBar)
continue;

Expand Down Expand Up @@ -1310,15 +1397,26 @@ struct AntinucleiInJets {

// Analyze antiprotons in the Underlying Event (UE) using perpendicular cones
for (auto const& track : mcTracks) {
if (!passedTrackSelection(track) || std::fabs(track.dcaXY()) > maxDcaxy || std::fabs(track.dcaZ()) > maxDcaz)
continue;

if (!track.has_mcParticle())
continue;

const auto mcparticle = track.mcParticle();
if (track.sign() > 0 || mcparticle.pdgCode() != kProtonBar)
continue;

// Fill DCA templates
if (mcparticle.pdgCode() == kProtonBar && passedTrackSelection(track) && std::fabs(track.dcaZ()) < maxDcaz) {
if (mcparticle.isPhysicalPrimary()) {
registryMC.fill(HIST("antiproton_prim_dca_ue"), track.pt(), track.dcaXY());
} else {
registryMC.fill(HIST("antiproton_all_dca_ue"), track.pt(), track.dcaXY());
}
}

// Track Selection
if (!passedTrackSelection(track) || std::fabs(track.dcaXY()) > maxDcaxy || std::fabs(track.dcaZ()) > maxDcaz)
continue;

// Compute distance from UE cones
double deltaEtaUe1 = track.eta() - ueAxis1.Eta();
double deltaPhiUe1 = getDeltaPhi(track.phi(), ueAxis1.Phi());
Expand Down
Loading