MpdParticleIdentification.cxx 23 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531
  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. // Sergey Merts
  14. //
  15. //-----------------------------------------------------------
  16. // C/C++ Headers ----------------------
  17. #include <iostream>
  18. #include <TMath.h>
  19. #include <TF1.h>
  20. // This Class' Header ------------------
  21. #include "MpdParticleIdentification.h"
  22. using namespace std;
  23. using namespace TMath;
  24. // Class Member definitions -----------
  25. MpdParticleIdentification::MpdParticleIdentification() {
  26. }
  27. MpdParticleIdentification::~MpdParticleIdentification() {
  28. }
  29. //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) {
  30. //
  31. // /*Parameters for fitting*/
  32. // //AZ Float_t *ProtonPar, *PionPar, *KaonPar, *ElectronPar;
  33. // /************************/
  34. //
  35. // dedx *= 1000000; // converting to keV/cm
  36. // const Int_t Ntypes = 4; //order: pion, kaon, proton, electron
  37. // Float_t ProtonPar[9], PionPar[9], KaonPar[9], ElectronPar[9], resultProb[Ntypes]; //AZ
  38. // const Int_t Nintervals = 10;
  39. // Float_t sigmas[Ntypes][Nintervals] = {//array of sigmas for different particles and different intervals of momentum. They were got from 100000 BOX events
  40. // {0.113, 0.107, 0.100, 0.107, 0.097, 0.111, 0.122, 0.105, 0.115, 0.118},
  41. // {0.407, 0.276, 0.260, 0.159, 0.185, 0.154, 0.122, 0.105, 0.115, 0.118},
  42. // {0.507, 0.488, 0.394, 0.330, 0.268, 0.244, 0.260, 0.172, 0.193, 0.118},
  43. // {0.163, 0.160, 0.149, 0.159, 0.172, 0.205, 0.260, 0.172, 0.193, 0.156}
  44. // };
  45. //
  46. // Float_t sigma = 1.0 / TMath::Sqrt(nHits); //for gaussian
  47. // Float_t sum = 0.0; // for normalizing
  48. //
  49. // /*A priori coefficients for Bayes approach */
  50. // Float_t bayesAprioriCoefficients[Ntypes];
  51. // /*A "measured" probabilities */
  52. // Float_t gausProb[Ntypes];
  53. //
  54. // Float_t sigPi, sigPr, sigKa, sigEl; //sigmas for gaussians
  55. //
  56. // switch (method) {
  57. //
  58. // case 0: //bethe-bloch approximation with "standard" sigmas and equal bayesian coefficients
  59. //
  60. // //parameters were got from Bethe-Bloch approximation for 100000 BOX events
  61. // /*AZ
  62. // ProtonPar = new Float_t[4];
  63. // PionPar = new Float_t[4];
  64. // KaonPar = new Float_t[4];
  65. // ElectronPar = new Float_t[4];
  66. // */
  67. // ProtonPar[0] = -3.957;
  68. // ProtonPar[1] = 3.525;
  69. // ProtonPar[2] = 16.468;
  70. // ProtonPar[3] = 0.815;
  71. // ProtonPar[4] = 5.247;
  72. // PionPar[0] = -0.739;
  73. // PionPar[1] = 7.373;
  74. // PionPar[2] = 3904.790;
  75. // PionPar[3] = 0.274;
  76. // PionPar[4] = 5.497;
  77. // KaonPar[0] = -2.590;
  78. // KaonPar[1] = 4.918;
  79. // KaonPar[2] = 79.722;
  80. // KaonPar[3] = 0.357;
  81. // KaonPar[4] = 4.511;
  82. // ElectronPar[0] = -1.552;
  83. // ElectronPar[1] = 1.748;
  84. // ElectronPar[2] = 7.425;
  85. // ElectronPar[3] = 0.980;
  86. // ElectronPar[4] = 1.604;
  87. //
  88. // sigma = 0.07 * dedx;
  89. // sigPi = sigPr = sigKa = sigEl = sigma;
  90. // bayesAprioriCoefficients[0] = 1.0;
  91. // bayesAprioriCoefficients[1] = 1.0;
  92. // bayesAprioriCoefficients[2] = 1.0;
  93. // bayesAprioriCoefficients[3] = 1.0;
  94. //
  95. // gausProb[0] = Gaus(dedx, BetheBlochFunction(P, PionPar), sigPi, kTRUE);
  96. // gausProb[1] = Gaus(dedx, BetheBlochFunction(P, KaonPar), sigKa, kTRUE);
  97. // gausProb[2] = Gaus(dedx, BetheBlochFunction(P, ProtonPar), sigPr, kTRUE);
  98. // gausProb[3] = Gaus(dedx, BetheBlochFunction(P, ElectronPar), sigEl, kTRUE);
  99. //
  100. // sum = gausProb[0] + gausProb[1] + gausProb[2] + gausProb[3];
  101. // if (sum == 0.) return 1;
  102. // gausProb[0] /= sum;
  103. // gausProb[1] /= sum;
  104. // gausProb[2] /= sum;
  105. // gausProb[3] /= sum;
  106. // break;
  107. //
  108. // case 1: //bethe-bloch approximation with special different sigmas and byesian coefficients
  109. //
  110. // //parameters were got from Bethe-Bloch approximation for 100000 BOX events
  111. // /*AZ
  112. // ProtonPar = new Float_t[4];
  113. // PionPar = new Float_t[4];
  114. // KaonPar = new Float_t[4];
  115. // ElectronPar = new Float_t[4];
  116. // */
  117. // ProtonPar[0] = -3.957;
  118. // ProtonPar[1] = 3.525;
  119. // ProtonPar[2] = 16.468;
  120. // ProtonPar[3] = 0.815;
  121. // ProtonPar[4] = 5.247;
  122. // PionPar[0] = -0.739;
  123. // PionPar[1] = 7.373;
  124. // PionPar[2] = 3904.790;
  125. // PionPar[3] = 0.274;
  126. // PionPar[4] = 5.497;
  127. // KaonPar[0] = -2.590;
  128. // KaonPar[1] = 4.918;
  129. // KaonPar[2] = 79.722;
  130. // KaonPar[3] = 0.357;
  131. // KaonPar[4] = 4.511;
  132. // ElectronPar[0] = -1.552;
  133. // ElectronPar[1] = 1.748;
  134. // ElectronPar[2] = 7.425;
  135. // ElectronPar[3] = 0.980;
  136. // ElectronPar[4] = 1.604;
  137. //
  138. // if (P < 0.3) {
  139. // bayesAprioriCoefficients[0] = 0.28;
  140. // bayesAprioriCoefficients[1] = 0.25;
  141. // bayesAprioriCoefficients[2] = 0.26;
  142. // bayesAprioriCoefficients[3] = 0.21;
  143. // sigPi = sigmas[0][0];
  144. // sigKa = sigmas[1][0];
  145. // sigPr = sigmas[2][0];
  146. // sigEl = sigmas[3][0];
  147. // } else if ((P >= 0.3) && (P < 0.4)) {
  148. // bayesAprioriCoefficients[0] = 0.91;
  149. // bayesAprioriCoefficients[1] = 0.02;
  150. // bayesAprioriCoefficients[2] = 0.06;
  151. // bayesAprioriCoefficients[3] = 0.01;
  152. // sigPi = sigmas[0][1];
  153. // sigKa = sigmas[1][1];
  154. // sigPr = sigmas[2][1];
  155. // sigEl = sigmas[3][1];
  156. // } else if ((P >= 0.4) && (P < 0.5)) {
  157. // bayesAprioriCoefficients[0] = 0.70;
  158. // bayesAprioriCoefficients[1] = 0.11;
  159. // bayesAprioriCoefficients[2] = 0.13;
  160. // bayesAprioriCoefficients[3] = 0.06;
  161. // sigPi = sigmas[0][2];
  162. // sigKa = sigmas[1][2];
  163. // sigPr = sigmas[2][2];
  164. // sigEl = sigmas[3][2];
  165. // } else if ((P >= 0.5) && (P < 0.6)) {
  166. // bayesAprioriCoefficients[0] = 0.39;
  167. // bayesAprioriCoefficients[1] = 0.19;
  168. // bayesAprioriCoefficients[2] = 0.32;
  169. // bayesAprioriCoefficients[3] = 0.10;
  170. // sigPi = sigmas[0][3];
  171. // sigKa = sigmas[1][3];
  172. // sigPr = sigmas[2][3];
  173. // sigEl = sigmas[3][3];
  174. // } else if ((P >= 0.6) && (P < 0.7)) {
  175. // bayesAprioriCoefficients[0] = 0.41;
  176. // bayesAprioriCoefficients[1] = 0.17;
  177. // bayesAprioriCoefficients[2] = 0.37;
  178. // bayesAprioriCoefficients[3] = 0.5;
  179. // sigPi = sigmas[0][4];
  180. // sigKa = sigmas[1][4];
  181. // sigPr = sigmas[2][4];
  182. // sigEl = sigmas[3][4];
  183. // } else if ((P >= 0.7) && (P < 0.8)) {
  184. // bayesAprioriCoefficients[0] = 0.43;
  185. // bayesAprioriCoefficients[1] = 0.16;
  186. // bayesAprioriCoefficients[2] = 0.31;
  187. // bayesAprioriCoefficients[3] = 0.10;
  188. // sigPi = sigmas[0][5];
  189. // sigKa = sigmas[1][5];
  190. // sigPr = sigmas[2][5];
  191. // sigEl = sigmas[3][5];
  192. // } else if ((P >= 0.8) && (P < 0.9)) {
  193. // bayesAprioriCoefficients[0] = 0.44;
  194. // bayesAprioriCoefficients[1] = 0.11;
  195. // bayesAprioriCoefficients[2] = 0.43;
  196. // bayesAprioriCoefficients[3] = 0.02;
  197. // sigPi = sigmas[0][6];
  198. // sigKa = sigmas[1][6];
  199. // sigPr = sigmas[2][6];
  200. // sigEl = sigmas[3][6];
  201. // } else if ((P >= 0.9) && (P < 1.0)) {
  202. // bayesAprioriCoefficients[0] = 0.36;
  203. // bayesAprioriCoefficients[1] = 0.21;
  204. // bayesAprioriCoefficients[2] = 0.36;
  205. // bayesAprioriCoefficients[3] = 0.07;
  206. // sigPi = sigmas[0][7];
  207. // sigKa = sigmas[1][7];
  208. // sigPr = sigmas[2][7];
  209. // sigEl = sigmas[3][7];
  210. // } else if ((P >= 1.0) && (P < 1.2)) {
  211. // bayesAprioriCoefficients[0] = 0.32;
  212. // bayesAprioriCoefficients[1] = 0.32;
  213. // bayesAprioriCoefficients[2] = 0.32;
  214. // bayesAprioriCoefficients[3] = 0.04;
  215. // sigPi = sigmas[0][8];
  216. // sigKa = sigmas[1][8];
  217. // sigPr = sigmas[2][8];
  218. // sigEl = sigmas[3][8];
  219. // } else if (P >= 1.2) {
  220. // bayesAprioriCoefficients[0] = 0.34;
  221. // bayesAprioriCoefficients[1] = 0.27;
  222. // bayesAprioriCoefficients[2] = 0.31;
  223. // bayesAprioriCoefficients[3] = 0.08;
  224. // sigPi = sigmas[0][9];
  225. // sigKa = sigmas[1][9];
  226. // sigPr = sigmas[2][9];
  227. // sigEl = sigmas[3][9];
  228. // }
  229. //
  230. // gausProb[0] = Gaus(dedx, BetheBlochFunction(P, PionPar), sigPi, kTRUE);
  231. // gausProb[1] = Gaus(dedx, BetheBlochFunction(P, KaonPar), sigKa, kTRUE);
  232. // gausProb[2] = Gaus(dedx, BetheBlochFunction(P, ProtonPar), sigPr, kTRUE);
  233. // gausProb[3] = Gaus(dedx, BetheBlochFunction(P, ElectronPar), sigEl, kTRUE);
  234. //
  235. // sum = gausProb[0] + gausProb[1] + gausProb[2] + gausProb[3];
  236. // if (sum == 0.) return 1;
  237. // gausProb[0] /= sum;
  238. // gausProb[1] /= sum;
  239. // gausProb[2] /= sum;
  240. // gausProb[3] /= sum;
  241. // break;
  242. //
  243. // case 2: //bethe-bloch approximation with "standard" sigmas and different byesian coefficients
  244. //
  245. // //parameters were got from Bethe-Bloch approximation for 100000 BOX events
  246. // /*AZ
  247. // ProtonPar = new Float_t[4];
  248. // PionPar = new Float_t[4];
  249. // KaonPar = new Float_t[4];
  250. // ElectronPar = new Float_t[4];
  251. // */
  252. // ProtonPar[0] = -3.957;
  253. // ProtonPar[1] = 3.525;
  254. // ProtonPar[2] = 16.468;
  255. // ProtonPar[3] = 0.815;
  256. // ProtonPar[4] = 5.247;
  257. // PionPar[0] = -0.739;
  258. // PionPar[1] = 7.373;
  259. // PionPar[2] = 3904.790;
  260. // PionPar[3] = 0.274;
  261. // PionPar[4] = 5.497;
  262. // KaonPar[0] = -2.590;
  263. // KaonPar[1] = 4.918;
  264. // KaonPar[2] = 79.722;
  265. // KaonPar[3] = 0.357;
  266. // KaonPar[4] = 4.511;
  267. // ElectronPar[0] = -1.552;
  268. // ElectronPar[1] = 1.748;
  269. // ElectronPar[2] = 7.425;
  270. // ElectronPar[3] = 0.980;
  271. // ElectronPar[4] = 1.604;
  272. //
  273. // sigma = 0.07 * dedx;
  274. // sigPi = sigPr = sigKa = sigEl = sigma;
  275. //
  276. // if (P < 0.3) {
  277. // bayesAprioriCoefficients[0] = 0.28;
  278. // bayesAprioriCoefficients[1] = 0.25;
  279. // bayesAprioriCoefficients[2] = 0.26;
  280. // bayesAprioriCoefficients[3] = 0.21;
  281. // } else if ((P >= 0.3) && (P < 0.4)) {
  282. // bayesAprioriCoefficients[0] = 0.91;
  283. // bayesAprioriCoefficients[1] = 0.02;
  284. // bayesAprioriCoefficients[2] = 0.06;
  285. // bayesAprioriCoefficients[3] = 0.01;
  286. // } else if ((P >= 0.4) && (P < 0.5)) {
  287. // bayesAprioriCoefficients[0] = 0.70;
  288. // bayesAprioriCoefficients[1] = 0.11;
  289. // bayesAprioriCoefficients[2] = 0.13;
  290. // bayesAprioriCoefficients[3] = 0.06;
  291. // } else if ((P >= 0.5) && (P < 0.6)) {
  292. // bayesAprioriCoefficients[0] = 0.39;
  293. // bayesAprioriCoefficients[1] = 0.19;
  294. // bayesAprioriCoefficients[2] = 0.32;
  295. // bayesAprioriCoefficients[3] = 0.10;
  296. // } else if ((P >= 0.6) && (P < 0.7)) {
  297. // bayesAprioriCoefficients[0] = 0.41;
  298. // bayesAprioriCoefficients[1] = 0.17;
  299. // bayesAprioriCoefficients[2] = 0.37;
  300. // bayesAprioriCoefficients[3] = 0.5;
  301. // } else if ((P >= 0.7) && (P < 0.8)) {
  302. // bayesAprioriCoefficients[0] = 0.43;
  303. // bayesAprioriCoefficients[1] = 0.16;
  304. // bayesAprioriCoefficients[2] = 0.31;
  305. // bayesAprioriCoefficients[3] = 0.10;
  306. // } else if ((P >= 0.8) && (P < 0.9)) {
  307. // bayesAprioriCoefficients[0] = 0.44;
  308. // bayesAprioriCoefficients[1] = 0.11;
  309. // bayesAprioriCoefficients[2] = 0.43;
  310. // bayesAprioriCoefficients[3] = 0.02;
  311. // } else if ((P >= 0.9) && (P < 1.0)) {
  312. // bayesAprioriCoefficients[0] = 0.36;
  313. // bayesAprioriCoefficients[1] = 0.21;
  314. // bayesAprioriCoefficients[2] = 0.36;
  315. // bayesAprioriCoefficients[3] = 0.07;
  316. // } else if ((P >= 1.0) && (P < 1.2)) {
  317. // bayesAprioriCoefficients[0] = 0.32;
  318. // bayesAprioriCoefficients[1] = 0.32;
  319. // bayesAprioriCoefficients[2] = 0.32;
  320. // bayesAprioriCoefficients[3] = 0.04;
  321. // } else if (P >= 1.2) {
  322. // bayesAprioriCoefficients[0] = 0.34;
  323. // bayesAprioriCoefficients[1] = 0.27;
  324. // bayesAprioriCoefficients[2] = 0.31;
  325. // bayesAprioriCoefficients[3] = 0.08;
  326. // }
  327. //
  328. // gausProb[0] = Gaus(dedx, BetheBlochFunction(P, PionPar), sigPi, kTRUE);
  329. // gausProb[1] = Gaus(dedx, BetheBlochFunction(P, KaonPar), sigKa, kTRUE);
  330. // gausProb[2] = Gaus(dedx, BetheBlochFunction(P, ProtonPar), sigPr, kTRUE);
  331. // gausProb[3] = Gaus(dedx, BetheBlochFunction(P, ElectronPar), sigEl, kTRUE);
  332. //
  333. // sum = gausProb[0] + gausProb[1] + gausProb[2] + gausProb[3];
  334. // if (sum == 0.) return 1;
  335. // gausProb[0] /= sum;
  336. // gausProb[1] /= sum;
  337. // gausProb[2] /= sum;
  338. // gausProb[3] /= sum;
  339. // break;
  340. //
  341. // case 3: //parabolic approximation
  342. //
  343. // //parameters were got from parabolic approximation function for 2000 UrQMD events
  344. // ProtonPar[0] = 0.031;
  345. // ProtonPar[1] = -0.124;
  346. // ProtonPar[2] = 1.412;
  347. // PionPar[0] = 0.230;
  348. // PionPar[1] = 0.088;
  349. // PionPar[2] = 1.072;
  350. // KaonPar[0] = 0.676;
  351. // KaonPar[1] = 0.621;
  352. // KaonPar[2] = 0.831;
  353. // ElectronPar[0] = 0.000;
  354. // ElectronPar[1] = -0.018;
  355. // ElectronPar[2] = 2.055;
  356. //
  357. // Float_t invP = 1.0 / P;
  358. //
  359. // gausProb[0] = Gaus(dedx, ParabolicFunction(invP, PionPar), sigPi, kTRUE);
  360. // gausProb[1] = Gaus(dedx, ParabolicFunction(invP, KaonPar), sigKa, kTRUE);
  361. // gausProb[2] = Gaus(dedx, ParabolicFunction(invP, ProtonPar), sigPr, kTRUE);
  362. // gausProb[3] = Gaus(dedx, ParabolicFunction(invP, ElectronPar), sigEl, kTRUE);
  363. //
  364. // sum = gausProb[0] + gausProb[1] + gausProb[2] + gausProb[3];
  365. // if (sum == 0.) return 1;
  366. // gausProb[0] /= sum;
  367. // gausProb[1] /= sum;
  368. // gausProb[2] /= sum;
  369. // gausProb[3] /= sum;
  370. //
  371. // break;
  372. // }
  373. // //AZ Float_t *resultProb = new Float_t[Ntypes];
  374. // BayesFunction(gausProb, bayesAprioriCoefficients, resultProb, Ntypes);
  375. //
  376. // Ppi = resultProb[0];
  377. // Pk = resultProb[1];
  378. // Pp = resultProb[2];
  379. // Pe = resultProb[3];
  380. //
  381. // return 0;
  382. //}
  383. /*
  384. 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) {
  385. Float_t ProtonParHyper[2], PionParHyper[2], KaonParHyper[2], ElectronParHyper[2];
  386. ProtonParHyper[0] = 0.938 * 0.938;
  387. ProtonParHyper[1] = -0.065; //-0.111791
  388. PionParHyper[0] = 0.139 * 0.139;
  389. PionParHyper[1] = -0.042;
  390. KaonParHyper[0] = 0.494 * 0.494;
  391. KaonParHyper[1] = -0.0491; //-0.0717
  392. ElectronParHyper[0] = 0.000511 * 0.000511;
  393. ElectronParHyper[1] = -0.028;
  394. const Int_t Ntypes = 4; //order: pion, kaon, proton, electron
  395. const Int_t Nintervals = 72;
  396. 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,
  397. 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,
  398. 2.13,2.16,2.19,2.22, 3};
  399. if (Abs(beta) < 1e-6) return 1;
  400. const Float_t oneOverBeta = 1.0 / beta;
  401. Float_t sigmas[Ntypes][Nintervals] = {//array of sigmas for different particles and different intervals of momentum.
  402. {0.0693274, 0.0693274, 0.0544018, 0.0434532, 0.0383988, 0.0374991, 0.0346434, 0.03343, 0.0319856, 0.032, 0.0315147, 0.0310805, 0.0298654, 0.0303749, 0.0301869, 0.0296072, 0.0289299, 0.0285614, 0.0298113, 0.0287559, 0.0269426, 0.0276648, 0.0281308, 0.0289822, 0.0260458, 0.0274102, 0.0276939, 0.0264013, 0.0251029, 0.0270367, 0.0269619, 0.0278648, 0.0247407, 0.0273073, 0.0262839, 0.0284527, 0.024663, 0.0270564, 0.0232594, 0.0256739, 0.0257952, 0.0237649, 0.0239696, 0.0248299, 0.0283596, 0.0264846, 0.022326, 0.0208957, 0.0307023, 0.0284478, 0.0229504, 0.0221302, 0.0251733, 0.0258165, 0.0273607, 0.0246405, 0.0229405, 0.024879, 0.0270787, 0.0238393, 0.026333, 0.0112759, 0.0404941, 0.0184759, 0.0225569, 0.0327113, 0.0221111, 0.0187046, 0.0194951, 0.0189474, 0.0211018, 0.00812245},
  403. {0.213016, 0.213016, 0.213016, 0.213016, 0.213016, 0.0989433, 0.0834294, 0.0675234, 0.0595375, 0.058213, 0.0486858, 0.0448145, 0.0471928, 0.0466786, 0.0451925, 0.0373988, 0.0451043, 0.0343869, 0.039328, 0.038438, 0.0364276, 0.034408, 0.0344429, 0.03274, 0.0368628, 0.0303875, 0.0312515, 0.0331831, 0.0302307, 0.0310654, 0.0382999, 0.0311076, 0.0396424, 0.0350253, 0.0319718, 0.0367985, 0.0300645, 0.0263978, 0.0261712, 0.0314639, 0.0263756, 0.0338912, 0.0313582, 0.0200379, 0.0265785, 0.0300492, 0.0308032, 0.03746, 0.0220187, 0.0222599, 0.0244639, 0.0284196, 0.0286862, 0.0407272, 0.012562, 0.023188, 0.012329, 0.0297552, 0.0208175, 0.0426947, 0.0267256, 0.0157099, 0.0404566, 0.0162424, 0.0144609, 0.0304202, 0.019353, 0.0527709, 0.013652, 0.187007, 0.28284, 0.0161017},
  404. {0.220549, 0.220549, 0.220549, 0.220549, 0.220549, 0.220549, 0.220549, 0.220549, 0.220549, 0.094147, 0.0942735, 0.0958616, 0.0738113, 0.0776133, 0.0631406, 0.0543964, 0.0672916, 0.0562169, 0.0625832, 0.0549343, 0.0521649, 0.0535391, 0.0529288, 0.0466422, 0.0523035, 0.0472841, 0.0406543, 0.0460797, 0.0530522, 0.0439447, 0.0426407, 0.0495772, 0.0406714, 0.0341736, 0.0376888, 0.0472375, 0.0441835, 0.0363689, 0.0425563, 0.0421598, 0.0371204, 0.0462698, 0.0360139, 0.0384682, 0.0400049, 0.0387571, 0.0328137, 0.0384007, 0.0351909, 0.0594957, 0.0310895, 0.0357163, 0.0387537, 0.0358798, 0.0385559, 0.0277795, 0.0348187, 0.0269486, 0.0523929, 0.0282686, 0.0339755, 0.0238928, 0.0241012, 0.0337723, 0.0334852, 0.0248415, 0.0534019, 0.0270102, 0.0403495, 0.0284284, 0.0289664, 0.0240482},
  405. {0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01}
  406. };
  407. Int_t i_p = 0; // momentum interval
  408. Float_t sum = 0.0; // for normalizing
  409. Float_t gausProb[Ntypes] = {0.0, 0.0, 0.0, 0.0};
  410. switch (method) {
  411. case 0: //Hyperbolic (1/beta vs p)
  412. if (P > 3) i_p = 71;
  413. for (Int_t k = 0; k < Nintervals; k++) if (P >= P_int[k] && P < P_int[k + 1]) i_p = k;
  414. gausProb[0] = Gaus(oneOverBeta, HyperbolicFunction(P, PionParHyper), sigmas[0][i_p], kTRUE);
  415. gausProb[1] = Gaus(oneOverBeta, HyperbolicFunction(P, KaonParHyper), sigmas[1][i_p], kTRUE);
  416. gausProb[2] = Gaus(oneOverBeta, HyperbolicFunction(P, ProtonParHyper), sigmas[2][i_p], kTRUE);
  417. gausProb[3] = Gaus(oneOverBeta, HyperbolicFunction(P, ElectronParHyper), sigmas[3][i_p], kTRUE);
  418. sum = gausProb[0] + gausProb[1] + gausProb[2] + gausProb[3];
  419. if (sum == 0.) return 1;
  420. gausProb[0] /= sum;
  421. gausProb[1] /= sum;
  422. gausProb[2] /= sum;
  423. gausProb[3] /= sum;
  424. break;
  425. case 1:
  426. break;
  427. }
  428. Ppi = gausProb[0];
  429. Pk = gausProb[1];
  430. Pp = gausProb[2];
  431. Pe = gausProb[3];
  432. return 0;
  433. }
  434. */
  435. Int_t MpdParticleIdentification::GetCombinedProbs(Float_t *tofProbs, Float_t *tpcProbs, Float_t *resultProbs, Int_t N) {
  436. /*tofProbs - measured probabilities from TOF*/
  437. /*tpcProbs - measured probabilities from TPC*/
  438. /*resultProbs - combined probabilities */
  439. /*N - number of types {pion, kaon, proton, electron} */
  440. Float_t sumBayes = 0.0;
  441. for (Int_t i = 0; i < N; ++i) {
  442. sumBayes += tofProbs[i] * tpcProbs[i];
  443. }
  444. for (Int_t i = 0; i < N; ++i) {
  445. if ((tpcProbs[i] != 0.0) && (tofProbs[i] != 0.0)) {
  446. resultProbs[i] = tofProbs[i] * tpcProbs[i] / sumBayes;
  447. } else if (tpcProbs[i] == 0.0) {
  448. resultProbs[i] = tofProbs[i] / sumBayes;
  449. } else if (tofProbs[i] == 0.0) {
  450. resultProbs[i] = tpcProbs[i] / sumBayes;
  451. }
  452. }
  453. return 0;
  454. }
  455. //Float_t MpdParticleIdentification::BetheBlochFunction(Float_t x, Float_t *p) {
  456. // Float_t b = 1 / (x / Sqrt(1 + x * x));
  457. // return p[0] / Power(b, p[3]) * (p[1] - Power(b, p[3]) - Log(p[2] + Power(1 / x, p[4])));
  458. //}
  459. Float_t MpdParticleIdentification::ParabolicFunction(Float_t x, Float_t *p) {
  460. return p[0] * x * x + p[1] * x + p[2];
  461. }
  462. /*
  463. Float_t MpdParticleIdentification::HyperbolicFunction(Float_t x, Float_t *p) {
  464. return Sqrt(1 + p[0] / (x * x)) + p[1]; //p[0]=m^2, p[1]=additive const
  465. }
  466. */
  467. /*
  468. // Bayes function for probabilities calculation //
  469. Int_t MpdParticleIdentification::BayesFunction(Float_t *measProb, Float_t *aprioriProb, Float_t *bayesProb, Int_t N) {
  470. // measProb - measured probabilities from TPC/TOF/etc //
  471. // aprioriProb - a priori probabilities from some magic //
  472. // bayesProb - result bayes probabilities //
  473. // N - number of types {pion, kaon, proton, electron} //
  474. Float_t sumBayes = 0.0;
  475. for (Int_t i = 0; i < N; ++i) {
  476. sumBayes += measProb[i] * aprioriProb[i];
  477. }
  478. if (sumBayes == 0.) return 1;
  479. for (Int_t i = 0; i < N; ++i) {
  480. bayesProb[i] = measProb[i] * aprioriProb[i] / sumBayes;
  481. }
  482. return 0;
  483. }
  484. */
  485. ClassImp(MpdParticleIdentification);