Skip to content
Snippets Groups Projects
Commit 54c2b1c4 authored by Pavel Batyuk's avatar Pavel Batyuk
Browse files

PID TPC:: recalculated probability coefficients have been added (from G. Eyyubova, SINP, MSI).

parent cbda7250
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
input/MPDTPCPidResponseVhelle.root 0 → 100644
+ 0
− 0
View file @ 54c2b1c4
File added
mpddata/MpdFillDstTask.cxx
+ 238
− 241
View file @ 54c2b1c4
This diff is collapsed.
tpc/CMakeLists.txt
+ 1
− 0
View file @ 54c2b1c4
...@@ -50,6 +50,7 @@ MpdTpcHitProducer.cxx ...@@ -50,6 +50,7 @@ MpdTpcHitProducer.cxx
MpdTpcSectorGeo.cxx MpdTpcSectorGeo.cxx
MpdTpcDigitizerAZ.cxx MpdTpcDigitizerAZ.cxx
MpdTpcClusterFinderAZ.cxx MpdTpcClusterFinderAZ.cxx
MpdTPCpid.cxx
) )
# fill list of header files from list of source files # fill list of header files from list of source files
......
tpc/MpdParticleIdentification.cxx
+ 0
− 355
View file @ 54c2b1c4
...@@ -32,361 +32,6 @@ MpdParticleIdentification::MpdParticleIdentification() { ...@@ -32,361 +32,6 @@ MpdParticleIdentification::MpdParticleIdentification() {
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) { 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) const Float_t m2_proton = 0.938 * 0.938; //square of proton mass (GeV)
......
tpc/MpdParticleIdentification.h
+ 6
− 3
View file @ 54c2b1c4
...@@ -19,15 +19,18 @@ ...@@ -19,15 +19,18 @@
// Base Class Headers ---------------- // Base Class Headers ----------------
#include "TSystem.h" #include "TSystem.h"
#include "MpdTPCpid.h"
class MpdParticleIdentification {
class MpdParticleIdentification : public MpdTPCpid {
public: public:
// Constructors/Destructors --------- // Constructors/Destructors ---------
MpdParticleIdentification(); MpdParticleIdentification();
~MpdParticleIdentification(); virtual ~MpdParticleIdentification();
Int_t 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); // Int_t 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);
// The function is realized now in the MpdTPCpid-class with corrected probability coefficients
Int_t GetTofProbs(Float_t P, Float_t beta, Float_t& Ppi, Float_t& PK, Float_t& Pp, Float_t& Pe, Int_t method); Int_t GetTofProbs(Float_t P, Float_t beta, Float_t& Ppi, Float_t& PK, Float_t& Pp, Float_t& Pe, Int_t method);
Int_t GetCombinedProbs(Float_t *tofProbs, Float_t *tpcProbs, Float_t *resultProbs, Int_t N); Int_t GetCombinedProbs(Float_t *tofProbs, Float_t *tpcProbs, Float_t *resultProbs, Int_t N);
......
tpc/MpdTPCpid.cxx 0 → 100644
+ 285
− 0
View file @ 54c2b1c4
//-----------------------------------------------------------
// File and Version Information:
// $Id$
//
// Description:
// Implementation of class ParticleIdentification
// see ParticleIdentification.h for details
//
// Environment:
// Software developed for MPD at NICA.
//
// Author List:
// Gyulnara Eyyubova
//
//-----------------------------------------------------------
// C/C++ Headers ----------------------
#include <iostream>
#include <TMath.h>
#include <TF1.h>
#include <TH1.h>
#include <TH2.h>
// This Class' Header ------------------
#include "MpdTPCpid.h"
using namespace std;
using namespace TMath;
// Class Member definitions -----------
MpdTPCpid::MpdTPCpid() :
ProtonPar(),
PionPar(),
KaonPar(),
ElectronPar(),
sigmasPi(),
sigmasPr(),
sigmasKa(),
sigmasEl(),
fCoefficients(NULL) {
ReadTPCResponse();
}
MpdTPCpid::~MpdTPCpid() {
fCoefficients->Delete();
}
Int_t MpdTPCpid::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) {
const Int_t Ntypes = 4; //order: pion, kaon, proton, electron
Float_t resultProb[Ntypes];
const Int_t Nintervals = 10;
Double_t P_int[Nintervals + 1] = {0.1, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.2, 3};
//default sigmas
Float_t sigma = 0.07 * dedx;
Float_t sum = 0.0; // for normalizing
/*A priori coefficients for Bayes approach */
Float_t bayesAprioriCoefficients[Ntypes];
/*A "measured" probabilities */
Float_t gausProb[Ntypes];
for (Int_t i = 0; i < Ntypes; i++) gausProb[i] = 0;
Int_t i_p = 0; // momentum interval
switch (method) {
case 0: //bethe-bloch approximation equal bayesian coefficients
if (P > 3) i_p = 9;
for (Int_t k = 1; k < Nintervals; k++) if (P >= P_int[k] && P < P_int[k + 1]) i_p = k;
bayesAprioriCoefficients[0] = 1.0;
bayesAprioriCoefficients[1] = 1.0;
bayesAprioriCoefficients[2] = 1.0;
bayesAprioriCoefficients[3] = 1.0;
gausProb[0] = Gaus(dedx, BetheBlochFunction(P, PionPar), sigmasPi[i_p], kTRUE); // kTRUE = normilized by sqrt(2*Pi)*sigma.
gausProb[1] = Gaus(dedx, BetheBlochFunction(P, KaonPar), sigmasKa[i_p], kTRUE);
gausProb[2] = Gaus(dedx, BetheBlochFunction(P, ProtonPar), sigmasPr[i_p], kTRUE);
gausProb[3] = Gaus(dedx, BetheBlochFunction(P, ElectronPar), sigmasEl[i_p], 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 byesian coefficients
if (P > 3) i_p = 9;
for (Int_t k = 1; k < Nintervals; k++) if (P >= P_int[k] && P < P_int[k + 1]) i_p = k;
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), sigmasPi[i_p], kTRUE);
gausProb[1] = Gaus(dedx, BetheBlochFunction(P, KaonPar), sigmasKa[i_p], kTRUE);
gausProb[2] = Gaus(dedx, BetheBlochFunction(P, ProtonPar), sigmasPr[i_p], kTRUE);
gausProb[3] = Gaus(dedx, BetheBlochFunction(P, ElectronPar), sigmasEl[i_p], 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;
}
BayesFunction(gausProb, bayesAprioriCoefficients, resultProb, Ntypes);
Ppi = resultProb[0];
Pk = resultProb[1];
Pp = resultProb[2];
Pe = resultProb[3];
return 0;
}
Float_t MpdTPCpid::BetheBlochFunction(Float_t x, Float_t *p) {
Float_t b = 1 / (x / Sqrt(1 + x * x));
return p[0] / Power(b, p[3]) * (p[1] - Power(b, p[3]) - Log(p[2] + Power(1 / x, p[4])));
}
/*Bayes function for probabilities calculation*/
Int_t MpdTPCpid::BayesFunction(Float_t *measProb, Float_t *aprioriProb, Float_t *bayesProb, Int_t N) {
/*measProb - measured probabilities from TPC/TOF/etc*/
/*aprioriProb - a priori probabilities from some magic*/
/*bayesProb - result bayes probabilities */
/*N - number of types {pion, kaon, proton, electron} */
Float_t sumBayes = 0.0;
for (Int_t i = 0; i < N; ++i) {
sumBayes += measProb[i] * aprioriProb[i];
}
if (sumBayes == 0.) return 1;
for (Int_t i = 0; i < N; ++i) {
bayesProb[i] = measProb[i] * aprioriProb[i] / sumBayes;
}
return 0;
}
void MpdTPCpid::ReadTPCResponse() {
TString c;
/* Accessing the DeDx response, obtained with VHLEE generator */
fCoefficients = new TFile("$VMCWORKDIR/input/MPDTPCPidResponseVhelle.root");
if (fCoefficients->IsZombie()) {
printf("File MPDTPCPidResponseVhelle.root does not exist.\n");
}
TH2D *hdedx_pi = (TH2D*) fCoefficients->Get("dedx_pi");
TH2D *hdedx_pr = (TH2D*) fCoefficients->Get("dedx_pr");
TH2D *hdedx_ka = (TH2D*) fCoefficients->Get("dedx_ka");
TH2D *hdedx_el = (TH2D*) fCoefficients->Get("dedx_el");
TF1 *fBetheBlochPion = hdedx_pi->GetFunction("BetheBlochALEPH");
PionPar[0] = fBetheBlochPion->GetParameter(0);
PionPar[1] = fBetheBlochPion->GetParameter(1);
PionPar[2] = fBetheBlochPion->GetParameter(2);
PionPar[3] = fBetheBlochPion->GetParameter(3);
PionPar[4] = fBetheBlochPion->GetParameter(4);
TF1 *fBetheBlochProton = hdedx_pr->GetFunction("BetheBlochALEPH");
ProtonPar[0] = fBetheBlochProton->GetParameter(0);
ProtonPar[1] = fBetheBlochProton->GetParameter(1);
ProtonPar[2] = fBetheBlochProton->GetParameter(2);
ProtonPar[3] = fBetheBlochProton->GetParameter(3);
ProtonPar[4] = fBetheBlochProton->GetParameter(4);
TF1 *fBetheBlochKaon = hdedx_ka->GetFunction("BetheBlochALEPH");
KaonPar[0] = fBetheBlochKaon->GetParameter(0);
KaonPar[1] = fBetheBlochKaon->GetParameter(1);
KaonPar[2] = fBetheBlochKaon->GetParameter(2);
KaonPar[3] = fBetheBlochKaon->GetParameter(3);
KaonPar[4] = fBetheBlochKaon->GetParameter(4);
TF1 *fBetheBlochElectron = hdedx_el->GetFunction("BetheBlochALEPH");
ElectronPar[0] = fBetheBlochElectron->GetParameter(0);
ElectronPar[1] = fBetheBlochElectron->GetParameter(1);
ElectronPar[2] = fBetheBlochElectron->GetParameter(2);
ElectronPar[3] = fBetheBlochElectron->GetParameter(3);
ElectronPar[4] = fBetheBlochElectron->GetParameter(4);
const Int_t Nintervals = 10;
Double_t P_int[Nintervals + 1] = {0.1, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.2, 3};
for (Int_t k = 0; k < Nintervals; k++) {
Int_t p_low_bin = hdedx_pi->GetXaxis()->FindBin(P_int[k]);
Int_t p_up_bin = hdedx_pi->GetXaxis()->FindBin(P_int[k + 1]);
c = Form("hpi%i", k);
TH1D *hpi = (TH1D*) hdedx_pi->ProjectionY(c, p_low_bin, p_up_bin);
hpi->Fit("gaus");
sigmasPi[k] = hpi->GetFunction("gaus")->GetParameter(2);
c = Form("hpr%i", k);
TH1D *hpr = (TH1D*) hdedx_pr->ProjectionY("hpr", p_low_bin, p_up_bin);
hpr->Fit("gaus");
sigmasPr[k] = hpr->GetFunction("gaus")->GetParameter(2);
c = Form("hka%i", k);
TH1D *hka = (TH1D*) hdedx_ka->ProjectionY("hka", p_low_bin, p_up_bin);
hka->Fit("gaus");
sigmasKa[k] = hka->GetFunction("gaus")->GetParameter(2);
c = Form("hel%i", k);
TH1D *hel = (TH1D*) hdedx_el->ProjectionY("hel", p_low_bin, p_up_bin);
hel->Fit("gaus");
sigmasEl[k] = hel->GetFunction("gaus")->GetParameter(2);
}
// f->Close();
}
ClassImp(MpdTPCpid);
tpc/MpdTPCpid.h 0 → 100644
+ 56
− 0
View file @ 54c2b1c4
//-----------------------------------------------------------
// File and Version Information:
// $Id$
//
// Description:
// class for particle identification
//
//
// Environment:
// Software developed for MPD at NICA.
//
// Author List:
// Gyulnara Eyyubova
//
//-----------------------------------------------------------
#ifndef MpdTPCpid_HH
#define MpdTPCpid_HH
// Base Class Headers ----------------
#include "TSystem.h"
#include <TGraph.h>
#include <TFile.h>
class MpdTPCpid {
public:
// Constructors/Destructors ---------
MpdTPCpid();
virtual ~MpdTPCpid();
Int_t 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);
Float_t BetheBlochFunction(Float_t x, Float_t *p);
Int_t BayesFunction(Float_t *measProb, Float_t *aprioriProb, Float_t *bayesProb, Int_t N);
void ReadTPCResponse();
private:
Float_t ProtonPar[5];
Float_t PionPar[5];
Float_t KaonPar[5];
Float_t ElectronPar[5];
//const Int_t Nintervals = 10;
Float_t sigmasPi[10];
Float_t sigmasPr[10];
Float_t sigmasKa[10];
Float_t sigmasEl[10];
TFile* fCoefficients;
public:
ClassDef(MpdTPCpid, 1)
};
#endif
tpc/tpcLinkDef.h
+ 1
− 0
View file @ 54c2b1c4
...@@ -27,6 +27,7 @@ ...@@ -27,6 +27,7 @@
#pragma link C++ class MpdTpcSectorGeo+; #pragma link C++ class MpdTpcSectorGeo+;
#pragma link C++ class MpdTpcDigitizerAZ+; #pragma link C++ class MpdTpcDigitizerAZ+;
#pragma link C++ class MpdTpcClusterFinderAZ+; #pragma link C++ class MpdTpcClusterFinderAZ+;
#pragma link C++ class MpdTPCpid+;
#endif #endif
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment