MpdLAQGSMGenerator.cxx 26 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981
  1. // -------------------------------------------------------------------------
  2. // ----- MpdLAQGSMGenerator source file -----
  3. // -------------------------------------------------------------------------
  4. /** MpdLAQGSMGenerator
  5. *@author Elena Litvinenko <litvin@nf.jinr.ru>
  6. *@version 15.02.2016
  7. *
  8. ** The MpdLAQGSMGenerator uses the ASCII input
  9. ** provided by K.Gudima LAQGSM event generator.
  10. **/
  11. #include "MpdLAQGSMGenerator.h"
  12. #include "MpdHypYPtGenerator.h" //AZ
  13. #include "MpdStack.h"
  14. #include "FairMCEventHeader.h"
  15. #include "FairPrimaryGenerator.h"
  16. #include "FairIon.h"
  17. #include "FairRunSim.h"
  18. #include "FairParticle.h"
  19. #include "TDatabasePDG.h"
  20. #include "TParticlePDG.h"
  21. #include "TParticle.h"
  22. #include "TSystem.h"
  23. #include "TParticle.h"
  24. #include "TVirtualMC.h"
  25. #include <iostream>
  26. using std::cout;
  27. using std::endl;
  28. using std::map;
  29. const Double_t kProtonMass = 0.938272013;
  30. const Double_t kNeutronMass = 0.93957;
  31. // ----- Default constructor ------------------------------------------
  32. MpdLAQGSMGenerator::MpdLAQGSMGenerator() :
  33. fInputFile(NULL),
  34. fGZInputFile(NULL)
  35. {}
  36. // ------------------------------------------------------------------------
  37. // ----- Standard constructor -----------------------------------------
  38. MpdLAQGSMGenerator::MpdLAQGSMGenerator(const char* fileName, const Bool_t use_collider_system, Int_t QGSM_format_ID,Int_t Max_Event_Number ) :
  39. FairGenerator(),
  40. fIonMap(),
  41. fInputFile(NULL),
  42. fGZInputFile(NULL)
  43. {
  44. //AZ memset(la_tab,0,sizeof(la_tab));
  45. fPDG=TDatabasePDG::Instance();
  46. fFileName = fileName;
  47. cout << "-I- MpdLAQGSMGenerator: Opening input file " << fileName << endl;
  48. TString sFileName = fileName;
  49. if (sFileName.Contains(".gz"))
  50. fGZ_input=1;
  51. else
  52. fGZ_input=0;
  53. if (fGZ_input)
  54. fGZInputFile = gzopen(fFileName, "rb");
  55. else
  56. fInputFile = fopen(fFileName, "r");
  57. if (( ! fInputFile ) && ( ! fGZInputFile ))
  58. Fatal("MpdLAQGSMGenerator","Cannot open input file.");
  59. fUseColliderSystem = use_collider_system;
  60. fQGSM_format_ID=QGSM_format_ID;
  61. Init(); //AZ
  62. cout << "-I- MpdLAQGSMGenerator: Looking for ions..." << endl;
  63. //AZ Int_t nIons = RegisterIons();
  64. // Int_t nIons = RegisterIons1();
  65. // Int_t nIons = RegisterIons(); // MG
  66. Int_t nIons = RegisterIons(Max_Event_Number); // EL
  67. cout << "-I- MpdLAQGSMGenerator: " << nIons << " ions registered."
  68. << endl;
  69. CloseInput();
  70. cout << "-I- MpdLAQGSMGenerator: Reopening input file " << fileName
  71. << endl;
  72. if (fGZ_input)
  73. fGZInputFile= gzopen(fFileName, "rb");
  74. else
  75. fInputFile = fopen(fFileName, "r");
  76. }
  77. // ------------------------------------------------------------------------
  78. // ----- Destructor ---------------------------------------------------
  79. MpdLAQGSMGenerator::~MpdLAQGSMGenerator()
  80. {
  81. CloseInput();
  82. Int_t nPart = fLa_tab.size();
  83. for (Int_t i = 0; i < nPart; ++i) delete fLa_tab[i];
  84. fLa_tab.clear();
  85. }
  86. // ------------------------------------------------------------------------
  87. // --------------Public method Init (fill list of known light particles)---
  88. void MpdLAQGSMGenerator::Init (const char *light_particles_filename)
  89. {
  90. TString fname, dd;
  91. if (!light_particles_filename) {
  92. dd = getenv("VMCWORKDIR");
  93. if (gSystem->OpenDirectory(dd+"/mpdgenerators"))
  94. dd += "/mpdgenerators";
  95. else {
  96. if (gSystem->OpenDirectory(dd+"/generators"))
  97. dd += "/generators";
  98. else
  99. dd = gSystem->WorkingDirectory();
  100. }
  101. fname = dd+ "/MpdLAQGSM_light_particles.dat";
  102. }
  103. else
  104. fname = light_particles_filename;
  105. Int_t i=0, i_val, i_max=84;
  106. Float_t a_val;
  107. char ss[300];
  108. FILE *fInputFile_light = fopen(fname.Data(), "r");
  109. //AZ while (!(feof(fInputFile_light))&&(i<i_max)) {
  110. while (!(feof(fInputFile_light))) {
  111. la_tab_t *lt = new la_tab_t;
  112. /*AZ
  113. fscanf(fInputFile_light, "%d", &i_val); la_tab[i].pdg=i_val;
  114. fscanf(fInputFile_light, "%d", &i_val); la_tab[i].Z=i_val;
  115. fscanf(fInputFile_light, "%d", &i_val); la_tab[i].lepton=i_val;
  116. fscanf(fInputFile_light, "%d", &i_val); la_tab[i].strange=i_val;
  117. fscanf(fInputFile_light, "%d", &i_val); la_tab[i].A=i_val;
  118. fscanf(fInputFile_light, "%f", &a_val); la_tab[i].mass=a_val;
  119. fscanf(fInputFile_light, "%s", ss); strncpy(la_tab[i].name,ss,10);
  120. i++;
  121. */
  122. i = fscanf(fInputFile_light, "%d", &i_val); lt->pdg=i_val;
  123. i = fscanf(fInputFile_light, "%d", &i_val); lt->Z=i_val;
  124. i = fscanf(fInputFile_light, "%d", &i_val); lt->lepton=i_val;
  125. i = fscanf(fInputFile_light, "%d", &i_val); lt->strange=i_val;
  126. i = fscanf(fInputFile_light, "%d", &i_val); lt->A=i_val;
  127. i = fscanf(fInputFile_light, "%f", &a_val); lt->mass=a_val;
  128. i = fscanf(fInputFile_light, "%s", ss); strncpy(lt->name,ss,10);
  129. fLa_tab.push_back(lt);
  130. if(!i) cerr<<__FILE__<<__func__<<" something went wrong in reading input.";
  131. }
  132. fclose(fInputFile_light);
  133. }
  134. // ------------------------------------------------------------------------
  135. Bool_t MpdLAQGSMGenerator::general_fgets (char *ss, Int_t nn, FILE* p)
  136. {
  137. Bool_t finished=0;
  138. if (fGZ_input) {
  139. gzgets(fGZInputFile,ss,nn);
  140. finished = (gzeof(fGZInputFile));
  141. }
  142. else {
  143. if(!fgets(ss,nn,fInputFile))
  144. cerr<<__FILE__<<__func__<<" something went wrong in reading input.";
  145. finished = ( feof(fInputFile) );
  146. }
  147. return finished;
  148. }
  149. Bool_t MpdLAQGSMGenerator::GetEventHeader(char *ss)
  150. {
  151. Bool_t finished=0;
  152. // if (fGZ_input) {
  153. // gzgets(fGZInputFile,ss,250);
  154. // finished = (gzeof(fGZInputFile));
  155. // }
  156. // else {
  157. // fgets(ss,250,fInputFile);
  158. // finished = ( feof(fInputFile) );
  159. // }
  160. finished = general_fgets (ss,250,fInputFile);
  161. // If end of input file is reached : close it and abort run
  162. if ( finished ) {
  163. cout << "-I- MpdLAQGSMGenerator: End of input file reached " << endl;
  164. CloseInput();
  165. return kFALSE;
  166. }
  167. Bool_t wrong_file=kFALSE;
  168. TString tss(ss);
  169. if (tss.Contains("QGSM")) { // 0
  170. general_fgets(ss,250,fInputFile);
  171. tss=ss;
  172. Int_t lines=0;
  173. while (!(tss.Contains("particles, b,bx,by"))) {
  174. if (tss.Contains("QGSM format ID")) {
  175. sscanf(ss,"QGSM format ID=%d",&fQGSM_format_ID);
  176. cout << "QGSM format ID " << fQGSM_format_ID << endl;
  177. if (fQGSM_format_ID==2) { // correction of incorrect format_ID in some data files
  178. fpos_t file_loc;
  179. z_off_t gz_file_loc;
  180. if (fGZ_input)
  181. gz_file_loc = gztell(fGZInputFile);
  182. else
  183. fgetpos(fInputFile, &file_loc);
  184. for (int k=0;k<5;k++) general_fgets(ss,250,fInputFile);
  185. if (strlen(ss)>90)
  186. fQGSM_format_ID=3;
  187. cout << "QGSM format ID now " << fQGSM_format_ID << endl;
  188. if (fGZ_input)
  189. gzseek(fGZInputFile,gz_file_loc,SEEK_SET);
  190. else
  191. fsetpos(fInputFile, &file_loc);
  192. }
  193. }
  194. general_fgets(ss,250,fInputFile);
  195. tss=ss;
  196. lines++;
  197. if ((fQGSM_format_ID>=2)&&(lines>=4))
  198. return kTRUE;
  199. if (lines>25) {
  200. wrong_file=kTRUE;
  201. break;
  202. }
  203. }
  204. }
  205. else {
  206. if (fQGSM_format_ID>=2)
  207. return kTRUE;
  208. if (!tss.Contains("particles, b,bx,by")) {
  209. general_fgets(ss,250,fInputFile);
  210. tss=ss;
  211. }
  212. }
  213. if (wrong_file) {
  214. cout << "-E- MpdLAQGSMGenerator: Wrong file header! " << endl;
  215. CloseInput();
  216. return kFALSE;
  217. }
  218. if (!tss.Contains("event")) { // 3
  219. if (fQGSM_format_ID>=2)
  220. return kTRUE;
  221. else {
  222. cout << "-E- MpdLAQGSMGenerator: Wrong event header!" << endl;
  223. CloseInput();
  224. return kFALSE;
  225. }
  226. }
  227. char tmp[250];
  228. switch (fQGSM_format_ID) {
  229. case 0:
  230. case 1:
  231. general_fgets(tmp,250,fInputFile);
  232. break;
  233. case 2:
  234. case 3:
  235. general_fgets(ss,250,fInputFile);
  236. general_fgets(ss,250,fInputFile);
  237. break;
  238. case 4:
  239. general_fgets(ss,250,fInputFile);
  240. general_fgets(ss,250,fInputFile);
  241. break;
  242. default:
  243. break;
  244. }
  245. return kTRUE;
  246. }
  247. // ------------------------------------------------------------------------
  248. //AZ ----- Public method SkipEvents --------------------------------------
  249. Bool_t MpdLAQGSMGenerator::SkipEvents(Int_t count) {
  250. if (count<=0) return kTRUE;
  251. // Check for input file
  252. if (fGZ_input) {
  253. if ( ! fGZInputFile ) {
  254. cout << "-E- MpdLAQGSMGenerator: Input file not open!" << endl;
  255. return kFALSE;
  256. }
  257. }
  258. else {
  259. if ( ! fInputFile ) {
  260. cout << "-E- MpdLAQGSMGenerator: Input file not open!" << endl;
  261. return kFALSE;
  262. }
  263. }
  264. Int_t eventId = 0, nTracks = 0;
  265. for(Int_t ii=0; ii<count;ii++) {
  266. char ss[252];
  267. if (!GetEventHeader(ss))
  268. return kFALSE;
  269. sscanf(ss," %d %d", &eventId, &nTracks);
  270. //cout << ii << " " << eventId << endl;
  271. for (Int_t itrack=0; itrack<nTracks; itrack++) {
  272. general_fgets(ss,250,fInputFile);
  273. }
  274. }
  275. return kTRUE;
  276. }
  277. // ------------------------------------------------------------------------
  278. // ----- Public method ReadEvent --------------------------------------
  279. Bool_t MpdLAQGSMGenerator::ReadEvent(FairPrimaryGenerator* primGen) {
  280. // Check for input file
  281. if (fGZ_input) {
  282. if ( ! fGZInputFile ) {
  283. cout << "-E- MpdLAQGSMGenerator: Input file not open!" << endl;
  284. return kFALSE;
  285. }
  286. }
  287. else {
  288. if ( ! fInputFile ) {
  289. cout << "-E- MpdLAQGSMGenerator: Input file not open!" << endl;
  290. return kFALSE;
  291. }
  292. }
  293. // ---> Check for primary generator
  294. if ( ! primGen ) {
  295. cout << "-E- MpdLAQGSMGenerator::ReadEvent: "
  296. << "No PrimaryGenerator!" << endl;
  297. return kFALSE;
  298. }
  299. //AZ if (la_tab[0].pdg<=0)
  300. if (fLa_tab[0]->pdg<=0)
  301. Init();
  302. // Define event variable to be read from file
  303. Int_t eventId = 0;
  304. Int_t nTracks = 0;
  305. Float_t b = 0., bx=0, by=0;
  306. // Define track variables to be read from file
  307. Int_t
  308. iCharge = 0,
  309. iLeptonic = 0,
  310. iStrange = 0,
  311. iBarionic=0, iCode=0, io1=0, io2=0, io3=0 //,
  312. //A=0,
  313. //Z=0
  314. ;
  315. Float_t
  316. px = 0.,
  317. py = 0.,
  318. pz = 0.,
  319. pz1 = 0,
  320. pza1 = 0,
  321. // Ekin = 0,
  322. // Ekin1 = 0,
  323. mass = 0, pol0 = 0, pol1 = 0;
  324. Int_t PDG=0;
  325. char ss[252];
  326. char ionName[30];
  327. if (!GetEventHeader(ss))
  328. return kFALSE;
  329. //Double_t projA,projZ,projEex, targA, targZ, targEex;
  330. Float_t projA,projZ,projEex, targA, targZ, targEex;
  331. //cout <<"ss " <<ss <<endl;
  332. int bpos;
  333. sscanf(ss," %d %d %n", &eventId, &nTracks, &bpos);
  334. //cout <<"eventId " <<eventId <<" " <<nTracks <<" " <<fQGSM_format_ID <<endl;
  335. if (!fQGSM_format_ID)
  336. sscanf(&(ss[bpos]),"%g %g %g %g %g %g %g %g %g", &b, &bx, &by, &projA,&projZ,&projEex, &targA, &targZ, &targEex );
  337. else
  338. sscanf(&(ss[bpos]),"%g %g %g", &b, &bx, &by);
  339. cout << "-I- MpdLAQGSMGenerator::ReadEvent: Event " << eventId
  340. << ", b = " << b << " fm, multiplicity " << nTracks << endl;
  341. TVector2 bb; // MG
  342. bb.Set(bx,by); // MG
  343. // Set event id and impact parameter in MCEvent if not yet done
  344. FairMCEventHeader* event = primGen->GetEvent();
  345. if ( event && (! event->IsSet()) ) {
  346. event->SetEventID(eventId-1);
  347. event->SetB(b);
  348. event->SetRotZ(bb.Phi()); // MG
  349. event->MarkSet(kTRUE);
  350. }
  351. /*
  352. //AZ - Loop over tracks to find hypernuclei
  353. fpos_t position;
  354. fgetpos (fInputFile, &position); // save position
  355. Int_t iHyperOk = 0;
  356. for (Int_t itrack = 0; itrack < nTracks; itrack++) {
  357. fgets(ss,250,fInputFile);
  358. sscanf(ss," %d %d %d %d %d", &iCharge,&iLeptonic,&iStrange,&iBarionic, &iCode);
  359. if (iStrange == -1 && iCharge == 1 && iBarionic == 3) {
  360. // select only H3L
  361. iHyperOk = 1;
  362. break;
  363. }
  364. // Select H4LL
  365. //if (iStrange == -2 && iCharge == 1 && iBarionic == 4) {
  366. //iHyperOk = 1;
  367. //break;
  368. //}
  369. }
  370. if (iHyperOk) fsetpos (fInputFile, &position);
  371. //else nTracks = 0; // do not store any track
  372. else {
  373. // Add particle manually
  374. TParticlePDG* part = fPDG->GetParticle("H3L");
  375. Int_t pdgType = part->PdgCode();
  376. static MpdHypYPtGenerator *hypGen = 0x0;
  377. if (hypGen == 0x0) {
  378. hypGen = new MpdHypYPtGenerator(pdgType,1);
  379. hypGen->SetDistributionPt(0.163); // temperature - H3L: NICA @ 5 GeV
  380. hypGen->SetDistributionY(0.0, 0.373); // y0, Sigma_y - H3L: NICA @ 5 GeV
  381. hypGen->Init();
  382. }
  383. hypGen->ReadEvent(primGen);
  384. fsetpos (fInputFile, &position);
  385. }
  386. // AZ
  387. */
  388. // Loop over tracks in the current event
  389. for (Int_t itrack=0; itrack<nTracks; itrack++) {
  390. general_fgets(ss,250,fInputFile);
  391. if (fQGSM_format_ID<3)
  392. sscanf(ss," %d %d %d %d %d", &iCharge,&iLeptonic,&iStrange,&iBarionic, &iCode);
  393. else{
  394. if(fQGSM_format_ID==4){
  395. //cout << ss << endl;
  396. sscanf(ss," %d %d %d %d %d %g %g %g %g %g %g", &iCharge,&iLeptonic,&iStrange,&iBarionic, &PDG, &px, &py, &pz, &pz1, &pza1, &mass);
  397. //printf("!!!!!!!!!!!!!!!!!!!!!reading: %d %d %d %d %d %g %g %g %g %g %g\n", iCharge,iLeptonic,iStrange,iBarionic, PDG, px, py, pz, pz1, pza1, mass);
  398. //primGen->AddTrack(PDG, px, py, pz, 0., 0., 0.);
  399. if (fUseColliderSystem)
  400. primGen->AddTrack(PDG, px, py, pz, 0., 0., 0.);
  401. else
  402. primGen->AddTrack(PDG, px, py, pz1, 0., 0., 0.);
  403. }
  404. else
  405. sscanf(ss," %d %d %d %d %d %d %d", &iCharge,&iLeptonic,&iStrange,&iBarionic, &io1, &io2, &io3);
  406. }
  407. if(fQGSM_format_ID!=4){
  408. Int_t p_num = 21;
  409. if (fQGSM_format_ID==3)
  410. p_num=36;
  411. if (!fQGSM_format_ID)
  412. sscanf(&(ss[p_num]),"%g%g%g%g%g",&px,&py,&pz,&pz1,&mass);
  413. else if (strlen(ss) > 105)
  414. sscanf(&(ss[p_num]),"%g%g%g%g%g%g%g",&px,&py,&pz,&pz1,&pol0,&pol1,&mass); //AZ - polarization data
  415. else
  416. sscanf(&(ss[p_num]),"%g%g%g%g%g%g",&px,&py,&pz,&pz1,&pza1,&mass);
  417. if (FindParticle (iCharge, iStrange, iLeptonic, iBarionic, mass, PDG, ionName)) {
  418. if (!fPDG->GetParticle(PDG)) {
  419. cout << "-W- MpdLAQGSMGenerator::ReadEvent: Cannot find " << PDG << " "
  420. << ionName << " in database!" << endl;
  421. continue;
  422. }
  423. // Give track to PrimaryGenerator
  424. if (fUseColliderSystem)
  425. primGen->AddTrack(PDG, px, py, pz, 0., 0., 0.);
  426. else
  427. primGen->AddTrack(PDG, px, py, pz1, 0., 0., 0.);
  428. //AZ - add polarization for light hyperons (pdg == 31** and 32**)
  429. //*
  430. //if (PDG == 3122) {
  431. if (TMath::Abs(PDG/100) == 31 || TMath::Abs(PDG/100) == 32) {
  432. Int_t nTr = gMC->GetStack()->GetNtrack();
  433. TParticle *part = ((MpdStack*)gMC->GetStack())->GetParticle(nTr-1);
  434. Double_t pol[3] = {pol1,0,0}; // polarization after rescattering
  435. pol[1] = TMath::Sqrt (1.0 - pol[0] * pol[0]);
  436. part->SetPolarisation(pol[0], pol[1], pol[2]);
  437. }
  438. //*/
  439. }//if (FindParticle (iCharge
  440. }//if(fQGSM_format_ID!=4)
  441. }
  442. return kTRUE;
  443. }
  444. // ------------------------------------------------------------------------
  445. // ----- Private method CloseInput ------------------------------------
  446. void MpdLAQGSMGenerator::CloseInput() {
  447. if ( fInputFile ) {
  448. if (!fGZ_input)
  449. fclose(fInputFile);
  450. fInputFile = NULL;
  451. }
  452. if ( fGZInputFile ) {
  453. if (fGZ_input)
  454. gzclose(fGZInputFile);
  455. fGZInputFile = NULL;
  456. }
  457. }
  458. // ------------------------------------------------------------------------
  459. Int_t MpdLAQGSMGenerator::CreatePdgCode(Int_t Z, Int_t A, Int_t Strange,Int_t user)
  460. {
  461. // 10LZZZAAAI
  462. const Int_t kion =1000000000;
  463. const Int_t kion_strange = kion+10000000;
  464. const Int_t kspecial =50000000;
  465. if (Strange) {
  466. Int_t ss=abs(Strange);
  467. if (Z) {
  468. if (ss==1)
  469. return (kion_strange+Z*10000+A*10);
  470. else
  471. return (kion+10000000*ss+Z*10000+A*10);
  472. }
  473. else
  474. return (kspecial + A*10);
  475. }
  476. else {
  477. return (kion+Z*10000+A*10+user);
  478. }
  479. }
  480. Bool_t MpdLAQGSMGenerator::general_feof (void *p)
  481. {
  482. Bool_t finished=0;
  483. if (fGZ_input)
  484. finished = (gzeof(fGZInputFile));
  485. else
  486. finished = ( feof(fInputFile) );
  487. return finished;
  488. }
  489. // ----- Private method RegisterIons ----------------------------------
  490. Int_t MpdLAQGSMGenerator::RegisterIons(Int_t Max_Event_Number) {
  491. Int_t nIons = 0;
  492. // Int_t eventId, nTracks, iPid, iMass, iCharge;
  493. // Double_t pBeam, b, px, py, pz;
  494. Int_t eventId = 0;
  495. Int_t nTracks = 0;
  496. // Float_t b = 0., bx=0, by=0;
  497. // Define track variables to be read from file
  498. Int_t
  499. iCharge = 0,
  500. iLeptonic = 0,
  501. iStrange = 0,
  502. iBarionic=0,
  503. iCode =0;
  504. Float_t
  505. px = 0.,
  506. py = 0.,
  507. pz = 0.,
  508. pz1 = 0,
  509. pza1 = 0,
  510. // Ekin = 0,
  511. // Ekin1 = 0,
  512. mass = 0,
  513. excEnergy=0;
  514. char
  515. ss[252],
  516. buf_ionName[30];
  517. Int_t
  518. A=0, Z=0, Q=0;
  519. fIonMap.clear();
  520. if (!GetEventHeader(ss))
  521. return 0;
  522. Bool_t geant4_bug_case=1;
  523. Int_t PDG;
  524. TString ionName;
  525. while ( (!general_feof(fInputFile)) && ((!Max_Event_Number)||(eventId<Max_Event_Number)) ) {
  526. sscanf(ss," %d %d", &eventId, &nTracks);
  527. if ( general_feof(fInputFile) ) continue;
  528. for (Int_t iTrack=0; iTrack<nTracks; iTrack++) {
  529. general_fgets(ss,250,fInputFile);
  530. sscanf(ss," %d %d %d %d %d", &iCharge,&iLeptonic,&iStrange,&iBarionic,&iCode);
  531. // if (!fQGSM_format_ID)
  532. // sscanf(&(ss[21]),"%g%g%g%g%g",&px,&py,&pz,&pz1,&mass);
  533. // else
  534. // sscanf(&(ss[21]),"%g%g%g%g%g%g",&px,&py,&pz,&pz1,&pza1,&mass);
  535. Int_t p_num = 21;
  536. if (fQGSM_format_ID>2)
  537. p_num=36;
  538. if (!fQGSM_format_ID) {
  539. sscanf(&(ss[p_num]),"%g%g%g%g%g",&px,&py,&pz,&pz1,&mass);
  540. //if(iCharge>20)
  541. // cout << "format 1 " << mass << " " << fQGSM_format_ID << endl;
  542. }
  543. else{
  544. sscanf(&(ss[p_num]),"%g%g%g%g%g%g",&px,&py,&pz,&pz1,&pza1,&mass);
  545. //if(iCharge>20)
  546. // cout << "format 2 " << px << " " << py << " " << pz << " " << pz1 << " " << pza1 << " "<< mass << " " << fQGSM_format_ID << endl;
  547. }
  548. if (iBarionic>1){ // ion
  549. A = iBarionic;
  550. Z = iCharge;
  551. Q = iCharge;
  552. if (A>9)
  553. geant4_bug_case=0;
  554. //if(iCharge>20)
  555. // cout << "??????????????????1 " << iCharge << " " << iStrange << " " << iLeptonic << " " << iBarionic<< " " << mass << endl;
  556. if (FindParticle (iCharge, iStrange, iLeptonic, iBarionic, mass, PDG, buf_ionName)) {
  557. //if(iCharge>20)
  558. // cout << "??????????????????2 " << iCharge << " " << iStrange << " " << iLeptonic << " " << iBarionic<< " " << mass << endl;
  559. if (fPDG->GetParticle(PDG))
  560. continue;
  561. // cout << "!!!!!!!!!!!!!!!!!!!! " << PDG << endl;
  562. ionName = buf_ionName;
  563. if (!iStrange) { // normal ion
  564. if(Z>2 || (Z==2 && A>4)) { // MG
  565. if ( fIonMap.find(ionName) == fIonMap.end() ) { // new ion
  566. // excEnergy = fabs(mass - kProtonMass*Z -kNeutronMass*(iBarionic-Z));
  567. FairIon* ion = new FairIon(ionName, Z, A, Q) ; // , excEnergy,mass); // from MG
  568. fIonMap[ionName] = ion;
  569. nIons++;
  570. }
  571. }
  572. }
  573. else { // hyper
  574. if (CreateNucleus (Z, mass, PDG, buf_ionName)) {
  575. // FairParticle *pa=new FairParticle(buf_ionName,Z, iBarionic, abs(iStrange), mass,Z,0,2.6e-10);
  576. /*
  577. TParticle *tpa= new TParticle();
  578. if (tpa) {
  579. tpa->SetPdgCode(PDG);
  580. FairParticle *pa=new FairParticle(PDG,tpa);
  581. if (pa)
  582. (FairRunSim::Instance())->AddNewParticle(pa);
  583. }
  584. */
  585. if (Z) { //AZ
  586. FairParticle *pa = new FairParticle(PDG, buf_ionName, kPTHadron, mass, Z, 2.6e-10, "hadron");
  587. if (pa) FairRunSim::Instance()->AddNewParticle(pa);
  588. } else {
  589. FairParticle *pa = new FairParticle(PDG, buf_ionName, kPTNeutron, mass, Z, 2.6e-10, "special");
  590. if (pa) FairRunSim::Instance()->AddNewParticle(pa);
  591. } //AZ
  592. }
  593. }
  594. }
  595. } // ion
  596. } // track loop
  597. if (!GetEventHeader(ss))
  598. break;
  599. } // event loop
  600. if (geant4_bug_case) { // artificial heavy ion
  601. mass = 9.98;
  602. excEnergy = fabs(mass - (kProtonMass+kNeutronMass)*5);
  603. TString ionName_dummy="Ion_10_5";
  604. if ( fIonMap.find(ionName_dummy) == fIonMap.end() ) {
  605. FairIon* ion = new FairIon(ionName_dummy, 5,10, 5, excEnergy,mass);
  606. fIonMap[ionName_dummy] = ion;
  607. nIons++;
  608. }
  609. }
  610. FairRunSim* run = FairRunSim::Instance();
  611. map<TString, FairIon*>::iterator mapIt;
  612. for (mapIt=fIonMap.begin(); mapIt!=fIonMap.end(); ++mapIt) {
  613. FairIon* ion = (*mapIt).second;
  614. run->AddNewIon(ion);
  615. }
  616. return nIons;
  617. }
  618. // ----- Private method RegisterIons1 ----------------------------------
  619. Int_t MpdLAQGSMGenerator::RegisterIons1()
  620. {
  621. // Read file with particle data
  622. Int_t nIons = 0, nPart = fLa_tab.size();
  623. fIonMap.clear();
  624. Bool_t geant4_bug_case = 1;
  625. Int_t PDG;
  626. TString ionName;
  627. Float_t mass = 0.0, excEnergy = 0.0;
  628. char buf_ionName[30];
  629. for (Int_t i = 0; i < nPart; ++i) {
  630. Int_t Z = fLa_tab[i]->Z, iCharge = Z, Q = Z;
  631. Int_t iLeptonic = fLa_tab[i]->lepton;
  632. Int_t iStrange = fLa_tab[i]->strange;
  633. Int_t A = fLa_tab[i]->A, iBarionic = A;
  634. mass = fLa_tab[i]->mass;
  635. if (iBarionic > 1) {
  636. // ion
  637. //if (A > 9) geant4_bug_case = 0;
  638. if (FindParticle (iCharge, iStrange, iLeptonic, iBarionic, mass, PDG, buf_ionName)) {
  639. if (fPDG->GetParticle(PDG)) continue;
  640. ionName = buf_ionName;
  641. if (!iStrange) {
  642. // normal ion
  643. if ( fIonMap.find(ionName) == fIonMap.end() ) { // new ion
  644. excEnergy = fabs(mass - kProtonMass*Z -kNeutronMass*(iBarionic-Z));
  645. FairIon* ion = new FairIon(ionName, Z, A, Q, excEnergy,mass);
  646. fIonMap[ionName] = ion;
  647. nIons++;
  648. }
  649. } else {
  650. // hyper
  651. if (CreateNucleus (Z, mass, PDG, buf_ionName)) {
  652. // FairParticle *pa=new FairParticle(buf_ionName,Z, iBarionic, abs(iStrange), mass,Z,0,2.6e-10);
  653. if (Z) {
  654. FairParticle *pa = new FairParticle(PDG, buf_ionName, kPTHadron, mass, Z, 2.6e-10, "hadron");
  655. if (pa) FairRunSim::Instance()->AddNewParticle(pa);
  656. } else {
  657. FairParticle *pa = new FairParticle(PDG, buf_ionName, kPTNeutron, mass, Z, 2.6e-10, "special");
  658. if (pa) FairRunSim::Instance()->AddNewParticle(pa);
  659. }
  660. }
  661. }
  662. }
  663. } // if (iBarionic > 1)
  664. } // for (Int_t i = 0; i < nPart;
  665. if (geant4_bug_case) { // artificial heavy ion
  666. mass = 9.98;
  667. excEnergy = fabs(mass - (kProtonMass+kNeutronMass)*5);
  668. TString ionName_dummy="Ion_10_5";
  669. if ( fIonMap.find(ionName_dummy) == fIonMap.end() ) {
  670. FairIon* ion = new FairIon(ionName_dummy, 5,10, 5, excEnergy,mass);
  671. fIonMap[ionName_dummy] = ion;
  672. nIons++;
  673. }
  674. }
  675. FairRunSim* run = FairRunSim::Instance();
  676. map<TString, FairIon*>::iterator mapIt;
  677. for (mapIt=fIonMap.begin(); mapIt!=fIonMap.end(); ++mapIt) {
  678. FairIon* ion = (*mapIt).second;
  679. run->AddNewIon(ion);
  680. }
  681. return nIons;
  682. }
  683. // ------------------------------------------------------------------------
  684. Bool_t MpdLAQGSMGenerator::FindParticle (Int_t Z, Int_t strange, Int_t lepton, Int_t A, Float_t mass, Int_t &PDG, char name[11])
  685. {
  686. Int_t
  687. n_mass,
  688. result=0;
  689. if (mass<6) { // use table
  690. Int_t
  691. l_min=-1, // line number first
  692. l_max=-1; // line number last
  693. n_mass = (Int_t) (mass/0.1);
  694. switch (n_mass) {
  695. case 1: l_min=8; l_max=12; // MU+,MU-,PI0,PI+,PI-
  696. break;
  697. case 9: l_min=28; l_max=32; // P,AP,AN,N,ETAP
  698. break;
  699. case 4: l_min=13; l_max=18; // K+,K-,AK0,KL,KS,K0
  700. break;
  701. case 18: l_min=66; l_max=66; // d
  702. break;
  703. case 5: l_min=19; l_max=19; // ETA
  704. break;
  705. case 11: l_min=34; l_max=41; // AL,L,S+,AS+,AS0,S0,AS-,S-
  706. break;
  707. case 28: l_min=68; l_max=69; // He3,t
  708. break;
  709. case 37: l_min=71; l_max=71; // He4
  710. break;
  711. case 13: l_min=50; l_max=59; // AXI0,XI0,AXI-,XI-,AS*0,S*0,S*+,AS*+,S*-,AS*-
  712. break;
  713. case 29: l_min=70; l_max=70; // H3L
  714. break;
  715. case 22:
  716. if (!Z)
  717. { l_min=67; l_max=67;} // LL
  718. else
  719. { l_min=81; l_max=81;} // KsiNP
  720. break;
  721. case 39: l_min=72; l_max=73; // He4L,H4L
  722. break;
  723. case 48: l_min=75; l_max=75; // He5L
  724. break;
  725. case 16: l_min=64; l_max=65; // AOM-,OM-
  726. break;
  727. case 41: l_min=74; l_max=74; // H4LL
  728. break;
  729. case 50: l_min=76; l_max=77; // He5LL,H5LL
  730. break;
  731. case 59: l_min=78; l_max=78; // He6LL
  732. break;
  733. case 0: l_min=1; l_max=7; // GM,ANUM,ANUE,NUE,NUM,E+,E-
  734. break;
  735. case 7: l_min=20; l_max=23; // RHO+,RHO-,RHO0,OMEG
  736. break;
  737. case 8: l_min=24; l_max=27; // K*+,K*-,AK*0,K*0
  738. break;
  739. case 10: l_min=33; l_max=33; // PHI
  740. break;
  741. case 12: l_min=42; l_max=49; // DL++,ADL++,DL+,DL-,ADL0,DL0,ADL-,ADL+
  742. break;
  743. case 15: l_min=60; l_max=63; // AXI*0,XI*0,AXI*-,XI*-
  744. break;
  745. case 20: l_min=79; l_max=80; // LN, LP
  746. break;
  747. case 24: l_min=82; l_max=83; // KsiL0, KsiL1
  748. break;
  749. case 26: l_min=84; l_max=84; // KsiKsi
  750. break;
  751. default: // unknown particle
  752. break;
  753. }
  754. if (l_min>0){
  755. for (int i=l_min-1; i<l_max; i++ ) {
  756. /*AZ
  757. if (Z!=la_tab[i].Z) continue;
  758. if (strange!=la_tab[i].strange) continue;
  759. if (lepton!=la_tab[i].lepton) continue;
  760. if (A!=la_tab[i].A) continue;
  761. if (mass<(la_tab[i].mass*0.99)) continue;
  762. if (mass>(la_tab[i].mass*1.01)) continue;
  763. PDG=la_tab[i].pdg;
  764. strncpy(name,la_tab[i].name,10);
  765. */
  766. if (Z != fLa_tab[i]->Z) continue;
  767. if (strange != fLa_tab[i]->strange) continue;
  768. if (lepton != fLa_tab[i]->lepton) continue;
  769. if (A != fLa_tab[i]->A) continue;
  770. if (mass < (fLa_tab[i]->mass*0.99)) continue;
  771. if (mass > (fLa_tab[i]->mass*1.01)) continue;
  772. PDG = fLa_tab[i]->pdg;
  773. strncpy(name,fLa_tab[i]->name,10);
  774. result=1;
  775. break;
  776. }
  777. }
  778. }
  779. else { // must be ion
  780. if (A>1) {
  781. PDG = CreatePdgCode(Z,A,strange,0 ); // 1); // MG
  782. sprintf(name,"Ion_%d_%d", A,Z);
  783. result=1;
  784. }
  785. }
  786. return (result);
  787. }
  788. // ------------------------------------------------------------------------
  789. Bool_t MpdLAQGSMGenerator::CreateNucleus (Int_t Z, Float_t mass, Int_t pdgCode, char pdgName[11])
  790. {
  791. // const Double_t kAu2Gev=0.9314943228; // constants from TGeant3::AddParticlesToPdgDataBase code
  792. const Double_t khSlash = 1.0545726663e-27; // constants from TGeant3::AddParticlesToPdgDataBase code
  793. const Double_t kErg2Gev = 1/1.6021773349e-3;
  794. const Double_t khShGev = khSlash*kErg2Gev;
  795. const Double_t kYear2Sec = 3600*24*365.25;
  796. Double_t width = khShGev/(12.33*kYear2Sec); // ? EL
  797. cout << "-I- MpdLAQGSMGenerator::CreateNucleus " << pdgName << " " << pdgCode << endl;
  798. TParticlePDG* p=0;
  799. if (Z)
  800. p = fPDG->AddParticle(pdgName,pdgName,mass,kFALSE,
  801. width,Z*3,"nucleus",pdgCode,0,kPTHadron);
  802. else
  803. p = fPDG->AddParticle(pdgName,pdgName,mass,kFALSE,
  804. 6.58211889e-25/2.6e-10,0,
  805. // "special",pdgCode,0,kPTHadron);
  806. "special",pdgCode,0,kPTNeutron); //AZ!!!
  807. if (p)
  808. p->Print();
  809. return (p);
  810. }
  811. // ------------------------------------------------------------------------
  812. ClassImp(MpdLAQGSMGenerator)