Skip to content

Commit 4261565

Browse files
alcalivaalcaliva
authored
[PWGLF] Extend histograms with subsample dimension; random subsample … (#13927)
1 parent c585c36 commit 4261565

File tree

1 file changed

+75
-33
lines changed

1 file changed

+75
-33
lines changed

‎PWGLF/Tasks/Nuspex/antinucleiInJets.cxx‎

Lines changed: 75 additions & 33 deletions
Original file line numberDiff line numberDiff line change
@@ -51,7 +51,7 @@
5151
#include "TGrid.h"
5252
#include <TList.h>
5353
#include <TPDGCode.h>
54-
#include <TRandom.h>
54+
#include <TRandom3.h>
5555
#include <TVector2.h>
5656
#include <TVector3.h>
5757

@@ -64,6 +64,7 @@
6464
#include <fastjet/tools/JetMedianBackgroundEstimator.hh>
6565
#include <fastjet/tools/Subtractor.hh>
6666

67+
#include <chrono>
6768
#include <cmath>
6869
#include <memory>
6970
#include <random>
@@ -96,6 +97,9 @@ struct AntinucleiInJets {
9697
HistogramRegistry registryMult{"registryMult", {}, OutputObjHandlingPolicy::AnalysisObject, true, true};
9798
HistogramRegistry registryCorr{"registryCorr", {}, OutputObjHandlingPolicy::AnalysisObject, true, true};
9899

100+
// Random generator for subsample assignment
101+
TRandom3 mRand;
102+
99103
// Event selection criteria
100104
Configurable<bool> rejectITSROFBorder{"rejectITSROFBorder", true, "Reject events near the ITS ROF border"};
101105
Configurable<bool> rejectTFBorder{"rejectTFBorder", true, "Reject events near the TF border"};
@@ -205,6 +209,10 @@ struct AntinucleiInJets {
205209
itsResponse.setMCDefaultParameters();
206210
}
207211

212+
// Initialize random seed using high-resolution clock to ensure unique sequences across parallel Grid jobs
213+
auto time_seed = std::chrono::high_resolution_clock::now().time_since_epoch().count();
214+
mRand.SetSeed(time_seed);
215+
208216
// Load reweighting histograms from CCDB if antinuclei efficiency processing is enabled
209217
if (doprocessAntinucleiEfficiency || doprocessJetsMCgen || doprocessJetsMCrec) {
210218
ccdb->setURL(urlToCcdb.value);
@@ -252,6 +260,9 @@ struct AntinucleiInJets {
252260
// Event counters
253261
registryData.add("number_of_events_data", "number of events in data", HistType::kTH1F, {{20, 0, 20, "counter"}});
254262

263+
// Configuration
264+
registryData.add("settingData", "settingData", HistType::kTH2F, {{100, 0.0, 50.0, "min #it{p}^{jet}_{T} (GeV/#it{c})"}, {20, 0.0, 1.0, "#it{R}_{jet}"}});
265+
255266
// Jet effective area over piR^2
256267
registryData.add("jetEffectiveAreaOverPiR2", "jet effective area / piR^2", HistType::kTH1F, {{2000, 0, 2, "Area/#piR^{2}"}});
257268

@@ -308,6 +319,7 @@ struct AntinucleiInJets {
308319
// Generated spectra of antiprotons
309320
registryMC.add("antiproton_gen_jet", "antiproton_gen_jet", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}});
310321
registryMC.add("antiproton_gen_ue", "antiproton_gen_ue", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}});
322+
registryMC.add("antiproton_gen_full", "antiproton_gen_full", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}});
311323

