MpdPidQA.cxx 59 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921
  1. #include "MpdPidQA.h"
  2. MpdPidQA::MpdPidQA() : MpdPid() {
  3. /// Default constructor
  4. fSumBetheBlochHists = nullptr;
  5. fChBetheBlochHists = nullptr;
  6. fm2LightHist = nullptr;
  7. fm2HeavyHist = nullptr;
  8. Init("pikapr");
  9. }
  10. MpdPidQA::MpdPidQA(Double_t sigmaTof, Double_t sigmaEloss, Double_t sqrts,
  11. Double_t EnLossCoef, TString Generator, TString Tracking, TString NSigPart)
  12. : MpdPid(sigmaTof, sigmaEloss, sqrts, EnLossCoef, Generator, Tracking, NSigPart)
  13. {
  14. Init(NSigPart);
  15. }
  16. MpdPidQA::~MpdPidQA() {
  17. for ( auto it = fEnLossMap.begin(); it != fEnLossMap.end(); ++it ) {
  18. for ( Int_t i = 0; i < MpdPidUtils::nQAHists; i++ ) delete it->second.dEdXPart[i]; delete it->second.BetheBlochHist;
  19. }
  20. for ( auto it = fMSquaredMap.begin(); it != fMSquaredMap.end(); ++it ) delete it->second;
  21. for ( auto it = fAbundanceMap.begin(); it != fAbundanceMap.end(); ++it ) {
  22. for ( TH1D* hYield : it->second ) {
  23. delete hYield;
  24. }
  25. }
  26. for ( auto it = fEffContMap.begin(); it != fEffContMap.end(); ++it ) {
  27. for ( vecTH1Dptrs vecEC : it->second ) {
  28. for ( TH1D* hEC : vecEC ) {
  29. delete hEC;
  30. }
  31. }
  32. }
  33. delete fSumBetheBlochHists;
  34. delete fChBetheBlochHists;
  35. delete fm2LightHist;
  36. delete fm2HeavyHist;
  37. }
  38. void MpdPidQA::Init(TString Particles) {
  39. cout << "MpdPidQA::Init().." << endl;
  40. Double_t PMIN = 0.0, PMAX = 5.0;
  41. Double_t XMIN[MpdPidUtils::nQAHists][MpdPidUtils::kNSpecies], XMAX[MpdPidUtils::nQAHists][MpdPidUtils::kNSpecies]; TString histName;
  42. const Int_t nbins = 200, nbins2 = 100;
  43. for (Int_t i = 0; i < MpdPidUtils::nQAHists; i++) { for (Int_t j = 0; j < MpdPidUtils::kNSpecies; j++) { XMIN[i][j] = -1.; XMAX[i][j] = -1.; } }
  44. TString ECHistNames[4] = { "All", "Id", "IdRight", "IdWrong" };
  45. nSigPart = Particles;
  46. Double_t MinEnLoss = 0.0, MaxEnLoss = 0.0, MaxEnLossLoc;
  47. /// filling out fPartTypeMap
  48. fPartTypeMap.insert( pair<Int_t,MpdPidUtils::ePartType>(11,MpdPidUtils::kElectron) );
  49. fPartTypeMap.insert( pair<Int_t,MpdPidUtils::ePartType>(-11,MpdPidUtils::kElectron) );
  50. fPartTypeMap.insert( pair<Int_t,MpdPidUtils::ePartType>(13,MpdPidUtils::kMuon) );
  51. fPartTypeMap.insert( pair<Int_t,MpdPidUtils::ePartType>(-13,MpdPidUtils::kMuon) );
  52. fPartTypeMap.insert( pair<Int_t,MpdPidUtils::ePartType>(211,MpdPidUtils::kPion) );
  53. fPartTypeMap.insert( pair<Int_t,MpdPidUtils::ePartType>(-211,MpdPidUtils::kPion) );
  54. fPartTypeMap.insert( pair<Int_t,MpdPidUtils::ePartType>(321,MpdPidUtils::kKaon) );
  55. fPartTypeMap.insert( pair<Int_t,MpdPidUtils::ePartType>(-321,MpdPidUtils::kKaon) );
  56. fPartTypeMap.insert( pair<Int_t,MpdPidUtils::ePartType>(2212,MpdPidUtils::kProton) );
  57. fPartTypeMap.insert( pair<Int_t,MpdPidUtils::ePartType>(-2212,MpdPidUtils::kProton) );
  58. fPartTypeMap.insert( pair<Int_t,MpdPidUtils::ePartType>(PDG_DEUTERON,MpdPidUtils::kDeuteron) );
  59. fPartTypeMap.insert( pair<Int_t,MpdPidUtils::ePartType>(PDG_TRITON,MpdPidUtils::kTriton) );
  60. fPartTypeMap.insert( pair<Int_t,MpdPidUtils::ePartType>(PDG_HE3,MpdPidUtils::kHe3) );
  61. fPartTypeMap.insert( pair<Int_t,MpdPidUtils::ePartType>(PDG_HE4,MpdPidUtils::kHe4) );
  62. if (Particles.Contains("el")) {
  63. /// dE/dx
  64. if ( ( fTrackingState == MpdPidUtils::kCF ) || ( fTrackingState == MpdPidUtils::kCFHM ) ) {
  65. MaxEnLossLoc = 5500.;
  66. for ( Int_t i = 0; i < 5; i++ ) { XMIN[i][MpdPidUtils::kElectron] = 2000.; XMAX[i][MpdPidUtils::kElectron] = 5000.; }
  67. for ( Int_t i = 5; i < MpdPidUtils::nQAHists; i++ ) { XMIN[i][MpdPidUtils::kElectron] = 2500.; XMAX[i][MpdPidUtils::kElectron] = 5500.; }
  68. } else {
  69. MaxEnLossLoc = 3.e-06;
  70. for ( Int_t i = 0; i < MpdPidUtils::nQAHists; i++ ) { XMIN[i][MpdPidUtils::kElectron] = 1.e-06; XMAX[i][MpdPidUtils::kElectron] = 3.e-06; }
  71. }
  72. MpdPidUtils::dEdXStruct_t SEl;
  73. for ( Int_t iHist = 0; iHist < MpdPidUtils::nQAHists; iHist++ ) {
  74. histName = "El_" + TString::Itoa(iHist, 10);
  75. SEl.dEdXPart[iHist] = new TH1D(histName, "", nbins, XMIN[iHist][MpdPidUtils::kElectron], XMAX[iHist][MpdPidUtils::kElectron]);
  76. }
  77. SEl.BetheBlochHist = new TH2D("hElBB_QA", "", 300, 0.05, 3.0, 300, MinEnLoss, MaxEnLossLoc);
  78. SEl.ibeg = 0; SEl.iend = 39;
  79. fEnLossMap.insert( pair<MpdPidUtils::ePartType,MpdPidUtils::dEdXStruct_t>(MpdPidUtils::kElectron,SEl) );
  80. /// m-squared
  81. TH2D *mass2El = new TH2D("mass2El", "", nbins, 0., 4., nbins, -0.7, 0.7);
  82. fMSquaredMap.insert( pair<MpdPidUtils::ePartType,TH2D*>(MpdPidUtils::kElectron,mass2El) );
  83. /// Abundance
  84. vecTH1Dptrs vecYield;
  85. TH1D *hYield_elpos = new TH1D("hYield_elpos","",600,0.,3.); vecYield.push_back(hYield_elpos);
  86. TH1D *hYield_elneg = new TH1D("hYield_elneg","",600,0.,3.); vecYield.push_back(hYield_elneg);
  87. fAbundanceMap.insert( pair<MpdPidUtils::ePartType,vecTH1Dptrs>(MpdPidUtils::kElectron,vecYield) );
  88. /// Efficiency & Contamination
  89. TH1D *hEC_pos[4], *hEC_neg[4];
  90. vecTH1Dptrs vecEC_pos, vecEC_neg;
  91. std::vector<vecTH1Dptrs> vecvecEC;
  92. for ( Int_t iHist = 0; iHist < 4; iHist++ ) {
  93. hEC_pos[iHist] = new TH1D( ECHistNames[iHist] + "_poschar_" + MpdPidUtils::cParticleShortName[MpdPidUtils::kElectron],"", 50, PMIN, PMAX);
  94. hEC_pos[iHist]->Sumw2();
  95. vecEC_pos.push_back( hEC_pos[iHist] );
  96. hEC_neg[iHist] = new TH1D( ECHistNames[iHist] + "_negchar_" + MpdPidUtils::cParticleShortName[MpdPidUtils::kElectron], "", 50, PMIN, PMAX);
  97. hEC_neg[iHist]->Sumw2();
  98. vecEC_neg.push_back( hEC_neg[iHist] );
  99. }
  100. vecvecEC.push_back( vecEC_pos );
  101. vecvecEC.push_back( vecEC_neg );
  102. fEffContMap.insert( pair<MpdPidUtils::ePartType,std::vector<vecTH1Dptrs>>(MpdPidUtils::kElectron,vecvecEC) );
  103. }
  104. if (Particles.Contains("mu"))
  105. {
  106. /// dE/dx
  107. if ( ( fTrackingState == MpdPidUtils::kCF ) || ( fTrackingState == MpdPidUtils::kCFHM ) ) {
  108. MaxEnLossLoc = 15000.;
  109. XMIN[0][MpdPidUtils::kMuon] = 2000.; XMIN[1][MpdPidUtils::kMuon] = 1000.; XMAX[0][MpdPidUtils::kMuon] = 15000.; XMAX[1][MpdPidUtils::kMuon] = 8000.;
  110. for ( Int_t i = 2; i < MpdPidUtils::nQAHists; i++ ) { XMIN[i][MpdPidUtils::kMuon] = 1000.; XMAX[i][MpdPidUtils::kMuon] = 5000.; }
  111. } else {
  112. MaxEnLossLoc = 15.e-06;
  113. XMAX[0][MpdPidUtils::kMuon] = 15.e-06; XMAX[1][MpdPidUtils::kMuon] = 5.e-06;
  114. for ( Int_t i = 0; i < MpdPidUtils::nQAHists; i++ ) XMIN[i][MpdPidUtils::kMuon] = 1.e-06;
  115. for ( Int_t i = 2; i < MpdPidUtils::nQAHists; i++ ) XMAX[i][MpdPidUtils::kMuon] = 3.e-06;
  116. }
  117. MpdPidUtils::dEdXStruct_t SMu;
  118. for ( Int_t iHist = 0; iHist < MpdPidUtils::nQAHists; iHist++ ) {
  119. histName = "Mu_" + TString::Itoa(iHist, 10);
  120. SMu.dEdXPart[iHist] = new TH1D(histName, "", nbins, XMIN[iHist][MpdPidUtils::kMuon], XMAX[iHist][MpdPidUtils::kMuon]);
  121. }
  122. SMu.BetheBlochHist = new TH2D("hMuBB_QA", "", 300, 0.05, 3.0, 300, MinEnLoss, MaxEnLossLoc);
  123. SMu.ibeg = 0; SMu.iend = 39;
  124. fEnLossMap.insert( pair<MpdPidUtils::ePartType,MpdPidUtils::dEdXStruct_t>(MpdPidUtils::kMuon,SMu) );
  125. /// m-squared
  126. TH2D *mass2Mu = new TH2D("mass2Mu", "", nbins, 0., 4., nbins, -0.8, 0.8);
  127. fMSquaredMap.insert( pair<MpdPidUtils::ePartType,TH2D*>(MpdPidUtils::kMuon,mass2Mu) );
  128. /// Abundance
  129. vecTH1Dptrs vecYield;
  130. TH1D *hYield_mupos = new TH1D("hYield_mupos","",600,0.,3.); vecYield.push_back(hYield_mupos);
  131. TH1D *hYield_muneg = new TH1D("hYield_muneg","",600,0.,3.); vecYield.push_back(hYield_muneg);
  132. fAbundanceMap.insert( pair<MpdPidUtils::ePartType,vecTH1Dptrs>(MpdPidUtils::kMuon,vecYield) );
  133. /// Efficiency & Contamination
  134. TH1D *hEC_pos[4], *hEC_neg[4];
  135. vecTH1Dptrs vecEC_pos, vecEC_neg;
  136. std::vector<vecTH1Dptrs> vecvecEC;
  137. for ( Int_t iHist = 0; iHist < 4; iHist++ ) {
  138. hEC_pos[iHist] = new TH1D( ECHistNames[iHist] + "_poschar_" + MpdPidUtils::cParticleShortName[MpdPidUtils::kMuon],"", 50, PMIN, PMAX);
  139. hEC_pos[iHist]->Sumw2();
  140. vecEC_pos.push_back( hEC_pos[iHist] );
  141. hEC_neg[iHist] = new TH1D( ECHistNames[iHist] + "_negchar_" + MpdPidUtils::cParticleShortName[MpdPidUtils::kMuon], "", 50, PMIN, PMAX);
  142. hEC_neg[iHist]->Sumw2();
  143. vecEC_neg.push_back( hEC_neg[iHist] );
  144. }
  145. vecvecEC.push_back( vecEC_pos );
  146. vecvecEC.push_back( vecEC_neg );
  147. fEffContMap.insert( pair<MpdPidUtils::ePartType,std::vector<vecTH1Dptrs>>(MpdPidUtils::kMuon,vecvecEC) );
  148. }
  149. if (Particles.Contains("pi")) {
  150. /// dE/dx
  151. if ( fTrackingState == MpdPidUtils::kCFHM ) {
  152. MaxEnLossLoc = 14.;
  153. for ( Int_t i = 0; i < 9; i++ ) XMIN[i][MpdPidUtils::kPion] = 0.7;
  154. for ( Int_t i = 9; i < 17; i++ ) XMIN[i][MpdPidUtils::kPion] = 0.8;
  155. for ( Int_t i = 17; i < 30; i++ ) XMIN[i][MpdPidUtils::kPion] = 0.9;
  156. for ( Int_t i = 30; i < MpdPidUtils::nQAHists; i++ ) XMIN[i][MpdPidUtils::kPion] = 1.0;
  157. XMAX[0][MpdPidUtils::kPion] = 14.; XMAX[1][MpdPidUtils::kPion] = 6.0; XMAX[2][MpdPidUtils::kPion] = 4.0; XMAX[3][MpdPidUtils::kPion] = 4.0;
  158. XMAX[4][MpdPidUtils::kPion] = 3.0; XMAX[5][MpdPidUtils::kPion] = 2.5; XMAX[6][MpdPidUtils::kPion] = 2.5; XMAX[16][MpdPidUtils::kPion] = 2.2;
  159. for ( Int_t i = 7; i < 13; i++ ) XMAX[i][MpdPidUtils::kPion] = 2.4;
  160. for ( Int_t i = 13; i < 16; i++ ) XMAX[i][MpdPidUtils::kPion] = 2.3;
  161. for ( Int_t i = 17; i < 35; i++ ) XMAX[i][MpdPidUtils::kPion] = 2.1;
  162. for ( Int_t i = 35; i < 38; i++ ) XMAX[i][MpdPidUtils::kPion] = 2.2;
  163. for ( Int_t i = 38; i < MpdPidUtils::nQAHists; i++ ) XMAX[i][MpdPidUtils::kPion] = 2.3;
  164. } else if ( fTrackingState == MpdPidUtils::kCF ) {
  165. MaxEnLossLoc = 25000.;
  166. XMIN[0][MpdPidUtils::kPion] = 1000.; XMAX[0][MpdPidUtils::kPion] = 25000.; XMAX[1][MpdPidUtils::kPion] = 10000.; XMAX[2][MpdPidUtils::kPion] = 8000.;
  167. for ( Int_t i = 1; i < MpdPidUtils::nQAHists; i++ ) XMIN[i][MpdPidUtils::kPion] = 2000.;
  168. for ( Int_t i = 3; i < MpdPidUtils::nQAHists; i++ ) XMAX[i][MpdPidUtils::kPion] = 6000.;
  169. } else {
  170. MaxEnLossLoc = 15.e-06;
  171. XMAX[0][MpdPidUtils::kPion] = 15.e-06; XMAX[1][MpdPidUtils::kPion] = 6.e-06; XMAX[2][MpdPidUtils::kPion] = 3.5e-06; XMAX[3][MpdPidUtils::kPion] = 3.e-06;
  172. XMAX[4][MpdPidUtils::kPion] = 2.5e-06; XMAX[5][MpdPidUtils::kPion] = 2.5e-06; XMAX[6][MpdPidUtils::kPion] = 2.2e-06;
  173. for ( Int_t i = 0; i < MpdPidUtils::nQAHists; i++ ) XMIN[i][MpdPidUtils::kPion] = 1.e-06;
  174. for ( Int_t i = 7; i < MpdPidUtils::nQAHists; i++ ) XMAX[i][MpdPidUtils::kPion] = 2.1e-06;
  175. }
  176. MpdPidUtils::dEdXStruct_t SPi;
  177. for ( Int_t iHist = 0; iHist < MpdPidUtils::nQAHists; iHist++ ) {
  178. histName = "Pi_" + TString::Itoa(iHist, 10);
  179. SPi.dEdXPart[iHist] = new TH1D(histName, "", nbins, XMIN[iHist][MpdPidUtils::kPion], XMAX[iHist][MpdPidUtils::kPion]);
  180. }
  181. SPi.BetheBlochHist = new TH2D("hPiBB_QA", "", 300, 0.05, 3.0, 300, MinEnLoss, MaxEnLossLoc);
  182. SPi.ibeg = 0; SPi.iend = 39;
  183. fEnLossMap.insert( pair<MpdPidUtils::ePartType,MpdPidUtils::dEdXStruct_t>(MpdPidUtils::kPion,SPi) );
  184. /// m-squared
  185. TH2D *mass2Pi = new TH2D("mass2Pi", "", nbins, 0., 4., nbins, -0.5, 1.5);
  186. fMSquaredMap.insert( pair<MpdPidUtils::ePartType,TH2D*>(MpdPidUtils::kPion,mass2Pi) );
  187. /// Abundance
  188. vecTH1Dptrs vecYield;
  189. TH1D *hYield_pipos = new TH1D("hYield_pipos","",600,0.,3.); vecYield.push_back(hYield_pipos);
  190. TH1D *hYield_pineg = new TH1D("hYield_pineg","",600,0.,3.); vecYield.push_back(hYield_pineg);
  191. fAbundanceMap.insert( pair<MpdPidUtils::ePartType,vecTH1Dptrs>(MpdPidUtils::kPion,vecYield) );
  192. /// Efficiency & Contamination
  193. TH1D *hEC_pos[4], *hEC_neg[4];
  194. vecTH1Dptrs vecEC_pos, vecEC_neg;
  195. std::vector<vecTH1Dptrs> vecvecEC;
  196. for ( Int_t iHist = 0; iHist < 4; iHist++ ) {
  197. hEC_pos[iHist] = new TH1D( ECHistNames[iHist] + "_poschar_" + MpdPidUtils::cParticleShortName[MpdPidUtils::kPion],"", 50, PMIN, PMAX);
  198. hEC_pos[iHist]->Sumw2();
  199. vecEC_pos.push_back( hEC_pos[iHist] );
  200. hEC_neg[iHist] = new TH1D( ECHistNames[iHist] + "_negchar_" + MpdPidUtils::cParticleShortName[MpdPidUtils::kPion], "", 50, PMIN, PMAX);
  201. hEC_neg[iHist]->Sumw2();
  202. vecEC_neg.push_back( hEC_neg[iHist] );
  203. }
  204. vecvecEC.push_back( vecEC_pos );
  205. vecvecEC.push_back( vecEC_neg );
  206. fEffContMap.insert( pair<MpdPidUtils::ePartType,std::vector<vecTH1Dptrs>>(MpdPidUtils::kPion,vecvecEC) );
  207. }
  208. if (Particles.Contains("ka")) {
  209. /// dE/dx
  210. if ( fTrackingState == MpdPidUtils::kCFHM ) {
  211. MaxEnLossLoc = 70;
  212. XMIN[0][MpdPidUtils::kKaon] = 10.; XMIN[1][MpdPidUtils::kKaon] = 6.; XMIN[2][MpdPidUtils::kKaon] = 5.; XMIN[3][MpdPidUtils::kKaon] = 3.;
  213. XMIN[4][MpdPidUtils::kKaon] = 2.; XMIN[5][MpdPidUtils::kKaon] = 2.; XMIN[21][MpdPidUtils::kKaon] = 0.9; XMIN[22][MpdPidUtils::kKaon] = 0.9;
  214. for ( Int_t i = 6; i < 21; i++ ) XMIN[i][MpdPidUtils::kKaon] = 1.;
  215. for ( Int_t i = 23; i < MpdPidUtils::nQAHists; i++ ) XMIN[i][MpdPidUtils::kKaon] = 0.8;
  216. XMAX[0][MpdPidUtils::kKaon] = 70.; XMAX[1][MpdPidUtils::kKaon] = 40.; XMAX[2][MpdPidUtils::kKaon] = 20.; XMAX[3][MpdPidUtils::kKaon] = 12.;
  217. XMAX[4][MpdPidUtils::kKaon] = 8.; XMAX[5][MpdPidUtils::kKaon] = 6.; XMAX[6][MpdPidUtils::kKaon] = 5.5; XMAX[7][MpdPidUtils::kKaon] = 4.5;
  218. XMAX[8][MpdPidUtils::kKaon] = 4.; XMAX[9][MpdPidUtils::kKaon] = 3.5; XMAX[10][MpdPidUtils::kKaon] = 3.5; XMAX[11][MpdPidUtils::kKaon] = 3.;
  219. XMAX[12][MpdPidUtils::kKaon] = 2.8; XMAX[13][MpdPidUtils::kKaon] = 2.6; XMAX[14][MpdPidUtils::kKaon] = 2.4; XMAX[15][MpdPidUtils::kKaon] = 2.3;
  220. XMAX[16][MpdPidUtils::kKaon] = 2.2; XMAX[17][MpdPidUtils::kKaon] = 2.2; XMAX[18][MpdPidUtils::kKaon] = 2.1;
  221. for ( Int_t i = 19; i < 22; i++ ) XMAX[i][MpdPidUtils::kKaon] = 2.;
  222. for ( Int_t i = 22; i < 25; i++ ) XMAX[i][MpdPidUtils::kKaon] = 1.9;
  223. for ( Int_t i = 25; i < 38; i++ ) XMAX[i][MpdPidUtils::kKaon] = 1.8;
  224. for ( Int_t i = 38; i < MpdPidUtils::nQAHists; i++ ) XMAX[i][MpdPidUtils::kKaon] = 1.7;
  225. } else if ( fTrackingState == MpdPidUtils::kCF ) {
  226. MaxEnLossLoc = 200000.;
  227. XMIN[0][MpdPidUtils::kKaon] = 10000.; XMIN[1][MpdPidUtils::kKaon] = 12000.; XMIN[2][MpdPidUtils::kKaon] = 8000.; XMIN[3][MpdPidUtils::kKaon] = 5000.; XMIN[4][MpdPidUtils::kKaon] = 3000.;
  228. for ( Int_t i = 5; i < 14; i++ ) XMIN[i][MpdPidUtils::kKaon] = 2000.; for ( Int_t i = 14; i < MpdPidUtils::nQAHists; i++ ) XMIN[i][MpdPidUtils::kKaon] = 1000.;
  229. XMAX[0][MpdPidUtils::kKaon] = 100000.; XMAX[1][MpdPidUtils::kKaon] = 60000.; XMAX[2][MpdPidUtils::kKaon] = 40000.; XMAX[3][MpdPidUtils::kKaon] = 22000.; XMAX[4][MpdPidUtils::kKaon] = 17000.;
  230. XMAX[5][MpdPidUtils::kKaon] = 13000.; XMAX[6][MpdPidUtils::kKaon] = 13000.; XMAX[7][MpdPidUtils::kKaon] = 10000.; XMAX[8][MpdPidUtils::kKaon] = 10000.; XMAX[9][MpdPidUtils::kKaon] = 9000.;
  231. XMAX[10][MpdPidUtils::kKaon] = 9000.; XMAX[11][MpdPidUtils::kKaon] = 7000.; XMAX[12][MpdPidUtils::kKaon] = 7000.; XMAX[13][MpdPidUtils::kKaon] = 6000.; XMAX[14][MpdPidUtils::kKaon] = 6000.;
  232. for ( Int_t i = 15; i < MpdPidUtils::nQAHists; i++ ) XMAX[i][MpdPidUtils::kKaon] = 5000.;
  233. } else {
  234. MaxEnLossLoc = 100.e-06;
  235. XMIN[0][MpdPidUtils::kKaon] = 10.e-06; XMIN[1][MpdPidUtils::kKaon] = 10.e-06; XMIN[2][MpdPidUtils::kKaon] = 3.e-06; XMIN[3][MpdPidUtils::kKaon] = 3.e-06; XMIN[4][MpdPidUtils::kKaon] = 3.e-06;
  236. XMIN[5][MpdPidUtils::kKaon] = 2.e-06; XMIN[6][MpdPidUtils::kKaon] = 2.e-06; XMIN[7][MpdPidUtils::kKaon] = 1.5e-06; XMIN[8][MpdPidUtils::kKaon] = 1.5e-06; XMIN[9][MpdPidUtils::kKaon] = 1.5e-06;
  237. XMIN[10][MpdPidUtils::kKaon] = 1.5e-06; for ( Int_t i = 11; i < MpdPidUtils::nQAHists; i++ ) XMIN[i][MpdPidUtils::kKaon] = 1.e-06;
  238. XMAX[0][MpdPidUtils::kKaon] = 100.e-06; XMAX[1][MpdPidUtils::kKaon] = 50.e-06; XMAX[2][MpdPidUtils::kKaon] = 25.e-06; XMAX[3][MpdPidUtils::kKaon] = 12.e-06; XMAX[4][MpdPidUtils::kKaon] = 8.e-06;
  239. XMAX[5][MpdPidUtils::kKaon] = 7.e-06; XMAX[6][MpdPidUtils::kKaon] = 5.e-06; XMAX[7][MpdPidUtils::kKaon] = 4.5e-06; XMAX[8][MpdPidUtils::kKaon] = 4.e-06; XMAX[9][MpdPidUtils::kKaon] = 3.5e-06;
  240. XMAX[10][MpdPidUtils::kKaon] = 3.e-06; XMAX[11][MpdPidUtils::kKaon] = 3.e-06; XMAX[12][MpdPidUtils::kKaon] = 3.e-06; XMAX[13][MpdPidUtils::kKaon] = 3.e-06; XMAX[14][MpdPidUtils::kKaon] = 3.e-06;
  241. for ( Int_t i = 15; i < MpdPidUtils::nQAHists; i++ ) XMAX[i][MpdPidUtils::kKaon] = 2.5e-06;
  242. }
  243. MpdPidUtils::dEdXStruct_t SKa;
  244. for ( Int_t iHist = 0; iHist < MpdPidUtils::nQAHists; iHist++ ) {
  245. histName = "Ka_" + TString::Itoa(iHist, 10);
  246. SKa.dEdXPart[iHist] = new TH1D(histName, "", nbins, XMIN[iHist][MpdPidUtils::kKaon], XMAX[iHist][MpdPidUtils::kKaon]);
  247. }
  248. SKa.BetheBlochHist = new TH2D("hKaBB_QA", "", 300, 0.05, 3.0, 300, MinEnLoss, MaxEnLossLoc);
  249. SKa.ibeg = 0; SKa.iend = 39;
  250. fEnLossMap.insert( pair<MpdPidUtils::ePartType,MpdPidUtils::dEdXStruct_t>(MpdPidUtils::kKaon,SKa) );
  251. /// m-squared
  252. TH2D *mass2Ka = new TH2D("mass2Ka", "", nbins, 0., 4., nbins, -0.5, 1.5);
  253. fMSquaredMap.insert( pair<MpdPidUtils::ePartType,TH2D*>(MpdPidUtils::kKaon,mass2Ka) );
  254. /// Abundance
  255. vecTH1Dptrs vecYield;
  256. TH1D *hYield_kapos = new TH1D("hYield_kapos","",600,0.,3.); vecYield.push_back(hYield_kapos);
  257. TH1D *hYield_kaneg = new TH1D("hYield_kaneg","",600,0.,3.); vecYield.push_back(hYield_kaneg);
  258. fAbundanceMap.insert( pair<MpdPidUtils::ePartType,vecTH1Dptrs>(MpdPidUtils::kKaon,vecYield) );
  259. /// Efficiency & Contamination
  260. TH1D *hEC_pos[4], *hEC_neg[4];
  261. vecTH1Dptrs vecEC_pos, vecEC_neg;
  262. std::vector<vecTH1Dptrs> vecvecEC;
  263. for ( Int_t iHist = 0; iHist < 4; iHist++ ) {
  264. hEC_pos[iHist] = new TH1D( ECHistNames[iHist] + "_poschar_" + MpdPidUtils::cParticleShortName[MpdPidUtils::kKaon],"", 50, PMIN, PMAX);
  265. hEC_pos[iHist]->Sumw2();
  266. vecEC_pos.push_back( hEC_pos[iHist] );
  267. hEC_neg[iHist] = new TH1D( ECHistNames[iHist] + "_negchar_" + MpdPidUtils::cParticleShortName[MpdPidUtils::kKaon], "", 50, PMIN, PMAX);
  268. hEC_neg[iHist]->Sumw2();
  269. vecEC_neg.push_back( hEC_neg[iHist] );
  270. }
  271. vecvecEC.push_back( vecEC_pos );
  272. vecvecEC.push_back( vecEC_neg );
  273. fEffContMap.insert( pair<MpdPidUtils::ePartType,std::vector<vecTH1Dptrs>>(MpdPidUtils::kKaon,vecvecEC) );
  274. }
  275. if (Particles.Contains("pr")) {
  276. /// dE/dx
  277. if ( fTrackingState == MpdPidUtils::kCFHM ) {
  278. MaxEnLossLoc = 90.;
  279. XMIN[0][MpdPidUtils::kProton] = 30.; XMIN[1][MpdPidUtils::kProton] = 12.; XMIN[2][MpdPidUtils::kProton] = 10.; XMIN[3][MpdPidUtils::kProton] = 10.; XMIN[4][MpdPidUtils::kProton] = 6.;
  280. XMIN[5][MpdPidUtils::kProton] = 6.; XMIN[6][MpdPidUtils::kProton] = 5.; XMIN[7][MpdPidUtils::kProton] = 4.; XMIN[8][MpdPidUtils::kProton] = 3.5; XMIN[9][MpdPidUtils::kProton] = 3.;
  281. XMIN[10][MpdPidUtils::kProton] = 2.7; XMIN[11][MpdPidUtils::kProton] = 2.5; XMIN[12][MpdPidUtils::kProton] = 2.2; XMIN[13][MpdPidUtils::kProton] = 2.; XMIN[14][MpdPidUtils::kProton] = 1.8;
  282. XMIN[15][MpdPidUtils::kProton] = 1.7; XMIN[16][MpdPidUtils::kProton] = 1.6; XMIN[17][MpdPidUtils::kProton] = 1.5; XMIN[18][MpdPidUtils::kProton] = 1.4; XMIN[19][MpdPidUtils::kProton] = 1.4;
  283. for ( Int_t i = 20; i < 23; i++ ) XMIN[i][MpdPidUtils::kProton] = 1.3;
  284. for ( Int_t i = 23; i < 26; i++ ) XMIN[i][MpdPidUtils::kProton] = 1.2;
  285. for ( Int_t i = 26; i < 31; i++ ) XMIN[i][MpdPidUtils::kProton] = 1.1;
  286. for ( Int_t i = 31; i < 35; i++ ) XMIN[i][MpdPidUtils::kProton] = 1.0;
  287. for ( Int_t i = 35; i < MpdPidUtils::nQAHists; i++ ) XMIN[i][MpdPidUtils::kProton] = 0.9;
  288. XMAX[0][MpdPidUtils::kProton] = 90.; XMAX[1][MpdPidUtils::kProton] = 80.; XMAX[2][MpdPidUtils::kProton] = 60.; XMAX[3][MpdPidUtils::kProton] = 37.; XMAX[4][MpdPidUtils::kProton] = 25.;
  289. XMAX[5][MpdPidUtils::kProton] = 16.; XMAX[6][MpdPidUtils::kProton] = 12.; XMAX[7][MpdPidUtils::kProton] = 10.; XMAX[8][MpdPidUtils::kProton] = 8.; XMAX[9][MpdPidUtils::kProton] = 7.;
  290. XMAX[10][MpdPidUtils::kProton] = 5.8; XMAX[11][MpdPidUtils::kProton] = 5.2; XMAX[12][MpdPidUtils::kProton] = 4.6; XMAX[13][MpdPidUtils::kProton] = 4.2; XMAX[14][MpdPidUtils::kProton] = 3.8;
  291. XMAX[15][MpdPidUtils::kProton] = 3.5; XMAX[16][MpdPidUtils::kProton] = 3.4; XMAX[17][MpdPidUtils::kProton] = 3.3; XMAX[18][MpdPidUtils::kProton] = 3.2; XMAX[19][MpdPidUtils::kProton] = 3.0;
  292. XMAX[20][MpdPidUtils::kProton] = 2.8; XMAX[21][MpdPidUtils::kProton] = 2.7; XMAX[22][MpdPidUtils::kProton] = 2.6; XMAX[23][MpdPidUtils::kProton] = 2.5; XMAX[24][MpdPidUtils::kProton] = 2.5;
  293. XMAX[25][MpdPidUtils::kProton] = 2.5; XMAX[26][MpdPidUtils::kProton] = 2.4; XMAX[27][MpdPidUtils::kProton] = 2.3; XMAX[28][MpdPidUtils::kProton] = 2.3; XMAX[29][MpdPidUtils::kProton] = 2.2;
  294. XMAX[30][MpdPidUtils::kProton] = 2.2; XMAX[31][MpdPidUtils::kProton] = 2.1; XMAX[32][MpdPidUtils::kProton] = 2.1;
  295. for ( Int_t i = 33; i < MpdPidUtils::nQAHists; i++ ) XMAX[i][MpdPidUtils::kProton] = 2.0;
  296. } else if ( fTrackingState == MpdPidUtils::kCF ) {
  297. MaxEnLossLoc = 400000.;
  298. XMIN[0][MpdPidUtils::kProton] = 100000.; XMIN[1][MpdPidUtils::kProton] = 50000.; XMIN[2][MpdPidUtils::kProton] = 20000.; XMIN[3][MpdPidUtils::kProton] = 15000.; XMIN[4][MpdPidUtils::kProton] = 10000.;
  299. XMIN[5][MpdPidUtils::kProton] = 10000.; XMIN[6][MpdPidUtils::kProton] = 8000.; XMIN[7][MpdPidUtils::kProton] = 8000.; XMIN[8][MpdPidUtils::kProton] = 8000.;
  300. for ( Int_t i = 9; i < 14; i++ ) XMIN[i][MpdPidUtils::kProton] = 5000.; XMIN[14][MpdPidUtils::kProton] = 2000.;
  301. for ( Int_t i = 15; i < 23; i++ ) XMIN[i][MpdPidUtils::kProton] = 3000.; for ( Int_t i = 23; i < MpdPidUtils::nQAHists; i++ ) XMIN[i][MpdPidUtils::kProton] = 2000.;
  302. XMAX[0][MpdPidUtils::kProton] = 800000.; XMAX[1][MpdPidUtils::kProton] = 250000.; XMAX[2][MpdPidUtils::kProton] = 180000.; XMAX[3][MpdPidUtils::kProton] = 90000.; XMAX[4][MpdPidUtils::kProton] = 50000.;
  303. XMAX[5][MpdPidUtils::kProton] = 30000.; XMAX[6][MpdPidUtils::kProton] = 25000.; XMAX[7][MpdPidUtils::kProton] = 19000.; XMAX[8][MpdPidUtils::kProton] = 17000.; XMAX[9][MpdPidUtils::kProton] = 15000.;
  304. XMAX[10][MpdPidUtils::kProton] = 15000.; XMAX[11][MpdPidUtils::kProton] = 15000.; XMAX[12][MpdPidUtils::kProton] = 12000.; XMAX[13][MpdPidUtils::kProton] = 12000.; XMAX[14][MpdPidUtils::kProton] = 12000.;
  305. XMAX[15][MpdPidUtils::kProton] = 9000.; XMAX[16][MpdPidUtils::kProton] = 9000.; XMAX[17][MpdPidUtils::kProton] = 8000.; XMAX[18][MpdPidUtils::kProton] = 8000.; XMAX[19][MpdPidUtils::kProton] = 8000.;
  306. XMAX[20][MpdPidUtils::kProton] = 7000.; XMAX[21][MpdPidUtils::kProton] = 7000.; XMAX[22][MpdPidUtils::kProton] = 7000.; for ( Int_t i = 23; i < MpdPidUtils::nQAHists; i++ ) XMAX[i][MpdPidUtils::kProton] = 6000.;
  307. } else {
  308. MaxEnLossLoc = 180.e-06;
  309. XMIN[0][MpdPidUtils::kProton] = 75.e-06; XMIN[1][MpdPidUtils::kProton] = 40.e-06; XMIN[2][MpdPidUtils::kProton] = 20.e-06; XMIN[3][MpdPidUtils::kProton] = 10.e-06; XMIN[4][MpdPidUtils::kProton] = 10.e-06;
  310. XMIN[5][MpdPidUtils::kProton] = 8.e-06; XMIN[6][MpdPidUtils::kProton] = 6.e-06; XMIN[7][MpdPidUtils::kProton] = 5.e-06; XMIN[8][MpdPidUtils::kProton] = 4.e-06; XMIN[9][MpdPidUtils::kProton] = 3.5e-06;
  311. XMIN[10][MpdPidUtils::kProton] = 3.e-06; XMIN[11][MpdPidUtils::kProton] = 3.e-06; XMIN[12][MpdPidUtils::kProton] = 2.5e-06; for ( Int_t i = 13; i < 17; i++ ) XMIN[i][MpdPidUtils::kProton] = 2.e-06;
  312. for ( Int_t i = 17; i < 23; i++ ) XMIN[i][MpdPidUtils::kProton] = 1.5e-06; for ( Int_t i = 23; i < MpdPidUtils::nQAHists; i++ ) XMIN[i][MpdPidUtils::kProton] = 1.e-06;
  313. XMAX[0][MpdPidUtils::kProton] = 180.e-06; XMAX[1][MpdPidUtils::kProton] = 100.e-06; XMAX[2][MpdPidUtils::kProton] = 70.e-06; XMAX[3][MpdPidUtils::kProton] = 35.e-06; XMAX[4][MpdPidUtils::kProton] = 25.e-06;
  314. XMAX[5][MpdPidUtils::kProton] = 17.e-06; XMAX[6][MpdPidUtils::kProton] = 13.e-06; XMAX[7][MpdPidUtils::kProton] = 11.e-06; XMAX[8][MpdPidUtils::kProton] = 9.e-06; XMAX[9][MpdPidUtils::kProton] = 7.5e-06;
  315. XMAX[10][MpdPidUtils::kProton] = 6.5e-06; XMAX[11][MpdPidUtils::kProton] = 6.e-06; XMAX[12][MpdPidUtils::kProton] = 5.e-06; XMAX[13][MpdPidUtils::kProton] = 5.e-06; XMAX[14][MpdPidUtils::kProton] = 4.5e-06;
  316. XMAX[15][MpdPidUtils::kProton] = 4.e-06; XMAX[16][MpdPidUtils::kProton] = 4.e-06; for ( Int_t i = 17; i < 21; i++ ) XMAX[i][MpdPidUtils::kProton] = 3.5e-06;
  317. for ( Int_t i = 21; i < 31; i++ ) XMAX[i][MpdPidUtils::kProton] = 3.e-06; for ( Int_t i = 31; i < MpdPidUtils::nQAHists; i++ ) XMAX[i][MpdPidUtils::kProton] = 2.5e-06;
  318. }
  319. MpdPidUtils::dEdXStruct_t SPr;
  320. for ( Int_t iHist = 0; iHist < MpdPidUtils::nQAHists; iHist++ ) {
  321. histName = "Pr_" + TString::Itoa(iHist, 10);
  322. SPr.dEdXPart[iHist] = new TH1D(histName, "", nbins, XMIN[iHist][MpdPidUtils::kProton], XMAX[iHist][MpdPidUtils::kProton]);
  323. }
  324. SPr.BetheBlochHist = new TH2D("hPrBB_QA", "", 300, 0.05, 3.0, 300, MinEnLoss, MaxEnLossLoc);
  325. SPr.ibeg = 0; SPr.iend = 39;
  326. fEnLossMap.insert( pair<MpdPidUtils::ePartType,MpdPidUtils::dEdXStruct_t>(MpdPidUtils::kProton,SPr) );
  327. /// m-squared
  328. TH2D *mass2Pr = new TH2D("mass2Pr", "", nbins, 0., 4., nbins, -0.5, 1.5);
  329. fMSquaredMap.insert( pair<MpdPidUtils::ePartType,TH2D*>(MpdPidUtils::kProton,mass2Pr) );
  330. /// Abundance
  331. vecTH1Dptrs vecYield;
  332. TH1D *hYield_prpos = new TH1D("hYield_prpos","",600,0.,3.); vecYield.push_back(hYield_prpos);
  333. TH1D *hYield_prneg = new TH1D("hYield_prneg","",600,0.,3.); vecYield.push_back(hYield_prneg);
  334. fAbundanceMap.insert( pair<MpdPidUtils::ePartType,vecTH1Dptrs>(MpdPidUtils::kProton,vecYield) );
  335. /// Efficiency & Contamination
  336. TH1D *hEC_pos[4], *hEC_neg[4];
  337. vecTH1Dptrs vecEC_pos, vecEC_neg;
  338. std::vector<vecTH1Dptrs> vecvecEC;
  339. for ( Int_t iHist = 0; iHist < 4; iHist++ ) {
  340. hEC_pos[iHist] = new TH1D( ECHistNames[iHist] + "_poschar_" + MpdPidUtils::cParticleShortName[MpdPidUtils::kProton],"", 50, PMIN, PMAX);
  341. hEC_pos[iHist]->Sumw2();
  342. vecEC_pos.push_back( hEC_pos[iHist] );
  343. hEC_neg[iHist] = new TH1D( ECHistNames[iHist] + "_negchar_" + MpdPidUtils::cParticleShortName[MpdPidUtils::kProton], "", 50, PMIN, PMAX);
  344. hEC_neg[iHist]->Sumw2();
  345. vecEC_neg.push_back( hEC_neg[iHist] );
  346. }
  347. vecvecEC.push_back( vecEC_pos );
  348. vecvecEC.push_back( vecEC_neg );
  349. fEffContMap.insert( pair<MpdPidUtils::ePartType,std::vector<vecTH1Dptrs>>(MpdPidUtils::kProton,vecvecEC) );
  350. }
  351. if (Particles.Contains("de")) {
  352. /// dE/dx
  353. if ( fTrackingState == MpdPidUtils::kCFHM ) {
  354. MaxEnLossLoc = 90.;
  355. XMIN[0][MpdPidUtils::kDeuteron] = 20.; XMIN[1][MpdPidUtils::kDeuteron] = 20.; XMIN[2][MpdPidUtils::kDeuteron] = 15.; XMIN[3][MpdPidUtils::kDeuteron] = 15.; XMIN[4][MpdPidUtils::kDeuteron] = 10.;
  356. XMIN[5][MpdPidUtils::kDeuteron] = 10.; XMIN[6][MpdPidUtils::kDeuteron] = 8.; XMIN[7][MpdPidUtils::kDeuteron] = 7.; XMIN[8][MpdPidUtils::kDeuteron] = 6.; XMIN[9][MpdPidUtils::kDeuteron] = 5.;
  357. XMIN[10][MpdPidUtils::kDeuteron] = 5.; XMIN[11][MpdPidUtils::kDeuteron] = 5.; XMIN[12][MpdPidUtils::kDeuteron] = 4.; XMIN[13][MpdPidUtils::kDeuteron] = 3.; XMIN[14][MpdPidUtils::kDeuteron] = 3.;
  358. XMIN[15][MpdPidUtils::kDeuteron] = 2.; XMIN[16][MpdPidUtils::kDeuteron] = 2.; XMIN[17][MpdPidUtils::kDeuteron] = 1.5; XMIN[18][MpdPidUtils::kDeuteron] = 1.5; XMIN[19][MpdPidUtils::kDeuteron] = 1.5;
  359. for ( Int_t i = 20; i < 30; i++ ) XMIN[i][MpdPidUtils::kDeuteron] = 1.0;
  360. for ( Int_t i = 30; i < 37; i++ ) XMIN[i][MpdPidUtils::kDeuteron] = 0.5;
  361. XMAX[0][MpdPidUtils::kDeuteron] = 90.; XMAX[1][MpdPidUtils::kDeuteron] = 70.; XMAX[2][MpdPidUtils::kDeuteron] = 55.; XMAX[3][MpdPidUtils::kDeuteron] = 45.; XMAX[4][MpdPidUtils::kDeuteron] = 40.;
  362. XMAX[5][MpdPidUtils::kDeuteron] = 32.; XMAX[6][MpdPidUtils::kDeuteron] = 26.; XMAX[7][MpdPidUtils::kDeuteron] = 22.; XMAX[8][MpdPidUtils::kDeuteron] = 19.; XMAX[9][MpdPidUtils::kDeuteron] = 16.;
  363. XMAX[10][MpdPidUtils::kDeuteron] = 14.; XMAX[11][MpdPidUtils::kDeuteron] = 12.; XMAX[12][MpdPidUtils::kDeuteron] = 11.; XMAX[13][MpdPidUtils::kDeuteron] = 10.; XMAX[14][MpdPidUtils::kDeuteron] = 9.;
  364. XMAX[15][MpdPidUtils::kDeuteron] = 9.; XMAX[16][MpdPidUtils::kDeuteron] = 8.; XMAX[17][MpdPidUtils::kDeuteron] = 7.5; XMAX[18][MpdPidUtils::kDeuteron] = 7.5; XMAX[19][MpdPidUtils::kDeuteron] = 7.;
  365. XMAX[20][MpdPidUtils::kDeuteron] = 6.5; XMAX[21][MpdPidUtils::kDeuteron] = 6.5; XMAX[22][MpdPidUtils::kDeuteron] = 6.; XMAX[23][MpdPidUtils::kDeuteron] = 5.5; XMAX[24][MpdPidUtils::kDeuteron] = 5.5;
  366. XMAX[25][MpdPidUtils::kDeuteron] = 5.; XMAX[26][MpdPidUtils::kDeuteron] = 5.; XMAX[27][MpdPidUtils::kDeuteron] = 5.; XMAX[28][MpdPidUtils::kDeuteron] = 4.5;
  367. for ( Int_t i = 29; i < 35; i++ ) XMAX[i][MpdPidUtils::kDeuteron] = 4.;
  368. XMAX[35][MpdPidUtils::kDeuteron] = 3.5; XMAX[36][MpdPidUtils::kDeuteron] = 3.5;
  369. } else if ( fTrackingState == MpdPidUtils::kCF ) {
  370. MaxEnLossLoc = 500000.;
  371. XMIN[0][MpdPidUtils::kDeuteron] = 50000.; XMIN[1][MpdPidUtils::kDeuteron] = 20000.; XMIN[2][MpdPidUtils::kDeuteron] = 10000.; XMIN[3][MpdPidUtils::kDeuteron] = 8000.;
  372. XMIN[4][MpdPidUtils::kDeuteron] = 8000.; XMIN[5][MpdPidUtils::kDeuteron] = 5000.; for ( Int_t i = 6; i < (MpdPidUtils::nQAHists - 2); i++ ) XMIN[i][MpdPidUtils::kDeuteron] = 2000.;
  373. XMAX[0][MpdPidUtils::kDeuteron] = 500000.; XMAX[1][MpdPidUtils::kDeuteron] = 400000.; XMAX[2][MpdPidUtils::kDeuteron] = 250000.; XMAX[3][MpdPidUtils::kDeuteron] = 200000.;
  374. XMAX[4][MpdPidUtils::kDeuteron] = 120000.; XMAX[5][MpdPidUtils::kDeuteron] = 80000.; XMAX[6][MpdPidUtils::kDeuteron] = 70000.; XMAX[7][MpdPidUtils::kDeuteron] = 60000.; XMAX[8][MpdPidUtils::kDeuteron] = 50000.;
  375. XMAX[9][MpdPidUtils::kDeuteron] = 40000.; XMAX[10][MpdPidUtils::kDeuteron] = 35000.; XMAX[11][MpdPidUtils::kDeuteron] = 30000.; XMAX[12][MpdPidUtils::kDeuteron] = 25000.; XMAX[13][MpdPidUtils::kDeuteron] = 20000.;
  376. XMAX[14][MpdPidUtils::kDeuteron] = 20000.; XMAX[15][MpdPidUtils::kDeuteron] = 18000.; XMAX[16][MpdPidUtils::kDeuteron] = 15000.; XMAX[17][MpdPidUtils::kDeuteron] = 15000.; XMAX[18][MpdPidUtils::kDeuteron] = 15000.;
  377. XMAX[19][MpdPidUtils::kDeuteron] = 13000.; XMAX[20][MpdPidUtils::kDeuteron] = 13000.; XMAX[21][MpdPidUtils::kDeuteron] = 13000.; XMAX[22][MpdPidUtils::kDeuteron] = 13000.; XMAX[23][MpdPidUtils::kDeuteron] = 11000.;
  378. XMAX[24][MpdPidUtils::kDeuteron] = 11000.; XMAX[25][MpdPidUtils::kDeuteron] = 10000.; XMAX[26][MpdPidUtils::kDeuteron] = 10000.; XMAX[27][MpdPidUtils::kDeuteron] = 10000.; XMAX[28][MpdPidUtils::kDeuteron] = 9000.;
  379. XMAX[29][MpdPidUtils::kDeuteron] = 8000.; XMAX[30][MpdPidUtils::kDeuteron] = 8000.; XMAX[31][MpdPidUtils::kDeuteron] = 7000.; XMAX[32][MpdPidUtils::kDeuteron] = 7000.; XMAX[33][MpdPidUtils::kDeuteron] = 7000.;
  380. XMAX[34][MpdPidUtils::kDeuteron] = 7000.; XMAX[35][MpdPidUtils::kDeuteron] = 6000.; XMAX[36][MpdPidUtils::kDeuteron] = 6000.; XMAX[37][MpdPidUtils::kDeuteron] = 6000.;
  381. } else {
  382. MaxEnLossLoc = 200.e-06;
  383. for ( Int_t i = 0; i < (MpdPidUtils::nQAHists - 2); i++ ) XMIN[i][MpdPidUtils::kDeuteron] = 1.e-06;
  384. XMAX[0][MpdPidUtils::kDeuteron] = 200.e-06; XMAX[1][MpdPidUtils::kDeuteron] = 200.e-06; XMAX[2][MpdPidUtils::kDeuteron] = 150.e-06; XMAX[3][MpdPidUtils::kDeuteron] = 140.e-06; XMAX[4][MpdPidUtils::kDeuteron] = 100.e-06;
  385. XMAX[5][MpdPidUtils::kDeuteron] = 80.e-06; XMAX[6][MpdPidUtils::kDeuteron] = 50.e-06; XMAX[7][MpdPidUtils::kDeuteron] = 40.e-06; XMAX[8][MpdPidUtils::kDeuteron] = 30.e-06; XMAX[9][MpdPidUtils::kDeuteron] = 20.e-06;
  386. XMAX[10][MpdPidUtils::kDeuteron] = 10.e-06; XMAX[11][MpdPidUtils::kDeuteron] = 8.e-06; for ( Int_t i = 12; i < (MpdPidUtils::nQAHists - 2); i++ ) XMAX[i][MpdPidUtils::kDeuteron] = 5.e-06;
  387. }
  388. MpdPidUtils::dEdXStruct_t SDe;
  389. for ( Int_t iHist = 0; iHist < MpdPidUtils::nQAHists; iHist++ ) {
  390. histName = "De_" + TString::Itoa(iHist, 10);
  391. SDe.dEdXPart[iHist] = new TH1D(histName, "", nbins, XMIN[iHist][MpdPidUtils::kDeuteron], XMAX[iHist][MpdPidUtils::kDeuteron]);
  392. }
  393. SDe.BetheBlochHist = new TH2D("hDeBB_QA", "", 300, 0.05, 3.0, 300, MinEnLoss, MaxEnLossLoc);
  394. SDe.ibeg = ( fTrackingState == MpdPidUtils::kCFHM ) ? 3 : 2; SDe.iend = 39;
  395. fEnLossMap.insert( pair<MpdPidUtils::ePartType,MpdPidUtils::dEdXStruct_t>(MpdPidUtils::kDeuteron,SDe) );
  396. /// m-squared
  397. TH2D *mass2De = new TH2D("mass2De", "", nbins, 0., 4., nbins, -0.5, 10.);
  398. fMSquaredMap.insert( pair<MpdPidUtils::ePartType,TH2D*>(MpdPidUtils::kDeuteron,mass2De) );
  399. /// Abundance
  400. vecTH1Dptrs vecYield;
  401. TH1D *hYield_de = new TH1D("hYield_de","",300,0.,3.); vecYield.push_back(hYield_de);
  402. fAbundanceMap.insert( pair<MpdPidUtils::ePartType,vecTH1Dptrs>(MpdPidUtils::kDeuteron,vecYield) );
  403. /// Efficiency & Contamination
  404. TH1D *hEC[4];
  405. vecTH1Dptrs vecEC;
  406. std::vector<vecTH1Dptrs> vecvecEC;
  407. for ( Int_t iHist = 0; iHist < 4; iHist++ ) {
  408. hEC[iHist] = new TH1D( ECHistNames[iHist] + "_" + MpdPidUtils::cParticleShortName[MpdPidUtils::kDeuteron],"", 50, PMIN, PMAX);
  409. hEC[iHist]->Sumw2();
  410. vecEC.push_back( hEC[iHist] );
  411. }
  412. vecvecEC.push_back( vecEC );
  413. fEffContMap.insert( pair<MpdPidUtils::ePartType,std::vector<vecTH1Dptrs>>(MpdPidUtils::kDeuteron,vecvecEC) );
  414. }
  415. if (Particles.Contains("tr")) {
  416. /// dE/dx
  417. if ( fTrackingState == MpdPidUtils::kCFHM ) {
  418. MaxEnLossLoc = 120.;
  419. for ( Int_t i = 0; i < 4; i++ ) XMIN[i][MpdPidUtils::kTriton] = 20.;
  420. for ( Int_t i = 4; i < 8; i++ ) XMIN[i][MpdPidUtils::kTriton] = 15.;
  421. XMIN[8][MpdPidUtils::kTriton] = 10.; XMIN[9][MpdPidUtils::kTriton] = 10.; XMIN[10][MpdPidUtils::kTriton] = 8.; XMIN[11][MpdPidUtils::kTriton] = 7.; XMIN[12][MpdPidUtils::kTriton] = 6.;
  422. XMIN[13][MpdPidUtils::kTriton] = 6.; XMIN[14][MpdPidUtils::kTriton] = 5.; XMIN[15][MpdPidUtils::kTriton] = 5.; XMIN[16][MpdPidUtils::kTriton] = 5.; XMIN[17][MpdPidUtils::kTriton] = 5.;
  423. XMIN[18][MpdPidUtils::kTriton] = 4.; XMIN[19][MpdPidUtils::kTriton] = 4.; XMIN[20][MpdPidUtils::kTriton] = 4.; XMIN[21][MpdPidUtils::kTriton] = 3.; XMIN[22][MpdPidUtils::kTriton] = 3.;
  424. XMIN[23][MpdPidUtils::kTriton] = 3.; XMIN[24][MpdPidUtils::kTriton] = 2.5; XMIN[25][MpdPidUtils::kTriton] = 2.5; XMIN[26][MpdPidUtils::kTriton] = 2.5; XMIN[27][MpdPidUtils::kTriton] = 2.5;
  425. for ( Int_t i = 28; i < 33; i++ ) XMIN[i][MpdPidUtils::kTriton] = 2.;
  426. for ( Int_t i = 33; i < 37; i++ ) XMIN[i][MpdPidUtils::kTriton] = 1.5;
  427. XMAX[0][MpdPidUtils::kTriton] = 120.; XMAX[1][MpdPidUtils::kTriton] = 120.; XMAX[2][MpdPidUtils::kTriton] = 100.; XMAX[3][MpdPidUtils::kTriton] = 80.; XMAX[4][MpdPidUtils::kTriton] = 70.;
  428. XMAX[5][MpdPidUtils::kTriton] = 60.; XMAX[6][MpdPidUtils::kTriton] = 50.; XMAX[7][MpdPidUtils::kTriton] = 45.; XMAX[8][MpdPidUtils::kTriton] = 40.; XMAX[9][MpdPidUtils::kTriton] = 35.;
  429. XMAX[10][MpdPidUtils::kTriton] = 30.; XMAX[11][MpdPidUtils::kTriton] = 26.; XMAX[12][MpdPidUtils::kTriton] = 24.; XMAX[13][MpdPidUtils::kTriton] = 21.; XMAX[14][MpdPidUtils::kTriton] = 19.;
  430. XMAX[15][MpdPidUtils::kTriton] = 17.; XMAX[16][MpdPidUtils::kTriton] = 15.; XMAX[17][MpdPidUtils::kTriton] = 13.; XMAX[18][MpdPidUtils::kTriton] = 13.; XMAX[19][MpdPidUtils::kTriton] = 12.;
  431. XMAX[20][MpdPidUtils::kTriton] = 11.; XMAX[21][MpdPidUtils::kTriton] = 11.; XMAX[22][MpdPidUtils::kTriton] = 10.; XMAX[23][MpdPidUtils::kTriton] = 9.; XMAX[24][MpdPidUtils::kTriton] = 8.5;
  432. XMAX[25][MpdPidUtils::kTriton] = 8.; XMAX[26][MpdPidUtils::kTriton] = 7.5; XMAX[27][MpdPidUtils::kTriton] = 7.; XMAX[28][MpdPidUtils::kTriton] = 7.; XMAX[29][MpdPidUtils::kTriton] = 6.5;
  433. XMAX[30][MpdPidUtils::kTriton] = 6.; XMAX[31][MpdPidUtils::kTriton] = 5.5; XMAX[32][MpdPidUtils::kTriton] = 5.; XMAX[33][MpdPidUtils::kTriton] = 4.5; XMAX[34][MpdPidUtils::kTriton] = 4.5;
  434. XMAX[35][MpdPidUtils::kTriton] = 4.; XMAX[36][MpdPidUtils::kTriton] = 4.;
  435. } else if ( fTrackingState == MpdPidUtils::kCF ) {
  436. MaxEnLossLoc = 500000.;
  437. XMIN[0][MpdPidUtils::kTriton] = 100000.; XMIN[1][MpdPidUtils::kTriton] = 50000.; XMIN[2][MpdPidUtils::kTriton] = 40000.; XMIN[3][MpdPidUtils::kTriton] = 30000.; XMIN[4][MpdPidUtils::kTriton] = 20000.;
  438. XMIN[5][MpdPidUtils::kTriton] = 10000.; XMIN[6][MpdPidUtils::kTriton] = 8000.; XMIN[7][MpdPidUtils::kTriton] = 5000.; for ( Int_t i = 8; i < 14; i++ ) XMIN[i][MpdPidUtils::kTriton] = 2000.;
  439. for ( Int_t i = 14; i < 20; i++ ) XMIN[i][MpdPidUtils::kTriton] = 5000.; for ( Int_t i = 20; i < 26; i++ ) XMIN[i][MpdPidUtils::kTriton] = 6000.;
  440. for ( Int_t i = 26; i < 29; i++ ) XMIN[i][MpdPidUtils::kTriton] = 5000.; for ( Int_t i = 29; i < 35; i++ ) XMIN[i][MpdPidUtils::kTriton] = 4000.;
  441. XMIN[35][MpdPidUtils::kTriton] = 3500.; XMIN[36][MpdPidUtils::kTriton] = 3500.; XMAX[0][MpdPidUtils::kTriton] = 500000.; XMAX[1][MpdPidUtils::kTriton] = 500000.; XMAX[2][MpdPidUtils::kTriton] = 350000.;
  442. XMAX[3][MpdPidUtils::kTriton] = 250000.; XMAX[4][MpdPidUtils::kTriton] = 200000.; XMAX[5][MpdPidUtils::kTriton] = 160000.; XMAX[6][MpdPidUtils::kTriton] = 140000.; XMAX[7][MpdPidUtils::kTriton] = 100000.;
  443. XMAX[8][MpdPidUtils::kTriton] = 90000.; XMAX[9][MpdPidUtils::kTriton] = 80000.; XMAX[10][MpdPidUtils::kTriton] = 60000.; XMAX[11][MpdPidUtils::kTriton] = 50000.; XMAX[12][MpdPidUtils::kTriton] = 45000.;
  444. XMAX[13][MpdPidUtils::kTriton] = 40000.; XMAX[14][MpdPidUtils::kTriton] = 40000.; XMAX[15][MpdPidUtils::kTriton] = 30000.; XMAX[16][MpdPidUtils::kTriton] = 30000.; XMAX[17][MpdPidUtils::kTriton] = 30000.;
  445. XMAX[18][MpdPidUtils::kTriton] = 25000.; XMAX[19][MpdPidUtils::kTriton] = 25000.; XMAX[20][MpdPidUtils::kTriton] = 22000.; XMAX[21][MpdPidUtils::kTriton] = 22000.; XMAX[22][MpdPidUtils::kTriton] = 20000.;
  446. XMAX[23][MpdPidUtils::kTriton] = 18000.; XMAX[24][MpdPidUtils::kTriton] = 18000.; XMAX[25][MpdPidUtils::kTriton] = 16000.; XMAX[26][MpdPidUtils::kTriton] = 16000.; XMAX[27][MpdPidUtils::kTriton] = 15000.;
  447. XMAX[28][MpdPidUtils::kTriton] = 14000.; XMAX[29][MpdPidUtils::kTriton] = 14000.; XMAX[30][MpdPidUtils::kTriton] = 13000.; XMAX[31][MpdPidUtils::kTriton] = 11000.; XMAX[32][MpdPidUtils::kTriton] = 11000.;
  448. XMAX[33][MpdPidUtils::kTriton] = 10000.; XMAX[34][MpdPidUtils::kTriton] = 10000.; XMAX[35][MpdPidUtils::kTriton] = 9000.; XMAX[36][MpdPidUtils::kTriton] = 8000.;
  449. } else {
  450. MaxEnLossLoc = 200.e-06;
  451. XMIN[0][MpdPidUtils::kTriton] = 3.e-06; XMIN[1][MpdPidUtils::kTriton] = 3.e-06; XMIN[2][MpdPidUtils::kTriton] = 2.e-06; for ( Int_t i = 0; i < (MpdPidUtils::nQAHists - 3); i++ ) XMIN[i][MpdPidUtils::kTriton] = 1.e-06;
  452. XMAX[0][MpdPidUtils::kTriton] = 200.e-06; XMAX[1][MpdPidUtils::kTriton] = 160.e-06; XMAX[2][MpdPidUtils::kTriton] = 140.e-06; XMAX[3][MpdPidUtils::kTriton] = 110.e-06; XMAX[4][MpdPidUtils::kTriton] = 80.e-06;
  453. XMAX[5][MpdPidUtils::kTriton] = 60.e-06; XMAX[6][MpdPidUtils::kTriton] = 40.e-06; for ( Int_t i = 7; i < (MpdPidUtils::nQAHists - 3); i++ ) XMAX[i][MpdPidUtils::kTriton] = 20.e-06;
  454. }
  455. MpdPidUtils::dEdXStruct_t STr;
  456. for ( Int_t iHist = 0; iHist < MpdPidUtils::nQAHists; iHist++ ) {
  457. histName = "Tr_" + TString::Itoa(iHist, 10);
  458. STr.dEdXPart[iHist] = new TH1D(histName, "", nbins, XMIN[iHist][MpdPidUtils::kTriton], XMAX[iHist][MpdPidUtils::kTriton]);
  459. }
  460. STr.BetheBlochHist = new TH2D("hTrBB_QA", "", 300, 0.05, 3.0, 300, MinEnLoss, MaxEnLossLoc);
  461. STr.ibeg = 3; STr.iend = 39;
  462. fEnLossMap.insert( pair<MpdPidUtils::ePartType,MpdPidUtils::dEdXStruct_t>(MpdPidUtils::kTriton,STr) );
  463. /// m-squared
  464. TH2D *mass2Tr = new TH2D("mass2Tr", "", nbins, 0., 4., nbins, -0.5, 15.);
  465. fMSquaredMap.insert( pair<MpdPidUtils::ePartType,TH2D*>(MpdPidUtils::kTriton,mass2Tr) );
  466. /// Abundance
  467. vecTH1Dptrs vecYield;
  468. TH1D *hYield_tr = new TH1D("hYield_tr","",300,0.,3.); vecYield.push_back(hYield_tr);
  469. fAbundanceMap.insert( pair<MpdPidUtils::ePartType,vecTH1Dptrs>(MpdPidUtils::kTriton,vecYield) );
  470. /// Efficiency & Contamination
  471. TH1D *hEC[4];
  472. vecTH1Dptrs vecEC;
  473. std::vector<vecTH1Dptrs> vecvecEC;
  474. for ( Int_t iHist = 0; iHist < 4; iHist++ ) {
  475. hEC[iHist] = new TH1D( ECHistNames[iHist] + "_" + MpdPidUtils::cParticleShortName[MpdPidUtils::kTriton],"", 50, PMIN, PMAX);
  476. hEC[iHist]->Sumw2();
  477. vecEC.push_back( hEC[iHist] );
  478. }
  479. vecvecEC.push_back( vecEC );
  480. fEffContMap.insert( pair<MpdPidUtils::ePartType,std::vector<vecTH1Dptrs>>(MpdPidUtils::kTriton,vecvecEC) );
  481. }
  482. if (Particles.Contains("he3")) {
  483. /// dE/dx
  484. if ( fTrackingState == MpdPidUtils::kCFHM ) {
  485. MaxEnLossLoc = 150.;
  486. for ( Int_t i = 0; i < 5; i++ ) XMIN[i][MpdPidUtils::kHe3] = 5.;
  487. XMIN[5][MpdPidUtils::kHe3] = 10.; XMIN[6][MpdPidUtils::kHe3] = 10.; XMIN[7][MpdPidUtils::kHe3] = 15.; XMIN[8][MpdPidUtils::kHe3] = 15.;
  488. for ( Int_t i = 9; i < 15; i++ ) XMIN[i][MpdPidUtils::kHe3] = 10.;
  489. for ( Int_t i = 15; i < 19; i++ ) XMIN[i][MpdPidUtils::kHe3] = 8.;
  490. for ( Int_t i = 19; i < 26; i++ ) XMIN[i][MpdPidUtils::kHe3] = 7.;
  491. XMIN[26][MpdPidUtils::kHe3] = 6.; XMIN[27][MpdPidUtils::kHe3] = 6.;
  492. XMAX[0][MpdPidUtils::kHe3] = 150.; XMAX[1][MpdPidUtils::kHe3] = 150.; XMAX[2][MpdPidUtils::kHe3] = 150.; XMAX[3][MpdPidUtils::kHe3] = 120.; XMAX[4][MpdPidUtils::kHe3] = 100.;
  493. XMAX[5][MpdPidUtils::kHe3] = 80.; XMAX[6][MpdPidUtils::kHe3] = 70.; XMAX[7][MpdPidUtils::kHe3] = 60.; XMAX[8][MpdPidUtils::kHe3] = 55.; XMAX[9][MpdPidUtils::kHe3] = 50.;
  494. XMAX[10][MpdPidUtils::kHe3] = 45.; XMAX[11][MpdPidUtils::kHe3] = 40.; XMAX[12][MpdPidUtils::kHe3] = 35.; XMAX[13][MpdPidUtils::kHe3] = 30.; XMAX[14][MpdPidUtils::kHe3] = 30.;
  495. XMAX[15][MpdPidUtils::kHe3] = 25.; XMAX[16][MpdPidUtils::kHe3] = 23.; XMAX[17][MpdPidUtils::kHe3] = 21.; XMAX[18][MpdPidUtils::kHe3] = 20.; XMAX[19][MpdPidUtils::kHe3] = 19.;
  496. XMAX[20][MpdPidUtils::kHe3] = 18.; XMAX[21][MpdPidUtils::kHe3] = 17.; XMAX[22][MpdPidUtils::kHe3] = 16.; XMAX[23][MpdPidUtils::kHe3] = 16.; XMAX[24][MpdPidUtils::kHe3] = 16.;
  497. XMAX[25][MpdPidUtils::kHe3] = 15.; XMAX[26][MpdPidUtils::kHe3] = 15.; XMAX[27][MpdPidUtils::kHe3] = 15.;
  498. } else if ( fTrackingState == MpdPidUtils::kCF ) {
  499. MaxEnLossLoc = 1800000.;
  500. XMIN[0][MpdPidUtils::kHe3] = 100000.; XMIN[1][MpdPidUtils::kHe3] = 100000.; XMIN[2][MpdPidUtils::kHe3] = 50000.; XMIN[3][MpdPidUtils::kHe3] = 30000.;
  501. for ( Int_t i = 4; i < (MpdPidUtils::nQAHists - 9); i++ ) XMIN[i][MpdPidUtils::kHe3] = 10000.;
  502. XMAX[0][MpdPidUtils::kHe3] = 1800000.; XMAX[1][MpdPidUtils::kHe3] = 1000000.; XMAX[2][MpdPidUtils::kHe3] = 800000.; XMAX[3][MpdPidUtils::kHe3] = 500000.; XMAX[4][MpdPidUtils::kHe3] = 350000.;
  503. XMAX[5][MpdPidUtils::kHe3] = 250000.; XMAX[6][MpdPidUtils::kHe3] = 200000.; XMAX[7][MpdPidUtils::kHe3] = 180000.; XMAX[8][MpdPidUtils::kHe3] = 150000.; XMAX[9][MpdPidUtils::kHe3] = 130000.;
  504. XMAX[10][MpdPidUtils::kHe3] = 110000.; XMAX[11][MpdPidUtils::kHe3] = 90000.; XMAX[12][MpdPidUtils::kHe3] = 80000.; XMAX[13][MpdPidUtils::kHe3] = 70000.; XMAX[14][MpdPidUtils::kHe3] = 60000.;
  505. XMAX[15][MpdPidUtils::kHe3] = 50000.; XMAX[16][MpdPidUtils::kHe3] = 45000.; XMAX[17][MpdPidUtils::kHe3] = 45000.; XMAX[18][MpdPidUtils::kHe3] = 40000.; XMAX[19][MpdPidUtils::kHe3] = 40000.;
  506. XMAX[20][MpdPidUtils::kHe3] = 35000.; XMAX[21][MpdPidUtils::kHe3] = 35000.; XMAX[22][MpdPidUtils::kHe3] = 30000.; XMAX[23][MpdPidUtils::kHe3] = 30000.; XMAX[24][MpdPidUtils::kHe3] = 28000.;
  507. XMAX[25][MpdPidUtils::kHe3] = 28000.; XMAX[26][MpdPidUtils::kHe3] = 26000.; XMAX[27][MpdPidUtils::kHe3] = 26000.; XMAX[28][MpdPidUtils::kHe3] = 30000.; XMAX[29][MpdPidUtils::kHe3] = 30000.; XMAX[30][MpdPidUtils::kHe3] = 30000.;
  508. } else {
  509. MaxEnLossLoc = 450.e-06;
  510. XMIN[0][MpdPidUtils::kHe3] = 3.e-06; XMIN[1][MpdPidUtils::kHe3] = 3.e-06; XMIN[2][MpdPidUtils::kHe3] = 3.e-06; XMIN[3][MpdPidUtils::kHe3] = 2.e-06; XMIN[4][MpdPidUtils::kHe3] = 2.e-06;
  511. for ( Int_t i = 5; i < 9; i++ ) XMIN[i][MpdPidUtils::kHe3] = 1.5e-06; for ( Int_t i = 9; i < (MpdPidUtils::nQAHists - 9); i++ ) XMIN[i][MpdPidUtils::kHe3] = 1.e-06;
  512. XMAX[0][MpdPidUtils::kHe3] = 450.e-06; XMAX[1][MpdPidUtils::kHe3] = 350.e-06; XMAX[2][MpdPidUtils::kHe3] = 300.e-06; XMAX[3][MpdPidUtils::kHe3] = 250.e-06; XMAX[4][MpdPidUtils::kHe3] = 200.e-06;
  513. XMAX[5][MpdPidUtils::kHe3] = 150.e-06; XMAX[6][MpdPidUtils::kHe3] = 100.e-06; XMAX[7][MpdPidUtils::kHe3] = 80.e-06; XMAX[8][MpdPidUtils::kHe3] = 70.e-06; XMAX[9][MpdPidUtils::kHe3] = 60.e-06;
  514. XMAX[10][MpdPidUtils::kHe3] = 50.e-06; XMAX[11][MpdPidUtils::kHe3] = 40.e-06; XMAX[12][MpdPidUtils::kHe3] = 35.e-06; for ( Int_t i = 13; i < (MpdPidUtils::nQAHists - 9); i++ ) XMAX[i][MpdPidUtils::kHe3] = XMAX[11][MpdPidUtils::kHe3] = 32.e-06;
  515. }
  516. MpdPidUtils::dEdXStruct_t SHe3;
  517. for ( Int_t iHist = 0; iHist < MpdPidUtils::nQAHists; iHist++ ) {
  518. histName = "He3_" + TString::Itoa(iHist, 10);
  519. SHe3.dEdXPart[iHist] = new TH1D(histName, "", nbins, XMIN[iHist][MpdPidUtils::kHe3], XMAX[iHist][MpdPidUtils::kHe3]);
  520. }
  521. SHe3.BetheBlochHist = new TH2D("hHe3BB_QA", "", 300, 0.05, 3.0, 300, MinEnLoss, MaxEnLossLoc);
  522. SHe3.ibeg = ( fTrackingState == MpdPidUtils::kCFHM ) ? 3 : 2; SHe3.iend = ( fTrackingState == MpdPidUtils::kCFHM ) ? 30 : 32;
  523. fEnLossMap.insert( pair<MpdPidUtils::ePartType,MpdPidUtils::dEdXStruct_t>(MpdPidUtils::kHe3,SHe3) );
  524. /// m-squared
  525. TH2D *mass2He3 = new TH2D("mass2He3", "", nbins, 0., 4., nbins, -0.5, 6.);
  526. fMSquaredMap.insert( pair<MpdPidUtils::ePartType,TH2D*>(MpdPidUtils::kHe3,mass2He3) );
  527. /// Abundance
  528. vecTH1Dptrs vecYield;
  529. TH1D *hYield_he3 = new TH1D("hYield_he3","",150,0.,1.5); vecYield.push_back(hYield_he3);
  530. fAbundanceMap.insert( pair<MpdPidUtils::ePartType,vecTH1Dptrs>(MpdPidUtils::kHe3,vecYield) );
  531. /// Efficiency & Contamination
  532. TH1D *hEC[4];
  533. vecTH1Dptrs vecEC;
  534. std::vector<vecTH1Dptrs> vecvecEC;
  535. for ( Int_t iHist = 0; iHist < 4; iHist++ ) {
  536. hEC[iHist] = new TH1D( ECHistNames[iHist] + "_" + MpdPidUtils::cParticleShortName[MpdPidUtils::kHe3],"", 50, PMIN, PMAX);
  537. hEC[iHist]->Sumw2();
  538. vecEC.push_back( hEC[iHist] );
  539. }
  540. vecvecEC.push_back( vecEC );
  541. fEffContMap.insert( pair<MpdPidUtils::ePartType,std::vector<vecTH1Dptrs>>(MpdPidUtils::kHe3,vecvecEC) );
  542. }
  543. if (Particles.Contains("he4")) {
  544. /// dE/dx
  545. if ( fTrackingState == MpdPidUtils::kCFHM ) {
  546. MaxEnLossLoc = 250.;
  547. for (Int_t i = 0; i < 4; i++) XMIN[i][8] = 5.;
  548. XMIN[4][8] = 20.; XMIN[5][8] = 20.; XMIN[6][8] = 25.; XMIN[7][8] = 25.; XMIN[8][8] = 25.;
  549. XMIN[9][8] = 20.; XMIN[10][8] = 20.; XMIN[11][8] = 20.; XMIN[12][8] = 18.; XMIN[13][8] = 15.;
  550. XMIN[14][8] = 15.; XMIN[15][8] = 14.; XMIN[16][8] = 14.; XMIN[17][8] = 13.; XMIN[18][8] = 12.;
  551. XMIN[19][8] = 11.; XMIN[26][8] = 9.; XMIN[27][8] = 8.;
  552. for (Int_t i = 20; i < 26; i++) XMIN[i][8] = 10.;
  553. XMAX[0][8] = 250.; XMAX[1][8] = 200.; XMAX[2][8] = 150.; XMAX[3][8] = 140.; XMAX[4][8] = 120.;
  554. XMAX[5][8] = 100.; XMAX[6][8] = 90.; XMAX[7][8] = 80.; XMAX[8][8] = 75.; XMAX[9][8] = 70.;
  555. XMAX[10][8] = 60.; XMAX[11][8] = 55.; XMAX[12][8] = 50.; XMAX[13][8] = 45.; XMAX[14][8] = 42.;
  556. XMAX[15][8] = 40.; XMAX[16][8] = 35.; XMAX[17][8] = 33.; XMAX[18][8] = 30.; XMAX[19][8] = 28.;
  557. XMAX[20][8] = 26.; XMAX[21][8] = 25.; XMAX[22][8] = 23.; XMAX[23][8] = 21.; XMAX[24][8] = 20.;
  558. XMAX[25][8] = 19.; XMAX[26][8] = 19.; XMAX[27][8] = 22.;
  559. } else if ( fTrackingState == MpdPidUtils::kCF ) {
  560. MaxEnLossLoc = 1800000.;
  561. XMIN[0][8] = 100000.; XMIN[1][8] = 80000.; XMIN[2][8] = 50000.; XMIN[3][8] = 30000.; XMIN[4][8] = 30000.;
  562. XMIN[5][8] = 20000.; for (Int_t i = 6; i < (MpdPidUtils::nQAHists - 12); i++) XMIN[i][8] = 10000.;
  563. XMAX[0][8] = 1800000.; XMAX[1][8] = 1200000.; XMAX[2][8] = 900000.; XMAX[3][8] = 700000.; XMAX[4][8] = 500000.;
  564. XMAX[5][8] = 350000.; XMAX[6][8] = 300000.; XMAX[7][8] = 250000.; XMAX[8][8] = 200000.; XMAX[9][8] = 180000.;
  565. XMAX[10][8] = 160000.; XMAX[11][8] = 140000.; XMAX[12][8] = 130000.; XMAX[13][8] = 120000.; XMAX[14][8] = 100000.;
  566. XMAX[15][8] = 90000.; XMAX[16][8] = 80000.; XMAX[17][8] = 70000.; XMAX[18][8] = 70000.; XMAX[19][8] = 60000.;
  567. XMAX[20][8] = 55000.; XMAX[21][8] = 50000.; XMAX[22][8] = 45000.; XMAX[23][8] = 45000.; XMAX[24][8] = 40000.;
  568. XMAX[25][8] = 40000.; XMAX[26][8] = 40000.; XMAX[27][8] = 40000.;
  569. } else {
  570. MaxEnLossLoc = 350.e-06;
  571. XMIN[0][8] = 3.e-06; XMIN[1][8] = 2.e-06; for (Int_t i = 2; i < (MpdPidUtils::nQAHists - 12); i++) XMIN[i][8] = 1.e-06;
  572. XMAX[0][8] = 350.e-06; XMAX[1][8] = 300.e-06; XMAX[2][8] = 250.e-06; XMAX[3][8] = 200.e-06; XMAX[4][8] = 150.e-06;
  573. XMAX[5][8] = 100.e-06; XMAX[6][8] = 80.e-06; XMAX[7][8] = 70.e-06; XMAX[8][8] = 60.e-06; XMAX[9][8] = 50.e-06;
  574. XMAX[10][8] = 40.e-06; for (Int_t i = 11; i < (MpdPidUtils::nQAHists - 12); i++) XMAX[i][8] = 32.e-06;
  575. }
  576. MpdPidUtils::dEdXStruct_t SHe4;
  577. for ( Int_t iHist = 0; iHist < MpdPidUtils::nQAHists; iHist++) {
  578. histName = "He4_" + TString::Itoa(iHist, 10);
  579. SHe4.dEdXPart[iHist] = new TH1D(histName, "", ( iHist < 27 ) ? nbins : 50, XMIN[iHist][8], XMAX[iHist][8]);
  580. }
  581. SHe4.BetheBlochHist = new TH2D("hHe4BB_QA", "", 300, 0.05, 3.0, 300, MinEnLoss, MaxEnLossLoc);
  582. SHe4.ibeg = 3; SHe4.iend = 30;
  583. fEnLossMap.insert( pair<MpdPidUtils::ePartType,MpdPidUtils::dEdXStruct_t>(MpdPidUtils::kHe4,SHe4) );
  584. /// m-squared
  585. TH2D *mass2He4 = new TH2D("mass2He4", "", nbins, 0., 4., nbins, -0.5, 8.);
  586. fMSquaredMap.insert( pair<MpdPidUtils::ePartType,TH2D*>(MpdPidUtils::kHe4,mass2He4) );
  587. /// Abundance
  588. vecTH1Dptrs vecYield;
  589. TH1D *hYield_he4 = new TH1D("hYield_he4","",150,0.,1.5); vecYield.push_back(hYield_he4);
  590. fAbundanceMap.insert( pair<MpdPidUtils::ePartType,vecTH1Dptrs>(MpdPidUtils::kHe4,vecYield) );
  591. /// Efficiency & Contamination
  592. TH1D *hEC[4];
  593. vecTH1Dptrs vecEC;
  594. std::vector<vecTH1Dptrs> vecvecEC;
  595. for ( Int_t iHist = 0; iHist < 4; iHist++ ) {
  596. hEC[iHist] = new TH1D( ECHistNames[iHist] + "_" + MpdPidUtils::cParticleShortName[MpdPidUtils::kHe4],"", 50, PMIN, PMAX);
  597. hEC[iHist]->Sumw2();
  598. vecEC.push_back( hEC[iHist] );
  599. }
  600. vecvecEC.push_back( vecEC );
  601. fEffContMap.insert( pair<MpdPidUtils::ePartType,std::vector<vecTH1Dptrs>>(MpdPidUtils::kHe4,vecvecEC) );
  602. }
  603. if ( fTrackingState == MpdPidUtils::kCFHM ) MaxEnLoss = 60.;
  604. else if ( fTrackingState == MpdPidUtils::kCF ) MaxEnLoss = 60000.;
  605. else MaxEnLoss = 60.e-06;
  606. fSumBetheBlochHists = new TH2D("fSumBetheBlochHists", "", 300, 0.05, 1.5, 300, MinEnLoss, MaxEnLoss);
  607. fChBetheBlochHists = new TH2D("fChBetheBlochHists", "", 600, -1.5, 1.5, 600, MinEnLoss, MaxEnLoss);
  608. fm2LightHist = new TH2D("fm2LightHist", "", 200, 0.0, 3.0, 200, -0.5, 1.5);
  609. fm2HeavyHist = new TH2D("fm2HeavyHist", "", 200, 0.0, 5.0, 200, 1.0, 12.0);
  610. for (Int_t i = 0; i < (MpdPidUtils::nQAHists - 8); i++) Xlow[i] = 0.05 * (1. + i);
  611. Xlow[32] = 1.7; Xlow[33] = 1.8; Xlow[34] = 1.9; Xlow[35] = 2.0; Xlow[36] = 2.1; Xlow[37] = 2.2; Xlow[38] = 2.4; Xlow[39] = 2.6;
  612. for (Int_t i = 0; i < (MpdPidUtils::nQAHists - 9); i++) Xhigh[i] = 0.05 + 0.05 * (1. + i);
  613. Xhigh[31] = 1.7; Xhigh[32] = 1.8; Xhigh[33] = 1.9; Xhigh[34] = 2.0; Xhigh[35] = 2.1; Xhigh[36] = 2.2; Xhigh[37] = 2.4; Xhigh[38] = 2.6; Xhigh[39] = 3.0;
  614. for (Int_t i=0; i<40; i++) X[i] = (Xlow[i] + Xhigh[i]) / 2;
  615. }
  616. MpdPidUtils::ePartType MpdPidQA::GetPartType(Int_t pdg) {
  617. map <Int_t,MpdPidUtils::ePartType>::iterator ret = fPartTypeMap.find(pdg);
  618. if ( ret != fPartTypeMap.end() ) return ret->second;
  619. else return MpdPidUtils::kUnknown;
  620. }
  621. void MpdPidQA::FillDedxHists(Double_t p, Double_t dedx, Int_t pdg) {
  622. Double_t sign = pdg > 0 ? 1.0 : -1.0;
  623. auto ret = fEnLossMap.find( GetPartType(pdg) );
  624. if ( ret != fEnLossMap.end() ) {
  625. fSumBetheBlochHists->Fill(p, dedx);
  626. fChBetheBlochHists->Fill(sign*p, dedx);
  627. ret->second.BetheBlochHist->Fill(p, dedx);
  628. for ( Int_t i = 0; i < MpdPidUtils::nQAHists; i++ ) {
  629. if ( (p > Xlow[i]) && (p <= Xhigh[i]) ) {
  630. if ( (i < ret->second.ibeg) || (i > ret->second.iend) ) continue;
  631. ret->second.dEdXPart[i - ret->second.ibeg]->Fill(dedx);
  632. }
  633. }
  634. }
  635. }
  636. void MpdPidQA::Fillm2Hists(Double_t p, Double_t m2, Int_t pdg) {
  637. Int_t pdgc = TMath::Abs(pdg);
  638. auto ret = fMSquaredMap.find( GetPartType(pdg) );
  639. if (ret != fMSquaredMap.end()) ret->second->Fill(p, m2);
  640. if ( ( pdgc == 211 ) || ( pdgc == 321 ) || ( pdgc == 2212 ) ) fm2LightHist->Fill(p, m2);
  641. if ( ( pdgc == PDG_DEUTERON ) || ( pdgc == PDG_TRITON ) || ( pdgc == PDG_HE3 ) || ( pdgc == PDG_HE4 ) ) fm2HeavyHist->Fill(p, m2);
  642. }
  643. void MpdPidQA::FillAmplHists(Double_t p, Int_t pdg) {
  644. if (p <= 5.0) {
  645. MpdPidUtils::ePartCharge ech;
  646. if ( ( TMath::Abs(pdg) == 11 ) || ( TMath::Abs(pdg) == 13 ) ) ech = pdg > 0 ? MpdPidUtils::kNeg : MpdPidUtils::kPos;
  647. else ech = pdg > 0 ? MpdPidUtils::kPos : MpdPidUtils::kNeg;
  648. auto ret = fAbundanceMap.find( GetPartType(pdg) );
  649. if ( ret != fAbundanceMap.end() ) ret->second[ech]->Fill(p);
  650. }
  651. }
  652. Bool_t MpdPidQA::FillEffContHists(Double_t p, Double_t dedx, Int_t charge, Int_t pdg, Double_t fProbCut) {
  653. FillEffDenominator(p, pdg);
  654. if ( !FillProbs(p, dedx, charge) ) return kFALSE;
  655. Int_t pidpdg = GetMaxProb();
  656. Double_t maxprob = GetProb( GetPartType(pidpdg) );
  657. if ( maxprob < fProbCut ) return kFALSE;
  658. MpdPidUtils::ePartCharge ech = charge > 0 ? MpdPidUtils::kPos : MpdPidUtils::kNeg;
  659. FillEffContHists(p, pidpdg, pdg, ech);
  660. return kTRUE;
  661. }
  662. Bool_t MpdPidQA::FillEffContHists(Double_t p, Double_t dedx, Double_t m2, Int_t charge, Int_t pdg, Double_t fProbCut) {
  663. FillEffDenominator(p, pdg);
  664. if ( !FillProbs(p, dedx, m2, charge) ) return kFALSE;
  665. Int_t pidpdg = GetMaxProb();
  666. Double_t maxprob = GetProb( GetPartType(pidpdg) );
  667. if ( maxprob < fProbCut ) return kFALSE;
  668. MpdPidUtils::ePartCharge ech = charge > 0 ? MpdPidUtils::kPos : MpdPidUtils::kNeg;
  669. FillEffContHists(p, pidpdg, pdg, ech);
  670. return kTRUE;
  671. }
  672. Bool_t MpdPidQA::FillEffContHists(Double_t p, Int_t pidpdg, Int_t pdg, MpdPidUtils::ePartCharge ech) {
  673. auto ret = fEffContMap.find( GetPartType(pidpdg) );
  674. if ( ret != fEffContMap.end() ) {
  675. ret->second[ech][1]->Fill(p);
  676. if ( pidpdg == pdg ) ret->second[ech][2]->Fill(p);
  677. else ret->second[ech][3]->Fill(p);
  678. }
  679. return kTRUE;
  680. }
  681. void MpdPidQA::FillEffDenominator(Double_t p, Int_t pdg) {
  682. auto ret = fEffContMap.find( GetPartType(pdg) );
  683. MpdPidUtils::ePartCharge ech;
  684. if ( ( TMath::Abs(pdg) == 11 ) || ( TMath::Abs(pdg) == 13 ) ) ech = pdg > 0 ? MpdPidUtils::kNeg : MpdPidUtils::kPos;
  685. else ech = pdg > 0 ? MpdPidUtils::kPos : MpdPidUtils::kNeg;
  686. if ( ret != fEffContMap.end() ) ret->second[ech][0]->Fill(p);
  687. }
  688. void MpdPidQA::GetDedxQA(TString dir) {
  689. Double_t XFUNCMIN, XFUNCMAX, mom, intF_MultX, intF;
  690. TString FName, FMXName;
  691. Int_t ibegloc, iendloc, nQAHistsloc, inam = 0;
  692. map <MpdPidUtils::ePartType,TGraphAsymmErrors*> graphs;
  693. map <MpdPidUtils::ePartType,vecTF1ptrs>::iterator it_parBB;
  694. vecTF1ptrs vecParBB, vecParBBMultX;
  695. for ( auto it = fEnLossMap.begin(); it != fEnLossMap.end(); ++it ) {
  696. ibegloc = it->second.ibeg; iendloc = it->second.iend; nQAHistsloc = iendloc - ibegloc + 1;
  697. it_parBB = fdEdxBBMap.find( it->first );
  698. vecParBB = it_parBB->second;
  699. for ( TF1* BBF : vecParBB ) {
  700. FName = BBF->GetName(); FName.Append("_MultX");
  701. TF1* parBBMultX = new TF1( FName,[&](Double_t *x, Double_t *p){ return p[0] * x[0] * BBF->Eval(x[0]); },BBF->GetXmin(),BBF->GetXmax(),1);
  702. parBBMultX->SetParameter(0, 1.0);
  703. vecParBBMultX.push_back(parBBMultX);
  704. }
  705. Double_t *Xgraph = new Double_t[nQAHistsloc]; Double_t *Xerrgraph = new Double_t[nQAHistsloc];
  706. Double_t *Ygraph = new Double_t[nQAHistsloc]; Double_t *sigma = new Double_t[nQAHistsloc];
  707. Double_t *sigma1 = new Double_t[nQAHistsloc]; Double_t *sigma2 = new Double_t[nQAHistsloc];
  708. for ( Int_t k = 0; k < nQAHistsloc; k++ ) { Xgraph[k] = X[k + ibegloc]; Xerrgraph[k] = 0.; }
  709. for ( Int_t i = 0; i < nQAHistsloc; i++ ) {
  710. XFUNCMIN = it->second.dEdXPart[i]->GetXaxis()->GetXmin(), XFUNCMAX = it->second.dEdXPart[i]->GetXaxis()->GetXmax();
  711. TF1* Gaus = new TF1("Gaus", "gaus", XFUNCMIN, XFUNCMAX);
  712. it->second.dEdXPart[i]->Fit(Gaus, "Q0R");
  713. TF1* AGaus = new TF1("AGaus", this, &MpdPid::AsymGaus, XFUNCMIN, XFUNCMAX, 4, "MpdPid", "AsymGaus");
  714. AGaus->SetParameters(Gaus->GetParameter(0), Gaus->GetParameter(1), Gaus->GetParameter(2), 0.01);
  715. it->second.dEdXPart[i]->Fit("AGaus","Q0RW");
  716. TF1 *Novosib = new TF1("Novosib", this, &MpdPidQA::Novosibirsk, XFUNCMIN, XFUNCMAX, 4, "MpdPidQA", "Novosibirsk");
  717. Novosib->SetParameters(Gaus->GetParameter(0), 0.01, Gaus->GetParameter(2), Gaus->GetParameter(1));
  718. it->second.dEdXPart[i]->Fit("Novosib","Q0RW");
  719. intF_MultX = 0.0; intF = 0.0;
  720. for ( TF1* BBF : vecParBB ) {
  721. /// Xlow and Xhigh are within one Bethe-Bloch function
  722. if ( ( BBF->GetXmin() <= Xlow[i+ibegloc] ) && ( BBF->GetXmax() > Xlow[i+ibegloc] ) && ( BBF->GetXmin() <= Xhigh[i+ibegloc] ) && ( BBF->GetXmax() > Xhigh[i+ibegloc] ) ) {
  723. for ( TF1* BBF_MultX : vecParBBMultX ) {
  724. FName = BBF->GetName(); FName.Append("_MultX");
  725. FMXName = BBF_MultX->GetName();
  726. if ( FMXName == FName ) {
  727. intF_MultX = BBF_MultX->Integral( Xlow[i+ibegloc], Xhigh[i+ibegloc] );
  728. intF = BBF->Integral( Xlow[i+ibegloc], Xhigh[i+ibegloc] );
  729. mom = intF_MultX / intF;
  730. }
  731. }
  732. /// Xlow and Xhigh are within two functions
  733. } else if ( ( BBF->GetXmin() <= Xlow[i+ibegloc] ) && ( BBF->GetXmax() > Xlow[i+ibegloc] ) ) {
  734. intF += BBF->Integral( Xlow[i+ibegloc], BBF->GetXmax() );
  735. for ( TF1* BBF_MultX : vecParBBMultX ) {
  736. FName = BBF->GetName(); FName.Append("_MultX");
  737. FMXName = BBF_MultX->GetName();
  738. if ( FMXName == FName ) {
  739. intF_MultX += BBF_MultX->Integral( Xlow[i+ibegloc], BBF_MultX->GetXmax() );
  740. }
  741. }
  742. } else if ( ( BBF->GetXmin() > Xlow[i+ibegloc] ) && ( BBF->GetXmax() < Xhigh[i+ibegloc] ) ) {
  743. if ( intF != 0.0 ) {
  744. intF += BBF->Integral( BBF->GetXmin(), BBF->GetXmax() );
  745. for ( TF1* BBF_MultX : vecParBBMultX ) {
  746. FName = BBF->GetName(); FName.Append("_MultX");
  747. FMXName = BBF_MultX->GetName();
  748. if ( FMXName == FName ) {
  749. intF_MultX += BBF_MultX->Integral( BBF->GetXmin(), BBF->GetXmax() );
  750. }
  751. }
  752. }
  753. } else if ( ( BBF->GetXmin() <= Xhigh[i+ibegloc] ) && ( BBF->GetXmax() > Xhigh[i+ibegloc] ) ) {
  754. if ( intF != 0.0 ) {
  755. intF += BBF->Integral( BBF->GetXmin(), Xhigh[i+ibegloc] );
  756. for ( TF1* BBF_MultX : vecParBBMultX ) {
  757. FName = BBF->GetName(); FName.Append("_MultX");
  758. FMXName = BBF_MultX->GetName();
  759. if ( FMXName == FName ) {
  760. intF_MultX += BBF_MultX->Integral( BBF->GetXmin(), Xhigh[i+ibegloc] );
  761. mom = intF_MultX / intF;
  762. }
  763. }
  764. }
  765. } else mom = -999.0;
  766. }
  767. Ygraph[i] = ( fTrackingState == MpdPidUtils::kCFHM ) ? AGaus->GetParameter(1) : Novosib->GetParameter(3);
  768. Ygraph[i] /= GetDedxParam(mom,it->first);
  769. sigma1[i] = AGaus->GetParameter(2); sigma2[i] = (1. + AGaus->GetParameter(3)) * AGaus->GetParameter(2);
  770. sigma1[i] /= GetDedxParam(mom,it->first); sigma2[i] /= GetDedxParam(mom,it->first);
  771. delete Gaus; delete AGaus; delete Novosib;
  772. }
  773. TGraphAsymmErrors *gr = new TGraphAsymmErrors(nQAHistsloc,Xgraph,Ygraph,Xerrgraph,Xerrgraph,sigma1,sigma2);
  774. FName = "Graph_"; FName.Append( MpdPidUtils::cParticleShortName[it->first] );
  775. gr->SetName(FName);
  776. graphs.insert( pair<MpdPidUtils::ePartType,TGraphAsymmErrors*>(it->first,gr) );
  777. delete Xgraph; delete Ygraph; delete Xerrgraph; delete sigma;
  778. }
  779. if ( !dir.EndsWith("/") ) dir.Append("/");
  780. TFile *fFile = new TFile(dir + "dEdxHists.root", "RECREATE");
  781. for ( auto it = fEnLossMap.begin(); it != fEnLossMap.end(); ++it ) {
  782. for ( Int_t i = 0; i < ( it->second.iend - it->second.ibeg + 1 ); i++ ) {
  783. it->second.dEdXPart[i]->Write();
  784. }
  785. it->second.BetheBlochHist->Write();
  786. }
  787. for ( auto it = graphs.begin(); it != graphs.end(); ++it ) {
  788. it->second->Write();
  789. }
  790. fSumBetheBlochHists->Write();
  791. fChBetheBlochHists->Write();
  792. fFile->Close();
  793. }
  794. void MpdPidQA::Getm2QA(TString dir) {
  795. if ( !dir.EndsWith("/") ) dir.Append("/");
  796. TFile *fFile = new TFile(dir + "m2Hists.root", "RECREATE");
  797. for ( auto it = fMSquaredMap.begin(); it != fMSquaredMap.end(); ++it ) {
  798. it->second->Write();
  799. }
  800. fm2LightHist->Write();
  801. fm2HeavyHist->Write();
  802. fFile->Close();
  803. }
  804. void MpdPidQA::GetAmplQA(TString dir) {
  805. if ( !dir.EndsWith("/") ) dir.Append("/");
  806. TFile *fFile = new TFile(dir + "ParticleYields.root", "RECREATE");
  807. for ( auto it = fAbundanceMap.begin(); it != fAbundanceMap.end(); ++it ) {
  808. for ( TH1D* hYield : it->second ) {
  809. hYield->Write();
  810. }
  811. }
  812. fFile->Close();
  813. }
  814. void MpdPidQA::GetEffContQA(TString dir) {
  815. if ( !dir.EndsWith("/") ) dir.Append("/");
  816. TFile *fFile = new TFile(dir + "EffCont.root", "RECREATE");
  817. for ( auto it = fEffContMap.begin(); it != fEffContMap.end(); ++it ) {
  818. for ( vecTH1Dptrs vecEC : it->second ) {
  819. vecEC[2]->Divide( vecEC[0] );
  820. vecEC[3]->Divide( vecEC[1] );
  821. vecEC[2]->Write();
  822. vecEC[3]->Write();
  823. }
  824. }
  825. fFile->Close();
  826. }
  827. Double_t MpdPidQA::Novosibirsk(Double_t *x, Double_t *par) {
  828. Double_t tail = par[1];
  829. Double_t width = par[2];
  830. Double_t peak = par[3];
  831. if (TMath::Abs(tail) < 1.e-7)
  832. return par[0]*TMath::Exp( -0.5 * TMath::Power( ( (x[0] - peak) / width ), 2 ));
  833. Double_t arg = 1.0 - ( x[0] - peak ) * tail / width;
  834. if (arg < 1.e-6) return 0.0; ///Argument of logaritem negative. Real continuation -> function equals zero
  835. Double_t log = TMath::Log(arg);
  836. static const Double_t xi = 2.3548200450309494; /// 2 Sqrt( Ln(4) )
  837. Double_t width_zero = ( 2.0 / xi ) * TMath::ASinH( tail * xi * 0.5 );
  838. Double_t width_zero2 = width_zero * width_zero;
  839. Double_t exponent = ( -0.5 / (width_zero2) * log * log ) - ( width_zero2 * 0.5 );
  840. return par[0]*TMath::Exp(exponent);
  841. }
  842. ClassImp(MpdPidQA);