MpdUnigenGenerator.cxx 5.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194
  1. #include "MpdUnigenGenerator.h"
  2. MpdUnigenGenerator::MpdUnigenGenerator() :
  3. FairGenerator(),
  4. fEventNumber(0),
  5. fInFile(nullptr),
  6. fInTree(nullptr),
  7. fEvent(nullptr),
  8. fParticle(nullptr),
  9. fEventPlaneSet(kFALSE),
  10. fSpectatorsON(kFALSE),
  11. fPhiMin(0.),
  12. fPhiMax(0.) {
  13. }
  14. MpdUnigenGenerator::MpdUnigenGenerator(TString fileName, Bool_t isSpectator) :
  15. FairGenerator(),
  16. fEventNumber(0),
  17. fInFile(nullptr),
  18. fInTree(nullptr),
  19. fEvent(nullptr),
  20. fParticle(nullptr),
  21. fEventPlaneSet(kFALSE),
  22. fPhiMin(0.),
  23. fPhiMax(0.) {
  24. std::cout << "-I MpdUnigenGenerator: Opening input file " << fileName.Data() << std::endl;
  25. fInFile = new TFile(fileName,"read");
  26. if (!fInFile){
  27. Fatal("MpdUnigenGenerator", "Cannot open input file.");
  28. exit(1);
  29. }
  30. fInTree = (TTree*) fInFile->Get("events");
  31. if (!fInTree){
  32. Fatal("MpdUnigenGenerator", "Cannot open TTree from the file.");
  33. exit(1);
  34. }
  35. fSpectatorsON = isSpectator;
  36. if (fSpectatorsON)
  37. std::cout << "-I MpdUnigenGenerator: Spectators ON" << std::endl;
  38. Long64_t nEntries = fInTree->GetEntries();
  39. std::cout << "-I MpdUnigenGenerator: Number of entries in tree: " << nEntries << std::endl;
  40. fInTree->SetBranchAddress("event", &fEvent);
  41. if (fSpectatorsON){
  42. RegisterIons();
  43. }
  44. }
  45. Bool_t MpdUnigenGenerator::ReadEvent(FairPrimaryGenerator* primGen){
  46. // Check for input file
  47. if (!fInFile) {
  48. cout << "-E- MpdUnigenGenerator: Input file not open! " << endl;
  49. return kFALSE;
  50. }
  51. // Check for primary generator
  52. if (!primGen) {
  53. cout << "-E- MpdUnigenGenerator::ReadEvent: "
  54. << "No PrimaryGenerator!" << endl;
  55. return kFALSE;
  56. }
  57. fInTree->GetEntry(fEventNumber);
  58. std::cout << "-I- MpdUnigenGenerator: Event " << fEventNumber << std::endl;
  59. Double_t phi = fEvent->GetPhi();
  60. Double_t dphi = 0.;
  61. // ---> Generate rotation angle
  62. if (fEventPlaneSet) {
  63. gRandom->SetSeed(0);
  64. dphi = gRandom->Uniform(fPhiMin, fPhiMax);
  65. phi += dphi;
  66. }
  67. FairMCEventHeader* header = primGen->GetEvent();
  68. if (header && (!header->IsSet())) {
  69. header->SetEventID(fEvent->GetEventNr());
  70. header->SetB(fEvent->GetB());
  71. header->MarkSet(kTRUE);
  72. header->SetRotZ(phi);
  73. }
  74. UInt_t nTracks = fEvent->GetNpa();
  75. Int_t ionPdg;
  76. for (Int_t iTrack=0; iTrack<nTracks; iTrack++){
  77. fParticle = fEvent->GetParticle(iTrack);
  78. if (!fParticle)
  79. continue;
  80. Double_t px = fParticle->Px();
  81. Double_t py = fParticle->Py();
  82. if (fEventPlaneSet) {
  83. Double_t pt = TMath::Sqrt(px * px + py * py);
  84. Double_t azim = TMath::ATan2(py, px);
  85. azim += dphi;
  86. px = pt * TMath::Cos(azim);
  87. py = pt * TMath::Sin(azim);
  88. }
  89. if (fParticle->GetPdg() < 1e9) {
  90. primGen->AddTrack(fParticle->GetPdg(), px, py, fParticle->Pz(), 0., 0., 0.);
  91. } else {
  92. // Since hyper-nuclei are not (yet) supported by FairRoot, their PDG
  93. // is replaced by that of the non-strange analogue.
  94. ionPdg = fParticle->GetPdg();
  95. if (GetIonLambdas(fParticle->GetPdg())) {
  96. ionPdg = fParticle->GetPdg() - GetIonLambdas(fParticle->GetPdg()) * kPdgLambda;
  97. std::cout << "-W- MpdUnigenGenerator: Replacing hypernucleus (PDG " << fParticle->GetPdg() << ") by PDG " << ionPdg << std::endl;
  98. }
  99. // Charged ions can be registered
  100. if (GetIonCharge(ionPdg)) {
  101. primGen->AddTrack(ionPdg, px, py, fParticle->Pz(), 0., 0., 0.);
  102. } else {
  103. // Neutral ions are not supported by GEANT4.
  104. // They are thus decomposed into neutrons (as an approximation)
  105. Int_t mass = GetIonMass(ionPdg);
  106. for (Int_t iNeutron = 0; iNeutron < mass; iNeutron++) {
  107. Double_t pxNeutron = px/(Double_t)mass;
  108. Double_t pyNeutron = py/(Double_t)mass;
  109. Double_t pzNeutron = fParticle->Pz()/(Double_t)mass;
  110. Double_t eNeutron = fParticle->E()/(Double_t)mass;
  111. primGen->AddTrack(2112, pxNeutron, pyNeutron, pzNeutron, 0., 0., 0.);
  112. }
  113. std::cout << "-W- MpdUnigenGenerator: Neutral fragment (PDG " << ionPdg << ") is split into " << mass << " neutrons." << std::endl;
  114. } // Neutral ion ?
  115. } // ion ?
  116. }
  117. fEventNumber++;
  118. return kTRUE;
  119. }
  120. MpdUnigenGenerator::~MpdUnigenGenerator(){
  121. fInFile->Close();
  122. delete fInTree;
  123. delete fEvent;
  124. delete fParticle;
  125. }
  126. Int_t MpdUnigenGenerator::RegisterIons(void) {
  127. std::cout << "-I- MpdUnigenGenerator: Registering ions..." << std::endl;
  128. UParticle* particle {nullptr};
  129. Int_t nIons {0};
  130. fIonMap.clear();
  131. for (Int_t iEvent=0; iEvent<fInTree->GetEntries(); iEvent++) {
  132. fInTree->GetEntry(iEvent);
  133. for (Int_t iParticle=0; iParticle<fEvent->GetNpa(); iParticle++) {
  134. particle = fEvent->GetParticle(iParticle);
  135. Int_t pdgCode = particle->GetPdg();
  136. if (pdgCode > 1e9) { // ion
  137. // For ions the PDG code is +-10LZZZAAAI,
  138. // where A = n_Lambda + n_proton + n_neutron, Z = n_proton, L = n_Lambda
  139. // and I = 0 - ground state, I>0 - excitations
  140. const Int_t nLambda = GetIonLambdas(pdgCode);
  141. const Int_t chrg = GetIonCharge(pdgCode);
  142. const Int_t mass = GetIonMass(pdgCode);
  143. // Neutral ions are skipped (not present in G4IonTable)
  144. if (chrg == 0) continue;
  145. // Add ion to ion map
  146. TString ionName = Form("Ion_%d_%d_%d", mass, chrg, nLambda);
  147. if (fIonMap.find(ionName) == fIonMap.end()) { // new ion
  148. fIonMap[ionName] = new FairIon(ionName, chrg, mass, chrg);
  149. std::cout << "-I- MpdUnigenGenerator: Added new ion with PDG " << pdgCode
  150. << " (Charge " << chrg << ", Mass " << mass << ", nLambda " << nLambda
  151. << ") from event " << iEvent << " at index " << iParticle << std::endl;
  152. nIons++;
  153. } //new ion ?
  154. } // ion ?
  155. } // particle loop
  156. } // event loop
  157. // Adding new ions to FairRunSim instance
  158. for (const auto& entry : fIonMap) {
  159. FairRunSim::Instance()->AddNewIon(entry.second);
  160. }
  161. return fIonMap.size();
  162. }
  163. ClassImp(MpdUnigenGenerator);