Skip to content

[PWGDQ] Implemented a postcalibration scheme for electrons in Pb-Pb collisions, and some cleanup #12679

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
lucamicheletti93 merged 6 commits into from
Aug 22, 2025
Merged

[PWGDQ] Implemented a postcalibration scheme for electrons in Pb-Pb collisions, and some cleanup #12679

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
581ad1b
Adding a guard if constexpr for TOF n-sigma's in the VarManager
Jun 4, 2025
519d5d6
Removed TPC randomized pid quantities from VarManager
Aug 13, 2025
dbeda0c
Implemented a new TPC pid postcalibration type for PbPb collisions
Aug 13, 2025
c3ed238
fix a few compilation errors
Aug 21, 2025
4e1b512
Please consider the following formatting changes
alibuild Aug 21, 2025
3cb240c
Merge pull request #25 from alibuild/alibot-cleanup-12679
iarsene Aug 21, 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
49 changes: 0 additions & 49 deletions PWGDQ/Core/CutsLibrary.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@

#include <RtypesCore.h>

#include <iostream>

Check failure on line 23 in PWGDQ/Core/CutsLibrary.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[include-iostream] 

Do not include iostream. Use O2 logging instead.
#include <set>
#include <string>
#include <vector>
Expand Down Expand Up @@ -1169,7 +1169,7 @@
return cut;
}

for (int iCut = 0; iCut < 10; iCut++) {

Check failure on line 1172 in PWGDQ/Core/CutsLibrary.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number] 

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
if (!nameStr.compare(Form("jpsiEleSel%d_ionut", iCut))) {
cut->AddCut(GetAnalysisCut("kineJpsiEle_ionut"));
cut->AddCut(GetAnalysisCut("dcaCut1_ionut"));
Expand Down Expand Up @@ -1469,41 +1469,13 @@
return cut;
}

for (int i = 1; i <= 8; i++) {

Check failure on line 1472 in PWGDQ/Core/CutsLibrary.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number] 

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
if (!nameStr.compare(Form("dalitzSelected%d", i))) {
cut->AddCut(GetAnalysisCut(Form("dalitzLeg%d", i)));
return cut;
}
}

//---------------------------------------------------------------------------------------
// NOTE: Below there are several TPC pid cuts used for studies of the dE/dx degradation
// and its impact on the high lumi pp quarkonia triggers
// To be removed when not needed anymore
if (!nameStr.compare("jpsiPID1Randomized")) {
cut->AddCut(GetAnalysisCut("jpsiStandardKine")); // standard kine cuts usually are applied via Filter in the task
cut->AddCut(GetAnalysisCut("electronStandardQuality"));
cut->AddCut(GetAnalysisCut("standardPrimaryTrack"));
cut->AddCut(GetAnalysisCut("electronPID1randomized"));
return cut;
}

if (!nameStr.compare("jpsiPID2Randomized")) {
cut->AddCut(GetAnalysisCut("jpsiStandardKine"));
cut->AddCut(GetAnalysisCut("electronStandardQuality"));
cut->AddCut(GetAnalysisCut("standardPrimaryTrack"));
cut->AddCut(GetAnalysisCut("electronPID2randomized"));
return cut;
}

if (!nameStr.compare("jpsiPIDnsigmaRandomized")) {
cut->AddCut(GetAnalysisCut("jpsiStandardKine"));
cut->AddCut(GetAnalysisCut("electronStandardQuality"));
cut->AddCut(GetAnalysisCut("standardPrimaryTrack"));
cut->AddCut(GetAnalysisCut("electronPIDnsigmaRandomized"));
return cut;
}

if (!nameStr.compare("jpsiPIDworseRes")) {
cut->AddCut(GetAnalysisCut("jpsiStandardKine"));
cut->AddCut(GetAnalysisCut("electronStandardQuality"));
Expand Down Expand Up @@ -1971,7 +1943,7 @@
return cut;
}