312324
// Normalization histogram
313325
registryMC.add("antiproton_deltay_deltaphi_jet", "antiproton_deltay_deltaphi_jet", HistType::kTH2F, {{2000, -1.0, 1.0, "#Delta#it{y}"}, {2000, 0.0, 2.0, "#Delta#phi"}});
@@ -325,11 +337,16 @@ struct AntinucleiInJets {
325337
registryMC.add("recEvents", "number of reconstructed events in mc", HistType::kTH1F, {{20, 0, 20, "counter"}});
326338
registryMC.add("recJets", "number of reconstructed jets", HistType::kTH1F, {{10, 0, 10, "counter"}});
327339

340+
// Configuration
341+
registryMC.add("settingMC", "settingMC", HistType::kTH2F, {{100, 0.0, 50.0, "min #it{p}^{jet}_{T} (GeV/#it{c})"}, {20, 0.0, 1.0, "#it{R}_{jet}"}});
342+
328343
// Reconstructed spectra of antiprotons
329344
registryMC.add("antiproton_rec_tpc_jet", "antiproton_rec_tpc_jet", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}});
330345
registryMC.add("antiproton_rec_tof_jet", "antiproton_rec_tof_jet", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}});
331346
registryMC.add("antiproton_rec_tpc_ue", "antiproton_rec_tpc_ue", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}});
332347
registryMC.add("antiproton_rec_tof_ue", "antiproton_rec_tof_ue", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}});
348+
registryMC.add("antiproton_rec_tpc_full", "antiproton_rec_tpc_full", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}});
349+
registryMC.add("antiproton_rec_tof_full", "antiproton_rec_tof_full", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}});
333350

