Skip to content

Commit d7e0aab

Browse files
alcalivaalcaliva
authored
[PWGLF] added DCA templates and function to compute multiplicity in different event classes (#11768)
1 parent 5691e46 commit d7e0aab

File tree

1 file changed

+109
-11
lines changed

1 file changed

+109
-11
lines changed

‎PWGLF/Tasks/Nuspex/antinucleiInJets.cxx‎

Lines changed: 109 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -77,6 +77,7 @@ using std::array;
7777

7878
using SelectedCollisions = soa::Join<aod::Collisions, aod::EvSels>;
7979
using SimCollisions = soa::Join<aod::Collisions, aod::EvSels, aod::McCollisionLabels>;
80+
using GenCollisions = aod::McCollisions;
8081

8182
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>;
8283

@@ -88,9 +89,11 @@ struct AntinucleiInJets {
8889
HistogramRegistry registryData{"registryData", {}, OutputObjHandlingPolicy::AnalysisObject, true, true};
8990
HistogramRegistry registryMC{"registryMC", {}, OutputObjHandlingPolicy::AnalysisObject, true, true};
9091
HistogramRegistry registryQC{"registryQC", {}, OutputObjHandlingPolicy::AnalysisObject, true, true};
92+
HistogramRegistry registryMult{"registryMult", {}, OutputObjHandlingPolicy::AnalysisObject, true, true};
9193

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

185+
if (doprocessMultEvents) {
186+
registryMult.add("multiplicityEvtsPtLeading", "multiplicityEvtsPtLeading", HistType::kTH1F, {{1000, 0, 1000, "#it{N}_{ch}"}});
187+
registryMult.add("multiplicityEvtsWithJet", "multiplicityEvtsWithJet", HistType::kTH1F, {{1000, 0, 1000, "#it{N}_{ch}"}});
188+
}
189+
182190
// data histograms
183191
if (doprocessData) {
184192

@@ -191,8 +199,8 @@ struct AntinucleiInJets {
191199
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}"}});
192200
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}"}});
193201
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}"}});
194-
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)"}});
195-
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)"}});
202+
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)"}});
203+
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)"}});
196204

197205
// antideuterons
198206
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}"}});
@@ -271,6 +279,12 @@ struct AntinucleiInJets {
271279
registryMC.add("antiproton_prim_ue", "antiproton_prim_ue", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}});
272280
registryMC.add("antiproton_incl_ue", "antiproton_incl_ue", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}});
273281