for (unsigned int i = 0; i < 30; i++) {

Check failure on line 1946 in PWGDQ/Core/CutsLibrary.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number] 

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
if (!nameStr.compare(Form("ElSelCutVar%s%i", vecPIDcase.at(icase).Data(), i))) {
cut->AddCut(GetAnalysisCut("lmeeStandardKine"));
cut->AddCut(GetAnalysisCut(Form("lmeeCutVarTrackCuts%i", i)));
Expand Down Expand Up @@ -2741,7 +2713,7 @@
return cut;
}

for (int i = 1; i <= 8; i++) {

Check failure on line 2716 in PWGDQ/Core/CutsLibrary.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number] 

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
if (!nameStr.compare(Form("lmee%s_pp502TeV_PID%s_UsePrefilter%d", vecTypetrackWithPID.at(jcase).Data(), vecPIDcase.at(icase).Data(), i))) {
cut->AddCut(GetAnalysisCut(Form("notDalitzLeg%d", i)));
cut->AddCut(GetAnalysisCut("lmeeStandardKine"));
Expand Down Expand Up @@ -2798,7 +2770,7 @@
return cut;
}

for (int i = 1; i <= 8; i++) {

Check failure on line 2773 in PWGDQ/Core/CutsLibrary.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number] 

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
if (!nameStr.compare(Form("lmee%s_eNSigmaRun3%s_UsePrefilter%d", vecTypetrackWithPID.at(jcase).Data(), vecPIDcase.at(icase).Data(), i))) {
cut->AddCut(GetAnalysisCut(Form("notDalitzLeg%d", i)));
cut->AddCut(GetAnalysisCut("lmeeStandardKine"));
Expand Down Expand Up @@ -4493,7 +4465,7 @@
cut->AddCut(VarManager::kITSncls, 6.5, 7.5);
cut->AddCut(VarManager::kTPCnclsCR, 80.0, 161.);
cut->AddCut(VarManager::kTPCncls, 90.0, 170.);
} else if (icase == 2) {

Check failure on line 4468 in PWGDQ/Core/CutsLibrary.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number] 

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
cut->AddCut(VarManager::kIsSPDfirst, 0.5, 1.5);
cut->AddCut(VarManager::kITSchi2, 0.0, 5.0);
cut->AddCut(VarManager::kITSncls, 4.5, 7.5);
Expand Down Expand Up @@ -5136,13 +5108,6 @@
return cut;
}

if (!nameStr.compare("electronPID1randomized")) {
cutLow1->SetParameters(130., -40.0);
cut->AddCut(VarManager::kTPCsignalRandomized, 70., 100.);
cut->AddCut(VarManager::kTPCsignalRandomized, cutLow1, 100.0, false, VarManager::kPin, 0.5, 3.0);
return cut;
}

if (!nameStr.compare("electronPID2")) {
cutLow1->SetParameters(130., -40.0);
cut->AddCut(VarManager::kTPCsignal, 73., 100.);
Expand All @@ -5157,13 +5122,6 @@
return cut;
}

if (!nameStr.compare("electronPID2randomized")) {
cutLow1->SetParameters(130., -40.0);
cut->AddCut(VarManager::kTPCsignalRandomized, 73., 100.);
cut->AddCut(VarManager::kTPCsignalRandomized, cutLow1, 100.0, false, VarManager::kPin, 0.5, 3.0);
return cut;
}

if (!nameStr.compare("electronPIDnsigma")) {
cut->AddCut(VarManager::kTPCnSigmaEl, -3.0, 3.0);
cut->AddCut(VarManager::kTPCnSigmaPr, 3.0, 3000.0);
Expand Down Expand Up @@ -5645,13 +5603,6 @@
return cut;
}

if (!nameStr.compare("electronPIDnsigmaRandomized")) {
cut->AddCut(VarManager::kTPCnSigmaElRandomized, -3.0, 3.0);
cut->AddCut(VarManager::kTPCnSigmaPrRandomized, 3.0, 3000.0);
cut->AddCut(VarManager::kTPCnSigmaPiRandomized, 3.0, 3000.0);
return cut;
}

