MpdTOFpid.cxx 12 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358
  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 "MpdTOFpid.h"
  27. using namespace std;
  28. using namespace TMath;
  29. // Class Member definitions -----------
  30. MpdTOFpid::MpdTOFpid() :
  31. fnBeta(1.),
  32. ProtonParHyper(),
  33. PionParHyper(),
  34. KaonParHyper(),
  35. ElectronParHyper(),
  36. sigmasTofPi(),
  37. sigmasTofPr(),
  38. sigmasTofKa(),
  39. sigmasTofEl()
  40. {
  41. ReadTOFResponse();
  42. }
  43. MpdTOFpid::~MpdTOFpid() {
  44. }
  45. Int_t MpdTOFpid::GetTofProbs(Float_t P, Float_t beta, 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. if (Abs(beta) < 1e-6) return 1;
  53. const Float_t oneOverBeta = 1.0 / beta;
  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: //Hyperbolic (1/beta vs p) with bayesian coefficients
  63. bayesAprioriCoefficients[0] = 1.0;
  64. bayesAprioriCoefficients[1] = 1.0;
  65. bayesAprioriCoefficients[2] = 1.0;
  66. bayesAprioriCoefficients[3] = 1.0;
  67. if (P > 3) i_p = 71;
  68. for (Int_t k = 0; k < Nintervals; k++) if (P >= P_int[k] && P < P_int[k + 1]) i_p = k;
  69. gausProb[0] = Gaus(oneOverBeta, HyperbolicFunction(P, PionParHyper), sigmasTofPi[i_p], kTRUE);
  70. gausProb[1] = Gaus(oneOverBeta, HyperbolicFunction(P, KaonParHyper), sigmasTofKa[i_p], kTRUE);
  71. gausProb[2] = Gaus(oneOverBeta, HyperbolicFunction(P, ProtonParHyper), sigmasTofPr[i_p], kTRUE);
  72. gausProb[3] = Gaus(oneOverBeta, HyperbolicFunction(P, ElectronParHyper), sigmasTofEl[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: //Hyperbolic (1/beta vs p) 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(oneOverBeta, HyperbolicFunction(P, PionParHyper), sigmasTofPi[i_p], kTRUE);
  135. gausProb[1] = Gaus(oneOverBeta, HyperbolicFunction(P, KaonParHyper), sigmasTofKa[i_p], kTRUE);
  136. gausProb[2] = Gaus(oneOverBeta, HyperbolicFunction(P, ProtonParHyper), sigmasTofPr[i_p], kTRUE);
  137. gausProb[3] = Gaus(oneOverBeta, HyperbolicFunction(P, ElectronParHyper), sigmasTofEl[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(oneOverBeta - HyperbolicFunction(P, PionParHyper))/sigmasTofPi[i_p];
  154. delta[1]=TMath::Abs(oneOverBeta - HyperbolicFunction(P, KaonParHyper))/sigmasTofKa[i_p];
  155. delta[2]=TMath::Abs(oneOverBeta - HyperbolicFunction(P, ProtonParHyper))/sigmasTofPr[i_p];
  156. delta[3]=TMath::Abs(oneOverBeta - HyperbolicFunction(P, ElectronParHyper))/sigmasTofEl[i_p];
  157. for(Int_t ii=0; ii<4; ii++) {
  158. if(delta[ii]<fnBeta) 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. break;
  168. }
  169. if(method!=0 && method!=1 && method!=2) {
  170. printf("No TOF method matches.\n");
  171. }
  172. else {
  173. Ppi=0.; Pk=0.; Pp=0.; Pe=0.;
  174. if(BayesFunction(gausProb, bayesAprioriCoefficients, resultProb, Ntypes)!=1){
  175. Ppi = resultProb[0];
  176. Pk = resultProb[1];
  177. Pp = resultProb[2];
  178. Pe = resultProb[3];
  179. }
  180. }
  181. return 0;
  182. }
  183. Float_t MpdTOFpid::HyperbolicFunction(Float_t x, Float_t *p) {
  184. return Sqrt(1 + p[0]*p[0] / (x * x)) + p[1]; //p[0]=m, p[1]=additive const
  185. }
  186. /*Bayes function for probabilities calculation*/
  187. Int_t MpdTOFpid::BayesFunction(Float_t *measProb, Float_t *aprioriProb, Float_t *bayesProb, Int_t N) {
  188. /*measProb - measured probabilities from TOF/TOF/etc*/
  189. /*aprioriProb - a priori probabilities from some magic*/
  190. /*bayesProb - result bayes probabilities */
  191. /*N - number of types {pion, kaon, proton, electron} */
  192. Float_t sumBayes = 0.0;
  193. for (Int_t i = 0; i < N; ++i) {
  194. sumBayes += measProb[i] * aprioriProb[i];
  195. }
  196. if (sumBayes == 0.) return 1;
  197. for (Int_t i = 0; i < N; ++i) {
  198. bayesProb[i] = measProb[i] * aprioriProb[i] / sumBayes;
  199. }
  200. return 0;
  201. }
  202. void MpdTOFpid::ReadTOFResponse() {
  203. std::string s(100,' ');
  204. /* Accessing the oneOverBeta response */
  205. TString dir = getenv("VMCWORKDIR");
  206. TString BetaFile = dir + "/input/MPDTOFPidResponseVhelle.txt";
  207. ifstream input;
  208. input.open(BetaFile);
  209. if (!input.is_open()) {
  210. cerr<<"PID: READ TOF RESPONSE: Cannot open file "<<BetaFile.Data()<<endl;
  211. }
  212. //const Int_t Ntypes = 4; //order: pion, kaon, proton, electron
  213. if(input) {
  214. getline(input, s);
  215. for(Int_t j=0; j<2; j++) input>>PionParHyper[j];
  216. getline(input,s);
  217. getline(input,s);
  218. for(Int_t j=0; j<2; j++) input>>KaonParHyper[j];
  219. getline(input,s);
  220. getline(input,s);
  221. for(Int_t j=0; j<2; j++) input>>ProtonParHyper[j];
  222. getline(input,s);
  223. getline(input,s);
  224. for(Int_t j=0; j<2; j++) input>>ElectronParHyper[j];
  225. getline(input, s);
  226. getline(input,s);
  227. for(Int_t j=0; j<72; j++) input>>sigmasTofPi[j];
  228. getline(input, s);
  229. getline(input,s);
  230. for(Int_t j=0; j<72; j++) input>>sigmasTofKa[j];
  231. getline(input, s);
  232. getline(input,s);
  233. for(Int_t j=0; j<72; j++) input>>sigmasTofPr[j];
  234. getline(input, s);
  235. getline(input,s);
  236. for(Int_t j=0; j<72; j++) input>>sigmasTofEl[j];
  237. input.close();
  238. }
  239. // Draw response
  240. /*
  241. const Int_t Nintervals = 72;
  242. 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,
  243. 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,
  244. 2.13,2.16,2.19,2.22, 3};
  245. Float_t P[72], P_err[72], mean_pi[72], mean_pr[72], mean_ka[72], mean_el[72];
  246. for(Int_t i=0; i<72; i++){
  247. P[i]=P_int[i]+0.5*(P_int[i+1]-P_int[i]);
  248. P_err[i]=0;
  249. mean_pi[i]=HyperbolicFunction(P[i], PionParHyper);
  250. mean_ka[i]=HyperbolicFunction(P[i], KaonParHyper);
  251. mean_pr[i]=HyperbolicFunction(P[i], ProtonParHyper);
  252. mean_el[i]=HyperbolicFunction(P[i], ElectronParHyper);
  253. }
  254. TGraphErrors *pi= new TGraphErrors(72, P, mean_pi, P_err, sigmasTofPi);
  255. pi->SetMarkerColor(2);
  256. pi->SetLineColor(2);
  257. pi->SetMarkerStyle(20);
  258. pi->SetMarkerSize(0.5);
  259. pi->SetFillColor(1);
  260. pi->SetFillStyle(3002);
  261. TGraphErrors *ka= new TGraphErrors(72, P, mean_ka, P_err, sigmasTofKa);
  262. ka->SetMarkerColor(2);
  263. ka->SetLineColor(2);
  264. ka->SetMarkerStyle(20);
  265. ka->SetMarkerSize(0.5);
  266. ka->SetFillColor(1);
  267. ka->SetFillStyle(3002);
  268. TGraphErrors *pr= new TGraphErrors(72, P, mean_pr, P_err, sigmasTofPr);
  269. pr->SetMarkerColor(2);
  270. pr->SetLineColor(2);
  271. pr->SetMarkerStyle(20);
  272. pr->SetMarkerSize(0.5);
  273. pr->SetFillColor(1);
  274. pr->SetFillStyle(3002);
  275. TGraphErrors *el= new TGraphErrors(72, P, mean_el, P_err, sigmasTofEl);
  276. el->SetMarkerColor(2);
  277. el->SetLineColor(2);
  278. el->SetMarkerStyle(20);
  279. el->SetMarkerSize(0.5);
  280. el->SetFillColor(1);
  281. el->SetFillStyle(3002);
  282. pr->Draw("AC3");
  283. ka->Draw("3Csame");
  284. pi->Draw("3Csame");
  285. el->Draw("3Csame");
  286. // pr->Draw("LXsame");
  287. // ka->Draw("LXsame");
  288. // pi->Draw("LXsame");
  289. // el->Draw("LXsame");
  290. */
  291. }
  292. ClassImp(MpdTOFpid);