MpdTPythia8Generator.cxx 6.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177
  1. /**
  2. * MpdTPythia8Generator.cxx
  3. * *@author A.Zinchenko <Alexander.Zinchenko@jinr.ru>
  4. * The MpdTPythia8Generator produces particles from TPythia8 output tree.
  5. **/
  6. #include "MpdTPythia8Generator.h"
  7. #include "FairMCEventHeader.h"
  8. #include "FairPrimaryGenerator.h"
  9. #include "TClonesArray.h"
  10. #include "TDatabasePDG.h"
  11. #include "TFile.h"
  12. #include "TLorentzVector.h"
  13. #include "TMath.h"
  14. //#include "TMCParticle.h"
  15. #include "TParticle.h"
  16. #include "TPythia8.h"
  17. #include "TRandom.h"
  18. #include "TSystem.h"
  19. #include "TTree.h"
  20. #include "TVirtualMC.h"
  21. #include <fstream>
  22. #include <iostream>
  23. #include <iomanip>
  24. // #define debug_hsd
  25. using namespace std;
  26. // ------------------------------------------------------------------------
  27. MpdTPythia8Generator::MpdTPythia8Generator() : FairGenerator(),
  28. fInputFile(NULL), fTree(NULL), fParticles(NULL), fNevents(0), fEdit(0)
  29. {
  30. // Constructor
  31. }
  32. // ------------------------------------------------------------------------
  33. MpdTPythia8Generator::MpdTPythia8Generator(TString fileName, TString treeName, TString branchName)
  34. : FairGenerator(),
  35. fParticles(new TClonesArray("TParticle",100)),
  36. //fParticles(new TClonesArray("TMCParticle",100)),
  37. fNevents(0),
  38. fEdit(0)
  39. {
  40. // Constructor
  41. //gSystem->Setenv("PYTHIA8DATA", "/home/zinchenk/pythia8/pythia8215/share/Pythia8/xmldoc");
  42. //TPythia8 pyt8; // to update PDG data base
  43. AddParticlesToPdgDataBase();
  44. fInputFile = new TFile(fileName);
  45. fTree = (TTree*) fInputFile->Get(treeName.Data());
  46. fTree->SetBranchAddress(branchName.Data(),&fParticles);
  47. cout << "\n ***** MpdPythia8Generator ***** " << endl;
  48. cout << " Number of events in input file: " << fileName << " " << fTree->GetEntries() << endl;
  49. cout << " ***** ***** \n" << endl;
  50. }
  51. // ------------------------------------------------------------------------
  52. MpdTPythia8Generator::~MpdTPythia8Generator()
  53. {
  54. // Destructor
  55. fParticles->Delete();
  56. fInputFile->Close();
  57. delete fInputFile;
  58. }
  59. //-------------------------------------------------------------------------
  60. Bool_t MpdTPythia8Generator::ReadEvent(FairPrimaryGenerator* primGen)
  61. {
  62. // Generate event
  63. static Int_t first = 1;
  64. if (first) {
  65. first = 0;
  66. //cout << " histo " << fHist << " " << fHist->GetNbinsX() << " " << fHist->GetNbinsY() << endl;
  67. //if (TDatabasePDG::Instance()->GetParticle(fPdgCode) == NULL) {
  68. // Define particle with fPdgCode (take it from DCM-QGSM input file)
  69. // To be done ...
  70. //}
  71. }
  72. //
  73. fTree->GetEntry(fNevents++);
  74. Int_t nPart = fParticles->GetEntriesFast(), iFirstIndx = 0;
  75. cout << " Number of particles in event " << fNevents-1 << ": " << nPart << endl;
  76. //TMCParticle *part = NULL;
  77. TParticle *part = NULL;
  78. Bool_t track = kTRUE;
  79. for (Int_t ipart = 0; ipart < nPart; ++ipart) {
  80. //part = (TMCParticle*) fParticles->UncheckedAt(ipart);
  81. part = (TParticle*) fParticles->UncheckedAt(ipart);
  82. //if (fEdit == 0 && part->GetKS() == 21) { ++iFirstIndx; continue; }
  83. //track = (part->GetKS() == 1) ? kTRUE : kFALSE;
  84. track = (part->GetStatusCode() == 1) ? kTRUE : kFALSE;
  85. //FairMCEventHeader* eventHeader = primGen->GetEvent();
  86. /* Set event impact parameter in MCEvent if not yet done */
  87. /*
  88. if (eventHeader && (!eventHeader->IsSet()) ) {
  89. eventHeader->SetB(fB);
  90. eventHeader->MarkSet(kTRUE);
  91. }
  92. */
  93. //Int_t parent = 0; //TMath::Max (part->GetParent() - 1 - iFirstIndx, -1);
  94. //Int_t parent = TMath::Max (part->GetFirstMother() - 1 - iFirstIndx, -1);
  95. Int_t parent = part->GetFirstMother();
  96. //cout << ipart << " " << part->GetStatusCode() << " " << part->GetPdgCode() << " " << part->GetFirstMother() << " " << parent << endl;
  97. //if (parent == -1) continue; // do not add beam protons - 5.12.2017
  98. primGen->AddTrack(part->GetPdgCode(), part->Px(), part->Py(), part->Pz(),
  99. part->Vx()*0.1, part->Vy()*0.1, part->Vz()*0.1, parent, track);
  100. //part->Vx(), part->Vy(), part->Vz(), -TMath::Abs(parent)-1, track); //AZ 5.12.2017
  101. }
  102. return kTRUE;
  103. }
  104. // ------------------------------------------------------------------------
  105. void MpdTPythia8Generator::AddParticlesToPdgDataBase()
  106. {
  107. // Add some pythia specific particle code to the data base - from TPythia8
  108. TDatabasePDG *pdgDB = TDatabasePDG::Instance();
  109. pdgDB->AddParticle("psi(3770)","psi(3770)", 3.77315, kTRUE,
  110. 0, 0, "psi(3770)", 30443);
  111. pdgDB->AddParticle("f_0(980)","f_0(980)", 1.00000, kTRUE,
  112. 0, 3, "f_0 meson", 9010221);
  113. pdgDB->AddParticle("Pomeron","Pomeron", 0.00000, kTRUE,
  114. 0, 0, "Pomeron", 990);
  115. pdgDB->AddParticle("p_diffr+","p_diffr+", 0.0, kTRUE,
  116. 0, 0, "p_diffr+", 9902210);
  117. pdgDB->AddParticle("J/psi[3S1(8)]","J/psi[3S1(8)]", 3.29692, kTRUE,
  118. 0, 0, "Ccbar multiplet", 9940003);
  119. pdgDB->AddParticle("chi_2c[3S1(8)]","chi_2c[3S1(8)]", 3.75620, kTRUE,
  120. 0, 0, "Ccbar multiplet", 9940005);
  121. pdgDB->AddParticle("chi_0c[3S1(8)]","chi_0c[3S1(8)]", 3.61475, kTRUE,
  122. 0, 0, "Ccbar multiplet", 9940011);
  123. pdgDB->AddParticle("chi_1c[3S1(8)]","chi_1c[3S1(8)]", 3.71066, kTRUE,
  124. 0, 0, "Ccbar multiplet", 9940023);
  125. pdgDB->AddParticle("psi(2S)[3S1(8)]","psi(2S)[3S1(8)]", 3.88611, kTRUE,
  126. 0, 0, "QCD diffr. state", 9940103);
  127. pdgDB->AddParticle("J/psi[1S0(8)]","J/psi[1S0(8)]", 3.29692, kTRUE,
  128. 0, 0, "Ccbar multiplet", 9941003);
  129. pdgDB->AddParticle("psi(2S)[1S0(8)]","psi(2S)[1S0(8)]", 3.88611, kTRUE,
  130. 0, 0, "Ccbar multiplet", 9941103);
  131. pdgDB->AddParticle("J/psi[3PJ(8)]","J/psi[3PJ(8)]", 3.29692, kTRUE,
  132. 0, 0, "Ccbar multiplet", 9942003);
  133. pdgDB->AddParticle("psi(3770)[3PJ(8)]","psi(3770)[3PJ(8)]", 3.97315, kTRUE,
  134. 0, 0, "Ccbar multiplet", 9942033);
  135. pdgDB->AddParticle("psi(2S)[3PJ(8)]","psi(2S)[3PJ(8)]", 3.88611, kTRUE,
  136. 0, 0, "Ccbar multiplet", 9942103);
  137. pdgDB->AddParticle("Upsilon[3S1(8)]","Upsilon[3S1(8)]", 9.6603, kTRUE,
  138. 0, 0, "Bbbar multiplet", 9950003);
  139. pdgDB->AddParticle("chi_2b[3S1(8)]","chi_2b[3S1(8)]", 10.1122, kTRUE,
  140. 0, 0, "Bbbar multiplet", 9950005);
  141. pdgDB->AddParticle("Upsilon[1S0(8)]","Upsilon[1S0(8)]", 9.6603, kTRUE,
  142. 0, 0, "Bbbar multiplet", 9951003);
  143. pdgDB->AddParticle("Upsilon(2S)[1S0(8)]","Upsilon(2S)[1S0(8)]", 10.22326, kTRUE,
  144. 0, 0, "Bbbar multiplet", 9951103);
  145. pdgDB->AddParticle("Upsilon(2S)[3PJ(8)]","Upsilon(2S)[3PJ(8)]", 10.22326, kTRUE,
  146. 0, 0, "Bbbar multiplet", 9952103);
  147. }
  148. // ------------------------------------------------------------------------