MpdPHQMDGenerator.cxx 12 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277
  1. /**
  2. *@class MpdPHQMDGenerator
  3. *@author V.Kireyeu, V. Voronyuk
  4. **/
  5. #include <stdio.h>
  6. #include <stdlib.h>
  7. #include <TMath.h>
  8. #include <TMCParticleType.h>
  9. #include <TDatabasePDG.h>
  10. #include <TParticlePDG.h>
  11. #include <TLorentzVector.h>
  12. #include "FairMCEventHeader.h"
  13. #include "MpdMCEventHeader.h"
  14. #include "MpdPHQMDGenerator.h"
  15. // ---------------------------------------------------------------------
  16. MpdPHQMDGenerator::MpdPHQMDGenerator()
  17. {
  18. /* It is better to leave empty */
  19. };
  20. // ---------------------------------------------------------------------
  21. MpdPHQMDGenerator::MpdPHQMDGenerator(const char *filename="phsd.dat.gz", const char *f79name="fort.79.gz")
  22. {
  23. fgzFile = gzopen(filename,"rb"); // zlib
  24. fgzFile79 = gzopen(f79name,"rb"); // zlib
  25. if (!fgzFile) {printf("-E- MpdPHQMDGenerator: can not open file: %s\n",filename); exit(1);}
  26. printf("-I- MpdPHQMDGenerator: open %s\n",filename);
  27. if (!fgzFile79) {printf("-E- MpdPHQMDGenerator: can not open file: %s\n",f79name); exit(1);}
  28. printf("-I- MpdPHQMDGenerator: open %s\n",f79name);
  29. fPsiRP=0.;
  30. fisRP=kTRUE; // by default RP is random
  31. frandom = new TRandom2();
  32. frandom->SetSeed(0);
  33. TDatabasePDG *fPDG = TDatabasePDG::Instance();
  34. TParticlePDG* particle = 0;
  35. kProtonMass = fPDG -> GetParticle(2212) -> Mass();
  36. kNeutronMass = fPDG -> GetParticle(2112) -> Mass();
  37. kLambdaMass = fPDG -> GetParticle(3122) -> Mass();
  38. // Light nuclei
  39. particle = fPDG->AddParticle( "d", "d", 1.8761199, 0, 1.6916047e-33, 1*3, "nucleus", 1000010020, 0, kPTHadron);
  40. particle = fPDG->AddParticle( "t", "t", 2.8094298, 0, 1.6916047e-33, 1*3, "nucleus", 1000010030, 0, kPTHadron);
  41. particle = fPDG->AddParticle( "He3", "He3", 2.8094101, 0, 1.6916047e-33, 2*3, "nucleus", 1000020030, 0, kPTHadron);
  42. particle = fPDG->AddParticle( "He4", "He4", 3.7283999, 0, 1.6916047e-33, 2*3, "nucleus", 1000020040, 0, kPTHadron);
  43. particle = fPDG->AddParticle( "He5", "He5", 5.01222363, 0, 1.6916047e-33, 2*3, "nucleus", 1000020050, 0, kPTHadron);
  44. particle = fPDG->AddParticle( "He6", "He6", 5.6056, 0, 1.6916047e-33, 2*3, "nucleus", 1000020060, 0, kPTHadron);
  45. particle = fPDG->AddParticle( "Li6", "Li6", 5.6016, 0, 1.6916047e-33, 3*3, "nucleus", 1000030060, 0, kPTHadron);
  46. particle = fPDG->AddParticle( "Li7", "Li7", 6.5339, 0, 1.6916047e-33, 3*3, "nucleus", 1000030070, 0, kPTHadron);
  47. particle = fPDG->AddParticle( "Li8", "Li8", 7.4714, 0, 1.6916047e-33, 3*3, "nucleus", 1000030080, 0, kPTHadron);
  48. particle = fPDG->AddParticle( "Be7", "Be7", 6.5342, 0, 1.6916047e-33, 4*3, "nucleus", 1000040070, 0, kPTHadron);
  49. particle = fPDG->AddParticle( "Be8", "Be8", 7.4569, 0, 1.6916047e-33, 4*3, "nucleus", 1000040080, 0, kPTHadron);
  50. particle = fPDG->AddParticle( "Be9", "Be9", 8.3959, 0, 1.6916047e-33, 4*3, "nucleus", 1000040090, 0, kPTHadron);
  51. particle = fPDG->AddParticle( "B8", "B8", 7.4724, 0, 1.6916047e-33, 5*3, "nucleus", 1000050080, 0, kPTHadron);
  52. particle = fPDG->AddParticle( "B9", "B9", 8.3959, 0, 1.6916047e-33, 5*3, "nucleus", 1000050090, 0, kPTHadron);
  53. // Hypernuclei desribed in UserDecay.C
  54. };
  55. // ---------------------------------------------------------------------
  56. MpdPHQMDGenerator::~MpdPHQMDGenerator()
  57. {
  58. if (fgzFile) {gzclose(fgzFile); fgzFile=NULL;}
  59. if (fgzFile79) {gzclose(fgzFile79); fgzFile79=NULL;}
  60. delete frandom;
  61. };
  62. // ---------------------------------------------------------------------
  63. Bool_t MpdPHQMDGenerator::ReadEvent(FairPrimaryGenerator *primGen)
  64. {
  65. if (!ReadHeader())
  66. {
  67. printf("-I- MpdPHQMDGenerator: End of file reached.\n");
  68. return kFALSE; // end of file
  69. }
  70. /* set random Reaction Plane angle */
  71. if (fisRP) {fPsiRP=frandom->Uniform(2.0*TMath::Pi());}
  72. //printf("-I- MpdPHQMDGenerator: RP angle = %e\n",fPsiRP);
  73. /* Set event impact parameter in MCEventHeader if not yet done */
  74. FairMCEventHeader *eventHeader = primGen->GetEvent();
  75. if (eventHeader && (!eventHeader->IsSet()))
  76. {
  77. eventHeader->SetB(fb);
  78. eventHeader->MarkSet(kTRUE);
  79. eventHeader->SetRotZ(fPsiRP);
  80. }
  81. Int_t ipdg;
  82. Float_t px,py,pz;
  83. Float_t mass, amass;
  84. Int_t nZ, nN, nL, nS;
  85. /* read tracks: fort.79 - baryons + clusters */
  86. for(Int_t i=1; i<=fntr79; i++)
  87. {
  88. if(gzgets(fgzFile79,fbuffer79,256)==Z_NULL) {printf("-E- MpdPHQMDGenerator: unexpected end of baryons file %d/%d\n",i,fntr79); exit(1);}
  89. int res79=sscanf(fbuffer79,"%*d %d %d %e %e %e %*e %*e %*e %*e %e %e %d %d", &nZ, &nN, &px, &py, &pz, &mass, &amass, &nL, &nS);
  90. if (res79!=9) {printf("-E- MpdPHQMDGenerator: selftest error in track, scan %d of 9\n",res79); exit(1);}
  91. ipdg = BaryonPDG(nN, nZ, nL, nS, amass);
  92. Float_t pAbs = 0.0, tmpZ, st, phi;
  93. TVector3 velocity;
  94. TLorentzVector pMother, tmpMom;
  95. if(nZ > 1 && nN == 0 && nL == 0 && nS == 0){ // bounded protons case
  96. pAbs = 0.5 * sqrt(4 * pow(kProtonMass, 2*nZ)) / amass;
  97. pMother.SetPxPyPzE(px, py, pz, sqrt(pow(px,2) + pow(py,2) + pow(pz,2) +pow(amass,2)));
  98. for(int bound = 0; bound < nZ; ++bound){
  99. tmpZ = 1. - 2.* (gRandom->Rndm());
  100. st = sqrt((1.-tmpZ) * (1.+tmpZ)) * pAbs;
  101. phi = 2. * TMath::Pi() * (gRandom->Rndm());
  102. tmpMom.SetXYZM(st*cos(phi), st*sin(phi), tmpZ*pAbs, kProtonMass);
  103. tmpMom.Boost(velocity);
  104. primGen->AddTrack(2212, tmpMom.Px(), tmpMom.Py(), tmpMom.Pz(), 0., 0., 0.);
  105. }
  106. }
  107. else if (nZ == 0 && nN > 1 && nL == 0 && nS == 0){ // bounded neutrons case
  108. pAbs = 0.5 * sqrt(4 * pow(kNeutronMass, 2*nN)) / amass;
  109. pMother.SetPxPyPzE(px, py, pz, sqrt(pow(px,2) + pow(py,2) + pow(pz,2) +pow(amass,2)));
  110. for(int bound = 0; bound < nN; ++bound){
  111. tmpZ = 1. - 2.* (gRandom->Rndm());
  112. st = sqrt((1.-tmpZ) * (1.+tmpZ)) * pAbs;
  113. phi = 2. * TMath::Pi() * (gRandom->Rndm());
  114. tmpMom.SetXYZM(st*cos(phi), st*sin(phi), tmpZ*pAbs, kNeutronMass);
  115. tmpMom.Boost(velocity);
  116. primGen->AddTrack(2112, tmpMom.Px(), tmpMom.Py(), tmpMom.Pz(), 0., 0., 0.);
  117. }
  118. }
  119. else if (nZ == 0 && nN > 1 && nL >= 1 && nS == 0){ // bounded lambdas + neutrons case
  120. pAbs = 0.5 * sqrt(4 * pow(kNeutronMass, 2*nN) * pow(kLambdaMass, 2*nL)) / amass;
  121. pMother.SetPxPyPzE(px, py, pz, sqrt(pow(px,2) + pow(py,2) + pow(pz,2) +pow(amass,2)));
  122. for(int bound = 0; bound < nN; ++bound){ // free neutrons
  123. tmpZ = 1. - 2.* (gRandom->Rndm());
  124. st = sqrt((1.-tmpZ) * (1.+tmpZ)) * pAbs;
  125. phi = 2. * TMath::Pi() * (gRandom->Rndm());
  126. tmpMom.SetXYZM(st*cos(phi), st*sin(phi), tmpZ*pAbs, kNeutronMass);
  127. tmpMom.Boost(velocity);
  128. primGen->AddTrack(2112, tmpMom.Px(), tmpMom.Py(), tmpMom.Pz(), 0., 0., 0.);
  129. }
  130. for(int bound = 0; bound < nL; ++bound){ // free lambdas
  131. tmpZ = 1. - 2.* (gRandom->Rndm());
  132. st = sqrt((1.-tmpZ) * (1.+tmpZ)) * pAbs;
  133. phi = 2. * TMath::Pi() * (gRandom->Rndm());
  134. tmpMom.SetXYZM(st*cos(phi), st*sin(phi), tmpZ*pAbs, kLambdaMass);
  135. tmpMom.Boost(velocity);
  136. primGen->AddTrack(3122, tmpMom.Px(), tmpMom.Py(), tmpMom.Pz(), 0., 0., 0.);
  137. }
  138. }
  139. else{ // normal case
  140. primGen->AddTrack(ipdg, px, py, pz, 0., 0., 0.);
  141. }
  142. }
  143. /* read tracks: phsd.dat */
  144. for(Int_t i=1; i<=fntr; i++)
  145. {
  146. /* read track */
  147. if(gzgets(fgzFile,fbuffer,256)==Z_NULL) {printf("-E- MpdPHQMDGenerator: unexpected end of file %d/%d\n",i,fntr); exit(1);}
  148. /* scan values */
  149. int res=sscanf(fbuffer,"%d %*d %e %e %e",&ipdg,&px,&py,&pz);
  150. if (res!=4) {printf("-E- MpdPHQMDGenerator: selftest error in track, scan %d of 4\n",res); exit(1);}
  151. if (ipdg==0) {printf("-W- MpdPHQMDGenerator: particle with pdg=0\n"); continue;}
  152. if (ipdg > 999) continue; // only mesons should be readed
  153. // replacement of heavy particles for Geant4
  154. // this decays will be improved in the next versions of PHQMD
  155. if (ipdg==+413) ipdg=+421; // D*+ -> D0
  156. if (ipdg==-413) ipdg=-421; // D*- -> D0_bar
  157. if (ipdg==+423) ipdg=+421; // D*0 -> D0
  158. if (ipdg==-423) ipdg=-421; // D*0_bar -> D0_bar
  159. if (ipdg==+4214) ipdg=+4122; // Sigma*_c+ -> Lambda_c+
  160. if (ipdg==-4214) ipdg=-4122; // Sigma*_c- -> Lambda_c-
  161. if (ipdg==+4114) ipdg=+4122; // Sigma*_c0 -> Lambda_c+
  162. if (ipdg==-4114) ipdg=-4122; // Sigma*_c0 -> Lambda_c-
  163. if (ipdg==+4224) ipdg=+4122; // Sigma*_c++ -> Lambda_c+
  164. if (ipdg==-4224) ipdg=-4122; // Sigma*_c-- -> Lambda_c-
  165. /* rotate RP angle */
  166. if (fPsiRP!=0.)
  167. {
  168. Double_t cosPsi = TMath::Cos(fPsiRP);
  169. Double_t sinPsi = TMath::Sin(fPsiRP);
  170. Double_t px1=px*cosPsi-py*sinPsi;
  171. Double_t py1=px*sinPsi+py*cosPsi;
  172. px=px1; py=py1;
  173. }
  174. /* add track to simulation */
  175. primGen->AddTrack(ipdg, px, py, pz, 0., 0., 0.);
  176. }
  177. return kTRUE;
  178. };
  179. // ---------------------------------------------------------------------
  180. // ---------------------------------------------------------------------
  181. Int_t MpdPHQMDGenerator::BaryonPDG(Int_t N, Int_t Z, Int_t L, Int_t S, Float_t mass)
  182. {
  183. Int_t I = 0;
  184. Int_t A = N + Z + L + S;
  185. Int_t LS = L + S;
  186. Int_t PDG = 999; // some default value
  187. if (N == 1 && Z == 0 && L == 0 && S == 0) PDG = 2112; // neutron
  188. else if (N == 0 && Z == 1 && L == 0 && S == 0) PDG = 2212; // proton
  189. else if (N == 0 && Z == 0 && L == 1 && S == 0) PDG = 3122; // lambda0
  190. else if (N == 0 && Z == 1 && L == 0 && S == 1 && abs(mass - kSigmaPMass) < 0.0001) PDG = 3222; // sigma+
  191. else if (N == 0 && Z == -1 && L == 0 && S == 1 && abs(mass - kSigmaPMass) < 0.0001) PDG = -3222; // sigma+
  192. else if (N == 0 && Z == 1 && L == 0 && S == 1 && abs(mass - kSigmaNMass) < 0.0001) PDG = 3212; // sigma0
  193. else if (N == 0 && Z == -1 && L == 0 && S == 1 && abs(mass - kSigmaNMass) < 0.0001) PDG = -3212; // sigma0
  194. else if (N == 0 && Z == 1 && L == 0 && S == 1 && abs(mass - kSigmaMMass) < 0.0001) PDG = 3112; // sigma-
  195. else if (N == 0 && Z == -1 && L == 0 && S == 1 && abs(mass - kSigmaMMass) < 0.0001) PDG = 3112; // sigma-
  196. else PDG = 10*10E7 + LS*10E6 + Z*10E3 + A*10 + I;
  197. return PDG;
  198. };
  199. // ---------------------------------------------------------------------
  200. Bool_t MpdPHQMDGenerator::ReadHeader()
  201. {
  202. /* read header */
  203. gzgets(fgzFile,fbuffer,256); // 1st line
  204. if (gzeof(fgzFile)) return kFALSE; // end of file reached
  205. int res=sscanf(fbuffer,"%d %*d %*d %e",&fntr,&fb);
  206. if (res!=2) {printf("-E- MpdPHQMDGenerator: selftest error in header, scan %d of 2\n",res); exit(1);}
  207. gzgets(fgzFile,fbuffer,256); // 2nd line
  208. if (gzeof(fgzFile)) {printf("-E- MpdPHQMDGenerator: unexpected end of file\n"); exit(1);}
  209. gzgets(fgzFile79,fbuffer79,256); // 1st line - unused
  210. gzgets(fgzFile79,fbuffer79,256); // 2nd line - unused
  211. gzgets(fgzFile79,fbuffer79,256); // event header
  212. if (gzeof(fgzFile79)) return kFALSE; // end of file reached
  213. int res79=sscanf(fbuffer79,"%d",&fntr79);
  214. if (res79!=1) {printf("-E- MpdPHQMDGenerator: selftest error in header, scan %d of 1\n",res79); exit(1);}
  215. if (gzeof(fgzFile79)) {printf("-E- MpdPHQMDGenerator: unexpected end of baryons file\n"); exit(1);}
  216. return kTRUE;
  217. };
  218. // ---------------------------------------------------------------------
  219. void MpdPHQMDGenerator::SkipTrack()
  220. {
  221. gzgets(fgzFile,fbuffer,256);
  222. if (gzeof(fgzFile)) {printf("-E- MpdPHQMDGenerator: unexpected end of file\n"); exit(1);}
  223. };
  224. // ---------------------------------------------------------------------
  225. Bool_t MpdPHQMDGenerator::SkipEvents(int n)
  226. {
  227. for (Int_t k=1; k<=n; ++k)
  228. {
  229. if (!ReadHeader()) return kFALSE; // end of file
  230. for (Int_t i=1; i<=fntr; ++i) SkipTrack();
  231. }
  232. return kTRUE;
  233. };
  234. // ---------------------------------------------------------------------
  235. ClassImp(MpdPHQMDGenerator);