From a5f028ffc5f1547f17093ce162fdcdc7542eee74 Mon Sep 17 00:00:00 2001
From: Pavel Batyuk <pavel.batyuk@jinr.ru>
Date: Mon, 5 Sep 2016 14:04:48 +0300
Subject: [PATCH] PID TPC:: recalculated probability coefficients have been
added (from G. Eyyubova, SINP, MSI).
---
tpc/MpdParticleIdentification.cxx | 355 ++++++++++++++++++++++++++++++
1 file changed, 355 insertions(+)
diff --git a/tpc/MpdParticleIdentification.cxx b/tpc/MpdParticleIdentification.cxx
index 554a64a..74b0836 100644
--- a/tpc/MpdParticleIdentification.cxx
+++ b/tpc/MpdParticleIdentification.cxx
@@ -32,6 +32,361 @@ MpdParticleIdentification::MpdParticleIdentification() {
MpdParticleIdentification::~MpdParticleIdentification() {
}
+//Int_t MpdParticleIdentification::GetTpcProbs(Float_t P, Float_t dedx, Int_t nHits, Float_t& Ppi, Float_t& Pk, Float_t& Pp, Float_t& Pe, Int_t method) {
+//
+// /*Parameters for fitting*/
+// //AZ Float_t *ProtonPar, *PionPar, *KaonPar, *ElectronPar;
+// /************************/
+//
+// dedx *= 1000000; // converting to keV/cm
+// const Int_t Ntypes = 4; //order: pion, kaon, proton, electron
+// Float_t ProtonPar[9], PionPar[9], KaonPar[9], ElectronPar[9], resultProb[Ntypes]; //AZ
+// const Int_t Nintervals = 10;
+// Float_t sigmas[Ntypes][Nintervals] = {//array of sigmas for different particles and different intervals of momentum. They were got from 100000 BOX events
+// {0.113, 0.107, 0.100, 0.107, 0.097, 0.111, 0.122, 0.105, 0.115, 0.118},
+// {0.407, 0.276, 0.260, 0.159, 0.185, 0.154, 0.122, 0.105, 0.115, 0.118},
+// {0.507, 0.488, 0.394, 0.330, 0.268, 0.244, 0.260, 0.172, 0.193, 0.118},
+// {0.163, 0.160, 0.149, 0.159, 0.172, 0.205, 0.260, 0.172, 0.193, 0.156}
+// };
+//
+// Float_t sigma = 1.0 / TMath::Sqrt(nHits); //for gaussian
+// Float_t sum = 0.0; // for normalizing
+//
+// /*A priori coefficients for Bayes approach */
+// Float_t bayesAprioriCoefficients[Ntypes];
+// /*A "measured" probabilities */
+// Float_t gausProb[Ntypes];
+//
+// Float_t sigPi, sigPr, sigKa, sigEl; //sigmas for gaussians
+//
+// switch (method) {
+//
+// case 0: //bethe-bloch approximation with "standard" sigmas and equal bayesian coefficients
+//
+// //parameters were got from Bethe-Bloch approximation for 100000 BOX events
+// /*AZ
+// ProtonPar = new Float_t[4];
+// PionPar = new Float_t[4];
+// KaonPar = new Float_t[4];
+// ElectronPar = new Float_t[4];
+// */
+// ProtonPar[0] = -3.957;
+// ProtonPar[1] = 3.525;
+// ProtonPar[2] = 16.468;
+// ProtonPar[3] = 0.815;
+// ProtonPar[4] = 5.247;
+// PionPar[0] = -0.739;
+// PionPar[1] = 7.373;
+// PionPar[2] = 3904.790;
+// PionPar[3] = 0.274;
+// PionPar[4] = 5.497;
+// KaonPar[0] = -2.590;
+// KaonPar[1] = 4.918;
+// KaonPar[2] = 79.722;
+// KaonPar[3] = 0.357;
+// KaonPar[4] = 4.511;
+// ElectronPar[0] = -1.552;
+// ElectronPar[1] = 1.748;
+// ElectronPar[2] = 7.425;
+// ElectronPar[3] = 0.980;
+// ElectronPar[4] = 1.604;
+//
+// sigma = 0.07 * dedx;
+// sigPi = sigPr = sigKa = sigEl = sigma;
+// bayesAprioriCoefficients[0] = 1.0;
+// bayesAprioriCoefficients[1] = 1.0;
+// bayesAprioriCoefficients[2] = 1.0;
+// bayesAprioriCoefficients[3] = 1.0;
+//
+// gausProb[0] = Gaus(dedx, BetheBlochFunction(P, PionPar), sigPi, kTRUE);
+// gausProb[1] = Gaus(dedx, BetheBlochFunction(P, KaonPar), sigKa, kTRUE);
+// gausProb[2] = Gaus(dedx, BetheBlochFunction(P, ProtonPar), sigPr, kTRUE);
+// gausProb[3] = Gaus(dedx, BetheBlochFunction(P, ElectronPar), sigEl, kTRUE);
+//
+// sum = gausProb[0] + gausProb[1] + gausProb[2] + gausProb[3];
+// if (sum == 0.) return 1;
+// gausProb[0] /= sum;
+// gausProb[1] /= sum;
+// gausProb[2] /= sum;
+// gausProb[3] /= sum;
+// break;
+//
+// case 1: //bethe-bloch approximation with special different sigmas and byesian coefficients
+//
+// //parameters were got from Bethe-Bloch approximation for 100000 BOX events
+// /*AZ
+// ProtonPar = new Float_t[4];
+// PionPar = new Float_t[4];
+// KaonPar = new Float_t[4];
+// ElectronPar = new Float_t[4];
+// */
+// ProtonPar[0] = -3.957;
+// ProtonPar[1] = 3.525;
+// ProtonPar[2] = 16.468;
+// ProtonPar[3] = 0.815;
+// ProtonPar[4] = 5.247;
+// PionPar[0] = -0.739;
+// PionPar[1] = 7.373;
+// PionPar[2] = 3904.790;
+// PionPar[3] = 0.274;
+// PionPar[4] = 5.497;
+// KaonPar[0] = -2.590;
+// KaonPar[1] = 4.918;
+// KaonPar[2] = 79.722;
+// KaonPar[3] = 0.357;
+// KaonPar[4] = 4.511;
+// ElectronPar[0] = -1.552;
+// ElectronPar[1] = 1.748;
+// ElectronPar[2] = 7.425;
+// ElectronPar[3] = 0.980;
+// ElectronPar[4] = 1.604;
+//
+// if (P < 0.3) {
+// bayesAprioriCoefficients[0] = 0.28;
+// bayesAprioriCoefficients[1] = 0.25;
+// bayesAprioriCoefficients[2] = 0.26;
+// bayesAprioriCoefficients[3] = 0.21;
+// sigPi = sigmas[0][0];
+// sigKa = sigmas[1][0];
+// sigPr = sigmas[2][0];
+// sigEl = sigmas[3][0];
+// } else if ((P >= 0.3) && (P < 0.4)) {
+// bayesAprioriCoefficients[0] = 0.91;
+// bayesAprioriCoefficients[1] = 0.02;
+// bayesAprioriCoefficients[2] = 0.06;
+// bayesAprioriCoefficients[3] = 0.01;
+// sigPi = sigmas[0][1];
+// sigKa = sigmas[1][1];
+// sigPr = sigmas[2][1];
+// sigEl = sigmas[3][1];
+// } else if ((P >= 0.4) && (P < 0.5)) {
+// bayesAprioriCoefficients[0] = 0.70;
+// bayesAprioriCoefficients[1] = 0.11;
+// bayesAprioriCoefficients[2] = 0.13;
+// bayesAprioriCoefficients[3] = 0.06;
+// sigPi = sigmas[0][2];
+// sigKa = sigmas[1][2];
+// sigPr = sigmas[2][2];
+// sigEl = sigmas[3][2];
+// } else if ((P >= 0.5) && (P < 0.6)) {
+// bayesAprioriCoefficients[0] = 0.39;
+// bayesAprioriCoefficients[1] = 0.19;
+// bayesAprioriCoefficients[2] = 0.32;
+// bayesAprioriCoefficients[3] = 0.10;
+// sigPi = sigmas[0][3];
+// sigKa = sigmas[1][3];
+// sigPr = sigmas[2][3];
+// sigEl = sigmas[3][3];
+// } else if ((P >= 0.6) && (P < 0.7)) {
+// bayesAprioriCoefficients[0] = 0.41;
+// bayesAprioriCoefficients[1] = 0.17;
+// bayesAprioriCoefficients[2] = 0.37;
+// bayesAprioriCoefficients[3] = 0.5;
+// sigPi = sigmas[0][4];
+// sigKa = sigmas[1][4];
+// sigPr = sigmas[2][4];
+// sigEl = sigmas[3][4];
+// } else if ((P >= 0.7) && (P < 0.8)) {
+// bayesAprioriCoefficients[0] = 0.43;
+// bayesAprioriCoefficients[1] = 0.16;
+// bayesAprioriCoefficients[2] = 0.31;
+// bayesAprioriCoefficients[3] = 0.10;
+// sigPi = sigmas[0][5];
+// sigKa = sigmas[1][5];
+// sigPr = sigmas[2][5];
+// sigEl = sigmas[3][5];
+// } else if ((P >= 0.8) && (P < 0.9)) {
+// bayesAprioriCoefficients[0] = 0.44;
+// bayesAprioriCoefficients[1] = 0.11;
+// bayesAprioriCoefficients[2] = 0.43;
+// bayesAprioriCoefficients[3] = 0.02;
+// sigPi = sigmas[0][6];
+// sigKa = sigmas[1][6];
+// sigPr = sigmas[2][6];
+// sigEl = sigmas[3][6];
+// } else if ((P >= 0.9) && (P < 1.0)) {
+// bayesAprioriCoefficients[0] = 0.36;
+// bayesAprioriCoefficients[1] = 0.21;
+// bayesAprioriCoefficients[2] = 0.36;
+// bayesAprioriCoefficients[3] = 0.07;
+// sigPi = sigmas[0][7];
+// sigKa = sigmas[1][7];
+// sigPr = sigmas[2][7];
+// sigEl = sigmas[3][7];
+// } else if ((P >= 1.0) && (P < 1.2)) {
+// bayesAprioriCoefficients[0] = 0.32;
+// bayesAprioriCoefficients[1] = 0.32;
+// bayesAprioriCoefficients[2] = 0.32;
+// bayesAprioriCoefficients[3] = 0.04;
+// sigPi = sigmas[0][8];
+// sigKa = sigmas[1][8];
+// sigPr = sigmas[2][8];
+// sigEl = sigmas[3][8];
+// } else if (P >= 1.2) {
+// bayesAprioriCoefficients[0] = 0.34;
+// bayesAprioriCoefficients[1] = 0.27;
+// bayesAprioriCoefficients[2] = 0.31;
+// bayesAprioriCoefficients[3] = 0.08;
+// sigPi = sigmas[0][9];
+// sigKa = sigmas[1][9];
+// sigPr = sigmas[2][9];
+// sigEl = sigmas[3][9];
+// }
+//
+// gausProb[0] = Gaus(dedx, BetheBlochFunction(P, PionPar), sigPi, kTRUE);
+// gausProb[1] = Gaus(dedx, BetheBlochFunction(P, KaonPar), sigKa, kTRUE);
+// gausProb[2] = Gaus(dedx, BetheBlochFunction(P, ProtonPar), sigPr, kTRUE);
+// gausProb[3] = Gaus(dedx, BetheBlochFunction(P, ElectronPar), sigEl, kTRUE);
+//
+// sum = gausProb[0] + gausProb[1] + gausProb[2] + gausProb[3];
+// if (sum == 0.) return 1;
+// gausProb[0] /= sum;
+// gausProb[1] /= sum;
+// gausProb[2] /= sum;
+// gausProb[3] /= sum;
+// break;
+//
+// case 2: //bethe-bloch approximation with "standard" sigmas and different byesian coefficients
+//
+// //parameters were got from Bethe-Bloch approximation for 100000 BOX events
+// /*AZ
+// ProtonPar = new Float_t[4];
+// PionPar = new Float_t[4];
+// KaonPar = new Float_t[4];
+// ElectronPar = new Float_t[4];
+// */
+// ProtonPar[0] = -3.957;
+// ProtonPar[1] = 3.525;
+// ProtonPar[2] = 16.468;
+// ProtonPar[3] = 0.815;
+// ProtonPar[4] = 5.247;
+// PionPar[0] = -0.739;
+// PionPar[1] = 7.373;
+// PionPar[2] = 3904.790;
+// PionPar[3] = 0.274;
+// PionPar[4] = 5.497;
+// KaonPar[0] = -2.590;
+// KaonPar[1] = 4.918;
+// KaonPar[2] = 79.722;
+// KaonPar[3] = 0.357;
+// KaonPar[4] = 4.511;
+// ElectronPar[0] = -1.552;
+// ElectronPar[1] = 1.748;
+// ElectronPar[2] = 7.425;
+// ElectronPar[3] = 0.980;
+// ElectronPar[4] = 1.604;
+//
+// sigma = 0.07 * dedx;
+// sigPi = sigPr = sigKa = sigEl = sigma;
+//
+// if (P < 0.3) {
+// bayesAprioriCoefficients[0] = 0.28;
+// bayesAprioriCoefficients[1] = 0.25;
+// bayesAprioriCoefficients[2] = 0.26;
+// bayesAprioriCoefficients[3] = 0.21;
+// } else if ((P >= 0.3) && (P < 0.4)) {
+// bayesAprioriCoefficients[0] = 0.91;
+// bayesAprioriCoefficients[1] = 0.02;
+// bayesAprioriCoefficients[2] = 0.06;
+// bayesAprioriCoefficients[3] = 0.01;
+// } else if ((P >= 0.4) && (P < 0.5)) {
+// bayesAprioriCoefficients[0] = 0.70;
+// bayesAprioriCoefficients[1] = 0.11;
+// bayesAprioriCoefficients[2] = 0.13;
+// bayesAprioriCoefficients[3] = 0.06;
+// } else if ((P >= 0.5) && (P < 0.6)) {
+// bayesAprioriCoefficients[0] = 0.39;
+// bayesAprioriCoefficients[1] = 0.19;
+// bayesAprioriCoefficients[2] = 0.32;
+// bayesAprioriCoefficients[3] = 0.10;
+// } else if ((P >= 0.6) && (P < 0.7)) {
+// bayesAprioriCoefficients[0] = 0.41;
+// bayesAprioriCoefficients[1] = 0.17;
+// bayesAprioriCoefficients[2] = 0.37;
+// bayesAprioriCoefficients[3] = 0.5;
+// } else if ((P >= 0.7) && (P < 0.8)) {
+// bayesAprioriCoefficients[0] = 0.43;
+// bayesAprioriCoefficients[1] = 0.16;
+// bayesAprioriCoefficients[2] = 0.31;
+// bayesAprioriCoefficients[3] = 0.10;
+// } else if ((P >= 0.8) && (P < 0.9)) {
+// bayesAprioriCoefficients[0] = 0.44;
+// bayesAprioriCoefficients[1] = 0.11;
+// bayesAprioriCoefficients[2] = 0.43;
+// bayesAprioriCoefficients[3] = 0.02;
+// } else if ((P >= 0.9) && (P < 1.0)) {
+// bayesAprioriCoefficients[0] = 0.36;
+// bayesAprioriCoefficients[1] = 0.21;
+// bayesAprioriCoefficients[2] = 0.36;
+// bayesAprioriCoefficients[3] = 0.07;
+// } else if ((P >= 1.0) && (P < 1.2)) {
+// bayesAprioriCoefficients[0] = 0.32;
+// bayesAprioriCoefficients[1] = 0.32;
+// bayesAprioriCoefficients[2] = 0.32;
+// bayesAprioriCoefficients[3] = 0.04;
+// } else if (P >= 1.2) {
+// bayesAprioriCoefficients[0] = 0.34;
+// bayesAprioriCoefficients[1] = 0.27;
+// bayesAprioriCoefficients[2] = 0.31;
+// bayesAprioriCoefficients[3] = 0.08;
+// }
+//
+// gausProb[0] = Gaus(dedx, BetheBlochFunction(P, PionPar), sigPi, kTRUE);
+// gausProb[1] = Gaus(dedx, BetheBlochFunction(P, KaonPar), sigKa, kTRUE);
+// gausProb[2] = Gaus(dedx, BetheBlochFunction(P, ProtonPar), sigPr, kTRUE);
+// gausProb[3] = Gaus(dedx, BetheBlochFunction(P, ElectronPar), sigEl, kTRUE);
+//
+// sum = gausProb[0] + gausProb[1] + gausProb[2] + gausProb[3];
+// if (sum == 0.) return 1;
+// gausProb[0] /= sum;
+// gausProb[1] /= sum;
+// gausProb[2] /= sum;
+// gausProb[3] /= sum;
+// break;
+//
+// case 3: //parabolic approximation
+//
+// //parameters were got from parabolic approximation function for 2000 UrQMD events
+// ProtonPar[0] = 0.031;
+// ProtonPar[1] = -0.124;
+// ProtonPar[2] = 1.412;
+// PionPar[0] = 0.230;
+// PionPar[1] = 0.088;
+// PionPar[2] = 1.072;
+// KaonPar[0] = 0.676;
+// KaonPar[1] = 0.621;
+// KaonPar[2] = 0.831;
+// ElectronPar[0] = 0.000;
+// ElectronPar[1] = -0.018;
+// ElectronPar[2] = 2.055;
+//
+// Float_t invP = 1.0 / P;
+//
+// gausProb[0] = Gaus(dedx, ParabolicFunction(invP, PionPar), sigPi, kTRUE);
+// gausProb[1] = Gaus(dedx, ParabolicFunction(invP, KaonPar), sigKa, kTRUE);
+// gausProb[2] = Gaus(dedx, ParabolicFunction(invP, ProtonPar), sigPr, kTRUE);
+// gausProb[3] = Gaus(dedx, ParabolicFunction(invP, ElectronPar), sigEl, kTRUE);
+//
+// sum = gausProb[0] + gausProb[1] + gausProb[2] + gausProb[3];
+// if (sum == 0.) return 1;
+// gausProb[0] /= sum;
+// gausProb[1] /= sum;
+// gausProb[2] /= sum;
+// gausProb[3] /= sum;
+//
+// break;
+// }
+// //AZ Float_t *resultProb = new Float_t[Ntypes];
+// BayesFunction(gausProb, bayesAprioriCoefficients, resultProb, Ntypes);
+//
+// Ppi = resultProb[0];
+// Pk = resultProb[1];
+// Pp = resultProb[2];
+// Pe = resultProb[3];
+//
+// return 0;
+//}
+
Int_t MpdParticleIdentification::GetTofProbs(Float_t P, Float_t beta, Float_t& Ppi, Float_t& Pk, Float_t& Pp, Float_t& Pe, Int_t method) {
const Float_t m2_proton = 0.938 * 0.938; //square of proton mass (GeV)
--
GitLab