Mpd3fdGenerator.cxx 6.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209
  1. // -------------------------------------------------------------------------
  2. // ----- Mpd3fdGenerator source file -----
  3. // ----- Created 24/06/04 by V. Friese -----
  4. // -------------------------------------------------------------------------
  5. #include "Mpd3fdGenerator.h"
  6. #include "FairPrimaryGenerator.h"
  7. #include "FairMCEventHeader.h"
  8. #include "constants.h"
  9. #include "TMCProcess.h"
  10. #include "TObjArray.h"
  11. #include "TPDGCode.h"
  12. #include "TParticle.h"
  13. #include "TRandom.h"
  14. #include "TString.h"
  15. #include "TVirtualMCStack.h"
  16. #include "TLorentzVector.h"
  17. #include "TDatabasePDG.h"
  18. #include "TParticlePDG.h"
  19. #include "TMath.h"
  20. #include "TString.h"
  21. #include "TRegexp.h"
  22. #include <iostream>
  23. #include <cstring>
  24. #include <stdio.h>
  25. using namespace std;
  26. using namespace TMath;
  27. // ----- Default constructor ------------------------------------------
  28. Mpd3fdGenerator::Mpd3fdGenerator()
  29. : FairGenerator(),
  30. fInputFile(NULL),
  31. fFileName("") {
  32. }
  33. // ------------------------------------------------------------------------
  34. // ----- Standard constructor -----------------------------------------
  35. Mpd3fdGenerator::Mpd3fdGenerator(TString fileName)
  36. : FairGenerator(),
  37. fInputFile(NULL),
  38. fFileName(fileName) {
  39. // fFileName = fileName;
  40. cout << "-I Mpd3fdGenerator: Opening input file " << fFileName << endl;
  41. fInputFile = new TFile(fFileName.Data());
  42. if (!fInputFile) {
  43. Fatal("Mpd3fdGenerator", "Cannot open input file.");
  44. exit(1);
  45. }
  46. fDstTree = new TChain("out");
  47. fDstTree->Add(fFileName);
  48. // Deactivate All
  49. fDstTree->SetBranchStatus("*",0);
  50. // Activate only selected branches
  51. fDstTree->SetBranchStatus("npart",1);
  52. fDstTree->SetBranchStatus("px",1);
  53. fDstTree->SetBranchStatus("py",1);
  54. fDstTree->SetBranchStatus("pz",1);
  55. //fDstTree->SetBranchStatus("x",1);
  56. //fDstTree->SetBranchStatus("y",1);
  57. //fDstTree->SetBranchStatus("z",1);
  58. fDstTree->SetBranchStatus("E",1);
  59. fDstTree->SetBranchStatus("id",1);
  60. fDstTree->SetBranchAddress("px", fPx);
  61. fDstTree->SetBranchAddress("py", fPy);
  62. fDstTree->SetBranchAddress("pz", fPz);
  63. //fDstTree->SetBranchAddress("x", fX); // [fm]
  64. //fDstTree->SetBranchAddress("y", fY); // [fm]
  65. //fDstTree->SetBranchAddress("z", fZ); // [fm]
  66. fDstTree->SetBranchAddress("E", fE);
  67. fDstTree->SetBranchAddress("npart", &fNpart);
  68. fDstTree->SetBranchAddress("id", fPID);
  69. fEventNumber = 0;
  70. fPsiRP=0.;
  71. fisRP=kTRUE; // by default RP is random
  72. frandom = new TRandom2();
  73. frandom->SetSeed(0);
  74. fPNCorr = 0.; // by default no correction
  75. }
  76. // ------------------------------------------------------------------------
  77. // ----- Destructor ---------------------------------------------------
  78. Mpd3fdGenerator::~Mpd3fdGenerator() {
  79. delete fInputFile;
  80. delete fDstTree;
  81. delete frandom;
  82. }
  83. // ------------------------------------------------------------------------
  84. // ----- Public method ReadEvent --------------------------------------
  85. Bool_t Mpd3fdGenerator::ReadEvent(FairPrimaryGenerator* primGen) {
  86. // ---> Check for input file
  87. if (!fInputFile) {
  88. cout << "-E Mpd3fdGenerator: Input file not open! " << endl;
  89. return kFALSE;
  90. }
  91. // ---> Check for primary generator
  92. if (!primGen) {
  93. cout << "-E- Mpd3fdGenerator::ReadEvent: "
  94. << "No PrimaryGenerator!" << endl;
  95. return kFALSE;
  96. }
  97. fDstTree->GetEntry(fEventNumber);
  98. //Int_t events = fDstTree->GetEntries();
  99. // ---> Define event variables to be read from file
  100. Int_t energ, impact;
  101. const TRegexp elb("_elb[0-9]*gev_");
  102. TSubString subelb = fFileName(elb); // select substring by regexp "elab"
  103. sscanf(subelb.Data(), "_elb%dgev_", &energ);
  104. const TRegexp imp("_[0-9]*fm_");
  105. TSubString subimp = fFileName(imp); // select substring by regexp "impact"
  106. sscanf(subelb.Data(), "_%dfm_", &impact);
  107. Float_t b = (Float_t) impact;
  108. /*
  109. // ---> Calculate beta and gamma for Lorentztransformation
  110. TDatabasePDG* pdgDB = TDatabasePDG::Instance();
  111. TParticlePDG* kProton = pdgDB->GetParticle(2212);
  112. Double_t kProtonMass = kProton->Mass();
  113. Double_t eBeam = energ + kProtonMass;
  114. Double_t pBeam = TMath::Sqrt(eBeam * eBeam - kProtonMass * kProtonMass);
  115. Double_t betaCM = pBeam / (eBeam + kProtonMass);
  116. Double_t gammaCM = TMath::Sqrt(1. / (1. - betaCM * betaCM));
  117. */
  118. cout << "-I Mpd3fdGenerator: Event " << fEventNumber << ", b = " << b
  119. << " fm, multiplicity " << fNpart << ", Elab: " << energ << endl;
  120. if (fNpart>kBatyukConst) {cout <<"-E- 3fdGenerator SELFCHECK ERROR"<< endl; exit(1);}
  121. /* set random Reaction Plane angle */
  122. if (fisRP) {fPsiRP=frandom->Uniform(2.0*TMath::Pi());}
  123. // Set event id and impact parameter in MCEvent if not yet done
  124. FairMCEventHeader* event = primGen->GetEvent();
  125. if (event && (!event->IsSet())) {
  126. //event->SetEventID(evnr);
  127. event->SetB(b);
  128. event->SetRotZ(fPsiRP);
  129. event->MarkSet(kTRUE);
  130. }
  131. // ---> Loop over tracks in the current event
  132. for (int itrack = 0; itrack < fNpart; itrack++) {
  133. Double_t px = fPx[itrack];
  134. Double_t py = fPy[itrack];
  135. Double_t pz = fPz[itrack];
  136. /*
  137. // Lorentztransformation to lab
  138. Double_t e = fE[itrack];
  139. Double_t mass = Sqrt(e * e - px * px - py * py - pz * pz);
  140. if (gCoordinateSystem == sysLaboratory) {
  141. pz = gammaCM * (pz + betaCM * e);
  142. }
  143. Double_t ee = sqrt(mass * mass + px * px + py * py + pz * pz);
  144. */
  145. /* in Theseus 2018-03-17-bc2a06d <neutrons>=<protons> */
  146. Int_t pid= fPID[itrack];
  147. if (fPNCorr!=0. && (pid==2212 || pid==2112)) {
  148. Float_t xx=frandom->Rndm();
  149. if (xx>fPNCorr) { pid=2112; } else { pid=2212; };
  150. }
  151. /* rotate RP angle */
  152. if (fPsiRP!=0.)
  153. {
  154. Double_t cosPsi = TMath::Cos(fPsiRP);
  155. Double_t sinPsi = TMath::Sin(fPsiRP);
  156. Double_t px1=px*cosPsi-py*sinPsi;
  157. Double_t py1=px*sinPsi+py*cosPsi;
  158. px=px1; py=py1;
  159. }
  160. // Give track to PrimaryGenerator
  161. primGen->AddTrack(pid, px,py,pz, 0.0, 0.0, 0.0); // pdg, mom [GeV/c], vertex [cm]
  162. }
  163. fEventNumber++;
  164. return kTRUE;
  165. }
  166. // ------------------------------------------------------------------------
  167. ClassImp(Mpd3fdGenerator);