MpdRoInvMassTask.cxx 4.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166
  1. // skelet
  2. #ifndef ROOT_MpdRoInvMassTask
  3. #include "MpdRoInvMassTask.h"
  4. #endif
  5. #include <TLorentzVector.h>
  6. #include "FairRootManager.h"
  7. #include "MpdMCTrack.h"
  8. #include "MpdTrack.h"
  9. #include <iostream>
  10. using std::cout;
  11. using std::endl;
  12. // ----- Default constructor ---------------------------------------
  13. MpdRoInvMassTask::MpdRoInvMassTask():
  14. FairTask(),
  15. fEventCounter(0)
  16. {}
  17. // ----- constructor with names ------------------------------------
  18. MpdRoInvMassTask::MpdRoInvMassTask(const char *name, const char *title):
  19. FairTask(name),
  20. fEventCounter(0)
  21. {}
  22. // ----- Destructor -----------------------------------------------
  23. MpdRoInvMassTask::~MpdRoInvMassTask()
  24. {}
  25. // -------------------------------------------------------------------
  26. InitStatus MpdRoInvMassTask::Init()
  27. {
  28. cout<<"\n-I- [MpdRoInvMassTask::Init] " <<endl;
  29. FairRootManager *manager = FairRootManager::Instance();
  30. fMCTracks = (TClonesArray *) manager->GetObject("MCTrack");
  31. fDstEvent = (MpdEvent *) manager->GetObject("MPDEvent.");
  32. Register();
  33. fRoInvMass = new TH1F("RoInvMass","Invariant Mass: (#pi_{+}#pi_{-}) 90%CL",200,0.,2.);
  34. fRoInvMass->GetXaxis()->SetTitle("M_{inv}");
  35. fRoInvMassMC = new TH1F("RoInvMassMC","Invariant Mass (MC): (#pi_{+}#pi_{-}) 90%CL",200,0.,2.);
  36. fRoInvMassMC->GetXaxis()->SetTitle("M_{inv}");
  37. fRoInvMassMC->SetLineColor(kRed);
  38. fPDG = TDatabasePDG::Instance();
  39. return kSUCCESS;
  40. }
  41. // -------------------------------------------------------------------
  42. void MpdRoInvMassTask::Exec(Option_t * option)
  43. {
  44. fEventCounter++;
  45. cout<<"-I- [MpdRoInvMassTask::Exec] " << "{" << fEventCounter << "}" <<endl;
  46. //MpdEvent *mpdEvent = fDstEvent; // one event by one event
  47. //event->Dump();
  48. TClonesArray *mpdTracks = fDstEvent->GetGlobalTracks();
  49. Int_t nTracks = mpdTracks->GetEntriesFast();
  50. cout << "N of Reconstructed tracks = " << nTracks << endl;
  51. //mpdTracks->Dump();
  52. /*
  53. * Let us calculate invariant mass of rho0 meson (rho0-> pi+ + pi-)
  54. */
  55. TLorentzVector lPos, lNeg;
  56. /* Events loop */
  57. for (Int_t i = 0; i < nTracks; i++)
  58. {
  59. MpdTrack *track1 = (MpdTrack *) mpdTracks->UncheckedAt(i);
  60. if ( !track1->GetTofFlag() ) continue; // no tof identification
  61. Float_t pion1 = track1->GetPidProbPion();
  62. if ( !(pion1 > 0.9) ) continue; // accept only pi+ 90% CL, reject other
  63. //track1->Dump();
  64. lPos.SetXYZM( track1->GetPx(), track1->GetPy(), track1->GetPz(), fPDG->GetParticle(211)->Mass() );
  65. for (Int_t j = 0; j<i; j++)
  66. {
  67. MpdTrack *track2 = (MpdTrack *) mpdTracks->UncheckedAt(j);
  68. if ( !track2->GetTofFlag() ) continue; // no tof identification
  69. Float_t pion2 = track2->GetPidProbPion();
  70. if( !(pion2 < 0.9) ) continue; // accept only pi- 90% CL, reject other
  71. //track2->Dump();
  72. lNeg.SetXYZM( track2->GetPx(), track2->GetPy(), track2->GetPz(), fPDG->GetParticle(-211)->Mass() );
  73. Float_t minv = (lPos + lNeg).Mag();
  74. fRoInvMass->Fill(minv);
  75. }
  76. }
  77. /* End of events loop */
  78. /* Events loop. The same but for MC */
  79. Int_t nMCTracks = fMCTracks->GetEntriesFast();
  80. cout << "N of MC tracks = " << nMCTracks << endl;
  81. for (Int_t i = 0; i < nMCTracks; i++)
  82. {
  83. MpdMCTrack *MCtrack1 = (MpdMCTrack*) fMCTracks->UncheckedAt(i);
  84. //if( MCtrack1->GetMotherId() >0 ) continue; // start from initial vertex
  85. Int_t pdgCode1 = MCtrack1->GetPdgCode();
  86. if ( pdgCode1!=211 ) continue; // pi+
  87. MCtrack1->Get4Momentum(lPos);
  88. for (Int_t j = 0; j<i; j++)
  89. {
  90. MpdMCTrack *MCtrack2 = (MpdMCTrack*) fMCTracks->UncheckedAt(j);
  91. //if( MCtrack2->GetMotherId() >0 ) continue; // start from initial vertex
  92. Int_t pdgCode2 = MCtrack2->GetPdgCode();
  93. if ( pdgCode2!=-211 ) continue; // pi-
  94. MCtrack2->Get4Momentum(lNeg);
  95. Float_t minv = (lPos + lNeg).Mag();
  96. fRoInvMassMC->Fill(minv); // with probability weight
  97. }
  98. }
  99. /* End of MC events loop */
  100. }
  101. // -------------------------------------------------------------------
  102. void MpdRoInvMassTask::Reset()
  103. {}
  104. // -------------------------------------------------------------------
  105. void MpdRoInvMassTask::Finish()
  106. {
  107. //cout<<"\n-I- [MpdRoInvMassTask::Finish] "<< endl;
  108. //TFile fileOut("MinvRho0.root","recreate");
  109. Double_t scale = 1./ (Double_t) fEventCounter;
  110. fRoInvMass->Scale(scale); // scale to total number of events
  111. fRoInvMassMC->Scale(scale); // scale to total number of events
  112. fRoInvMass->Write("");
  113. fRoInvMassMC->Write("");
  114. }
  115. // -------------------------------------------------------------------
  116. void MpdRoInvMassTask::Register()
  117. {
  118. //FairRootManager::Instance()->Register("MCTrack", "MC", fMCTracks, kTRUE);
  119. }
  120. // -------------------------------------------------------------------
  121. ClassImp(MpdRoInvMassTask);