Skip to content

Commit c7cc561

Browse files
iarseneIonut Cristian Arsenealibuild
authored
[PWGDQ] Implemented a postcalibration scheme for electrons in Pb-Pb collisions, and some cleanup (#12679)
Co-authored-by: Ionut Cristian Arsene <iarsene@cern.ch> Co-authored-by: ALICE Action Bot <alibuild@cern.ch>
1 parent 8e59519 commit c7cc561

File tree

5 files changed

+193
-146
lines changed

5 files changed

+193
-146
lines changed

‎PWGDQ/Core/CutsLibrary.cxx‎

Lines changed: 0 additions & 49 deletions
Original file line numberDiff line numberDiff line change
@@ -1476,34 +1476,6 @@ AnalysisCompositeCut* o2::aod::dqcuts::GetCompositeCut(const char* cutName)
14761476
}
14771477
}
14781478

1479-
//---------------------------------------------------------------------------------------
1480-
// NOTE: Below there are several TPC pid cuts used for studies of the dE/dx degradation
1481-
// and its impact on the high lumi pp quarkonia triggers
1482-
// To be removed when not needed anymore
1483-
if (!nameStr.compare("jpsiPID1Randomized")) {
1484-
cut->AddCut(GetAnalysisCut("jpsiStandardKine")); // standard kine cuts usually are applied via Filter in the task
1485-
cut->AddCut(GetAnalysisCut("electronStandardQuality"));
1486-
cut->AddCut(GetAnalysisCut("standardPrimaryTrack"));
1487-
cut->AddCut(GetAnalysisCut("electronPID1randomized"));
1488-
return cut;
1489-
}
1490-
1491-
if (!nameStr.compare("jpsiPID2Randomized")) {
1492-
cut->AddCut(GetAnalysisCut("jpsiStandardKine"));
1493-
cut->AddCut(GetAnalysisCut("electronStandardQuality"));
1494-
cut->AddCut(GetAnalysisCut("standardPrimaryTrack"));
1495-
cut->AddCut(GetAnalysisCut("electronPID2randomized"));
1496-
return cut;
1497-
}
1498-
1499-
if (!nameStr.compare("jpsiPIDnsigmaRandomized")) {
1500-
cut->AddCut(GetAnalysisCut("jpsiStandardKine"));
1501-
cut->AddCut(GetAnalysisCut("electronStandardQuality"));
1502-
cut->AddCut(GetAnalysisCut("standardPrimaryTrack"));
1503-
cut->AddCut(GetAnalysisCut("electronPIDnsigmaRandomized"));
1504-
return cut;
1505-
}
1506-
15071479
if (!nameStr.compare("jpsiPIDworseRes")) {
15081480
cut->AddCut(GetAnalysisCut("jpsiStandardKine"));
15091481
cut->AddCut(GetAnalysisCut("electronStandardQuality"));
@@ -5136,13 +5108,6 @@ AnalysisCut* o2::aod::dqcuts::GetAnalysisCut(const char* cutName)
51365108
return cut;
51375109
}
51385110

5139-
if (!nameStr.compare("electronPID1randomized")) {
5140-
cutLow1->SetParameters(130., -40.0);
5141-
cut->AddCut(VarManager::kTPCsignalRandomized, 70., 100.);
5142-
cut->AddCut(VarManager::kTPCsignalRandomized, cutLow1, 100.0, false, VarManager::kPin, 0.5, 3.0);
5143-
return cut;
5144-
}
5145-
51465111
if (!nameStr.compare("electronPID2")) {
51475112
cutLow1->SetParameters(130., -40.0);
51485113
cut->AddCut(VarManager::kTPCsignal, 73., 100.);
@@ -5157,13 +5122,6 @@ AnalysisCut* o2::aod::dqcuts::GetAnalysisCut(const char* cutName)
51575122
return cut;
51585123
}
51595124

5160-
if (!nameStr.compare("electronPID2randomized")) {
5161-
cutLow1->SetParameters(130., -40.0);
5162-
cut->AddCut(VarManager::kTPCsignalRandomized, 73., 100.);
5163-
cut->AddCut(VarManager::kTPCsignalRandomized, cutLow1, 100.0, false, VarManager::kPin, 0.5, 3.0);
5164-
return cut;
5165-
}
5166-
51675125
if (!nameStr.compare("electronPIDnsigma")) {
51685126
cut->AddCut(VarManager::kTPCnSigmaEl, -3.0, 3.0);
51695127
cut->AddCut(VarManager::kTPCnSigmaPr, 3.0, 3000.0);
@@ -5645,13 +5603,6 @@ AnalysisCut* o2::aod::dqcuts::GetAnalysisCut(const char* cutName)
56455603
return cut;
56465604
}
56475605