if (!nameStr.compare("electronPIDworseRes")) {
cut->AddCut(VarManager::kTPCnSigmaEl, -3.0, 3.0);
cut->AddCut(VarManager::kTPCnSigmaPr, 3.0 * 0.8, 3000.0); // emulates a 20% degradation in PID resolution
Expand Down Expand Up @@ -7258,11 +7209,11 @@
labelsFlatBin.push_back(Form("%s_cent%.0f_%.0f_pt%.1f_%.1f", cent.c_str(), centMin, centMax, ptMin, ptMax));
LOG(info) << "Added cut for " << Form("%s_cent%.0f_%.0f_pt%.1f_%.1f", cent.c_str(), centMin, centMax, ptMin, ptMax) << " with cuts: [";
for (size_t i = 0; i < binCuts.size(); ++i) {
std::cout << binCuts[i];

Check failure on line 7212 in PWGDQ/Core/CutsLibrary.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[logging] 

Use O2 logging (LOG, LOGF, LOGP).
if (i != binCuts.size() - 1)
std::cout << ", ";

Check failure on line 7214 in PWGDQ/Core/CutsLibrary.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[logging] 

Use O2 logging (LOG, LOGF, LOGP).
}
std::cout << "] and direction: " << (exclude ? "CutGreater" : "CutSmaller") << std::endl;

Check failure on line 7216 in PWGDQ/Core/CutsLibrary.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[logging] 

Use O2 logging (LOG, LOGF, LOGP).
}
}

Expand Down
152 changes: 144 additions & 8 deletions PWGDQ/Core/VarManager.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -45,6 +45,8 @@ o2::vertexing::FwdDCAFitterN<3> VarManager::fgFitterThreeProngFwd;
o2::globaltracking::MatchGlobalFwd VarManager::mMatching;
std::map<VarManager::CalibObjects, TObject*> VarManager::fgCalibs;
bool VarManager::fgRunTPCPostCalibration[4] = {false, false, false, false};
int VarManager::fgCalibrationType = 0; // 0 - no calibration, 1 - calibration vs (TPCncls,pIN,eta) typically for pp, 2 - calibration vs (eta,nPV,nLong,tLong) typically for PbPb
bool VarManager::fgUseInterpolatedCalibration = true; // use interpolated calibration histograms (default: true)

//__________________________________________________________________
VarManager::VarManager() : TObject()
Expand Down Expand Up @@ -209,6 +211,148 @@ float VarManager::calculateCosPA(KFParticle kfp, KFParticle PV)
{
return cpaFromKF(kfp, PV);
}

