MpdSmashGenerator.cxx 5.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206
  1. // -------------------------------------------------------------------------
  2. // ----- MpdSmashGenerator source file -----
  3. // ----- Adopted by I.Altsybeev 17/07/2020 -----
  4. // ----- SMASH v.>=1.8 with ROOT output is required -----
  5. // -------------------------------------------------------------------------
  6. #include "MpdSmashGenerator.h"
  7. #include "FairPrimaryGenerator.h"
  8. #include "FairMCEventHeader.h"
  9. #include "constants.h"
  10. #include "TRandom.h"
  11. #include "TMath.h"
  12. #include "TString.h"
  13. //#include "TTree.h"
  14. //#include "TChain.h"
  15. #include "TFile.h"
  16. #include "TMCProcess.h"
  17. #include "TObjArray.h"
  18. #include "TPDGCode.h"
  19. #include "TParticle.h"
  20. #include "TVirtualMCStack.h"
  21. #include "TLorentzVector.h"
  22. //#include "TDatabasePDG.h"
  23. //#include "TParticlePDG.h"
  24. #include <iostream>
  25. #include <cstring>
  26. #include <stdio.h>
  27. using namespace std;
  28. using namespace TMath;
  29. // ----- Default constructor ------------------------------------------
  30. MpdSmashGenerator::MpdSmashGenerator()
  31. : FairGenerator()
  32. , fInputFile(NULL)
  33. , fFileName("")
  34. , fImpPar (0)
  35. , fNpart (0)
  36. , fEmptyEv (false)
  37. , fIsRandomRP(false)
  38. , fPsiRP (0)
  39. , fRandom (0)
  40. {
  41. }
  42. // ------------------------------------------------------------------------
  43. // ----- Standard constructor -----------------------------------------
  44. MpdSmashGenerator::MpdSmashGenerator(TString fileName)
  45. : FairGenerator()
  46. , fInputFile(NULL)
  47. , fFileName(fileName)
  48. , fImpPar (0)
  49. , fNpart (0)
  50. , fEmptyEv (false)
  51. , fIsRandomRP(false)
  52. , fPsiRP (0)
  53. , fRandom (0)
  54. {
  55. cout << "-I MpdSmashGenerator: Opening input file " << fFileName << endl;
  56. fInputFile = new TFile(fFileName.Data());
  57. if (!fInputFile) {
  58. Fatal("MpdSmashGenerator", "Cannot open input file.");
  59. exit(1);
  60. }
  61. fInputTree = (TTree*)fInputFile->Get("particles");
  62. // fInputTree = new TChain("out");
  63. // fInputTree->Add(fFileName);
  64. cout << "-I MpdSmashGenerator: Number of entries in tree: " << fInputTree->GetEntries() << endl;
  65. fInputTree->SetBranchAddress( "impact_b", &fImpPar );
  66. fInputTree->SetBranchAddress( "npart", &fNpart );
  67. fInputTree->SetBranchAddress( "empty_event", &fEmptyEv );
  68. fInputTree->SetBranchAddress( "px", fPx );
  69. fInputTree->SetBranchAddress( "py", fPy );
  70. fInputTree->SetBranchAddress( "pz", fPz );
  71. fInputTree->SetBranchAddress( "pdgcode", fPID );
  72. // fInputTree->SetBranchAddress( "charge", fCharge );
  73. fEventNumber = 0;
  74. fPsiRP = 0.;
  75. fIsRandomRP = kTRUE; // by default RP is random
  76. fRandom = new TRandom3();
  77. fRandom->SetSeed(0);
  78. // fPNCorr = 0.; // by default no correction
  79. }
  80. // ------------------------------------------------------------------------
  81. // ----- Destructor ---------------------------------------------------
  82. MpdSmashGenerator::~MpdSmashGenerator()
  83. {
  84. delete fInputFile;
  85. delete fRandom;
  86. }
  87. // ------------------------------------------------------------------------
  88. // ----- Public method ReadEvent --------------------------------------
  89. Bool_t MpdSmashGenerator::ReadEvent(FairPrimaryGenerator* primGen)
  90. {
  91. // Check for input file
  92. if (!fInputFile) {
  93. cout << "-E MpdSmashGenerator: Input file not open! " << endl;
  94. return kFALSE;
  95. }
  96. // Check for primary generator
  97. if (!primGen) {
  98. cout << "-E- MpdSmashGenerator::ReadEvent: "
  99. << "No PrimaryGenerator!" << endl;
  100. return kFALSE;
  101. }
  102. // Skip empty SMASH events:
  103. // bool hasInteraction = false;
  104. // while( !hasInteraction )
  105. // {
  106. // fInputTree->GetEntry(fEventNumber);
  107. // if ( fEmptyEv )
  108. // {
  109. // fEventNumber++;
  110. // continue;
  111. // }
  112. // hasInteraction = true;
  113. // }
  114. fInputTree->GetEntry(fEventNumber);
  115. cout << "-I MpdSmashGenerator: Event " << fEventNumber << ", b = " << fImpPar
  116. << " fm, multiplicity " << fNpart << endl;// << ", Elab: " << energ << endl;
  117. if ( fNpart > MAX_N_PART )
  118. {
  119. cout <<"-E- MpdSmashGenerator: ERROR! fNpart>MAX_N_PART (fNpart=" << fNpart << ")" << endl;
  120. exit(1);
  121. }
  122. /* set random Reaction Plane angle */
  123. fPsiRP = 0;
  124. if( fIsRandomRP )
  125. fPsiRP = fRandom->Uniform(2.0*TMath::Pi());
  126. // Set event id and impact parameter in MCEvent if not yet done
  127. FairMCEventHeader* event = primGen->GetEvent();
  128. if (event && (!event->IsSet()))
  129. {
  130. //event->SetEventID(evnr);
  131. event->SetB(fImpPar);
  132. event->SetRotZ(fPsiRP);
  133. event->MarkSet(kTRUE);
  134. }
  135. // Loop over tracks in the current event
  136. for (int itrack = 0; itrack < fNpart; itrack++)
  137. {
  138. Double_t px = fPx[itrack];
  139. Double_t py = fPy[itrack];
  140. Double_t pz = fPz[itrack];
  141. Int_t pid = fPID[itrack];
  142. /* rotate RP angle */
  143. if (fPsiRP!=0.)
  144. {
  145. Double_t cosPsi = TMath::Cos(fPsiRP);
  146. Double_t sinPsi = TMath::Sin(fPsiRP);
  147. Double_t px1 = px*cosPsi-py*sinPsi;
  148. Double_t py1 = px*sinPsi+py*cosPsi;
  149. px = px1;
  150. py = py1;
  151. }
  152. // Give track to PrimaryGenerator
  153. primGen->AddTrack( pid, px,py,pz, 0.0, 0.0, 0.0); // pdg, mom [GeV/c], vertex [cm]
  154. }
  155. fEventNumber++;
  156. return kTRUE;
  157. }
  158. // ------------------------------------------------------------------------
  159. ClassImp(MpdSmashGenerator);