MpdPHSDGenerator.cxx 5.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162
  1. /**
  2. *@class MpdPHSDGenerator
  3. *@author V.Voronyuk <vadimv@jinr.ru>
  4. * MpdPHSDGenerator reads newer (instead of default) output (phsd.dat/phsd.dat.gz)
  5. * of HSD/PHSD transport model.
  6. **/
  7. #include <stdio.h>
  8. #include <stdlib.h>
  9. #include <TMath.h>
  10. #include <TParticle.h>
  11. #include <TVirtualMC.h>
  12. #include "FairMCEventHeader.h"
  13. #include "MpdMCEventHeader.h"
  14. #include "MpdPHSDGenerator.h"
  15. #include "MpdStack.h"
  16. // ---------------------------------------------------------------------
  17. MpdPHSDGenerator::MpdPHSDGenerator()
  18. {
  19. /* It is better to leave empty */
  20. };
  21. // ---------------------------------------------------------------------
  22. MpdPHSDGenerator::MpdPHSDGenerator(const char *filename="phsd.dat.gz")
  23. {
  24. fgzFile = gzopen(filename,"rb"); // zlib
  25. if (!fgzFile) {printf("-E- MpdPHSDGenerator: can not open file: %s\n",filename); exit(1);}
  26. printf("-I- MpdPHSDGenerator: open %s\n",filename);
  27. fPsiRP=0.;
  28. fisRP=kTRUE; // by default RP is random
  29. fHPol=kFALSE; // by default hyperons are not polarized
  30. frandom = new TRandom2();
  31. frandom->SetSeed(0);
  32. };
  33. // ---------------------------------------------------------------------
  34. MpdPHSDGenerator::~MpdPHSDGenerator()
  35. {
  36. if (fgzFile) {gzclose(fgzFile); fgzFile=NULL;}
  37. delete frandom;
  38. };
  39. // ---------------------------------------------------------------------
  40. Bool_t MpdPHSDGenerator::ReadEvent(FairPrimaryGenerator *primGen)
  41. {
  42. if (!ReadHeader())
  43. {
  44. printf("-I- MpdPHSDGenerator: End of file reached.\n");
  45. return kFALSE; // end of file
  46. }
  47. /* set random Reaction Plane angle */
  48. if (fisRP) {fPsiRP=frandom->Uniform(2.0*TMath::Pi());}
  49. //printf("-I- MpdPHSDGenerator: RP angle = %e\n",fPsiRP);
  50. /* Set event impact parameter in MCEventHeader if not yet done */
  51. FairMCEventHeader *eventHeader = primGen->GetEvent();
  52. if (eventHeader && (!eventHeader->IsSet()))
  53. {
  54. eventHeader->SetB(fb);
  55. eventHeader->MarkSet(kTRUE);
  56. eventHeader->SetRotZ(fPsiRP);
  57. }
  58. /* read tracks */
  59. for(Int_t i=1; i<=fntr; i++)
  60. {
  61. Int_t ipdg; Float_t px,py,pz;
  62. Float_t pol[3] = {0.,0.,0.}; // polarization
  63. /* read track */
  64. if(gzgets(fgzFile,fbuffer,256)==Z_NULL) {printf("-E- MpdPHSDGenerator: unexpected end of file %d/%d\n",i,fntr); exit(1);}
  65. /* scan values */
  66. int res=sscanf(fbuffer,"%d %*d %e %e %e",&ipdg,&px,&py,&pz);
  67. if (res!=4) {printf("-E- MpdPHSDGenerator: selftest error in track, scan %d of 4\n",res); exit(1);}
  68. if (ipdg==0) {printf("-W- MpdPHSDGenerator: particle with pdg=0\n"); continue;}
  69. // read polarizarion -- extention of standart PHSD output
  70. if (fHPol==kTRUE && (TMath::Abs(ipdg/100) == 31 || TMath::Abs(ipdg/100) == 32))
  71. {
  72. res=sscanf(fbuffer,"%*d %*d %*e %*e %*e %*e %*d %*d %*d %e %e %e",&pol[0],&pol[1],&pol[2]);
  73. if (res!=3) {printf("-E- MpdPHSDGenerator: no Hyperon polarization info %d \n",res); exit(1);}
  74. }
  75. // replacement of heavy particles for Geant4
  76. // this decays will be improved in the next versions of PHSD
  77. if (ipdg==+413) ipdg=+421; // D*+ -> D0
  78. if (ipdg==-413) ipdg=-421; // D*- -> D0_bar
  79. if (ipdg==+423) ipdg=+421; // D*0 -> D0
  80. if (ipdg==-423) ipdg=-421; // D*0_bar -> D0_bar
  81. if (ipdg==+4214) ipdg=+4122; // Sigma*_c+ -> Lambda_c+
  82. if (ipdg==-4214) ipdg=-4122; // Sigma*_c- -> Lambda_c-
  83. if (ipdg==+4114) ipdg=+4122; // Sigma*_c0 -> Lambda_c+
  84. if (ipdg==-4114) ipdg=-4122; // Sigma*_c0 -> Lambda_c-
  85. if (ipdg==+4224) ipdg=+4122; // Sigma*_c++ -> Lambda_c+
  86. if (ipdg==-4224) ipdg=-4122; // Sigma*_c-- -> Lambda_c-
  87. /* rotate RP angle */
  88. if (fPsiRP!=0.)
  89. {
  90. Double_t cosPsi = TMath::Cos(fPsiRP);
  91. Double_t sinPsi = TMath::Sin(fPsiRP);
  92. Double_t px1=px*cosPsi-py*sinPsi;
  93. Double_t py1=px*sinPsi+py*cosPsi;
  94. px=px1; py=py1;
  95. }
  96. /* add track to simulation */
  97. primGen->AddTrack(ipdg, px, py, pz, 0., 0., 0.);
  98. // add polarization for light hyperons (pdg == 31** and 32**)
  99. if (fHPol==kTRUE && (TMath::Abs(ipdg/100) == 31 || TMath::Abs(ipdg/100) == 32))
  100. {
  101. Int_t nTr = gMC->GetStack()->GetNtrack();
  102. TParticle *part = ((MpdStack*)gMC->GetStack())->GetParticle(nTr-1);
  103. if (TMath::Abs(pol[0]) > 1. || TMath::Abs(pol[1]) > 1. || TMath::Abs(pol[2]) > 1.)
  104. {
  105. part->SetPolarisation(0., 0., 0.);
  106. std::cout << "Component of the polarization vector is higher than 1!" << std::endl;
  107. }else
  108. {
  109. part->SetPolarisation(pol[0], pol[1], pol[2]);
  110. Float_t weight_pol = TMath::Sqrt(pol[0]*pol[0]+pol[1]*pol[1]+pol[2]*pol[2]);
  111. part->SetWeight(weight_pol);
  112. }
  113. }
  114. }
  115. return kTRUE;
  116. };
  117. // ---------------------------------------------------------------------
  118. Bool_t MpdPHSDGenerator::ReadHeader()
  119. {
  120. /* read header */
  121. gzgets(fgzFile,fbuffer,256); // 1st line
  122. if (gzeof(fgzFile)) return kFALSE; // end of file reached
  123. int res=sscanf(fbuffer,"%d %*d %*d %e",&fntr,&fb);
  124. if (res!=2) {printf("-E- MpdPHSDGenerator: selftest error in header, scan %d of 2\n",res); exit(1);}
  125. gzgets(fgzFile,fbuffer,256); // 2nd line
  126. if (gzeof(fgzFile)) {printf("-E- MpdPHSDGenerator: unexpected end of file\n"); exit(1);}
  127. return kTRUE;
  128. };
  129. // ---------------------------------------------------------------------
  130. void MpdPHSDGenerator::SkipTrack()
  131. {
  132. gzgets(fgzFile,fbuffer,256);
  133. if (gzeof(fgzFile)) {printf("-E- MpdPHSDGenerator: unexpected end of file\n"); exit(1);}
  134. };
  135. // ---------------------------------------------------------------------
  136. Bool_t MpdPHSDGenerator::SkipEvents(int n)
  137. {
  138. for (Int_t k=1; k<=n; ++k)
  139. {
  140. if (!ReadHeader()) return kFALSE; // end of file
  141. for (Int_t i=1; i<=fntr; ++i) SkipTrack();
  142. }
  143. return kTRUE;
  144. };
  145. // ---------------------------------------------------------------------
  146. ClassImp(MpdPHSDGenerator);