MpdPlutoGenerator.cxx 4.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164
  1. // -------------------------------------------------------------------------
  2. // ----- MpdPlutoGenerator source file -----
  3. // ----- Created 13/07/04 by V. Friese / D.Bertini -----
  4. // ----- Modified from FairPlutoGenerator for MPD by V. Zhezher -----
  5. // -------------------------------------------------------------------------
  6. #include "MpdPlutoGenerator.h"
  7. #include "FairPrimaryGenerator.h"
  8. #include "TClonesArray.h"
  9. #include "TDatabasePDG.h"
  10. #include "TFile.h"
  11. #include "TLorentzVector.h"
  12. #include "TTree.h"
  13. #include "TVector3.h"
  14. #include "PParticle.h"
  15. #include <iostream>
  16. const Double_t kProtonMass = 0.938271998;
  17. const Double_t kElectronMass = 0.00051098892;
  18. // ----- Default constructor ------------------------------------------
  19. MpdPlutoGenerator::MpdPlutoGenerator()
  20. :FairGenerator(),
  21. iEvent(0),
  22. fEkin(0),
  23. fFileName(NULL),
  24. fInputFile(NULL),
  25. fInputTree(NULL),
  26. fParticles(NULL)
  27. {
  28. }
  29. // ------------------------------------------------------------------------
  30. // ----- Standard constructor -----------------------------------------
  31. MpdPlutoGenerator::MpdPlutoGenerator(const Char_t* fileName, Double_t ekin = 25.0)
  32. :FairGenerator(),
  33. iEvent(0),
  34. fFileName(fileName),
  35. fInputFile(new TFile(fileName)),
  36. fInputTree(NULL),
  37. fParticles(new TClonesArray("PParticle",100))
  38. {
  39. fEkin = ekin;
  40. fInputTree = (TTree*) fInputFile->Get("data");
  41. fInputTree->SetBranchAddress("Particles", &fParticles);
  42. }
  43. // ------------------------------------------------------------------------
  44. // ----- Destructor ---------------------------------------------------
  45. MpdPlutoGenerator::~MpdPlutoGenerator() {
  46. CloseInput();
  47. }
  48. // ------------------------------------------------------------------------
  49. // ----- Public method ReadEvent --------------------------------------
  50. Bool_t MpdPlutoGenerator::ReadEvent(FairPrimaryGenerator* primGen) {
  51. // Check for input file
  52. if ( ! fInputFile ) {
  53. cout << "-E MpdPlutoGenerator: Input file nor open!" << endl;
  54. return kFALSE;
  55. }
  56. // Check for number of events in input file
  57. if ( iEvent > fInputTree->GetEntries() ) {
  58. cout << "-E MpdPlutoGenerator: No more events in input file!" << endl;
  59. CloseInput();
  60. return kFALSE;
  61. }
  62. TFile *g=gFile;
  63. fInputFile->cd();
  64. fInputTree->GetEntry(iEvent++);
  65. g->cd();
  66. // Get PDG database
  67. TDatabasePDG* dataBase = TDatabasePDG::Instance();
  68. // Get number of particle in TClonesrray
  69. Int_t nParts = fParticles->GetEntriesFast();
  70. // Loop over particles in TClonesArray
  71. for (Int_t iPart=0; iPart < nParts; iPart++) {
  72. PParticle* part = (PParticle*) fParticles->At(iPart);
  73. Int_t pdgType = dataBase->ConvertGeant3ToPdg( part->ID() );
  74. // Check if particle type is known to database
  75. if ( ! pdgType ) {
  76. cout << "-W MpdPlutoGenerator: Unknown type " << part->ID()
  77. << ", skipping particle." << endl;
  78. continue;
  79. }
  80. TLorentzVector mom = part->Vect4();
  81. Double_t px = mom.Px();
  82. Double_t py = mom.Py();
  83. Double_t pz = mom.Pz();
  84. TVector3 vertex = part->getVertex();
  85. Double_t vx = vertex.x();
  86. Double_t vy = vertex.y();
  87. Double_t vz = vertex.z();
  88. // ---> Calculate beta and gamma for Lorentztransformation
  89. Double_t eBeam = fEkin + kProtonMass;
  90. Double_t pBeam = TMath::Sqrt(eBeam*eBeam - kProtonMass*kProtonMass);
  91. Double_t betaCM = pBeam / (eBeam + kProtonMass);
  92. Double_t gammaCM = TMath::Sqrt( 1. / ( 1. - betaCM*betaCM) );
  93. //AZ Double_t mass = kElectronMass;
  94. Double_t mass = dataBase->GetParticle(pdgType)->Mass();
  95. Double_t e = sqrt( mass*mass + px*px + py*py + pz*pz );
  96. pz = gammaCM * ( pz - betaCM * e );
  97. // Give track to PrimaryGenerator
  98. //primGen->AddTrack(pdgType, px, py, pz, vx, vy, vz);
  99. //AZ - keep track of parents
  100. Int_t parentIndx = part->GetParentIndex();
  101. primGen->AddTrack(pdgType, px, py, pz, vx, vy, vz, -parentIndx-10);
  102. //AZ
  103. } // Loop over particle in event
  104. return kTRUE;
  105. }
  106. // ------------------------------------------------------------------------
  107. // ----- Public method SkipEvents --------------------------------------
  108. Bool_t MpdPlutoGenerator::SkipEvents(Int_t count) {
  109. if (count<=0) return kTRUE;
  110. iEvent = count;
  111. std::cout << "-I MpdPlutoGenerator: Skipped " << count << " Events!" << std::endl;
  112. return kTRUE;
  113. }
  114. // ----- Private method CloseInput ------------------------------------
  115. void MpdPlutoGenerator::CloseInput() {
  116. if ( fInputFile ) {
  117. cout << "-I MpdPlutoGenerator: Closing input file " << fFileName
  118. << endl;
  119. fInputFile->Close();
  120. delete fInputFile;
  121. }
  122. fInputFile = NULL;
  123. }
  124. // ------------------------------------------------------------------------
  125. ClassImp(MpdPlutoGenerator)