334351
// Fraction of primary antiprotons
335352
registryMC.add("antiproton_prim_jet", "antiproton_prim_jet", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}});
@@ -455,43 +472,44 @@ struct AntinucleiInJets {
455472
const AxisSpec nBarD2Axis{100, 0.0, 100.0, "N_{#bar{d}}^{i} #times N_{#bar{d}}^{j}"};
456473
const AxisSpec nBarP2Axis{100, 0.0, 100.0, "N_{#bar{p}}^{i} #times N_{#bar{p}}^{j}"};
457474
const AxisSpec nBarDnBarPAxis{100, 0.0, 100.0, "N_{#bar{d}}^{i} #times N_{#bar{p}}^{j}"};
475+
const AxisSpec subsampleAxis{20, 0, 20, "Subsample Index"};
458476

459477
// Event counter
460478
registryCorr.add("eventCounter", "number of events", HistType::kTH1F, {{20, 0, 20, "counter"}});
461-
registryCorr.add("eventCounter_centrality_fullEvent", "Number of events per centrality (Full Event)", HistType::kTH1F, {multiplicityAxis});
462-
registryCorr.add("eventCounter_centrality_jet", "Number of events per centrality (Jet)", HistType::kTH1F, {multiplicityAxis});
463-
registryCorr.add("eventCounter_centrality_ue", "Number of events per centrality (Underlying Event)", HistType::kTH1F, {multiplicityAxis});
479+
registryCorr.add("eventCounter_centrality_fullEvent", "Number of events per centrality (Full Event)", HistType::kTH2F, {multiplicityAxis, subsampleAxis});
480+
// registryCorr.add("eventCounter_centrality_jet", "Number of events per centrality (Jet)", HistType::kTH1F, {multiplicityAxis});
481+
// registryCorr.add("eventCounter_centrality_ue", "Number of events per centrality (Underlying Event)", HistType::kTH1F, {multiplicityAxis});
464482

465483
// Correlation histograms: antiproton vs. antideuteron number vs. event multiplicity
466-
registryCorr.add("rho_jet", "rho_jet", HistType::kTH3F, {nAntideuteronsAxis, nAntiprotonsAxis, multiplicityAxis});
467-
registryCorr.add("rho_ue", "rho_ue", HistType::kTH3F, {nAntideuteronsAxis, nAntiprotonsAxis, multiplicityAxis});
468-
registryCorr.add("rho_fullEvent", "rho_fullEvent", HistType::kTH3F, {nAntideuteronsAxis, nAntiprotonsAxis, multiplicityAxis});
484+
// registryCorr.add("rho_jet", "rho_jet", HistType::kTH3F, {nAntideuteronsAxis, nAntiprotonsAxis, multiplicityAxis});
485+
// registryCorr.add("rho_ue", "rho_ue", HistType::kTH3F, {nAntideuteronsAxis, nAntiprotonsAxis, multiplicityAxis});
486+
registryCorr.add("rho_fullEvent", "rho_fullEvent", HistType::kTHnSparseD, {nAntideuteronsAxis, nAntiprotonsAxis, multiplicityAxis, subsampleAxis});
469487

470488
// Correlation histograms: net antiproton vs. net antideuteron numbers
471-
registryCorr.add("rho_netP_netD_jet", "rho_netP_netD_jet", HistType::kTH2F, {nAntideuteronsAxis, nAntiprotonsAxis});
472-
registryCorr.add("rho_netP_netD_ue", "rho_netP_netD_ue", HistType::kTH2F, {nAntideuteronsAxis, nAntiprotonsAxis});
473-
registryCorr.add("rho_netP_netD_fullEvent", "rho_netP_netD_fullEvent", HistType::kTH2F, {nAntideuteronsAxis, nAntiprotonsAxis});
489+
// registryCorr.add("rho_netP_netD_jet", "rho_netP_netD_jet", HistType::kTH2F, {nAntideuteronsAxis, nAntiprotonsAxis});
490+
// registryCorr.add("rho_netP_netD_ue", "rho_netP_netD_ue", HistType::kTH2F, {nAntideuteronsAxis, nAntiprotonsAxis});
491+
registryCorr.add("rho_netP_netD_fullEvent", "rho_netP_netD_fullEvent", HistType::kTH3F, {nAntideuteronsAxis, nAntiprotonsAxis, subsampleAxis});
474492

475493
// Efficiency histograms jet
476-
registryCorr.add("q1d_jet", "q1d_jet", HistType::kTH3F, {nAntideuteronsAxis, ptPerNucleonAxis, multiplicityAxis});
477-
registryCorr.add("q1p_jet", "q1p_jet", HistType::kTH3F, {nAntiprotonsAxis, ptPerNucleonAxis, multiplicityAxis});
478-
registryCorr.add("q1d_square_jet", "q1d_square_jet", HistType::kTHnSparseD, {ptPerNucleonAxis, ptPerNucleonAxis, nBarD2Axis, multiplicityAxis});
479-
registryCorr.add("q1p_square_jet", "q1p_square_jet", HistType::kTHnSparseD, {ptPerNucleonAxis, ptPerNucleonAxis, nBarP2Axis, multiplicityAxis});
480-
registryCorr.add("q1d_q1p_jet", "q1d_q1p_jet", HistType::kTHnSparseD, {ptPerNucleonAxis, ptPerNucleonAxis, nBarDnBarPAxis, multiplicityAxis});
494+
// registryCorr.add("q1d_jet", "q1d_jet", HistType::kTH3F, {nAntideuteronsAxis, ptPerNucleonAxis, multiplicityAxis});
495+
// registryCorr.add("q1p_jet", "q1p_jet", HistType::kTH3F, {nAntiprotonsAxis, ptPerNucleonAxis, multiplicityAxis});
496+
// registryCorr.add("q1d_square_jet", "q1d_square_jet", HistType::kTHnSparseD, {ptPerNucleonAxis, ptPerNucleonAxis, nBarD2Axis, multiplicityAxis});
497+
// registryCorr.add("q1p_square_jet", "q1p_square_jet", HistType::kTHnSparseD, {ptPerNucleonAxis, ptPerNucleonAxis, nBarP2Axis, multiplicityAxis});
498+
// registryCorr.add("q1d_q1p_jet", "q1d_q1p_jet", HistType::kTHnSparseD, {ptPerNucleonAxis, ptPerNucleonAxis, nBarDnBarPAxis, multiplicityAxis});
481499

482500
// Efficiency histograms UE
483-
registryCorr.add("q1d_ue", "q1d_ue", HistType::kTH3F, {nAntideuteronsAxis, ptPerNucleonAxis, multiplicityAxis});
484-
registryCorr.add("q1p_ue", "q1p_ue", HistType::kTH3F, {nAntiprotonsAxis, ptPerNucleonAxis, multiplicityAxis});
485-
registryCorr.add("q1d_square_ue", "q1d_square_ue", HistType::kTHnSparseD, {ptPerNucleonAxis, ptPerNucleonAxis, nBarD2Axis, multiplicityAxis});
486-
registryCorr.add("q1p_square_ue", "q1p_square_ue", HistType::kTHnSparseD, {ptPerNucleonAxis, ptPerNucleonAxis, nBarP2Axis, multiplicityAxis});
487-
registryCorr.add("q1d_q1p_ue", "q1d_q1p_ue", HistType::kTHnSparseD, {ptPerNucleonAxis, ptPerNucleonAxis, nBarDnBarPAxis, multiplicityAxis});
501+
// registryCorr.add("q1d_ue", "q1d_ue", HistType::kTH3F, {nAntideuteronsAxis, ptPerNucleonAxis, multiplicityAxis});
502+
// registryCorr.add("q1p_ue", "q1p_ue", HistType::kTH3F, {nAntiprotonsAxis, ptPerNucleonAxis, multiplicityAxis});
503+
// registryCorr.add("q1d_square_ue", "q1d_square_ue", HistType::kTHnSparseD, {ptPerNucleonAxis, ptPerNucleonAxis, nBarD2Axis, multiplicityAxis});
504+
// registryCorr.add("q1p_square_ue", "q1p_square_ue", HistType::kTHnSparseD, {ptPerNucleonAxis, ptPerNucleonAxis, nBarP2Axis, multiplicityAxis});
505+
// registryCorr.add("q1d_q1p_ue", "q1d_q1p_ue", HistType::kTHnSparseD, {ptPerNucleonAxis, ptPerNucleonAxis, nBarDnBarPAxis, multiplicityAxis});
488506

489507
// Efficiency histograms full event
490-
registryCorr.add("q1d_fullEvent", "q1d_fullEvent", HistType::kTH3F, {nAntideuteronsAxis, ptPerNucleonAxis, multiplicityAxis});
491-
registryCorr.add("q1p_fullEvent", "q1p_fullEvent", HistType::kTH3F, {nAntiprotonsAxis, ptPerNucleonAxis, multiplicityAxis});
492-
registryCorr.add("q1d_square_fullEvent", "q1d_square_fullEvent", HistType::kTHnSparseD, {ptPerNucleonAxis, ptPerNucleonAxis, nBarD2Axis, multiplicityAxis});
493-
registryCorr.add("q1p_square_fullEvent", "q1p_square_fullEvent", HistType::kTHnSparseD, {ptPerNucleonAxis, ptPerNucleonAxis, nBarP2Axis, multiplicityAxis});
494-
registryCorr.add("q1d_q1p_fullEvent", "q1d_q1p_fullEvent", HistType::kTHnSparseD, {ptPerNucleonAxis, ptPerNucleonAxis, nBarDnBarPAxis, multiplicityAxis});
508+
registryCorr.add("q1d_fullEvent", "q1d_fullEvent", HistType::kTHnSparseD, {nAntideuteronsAxis, ptPerNucleonAxis, multiplicityAxis, subsampleAxis});
509+
registryCorr.add("q1p_fullEvent", "q1p_fullEvent", HistType::kTHnSparseD, {nAntiprotonsAxis, ptPerNucleonAxis, multiplicityAxis, subsampleAxis});
510+
registryCorr.add("q1d_square_fullEvent", "q1d_square_fullEvent", HistType::kTHnSparseD, {ptPerNucleonAxis, ptPerNucleonAxis, nBarD2Axis, multiplicityAxis, subsampleAxis});
511+
registryCorr.add("q1p_square_fullEvent", "q1p_square_fullEvent", HistType::kTHnSparseD, {ptPerNucleonAxis, ptPerNucleonAxis, nBarP2Axis, multiplicityAxis, subsampleAxis});
512+
registryCorr.add("q1d_q1p_fullEvent", "q1d_q1p_fullEvent", HistType::kTHnSparseD, {ptPerNucleonAxis, ptPerNucleonAxis, nBarDnBarPAxis, multiplicityAxis, subsampleAxis});
495513
}
496514
}
497515