282+
// DCA Templates
283+
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)"}});
284+
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)"}});
285+
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)"}});
286+
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)"}});
287+
274288
// antiproton reweighting
275289
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}"}});
276290
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}"}});
@@ -597,7 +611,8 @@ struct AntinucleiInJets {
597611
// jet pt must be larger than threshold
598612
auto jetForSub = jet;
599613
fastjet::PseudoJet jetMinusBkg = backgroundSub.doRhoAreaSub(jetForSub, rhoPerp, rhoMPerp);
600-
if (getCorrectedPt(jetMinusBkg.pt(), responseMatrix) < minJetPt)
614+
// if (getCorrectedPt(jetMinusBkg.pt(), responseMatrix) < minJetPt)
615+
if (jetMinusBkg.pt() < minJetPt)
601616
continue;
602617
isAtLeastOneJetSelected = true;
603618

@@ -768,6 +783,68 @@ struct AntinucleiInJets {
768783
}
769784
PROCESS_SWITCH(AntinucleiInJets, processData, "Process Data", true);
770785

786+
void processMultEvents(SelectedCollisions::iterator const& collision, FullNucleiTracks const& tracks)
787+
{
788+
// event selection
789+
if (!collision.sel8() || std::fabs(collision.posZ()) > zVtx)
790+
return;
791+
792+
// Leading Track
793+
double ptMax(0.0);
794+
795+
// loop over reconstructed tracks
796+
int id(-1);
797+
std::vector<fastjet::PseudoJet> fjParticles;
798+
for (auto const& track : tracks) {
799+
id++;
800+
if (!passedTrackSelectionForJetReconstruction(track))
801+
continue;
802+
if (track.pt() > ptMax) {
803+
ptMax = track.pt();
804+
}
805+
806+
// 4-momentum representation of a particle
807+
fastjet::PseudoJet fourMomentum(track.px(), track.py(), track.pz(), track.energy(MassPionCharged));
808+
fourMomentum.set_user_index(id);
809+
fjParticles.emplace_back(fourMomentum);
810+
}
811+
// reject empty events
812+
if (fjParticles.empty()) {
813+
return;
814+
}
815+
816+
if (ptMax > ptLeadingMin) {
817+
registryMult.fill(HIST("multiplicityEvtsPtLeading"), fjParticles.size());
818+
}
819+
820+
// cluster particles using the anti-kt algorithm
821+
fastjet::JetDefinition jetDef(fastjet::antikt_algorithm, rJet);
822+
fastjet::AreaDefinition areaDef(fastjet::active_area, fastjet::GhostedAreaSpec(1.0));
823+
fastjet::ClusterSequenceArea cs(fjParticles, jetDef, areaDef);
824+
std::vector<fastjet::PseudoJet> jets = fastjet::sorted_by_pt(cs.inclusive_jets());
825+
auto [rhoPerp, rhoMPerp] = backgroundSub.estimateRhoPerpCone(fjParticles, jets);
826+
827+
// loop over reconstructed jets
828+
bool isAtLeastOneJetSelected = false;
829+
for (const auto& jet : jets) {
830+
831+
// jet must be fully contained in the acceptance
832+
if ((std::fabs(jet.eta()) + rJet) > (maxEta - deltaEtaEdge))
833+
continue;
834+
835+
// jet pt must be larger than threshold
836+
auto jetForSub = jet;
837+
fastjet::PseudoJet jetMinusBkg = backgroundSub.doRhoAreaSub(jetForSub, rhoPerp, rhoMPerp);
838+
if (jetMinusBkg.pt() < minJetPt)
839+
continue;
840+
isAtLeastOneJetSelected = true;
841+
}
842+
if (isAtLeastOneJetSelected) {
843+
registryMult.fill(HIST("multiplicityEvtsWithJet"), fjParticles.size());
844+
}
845+
}
846+
PROCESS_SWITCH(AntinucleiInJets, processMultEvents, "Process Mult Events", true);
847+
771848
// Process QC
772849
void processQC(SelectedCollisions::iterator const& collision, FullNucleiTracks const& tracks)
773850
{
@@ -1087,13 +1164,13 @@ struct AntinucleiInJets {
10871164
}
10881165
PROCESS_SWITCH(AntinucleiInJets, processEfficiency, "process efficiency", false);
10891166

1090-
void processJetsMCgen(SimCollisions const& collisions, aod::McParticles const& mcParticles)
1167+
void processJetsMCgen(GenCollisions const& collisions, aod::McParticles const& mcParticles)
10911168
{
10921169
// Loop over all simulated collision events
10931170
for (const auto& collision : collisions) {
10941171

1095-
// Apply event selection: require sel8 and vertex position within the allowed z range
1096-
if (!collision.sel8() || std::fabs(collision.posZ()) > zVtx)
1172+
// Apply event selection: require vertex position within the allowed z range
1173+
if (std::fabs(collision.posZ()) > zVtx)
10971174
continue;
10981175

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

12581335
// Apply jet pT threshold
1259-
if (getCorrectedPt(jetMinusBkg.pt(), responseMatrix) < minJetPt)
1336+
// if (getCorrectedPt(jetMinusBkg.pt(), responseMatrix) < minJetPt)
1337+
if (jetMinusBkg.pt() < minJetPt)
12601338
continue;
12611339

12621340
// Set up two perpendicular cone axes for underlying event estimation
@@ -1271,11 +1349,20 @@ struct AntinucleiInJets {
12711349
auto const& track = mcTracks.iteratorAt(particle.user_index());
12721350
if (!track.has_mcParticle())
12731351
continue;
1352+
const auto mcparticle = track.mcParticle();
1353+
1354+
// Fill DCA templates
1355+
if (mcparticle.pdgCode() == kProtonBar && passedTrackSelection(track) && std::fabs(track.dcaZ()) < maxDcaz) {
1356+
if (mcparticle.isPhysicalPrimary()) {
1357+
registryMC.fill(HIST("antiproton_prim_dca_jet"), track.pt(), track.dcaXY());
1358+
} else {
1359+
registryMC.fill(HIST("antiproton_all_dca_jet"), track.pt(), track.dcaXY());
1360+
}
1361+
}
12741362

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

@@ -1310,15 +1397,26 @@ struct AntinucleiInJets {
13101397

13111398
// Analyze antiprotons in the Underlying Event (UE) using perpendicular cones
13121399
for (auto const& track : mcTracks) {
1313-
if (!passedTrackSelection(track) || std::fabs(track.dcaXY()) > maxDcaxy || std::fabs(track.dcaZ()) > maxDcaz)
1314-
continue;
1400+
13151401
if (!track.has_mcParticle())
13161402
continue;
1317-
13181403
const auto mcparticle = track.mcParticle();
13191404
if (track.sign() > 0 || mcparticle.pdgCode() != kProtonBar)
13201405
continue;
13211406

1407+
// Fill DCA templates
1408+
if (mcparticle.pdgCode() == kProtonBar && passedTrackSelection(track) && std::fabs(track.dcaZ()) < maxDcaz) {
1409+
if (mcparticle.isPhysicalPrimary()) {
1410+
registryMC.fill(HIST("antiproton_prim_dca_ue"), track.pt(), track.dcaXY());
1411+
} else {
1412+
registryMC.fill(HIST("antiproton_all_dca_ue"), track.pt(), track.dcaXY());
1413+
}
1414+
}
1415+
1416+
// Track Selection
1417+
if (!passedTrackSelection(track) || std::fabs(track.dcaXY()) > maxDcaxy || std::fabs(track.dcaZ()) > maxDcaz)
1418+
continue;
1419+
13221420
// Compute distance from UE cones
13231421
double deltaEtaUe1 = track.eta() - ueAxis1.Eta();
13241422
double deltaPhiUe1 = getDeltaPhi(track.phi(), ueAxis1.Phi());

0 commit comments

Comments
 (0)