MpdUrqmdGenerator.cxx 13 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432
  1. // -------------------------------------------------------------------------
  2. // ----- MpdUrqmdGenerator source file -----
  3. // ----- Created 24/06/04 by V. Friese -----
  4. // -------------------------------------------------------------------------
  5. #include "MpdUrqmdGenerator.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 <iostream>
  20. #include <cstring>
  21. #include <stdio.h>
  22. using namespace std;
  23. using namespace TMath;
  24. //const Double_t kProtonMass = 0.938271998;
  25. // ----- Default constructor ------------------------------------------
  26. MpdUrqmdGenerator::MpdUrqmdGenerator()
  27. : FairGenerator(),
  28. fInputFile(NULL),
  29. fParticleTable(),
  30. fFileName(NULL),
  31. fPhiMin(0.),
  32. fPhiMax(0.),
  33. fEventPlaneSet(kFALSE) {
  34. }
  35. // ------------------------------------------------------------------------
  36. // ----- Standard constructor -----------------------------------------
  37. MpdUrqmdGenerator::MpdUrqmdGenerator(const char* fileName)
  38. : FairGenerator(),
  39. fInputFile(NULL),
  40. fParticleTable(),
  41. fFileName(fileName),
  42. fPhiMin(0.),
  43. fPhiMax(0.),
  44. fEventPlaneSet(kFALSE) {
  45. // fFileName = fileName;
  46. cout << "-I MpdUrqmdGenerator: Opening input file " << fFileName << endl;
  47. fInputFile = gzopen(fFileName, "rb");
  48. if (!fInputFile) {
  49. Fatal("MpdUrqmdGenerator", "Cannot open input file.");
  50. exit(1);
  51. }
  52. ReadConversionTable();
  53. }
  54. // ------------------------------------------------------------------------
  55. // ----- Destructor ---------------------------------------------------
  56. MpdUrqmdGenerator::~MpdUrqmdGenerator() {
  57. // cout<<"Enter Destructor of MpdUrqmdGenerator"<<endl;
  58. if (fInputFile) {
  59. #ifdef GZIP_SUPPORT
  60. gzclose(fInputFile);
  61. #else
  62. fclose(fInputFile);
  63. #endif
  64. fInputFile = NULL;
  65. }
  66. fParticleTable.clear();
  67. // cout<<"Leave Destructor of MpdUrqmdGenerator"<<endl;
  68. }
  69. // ------------------------------------------------------------------------
  70. // ----- Public method ReadEvent --------------------------------------
  71. Bool_t MpdUrqmdGenerator::ReadEvent(FairPrimaryGenerator* primGen) {
  72. // ---> Check for input file
  73. if (!fInputFile) {
  74. cout << "-E MpdUrqmdGenerator: Input file not open! " << endl;
  75. return kFALSE;
  76. }
  77. // ---> Check for primary generator
  78. if (!primGen) {
  79. cout << "-E- MpdUrqmdGenerator::ReadEvent: "
  80. << "No PrimaryGenerator!" << endl;
  81. return kFALSE;
  82. }
  83. // ---> Define event variables to be read from file
  84. int evnr = 0, ntracks = 0, aProj = 0, zProj = 0, aTarg = 0, zTarg = 0;
  85. float b = 0., ekin = 0.;
  86. int ityp = 0, i3 = 0, ichg = 0, pid = 0;
  87. float ppx = 0., ppy = 0., ppz = 0., m = 0.;
  88. // ---> Read and check first event header line from input file
  89. char read[200];
  90. #ifdef GZIP_SUPPORT
  91. gzgets(fInputFile, read, 200); // line 1
  92. #else
  93. fgets(read, 200, fInputFile);
  94. #endif
  95. Int_t urqmdVersion = 0;
  96. sscanf(read, "UQMD version: %d 1000 %d output_file 14", &urqmdVersion, &urqmdVersion);
  97. cout << "URQMD VERSION USED = " << urqmdVersion << endl;
  98. #ifdef GZIP_SUPPORT
  99. if (gzeof(fInputFile)) {
  100. cout << "-I MpdUrqmdGenerator : End of input file reached." << endl;
  101. gzclose(fInputFile);
  102. #else
  103. if ( feof(fInputFile) ) {
  104. cout << "-I MpdUrqmdGenerator : End of input file reached." << endl;
  105. fclose(fInputFile);
  106. #endif
  107. fInputFile = NULL;
  108. return kFALSE;
  109. }
  110. if (read[0] != 'U') {
  111. cout << "-E MpdUrqmdGenerator: Wrong event header" << endl;
  112. return kFALSE;
  113. }
  114. // ---> Read rest of event header
  115. #ifdef GZIP_SUPPORT
  116. gzgets(fInputFile, read, 200); // line 2
  117. sscanf(read, "projectile: (mass, char) %d %d target: (mass, char) %d %d",
  118. &aProj, &zProj, &aTarg, &zTarg); // line 2
  119. gzgets(fInputFile, read, 200); // line 3
  120. gzgets(fInputFile, read, 36); // line 4
  121. gzgets(fInputFile, read, 200); // line 4
  122. sscanf(read, "%f", &b); // line 4
  123. gzgets(fInputFile, read, 39); // line 5
  124. gzgets(fInputFile, read, 200); // line 5
  125. sscanf(read, "%e", &ekin); // line 5
  126. gzgets(fInputFile, read, 7); // line 6
  127. gzgets(fInputFile, read, 200); // line 6
  128. sscanf(read, "%d", &evnr); // line 6
  129. for (Int_t iline = 0; iline < ((urqmdVersion == 30400) ? 11 : 8); iline++)
  130. gzgets(fInputFile, read, 200);
  131. gzgets(fInputFile, read, 200); // line 18
  132. sscanf(read, "%d", &ntracks); // line 18
  133. gzgets(fInputFile, read, 200); // line 19
  134. #else
  135. fgets(read, 26, fInputFile);
  136. fscanf(fInputFile, "%d", &aProj);
  137. fscanf(fInputFile, "%d", &zProj);
  138. fgets(read, 25, fInputFile);
  139. fscanf(fInputFile, "%d", &aTarg);
  140. fscanf(fInputFile, "%d", &zTarg);
  141. fgets(read, 200, fInputFile);
  142. fgets(read, 200, fInputFile);
  143. fgets(read, 36, fInputFile);
  144. fscanf(fInputFile, "%f", &b);
  145. fgets(read, 200, fInputFile);
  146. fgets(read, 39, fInputFile);
  147. fscanf(fInputFile, "%e", &ekin);
  148. fgets(read, 200, fInputFile);
  149. fgets(read, 7, fInputFile);
  150. fscanf(fInputFile, "%d", &evnr);
  151. fgets(read, 200, fInputFile);
  152. for (int iline=0; iline<8; iline++) { fgets(read, 200,fInputFile); }
  153. fscanf(fInputFile, "%d", &ntracks);
  154. fgets(read, 200, fInputFile);
  155. fgets(read, 200, fInputFile);
  156. #endif
  157. // ---> Calculate beta and gamma for Lorentztransformation
  158. TDatabasePDG* pdgDB = TDatabasePDG::Instance();
  159. TParticlePDG* kProton = pdgDB->GetParticle(2212);
  160. Double_t kProtonMass = kProton->Mass();
  161. Double_t eBeam = ekin + kProtonMass;
  162. Double_t pBeam = TMath::Sqrt(eBeam * eBeam - kProtonMass * kProtonMass);
  163. Double_t betaCM = pBeam / (eBeam + kProtonMass);
  164. Double_t gammaCM = TMath::Sqrt(1. / (1. - betaCM * betaCM));
  165. cout << "-I MpdUrqmdGenerator: Event " << evnr << ", b = " << b
  166. << " fm, multiplicity " << ntracks << ", ekin: " << ekin << endl;
  167. Double_t phi = 0.;
  168. // ---> Generate rotation angle
  169. if (fEventPlaneSet) {
  170. phi = gRandom->Uniform(fPhiMin, fPhiMax);
  171. }
  172. // Set event id and impact parameter in MCEvent if not yet done
  173. FairMCEventHeader* event = primGen->GetEvent();
  174. if (event && (!event->IsSet())) {
  175. event->SetEventID(evnr);
  176. event->SetB(b);
  177. event->MarkSet(kTRUE);
  178. event->SetRotZ(phi);
  179. }
  180. // ---> Loop over tracks in the current event
  181. for (int itrack = 0; itrack < ntracks; itrack++) {
  182. // Read momentum and PID from file
  183. #ifdef GZIP_SUPPORT
  184. gzgets(fInputFile, read, 81);
  185. gzgets(fInputFile, read, 200);
  186. sscanf(read, "%e %e %e %e %d %d %d", &ppx, &ppy, &ppz, &m, &ityp, &i3, &ichg);
  187. #else
  188. fgets(read, 81, fInputFile);
  189. fscanf(fInputFile, "%e", &ppx);
  190. fscanf(fInputFile, "%e", &ppy);
  191. fscanf(fInputFile, "%e", &ppz);
  192. fscanf(fInputFile, "%e", &m);
  193. fscanf(fInputFile, "%d", &ityp);
  194. fscanf(fInputFile, "%d", &i3);
  195. fscanf(fInputFile, "%d", &ichg);
  196. fgets(read, 200, fInputFile);
  197. #endif
  198. // Convert UrQMD type and charge to unique pid identifier
  199. if (ityp >= 0) {
  200. pid = 1000 * (ichg + 2) + ityp;
  201. } else {
  202. pid = -1000 * (ichg + 2) + ityp;
  203. }
  204. // Convert Unique PID into PDG particle code
  205. if (fParticleTable.find(pid) == fParticleTable.end()) {
  206. cout << "-W MpdUrqmdGenerator: PID " << ityp << " charge "
  207. << ichg << " not found in table (" << pid << ")" << endl;
  208. continue;
  209. }
  210. Int_t pdgID = fParticleTable[pid];
  211. // Lorentztransformation to lab
  212. Double_t mass = Double_t(m);
  213. Double_t px = Double_t(ppx);
  214. Double_t py = Double_t(ppy);
  215. Double_t pz = Double_t(ppz);
  216. Double_t e = sqrt(mass * mass + px * px + py * py + pz * pz);
  217. if (gCoordinateSystem == sysLaboratory)
  218. pz = gammaCM * (pz + betaCM * e);
  219. Double_t ee = sqrt(mass * mass + px * px + py * py + pz * pz);
  220. if (fEventPlaneSet) {
  221. Double_t pt = Sqrt(px * px + py * py);
  222. Double_t azim = ATan2(py, px);
  223. azim += phi;
  224. px = pt * Cos(azim);
  225. py = pt * Sin(azim);
  226. }
  227. TLorentzVector pp;
  228. pp.SetPx(px);
  229. pp.SetPy(py);
  230. pp.SetPz(pz);
  231. pp.SetE(ee);
  232. // Give track to PrimaryGenerator
  233. primGen->AddTrack(pdgID, px, py, pz, 0., 0., 0.);
  234. }
  235. return kTRUE;
  236. }
  237. // ------------------------------------------------------------------------
  238. // ----- Public method ReadEvent --------------------------------------
  239. Bool_t MpdUrqmdGenerator::SkipEvents(Int_t count) {
  240. if (count <= 0) {
  241. return kTRUE;
  242. }
  243. for (Int_t ii = 0; ii < count; ii++) {
  244. // ---> Check for input file
  245. if (!fInputFile) {
  246. cout << "-E MpdUrqmdGenerator: Input file not open! " << endl;
  247. return kFALSE;
  248. }
  249. // ---> Define event variables to be read from file
  250. int evnr = 0, ntracks = 0, aProj = 0, zProj = 0, aTarg = 0, zTarg = 0;
  251. float b = 0., ekin = 0.;
  252. // ---> Read and check first event header line from input file
  253. char read[200];
  254. #ifdef GZIP_SUPPORT
  255. gzgets(fInputFile, read, 200); // line 1
  256. #else
  257. fgets(read, 200, fInputFile);
  258. #endif
  259. Int_t urqmdVersion = 0;
  260. sscanf(read, "UQMD version: %d 1000 %d output_file 14", &urqmdVersion, &urqmdVersion);
  261. cout << "URQMD VERSION USED = " << urqmdVersion << endl;
  262. #ifdef GZIP_SUPPORT
  263. if (gzeof(fInputFile)) {
  264. cout << "-I MpdUrqmdGenerator : End of input file reached." << endl;
  265. gzclose(fInputFile);
  266. #else
  267. if ( feof(fInputFile) ) {
  268. cout << "-I MpdUrqmdGenerator : End of input file reached." << endl;
  269. fclose(fInputFile);
  270. #endif
  271. fInputFile = NULL;
  272. return kFALSE;
  273. }
  274. if (read[0] != 'U') {
  275. cout << "-E MpdUrqmdGenerator: Wrong event header" << endl;
  276. return kFALSE;
  277. }
  278. // ---> Read rest of event header
  279. #ifdef GZIP_SUPPORT
  280. gzgets(fInputFile, read, 200); // line 2
  281. sscanf(read, "projectile: (mass, char) %d %d target: (mass, char) %d %d",
  282. &aProj, &zProj, &aTarg, &zTarg); // line 2
  283. gzgets(fInputFile, read, 200); // line 3
  284. gzgets(fInputFile, read, 36); // line 4
  285. gzgets(fInputFile, read, 200); // line 4
  286. sscanf(read, "%f", &b); // line 4
  287. gzgets(fInputFile, read, 39); // line 5
  288. gzgets(fInputFile, read, 200); // line 5
  289. sscanf(read, "%e", &ekin); // line 5
  290. gzgets(fInputFile, read, 7); // line 6
  291. gzgets(fInputFile, read, 200); // line 6
  292. sscanf(read, "%d", &evnr); // line 6
  293. for (Int_t iline = 0; iline < ((urqmdVersion == 30400) ? 11 : 8); iline++)
  294. gzgets(fInputFile, read, 200);
  295. gzgets(fInputFile, read, 200); // line 15
  296. sscanf(read, "%d", &ntracks); // line 15
  297. gzgets(fInputFile, read, 200); // line 16
  298. #else
  299. fgets(read, 26, fInputFile);
  300. fscanf(fInputFile, "%d", &aProj);
  301. fscanf(fInputFile, "%d", &zProj);
  302. fgets(read, 25, fInputFile);
  303. fscanf(fInputFile, "%d", &aTarg);
  304. fscanf(fInputFile, "%d", &zTarg);
  305. fgets(read, 200, fInputFile);
  306. fgets(read, 200, fInputFile);
  307. fgets(read, 36, fInputFile);
  308. fscanf(fInputFile, "%f", &b);
  309. fgets(read, 200, fInputFile);
  310. fgets(read, 39, fInputFile);
  311. fscanf(fInputFile, "%e", &ekin);
  312. fgets(read, 200, fInputFile);
  313. fgets(read, 7, fInputFile);
  314. fscanf(fInputFile, "%d", &evnr);
  315. fgets(read, 200, fInputFile);
  316. for (int iline=0; iline<8; iline++) { fgets(read, 200,fInputFile); }
  317. fscanf(fInputFile, "%d", &ntracks);
  318. fgets(read, 200, fInputFile);
  319. fgets(read, 200, fInputFile);
  320. #endif
  321. cout << "-I MpdUrqmdGenerator: Event " << evnr << " skipped!" << endl;
  322. // ---> Loop over tracks in the current event
  323. for (int itrack = 0; itrack < ntracks; itrack++) {
  324. // Read momentum and PID from file
  325. #ifdef GZIP_SUPPORT
  326. gzgets(fInputFile, read, 200);
  327. #else
  328. fgets(read, 81, fInputFile);
  329. fgets(read, 200, fInputFile);
  330. #endif
  331. }
  332. }
  333. return kTRUE;
  334. }
  335. // ------------------------------------------------------------------------
  336. // ----- Private method ReadConverisonTable ---------------------------
  337. void MpdUrqmdGenerator::ReadConversionTable() {
  338. TString work = getenv("VMCWORKDIR");
  339. TString fileName = work + "/input/urqmd_pdg.dat";
  340. ifstream* pdgconv = new ifstream(fileName.Data());
  341. Int_t index = 0;
  342. Int_t pdgId = 0;
  343. while (!pdgconv->eof()) {
  344. index = pdgId = 0;
  345. *pdgconv >> index >> pdgId;
  346. fParticleTable[index] = pdgId;
  347. }
  348. pdgconv->close();
  349. delete pdgconv;
  350. cout << "-I MpdUrqmdGenerator: Particle table for conversion from "
  351. << "UrQMD loaded" << endl;
  352. }
  353. // ------------------------------------------------------------------------
  354. void MpdUrqmdGenerator::SetEventPlane(Double_t phiMin, Double_t phiMax) {
  355. fPhiMin = phiMin;
  356. fPhiMax = phiMax;
  357. fEventPlaneSet = kTRUE;
  358. }
  359. ClassImp(MpdUrqmdGenerator);