MpdMcDstGenerator.cxx 3.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117
  1. #include "MpdMcDstGenerator.h"
  2. TClonesArray** McDst::mcArrays = nullptr;
  3. MpdMcDstGenerator::MpdMcDstGenerator() :
  4. fEventNumber(0),
  5. myReader(nullptr),
  6. fGenTracks(nullptr),
  7. fEventPlaneSet(kFALSE),
  8. fPhiMin(0.),
  9. fPhiMax(0.) {
  10. }
  11. MpdMcDstGenerator::MpdMcDstGenerator(TString fileName) :
  12. fEventNumber(0),
  13. myReader(nullptr),
  14. fGenTracks(nullptr),
  15. fEventPlaneSet(kFALSE),
  16. fPhiMin(0.),
  17. fPhiMax(0.) {
  18. cout << "-I MpdMcDstGenerator: Opening input file " << fileName << endl;
  19. myReader = new McDstReader(fileName.Data());
  20. myReader->Init();
  21. myReader->setStatus("*", 0);
  22. myReader->setStatus("Event", 1);
  23. myReader->setStatus("Particle", 1);
  24. if (!myReader->chain()) {
  25. cout << "No chain has been found." << endl;
  26. return;
  27. }
  28. Long64_t eventsInTree = myReader->tree()->GetEntries();
  29. cout << "eventsInTree: " << eventsInTree << endl;
  30. Long64_t events2read = myReader->chain()->GetEntries();
  31. cout << "Number of events to read: " << events2read << endl;
  32. // Create structure to store generator tracks ...
  33. FairRunSim* runSim = FairRunSim::Instance();
  34. MpdGenTrackTask* tracksInfo = new MpdGenTrackTask();
  35. runSim->AddTask(tracksInfo);
  36. fGenTracks = tracksInfo->GetTracksInfo();
  37. }
  38. Bool_t MpdMcDstGenerator::ReadEvent(FairPrimaryGenerator* primGen) {
  39. fGenTracks->Delete(); // Clear array with gen. tracks ...
  40. myReader->chain()->GetEntry(fEventNumber);
  41. McDst* dst = myReader->mcDst();
  42. McEvent* event = dst->event();
  43. if (!event) {
  44. cout << "Something went wrong when getting event, stopping here ..." << endl;
  45. return kFALSE;
  46. }
  47. Double_t phi = 0.;
  48. // ---> Generate rotation angle
  49. if (fEventPlaneSet) {
  50. gRandom->SetSeed(0);
  51. phi = gRandom->Uniform(fPhiMin, fPhiMax);
  52. }
  53. FairMCEventHeader* header = primGen->GetEvent();
  54. if (header && (!header->IsSet())) {
  55. header->SetEventID(event->eventNr());
  56. header->SetNPrim(event->npart());
  57. header->SetB(event->impact());
  58. header->MarkSet(kTRUE);
  59. header->SetRotZ(phi);
  60. }
  61. UInt_t nTracks = dst->numberOfParticles();
  62. for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
  63. McParticle* particle = dst->particle(iTrack);
  64. if (!particle)
  65. continue;
  66. Double_t px = particle->px();
  67. Double_t py = particle->py();
  68. if (fEventPlaneSet) {
  69. Double_t pt = Sqrt(px * px + py * py);
  70. Double_t azim = ATan2(py, px);
  71. azim += phi;
  72. px = pt * Cos(azim);
  73. py = pt * Sin(azim);
  74. }
  75. primGen->AddTrack(particle->pdg(), px, py, particle->pz(), 0., 0., 0.);
  76. // Push a new track to the array of gen. tracks ...
  77. MpdGenTrack* track = new ((*fGenTracks)[fGenTracks->GetEntriesFast()]) MpdGenTrack();
  78. track->SetXYZT(px, py, particle->z(), particle->t());
  79. track->SetPxyz(px, py, particle->pz());
  80. track->SetPdg(particle->pdg());
  81. track->SetE(particle->e());
  82. }
  83. fEventNumber++;
  84. return kTRUE;
  85. }
  86. MpdMcDstGenerator::~MpdMcDstGenerator() {
  87. myReader->Finish();
  88. delete fGenTracks;
  89. }
  90. ClassImp(MpdMcDstGenerator);