5648-
if (!nameStr.compare("electronPIDnsigmaRandomized")) {
5649-
cut->AddCut(VarManager::kTPCnSigmaElRandomized, -3.0, 3.0);
5650-
cut->AddCut(VarManager::kTPCnSigmaPrRandomized, 3.0, 3000.0);
5651-
cut->AddCut(VarManager::kTPCnSigmaPiRandomized, 3.0, 3000.0);
5652-
return cut;
5653-
}
5654-
56555606
if (!nameStr.compare("electronPIDworseRes")) {
56565607
cut->AddCut(VarManager::kTPCnSigmaEl, -3.0, 3.0);
56575608
cut->AddCut(VarManager::kTPCnSigmaPr, 3.0 * 0.8, 3000.0); // emulates a 20% degradation in PID resolution

‎PWGDQ/Core/VarManager.cxx‎

Lines changed: 144 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -45,6 +45,8 @@ o2::vertexing::FwdDCAFitterN<3> VarManager::fgFitterThreeProngFwd;
4545
o2::globaltracking::MatchGlobalFwd VarManager::mMatching;
4646
std::map<VarManager::CalibObjects, TObject*> VarManager::fgCalibs;
4747
bool VarManager::fgRunTPCPostCalibration[4] = {false, false, false, false};
48+
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
49+
bool VarManager::fgUseInterpolatedCalibration = true; // use interpolated calibration histograms (default: true)
4850

4951
//__________________________________________________________________
5052
VarManager::VarManager() : TObject()
@@ -209,6 +211,148 @@ float VarManager::calculateCosPA(KFParticle kfp, KFParticle PV)
209211
{
210212
return cpaFromKF(kfp, PV);
211213
}
214+
215+
//__________________________________________________________________
216+
double VarManager::ComputePIDcalibration(int species, double nSigmaValue)
217+
{
218+
// species: 0 - electron, 1 - pion, 2 - kaon, 3 - proton
219+
// Depending on the PID calibration type, we use different types of calibration histograms
220+
221+
if (fgCalibrationType == 1) {
222+
// get the calibration histograms
223+
CalibObjects calibMean, calibSigma;
224+
switch (species) {
225+
case 0:
226+
calibMean = kTPCElectronMean;
227+
calibSigma = kTPCElectronSigma;
228+
break;
229+
case 1:
230+
calibMean = kTPCPionMean;
231+
calibSigma = kTPCPionSigma;
232+
break;
233+
case 2:
234+
calibMean = kTPCKaonMean;
235+
calibSigma = kTPCKaonSigma;
236+
break;
237+
case 3:
238+
calibMean = kTPCProtonMean;
239+
calibSigma = kTPCProtonSigma;
240+
break;
241+
default:
242+
LOG(fatal) << "Invalid species for PID calibration: " << species;
243+
return -999.0; // Return zero if species is invalid
244+
};
245+
246+
TH3F* calibMeanHist = reinterpret_cast<TH3F*>(fgCalibs[calibMean]);
247+
TH3F* calibSigmaHist = reinterpret_cast<TH3F*>(fgCalibs[calibSigma]);
248+
if (!calibMeanHist || !calibSigmaHist) {
249+
LOG(fatal) << "Calibration histograms not found for species: " << species;
250+
return -999.0; // Return zero if histograms are not found
251+
}
252+
253+
// Get the bin indices for the calibration histograms
254+
int binTPCncls = calibMeanHist->GetXaxis()->FindBin(fgValues[kTPCncls]);
255+
binTPCncls = (binTPCncls == 0 ? 1 : binTPCncls);
256+
binTPCncls = (binTPCncls > calibMeanHist->GetXaxis()->GetNbins() ? calibMeanHist->GetXaxis()->GetNbins() : binTPCncls);
257+
int binPin = calibMeanHist->GetYaxis()->FindBin(fgValues[kPin]);
258+
binPin = (binPin == 0 ? 1 : binPin);
259+
binPin = (binPin > calibMeanHist->GetYaxis()->GetNbins() ? calibMeanHist->GetYaxis()->GetNbins() : binPin);
260+
int binEta = calibMeanHist->GetZaxis()->FindBin(fgValues[kEta]);
261+
binEta = (binEta == 0 ? 1 : binEta);
262+
binEta = (binEta > calibMeanHist->GetZaxis()->GetNbins() ? calibMeanHist->GetZaxis()->GetNbins() : binEta);
263+
264+
double mean = calibMeanHist->GetBinContent(binTPCncls, binPin, binEta);
265+
double sigma = calibSigmaHist->GetBinContent(binTPCncls, binPin, binEta);
266+
return (nSigmaValue - mean) / sigma; // Return the calibrated nSigma value
267+
} else if (fgCalibrationType == 2) {
268+
// get the calibration histograms
269+
CalibObjects calibMean, calibSigma, calibStatus;
270+
switch (species) {
271+
case 0:
272+
calibMean = kTPCElectronMean;
273+
calibSigma = kTPCElectronSigma;
274+
calibStatus = kTPCElectronStatus;
275+
break;
276+
case 1:
277+
calibMean = kTPCPionMean;
278+
calibSigma = kTPCPionSigma;
279+
calibStatus = kTPCPionStatus;
280+
break;
281+
case 2:
282+
calibMean = kTPCKaonMean;
283+
calibSigma = kTPCKaonSigma;
284+
calibStatus = kTPCKaonStatus;
285+
break;
286+
case 3:
287+
calibMean = kTPCProtonMean;
288+
calibSigma = kTPCProtonSigma;
289+
calibStatus = kTPCProtonStatus;
290+
break;
291+
default:
292+
LOG(fatal) << "Invalid species for PID calibration: " << species;
293+
return -999.0; // Return zero if species is invalid
294+
};
295+
296+
THnF* calibMeanHist = reinterpret_cast<THnF*>(fgCalibs[calibMean]);
297+
THnF* calibSigmaHist = reinterpret_cast<THnF*>(fgCalibs[calibSigma]);
298+
THnF* calibStatusHist = reinterpret_cast<THnF*>(fgCalibs[calibStatus]);
299+
if (!calibMeanHist || !calibSigmaHist || !calibStatusHist) {
300+
LOG(fatal) << "Calibration histograms not found for species: " << species;
301+
return -999.0; // Return zero if histograms are not found
302+
}
303+
304+
// Get the bin indices for the calibration histograms
305+
int binEta = calibMeanHist->GetAxis(0)->FindBin(fgValues[kEta]);
306+
binEta = (binEta == 0 ? 1 : binEta);
307+
binEta = (binEta > calibMeanHist->GetAxis(0)->GetNbins() ? calibMeanHist->GetAxis(0)->GetNbins() : binEta);
308+
int binNpv = calibMeanHist->GetAxis(1)->FindBin(fgValues[kVtxNcontribReal]);
309+
binNpv = (binNpv == 0 ? 1 : binNpv);
310+
binNpv = (binNpv > calibMeanHist->GetAxis(1)->GetNbins() ? calibMeanHist->GetAxis(1)->GetNbins() : binNpv);
311+
int binNlong = calibMeanHist->GetAxis(2)->FindBin(fgValues[kNTPCcontribLongA]);
312+
binNlong = (binNlong == 0 ? 1 : binNlong);
313+
binNlong = (binNlong > calibMeanHist->GetAxis(2)->GetNbins() ? calibMeanHist->GetAxis(2)->GetNbins() : binNlong);
314+
int binTlong = calibMeanHist->GetAxis(3)->FindBin(fgValues[kNTPCmedianTimeLongA]);
315+
binTlong = (binTlong == 0 ? 1 : binTlong);
316+
binTlong = (binTlong > calibMeanHist->GetAxis(3)->GetNbins() ? calibMeanHist->GetAxis(3)->GetNbins() : binTlong);
317+
318+
int bin[4] = {binEta, binNpv, binNlong, binTlong};
319+
int status = static_cast<int>(calibStatusHist->GetBinContent(bin));
320+
double mean = calibMeanHist->GetBinContent(bin);
321+
double sigma = calibSigmaHist->GetBinContent(bin);
322+
switch (status) {
323+
case 0:
324+
// good calibration, return the calibrated nSigma value
325+
return (nSigmaValue - mean) / sigma;
326+
break;
327+
case 1:
328+
// calibration not valid, return the original nSigma value
329+
return nSigmaValue;
330+
break;
331+
case 2: // calibration constant has poor stat uncertainty, consider the user option for what to do
332+
case 3:
333+
// calibration constants have been interpolated
334+
if (fgUseInterpolatedCalibration) {
335+
return (nSigmaValue - mean) / sigma;
336+
} else {
337+
// return the original nSigma value
338+
return nSigmaValue;
339+
}
340+
break;
341+
case 4:
342+
// calibration constants interpolation failed, return the original nSigma value
343+
return nSigmaValue;
344+
break;
345+
default:
346+
return nSigmaValue; // unknown status, return the original nSigma value
347+
break;
348+
};
349+
} else {
350+
// unknown calibration type, return the original nSigma value
351+
LOG(fatal) << "Unknown calibration type: " << fgCalibrationType;
352+
return nSigmaValue; // Return the original nSigma value
353+
}
354+
}
355+
212356
//__________________________________________________________________
213357
void VarManager::SetDefaultVarNames()
214358
{
@@ -1451,8 +1595,6 @@ void VarManager::SetDefaultVarNames()
14511595
fgVarNamesMap["kTPCnCRoverFindCls"] = kTPCnCRoverFindCls;
14521596
fgVarNamesMap["kTPCchi2"] = kTPCchi2;
14531597
fgVarNamesMap["kTPCsignal"] = kTPCsignal;
1454-
fgVarNamesMap["kTPCsignalRandomized"] = kTPCsignalRandomized;
1455-
fgVarNamesMap["kTPCsignalRandomizedDelta"] = kTPCsignalRandomizedDelta;
14561598
fgVarNamesMap["kPhiTPCOuter"] = kPhiTPCOuter;
14571599
fgVarNamesMap["kTrackIsInsideTPCModule"] = kTrackIsInsideTPCModule;
14581600
fgVarNamesMap["kTRDsignal"] = kTRDsignal;
@@ -1476,20 +1618,14 @@ void VarManager::SetDefaultVarNames()
14761618
fgVarNamesMap["kTrackCTglTgl"] = kTrackCTglTgl;
14771619
fgVarNamesMap["kTrackC1Pt21Pt2"] = kTrackC1Pt21Pt2;
14781620
fgVarNamesMap["kTPCnSigmaEl"] = kTPCnSigmaEl;
1479-
fgVarNamesMap["kTPCnSigmaElRandomized"] = kTPCnSigmaElRandomized;
1480-
fgVarNamesMap["kTPCnSigmaElRandomizedDelta"] = kTPCnSigmaElRandomizedDelta;
14811621
fgVarNamesMap["kTPCnSigmaMu"] = kTPCnSigmaMu;
14821622
fgVarNamesMap["kTPCnSigmaPi"] = kTPCnSigmaPi;
1483-
fgVarNamesMap["kTPCnSigmaPiRandomized"] = kTPCnSigmaPiRandomized;
1484-
fgVarNamesMap["kTPCnSigmaPiRandomizedDelta"] = kTPCnSigmaPiRandomizedDelta;
14851623
fgVarNamesMap["kTPCnSigmaKa"] = kTPCnSigmaKa;
14861624
fgVarNamesMap["kTPCnSigmaPr"] = kTPCnSigmaPr;
14871625
fgVarNamesMap["kTPCnSigmaEl_Corr"] = kTPCnSigmaEl_Corr;
14881626
fgVarNamesMap["kTPCnSigmaPi_Corr"] = kTPCnSigmaPi_Corr;
14891627
fgVarNamesMap["kTPCnSigmaKa_Corr"] = kTPCnSigmaKa_Corr;
14901628
fgVarNamesMap["kTPCnSigmaPr_Corr"] = kTPCnSigmaPr_Corr;
1491-
fgVarNamesMap["kTPCnSigmaPrRandomized"] = kTPCnSigmaPrRandomized;
1492-
fgVarNamesMap["kTPCnSigmaPrRandomizedDelta"] = kTPCnSigmaPrRandomizedDelta;
14931629
fgVarNamesMap["kTOFnSigmaEl"] = kTOFnSigmaEl;
14941630
fgVarNamesMap["kTOFnSigmaMu"] = kTOFnSigmaMu;
14951631
fgVarNamesMap["kTOFnSigmaPi"] = kTOFnSigmaPi;

0 commit comments

Comments
 (0)