MpdEPOSGenerator.cxx 4.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128
  1. /**
  2. *@class MpdEPOSGenerator
  3. *@author K.Shtejer <kshtejer@jinr.ru>
  4. * MpdEPOSGenerator reads output of EPOS transport model in HepMC format.
  5. * (crmc_epos199.hepmc / crmc_epos199.hepmc.gz)
  6. * This interfase assumes that the two header lines and the footer line
  7. * starting with "HepMC::" have been removed.
  8. *
  9. * Last updates: June 8, 2018
  10. *
  11. **/
  12. #include <stdio.h>
  13. #include <stdlib.h>
  14. #include <TMath.h>
  15. #include "FairMCEventHeader.h"
  16. #include "MpdMCEventHeader.h"
  17. #include "MpdEPOSGenerator.h"
  18. // ---------------------------------------------------------------------
  19. MpdEPOSGenerator::MpdEPOSGenerator()
  20. {
  21. // It is better to leave empty
  22. };
  23. // ---------------------------------------------------------------------
  24. MpdEPOSGenerator::MpdEPOSGenerator(const char *filename="crmc_epos199.hepmc.gz")
  25. {
  26. fgzFile = gzopen(filename,"rb"); // zlib
  27. if (!fgzFile) {printf("-E- MpdEPOSGenerator: can not open file: %s\n",filename); exit(1);}
  28. printf("-I- MpdEPOSGenerator: open %s\n",filename);
  29. fPsiRP=0.;
  30. fisRP=kTRUE; // by default RP is random
  31. frandom = new TRandom2();
  32. frandom->SetSeed(0);
  33. };
  34. // ---------------------------------------------------------------------
  35. MpdEPOSGenerator::~MpdEPOSGenerator()
  36. {
  37. if (fgzFile) {gzclose(fgzFile); fgzFile=NULL;}
  38. delete frandom;
  39. };
  40. // ---------------------------------------------------------------------
  41. Bool_t MpdEPOSGenerator::ReadEvent(FairPrimaryGenerator *primGen)
  42. {
  43. gzgets(fgzFile,fbuffer,256); // 1st line "E": event
  44. if (gzeof(fgzFile)) return kFALSE; // end of file reached
  45. res=sscanf(fbuffer, "%s %*d %*d %*e %*e %*e %*d %*d %d %*[^\n]", charId, &numVtx);
  46. if (strcmp(charId,"E") == 0){printf(" => %s \n",charId);}
  47. if (res!=2) {printf("-E- MpdEPOSGenerator: selftest error in header, scan %d of 2\n",res); exit(1);}
  48. gzgets(fgzFile,fbuffer,256); // 2st line "U": unities
  49. if (gzeof(fgzFile)) return kFALSE;
  50. res=sscanf(fbuffer, "%s %*[^\n]", charId);
  51. if (strcmp(charId,"U") == 0){printf(" => %s \n",charId);}
  52. if (res!=1) {printf("-E- MpdEPOSGenerator: selftest error in header, scan %d of 1\n",res); exit(1);}
  53. gzgets(fgzFile,fbuffer,256); // 3st line "C": cross section
  54. if (gzeof(fgzFile)) return kFALSE;
  55. res=sscanf(fbuffer, "%s %*[^\n]", charId);
  56. if (strcmp(charId,"C") == 0){printf(" => %s \n",charId);}
  57. if (res!=1) {printf("-E- MpdEPOSGenerator: selftest error in header, scan %d of 1\n",res); exit(1);}
  58. gzgets(fgzFile,fbuffer,256); // 4st line "H": heavy ion info
  59. if (gzeof(fgzFile)) return kFALSE;
  60. res=sscanf(fbuffer, "%s %*d %*d %*d %*d %*d %*d %*d %*d %*d %e %*[^\n]", charId, &fb);
  61. if (strcmp(charId,"H") == 0){printf(" => %s \n",charId);}
  62. if (res!=2) {printf("-E- MpdEPOSGenerator: selftest error in header, scan %d of 2\n",res); exit(1);}
  63. gzgets(fgzFile,fbuffer,256); // 5st line "F"
  64. if (gzeof(fgzFile)) return kFALSE;
  65. res=sscanf(fbuffer, "%s %*d %*[^\n]", charId);
  66. if (strcmp(charId,"F") == 0){printf(" => %s \n",charId);}
  67. if (res!=1) {printf("-E- MpdEPOSGenerator: selftest error in header, scan %d of 1\n",res); exit(1);}
  68. printf("-I- Info from EPOS-> Event-> numVtx, ImpParam(fm) : %d, %e \n",numVtx, fb);
  69. /* set random Reaction Plane angle */
  70. if (fisRP) {fPsiRP=frandom->Uniform(2.0*TMath::Pi());}
  71. //printf("-I- MpdPHSDGenerator: RP angle = %e\n",fPsiRP);
  72. /* Set event impact parameter in MCEventHeader if not yet done */
  73. FairMCEventHeader *eventHeader = primGen->GetEvent();
  74. if (eventHeader && (!eventHeader->IsSet()))
  75. {
  76. eventHeader->SetB(fb);
  77. eventHeader->MarkSet(kTRUE);
  78. eventHeader->SetRotZ(fPsiRP);
  79. }
  80. // read vertexes //
  81. for(Int_t j = 1; j <= numVtx; j++)
  82. {
  83. gzgets(fgzFile,fbuffer,256); // line "V": Vertex
  84. res=sscanf(fbuffer, "%s %*d %d %*e %*e %*e %*e %d %d %*[^\n]", charId, &ivtx, &inpart, &outpart);
  85. fntr = inpart + outpart;
  86. // read particles //
  87. for(Int_t i = 1; i <= fntr; i++)
  88. {
  89. gzgets(fgzFile,fbuffer,256); // line "P": Particle
  90. res=sscanf(fbuffer, "%s %d %d %e %e %e %*e %e %d %*[^\n]", charId, &itrk, &ipdg, &px, &py, &pz, &imass, &istatus);
  91. if (istatus == 1)
  92. {
  93. // add track to simulation //
  94. primGen->AddTrack(ipdg, px, py, pz, 0., 0., 0.);
  95. printf("\n -I- Info from EPOS-> Track -> itrk, istatus, ipdg, px, py, pz : %d %d %d %e %e %e", itrk, istatus, ipdg, px, py, pz);
  96. } else {
  97. printf("\n -I- Info from EPOS-> Skipped Track -> itrk, istatus : %d %d", itrk, istatus);
  98. }
  99. }
  100. if (gzeof(fgzFile)) return kFALSE;
  101. }
  102. return kTRUE;
  103. };
  104. // ---------------------------------------------------------------------
  105. ClassImp(MpdEPOSGenerator);