|
50 | 50 | #include "Math/VectorUtil.h" |
51 | 51 | #include "TGeoGlobalMagField.h" |
52 | 52 | #include "TH3F.h" |
| 53 | +#include "THn.h" |
53 | 54 | #include "TRandom.h" |
54 | 55 | #include <TObject.h> |
55 | 56 | #include <TString.h> |
@@ -1156,11 +1157,12 @@ class VarManager : public TObject |
1156 | 1157 | static void SetCalibrationType(int type, bool useInterpolation = true) |
1157 | 1158 | { |
1158 | 1159 | if (type < 0 || type > 2) { |
1159 | | - LOG(fatal) << "Invalid calibration type. Must be 0, 1, or 2."); |
| 1160 | + LOG(fatal) << "Invalid calibration type. Must be 0, 1, or 2."; |
1160 | 1161 | } |
1161 | 1162 | fgCalibrationType = type; |
1162 | 1163 | fgUseInterpolatedCalibration = useInterpolation; |
1163 | 1164 | } |
| 1165 | + static double ComputePIDcalibration(int species, double nSigmaValue); |
1164 | 1166 |
|
1165 | 1167 | static TObject* GetCalibrationObject(CalibObjects calib) |
1166 | 1168 | { |
@@ -1378,145 +1380,6 @@ o2::dataformats::GlobalFwdTrack VarManager::PropagateMuon(const T& muon, const C |
1378 | 1380 | return propmuon; |
1379 | 1381 | } |
1380 | 1382 |
|
1381 | | -double VarManager::ComputePIDcalibration(int species, double nSigmaValue) { |
1382 | | - // species: 0 - electron, 1 - pion, 2 - kaon, 3 - proton |
1383 | | - // Depending on the PID calibration type, we use different types of calibration histograms |
1384 | | - |
1385 | | - if (fgCalibrationType == 1) { |
1386 | | - // get the calibration histograms |
1387 | | - CalibObjects calibMean, calibSigma; |
1388 | | - switch (species) { |
1389 | | - case 0: |
1390 | | - calibMean = kTPCElectronMean; |
1391 | | - calibSigma = kTPCElectronSigma; |
1392 | | - break; |
1393 | | - case 1: |
1394 | | - calibMean = kTPCPionMean; |
1395 | | - calibSigma = kTPCPionSigma; |
1396 | | - break; |
1397 | | - case 2: |
1398 | | - calibMean = kTPCKaonMean; |
1399 | | - calibSigma = kTPCKaonSigma; |
1400 | | - break; |
1401 | | - case 3: |
1402 | | - calibMean = kTPCProtonMean; |
1403 | | - calibSigma = kTPCProtonSigma; |
1404 | | - break; |
1405 | | - default: |
1406 | | - LOG(fatal) << "Invalid species for PID calibration: " << species; |
1407 | | - return -999.0; // Return zero if species is invalid |
1408 | | - }; |
1409 | | - |
1410 | | - TH3F* calibMeanHist = reinterpret_cast<TH3F*>(fgCalibs[calibMean]); |
1411 | | - TH3F* calibSigmaHist = reinterpret_cast<TH3F*>(fgCalibs[calibSigma]); |
1412 | | - if (!calibMeanHist || !calibSigmaHist) { |
1413 | | - LOG(fatal) << "Calibration histograms not found for species: " << species; |
1414 | | - return -999.0; // Return zero if histograms are not found |
1415 | | - } |
1416 | | - |
1417 | | - // Get the bin indices for the calibration histograms |
1418 | | - int binTPCncls = calibMeanHist->GetXaxis()->FindBin(fgValues[kTPCncls]); |
1419 | | - binTPCncls = (binTPCncls == 0 ? 1 : binTPCncls); |
1420 | | - binTPCncls = (binTPCncls > calibMeanHist->GetXaxis()->GetNbins() ? calibMeanHist->GetXaxis()->GetNbins() : binTPCncls); |
1421 | | - int binPin = calibMeanHist->GetYaxis()->FindBin(fgValues[kPin]); |
1422 | | - binPin = (binPin == 0 ? 1 : binPin); |
1423 | | - binPin = (binPin > calibMeanHist->GetYaxis()->GetNbins() ? calibMeanHist->GetYaxis()->GetNbins() : binPin); |
1424 | | - int binEta = calibMeanHist->GetZaxis()->FindBin(fgValues[kEta]); |
1425 | | - binEta = (binEta == 0 ? 1 : binEta); |
1426 | | - binEta = (binEta > calibMeanHist->GetZaxis()->GetNbins() ? calibMeanHist->GetZaxis()->GetNbins() : binEta); |
1427 | | - |
1428 | | - double mean = calibMeanHist->GetBinContent(binTPCncls, binPin, binEta); |
1429 | | - double sigma = calibSigmaHist->GetBinContent(binTPCncls, binPin, binEta); |
1430 | | - return (nSigmaValue - mean) / sigma; // Return the calibrated nSigma value |
1431 | | - } |
1432 | | - else if (fgCalibrationType == 2) { |
1433 | | - // get the calibration histograms |
1434 | | - CalibObjects calibMean, calibSigma, calibStatus; |
1435 | | - switch (species) { |
1436 | | - case 0: |
1437 | | - calibMean = kTPCElectronMean; |
1438 | | - calibSigma = kTPCElectronSigma; |
1439 | | - calibStatus = kTPCElectronStatus; |
1440 | | - break; |
1441 | | - case 1: |
1442 | | - calibMean = kTPCPionMean; |
1443 | | - calibSigma = kTPCPionSigma; |
1444 | | - calibStatus = kTPCPionStatus; |
1445 | | - break; |
1446 | | - case 2: |
1447 | | - calibMean = kTPCKaonMean; |
1448 | | - calibSigma = kTPCKaonSigma; |
1449 | | - calibStatus = kTPCKaonStatus; |
1450 | | - break; |
1451 | | - case 3: |
1452 | | - calibMean = kTPCProtonMean; |
1453 | | - calibSigma = kTPCProtonSigma; |
1454 | | - calibStatus = kTPCProtonStatus; |
1455 | | - break; |
1456 | | - default: |
1457 | | - LOG(fatal) << "Invalid species for PID calibration: " << species; |
1458 | | - return -999.0; // Return zero if species is invalid |
1459 | | - }; |
1460 | | - |
1461 | | - THnF* calibMeanHist = reinterpret_cast<THnF*>(fgCalibs[calibMean]); |
1462 | | - THnF* calibSigmaHist = reinterpret_cast<THnF*>(fgCalibs[calibSigma]); |
1463 | | - THnF* calibStatusHist = reinterpret_cast<THnF*>(fgCalibs[calibStatus]); |
1464 | | - if (!calibMeanHist || !calibSigmaHist || !calibStatusHist) { |
1465 | | - LOG(fatal) << "Calibration histograms not found for species: " << species; |
1466 | | - return -999.0; // Return zero if histograms are not found |
1467 | | - } |
1468 | | - |
1469 | | - // Get the bin indices for the calibration histograms |
1470 | | - int binEta = calibMeanHist->GetAxis(0)->FindBin(fgValues[kEta]); |
1471 | | - binEta = (binEta == 0 ? 1 : binEta); |
1472 | | - binEta = (binEta > calibMeanHist->GetAxis(0)->GetNbins() ? calibMeanHist->GetAxis(0)->GetNbins() : binEta); |
1473 | | - int binNpv = calibMeanHist->GetAxis(1)->FindBin(fgValues[kVtxNcontribReal]); |
1474 | | - binNpv = (binNpv == 0 ? 1 : binNpv); |
1475 | | - binNpv = (binNpv > calibMeanHist->GetAxis(1)->GetNbins() ? calibMeanHist->GetAxis(1)->GetNbins() : binNpv); |
1476 | | - int binNlong = calibMeanHist->GetAxis(2)->FindBin(fgValues[kNTPCcontribLongA]); |
1477 | | - binNlong = (binNlong == 0 ? 1 : binNlong); |
1478 | | - binNlong = (binNlong > calibMeanHist->GetAxis(2)->GetNbins() ? calibMeanHist->GetAxis(2)->GetNbins() : binNlong); |
1479 | | - int binTlong = calibMeanHist->GetAxis(3)->FindBin(fgValues[kNTPCmedianTimeLongA]); |
1480 | | - binTlong = (binTlong == 0 ? 1 : binTlong); |
1481 | | - binTlong = (binTlong > calibMeanHist->GetAxis(3)->GetNbins() ? calibMeanHist->GetAxis(3)->GetNbins() : binTlong); |
1482 | | - |
1483 | | - int bin[4] = {binEta, binNpv, binNlong, binTlong}; |
1484 | | - int status = reinterpret_cast<int>(calibStatusHist->GetBinContent(bin)); |
1485 | | - double mean = calibMeanHist->GetBinContent(bin); |
1486 | | - double sigma = calibSigmaHist->GetBinContent(bin); |
1487 | | - switch (status) { |
1488 | | - case 0: |
1489 | | - // good calibration, return the calibrated nSigma value |
1490 | | - return (nSigmaValue - mean) / sigma; |
1491 | | - break; |
1492 | | - case 1: |
1493 | | - // calibration not valid, return the original nSigma value |
1494 | | - return nSigmaValue; |
1495 | | - break; |
1496 | | - case 2: // calibration constant has poor stat uncertainty, consider the user option for what to do |
1497 | | - case 3: |
1498 | | - // calibration constants have been interpolated |
1499 | | - if (fgUseInterpolatedCalibration) { |
1500 | | - return (nSigmaValue - mean) / sigma; |
1501 | | - } else { |
1502 | | - // return the original nSigma value |
1503 | | - return nSigmaValue; |
1504 | | - } |
1505 | | - break; |
1506 | | - case 4: |
1507 | | - // calibration constants interpolation failed, return the original nSigma value |
1508 | | - return nSigmaValue; |
1509 | | - break; |
1510 | | - default: |
1511 | | - return nSigmaValue; // unknown status, return the original nSigma value |
1512 | | - break; |
1513 | | - }; |
1514 | | - } else { |
1515 | | - // unknown calibration type, return the original nSigma value |
1516 | | - LOG(fatal) << "Unknown calibration type: " << fgCalibrationType; |
1517 | | - return nSigmaValue; // Return the original nSigma value |
1518 | | - } |
1519 | | -} |
1520 | 1383 |
|
1521 | 1384 | template <uint32_t fillMap, typename T, typename C> |
1522 | 1385 | void VarManager::FillMuonPDca(const T& muon, const C& collision, float* values) |
|
0 commit comments