MpdFemtoModelCorrFctnMomResolution.cxx 20 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620
  1. #include "MpdFemtoModelCorrFctnMomResolution.h"
  2. //_______________________
  3. MpdFemtoModelCorrFctnMomResolution::MpdFemtoModelCorrFctnMomResolution() :
  4. fManager(nullptr),
  5. fNumeratorTrue(nullptr),
  6. fNumeratorFake(nullptr),
  7. fDenominator(nullptr),
  8. fNumeratorTrueIdeal(nullptr),
  9. fNumeratorFakeIdeal(nullptr),
  10. fDenominatorIdeal(nullptr),
  11. fQgenQrec(nullptr),
  12. fKaonPDG(kFALSE),
  13. fFillkT(kFALSE) {
  14. // Default constructor
  15. fNumeratorTrue = new TH3D("ModelNumTrue", "ModelNumTrue", 50, 0.0, 0.5, 50, 0.0, 0.5, 50, 0.0, 0.5);
  16. fNumeratorFake = new TH3D("ModelNumFake", "ModelNumFake", 50, 0.0, 0.5, 50, 0.0, 0.5, 50, 0.0, 0.5);
  17. fDenominator = new TH3D("ModelDen", "ModelDen", 50, 0.0, 0.5, 50, 0.0, 0.5, 50, 0.0, 0.5);
  18. fNumeratorTrueIdeal = new TH3D("ModelNumTrueIdeal", "ModelNumTrueIdeal", 50, 0.0, 0.5, 50, 0.0, 0.5, 50, 0.0, 0.5);
  19. fNumeratorFakeIdeal = new TH3D("ModelNumFakeIdeal", "ModelNumFakeIdeal", 50, 0.0, 0.5, 50, 0.0, 0.5, 50, 0.0, 0.5);
  20. fDenominatorIdeal = new TH3D("ModelDenIdeal", "ModelDenIdeal", 50, 0.0, 0.5, 50, 0.0, 0.5, 50, 0.0, 0.5);
  21. fQgenQrec = new TH2D("QgenQrec", "QgenQrec", 50, 0.0, 0.5, 50, 0.0, 0.5);
  22. for (int i = 0; i < fNbbPairs; i++) {
  23. auto *name = Form("fkTdists[%i]", i);
  24. fkTdists[i] = new TH1D(name, name, 100, 0.0, 5.0);
  25. fkTdists[i]->Sumw2();
  26. }
  27. fNumeratorTrue->Sumw2();
  28. fNumeratorFake->Sumw2();
  29. fDenominator->Sumw2();
  30. fNumeratorTrueIdeal->Sumw2();
  31. fNumeratorFakeIdeal->Sumw2();
  32. fDenominatorIdeal->Sumw2();
  33. fQgenQrec->Sumw2();
  34. }
  35. //_______________________
  36. MpdFemtoModelCorrFctnMomResolution::MpdFemtoModelCorrFctnMomResolution(const char *title,
  37. Int_t aNbins,
  38. Double_t aQinvLo,
  39. Double_t aQinvHi) :
  40. fManager(nullptr),
  41. fNumeratorTrue(nullptr),
  42. fNumeratorFake(nullptr),
  43. fDenominator(nullptr),
  44. fNumeratorTrueIdeal(nullptr),
  45. fNumeratorFakeIdeal(nullptr),
  46. fDenominatorIdeal(nullptr),
  47. fQgenQrec(nullptr),
  48. fKaonPDG(kFALSE),
  49. fFillkT(kFALSE) {
  50. // Normal constructor
  51. char *buf;
  52. buf = Form("NumTrue%s", title);
  53. fNumeratorTrue = new TH3D(buf, buf, aNbins, aQinvLo, aQinvHi, aNbins, aQinvLo, aQinvHi, aNbins, aQinvLo, aQinvHi);
  54. buf = Form("NumFake%s", title);
  55. fNumeratorFake = new TH3D(buf, buf, aNbins, aQinvLo, aQinvHi, aNbins, aQinvLo, aQinvHi, aNbins, aQinvLo, aQinvHi);
  56. buf = Form("Den%s", title);
  57. fDenominator = new TH3D(buf, buf, aNbins, aQinvLo, aQinvHi, aNbins, aQinvLo, aQinvHi, aNbins, aQinvLo, aQinvHi);
  58. buf = Form("NumTrueIdeal%s", title);
  59. fNumeratorTrueIdeal = new TH3D(buf, buf, aNbins, aQinvLo, aQinvHi, aNbins, aQinvLo, aQinvHi, aNbins, aQinvLo, aQinvHi);
  60. buf = Form("NumFakeIdeal%s", title);
  61. fNumeratorFakeIdeal = new TH3D(buf, buf, aNbins, aQinvLo, aQinvHi, aNbins, aQinvLo, aQinvHi, aNbins, aQinvLo, aQinvHi);
  62. buf = Form("DenIdeal%s", title);
  63. fDenominatorIdeal = new TH3D(buf, buf, aNbins, aQinvLo, aQinvHi, aNbins, aQinvLo, aQinvHi, aNbins, aQinvLo, aQinvHi);
  64. buf = Form("QgenQrec%s", title);
  65. fQgenQrec = new TH2D(buf, buf, aNbins, aQinvLo, aQinvHi, aNbins, aQinvLo, aQinvHi);
  66. //test
  67. //fQgenQrec = new TH2D(buf,buf,aNbins,aQinvLo,aQinvHi,aNbins,-0.05,0.05);
  68. for (int i = 0; i < fNbbPairs; i++) {
  69. const char *name = Form("fkTdists[%i]_%s", i, title);
  70. fkTdists[i] = new TH1D(name, name, 100, 0.0, 5.0);
  71. fkTdists[i]->Sumw2();
  72. }
  73. fNumeratorTrue->Sumw2();
  74. fNumeratorFake->Sumw2();
  75. fDenominator->Sumw2();
  76. fNumeratorTrueIdeal->Sumw2();
  77. fNumeratorFakeIdeal->Sumw2();
  78. fDenominatorIdeal->Sumw2();
  79. fQgenQrec->Sumw2();
  80. }
  81. //_______________________
  82. MpdFemtoModelCorrFctnMomResolution::MpdFemtoModelCorrFctnMomResolution(const MpdFemtoModelCorrFctnMomResolution& aCorrFctn) :
  83. fManager(aCorrFctn.fManager),
  84. fNumeratorTrue(nullptr),
  85. fNumeratorFake(nullptr),
  86. fDenominator(nullptr),
  87. fNumeratorTrueIdeal(nullptr),
  88. fNumeratorFakeIdeal(nullptr),
  89. fDenominatorIdeal(nullptr),
  90. fQgenQrec(nullptr),
  91. fKaonPDG(aCorrFctn.fKaonPDG),
  92. fFillkT(aCorrFctn.fFillkT) {
  93. // Copy constructor
  94. fNumeratorTrue = new TH3D(*aCorrFctn.fNumeratorTrue);
  95. fNumeratorFake = new TH3D(*aCorrFctn.fNumeratorFake);
  96. fDenominator = new TH3D(*aCorrFctn.fDenominator);
  97. fNumeratorTrueIdeal = new TH3D(*aCorrFctn.fNumeratorTrueIdeal);
  98. fNumeratorFakeIdeal = new TH3D(*aCorrFctn.fNumeratorFakeIdeal);
  99. fDenominatorIdeal = new TH3D(*aCorrFctn.fDenominatorIdeal);
  100. fQgenQrec = new TH2D(*aCorrFctn.fQgenQrec);
  101. for (int i = 0; i < fNbbPairs; i++) {
  102. fkTdists[i] = new TH1D(*aCorrFctn.fkTdists[i]);
  103. }
  104. }
  105. //_______________________
  106. MpdFemtoModelCorrFctnMomResolution::~MpdFemtoModelCorrFctnMomResolution() {
  107. // Destructor
  108. delete fNumeratorTrue;
  109. delete fNumeratorFake;
  110. delete fDenominator;
  111. delete fNumeratorTrueIdeal;
  112. delete fNumeratorFakeIdeal;
  113. delete fDenominatorIdeal;
  114. delete fQgenQrec;
  115. for (int i = 0; i < fNbbPairs; i++) {
  116. delete fkTdists[i];
  117. }
  118. }
  119. //_______________________
  120. MpdFemtoModelCorrFctnMomResolution& MpdFemtoModelCorrFctnMomResolution::operator=(const MpdFemtoModelCorrFctnMomResolution& aCorrFctn) {
  121. // Assignment operator
  122. if (this == &aCorrFctn) {
  123. return *this;
  124. }
  125. MpdFemtoModelCorrFctnMomResolution::operator=(aCorrFctn);
  126. *fNumeratorTrue = *aCorrFctn.fNumeratorTrue;
  127. *fNumeratorFake = *aCorrFctn.fNumeratorFake;
  128. *fDenominator = *aCorrFctn.fDenominator;
  129. *fQgenQrec = *aCorrFctn.fQgenQrec;
  130. for (int i = 0; i < fNbbPairs; i++) {
  131. *fkTdists[i] = *aCorrFctn.fkTdists[i];
  132. }
  133. *fNumeratorTrueIdeal = *aCorrFctn.fNumeratorTrueIdeal;
  134. *fNumeratorFakeIdeal = *aCorrFctn.fNumeratorFakeIdeal;
  135. *fDenominatorIdeal = *aCorrFctn.fDenominatorIdeal;
  136. fManager = aCorrFctn.fManager;
  137. fKaonPDG = aCorrFctn.fKaonPDG;
  138. return *this;
  139. }
  140. //_______________________
  141. void MpdFemtoModelCorrFctnMomResolution::connectToManager(MpdFemtoModelManager *aManager) {
  142. fManager = aManager;
  143. }
  144. //_______________________
  145. MpdFemtoString MpdFemtoModelCorrFctnMomResolution::report() {
  146. // Construct the report
  147. TString report = "";
  148. return MpdFemtoString((const char *) report);
  149. }
  150. //_______________________
  151. void MpdFemtoModelCorrFctnMomResolution::addRealPair(MpdFemtoPair* aPair) {
  152. if (mPairCut && !mPairCut->pass(aPair)) {
  153. return;
  154. }
  155. // if (!fKaonPDG) {
  156. // cout<<" AliFemtoModelCorrFcn add real pair "<<endl;
  157. Double_t weight = fManager->weight(aPair);
  158. //cout<<" wight "<< weight<<endl;
  159. //cout<<"Qinv"<<aPair->QInv()<<endl;
  160. fNumeratorTrue->Fill(aPair->qOutCMS(), aPair->qSideCMS(), aPair->qLongCMS(), weight); //qOut, qSide, qLong
  161. //Double_t tQinvTrue = GetQinvTrue(aPair);
  162. fNumeratorTrueIdeal->Fill(GetQoutTrue(aPair), GetQsideTrue(aPair), GetQlongTrue(aPair), weight); //qOut_ideal, qSide_ideal, qLong_ideal
  163. //cout<<"Qinv true"<<tQinvTrue<<endl;
  164. // }//Special MC analysis for K selected by PDG code -->
  165. // else {
  166. // Double_t weight = fManager->GetWeight(aPair);
  167. // fNumeratorTrue->Fill(aPair->QOutCMS(), aPair->QSideCMS(), aPair->QLongCMS(), weight); //qOut, qSide, qLong
  168. // //Double_t tQinvTrue = GetQinvTrue(aPair);
  169. // fNumeratorTrueIdeal->Fill(GetQoutTrue(aPair), GetQsideTrue(aPair), GetQlongTrue(aPair), weight); //qOut_ideal, qSide_ideal, qLong_ideal
  170. // }
  171. }
  172. //_______________________
  173. void MpdFemtoModelCorrFctnMomResolution::addMixedPair(MpdFemtoPair* aPair) {
  174. if (mPairCut && !mPairCut->pass(aPair)) {
  175. return;
  176. }
  177. Double_t weight = fManager->weight(aPair);
  178. Double_t qinv = aPair->qInv();
  179. Double_t qinv_ideal = GetQinvTrue(aPair);
  180. // if (!fKaonPDG) {
  181. fNumeratorFake->Fill(aPair->qOutCMS(), aPair->qSideCMS(), aPair->qLongCMS(), weight); //qOut, qSide, qLong
  182. fDenominator->Fill(aPair->qOutCMS(), aPair->qSideCMS(), aPair->qLongCMS(), 1.0); //qOut, qSide, qLong
  183. fNumeratorFakeIdeal->Fill(GetQoutTrue(aPair), GetQsideTrue(aPair), GetQlongTrue(aPair), weight); //qOut_ideal, qSide_ideal, qLong_ideal
  184. fDenominatorIdeal->Fill(GetQoutTrue(aPair), GetQsideTrue(aPair), GetQlongTrue(aPair), 1.0); //qOut_ideal, qSide_ideal, qLong_ideal
  185. if (fFillkT) {
  186. int pairNumber = GetPairNumber(aPair);
  187. if (pairNumber >= 0) {
  188. if (fkTdists[pairNumber]) {
  189. fkTdists[pairNumber]->Fill(GetParentsKt(aPair));
  190. }
  191. }
  192. }
  193. fQgenQrec->Fill(qinv_ideal, qinv);
  194. // }//Special MC analysis for K selected by PDG code -->
  195. // else {
  196. // // AliFemtoTrack *inf1 = (AliFemtoTrack *) aPair->Track1()->Track();
  197. // // AliFemtoTrack *inf2 = (AliFemtoTrack *) aPair->Track2()->Track();
  198. // // Double_t pdg1 = ((AliFemtoModelHiddenInfo*)inf1->GetHiddenInfo())->GetPDGPid();
  199. // // Double_t pdg2 = ((AliFemtoModelHiddenInfo*)inf2->GetHiddenInfo())->GetPDGPid();
  200. // // if((aPair->KT())<0.5)cout<<" Corr Func pdg1 "<<pdg1<<" pdg2 "<<pdg2<<" qinv "<<aPair->QInv()<< " w "<<weight<<endl;
  201. // fNumeratorFake->Fill(aPair->QOutCMS(), aPair->QSideCMS(), aPair->QLongCMS(), weight); //qOut, qSide, qLong
  202. // fDenominator->Fill(aPair->QOutCMS(), aPair->QSideCMS(), aPair->QLongCMS(), 1.0); //qOut, qSide, qLong
  203. // //if(qinv_ideal>0) {
  204. // fNumeratorFakeIdeal->Fill(GetQoutTrue(aPair), GetQsideTrue(aPair), GetQlongTrue(aPair), weight); //qOut_ideal, qSide_ideal, qLong_ideal
  205. // fDenominatorIdeal->Fill(GetQoutTrue(aPair), GetQsideTrue(aPair), GetQlongTrue(aPair), 1.0); //qOut_ideal, qSide_ideal, qLong_ideal
  206. // fQgenQrec->Fill(qinv_ideal, qinv);
  207. // //}
  208. // //test
  209. // //if(tQinvTrue>0)fQgenQrec->Fill(tQinvTrue,tQinvTrue-aPair->QInv());
  210. // }
  211. }
  212. //void MpdFemtoModelCorrFctnMomResolution::SetSpecificPairCut(AliFemtoPairCut* aCut)
  213. //{
  214. // fPairCut = aCut;
  215. //}
  216. //_______________________
  217. Double_t MpdFemtoModelCorrFctnMomResolution::GetQinvTrue(MpdFemtoPair* aPair) {
  218. //if(!fKaonPDG) {
  219. MpdFemtoParticle* first = aPair->track1();
  220. MpdFemtoParticle* second = aPair->track2();
  221. if (!first || !second)
  222. return -1;
  223. MpdFemtoModelHiddenInfo* inf1 = (MpdFemtoModelHiddenInfo*) first->getHiddenInfo();
  224. MpdFemtoModelHiddenInfo* inf2 = (MpdFemtoModelHiddenInfo*) second->getHiddenInfo();
  225. if (!inf1 || !inf2)
  226. return -1;
  227. TVector3 mom1 = inf1->trueMomentum();
  228. Double_t mass1 = inf1->mass();
  229. Double_t ene1 = TMath::Sqrt(mom1.Mag2() + mass1 * mass1);
  230. TLorentzVector fm1(mom1.Px(), mom1.Py(), mom1.Pz(), ene1);
  231. TVector3 mom2 = inf2->trueMomentum();
  232. Double_t mass2 = inf2->mass();
  233. Double_t ene2 = TMath::Sqrt(mom2.Mag2() + mass2 * mass2);
  234. TLorentzVector fm2(mom2.Px(), mom2.Py(), mom2.Pz(), ene2);
  235. //std::cout<<" CFModel mass1 mass2 "<<am1<<" "<<am2<<std::endl;
  236. TLorentzVector tQinvTrueVec = fm1 - fm2;
  237. Double_t tQinvTrue = -1. * tQinvTrueVec.M();
  238. return tQinvTrue;
  239. //}
  240. //Special MC analysis for K selected by PDG code -->
  241. // else {
  242. // const AliFemtoTrack *inf1 = aPair->Track1()->Track(),
  243. // *inf2 = aPair->Track2()->Track();
  244. //
  245. // AliFemtoLorentzVector fm1;
  246. // AliFemtoThreeVector* temp = ((AliFemtoModelHiddenInfo*)inf1->GetHiddenInfo())->GetTrueMomentum();
  247. // fm1.SetVect(*temp);
  248. // Double_t am1 = ((AliFemtoModelHiddenInfo*)inf1->GetHiddenInfo())->GetMass();
  249. // Double_t am2 = ((AliFemtoModelHiddenInfo*)inf2->GetHiddenInfo())->GetMass();
  250. //
  251. // am1=0.493677;
  252. // am2=0.493677;
  253. //
  254. // //Double_t pdg1 = ((AliFemtoModelHiddenInfo*)inf1->GetHiddenInfo())->GetPDGPid();
  255. // //Double_t pdg2 = ((AliFemtoModelHiddenInfo*)inf2->GetHiddenInfo())->GetPDGPid();
  256. //
  257. // double ener = TMath::Sqrt(temp->Mag2()+am1*am1);
  258. // fm1.SetE(ener);
  259. //
  260. // AliFemtoLorentzVector fm2;
  261. // AliFemtoThreeVector* temp2 = ((AliFemtoModelHiddenInfo*)inf2->GetHiddenInfo())->GetTrueMomentum();
  262. // fm2.SetVect(*temp2);
  263. // ener = TMath::Sqrt(temp2->Mag2()+am2*am2);
  264. // fm2.SetE(ener);
  265. //
  266. //
  267. // AliFemtoLorentzVector tQinvTrueVec = (fm1-fm2);
  268. // Double_t tQinvTrue = -1.* tQinvTrueVec.m();
  269. //
  270. // // if(tQinvTrue<0 && am1!=0 && am2!=0)std::cout<<" CFModel Qinv mass1 mass2 "<<aPair->QInv()<<" Qinv_true "<<tQinvTrue<<" "<<am1<<" "<<am2<<" pdg1 "<<pdg1<<" pdg2 "<<pdg2<<std::endl;
  271. // // if(pdg1!=211 || pdg2!=211)std::cout<<" CFModel Qinv mass1 mass2 "<<aPair->QInv()<<" Qinv_true "<<tQinvTrue<<" "<<am1<<" "<<am2<<" pdg1 "<<pdg1<<" pdg2 "<<pdg2<<std::endl;
  272. //
  273. // if(am1==0 || am2==0)tQinvTrue=-10000;
  274. //
  275. // return tQinvTrue;
  276. // }
  277. }
  278. //_______________________
  279. Double_t MpdFemtoModelCorrFctnMomResolution::GetQoutTrue(MpdFemtoPair* aPair) {
  280. MpdFemtoParticle* first = aPair->track1();
  281. MpdFemtoParticle* second = aPair->track2();
  282. if (!first || !second)
  283. return -1;
  284. MpdFemtoModelHiddenInfo* inf1 = (MpdFemtoModelHiddenInfo*) first->getHiddenInfo();
  285. MpdFemtoModelHiddenInfo* inf2 = (MpdFemtoModelHiddenInfo*) second->getHiddenInfo();
  286. if (!inf1 || !inf2)
  287. return -1;
  288. TVector3 p1 = inf1->trueMomentum();
  289. TVector3 p2 = inf2->trueMomentum();
  290. Double_t dx = p1.Px() - p2.Px();
  291. Double_t dy = p1.Py() - p2.Py();
  292. Double_t px = p1.Px() + p2.Px();
  293. Double_t py = p1.Py() + p2.Py();
  294. Double_t qTkT = dx * px + dy * py;
  295. Double_t pT = TMath::Sqrt(px * px + py * py);
  296. // Double_t qOut;
  297. // if (!pT)
  298. // qOut = 0;
  299. // else
  300. // qOut = qTkT / pT;
  301. // if (!qOut)
  302. // qOut = -10000;
  303. //return qOut;
  304. return qTkT / pT;
  305. }
  306. Double_t MpdFemtoModelCorrFctnMomResolution::GetQsideTrue(MpdFemtoPair* aPair) {
  307. MpdFemtoParticle* first = aPair->track1();
  308. MpdFemtoParticle* second = aPair->track2();
  309. if (!first || !second)
  310. return -1;
  311. MpdFemtoModelHiddenInfo* inf1 = (MpdFemtoModelHiddenInfo*) first->getHiddenInfo();
  312. MpdFemtoModelHiddenInfo* inf2 = (MpdFemtoModelHiddenInfo*) second->getHiddenInfo();
  313. if (!inf1 || !inf2)
  314. return -1;
  315. TVector3 p1 = inf1->trueMomentum();
  316. TVector3 p2 = inf2->trueMomentum();
  317. Double_t x1 = p1.Px();
  318. Double_t x2 = p2.Px();
  319. Double_t y1 = p1.Py();
  320. Double_t y2 = p2.Py();
  321. Double_t px = p1.Px() + p2.Px();
  322. Double_t py = p1.Py() + p2.Py();
  323. Double_t qTkT = x2 * y1 - x1 * y2;
  324. Double_t pT = sqrt(px * px + py * py);
  325. // Double_t qSide;
  326. // if (!pT)
  327. // qSide = 0;
  328. // else
  329. // qSide = (2.0 * qTkT) / pT;
  330. //
  331. // if (!qSide)
  332. // qSide = -10000;
  333. //
  334. // return qSide;
  335. return (2. * qTkT) / pT;
  336. }
  337. Double_t MpdFemtoModelCorrFctnMomResolution::GetQlongTrue(MpdFemtoPair* aPair) {
  338. MpdFemtoParticle* first = aPair->track1();
  339. MpdFemtoParticle* second = aPair->track2();
  340. if (!first || !second)
  341. return -1;
  342. MpdFemtoModelHiddenInfo* inf1 = (MpdFemtoModelHiddenInfo*) first->getHiddenInfo();
  343. MpdFemtoModelHiddenInfo* inf2 = (MpdFemtoModelHiddenInfo*) second->getHiddenInfo();
  344. if (!inf1 || !inf2)
  345. return -1;
  346. TVector3 mom1 = inf1->trueMomentum();
  347. Double_t mass1 = inf1->mass();
  348. Double_t ene1 = TMath::Sqrt(mom1.Mag2() + mass1 * mass1);
  349. TLorentzVector fm1(mom1.Px(), mom1.Py(), mom1.Pz(), ene1);
  350. TVector3 mom2 = inf2->trueMomentum();
  351. Double_t mass2 = inf2->mass();
  352. Double_t ene2 = TMath::Sqrt(mom2.Mag2() + mass2 * mass2);
  353. TLorentzVector fm2(mom2.Px(), mom2.Py(), mom2.Pz(), ene2);
  354. //std::cout<<" CFModel mass1 mass2 "<<am1<<" "<<am2<<std::endl;
  355. Double_t dz = fm1.Z() - fm2.Z();
  356. Double_t zz = fm1.Z() + fm2.Z();
  357. Double_t dt = fm1.T() - fm2.T();
  358. Double_t tt = fm1.T() + fm2.T();
  359. Double_t beta = zz / tt;
  360. Double_t gamma = 1.0 / TMath::Sqrt((1. - beta)*(1. + beta));
  361. Double_t qLong = gamma * (dz - beta * dt);
  362. // if (!qLong)
  363. // qLong = -10000;
  364. return qLong;
  365. }
  366. //_______________________
  367. void MpdFemtoModelCorrFctnMomResolution::eventBegin(const MpdFemtoEvent* /* aEvent */) {
  368. /* Do nothing */
  369. }
  370. //_______________________
  371. void MpdFemtoModelCorrFctnMomResolution::eventEnd(const MpdFemtoEvent* /* aEvent */) {
  372. /* Do nothing */
  373. }
  374. //_______________________
  375. void MpdFemtoModelCorrFctnMomResolution::finish() {
  376. /* Do nothing */
  377. }
  378. //_______________________
  379. void MpdFemtoModelCorrFctnMomResolution::writeOutHistos() {
  380. // Write out data histos
  381. fQgenQrec->Write();
  382. fNumeratorTrue->Write();
  383. fNumeratorFake->Write();
  384. fDenominator->Write();
  385. if (fFillkT) {
  386. for (int i = 0; i < fNbbPairs; i++) {
  387. fkTdists[i]->Write();
  388. }
  389. }
  390. fNumeratorTrueIdeal->Write();
  391. fNumeratorFakeIdeal->Write();
  392. fDenominatorIdeal->Write();
  393. }
  394. //_________________________
  395. TList* MpdFemtoModelCorrFctnMomResolution::getOutputList() {
  396. // Prepare the list of objects to be written to the output
  397. TList *tOutputList = new TList();
  398. tOutputList->Add(fNumeratorTrue);
  399. tOutputList->Add(fNumeratorFake);
  400. tOutputList->Add(fDenominator);
  401. tOutputList->Add(fNumeratorTrueIdeal);
  402. tOutputList->Add(fNumeratorFakeIdeal);
  403. tOutputList->Add(fDenominatorIdeal);
  404. tOutputList->Add(fQgenQrec);
  405. if (fFillkT) {
  406. for (int i = 0; i < fNbbPairs; i++) {
  407. tOutputList->Add(fkTdists[i]);
  408. }
  409. }
  410. return tOutputList;
  411. }
  412. void MpdFemtoModelCorrFctnMomResolution::SetKaonPDG(Bool_t aSetKaonAna) {
  413. fKaonPDG = aSetKaonAna;
  414. }
  415. Double_t MpdFemtoModelCorrFctnMomResolution::GetParentsKt(MpdFemtoPair *pair) {
  416. // MpdFemtoParticle* first = new MpdFemtoParticle(*(pair->track1()));
  417. // MpdFemtoParticle* second = new MpdFemtoParticle(*(pair->track2()));
  418. //
  419. // if (!first) {
  420. // if (second)
  421. // delete second;
  422. // return -1;
  423. // }
  424. //
  425. // if (!second) {
  426. // if (first)
  427. // delete first;
  428. // return -1;
  429. // }
  430. //
  431. // MpdFemtoModelHiddenInfo* info1 = first->getHiddenInfo();
  432. // MpdFemtoModelHiddenInfo* info2 = second->getHiddenInfo();
  433. //
  434. // if (!info1 || !info2) {
  435. // if (first) delete first;
  436. // if (second) delete second;
  437. // return -1;
  438. // }
  439. // AliFemtoThreeVector* p1 = info1->GetMotherMomentum();
  440. // AliFemtoThreeVector* p2 = info2->GetMotherMomentum();
  441. //
  442. // if (!p1 || !p2) {
  443. // if (first) delete first;
  444. // if (second) delete second;
  445. // return -1;
  446. // }
  447. // double px = p1->x() + p2->x();
  448. // double py = p1->y() + p2->y();
  449. // double pT = sqrt(px * px + py * py);
  450. //
  451. // delete first;
  452. // delete second;
  453. //
  454. // return pT / 2.;
  455. return 0.;
  456. }
  457. //
  458. Int_t MpdFemtoModelCorrFctnMomResolution::GetPairNumber(MpdFemtoPair *pair) {
  459. // const AliFemtoModelHiddenInfo
  460. // *info1 = (AliFemtoModelHiddenInfo*) pair->Track1()->GetHiddenInfo(),
  461. // *info2 = (AliFemtoModelHiddenInfo*) pair->Track2()->GetHiddenInfo();
  462. //
  463. // if (!info1 || !info2) {
  464. // return -1;
  465. // }
  466. //
  467. // int pdg1 = TMath::Abs(info1->GetMotherPdgCode());
  468. // int pdg2 = TMath::Abs(info2->GetMotherPdgCode());
  469. //
  470. // if (pdg2 < pdg1) {
  471. // std::swap(pdg1, pdg2);
  472. // }
  473. //
  474. // if (pdg1 == 2212 && pdg2 == 2212) return 0; // pp
  475. // if (pdg1 == 2212 && pdg2 == 3122) return 1; // pΛ
  476. // if (pdg1 == 3122 && pdg2 == 3122) return 2; // ΛΛ
  477. // if (pdg1 == 2212 && pdg2 == 3222) return 3; // pΣ+
  478. // if (pdg1 == 3122 && pdg2 == 3222) return 4; // ΛΣ+
  479. // if (pdg1 == 3222 && pdg2 == 3222) return 5; // Σ+Σ+
  480. // if (pdg1 == 2212 && pdg2 == 3312) return 6; // pΞ-
  481. // if (pdg1 == 2212 && pdg2 == 3322) return 7; // pΞ0
  482. // if (pdg1 == 3122 && pdg2 == 3312) return 8; // ΛΞ-
  483. // if (pdg1 == 3122 && pdg2 == 3322) return 9; // ΛΞ0
  484. // if (pdg1 == 3222 && pdg2 == 3322) return 10; // Σ+Ξ0
  485. // if (pdg1 == 3222 && pdg2 == 3312) return 11; // Σ+Ξ-
  486. // if (pdg1 == 3212 && pdg2 == 3122) return 12; // Σ0Λ
  487. // if (pdg1 == 2212 && pdg2 == 3212) return 13; // pΣ0
  488. // if (pdg1 == 3212 && pdg2 == 3222) return 14; // Σ0Σ+
  489. // if (pdg1 == 3322 && pdg2 == 3322) return 15; // Ξ0Ξ0
  490. // if (pdg1 == 3312 && pdg2 == 3322) return 16; // Ξ-Ξ0
  491. // if (pdg1 == 3312 && pdg2 == 3312) return 17; // Ξ-Ξ-
  492. // if (pdg1 == 3212 && pdg2 == 3322) return 18; // Σ0Ξ0
  493. // if (pdg1 == 3212 && pdg2 == 3312) return 19; // Σ0Ξ-
  494. // if (pdg1 == 3212 && pdg2 == 3212) return 20; // Σ0Σ0
  495. //
  496. // return -1;
  497. return 0;
  498. }