//__________________________________________________________________
double VarManager::ComputePIDcalibration(int species, double nSigmaValue)
{
// species: 0 - electron, 1 - pion, 2 - kaon, 3 - proton
// Depending on the PID calibration type, we use different types of calibration histograms

if (fgCalibrationType == 1) {
// get the calibration histograms
CalibObjects calibMean, calibSigma;
switch (species) {
case 0:
calibMean = kTPCElectronMean;
calibSigma = kTPCElectronSigma;
break;
case 1:
calibMean = kTPCPionMean;
calibSigma = kTPCPionSigma;
break;
case 2:
calibMean = kTPCKaonMean;
calibSigma = kTPCKaonSigma;
break;
case 3:
calibMean = kTPCProtonMean;
calibSigma = kTPCProtonSigma;
break;
default:
LOG(fatal) << "Invalid species for PID calibration: " << species;
return -999.0; // Return zero if species is invalid
};

TH3F* calibMeanHist = reinterpret_cast<TH3F*>(fgCalibs[calibMean]);
TH3F* calibSigmaHist = reinterpret_cast<TH3F*>(fgCalibs[calibSigma]);
if (!calibMeanHist || !calibSigmaHist) {
LOG(fatal) << "Calibration histograms not found for species: " << species;
return -999.0; // Return zero if histograms are not found
}

// Get the bin indices for the calibration histograms
int binTPCncls = calibMeanHist->GetXaxis()->FindBin(fgValues[kTPCncls]);
binTPCncls = (binTPCncls == 0 ? 1 : binTPCncls);
binTPCncls = (binTPCncls > calibMeanHist->GetXaxis()->GetNbins() ? calibMeanHist->GetXaxis()->GetNbins() : binTPCncls);
int binPin = calibMeanHist->GetYaxis()->FindBin(fgValues[kPin]);
binPin = (binPin == 0 ? 1 : binPin);
binPin = (binPin > calibMeanHist->GetYaxis()->GetNbins() ? calibMeanHist->GetYaxis()->GetNbins() : binPin);
int binEta = calibMeanHist->GetZaxis()->FindBin(fgValues[kEta]);
binEta = (binEta == 0 ? 1 : binEta);
binEta = (binEta > calibMeanHist->GetZaxis()->GetNbins() ? calibMeanHist->GetZaxis()->GetNbins() : binEta);

double mean = calibMeanHist->GetBinContent(binTPCncls, binPin, binEta);
double sigma = calibSigmaHist->GetBinContent(binTPCncls, binPin, binEta);
return (nSigmaValue - mean) / sigma; // Return the calibrated nSigma value
} else if (fgCalibrationType == 2) {
// get the calibration histograms
CalibObjects calibMean, calibSigma, calibStatus;
switch (species) {
case 0:
calibMean = kTPCElectronMean;
calibSigma = kTPCElectronSigma;
calibStatus = kTPCElectronStatus;
break;
case 1:
calibMean = kTPCPionMean;
calibSigma = kTPCPionSigma;
calibStatus = kTPCPionStatus;
break;
case 2:
calibMean = kTPCKaonMean;
calibSigma = kTPCKaonSigma;
calibStatus = kTPCKaonStatus;
break;
case 3:
calibMean = kTPCProtonMean;
calibSigma = kTPCProtonSigma;
calibStatus = kTPCProtonStatus;
break;
default:
LOG(fatal) << "Invalid species for PID calibration: " << species;
return -999.0; // Return zero if species is invalid
};

THnF* calibMeanHist = reinterpret_cast<THnF*>(fgCalibs[calibMean]);
THnF* calibSigmaHist = reinterpret_cast<THnF*>(fgCalibs[calibSigma]);
THnF* calibStatusHist = reinterpret_cast<THnF*>(fgCalibs[calibStatus]);
if (!calibMeanHist || !calibSigmaHist || !calibStatusHist) {
LOG(fatal) << "Calibration histograms not found for species: " << species;
return -999.0; // Return zero if histograms are not found
}

// Get the bin indices for the calibration histograms
int binEta = calibMeanHist->GetAxis(0)->FindBin(fgValues[kEta]);
binEta = (binEta == 0 ? 1 : binEta);
binEta = (binEta > calibMeanHist->GetAxis(0)->GetNbins() ? calibMeanHist->GetAxis(0)->GetNbins() : binEta);
int binNpv = calibMeanHist->GetAxis(1)->FindBin(fgValues[kVtxNcontribReal]);
binNpv = (binNpv == 0 ? 1 : binNpv);
binNpv = (binNpv > calibMeanHist->GetAxis(1)->GetNbins() ? calibMeanHist->GetAxis(1)->GetNbins() : binNpv);
int binNlong = calibMeanHist->GetAxis(2)->FindBin(fgValues[kNTPCcontribLongA]);
binNlong = (binNlong == 0 ? 1 : binNlong);
binNlong = (binNlong > calibMeanHist->GetAxis(2)->GetNbins() ? calibMeanHist->GetAxis(2)->GetNbins() : binNlong);
int binTlong = calibMeanHist->GetAxis(3)->FindBin(fgValues[kNTPCmedianTimeLongA]);
binTlong = (binTlong == 0 ? 1 : binTlong);
binTlong = (binTlong > calibMeanHist->GetAxis(3)->GetNbins() ? calibMeanHist->GetAxis(3)->GetNbins() : binTlong);

int bin[4] = {binEta, binNpv, binNlong, binTlong};
int status = static_cast<int>(calibStatusHist->GetBinContent(bin));
double mean = calibMeanHist->GetBinContent(bin);
double sigma = calibSigmaHist->GetBinContent(bin);
switch (status) {
case 0:
// good calibration, return the calibrated nSigma value
return (nSigmaValue - mean) / sigma;
break;
case 1:
// calibration not valid, return the original nSigma value
return nSigmaValue;
break;
case 2: // calibration constant has poor stat uncertainty, consider the user option for what to do
case 3:
// calibration constants have been interpolated
if (fgUseInterpolatedCalibration) {
return (nSigmaValue - mean) / sigma;
} else {
// return the original nSigma value
return nSigmaValue;
}
break;
case 4:
// calibration constants interpolation failed, return the original nSigma value
return nSigmaValue;
break;
default:
return nSigmaValue; // unknown status, return the original nSigma value
break;
};
} else {
// unknown calibration type, return the original nSigma value
LOG(fatal) << "Unknown calibration type: " << fgCalibrationType;
return nSigmaValue; // Return the original nSigma value
}
}

