Skip to content

Commit 2bdde6c

Browse files
Francesca CasilloFrancesca Casillo
Francesca Casillo
authored and
Francesca Casillo
committed
[PWGLF] Added fullEvent histograms to processCoalescence and a new function (processSystCorr) for systematic studies
1 parent 8eb8a8e commit 2bdde6c

File tree

1 file changed

+244
-0
lines changed

1 file changed

+244
-0
lines changed

‎PWGLF/Tasks/Nuspex/antinucleiInJets.cxx‎

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

32093453
// Fill particle array to feed to Fastjet
32103454
for (const auto& part : chargedParticles) {

0 commit comments

Comments
 (0)