Skip to content

Commit f5ddf80

Browse files
njacazionjacazio
authored
Update onTheFlyRichPid.cxx
1 parent ce21fdb commit f5ddf80

File tree

1 file changed

+41
-62
lines changed

1 file changed

+41
-62
lines changed

‎ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx‎

Lines changed: 41 additions & 62 deletions
Original file line numberDiff line numberDiff line change
@@ -105,27 +105,6 @@ struct OnTheFlyRichPid {
105105
Configurable<bool> flagIncludeTrackAngularRes{"flagIncludeTrackAngularRes", true, "flag to include or exclude track time resolution"};
106106
Configurable<float> multiplicityEtaRange{"multiplicityEtaRange", 0.800000012, "eta range to compute the multiplicity"};
107107
Configurable<bool> flagRICHLoadDelphesLUTs{"flagRICHLoadDelphesLUTs", false, "flag to load Delphes LUTs for tracking correction (use recoTrack parameters if false)"};
108-
/*Configurable<float> bRichRefractiveIndexSector0AndNMinus1{"bRichRefractiveIndexSector0AndMinus1", 1.03, "barrel RICH refractive index sector 0 and N-1"};
109-
Configurable<float> bRichRefractiveIndexSector1AndNMinus2{"bRichRefractiveIndexSector1AndMinus2", 1.03, "barrel RICH refractive index sector 1 and N-2"};
110-
Configurable<float> bRichRefractiveIndexSector2AndNMinus3{"bRichRefractiveIndexSector2AndMinus3", 1.03, "barrel RICH refractive index sector 2 and N-3"};
111-
Configurable<float> bRichRefractiveIndexSector3AndNMinus4{"bRichRefractiveIndexSector3AndMinus4", 1.03, "barrel RICH refractive index sector 3 and N-4"};
112-
Configurable<float> bRichRefractiveIndexSector4AndNMinus5{"bRichRefractiveIndexSector4AndMinus5", 1.03, "barrel RICH refractive index sector 4 and N-5"};
113-
Configurable<float> bRichRefractiveIndexSector5AndNMinus6{"bRichRefractiveIndexSector5AndMinus6", 1.03, "barrel RICH refractive index sector 5 and N-6"};
114-
Configurable<float> bRichRefractiveIndexSector6AndNMinus7{"bRichRefractiveIndexSector6AndMinus7", 1.03, "barrel RICH refractive index sector 6 and N-7"};
115-
Configurable<float> bRichRefractiveIndexSector7AndNMinus8{"bRichRefractiveIndexSector7AndMinus8", 1.03, "barrel RICH refractive index sector 7 and N-8"};
116-
Configurable<float> bRichRefractiveIndexSector8AndNMinus9{"bRichRefractiveIndexSector8AndMinus9", 1.03, "barrel RICH refractive index sector 8 and N-9"};
117-
Configurable<float> bRichRefractiveIndexSector9AndNMinus10{"bRichRefractiveIndexSector9AndMinus10", 1.03, "barrel RICH refractive index sector 9 and N-10"};
118-
Configurable<float> bRichRefractiveIndexSector10AndNMinus11{"bRichRefractiveIndexSector10AndMinus11", 1.03, "barrel RICH refractive index sector 10 and N-11"};
119-
Configurable<float> bRichRefractiveIndexSector11AndNMinus12{"bRichRefractiveIndexSector11AndMinus12", 1.03, "barrel RICH refractive index sector 11 and N-12"};
120-
Configurable<float> bRichRefractiveIndexSector12AndNMinus13{"bRichRefractiveIndexSector12AndMinus13", 1.03, "barrel RICH refractive index sector 12 and N-13"};
121-
Configurable<float> bRichRefractiveIndexSector13AndNMinus14{"bRichRefractiveIndexSector13AndMinus14", 1.03, "barrel RICH refractive index sector 13 and N-14"};
122-
Configurable<float> bRichRefractiveIndexSector14AndNMinus15{"bRichRefractiveIndexSector14AndMinus15", 1.03, "barrel RICH refractive index sector 14 and N-15"};
123-
Configurable<float> bRichRefractiveIndexSector15AndNMinus16{"bRichRefractiveIndexSector15AndMinus16", 1.03, "barrel RICH refractive index sector 15 and N-16"};
124-
Configurable<float> bRichRefractiveIndexSector16AndNMinus17{"bRichRefractiveIndexSector16AndMinus17", 1.03, "barrel RICH refractive index sector 16 and N-17"};
125-
Configurable<float> bRichRefractiveIndexSector17AndNMinus18{"bRichRefractiveIndexSector17AndMinus18", 1.03, "barrel RICH refractive index sector 17 and N-18"};
126-
Configurable<float> bRichRefractiveIndexSector18AndNMinus19{"bRichRefractiveIndexSector18AndMinus19", 1.03, "barrel RICH refractive index sector 18 and N-19"};
127-
Configurable<float> bRichRefractiveIndexSector19AndNMinus20{"bRichRefractiveIndexSector19AndMinus20", 1.03, "barrel RICH refractive index sector 19 and N-20"};
128-
Configurable<float> bRichRefractiveIndexSector20AndNMinus21{"bRichRefractiveIndexSector20AndMinus21", 1.03, "barrel RICH refractive index sector 20 and N-21"};*/
129108
Configurable<float> bRichRefractiveIndexSector0{"bRichRefractiveIndexSector0", 1.03, "barrel RICH refractive index central(s)"}; // central(s)
130109
Configurable<float> bRichRefractiveIndexSector1{"bRichRefractiveIndexSector1", 1.03, "barrel RICH refractive index central(s)-1 and central(s)+1"}; // central(s)-1 and central(s)+1
131110
Configurable<float> bRichRefractiveIndexSector2{"bRichRefractiveIndexSector2", 1.03, "barrel RICH refractive index central(s)-2 and central(s)+2"}; // central(s)-2 and central(s)+2
@@ -543,21 +522,21 @@ struct OnTheFlyRichPid {
543522
float angularResolution(float eta)
544523
{
545524
// Vectors for sampling (USE ANALYTICAL EXTRAPOLATION FOR BETTER RESULTS)
546-
static constexpr float etaSampling[] = {-2.000000, -1.909740, -1.731184, -1.552999, -1.375325, -1.198342, -1.022276, -0.847390, -0.673976, -0.502324, -0.332683, -0.165221, 0.000000, 0.165221, 0.332683, 0.502324, 0.673976, 0.847390, 1.022276, 1.198342, 1.375325, 1.552999, 1.731184, 1.909740, 2.000000};
547-
static constexpr float resRingSamplingWithAbsWalls[] = {0.0009165, 0.000977, 0.001098, 0.001198, 0.001301, 0.001370, 0.001465, 0.001492, 0.001498, 0.001480, 0.001406, 0.001315, 0.001241, 0.001325, 0.001424, 0.001474, 0.001480, 0.001487, 0.001484, 0.001404, 0.001273, 0.001197, 0.001062, 0.000965, 0.0009165};
548-
static constexpr float resRingSamplingWithoutAbsWalls[] = {0.0009165, 0.000977, 0.001095, 0.001198, 0.001300, 0.001369, 0.001468, 0.001523, 0.001501, 0.001426, 0.001299, 0.001167, 0.001092, 0.001179, 0.001308, 0.001407, 0.001491, 0.001508, 0.001488, 0.001404, 0.001273, 0.001196, 0.001061, 0.000965, 0.0009165};
549-
static constexpr int sizeResVector = sizeof(etaSampling) / sizeof(etaSampling[0]);
525+
static constexpr float kEtaSampling[] = {-2.000000, -1.909740, -1.731184, -1.552999, -1.375325, -1.198342, -1.022276, -0.847390, -0.673976, -0.502324, -0.332683, -0.165221, 0.000000, 0.165221, 0.332683, 0.502324, 0.673976, 0.847390, 1.022276, 1.198342, 1.375325, 1.552999, 1.731184, 1.909740, 2.000000};
526+
static constexpr float kResRingSamplingWithAbsWalls[] = {0.0009165, 0.000977, 0.001098, 0.001198, 0.001301, 0.001370, 0.001465, 0.001492, 0.001498, 0.001480, 0.001406, 0.001315, 0.001241, 0.001325, 0.001424, 0.001474, 0.001480, 0.001487, 0.001484, 0.001404, 0.001273, 0.001197, 0.001062, 0.000965, 0.0009165};
527+
static constexpr float kResRingSamplingWithoutAbsWalls[] = {0.0009165, 0.000977, 0.001095, 0.001198, 0.001300, 0.001369, 0.001468, 0.001523, 0.001501, 0.001426, 0.001299, 0.001167, 0.001092, 0.001179, 0.001308, 0.001407, 0.001491, 0.001508, 0.001488, 0.001404, 0.001273, 0.001196, 0.001061, 0.000965, 0.0009165};
528+
static constexpr int sizeResVector = sizeof(kEtaSampling) / sizeof(kEtaSampling[0]);
550529
// Use binary search to find the lower and upper indices
551-
const int lowerIndex = std::lower_bound(etaSampling, etaSampling + sizeResVector, eta) - etaSampling - 1;
530+
const int lowerIndex = std::lower_bound(kEtaSampling, kEtaSampling + sizeResVector, eta) - kEtaSampling - 1;
552531
const int upperIndex = lowerIndex + 1;
553532
if (lowerIndex >= 0 && upperIndex < sizeResVector) {
554533
// Resolution interpolation
555534
if (bRichFlagAbsorbingWalls) {
556-
float interpolatedResRing = interpolate(eta, etaSampling[lowerIndex], etaSampling[upperIndex], resRingSamplingWithAbsWalls[lowerIndex], resRingSamplingWithAbsWalls[upperIndex]);
535+
float interpolatedResRing = interpolate(eta, kEtaSampling[lowerIndex], kEtaSampling[upperIndex], kResRingSamplingWithAbsWalls[lowerIndex], kResRingSamplingWithAbsWalls[upperIndex]);
557536
// std::cout << "Interpolated y value: " << interpolatedY << std::endl;
558537
return interpolatedResRing;
559538
} else {
560-
float interpolatedResRing = interpolate(eta, etaSampling[lowerIndex], etaSampling[upperIndex], resRingSamplingWithoutAbsWalls[lowerIndex], resRingSamplingWithoutAbsWalls[upperIndex]);
539+
float interpolatedResRing = interpolate(eta, kEtaSampling[lowerIndex], kEtaSampling[upperIndex], kResRingSamplingWithoutAbsWalls[lowerIndex], kResRingSamplingWithoutAbsWalls[upperIndex]);
561540
// std::cout << "Interpolated y value: " << interpolatedY << std::endl;
562541
return interpolatedResRing;
563542
}
@@ -630,16 +609,16 @@ struct OnTheFlyRichPid {
630609
const float cosPhi = std::cos(phiC);
631610
// const float ze = thicknessRad / (2. * zP);
632611
const float aZ = zP * cosThetaCherenkov - xP * sinThetaCherenkov * cosPhi;
633-
const float e3Z = std::sqrt(aZ * aZ + (nGas / n) * (nGas / n) - 1.);
612+
const float e3z = std::sqrt(aZ * aZ + (nGas / n) * (nGas / n) - 1.);
634613
const float Z = thicknessGas;
635-
const float alpha = e3Z / aZ;
636-
const float etac = e3Z * n * n;
614+
const float alpha = e3z / aZ;
615+
const float etac = e3z * n * n;
637616
const float k = thicknessRad / (2. * Z);
638617
const float m = 1. / (n * n);
639618
const float lambda = (1. + k * alpha * alpha * alpha) / (1. + k * alpha);
640619
// Derivative d(thetaCherenkov)/dx
641620
const float temp1 = etac / Z;
642-
const float temp2 = alpha * e3Z * cosPhi;
621+
const float temp2 = alpha * e3z * cosPhi;
643622
const float temp3 = xP * sinThetaCherenkov * sinPhi * sinPhi;
644623
const float dThetaX = temp1 * (temp2 - temp3);
645624
// Derivative d(thetaCherenkov)/dy
@@ -835,17 +814,17 @@ struct OnTheFlyRichPid {
835814
float nSigmaBarrelRich[kNspecies] = {kErrorValue, kErrorValue, kErrorValue, kErrorValue, kErrorValue};
836815
bool signalBarrelRich[kNspecies] = {false, false, false, false, false};
837816
float deltaThetaBarrelRich[kNspecies]; //, nSigmaBarrelRich[kNspecies];
838-
static constexpr int lpdgArray[kNspecies] = {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton};
839-
static constexpr float masses[kNspecies] = {o2::track::pid_constants::sMasses[o2::track::PID::Electron],
840-
o2::track::pid_constants::sMasses[o2::track::PID::Muon],
841-
o2::track::pid_constants::sMasses[o2::track::PID::Pion],
842-
o2::track::pid_constants::sMasses[o2::track::PID::Kaon],
843-
o2::track::pid_constants::sMasses[o2::track::PID::Proton]};
817+
static constexpr int kPdgArray[kNspecies] = {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton};
818+
static constexpr float kMasses[kNspecies] = {o2::track::pid_constants::sMasses[o2::track::PID::Electron],
819+
o2::track::pid_constants::sMasses[o2::track::PID::Muon],
820+
o2::track::pid_constants::sMasses[o2::track::PID::Pion],
821+
o2::track::pid_constants::sMasses[o2::track::PID::Kaon],
822+
o2::track::pid_constants::sMasses[o2::track::PID::Proton]};
844823

845824
for (int ii = 0; ii < kNspecies; ii++) { // Loop on the particle hypotheses
846825

847826
float hypothesisAngleBarrelRich = kErrorValue;
848-
const bool hypothesisAngleBarrelRichOk = cherenkovAngle(recoTrack.getP(), masses[ii], aerogelRindex[iSecor], hypothesisAngleBarrelRich);
827+
const bool hypothesisAngleBarrelRichOk = cherenkovAngle(recoTrack.getP(), kMasses[ii], aerogelRindex[iSecor], hypothesisAngleBarrelRich);
849828
signalBarrelRich[ii] = hypothesisAngleBarrelRichOk; // Particle is above the threshold and enough photons
850829

851830
// Evaluate total sigma (layer + tracking resolution)
@@ -855,48 +834,48 @@ struct OnTheFlyRichPid {
855834
double ptResolution = transverseMomentum * transverseMomentum * std::sqrt(recoTrack.getSigma1Pt2());
856835
double etaResolution = std::fabs(std::sin(2.0 * std::atan(std::exp(-recoTrack.getEta())))) * std::sqrt(recoTrack.getSigmaTgl2());
857836
if (flagRICHLoadDelphesLUTs) {
858-
ptResolution = mSmearer.getAbsPtRes(lpdgArray[ii], dNdEta, recoTrack.getEta(), transverseMomentum);
859-
etaResolution = mSmearer.getAbsEtaRes(lpdgArray[ii], dNdEta, recoTrack.getEta(), transverseMomentum);
837+
ptResolution = mSmearer.getAbsPtRes(kPdgArray[ii], dNdEta, recoTrack.getEta(), transverseMomentum);
838+
etaResolution = mSmearer.getAbsEtaRes(kPdgArray[ii], dNdEta, recoTrack.getEta(), transverseMomentum);
860839
}
861840
// cout << endl << "Pt resolution: " << ptResolution << ", Eta resolution: " << etaResolution << endl << endl;
862-
const float barrelTrackAngularReso = calculateTrackAngularResolutionAdvanced(recoTrack.getP() / std::cosh(recoTrack.getEta()), recoTrack.getEta(), ptResolution, etaResolution, masses[ii], aerogelRindex[iSecor]);
841+
const float barrelTrackAngularReso = calculateTrackAngularResolutionAdvanced(recoTrack.getP() / std::cosh(recoTrack.getEta()), recoTrack.getEta(), ptResolution, etaResolution, kMasses[ii], aerogelRindex[iSecor]);
863842
barrelTotalAngularReso = std::hypot(barrelRICHAngularResolution, barrelTrackAngularReso);
864843
if (doQAplots &&
865844
hypothesisAngleBarrelRich > kErrorValue + 1. &&
866845
measuredAngleBarrelRich > kErrorValue + 1. &&
867846
barrelRICHAngularResolution > kErrorValue + 1. &&
868847
flagReachesRadiator) {
869848
switch (mcParticle.pdgCode()) {
870-
case lpdgArray[kEl]: // Electron
871-
case -lpdgArray[kEl]: // Positron
849+
case kPdgArray[kEl]: // Electron
850+
case -kPdgArray[kEl]: // Positron
872851
if (ii == kEl) {
873852
histos.fill(HIST("h2dBarrelAngularResTrackElecVsP"), recoTrack.getP(), 1000.0 * barrelTrackAngularReso);
874853
histos.fill(HIST("h2dBarrelAngularResTotalElecVsP"), recoTrack.getP(), 1000.0 * barrelTotalAngularReso);
875854
}
876855
break;
877-
case lpdgArray[kMu]: // Muon
878-
case -lpdgArray[kMu]: // AntiMuon
856+
case kPdgArray[kMu]: // Muon
857+
case -kPdgArray[kMu]: // AntiMuon
879858
if (ii == kMu) {
880859
histos.fill(HIST("h2dBarrelAngularResTrackMuonVsP"), recoTrack.getP(), 1000.0 * barrelTrackAngularReso);
881860
histos.fill(HIST("h2dBarrelAngularResTotalMuonVsP"), recoTrack.getP(), 1000.0 * barrelTotalAngularReso);
882861
}
883862
break;
884-
case lpdgArray[kPi]: // Pion
885-
case -lpdgArray[kPi]: // AntiPion
863+
case kPdgArray[kPi]: // Pion
864+
case -kPdgArray[kPi]: // AntiPion
886865
if (ii == kPi) {
887866
histos.fill(HIST("h2dBarrelAngularResTrackPionVsP"), recoTrack.getP(), 1000.0 * barrelTrackAngularReso);
888867
histos.fill(HIST("h2dBarrelAngularResTotalPionVsP"), recoTrack.getP(), 1000.0 * barrelTotalAngularReso);
889868
}
890869
break;
891-
case lpdgArray[kKa]: // Kaon
892-
case -lpdgArray[kKa]: // AntiKaon
870+
case kPdgArray[kKa]: // Kaon
871+
case -kPdgArray[kKa]: // AntiKaon
893872
if (ii == kKa) {
894873
histos.fill(HIST("h2dBarrelAngularResTrackKaonVsP"), recoTrack.getP(), 1000.0 * barrelTrackAngularReso);
895874
histos.fill(HIST("h2dBarrelAngularResTotalKaonVsP"), recoTrack.getP(), 1000.0 * barrelTotalAngularReso);
896875
}
897876
break;
898-
case lpdgArray[kPr]: // Proton
899-
case -lpdgArray[kPr]: // AntiProton
877+
case kPdgArray[kPr]: // Proton
878+
case -kPdgArray[kPr]: // AntiProton
900879
if (ii == kPr) {
901880
histos.fill(HIST("h2dBarrelAngularResTrackProtVsP"), recoTrack.getP(), 1000.0 * barrelTrackAngularReso);
902881
histos.fill(HIST("h2dBarrelAngularResTotalProtVsP"), recoTrack.getP(), 1000.0 * barrelTotalAngularReso);
@@ -932,40 +911,40 @@ struct OnTheFlyRichPid {
932911
histos.fill(HIST("hSectorID"), iSecor);
933912

934913
switch (mcParticle.pdgCode()) {
935-
case lpdgArray[kEl]: // Electron
936-
case -lpdgArray[kEl]: // Positron
914+
case kPdgArray[kEl]: // Electron
915+
case -kPdgArray[kEl]: // Positron
937916
histos.fill(HIST("h2dBarrelNsigmaTrueElecVsElecHypothesis"), recoTrack.getP(), nSigmaBarrelRich[0]);
938917
histos.fill(HIST("h2dBarrelNsigmaTrueElecVsMuonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[1]);
939918
histos.fill(HIST("h2dBarrelNsigmaTrueElecVsPionHypothesis"), recoTrack.getP(), nSigmaBarrelRich[2]);
940919
histos.fill(HIST("h2dBarrelNsigmaTrueElecVsKaonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[3]);
941920
histos.fill(HIST("h2dBarrelNsigmaTrueElecVsProtHypothesis"), recoTrack.getP(), nSigmaBarrelRich[4]);
942921
break;
943-
case lpdgArray[kMu]: // Muon
944-
case -lpdgArray[kMu]: // AntiMuon
922+
case kPdgArray[kMu]: // Muon
923+
case -kPdgArray[kMu]: // AntiMuon
945924
histos.fill(HIST("h2dBarrelNsigmaTrueMuonVsElecHypothesis"), recoTrack.getP(), nSigmaBarrelRich[0]);
946925
histos.fill(HIST("h2dBarrelNsigmaTrueMuonVsMuonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[1]);
947926
histos.fill(HIST("h2dBarrelNsigmaTrueMuonVsPionHypothesis"), recoTrack.getP(), nSigmaBarrelRich[2]);
948927
histos.fill(HIST("h2dBarrelNsigmaTrueMuonVsKaonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[3]);
949928
histos.fill(HIST("h2dBarrelNsigmaTrueMuonVsProtHypothesis"), recoTrack.getP(), nSigmaBarrelRich[4]);
950929
break;
951-
case lpdgArray[kPi]: // Pion
952-
case -lpdgArray[kPi]: // AntiPion
930+
case kPdgArray[kPi]: // Pion
931+
case -kPdgArray[kPi]: // AntiPion
953932
histos.fill(HIST("h2dBarrelNsigmaTruePionVsElecHypothesis"), recoTrack.getP(), nSigmaBarrelRich[0]);
954933
histos.fill(HIST("h2dBarrelNsigmaTruePionVsMuonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[1]);
955934
histos.fill(HIST("h2dBarrelNsigmaTruePionVsPionHypothesis"), recoTrack.getP(), nSigmaBarrelRich[2]);
956935
histos.fill(HIST("h2dBarrelNsigmaTruePionVsKaonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[3]);
957936
histos.fill(HIST("h2dBarrelNsigmaTruePionVsProtHypothesis"), recoTrack.getP(), nSigmaBarrelRich[4]);
958937
break;
959-
case lpdgArray[kka]: // Kaon
960-
case -lpdgArray[kka]: // AntiKaon
938+
case kPdgArray[kka]: // Kaon
939+
case -kPdgArray[kka]: // AntiKaon
961940
histos.fill(HIST("h2dBarrelNsigmaTrueKaonVsElecHypothesis"), recoTrack.getP(), nSigmaBarrelRich[0]);
962941
histos.fill(HIST("h2dBarrelNsigmaTrueKaonVsMuonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[1]);
963942
histos.fill(HIST("h2dBarrelNsigmaTrueKaonVsPionHypothesis"), recoTrack.getP(), nSigmaBarrelRich[2]);
964943
histos.fill(HIST("h2dBarrelNsigmaTrueKaonVsKaonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[3]);
965944
histos.fill(HIST("h2dBarrelNsigmaTrueKaonVsProtHypothesis"), recoTrack.getP(), nSigmaBarrelRich[4]);
966945
break;
967-
case lpdgArray[kPr]: // Proton
968-
case -lpdgArray[kPr]: // AntiProton
946+
case kPdgArray[kPr]: // Proton
947+
case -kPdgArray[kPr]: // AntiProton
969948
histos.fill(HIST("h2dBarrelNsigmaTrueProtVsElecHypothesis"), recoTrack.getP(), nSigmaBarrelRich[0]);
970949
histos.fill(HIST("h2dBarrelNsigmaTrueProtVsMuonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[1]);
971950
histos.fill(HIST("h2dBarrelNsigmaTrueProtVsPionHypothesis"), recoTrack.getP(), nSigmaBarrelRich[2]);

0 commit comments

Comments
 (0)