//__________________________________________________________________
void VarManager::SetDefaultVarNames()
{
Expand Down Expand Up @@ -1451,8 +1595,6 @@ void VarManager::SetDefaultVarNames()
fgVarNamesMap["kTPCnCRoverFindCls"] = kTPCnCRoverFindCls;
fgVarNamesMap["kTPCchi2"] = kTPCchi2;
fgVarNamesMap["kTPCsignal"] = kTPCsignal;
fgVarNamesMap["kTPCsignalRandomized"] = kTPCsignalRandomized;
fgVarNamesMap["kTPCsignalRandomizedDelta"] = kTPCsignalRandomizedDelta;
fgVarNamesMap["kPhiTPCOuter"] = kPhiTPCOuter;
fgVarNamesMap["kTrackIsInsideTPCModule"] = kTrackIsInsideTPCModule;
fgVarNamesMap["kTRDsignal"] = kTRDsignal;
Expand All @@ -1476,20 +1618,14 @@ void VarManager::SetDefaultVarNames()
fgVarNamesMap["kTrackCTglTgl"] = kTrackCTglTgl;
fgVarNamesMap["kTrackC1Pt21Pt2"] = kTrackC1Pt21Pt2;
fgVarNamesMap["kTPCnSigmaEl"] = kTPCnSigmaEl;
fgVarNamesMap["kTPCnSigmaElRandomized"] = kTPCnSigmaElRandomized;
fgVarNamesMap["kTPCnSigmaElRandomizedDelta"] = kTPCnSigmaElRandomizedDelta;
fgVarNamesMap["kTPCnSigmaMu"] = kTPCnSigmaMu;
fgVarNamesMap["kTPCnSigmaPi"] = kTPCnSigmaPi;
fgVarNamesMap["kTPCnSigmaPiRandomized"] = kTPCnSigmaPiRandomized;
fgVarNamesMap["kTPCnSigmaPiRandomizedDelta"] = kTPCnSigmaPiRandomizedDelta;
fgVarNamesMap["kTPCnSigmaKa"] = kTPCnSigmaKa;
fgVarNamesMap["kTPCnSigmaPr"] = kTPCnSigmaPr;
fgVarNamesMap["kTPCnSigmaEl_Corr"] = kTPCnSigmaEl_Corr;
fgVarNamesMap["kTPCnSigmaPi_Corr"] = kTPCnSigmaPi_Corr;
fgVarNamesMap["kTPCnSigmaKa_Corr"] = kTPCnSigmaKa_Corr;
fgVarNamesMap["kTPCnSigmaPr_Corr"] = kTPCnSigmaPr_Corr;
fgVarNamesMap["kTPCnSigmaPrRandomized"] = kTPCnSigmaPrRandomized;
fgVarNamesMap["kTPCnSigmaPrRandomizedDelta"] = kTPCnSigmaPrRandomizedDelta;
fgVarNamesMap["kTOFnSigmaEl"] = kTOFnSigmaEl;
fgVarNamesMap["kTOFnSigmaMu"] = kTOFnSigmaMu;
fgVarNamesMap["kTOFnSigmaPi"] = kTOFnSigmaPi;
Expand Down
Loading
Loading