MpdTPCpid.cxx 12 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386
  1. //-----------------------------------------------------------
  2. // File and Version Information:
  3. // $Id$
  4. //
  5. // Description:
  6. // Implementation of class ParticleIdentification
  7. // see ParticleIdentification.h for details
  8. //
  9. // Environment:
  10. // Software developed for MPD at NICA.
  11. //
  12. // Author List:
  13. // Gyulnara Eyyubova
  14. //
  15. //-----------------------------------------------------------
  16. // C/C++ Headers ----------------------
  17. #include <iostream>
  18. #include <fstream>
  19. #include <stdio.h>
  20. #include <TMath.h>
  21. #include <TF1.h>
  22. #include <TH1.h>
  23. #include <TH2.h>
  24. #include <TGraphErrors.h>
  25. // This Class' Header ------------------
  26. #include "MpdTPCpid.h"
  27. using namespace std;
  28. using namespace TMath;
  29. // Class Member definitions -----------
  30. MpdTPCpid::MpdTPCpid() :
  31. fn(1.),
  32. ProtonPar(),
  33. PionPar(),
  34. KaonPar(),
  35. ElectronPar(),
  36. sigmasPi(),
  37. sigmasPr(),
  38. sigmasKa(),
  39. sigmasEl()
  40. {
  41. ReadTPCResponse();
  42. }
  43. MpdTPCpid::~MpdTPCpid() {
  44. }
  45. 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) {
  46. const Int_t Ntypes = 4; //order: pion, kaon, proton, electron
  47. Float_t resultProb[Ntypes];
  48. const Int_t Nintervals = 72;
  49. Double_t P_int[Nintervals + 1] = {0.09,0.12,0.15,0.18,0.21,0.24,0.27,0.3,0.33,0.36,0.39,0.42,0.45,0.48,0.51,0.54,0.57,0.6,0.63,0.66,0.69,0.72,0.75,0.78,0.81,0.84,0.87,0.9,0.93,0.96,0.99,1.02,1.05,1.08,
  50. 1.11,1.14,1.17,1.2,1.23,1.26,1.29,1.32,1.35,1.38,1.41,1.44,1.47,1.5,1.53,1.56,1.59,1.62,1.65,1.68,1.71,1.74,1.77,1.8,1.83,1.86,1.89,1.92,1.95,1.98,2.01,2.04,2.07,2.1,
  51. 2.13,2.16,2.19,2.22, 3};
  52. //default sigmas
  53. //Float_t sigma = 0.07 * dedx;
  54. Float_t sum = 0.0; // for normalizing
  55. /*A priori coefficients for Bayes approach */
  56. Float_t bayesAprioriCoefficients[Ntypes];
  57. /*A "measured" probabilities */
  58. Float_t gausProb[Ntypes];
  59. for (Int_t i = 0; i < Ntypes; i++) gausProb[i] = 0;
  60. Int_t i_p = 0; // momentum interval
  61. switch (method) {
  62. case 0: //bethe-bloch approximation equal bayesian coefficients
  63. if (P > 3) i_p = 71;
  64. for (Int_t k = 0; k < Nintervals; k++) if (P >= P_int[k] && P < P_int[k + 1]) i_p = k;
  65. bayesAprioriCoefficients[0] = 1.0;
  66. bayesAprioriCoefficients[1] = 1.0;
  67. bayesAprioriCoefficients[2] = 1.0;
  68. bayesAprioriCoefficients[3] = 1.0;
  69. gausProb[0] = Gaus(dedx, BetheBlochFunction(P, PionPar), sigmasPi[i_p], kTRUE); // kTRUE = normilized by sqrt(2*Pi)*sigma.
  70. gausProb[1] = Gaus(dedx, BetheBlochFunction(P, KaonPar), sigmasKa[i_p], kTRUE);
  71. gausProb[2] = Gaus(dedx, BetheBlochFunction(P, ProtonPar), sigmasPr[i_p], kTRUE);
  72. gausProb[3] = Gaus(dedx, BetheBlochFunction(P, ElectronPar), sigmasEl[i_p], kTRUE);
  73. sum = gausProb[0] + gausProb[1] + gausProb[2] + gausProb[3];
  74. if (sum == 0.) return 1;
  75. gausProb[0] /= sum;
  76. gausProb[1] /= sum;
  77. gausProb[2] /= sum;
  78. gausProb[3] /= sum;
  79. break;
  80. case 1: //bethe-bloch approximation with special byesian coefficients
  81. if (P > 3) i_p = 71;
  82. for (Int_t k = 0; k < Nintervals; k++) if (P >= P_int[k] && P < P_int[k + 1]) i_p = k;
  83. if (P < 0.3) {
  84. bayesAprioriCoefficients[0] = 0.28;
  85. bayesAprioriCoefficients[1] = 0.25;
  86. bayesAprioriCoefficients[2] = 0.26;
  87. bayesAprioriCoefficients[3] = 0.21;
  88. } else if ((P >= 0.3) && (P < 0.4)) {
  89. bayesAprioriCoefficients[0] = 0.91;
  90. bayesAprioriCoefficients[1] = 0.02;
  91. bayesAprioriCoefficients[2] = 0.06;
  92. bayesAprioriCoefficients[3] = 0.01;
  93. } else if ((P >= 0.4) && (P < 0.5)) {
  94. bayesAprioriCoefficients[0] = 0.70;
  95. bayesAprioriCoefficients[1] = 0.11;
  96. bayesAprioriCoefficients[2] = 0.13;
  97. bayesAprioriCoefficients[3] = 0.06;
  98. } else if ((P >= 0.5) && (P < 0.6)) {
  99. bayesAprioriCoefficients[0] = 0.39;
  100. bayesAprioriCoefficients[1] = 0.19;
  101. bayesAprioriCoefficients[2] = 0.32;
  102. bayesAprioriCoefficients[3] = 0.10;
  103. } else if ((P >= 0.6) && (P < 0.7)) {
  104. bayesAprioriCoefficients[0] = 0.41;
  105. bayesAprioriCoefficients[1] = 0.17;
  106. bayesAprioriCoefficients[2] = 0.37;
  107. bayesAprioriCoefficients[3] = 0.5;
  108. } else if ((P >= 0.7) && (P < 0.8)) {
  109. bayesAprioriCoefficients[0] = 0.43;
  110. bayesAprioriCoefficients[1] = 0.16;
  111. bayesAprioriCoefficients[2] = 0.31;
  112. bayesAprioriCoefficients[3] = 0.10;
  113. } else if ((P >= 0.8) && (P < 0.9)) {
  114. bayesAprioriCoefficients[0] = 0.44;
  115. bayesAprioriCoefficients[1] = 0.11;
  116. bayesAprioriCoefficients[2] = 0.43;
  117. bayesAprioriCoefficients[3] = 0.02;
  118. } else if ((P >= 0.9) && (P < 1.0)) {
  119. bayesAprioriCoefficients[0] = 0.36;
  120. bayesAprioriCoefficients[1] = 0.21;
  121. bayesAprioriCoefficients[2] = 0.36;
  122. bayesAprioriCoefficients[3] = 0.07;
  123. } else if ((P >= 1.0) && (P < 1.2)) {
  124. bayesAprioriCoefficients[0] = 0.32;
  125. bayesAprioriCoefficients[1] = 0.32;
  126. bayesAprioriCoefficients[2] = 0.32;
  127. bayesAprioriCoefficients[3] = 0.04;
  128. } else if (P >= 1.2) {
  129. bayesAprioriCoefficients[0] = 0.34;
  130. bayesAprioriCoefficients[1] = 0.27;
  131. bayesAprioriCoefficients[2] = 0.31;
  132. bayesAprioriCoefficients[3] = 0.08;
  133. }
  134. gausProb[0] = Gaus(dedx, BetheBlochFunction(P, PionPar), sigmasPi[i_p], kTRUE);
  135. gausProb[1] = Gaus(dedx, BetheBlochFunction(P, KaonPar), sigmasKa[i_p], kTRUE);
  136. gausProb[2] = Gaus(dedx, BetheBlochFunction(P, ProtonPar), sigmasPr[i_p], kTRUE);
  137. gausProb[3] = Gaus(dedx, BetheBlochFunction(P, ElectronPar), sigmasEl[i_p], kTRUE);
  138. sum = gausProb[0] + gausProb[1] + gausProb[2] + gausProb[3];
  139. if (sum == 0.) return 1;
  140. gausProb[0] /= sum;
  141. gausProb[1] /= sum;
  142. gausProb[2] /= sum;
  143. gausProb[3] /= sum;
  144. break;
  145. case 2: // n-sigma method,
  146. Float_t delta[Ntypes];
  147. bayesAprioriCoefficients[0] = 1.0;
  148. bayesAprioriCoefficients[1] = 1.0;
  149. bayesAprioriCoefficients[2] = 1.0;
  150. bayesAprioriCoefficients[3] = 1.0;
  151. if (P > 3) i_p = 71;
  152. for (Int_t k = 0; k < Nintervals; k++) if (P >= P_int[k] && P < P_int[k + 1]) i_p = k;
  153. delta[0]=TMath::Abs(dedx-BetheBlochFunction(P, PionPar))/sigmasPi[i_p];
  154. delta[1]=TMath::Abs(dedx-BetheBlochFunction(P, KaonPar))/sigmasKa[i_p];
  155. delta[2]=TMath::Abs(dedx-BetheBlochFunction(P, ProtonPar))/sigmasPr[i_p];
  156. delta[3]=TMath::Abs(dedx-BetheBlochFunction(P, ElectronPar))/sigmasEl[i_p];
  157. for(Int_t ii=0; ii<4; ii++) {
  158. if(delta[ii]<fn) gausProb[ii]=1;
  159. }
  160. Int_t c=0;
  161. for(Int_t ii=0; ii<4; ii++){
  162. c = gausProb[ii];
  163. for(Int_t jj=ii+1; jj<4; jj++) {
  164. if (c == 1 && gausProb[jj]==1) { gausProb[ii] = 0; gausProb[jj] = 0; }
  165. }
  166. }
  167. // cout<<"gausProb[ii] "<< gausProb[0]<<" "<<gausProb[1]<<" "<<gausProb[2]<<" "<<gausProb[3]<<endl;
  168. // cout << sigmasPi[i_p] << endl;
  169. break;
  170. }
  171. if(method!=0 && method!=1 && method!=2) {
  172. printf("No TPC method matches.\n");
  173. }
  174. else {
  175. Ppi=0.; Pk=0.; Pp=0.; Pe=0.;
  176. if(BayesFunction(gausProb, bayesAprioriCoefficients, resultProb, Ntypes)!=1){
  177. Ppi = resultProb[0];
  178. Pk = resultProb[1];
  179. Pp = resultProb[2];
  180. Pe = resultProb[3];
  181. }
  182. }
  183. return 0;
  184. }
  185. /*
  186. Float_t MpdTPCpid::BetheBlochFunction(Float_t x, Float_t *p) {
  187. //Float_t b = 1 / (x / Sqrt(1 + x * x));
  188. Float_t b = x / TMath::Sqrt(p[0]* p[0] + x * x);
  189. return p[1] / Power(b, p[4]) * (p[2] - Power(b, p[4]) - Log(p[3] + Power(1 / x, p[5])));
  190. }
  191. */
  192. Float_t MpdTPCpid::BetheBlochFunction(Float_t x, Float_t *par) {
  193. Float_t p=x;//kaon momentum
  194. Float_t mass=par[0];//Mass of particle
  195. Float_t beta=p/TMath::Sqrt(p*p+mass*mass);
  196. Float_t gamma=1/TMath::Sqrt(1-beta*beta);
  197. Float_t bg=beta*gamma;
  198. Float_t kp1=par[1];
  199. Float_t kp2=par[2];
  200. Float_t kp3=par[3];
  201. Float_t kp4=par[4];
  202. Float_t kp5=par[5];
  203. beta = bg/TMath::Sqrt(1.+ bg*bg);
  204. Float_t aa = TMath::Power(beta,kp4);
  205. Float_t bb = TMath::Power(1./bg,kp5);
  206. bb=TMath::Log(kp3+bb);
  207. return (kp2-aa-bb)*kp1/aa;
  208. }
  209. /*Bayes function for probabilities calculation*/
  210. Int_t MpdTPCpid::BayesFunction(Float_t *measProb, Float_t *aprioriProb, Float_t *bayesProb, Int_t N) {
  211. /*measProb - measured probabilities from TPC/TOF/etc*/
  212. /*aprioriProb - a priori probabilities from some magic*/
  213. /*bayesProb - result bayes probabilities */
  214. /*N - number of types {pion, kaon, proton, electron} */
  215. Float_t sumBayes = 0.0;
  216. for (Int_t i = 0; i < N; ++i) {
  217. sumBayes += measProb[i] * aprioriProb[i];
  218. }
  219. if (sumBayes == 0.) return 1;
  220. for (Int_t i = 0; i < N; ++i) {
  221. bayesProb[i] = measProb[i] * aprioriProb[i] / sumBayes;
  222. }
  223. return 0;
  224. }
  225. void MpdTPCpid::ReadTPCResponse() {
  226. std::string s(100,' ');
  227. /* Accessing the DeDx response */
  228. TString dir = getenv("VMCWORKDIR");
  229. TString dEdxFile = dir + "/input/MPDTPCPidResponseVhelle.txt";
  230. ifstream input;
  231. input.open(dEdxFile);
  232. if (!input.is_open()) {
  233. cerr<<"PID: READ TPC RESPONSE: Cannot open file "<<dEdxFile.Data()<<endl;
  234. }
  235. //const Int_t Ntypes = 4; //order: pion, kaon, proton, electron
  236. if(input) {
  237. getline(input, s);
  238. for(Int_t j=0; j<6; j++) input>>PionPar[j];
  239. getline(input,s);
  240. getline(input,s);
  241. for(Int_t j=0; j<6; j++) input>>KaonPar[j];
  242. getline(input,s);
  243. getline(input,s);
  244. for(Int_t j=0; j<6; j++) input>>ProtonPar[j];
  245. getline(input,s);
  246. getline(input,s);
  247. for(Int_t j=0; j<6; j++) input>>ElectronPar[j];
  248. getline(input, s);
  249. getline(input,s);
  250. for(Int_t j=0; j<72; j++) input>>sigmasPi[j];
  251. getline(input, s);
  252. getline(input,s);
  253. for(Int_t j=0; j<72; j++) input>>sigmasKa[j];
  254. getline(input, s);
  255. getline(input,s);
  256. for(Int_t j=0; j<72; j++) input>>sigmasPr[j];
  257. getline(input, s);
  258. getline(input,s);
  259. for(Int_t j=0; j<72; j++) input>>sigmasEl[j];
  260. input.close();
  261. }
  262. /*
  263. // Draw response
  264. const Int_t Nintervals = 72;
  265. Double_t P_int[Nintervals + 1] = {0.09,0.12,0.15,0.18,0.21,0.24,0.27,0.3,0.33,0.36,0.39,0.42,0.45,0.48,0.51,0.54,0.57,0.6,0.63,0.66,0.69,0.72,0.75,0.78,0.81,0.84,0.87,0.9,0.93,0.96,0.99,1.02,1.05,1.08,
  266. 1.11,1.14,1.17,1.2,1.23,1.26,1.29,1.32,1.35,1.38,1.41,1.44,1.47,1.5,1.53,1.56,1.59,1.62,1.65,1.68,1.71,1.74,1.77,1.8,1.83,1.86,1.89,1.92,1.95,1.98,2.01,2.04,2.07,2.1,
  267. 2.13,2.16,2.19,2.22, 3};
  268. Float_t P[72], P_err[72], mean_pi[72], mean_pr[72], mean_ka[72], mean_el[72];
  269. for(Int_t i=0; i<72; i++){
  270. P[i]=P_int[i]+0.5*(P_int[i+1]-P_int[i]);
  271. P_err[i]=0;
  272. mean_pi[i]=BetheBlochFunction(P[i], PionPar);
  273. mean_ka[i]=BetheBlochFunction(P[i], KaonPar);
  274. mean_pr[i]=BetheBlochFunction(P[i], ProtonPar);
  275. mean_el[i]=BetheBlochFunction(P[i], ElectronPar);
  276. }
  277. TGraphErrors *pi= new TGraphErrors(72, P, mean_pi, P_err, sigmasPi);
  278. pi->SetMarkerColor(2);
  279. pi->SetLineColor(2);
  280. pi->SetMarkerStyle(20);
  281. pi->SetMarkerSize(0.5);
  282. // pi->SetFillColor(1);
  283. // pi->SetFillStyle(3002);
  284. TGraphErrors *ka= new TGraphErrors(72, P, mean_ka, P_err, sigmasKa);
  285. ka->SetMarkerColor(2);
  286. ka->SetLineColor(2);
  287. ka->SetMarkerStyle(20);
  288. ka->SetMarkerSize(0.5);
  289. ka->SetFillColor(1);
  290. ka->SetFillStyle(3002);
  291. TGraphErrors *pr= new TGraphErrors(72, P, mean_pr, P_err, sigmasPr);
  292. pr->SetMarkerColor(2);
  293. pr->SetLineColor(2);
  294. pr->SetMarkerStyle(20);
  295. pr->SetMarkerSize(0.5);
  296. pr->SetFillColor(1);
  297. pr->SetFillStyle(3002);
  298. TGraphErrors *el= new TGraphErrors(72, P, mean_el, P_err, sigmasEl);
  299. el->SetMarkerColor(2);
  300. el->SetLineColor(2);
  301. el->SetMarkerStyle(20);
  302. el->SetMarkerSize(0.5);
  303. el->SetFillColor(1);
  304. el->SetFillStyle(3002);
  305. pr->Draw("C3same");
  306. ka->Draw("3Csame");
  307. pi->Draw("3Csame");
  308. el->Draw("3Csame");
  309. // pr->Draw("LXsame");
  310. // ka->Draw("LXsame");
  311. // pi->Draw("LXsame");
  312. // el->Draw("LXsame");
  313. */
  314. }
  315. ClassImp(MpdTPCpid);