Skip to content

Commit 31c5657

Browse files
FrancescaCasilloFrancesca Casillo
FrancescaCasillo
and
Francesca Casillo
authored
[PWGLF] Added fullEvent histograms to processCoalescence and a new function (processSystCorr) for systematic studies (#14379)
Co-authored-by: Francesca Casillo <francesca@MacBook-Air-di-Francesca-3.local>
1 parent 9439541 commit 31c5657

File tree

1 file changed

+281
-0
lines changed

1 file changed

+281
-0
lines changed

‎PWGLF/Tasks/Nuspex/antinucleiInJets.cxx‎

Lines changed: 281 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -476,8 +476,10 @@ struct AntinucleiInJets {
476476
// Coalescence
477477
if (doprocessCoalescence) {
478478
registryMC.add("genEventsCoalescence", "genEventsCoalescence", HistType::kTH1F, {{20, 0, 20, "counter"}});
479+
registryMC.add("antideuteron_coal_fullEvent", "antideuteron_coal_fullEvent", HistType::kTH1F, {{nbins, 2 * min, 2 * max, "#it{p}_{T} (GeV/#it{c})"}});
479480
registryMC.add("antideuteron_coal_jet", "antideuteron_coal_jet", HistType::kTH1F, {{nbins, 2 * min, 2 * max, "#it{p}_{T} (GeV/#it{c})"}});
480481
registryMC.add("antideuteron_coal_ue", "antideuteron_coal_ue", HistType::kTH1F, {{nbins, 2 * min, 2 * max, "#it{p}_{T} (GeV/#it{c})"}});
482+
registryMC.add("antiproton_coal_fullEvent", "antiproton_coal_fullEvent", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}});
481483
registryMC.add("antiproton_coal_jet", "antiproton_coal_jet", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}});
482484
registryMC.add("antiproton_coal_ue", "antiproton_coal_ue", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}});
483485
}
@@ -563,6 +565,34 @@ struct AntinucleiInJets {
563565
registryCorr.add("q1p_square_fullEvent", "q1p_square_fullEvent", HistType::kTHnSparseD, {ptPerNucleonAxis, ptPerNucleonAxis, nBarP2Axis, multiplicityAxis, subsampleAxis});
564566
registryCorr.add("q1d_q1p_fullEvent", "q1d_q1p_fullEvent", HistType::kTHnSparseD, {ptPerNucleonAxis, ptPerNucleonAxis, nBarDnBarPAxis, multiplicityAxis, subsampleAxis});
565567
}
568+
569+
// Systematic uncertainties on correlation analysis
570+
if (doprocessSystCorr) {
571+
572+
// Axes definitions for multidimensional histogram binning
573+
const AxisSpec multiplicityAxis{100, 0.0, 100.0, "multiplicity percentile"};
574+
const AxisSpec ptPerNucleonAxis{5, 0.4, 0.9, "{p}_{T}/A (GeV/#it{c})"};
575+
const AxisSpec nAntideuteronsAxis{10, 0.0, 10.0, "N_{#bar{d}}"};
576+
const AxisSpec nAntiprotonsAxis{10, 0.0, 10.0, "N_{#bar{p}}"};
577+
const AxisSpec nBarD2Axis{100, 0.0, 100.0, "N_{#bar{d}}^{i} #times N_{#bar{d}}^{j}"};
578+
const AxisSpec nBarP2Axis{100, 0.0, 100.0, "N_{#bar{p}}^{i} #times N_{#bar{p}}^{j}"};
579+
const AxisSpec nBarDnBarPAxis{100, 0.0, 100.0, "N_{#bar{d}}^{i} #times N_{#bar{p}}^{j}"};
580+
const AxisSpec systAxis{nSyst, 0, static_cast<double>(nSyst), "Systematic Variation Index"};
581+
582+
registryCorr.add("eventCounter_syst", "number of events syst", HistType::kTH1F, {{20, 0, 20, "counter"}});
583+
registryCorr.add("eventCounter_centrality_fullEvent_syst", "Number of events per centrality (Full Event) Syst", HistType::kTH2F, {multiplicityAxis, systAxis});
584+
585+
// Correlation histograms
586+
registryCorr.add("rho_fullEvent_syst", "rho_fullEvent_syst", HistType::kTHnSparseD, {nAntideuteronsAxis, nAntiprotonsAxis, multiplicityAxis, systAxis});
587+
registryCorr.add("rho_netP_netD_fullEvent_syst", "rho_netP_netD_fullEvent_syst", HistType::kTH3F, {nAntideuteronsAxis, nAntiprotonsAxis, systAxis});
588+
589+
// Efficiency histograms full event
590+
registryCorr.add("q1d_fullEvent_syst", "q1d_fullEvent_syst", HistType::kTHnSparseD, {nAntideuteronsAxis, ptPerNucleonAxis, multiplicityAxis, systAxis});
591+
registryCorr.add("q1p_fullEvent_syst", "q1p_fullEvent_syst", HistType::kTHnSparseD, {nAntiprotonsAxis, ptPerNucleonAxis, multiplicityAxis, systAxis});
592+
registryCorr.add("q1d_square_fullEvent_syst", "q1d_square_fullEvent_syst", HistType::kTHnSparseD, {ptPerNucleonAxis, ptPerNucleonAxis, nBarD2Axis, multiplicityAxis, systAxis});
593+
registryCorr.add("q1p_square_fullEvent_syst", "q1p_square_fullEvent_syst", HistType::kTHnSparseD, {ptPerNucleonAxis, ptPerNucleonAxis, nBarP2Axis, multiplicityAxis, systAxis});
594+
registryCorr.add("q1d_q1p_fullEvent_syst", "q1d_q1p_fullEvent_syst", HistType::kTHnSparseD, {ptPerNucleonAxis, ptPerNucleonAxis, nBarDnBarPAxis, multiplicityAxis, systAxis});
595+
}
566596
}
567597

568598
void getReweightingHistograms(o2::framework::Service<o2::ccdb::BasicCCDBManager> const& ccdbObj, TString filepath, TString antip, TString antilambda, TString antisigma, TString antixi, TString antiomega, TString jet, TString ue)
@@ -983,6 +1013,54 @@ struct AntinucleiInJets {
9831013
return (track.hasTOF() && std::abs(nsigmaTOF) < NsigmaMax);
9841014
}
9851015

1016+
// Selection of (anti)protons with systematic variations
1017+
template <typename ProtonTrack>
1018+
bool isProtonSyst(const ProtonTrack& track, double minSigTPC, double maxSigTPC, double minSigTOF, double maxSigTOF)
1019+
{
1020+
// Constant
1021+
static constexpr double PtThreshold = 0.6;
1022+
1023+
// PID variables and transverse momentum of the track
1024+
const double nsigmaTPC = track.tpcNSigmaPr();
1025+
const double nsigmaTOF = track.tofNSigmaPr();
1026+
const double pt = track.pt();
1027+
1028+
// Apply TPC PID cut (with systematic variations)
1029+
if (nsigmaTPC < minSigTPC || nsigmaTPC > maxSigTPC)
1030+
return false;
1031+
1032+
// Low-pt: TPC PID is sufficient
1033+
if (pt < PtThreshold)
1034+
return true;
1035+
1036+
// High-pt: require valid TOF match and pass TOF PID (with systematic variations)
1037+
return (track.hasTOF() && nsigmaTOF > minSigTOF && nsigmaTOF < maxSigTOF);
1038+
}
1039+
1040+
// Selection of (anti)deuterons with systematic variations
1041+
template <typename DeuteronTrack>
1042+
bool isDeuteronSyst(const DeuteronTrack& track, double minSigTPC, double maxSigTPC, double minSigTOF, double maxSigTOF)
1043+
{
1044+
// Constant
1045+
static constexpr double PtThreshold = 1.0;
1046+
1047+
// PID variables and transverse momentum of the track
1048+
const double nsigmaTPC = track.tpcNSigmaDe();
1049+
const double nsigmaTOF = track.tofNSigmaDe();
1050+
const double pt = track.pt();
1051+
1052+
// Apply TPC PID cut (with systematic variations)
1053+
if (nsigmaTPC < minSigTPC || nsigmaTPC > maxSigTPC)
1054+
return false;
1055+
1056+
// Low-pt: TPC PID is sufficient
1057+
if (pt < PtThreshold)
1058+
return true;
1059+
1060+
// High-pt: require valid TOF match and pass TOF PID (with systematic variations)
1061+
return (track.hasTOF() && nsigmaTOF > minSigTOF && nsigmaTOF < maxSigTOF);
1062+
}
1063+
9861064
// Process Data
9871065
void processData(SelectedCollisions::iterator const& collision, AntiNucleiTracks const& tracks, aod::BCsWithTimestamps const&)
9881066
{
@@ -3102,6 +3180,196 @@ struct AntinucleiInJets {
31023180
}
31033181
PROCESS_SWITCH(AntinucleiInJets, processCorr, "Process Correlation analysis", false);
31043182

3183+
// Process correlation analysis with systematic variations of analysis parameters
3184+
void processSystCorr(SelectedCollisions::iterator const& collision, AntiNucleiTracks const& tracks)
3185+
{
3186+
// cut settings (from processSystData)
3187+
static std::vector<double> maxDcaxySyst = {
3188+
0.071, 0.060, 0.066, 0.031, 0.052, 0.078, 0.045, 0.064, 0.036, 0.074,
3189+
0.079, 0.043, 0.067, 0.059, 0.032, 0.070, 0.048, 0.077, 0.062, 0.034,
3190+
0.057, 0.055, 0.073, 0.038, 0.050, 0.075, 0.041, 0.061, 0.033, 0.069,
3191+
0.035, 0.044, 0.076, 0.049, 0.037, 0.054, 0.072, 0.046, 0.058, 0.040,
3192+
0.068, 0.042, 0.056, 0.039, 0.047, 0.065, 0.051, 0.053, 0.063, 0.030};
3193+
3194+
static std::vector<double> maxDcazSyst = {
3195+
0.064, 0.047, 0.032, 0.076, 0.039, 0.058, 0.043, 0.069, 0.050, 0.035,
3196+
0.074, 0.061, 0.045, 0.033, 0.068, 0.055, 0.037, 0.071, 0.042, 0.053,
3197+
0.077, 0.038, 0.065, 0.049, 0.036, 0.059, 0.044, 0.067, 0.041, 0.034,
3198+
0.073, 0.052, 0.040, 0.063, 0.046, 0.031, 0.070, 0.054, 0.037, 0.062,
3199+
0.048, 0.035, 0.075, 0.051, 0.039, 0.066, 0.043, 0.060, 0.032, 0.056};
3200+
3201+
static std::vector<double> nSigmaItsMinSyst = {
3202+
-2.9, -2.8, -3.0, -3.4, -2.7, -3.3, -3.0, -3.1, -3.4, -3.1,
3203+
-3.0, -2.8, -3.2, -2.6, -2.7, -3.4, -2.9, -3.0, -3.0, -2.7,
3204+
-2.9, -3.3, -3.0, -3.1, -3.2, -3.0, -2.9, -2.7, -3.3, -3.0,
3205+
-2.8, -3.3, -2.7, -3.3, -2.8, -3.4, -2.8, -3.4, -2.9, -3.1,
3206+
-3.2, -2.6, -3.1, -2.9, -3.1, -2.8, -2.9, -3.3, -3.0, -2.8};
3207+
3208+
static std::vector<double> nSigmaItsMaxSyst = {
3209+
2.9, 2.8, 3.0, 3.4, 2.7, 3.3, 3.0, 3.1, 3.4, 3.1, 3.0, 2.8, 3.2,
3210+
2.6, 2.7, 3.4, 2.9, 3.0, 3.0, 2.7, 2.9, 3.3, 3.0, 3.1, 3.2, 3.0,
3211+
2.9, 2.7, 3.3, 3.0, 2.8, 3.3, 2.7, 3.3, 2.8, 3.4, 2.8, 3.4, 2.9,
3212+
3.1, 3.2, 2.6, 3.1, 2.9, 3.1, 2.8, 2.9, 3.3, 3.0, 2.8};
3213+
3214+
static std::vector<double> minNsigmaTpcSyst = {
3215+
-3.2, -2.9, -3.1, -2.9, -3.5, -2.6, -3.3, -3.0, -3.5, -2.7,
3216+
-3.0, -2.6, -3.3, -3.4, -2.8, -3.1, -2.6, -3.2, -3.1, -2.8,
3217+
-3.4, -2.7, -3.4, -2.9, -3.0, -2.5, -3.3, -2.8, -3.1, -2.7,
3218+
-3.4, -2.8, -3.3, -2.6, -3.1, -2.5, -3.4, -3.0, -3.2, -2.6,
3219+
-3.4, -2.8, -3.1, -2.6, -3.3, -2.7, -3.2, -2.7, -3.4, -2.9};
3220+
3221+
static std::vector<double> maxNsigmaTpcSyst = {
3222+
3.2, 2.9, 3.1, 2.9, 3.5, 2.6, 3.3, 3.0, 3.5, 2.7, 3.0, 2.6, 3.3,
3223+
3.4, 2.8, 3.1, 2.6, 3.2, 3.1, 2.8, 3.4, 2.7, 3.4, 2.9, 3.0, 2.5,
3224+
3.3, 2.8, 3.1, 2.7, 3.4, 2.8, 3.3, 2.6, 3.1, 2.5, 3.4, 3.0, 3.2,
3225+
2.6, 3.4, 2.8, 3.1, 2.6, 3.3, 2.7, 3.2, 2.7, 3.4, 2.9};
3226+
3227+
static std::vector<double> minNsigmaTofSyst = {
3228+
-3.2, -2.9, -3.1, -2.9, -3.5, -2.6, -3.3, -3.0, -3.5, -2.7,
3229+
-3.0, -2.6, -3.3, -3.4, -2.8, -3.1, -2.6, -3.2, -3.1, -2.8,
3230+
-3.4, -2.7, -3.4, -2.9, -3.0, -2.5, -3.3, -2.8, -3.1, -2.7,
3231+
-3.4, -2.8, -3.3, -2.6, -3.1, -2.5, -3.4, -3.0, -3.2, -2.6,
3232+
-3.4, -2.8, -3.1, -2.6, -3.3, -2.7, -3.2, -2.7, -3.4, -2.9};
3233+
3234+
static std::vector<double> maxNsigmaTofSyst = {
3235+
3.9, 3.6, 3.8, 3.2, 3.2, 3.5, 3.1, 3.8, 3.5, 3.4, 3.9, 3.8, 3.7,
3236+
3.0, 3.6, 3.1, 3.7, 3.4, 4.0, 3.0, 3.7, 3.3, 3.9, 3.1, 3.3, 3.5,
3237+
3.6, 3.2, 3.5, 3.3, 3.9, 3.0, 3.4, 3.2, 3.1, 3.9, 3.6, 3.1, 3.2,
3238+
4.0, 3.1, 3.7, 3.6, 3.1, 3.3, 3.5, 3.3, 3.4, 3.1, 3.8};
3239+
3240+
// Event counter: before event selection
3241+
registryCorr.fill(HIST("eventCounter_syst"), 0.5);
3242+
3243+
// Apply standard event selection (Same as processCorr)
3244+
if (!collision.sel8() || std::fabs(collision.posZ()) > zVtx)
3245+
return;
3246+
3247+
registryCorr.fill(HIST("eventCounter_syst"), 1.5);
3248+
3249+
if (rejectITSROFBorder && !collision.selection_bit(o2::aod::evsel::kNoITSROFrameBorder))
3250+
return;
3251+
registryCorr.fill(HIST("eventCounter_syst"), 2.5);
3252+
3253+
if (rejectTFBorder && !collision.selection_bit(o2::aod::evsel::kNoTimeFrameBorder))
3254+
return;
3255+
registryCorr.fill(HIST("eventCounter_syst"), 3.5);
3256+
3257+
if (requireVtxITSTPC && !collision.selection_bit(o2::aod::evsel::kIsVertexITSTPC))
3258+
return;
3259+
registryCorr.fill(HIST("eventCounter_syst"), 4.5);
3260+
3261+
if (rejectSameBunchPileup && !collision.selection_bit(o2::aod::evsel::kNoSameBunchPileup))
3262+
return;
3263+
registryCorr.fill(HIST("eventCounter_syst"), 5.5);
3264+
3265+
if (requireIsGoodZvtxFT0VsPV && !collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV))
3266+
return;
3267+
registryCorr.fill(HIST("eventCounter_syst"), 6.5);
3268+
3269+
if (requireIsVertexTOFmatched && !collision.selection_bit(o2::aod::evsel::kIsVertexTOFmatched))
3270+
return;
3271+
registryCorr.fill(HIST("eventCounter_syst"), 7.5);
3272+
3273+
const float multiplicity = collision.centFT0M();
3274+
3275+
// Bins setup
3276+
static const std::vector<double> ptOverAbins = {0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0};
3277+
const int nBins = ptOverAbins.size() - 1;
3278+
3279+
// Loop over systematic variations
3280+
for (int isyst = 0; isyst < nSyst; isyst++) {
3281+
3282+
// Fill event counter for this systematic
3283+
registryCorr.fill(HIST("eventCounter_centrality_fullEvent_syst"), multiplicity, static_cast<double>(isyst));
3284+
3285+
// Particle counters for this specific cut setting
3286+
std::vector<int> nAntiprotonFullEvent(nBins, 0);
3287+
std::vector<int> nAntideuteronFullEvent(nBins, 0);
3288+
int nTotProtonFullEvent(0);
3289+
int nTotDeuteronFullEvent(0);
3290+
int nTotAntiprotonFullEvent(0);
3291+
int nTotAntideuteronFullEvent(0);
3292+
3293+
// Loop over tracks
3294+
for (auto const& track : tracks) {
3295+
3296+
// Apply track selection (from processSystData)
3297+
if (!passedTrackSelectionSyst(track, isyst))
3298+
continue;
3299+
3300+
// Apply DCA selections (from processSystData)
3301+
if (std::fabs(track.dcaXY()) > maxDcaxySyst[isyst] || std::fabs(track.dcaZ()) > maxDcazSyst[isyst])
3302+
continue;
3303+
3304+
// Particle identification using the ITS cluster size (vary also PID ITS) (from processSystData)
3305+
bool passedItsPidProt(true), passedItsPidDeut(true);
3306+
double nSigmaITSprot = static_cast<double>(itsResponse.nSigmaITS<o2::track::PID::Proton>(track));
3307+
double nSigmaITSdeut = static_cast<double>(itsResponse.nSigmaITS<o2::track::PID::Deuteron>(track));
3308+
3309+
if (applyItsPid && track.pt() < ptMaxItsPidProt && (nSigmaITSprot < nSigmaItsMinSyst[isyst] || nSigmaITSprot > nSigmaItsMaxSyst[isyst])) {
3310+
passedItsPidProt = false;
3311+
}
3312+
if (applyItsPid && track.pt() < ptMaxItsPidDeut && (nSigmaITSdeut < nSigmaItsMinSyst[isyst] || nSigmaITSdeut > nSigmaItsMaxSyst[isyst])) {
3313+
passedItsPidDeut = false;
3314+
}
3315+
3316+
// Check Identity (Proton/Deuteron) using TPC/TOF
3317+
bool isProt = isProtonSyst(track, minNsigmaTpcSyst[isyst], maxNsigmaTpcSyst[isyst], minNsigmaTofSyst[isyst], maxNsigmaTofSyst[isyst]);
3318+
bool isDeut = isDeuteronSyst(track, minNsigmaTpcSyst[isyst], maxNsigmaTpcSyst[isyst], minNsigmaTofSyst[isyst], maxNsigmaTofSyst[isyst]);
3319+
3320+
// Kinematic range selection & counting
3321+
// (Anti)protons
3322+
if (isProt && passedItsPidProt) {
3323+
if (track.pt() >= ptOverAbins[0] && track.pt() < ptOverAbins[nBins]) {
3324+
if (track.sign() > 0) {
3325+
nTotProtonFullEvent++;
3326+
} else if (track.sign() < 0) {
3327+
nTotAntiprotonFullEvent++;
3328+
int ibin = findBin(ptOverAbins, track.pt());
3329+
nAntiprotonFullEvent[ibin]++;
3330+
}
3331+
}
3332+
}
3333+
3334+
// (Anti)deuterons
3335+
if (isDeut && passedItsPidDeut) {
3336+
double ptPerNucleon = 0.5 * track.pt();
3337+
if (ptPerNucleon >= ptOverAbins[0] && ptPerNucleon < ptOverAbins[nBins]) {
3338+
if (track.sign() > 0) {
3339+
nTotDeuteronFullEvent++;
3340+
} else if (track.sign() < 0) {
3341+
nTotAntideuteronFullEvent++;
3342+
int ibin = findBin(ptOverAbins, ptPerNucleon);
3343+
nAntideuteronFullEvent[ibin]++;
3344+
}
3345+
}
3346+
}
3347+
}
3348+
3349+
// Fill Correlation Histograms for systematic variations
3350+
int netProtonFullEvent = nTotProtonFullEvent - nTotAntiprotonFullEvent;
3351+
int netDeuteronFullEvent = nTotDeuteronFullEvent - nTotAntideuteronFullEvent;
3352+
3353+
registryCorr.fill(HIST("rho_fullEvent_syst"), nTotAntideuteronFullEvent, nTotAntiprotonFullEvent, multiplicity, static_cast<double>(isyst));
3354+
registryCorr.fill(HIST("rho_netP_netD_fullEvent_syst"), netDeuteronFullEvent, netProtonFullEvent, static_cast<double>(isyst));
3355+
3356+
// Fill efficiency histograms for systematic variations
3357+
for (int i = 0; i < nBins; i++) {
3358+
double ptAcenteri = 0.5 * (ptOverAbins[i] + ptOverAbins[i + 1]);
3359+
3360+
registryCorr.fill(HIST("q1d_fullEvent_syst"), nAntideuteronFullEvent[i], ptAcenteri, multiplicity, static_cast<double>(isyst));
3361+
registryCorr.fill(HIST("q1p_fullEvent_syst"), nAntiprotonFullEvent[i], ptAcenteri, multiplicity, static_cast<double>(isyst));
3362+
for (int j = 0; j < nBins; j++) {
3363+
double ptAcenterj = 0.5 * (ptOverAbins[j] + ptOverAbins[j + 1]);
3364+
registryCorr.fill(HIST("q1d_square_fullEvent_syst"), ptAcenteri, ptAcenterj, (nAntideuteronFullEvent[i] * nAntideuteronFullEvent[j]), multiplicity, static_cast<double>(isyst));
3365+
registryCorr.fill(HIST("q1p_square_fullEvent_syst"), ptAcenteri, ptAcenterj, (nAntiprotonFullEvent[i] * nAntiprotonFullEvent[j]), multiplicity, static_cast<double>(isyst));
3366+
registryCorr.fill(HIST("q1d_q1p_fullEvent_syst"), ptAcenteri, ptAcenterj, (nAntideuteronFullEvent[i] * nAntiprotonFullEvent[j]), multiplicity, static_cast<double>(isyst));
3367+
}
3368+
}
3369+
}
3370+
}
3371+
PROCESS_SWITCH(AntinucleiInJets, processSystCorr, "Process Correlation systematics", false);
3372+
31053373
// Process coalescence
31063374
void processCoalescence(GenCollisionsMc const& collisions, aod::McParticles const& mcParticles)
31073375
{
@@ -3206,6 +3474,19 @@ struct AntinucleiInJets {
32063474
}
32073475
}
32083476

3477+
for (const auto& part : chargedParticles) {
3478+
if (part.used)
3479+
continue;
3480+
3481+
// Fill histograms for Full Event
3482+
if (part.pdgCode == PDG_t::kProtonBar) {
3483+
registryMC.fill(HIST("antiproton_coal_fullEvent"), part.pt());
3484+
}
3485+
if (part.pdgCode == -o2::constants::physics::Pdg::kDeuteron) {
3486+
registryMC.fill(HIST("antideuteron_coal_fullEvent"), part.pt());
3487+
}
3488+
}
3489+
32093490
// Fill particle array to feed to Fastjet
32103491
for (const auto& part : chargedParticles) {
32113492

0 commit comments

Comments
 (0)