MpdDCMSMMGenerator.cxx 19 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507
  1. // -------------------------------------------------------------------------
  2. // ----- MpdDCMSMMGenerator source file -----
  3. // ----- Created 27-AUG-2019 by Igor Rufanov -----
  4. // -------------------------------------------------------------------------
  5. #include "MpdDCMSMMGenerator.h"
  6. #include "FairPrimaryGenerator.h"
  7. #include "FairMCEventHeader.h"
  8. #include "TRandom.h"
  9. #include "TDatabasePDG.h"
  10. #include "TParticlePDG.h"
  11. #include "FairIon.h"
  12. #include "FairParticle.h"
  13. #include "FairRunSim.h"
  14. using namespace std;
  15. using namespace TMath;
  16. // ----- Default constructor ------------------------------------------
  17. MpdDCMSMMGenerator::MpdDCMSMMGenerator()
  18. : FairGenerator(),
  19. fInputFile(NULL),
  20. fFileName(NULL){
  21. }
  22. // ------------------------------------------------------------------------
  23. // ----- Standard constructor -----------------------------------------
  24. MpdDCMSMMGenerator::MpdDCMSMMGenerator(const char* fileName)
  25. : FairGenerator(),
  26. fInputFile(NULL),
  27. fFileName(fileName),
  28. fSpectatorsON(kFALSE),
  29. fFixedTarget(kFALSE),
  30. fGammaCM(1.),
  31. fBetaCM(0.) {
  32. cout << "-I MpdDCMSMMGenerator: Opening input file " << fFileName << endl;
  33. #ifdef GZIP_SUPPORT
  34. fInputFile = gzopen(fFileName, "rb");
  35. #else
  36. fInputFile = fopen(fFileName, "r");
  37. #endif
  38. if (!fInputFile) {
  39. Fatal("MpdDCMSMMGenerator", "Cannot open input file.");
  40. exit(1);
  41. }
  42. // read file header
  43. char read[200];
  44. Int_t A1,Z1,A2,Z2;
  45. Double_t T0,sqS;
  46. for( Int_t i=0; i<3; i++) {
  47. #ifdef GZIP_SUPPORT
  48. char* ch= gzgets( fInputFile, read, 200);
  49. #else
  50. char* ch= fgets(read, 200, fInputFile);
  51. #endif
  52. cout<<"-I MpdDCMSMMGenerator:"<<read;
  53. if( i==0) {
  54. string str0(read);
  55. A1= stoi( str0.substr( str0.find("A1=")+3, 3));
  56. Z1= stoi( str0.substr( str0.find("Z1=")+3, 3));
  57. A2= stoi( str0.substr( str0.find("A2=")+3, 3));
  58. Z2= stoi( str0.substr( str0.find("Z2=")+3, 3));
  59. //Int_t A1= 0;
  60. } else if( i==1) {
  61. string str0(read);
  62. T0= stof( str0.substr( str0.find("T0=")+3, 8));
  63. sqS= stof( str0.substr( str0.find("sqrt(s)=")+8, 8));
  64. //Double_t mp=0.938272;
  65. Double_t mp=0.940; // to obtain equivalence of read "sqrt(s)=" and calculated below mCMS
  66. Double_t e=T0+mp;
  67. Double_t p= pow( e*e-mp*mp, 0.5);
  68. Double_t mCMS= pow( 2.*mp*(e+mp), 0.5); cout<<"mCMS="<<mCMS<<endl;
  69. fGammaCM=(e+mp)/mCMS;
  70. fBetaCM= pow( 1.-1./(fGammaCM*fGammaCM), 0.5);
  71. }
  72. }
  73. cout<<"-I MpdDCMSMMGenerator: A1="<<A1<<" Z1="<<Z1<<" A2="<<A2<<" Z2="<<Z2<<" T0="<<T0<<" sqS="<<sqS<<endl;
  74. if( fSpectatorsON) {
  75. Int_t n= RegisterIons();
  76. }
  77. }
  78. // ------------------------------------------------------------------------
  79. // ----- Destructor ---------------------------------------------------
  80. MpdDCMSMMGenerator::~MpdDCMSMMGenerator() {
  81. if (fInputFile) {
  82. #ifdef GZIP_SUPPORT
  83. gzclose(fInputFile);
  84. #else
  85. fclose(fInputFile);
  86. #endif
  87. fInputFile = NULL;
  88. }
  89. cout<<"Leave Destructor of MpdDCMSMMGenerator"<<endl;
  90. }
  91. // ------------------------------------------------------------------------
  92. // ----- Public method ReadEvent --------------------------------------
  93. Bool_t MpdDCMSMMGenerator::ReadEvent(FairPrimaryGenerator* primGen) {
  94. // ---> Check for input file
  95. // cout<<"MpdDCMSMMGenerator::ReadEvent -------------------------"<<endl;
  96. if (!fInputFile) {
  97. cout << "-E MpdDCMSMMGenerator: Input file not open! " << endl;
  98. return kFALSE;
  99. }
  100. // ---> Check for primary generator
  101. if (!primGen) {
  102. cout << "-E- MpdDCMSMMGenerator::ReadEvent: "
  103. << "No PrimaryGenerator!" << endl;
  104. return kFALSE;
  105. }
  106. char read[80];
  107. #ifdef GZIP_SUPPORT
  108. char* ch= gzgets( fInputFile, read, 80);
  109. #else
  110. char* ch= fgets(read, 80, fInputFile);
  111. #endif
  112. Int_t evnr=0; Float_t b, bimpX, bimpY;
  113. sscanf(read, "%d %f %f %f", &evnr, &b, &bimpX, &bimpY);
  114. #ifdef GZIP_SUPPORT
  115. if( gzeof(fInputFile) ) {
  116. #else
  117. if( feof(fInputFile) ) {
  118. #endif
  119. cout << "-I MpdDCMSMMGenerator : End of input file reached." << endl;
  120. const Bool_t ZeroSizeEvents=kFALSE;
  121. if( ZeroSizeEvents) {
  122. return kTRUE; // gives zero multiplicity events after EOF
  123. }
  124. else { // gives geant crash after EOF and one empty event in .root file
  125. #ifdef GZIP_SUPPORT
  126. gzclose(fInputFile);
  127. #else
  128. fclose(fInputFile);
  129. #endif
  130. fInputFile = NULL;
  131. return kFALSE;
  132. }
  133. }
  134. Float_t phi = atan2( bimpY, bimpX);
  135. // Set event id and impact parameter in MCEvent if not yet done
  136. FairMCEventHeader* event = primGen->GetEvent();
  137. if (event && (!event->IsSet())) {
  138. event->SetEventID(evnr);
  139. event->SetB(b);
  140. event->MarkSet(kTRUE);
  141. event->SetRotZ(phi);
  142. }
  143. TDatabasePDG* dbPDG = TDatabasePDG::Instance();
  144. Float_t px,py,pz;
  145. for( Int_t ibeam=0; ibeam<3; ibeam++) { // spectators pz+, spectators pz-, participants.
  146. Int_t np=0;
  147. #ifdef GZIP_SUPPORT
  148. ch= gzgets( fInputFile, read, 80);
  149. #else
  150. ch= fgets(read, 80, fInputFile);
  151. #endif
  152. sscanf(read, "%d", &np); //cout<<np<<" "<<endl;
  153. for( Int_t i=0; i<np; i++) {
  154. #ifdef GZIP_SUPPORT
  155. ch= gzgets( fInputFile, read, 80);
  156. #else
  157. ch= fgets(read, 80, fInputFile);
  158. #endif
  159. Int_t iN, iQ, iS=0;
  160. Float_t xxx=0.,mass;
  161. if( ibeam < 2) sscanf(read, "%d %d %f %f %f %f", &iN, &iQ, &xxx, &px,&py,&pz);
  162. else sscanf(read, "%d %d %d %f %f %f %f", &iN, &iQ, &iS, &px,&py,&pz, &mass);
  163. Int_t pdgID=0;
  164. if( ibeam==2) { // participants
  165. Double_t massFactor=1.;
  166. pdgID= FindPDGCodeParticipant( iN, iS, iQ, mass, massFactor);
  167. if( massFactor != 1.) { px*=massFactor; py*=massFactor; pz*=massFactor; }
  168. if( pdgID>1000030000) cout<<pdgID<<" "<<iN<<" "<<iS<<" "<<iQ<<" "<<mass<<endl;
  169. } else { // spectators
  170. if( fSpectatorsON) {
  171. Int_t dN=-999; // difference of number of nucleons iN-NTab between DCM-DCMSMM and registered ions
  172. pdgID= FindPDGCodeSpectator( iN, iQ, dN);
  173. }
  174. }
  175. if( pdgID) {
  176. if( fFixedTarget) {
  177. TParticlePDG* particle= dbPDG->GetParticle(pdgID);
  178. Double_t m= particle->Mass();
  179. Double_t e= pow( px*px+ py*py+ pz*pz+ m*m, 0.5 );
  180. Double_t pzF = fGammaCM * ( pz + fBetaCM * e); // Lorentz transformation to fixed target system
  181. Double_t eF= pow( px*px+ py*py+ pzF*pzF+ m*m, 0.5 );
  182. // cout<<fGammaCM<<" "<<fBetaCM<<" "<<pdgID<<" "<<m<<" "<<pz<<" "<<pzF<<" "<<eF-m<<endl;
  183. pz=pzF;
  184. }
  185. // Int_t Geant3ID= dbPDG->ConvertPdgToGeant3(pdgID);
  186. primGen->AddTrack( pdgID, px, py, pz, 0., 0., 0.);
  187. }
  188. else {
  189. if( ibeam==2 || (ibeam==2&&fSpectatorsON))
  190. cout<<"-I MpdDCMSMMGenerator : unknown particle N="<<iN<<" Q="<<iQ<<endl;
  191. }
  192. }
  193. }
  194. return kTRUE;
  195. }
  196. // ------------------------------------------------------------------------
  197. Bool_t MpdDCMSMMGenerator::SkipEvents(Int_t count) {
  198. if (count <= 0) return kTRUE;
  199. for (Int_t ii = 0; ii < count; ii++) {
  200. // ---> Check for input file
  201. if (!fInputFile) {
  202. cout << "-E MpdDCMSMMGenerator: Input file not open! " << endl;
  203. return kFALSE;
  204. }
  205. char read[80];
  206. #ifdef GZIP_SUPPORT
  207. char* ch= gzgets( fInputFile, read, 80);
  208. #else
  209. char* ch= fgets(read, 80, fInputFile);
  210. #endif
  211. Int_t evnr=0; Float_t b, bimpX, bimpY;
  212. sscanf(read, "%d %f %f %f", &evnr, &b, &bimpX, &bimpY);
  213. for( Int_t ibeam=0; ibeam<3; ibeam++) {
  214. Int_t np=0;
  215. #ifdef GZIP_SUPPORT
  216. ch= gzgets( fInputFile, read, 80);
  217. #else
  218. ch= fgets(read, 80, fInputFile);
  219. #endif
  220. sscanf(read, "%d", &np);
  221. for( Int_t i=0; i<np; i++) {
  222. #ifdef GZIP_SUPPORT
  223. ch= gzgets( fInputFile, read, 80);
  224. #else
  225. ch= fgets(read, 80, fInputFile);
  226. #endif
  227. Int_t iN, iQ, iS=0;
  228. Float_t xxx=0., mass, px,py,pz;
  229. if( ibeam < 2) sscanf(read, "%d %d %f %f %f %f", &iN, &iQ, &xxx, &px,&py,&pz); //cout<<np<<" "<<endl;
  230. else sscanf(read, "%d %d %d %f %f %f %f", &iN, &iQ, &iS, &px,&py,&pz, &mass);
  231. }
  232. }
  233. }
  234. return kTRUE;
  235. }
  236. // ------------------------------------------------------------------------
  237. Int_t MpdDCMSMMGenerator::FindPDGCodeParticipant( Int_t A, Int_t S, Int_t Z, Float_t mass, Double_t &massFactor) {
  238. Int_t k=0;
  239. Int_t aA= abs(A), aS= abs(S), aZ= abs(Z);
  240. if( aA == 0) { // mesons and gamma =====================================
  241. if( S == 0) { // not strange
  242. if( aZ == 0) { // neutral
  243. const Int_t n000=7; // gamma, pi0, eta, rho, omega, etaPri, phi
  244. Float_t M000M[n000]={ 0.0,0.13497,0.54745,0.7690,0.78265,0.95778,1.019455};
  245. Float_t M000L[n000]={ 0.0,0.13400,0.54700,0.7680,0.78200,0.95700,1.019390};
  246. Float_t M000U[n000]={ 0.0,0.13510,0.54810,0.7711,0.78310,0.95810,1.019610};
  247. Int_t C000[n000]= { 22, 111, 211, 113, 223, 331, 333};
  248. for( Int_t i=0; i<n000; i++) {
  249. if( mass < M000L[i] || mass > M000U[i]) continue;
  250. k= C000[i];
  251. break;
  252. }
  253. if( k==0 && mass>0.4979 && mass<0.4981 ) k=22; // rare particle in pairs in begining of list
  254. // invariant mass of these pairs (considering them as gamma) is close to mass of pion
  255. }
  256. else { // non-stranges charged mesons
  257. if( mass > 0.1389 && mass < 0.1401) k=211*Z; // pi+/-
  258. else if( mass > 0.134 && mass < 0.136) k=211*Z; // also pi+/-
  259. else if( mass == 0.776) k=213*Z; // rho+/-
  260. else if( mass > 0.0004 && mass < 0.0011) k=-11*Z; // e+/e-
  261. else if( mass > 0.1055 && mass < 0.1065) k=-13*Z; // mu+/mu-
  262. }
  263. }
  264. else { // A==0 && S!=0 strange mesons
  265. if( aZ == 0) { // neutral K0, K*0
  266. if( mass > 0.491 && mass < 0.4981) { // K0 S=1 311, AK0 S=-1 -311
  267. Double_t KSL= gRandom->Uniform(0.,1.);
  268. if( KSL<0.5) k= 130; // KL
  269. else k= 310; // KS
  270. } else if( mass == 0.8922) k=313*S; // K*0
  271. }
  272. if( aZ == 1) {
  273. if( mass > 0.491 && mass < 0.4951) k=321*Z; // K+
  274. else if( mass > 0.8815 && mass > 0.8895 ) k=323*Z; // K*+ (not occured)
  275. }
  276. }
  277. }
  278. else if ( aA==1 ) { // barions ============================================
  279. if( S == 0) { // Light baryons: n,p,Delta
  280. if( Z==0) {
  281. if( ( mass > 0.9389 && mass < 0.940) || mass < -1. ) k=2112*A; // n
  282. else if( mass == 1.232) k=2114*A; // Delta0 (not observed)
  283. }
  284. else if( aZ==1) {
  285. if( ( mass > 0.9379 && mass < 0.9401) || mass < -1. ) k=2212*A; // p
  286. else if( mass == 1.232) {
  287. if( A*Z == 1) k=2214*A; // Delta+ (not observed)
  288. else if( A*Z == -1) k=1114*A; // Delta- (not observed)
  289. }
  290. }
  291. else if( aZ==2) {
  292. if( mass == 1.232) k=2224*A; // Delta++ (not observed)
  293. }
  294. }
  295. else if( aS==1) { // Lambda, Sigma, Sigma*
  296. if( Z == 0) {
  297. if( mass > 1.1149 && mass < 1.1161) k= 3122*A; // Lambda
  298. else if( mass > 1.189 && mass < 1.1921 ) k= 3212*A; // Sigma0
  299. else if( mass == 1.382 ) k= 3214*A; // Sigma*0 (not observed)
  300. }
  301. else if( aZ == 1) { // Sigma-, Sigma+, Sigma*-, Sigma*+
  302. if( mass > 1.1889 && mass < 1.1921) k= 3222*A; // Sigma+
  303. else if( mass > 1.196 && mass < 1.1971) k= 3112*A; // Sigma-
  304. else if( mass == 1.3823) k= 3224*A; // Sigma*+ (not observed)
  305. else if( mass == 1.3875) k= 3114*A; // Sigma*- (not observed)
  306. }
  307. }
  308. else if( aS==2) { // Xi, Xi*
  309. if( Z == 0) { // Xi0, Xi*0
  310. if( mass > 1.3139 && mass < 1.3161 ) k= 3322*A; // Xi0
  311. else if( mass == 1.5318) k= 3324*A; // Xi*0
  312. }
  313. else if( aZ == 1) { // Xi-, Xi*-
  314. if( mass > 1.3139 && mass < 1.3221 ) k= 3312*A; // Xi-
  315. else if( mass == 1.535) k= 3314*A; // Xi*-
  316. }
  317. }
  318. else if( aS==3) { // Omega-
  319. if( aZ == 1) {
  320. if( mass > 1.6 && mass < 1.801 ) k= 3334*A; // Omega- 1.67245
  321. }
  322. }
  323. }
  324. else if ( A==2 ) { // d is only valid ion (particle). All over are unstable. ==
  325. if( S==0) {
  326. if( Z==1 && ( mass<-1 || ( mass>1.875 && mass<1.877))) k=1000010020; // d
  327. }
  328. else if( S==-1) { // unstable strange pairs => take the most heavy nucleon
  329. if( Z==-1) { k= 3112; massFactor=1.197436/mass;} // Sigma- n 2.13699 GeV (not occured)
  330. else if( Z==0 && mass>2.0545 &&mass<2.0555 ) { k= 3122; massFactor=1.11568/mass;} // Lambda n
  331. else if( Z==1 && mass>2.0535 &&mass<2.0545 ) { k= 3122; massFactor=1.11568/mass;} // Lambda p
  332. else if( Z==2 ) { k= 3222; massFactor=1.18937/mass;} // Sigma+ p 2.12764 GeV (not occured)
  333. }
  334. else if( S==-2) {
  335. if( Z==-1 && mass>2.2545 && mass<2.2615) { k= 3312; massFactor=1.32132/mass;} // Xi- n M=2.26088 ( also ghost at 2.255 GeV)
  336. else if( Z==0) {
  337. if( mass>2.2305 && mass<2.2315) { k= 3122; massFactor=1.11568/mass;} // Lambda Lambda M=2.23136
  338. else if( mass>2.2525 && mass<2.2555) { k= 3322; massFactor=1.19255/mass;} // Xi0 n M=2.25446
  339. else if( mass>2.2595 && mass<2.2605) { k= 3312; massFactor=1.197436/mass;} // Xi- p M=2.25959
  340. }
  341. else if( Z==1 && mass>2.2525 && mass<2.2535) { k= 3322; massFactor=1.3149/mass;} // Xi0 p M=2.25317 (132)
  342. }
  343. else if( S==-3) {
  344. if( Z==-1 && mass>2.4305 && mass<2.4375) { k= 3312; massFactor=1.197436/mass;} // Xi- Lambda M=2.43700 ( also ghost at 2.431 GeV) (37)
  345. else if( Z==0 && mass>2.4301 && mass<2.4315) { k=3322; massFactor=1.3149/mass;} // Xi0 Lambda M=2.43058 (36)
  346. }
  347. else if( S==-4 && Z==-1 && mass>2.6295 && mass<2.6365) { k=3312; massFactor=1.32132/mass;} // Xi- Xi0 M=2.63622 (1)
  348. }
  349. else if( A==3) {
  350. //
  351. if( S==0) {
  352. if( Z==1 && mass>2.8085 && mass<2.8095) k=1000010030; // t
  353. else if( Z==2 && mass>2.8085 && mass<2.8095) k=1000020030; // He3
  354. }
  355. else if( S==-1) {
  356. if( Z==1 && mass>2.9915 && mass<2.9925) k=1010010030; // H3L
  357. }
  358. }
  359. else if( A==4) {
  360. if( S==0) {
  361. if( Z==2 && mass>3.7275 && mass<3.7285) k=1000020040; // Alpha
  362. //else if( Z==1 && mass>3.7275 && mass<3.7285) k=1000020040; // H4
  363. }
  364. else if( S==-1) {
  365. if( Z==1 && mass>3.9245 && mass<3.9255) k=1010010040; // H4L
  366. if( Z==2 && mass>3.9245 && mass<3.9255) k=1010020040; // He4L
  367. }
  368. else if( S==-2) {
  369. if( Z==1 && mass>4.1065 && mass<4.1075 ) k=1020010040; // H4LL (6 per 100k MPD events)
  370. }
  371. }
  372. else if( A==5) {
  373. if( S==-2 && mass>5.0405 && mass<5.0415) k= 1020020050; // He5LL (0 per 100k MPD events)
  374. }
  375. //if( aA<=2 && k==0) {
  376. if( mass>-10 && k==0) {
  377. //if( A==0 && S==0 && Z==0 && fabs( mass)-0.498 < 0.001) {}
  378. //else
  379. cout<<"A="<<A<<" S="<<S<<" Z="<<Z<<" mass= "<<mass<<" Pdg="<<k<<endl;
  380. }
  381. //k= 1000030060;
  382. return k;
  383. }
  384. // ------------------------------------------------------------------------
  385. Int_t MpdDCMSMMGenerator::FindPDGCodeSpectator( Int_t N, Int_t Z, Int_t &dN) {
  386. Int_t k=0;
  387. dN=0;
  388. if( Z==0 && N==1) k=2112; // n
  389. else if( Z==1 && N==1 ) k=2212; // p
  390. else if( Z>=1 && Z<=fZMax ) {
  391. Int_t NTab;
  392. if( N<fN1[Z]) { NTab=fN1[Z]; dN=N-fN1[Z];}
  393. else if( N<=fN2[Z]) { NTab=N; dN=0;}
  394. else { NTab=fN2[Z]; dN=N-fN2[Z];}
  395. k=1000000000+ Z*10000+ NTab*10;
  396. }
  397. return k;
  398. }
  399. // ------------------------------------------------------------------------
  400. Int_t MpdDCMSMMGenerator::RegisterIons( void) {
  401. // see ./fairsoft/transport/geant3/TGeant3/TGeant3.cxx::AddParticlesToPdgDataBase()
  402. // or ./fairsoft/transport/geant4_vmc/source/physics/src/TG4ParticlesManager::DefineParticles()
  403. // Deuteron, Triton, Alpha, HE3, HyperTriton and their ANTI-particle are added there
  404. // also Special particles: "Cherenkov", "FeedbackPhoton".
  405. // He4L, H4L are added in mpdroot/gconfig/UserDecay.C
  406. struct gpions { Int_t N,Z; Char_t Name[6]; Double_t Mass; };
  407. const Int_t NSI=96;
  408. const struct gpions ions[NSI] = {
  409. { 2, 1, "d", 0.0}, { 3, 1, "t", 0.0},
  410. { 3, 2, "He3", 0.0}, { 4, 2, "He4", 0.0}, { 5, 2, "He5q", 4.6688534}, { 6, 2, "He6", 5.6065596},
  411. { 6, 3, "Li6", 5.60305}, {7, 3, "Li7", 6.53536}, { 8, 3, "Li8", 7.4728996}, { 9, 3, "Li9", 8.40840118},
  412. { 7, 4, "Be7", 6.53622}, { 8, 4, "Be8q", 7.4568945}, { 9, 4, "Be9", 8.39479},
  413. { 10, 5, "B10", 9.32699}, { 11, 5, "B11", 10.25510}, { 10, 6, "C10", 9.3306397},
  414. { 11, 6, "C11", 10.257085}, { 12, 6, "C12", 11.17793}, { 13, 6, "C13", 12.112548},
  415. { 14, 6, "C14", 13.043937}, { 14, 7, "N14", 13.04378}, { 16, 8, "O16", 14.89917},
  416. { 19, 9, "F19", 17.69690}, { 20, 10, "Ne20", 18.62284}, { 23, 11, "Na23", 21.41483},
  417. { 24, 12, "Mg24", 22.34193}, { 27, 13, "Al27", 25.13314}, { 28, 14, "Si28", 26.06034},
  418. { 31, 15, "P31", 28.85188}, { 32, 16, "S32", 29.78180}, { 35, 17, "Cl35", 32.57328},
  419. { 36, 18, "Ar36", 33.50356}, { 39, 19, "K39", 36.29447}, { 40, 20, "Ca40", 37.22492},
  420. { 45, 21, "Sc45", 41.87617}, { 48, 22, "Ti48", 44.66324}, { 51, 23, "V51", 47.45401},
  421. { 52, 24, "Cr52", 48.38228}, { 55, 25, "Mn55", 51.17447}, { 56, 26, "Fe56", 52.10307},
  422. { 59, 27, "Co59", 54.89593}, { 58, 28, "Ni58", 53.96644}, { 63, 29, "Cu63", 58.61856},
  423. { 64, 30, "Zn64", 59.54963}, { 69, 31, "Ga69", 64.2037653}, { 74, 32, "Ge74", 68.85715},
  424. { 75, 33, "As75", 69.78902528}, { 80, 34, "Se80", 74.44178}, { 81, 35, "Br81", 75.373047},
  425. { 84, 36, "Kr84", 78.16309}, { 85, 37, "Rb85", 79.09483}, { 88, 38, "Sr88", 81.88358},
  426. { 89, 39, "Y89", 82.81527}, { 90, 40, "Zr90", 83.74571}, { 93, 41, "Nb93", 86.541743},
  427. { 98, 42, "Mo98", 91.19832}, { 99, 43, "Tc99", 92.13059}, { 100, 44, "Ru100", 93.060191},
  428. { 103, 45, "Rh103", 93.9934966}, { 106, 46, "Pd106", 98.64997}, { 109, 47, "Ag109", 101.444134},
  429. { 114, 48, "Cd114", 106.10997}, { 117, 49, "In117", 108.89586}, { 120, 50, "Sn120", 111.68821},
  430. { 123, 51, "Sb123", 114.48455}, { 125, 52, "Te125", 116.3477}, { 127, 53, "I127", 118.210768},
  431. { 132, 54, "Xe132", 122.86796}, { 135, 55, "Cs135", 125.66412195}, { 138, 56, "Ba138", 128.45793},
  432. { 139, 57, "La139", 129.3904488}, { 140, 58, "Ce140", 130.32111}, { 141, 59, "Pr141", 131.2546475},
  433. { 145, 60, "Nd145", 134.9852}, { 149, 61, "Pm149", 138.716549}, { 152, 62, "Sm152", 141.51236},
  434. { 153, 63, "Eu153", 142.44522}, { 157, 64, "Gd157", 146.17374}, { 161, 65, "Tb161", 149.90308},
  435. { 164, 66, "Dy164", 152.69909}, { 165, 67, "Ho165", 153.63162}, { 168, 68, "Er168", 156.42801},
  436. { 171, 69, "Tm171", 159.2262758}, { 174, 70, "Yb174", 162.02245}, { 175, 71, "Lu175", 162.956297},
  437. { 178, 72, "Hf178", 165.7535059}, { 181, 73, "Ta181", 168.55199}, { 184, 74, "W184", 171.34924},
  438. { 185, 75, "Re185", 172.28258}, { 188, 76, "Os188", 175.079754}, { 191, 77, "Ir191", 177.878667},
  439. { 194, 78, "Pt194", 180.67513}, { 197, 79, "Au197", 183.47324}, { 202, 80, "Hg202", 188.13451},
  440. { 205, 81, "Tl205", 190.93247}, { 208, 82, "Pb208", 193.72907}};
  441. TDatabasePDG* dbPDG = TDatabasePDG::Instance();
  442. for( Int_t i=0; i<NSI; i++) {
  443. if( ions[i].Mass < 0.1) continue; // expected that it is already registered
  444. Int_t ID= 1000000000+ ions[i].Z*10000+ ions[i].N*10;
  445. if ( dbPDG->GetParticle(ID) ) continue;
  446. FairIon* ion = new FairIon( ions[i].Name, ions[i].Z, ions[i].N, ions[i].Z) ;
  447. FairRunSim::Instance()->AddNewIon(ion);
  448. }
  449. for( Int_t i=0; i<NSI; i++) {
  450. Int_t z= ions[i].Z; if( z<0 && z>fZMax) continue;
  451. fN2[z]=ions[i].N;
  452. }
  453. for( Int_t i=NSI-1; i>=0; i--) {
  454. Int_t z= ions[i].Z; if( z<0 && z>fZMax) continue;
  455. fN1[z]=ions[i].N;
  456. }
  457. return NSI;
  458. }
  459. ClassImp(MpdDCMSMMGenerator);