@@ -924,6 +942,7 @@ struct AntinucleiInJets {
924942
{
925943
// Event counter: before event selection
926944
registryData.fill(HIST("number_of_events_data"), 0.5);
945+
registryData.fill(HIST("settingData"), minJetPt.value, rJet.value);
927946

928947
// Retrieve the bunch crossing information with timestamps from the collision
929948
auto bc = collision.template bc_as<aod::BCsWithTimestamps>();
@@ -1801,6 +1820,7 @@ struct AntinucleiInJets {
18011820
if (particle.pdgCode() == PDG_t::kProtonBar) {
18021821
TVector3 pVec(particle.px(), particle.py(), particle.pz());
18031822
protonMomentum.emplace_back(pVec);
1823+
registryMC.fill(HIST("antiproton_gen_full"), particle.pt());
18041824
}
18051825

18061826
// 4-momentum representation of a particle
@@ -1950,6 +1970,9 @@ struct AntinucleiInJets {
19501970
// Loop over all reconstructed collisions
19511971
for (const auto& collision : collisions) {
19521972

1973+
// Configuration
1974+
registryMC.fill(HIST("settingMC"), minJetPt.value, rJet.value);
1975+
19531976
// Increment event counter
19541977
eventCounter++;
19551978

@@ -2013,6 +2036,19 @@ struct AntinucleiInJets {
20132036
// Store track index for antiproton tracks
20142037
if (passedTrackSelection(track) && track.sign() < 0 && mcparticle.pdgCode() == PDG_t::kProtonBar) {
20152038
antiprotonTrackIndex.emplace_back(id);
2039+
2040+
double nsigmaTPCPr = track.tpcNSigmaPr();
2041+
double nsigmaTOFPr = track.tofNSigmaPr();
2042+
double pt = track.pt();
2043+
double dcaxy = track.dcaXY();
2044+
double dcaz = track.dcaZ();
2045+
2046+
if (mcparticle.isPhysicalPrimary() && std::fabs(dcaxy) < maxDcaxy && std::fabs(dcaz) < maxDcaz && nsigmaTPCPr > minNsigmaTpc && nsigmaTPCPr < maxNsigmaTpc) {
2047+
registryMC.fill(HIST("antiproton_rec_tpc_full"), pt);
2048+
if (track.hasTOF() && nsigmaTOFPr > minNsigmaTof && nsigmaTOFPr < maxNsigmaTof) {
2049+
registryMC.fill(HIST("antiproton_rec_tof_full"), pt);
2050+
}
2051+
}
20162052
}
20172053

20182054
// Apply track selection for jet reconstruction
@@ -2649,11 +2685,14 @@ struct AntinucleiInJets {
26492685
return;
26502686
registryCorr.fill(HIST("eventCounter"), 7.5);
26512687

2688+
// Assign event to a random subsample (0-19)
2689+
double sampleId = mRand.Integer(20) + 0.5;
2690+
26522691
// Multiplicity percentile
26532692
const float multiplicity = collision.centFT0M();
26542693

26552694
// Fill event counter vs centrality (full Event region)
2656-
registryCorr.fill(HIST("eventCounter_centrality_fullEvent"), multiplicity);
2695+
registryCorr.fill(HIST("eventCounter_centrality_fullEvent"), multiplicity, sampleId);
26572696

26582697
// pt/A bins
26592698
std::vector<double> ptOverAbins = {0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0};
@@ -2732,23 +2771,25 @@ struct AntinucleiInJets {
27322771
// Fill correlation histograms
27332772
int netProtonFullEvent = nTotProtonFullEvent - nTotAntiprotonFullEvent;
27342773
int netDeuteronFullEvent = nTotDeuteronFullEvent - nTotAntideuteronFullEvent;
2735-
registryCorr.fill(HIST("rho_fullEvent"), nTotAntideuteronFullEvent, nTotAntiprotonFullEvent, multiplicity);
2736-
registryCorr.fill(HIST("rho_netP_netD_fullEvent"), netDeuteronFullEvent, netProtonFullEvent);
2774+
2775+
registryCorr.fill(HIST("rho_fullEvent"), nTotAntideuteronFullEvent, nTotAntiprotonFullEvent, multiplicity, sampleId);
2776+
registryCorr.fill(HIST("rho_netP_netD_fullEvent"), netDeuteronFullEvent, netProtonFullEvent, sampleId);
27372777

27382778
// Fill efficiency histograms
27392779
for (int i = 0; i < nBins; i++) {
27402780
double ptAcenteri = 0.5 * (ptOverAbins[i] + ptOverAbins[i + 1]);
27412781

2742-
registryCorr.fill(HIST("q1d_fullEvent"), nAntideuteronFullEvent[i], ptAcenteri, multiplicity);
2743-
registryCorr.fill(HIST("q1p_fullEvent"), nAntiprotonFullEvent[i], ptAcenteri, multiplicity);
2782+
registryCorr.fill(HIST("q1d_fullEvent"), nAntideuteronFullEvent[i], ptAcenteri, multiplicity, sampleId);
2783+
registryCorr.fill(HIST("q1p_fullEvent"), nAntiprotonFullEvent[i], ptAcenteri, multiplicity, sampleId);
27442784
for (int j = 0; j < nBins; j++) {
27452785
double ptAcenterj = 0.5 * (ptOverAbins[j] + ptOverAbins[j + 1]);
2746-
registryCorr.fill(HIST("q1d_square_fullEvent"), ptAcenteri, ptAcenterj, nAntideuteronFullEvent[i] * nAntideuteronFullEvent[j], multiplicity);
2747-
registryCorr.fill(HIST("q1p_square_fullEvent"), ptAcenteri, ptAcenterj, nAntiprotonFullEvent[i] * nAntiprotonFullEvent[j], multiplicity);
2748-
registryCorr.fill(HIST("q1d_q1p_fullEvent"), ptAcenteri, ptAcenterj, nAntideuteronFullEvent[i] * nAntiprotonFullEvent[j], multiplicity);
2786+
registryCorr.fill(HIST("q1d_square_fullEvent"), ptAcenteri, ptAcenterj, (nAntideuteronFullEvent[i] * nAntideuteronFullEvent[j]), multiplicity, sampleId);
2787+
registryCorr.fill(HIST("q1p_square_fullEvent"), ptAcenteri, ptAcenterj, (nAntiprotonFullEvent[i] * nAntiprotonFullEvent[j]), multiplicity, sampleId);
2788+
registryCorr.fill(HIST("q1d_q1p_fullEvent"), ptAcenteri, ptAcenterj, (nAntideuteronFullEvent[i] * nAntiprotonFullEvent[j]), multiplicity, sampleId);
27492789
}
27502790
}
27512791

2792+
/*
27522793
// Loop over reconstructed tracks (refactoring: this part can be incorporated above)
27532794
int id(-1);
27542795
std::vector<fastjet::PseudoJet> fjParticles;
@@ -3013,6 +3054,7 @@ struct AntinucleiInJets {
30133054
registryCorr.fill(HIST("eventCounter"), 9.5);
30143055
registryCorr.fill(HIST("eventCounter_centrality_jet"), multiplicity);
30153056
}
3057+
*/
30163058
}
30173059
PROCESS_SWITCH(AntinucleiInJets, processCorr, "Process Correlation analysis", false);
30183060
};

0 commit comments

Comments
 (0)