MpdHistoGenerator.cxx 3.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124
  1. /**
  2. * MpdHistoGenerator.cxx
  3. * *@author A.Zinchenko <Alexander.Zinchenko@jinr.ru>
  4. * The MpdHistoGenerator produces particles according to their Eta-Pt histogram
  5. **/
  6. #include "MpdHistoGenerator.h"
  7. #include "FairMCEventHeader.h"
  8. #include "FairPrimaryGenerator.h"
  9. #include "TDatabasePDG.h"
  10. #include "TFile.h"
  11. #include "TH2D.h"
  12. #include "TLorentzVector.h"
  13. #include "TMath.h"
  14. #include "TRandom.h"
  15. #include "TVirtualMC.h"
  16. #include <fstream>
  17. #include <iostream>
  18. #include <iomanip>
  19. // #define debug_hsd
  20. using namespace std;
  21. // ------------------------------------------------------------------------
  22. MpdHistoGenerator::MpdHistoGenerator() :
  23. FairGenerator(), fPdgCode(-211), fMult(1), fHist(NULL), fYield(-1.0),
  24. fFileName("")
  25. {
  26. // Constructor
  27. //gRandom->SetSeed(0); // change random seed
  28. // Get histograms Pt vs Eta (He3)
  29. //TFile f("he3_eta_pt.root");
  30. //f.ReadAll();
  31. //fHist = (TH2D*) f.FindObject("hPtEta");
  32. //fHist->SetDirectory(0);
  33. }
  34. // ------------------------------------------------------------------------
  35. MpdHistoGenerator::MpdHistoGenerator(Int_t pdgid, Int_t mult, Double_t yield) :
  36. FairGenerator(), fPdgCode(pdgid), fMult(mult), fHist(NULL), fYield(yield),
  37. fFileName("")
  38. {
  39. // Constructor
  40. //gRandom->SetSeed(0); // change random seed
  41. // Get histograms Pt vs Eta (He3)
  42. //TFile f("he3_eta_pt.root");
  43. //f.ReadAll();
  44. //fHist = (TH2D*) f.FindObject("hPtVsEta");
  45. //fHist->SetDirectory(0);
  46. }
  47. // ------------------------------------------------------------------------
  48. Bool_t MpdHistoGenerator::Init()
  49. {
  50. // Get histograms Pt vs Eta (He3)
  51. if (fFileName == "") fFileName = "genPtVsEta.root";
  52. TFile *f = TFile::Open(fFileName);
  53. f->ReadAll();
  54. fHist = (TH2D*) f->FindObject("hPtVsEta");
  55. fHist->SetDirectory(0);
  56. return kTRUE;
  57. }
  58. // ------------------------------------------------------------------------
  59. MpdHistoGenerator::~MpdHistoGenerator()
  60. {
  61. // Destructor
  62. }
  63. //-------------------------------------------------------------------------
  64. Bool_t MpdHistoGenerator::ReadEvent(FairPrimaryGenerator* primGen)
  65. {
  66. // Generate event
  67. static Int_t first = 1;
  68. if (first) {
  69. first = 0;
  70. //cout << " histo " << fHist << " " << fHist->GetNbinsX() << " " << fHist->GetNbinsY() << endl;
  71. if (TDatabasePDG::Instance()->GetParticle(fPdgCode) == NULL) {
  72. // Define particle with fPdgCode (take it from DCM-QGSM input file)
  73. // To be done ...
  74. }
  75. }
  76. if (fYield > 0 && gRandom->Rndm() > fYield) return kTRUE;
  77. Double_t pt, eta, phi;
  78. TVector3 part;
  79. for (Int_t i = 0; i < fMult; ++i) {
  80. phi = gRandom->Uniform(TMath::TwoPi());
  81. fHist->GetRandom2(eta, pt);
  82. part.SetPtEtaPhi(pt, eta, phi);
  83. //FairMCEventHeader* eventHeader = primGen->GetEvent();
  84. /* Set event impact parameter in MCEvent if not yet done */
  85. /*
  86. if (eventHeader && (!eventHeader->IsSet()) ) {
  87. eventHeader->SetB(fB);
  88. eventHeader->MarkSet(kTRUE);
  89. }
  90. */
  91. //primGen->AddTrack(1000020030, He3.Px(), He3.Py(), He3.Pz(), 0., 0., 0., -1);
  92. primGen->AddTrack(fPdgCode, part.Px(), part.Py(), part.Pz(), 0., 0., 0., -1);
  93. }
  94. return kTRUE;
  95. }
  96. // ------------------------------------------------------------------------