BmnRoot
Loading...
Searching...
No Matches
MpdDCMSMMGenerator.cxx
Go to the documentation of this file.
1// -------------------------------------------------------------------------
2// ----- MpdDCMSMMGenerator source file -----
3// ----- Created 27-AUG-2019 by Igor Rufanov -----
4// -------------------------------------------------------------------------
6
7#include "FairIon.h"
8#include "FairMCEventHeader.h"
9#include "FairPrimaryGenerator.h"
10#include "FairRunSim.h"
11#include "TDatabasePDG.h"
12#include "TParticlePDG.h"
13#include "TRandom.h"
14
15using namespace std;
16using namespace TMath;
17
18// ----- Default constructor ------------------------------------------
20 : FairGenerator()
21 , fInputFile(nullptr)
22 , fFileName("")
23 , fBoostType("CmsFixed")
24 , fUseLorentzBoost(kFALSE)
25 , fBoostBeta(0.0)
26 , fBoostGamma(0.0)
27 , fSpectatorsON(kTRUE)
28 , fInputFileVersion(1)
29{}
30// ------------------------------------------------------------------------
31
32// ----- Standard constructor -----------------------------------------
34 : FairGenerator()
35 , fInputFile(nullptr)
36 , fFileName(fileName)
37 , fBoostType("CmsFixed")
38 , fUseLorentzBoost(kFALSE)
39 , fBoostBeta(0.0)
40 , fBoostGamma(0.0)
41 , fSpectatorsON(kTRUE)
42 , fInputFileVersion(1)
43{
44 cout << "-I MpdDCMSMMGenerator: Opening input file " << fFileName.Data() << endl;
45#ifdef GZIP_SUPPORT
46 fInputFile = gzopen(fFileName.Data(), "rb");
47#else
48 fInputFile = fopen(fFileName.Data(), "r");
49#endif
50 if (!fInputFile) {
51 Fatal("MpdDCMSMMGenerator", "Cannot open input file.");
52 exit(1);
53 }
54 // read file header
55 char read[200];
56 Int_t A1, Z1, A2, Z2;
57 Double_t T0, sqS;
58 Int_t header_size = 1; // old and current formats have different header size
59
60 for (Int_t i = 0; i < header_size; i++) {
61#ifdef GZIP_SUPPORT
62 /*char* ch = */ gzgets(fInputFile, read, 200);
63#else
64 /*char* ch = */ fgets(read, 200, fInputFile);
65#endif
66 cout << "-I MpdDCMSMMGenerator:" << read;
67 if (i == 0) {
68 string str0(read);
69
70 if (str0.find("Results of DCM-SMM") != string::npos) { // new version
71 fInputFileVersion = 1;
72 header_size = 20;
73 } else if (str0.find("Results of QGSM") != string::npos) { // old version
74 fInputFileVersion = 0;
75 header_size = 3;
76 } else {
77 Fatal("MpdDCMSMMGenerator", "Wrong input file format.");
78 exit(1);
79 }
80
81 A1 = stoi(str0.substr(str0.find("A1=") + 3, 3));
82 Z1 = stoi(str0.substr(str0.find("Z1=") + 3, 3));
83 A2 = stoi(str0.substr(str0.find("A2=") + 3, 3));
84 Z2 = stoi(str0.substr(str0.find("Z2=") + 3, 3));
85 // Int_t A1= 0;
86 } else if (i == 1) {
87 string str0(read); // hypcoa-b1.f line 711: ss=sqrt(1.88*(T0+1.88)) => mp=0.940
88 T0 = stof(str0.substr(str0.find("T0=") + 3, 8));
89 sqS = stof(str0.substr(str0.find("sqrt(s)=") + 8, 8));
90
91 Double_t mp = 0.940; // to obtain equivalence of read "sqrt(s)=" and calculated below mCMS
92 Double_t e = T0 + mp;
93 Double_t p = pow(e * e - mp * mp, 0.5);
94 Double_t mCMS = pow(2. * mp * (e + mp), 0.5);
95 cout << "mCMS=" << mCMS << endl;
96
97 if (fInputFileVersion == 0) {
98 fUseLorentzBoost = kTRUE;
99 if (fBoostType == "None")
100 fUseLorentzBoost = kFALSE;
101 else if (fBoostType == "CmsFixed") { // from CMS to target (2-nd particle) rest system
102 fBoostGamma = (e + mp) / mCMS;
103 fBoostBeta = pow(1. - 1. / (fBoostGamma * fBoostGamma), 0.5);
104 } else if (fBoostType == "CmsFixedInverted") { // from CMS to beam ( 1-st particle) rest system
105 fBoostGamma = -(e + mp) / mCMS;
106 fBoostBeta = -pow(1. - 1. / (fBoostGamma * fBoostGamma), 0.5);
107 } else if (fBoostType == "FixedFixedInverted") { // from target rest system to beam RS
108 fBoostGamma = -e / mp;
109 fBoostBeta = -p / e;
110 }
111 }
112 }
113 }
114 cout << "-I MpdDCMSMMGenerator: A1=" << A1 << " Z1=" << Z1 << " A2=" << A2 << " Z2=" << Z2 << " T0=" << T0
115 << " sqS=" << sqS << endl;
116
117 if (fSpectatorsON) /*Int_t n = */
118 RegisterIons();
119}
120// ------------------------------------------------------------------------
121
122// ----- Destructor ---------------------------------------------------
124{
125 if (fInputFile) {
126#ifdef GZIP_SUPPORT
127 gzclose(fInputFile);
128#else
129 fclose(fInputFile);
130#endif
131 fInputFile = NULL;
132 }
133 // cout<<"Leave Destructor of MpdDCMSMMGenerator"<<endl;
134}
135// ------------------------------------------------------------------------
136
137// ----- Public method ReadEvent --------------------------------------
138Bool_t MpdDCMSMMGenerator::ReadEvent(FairPrimaryGenerator* primGen)
139{
140 // ---> Check for input file
141 // cout<<"MpdDCMSMMGenerator::ReadEvent -------------------------"<<endl;
142 if (!fInputFile) {
143 cout << "-E MpdDCMSMMGenerator: Input file not open! " << endl;
144 return kFALSE;
145 }
146
147 // ---> Check for primary generator
148 if (!primGen) {
149 cout << "-E- MpdDCMSMMGenerator::ReadEvent: "
150 << "No PrimaryGenerator!" << endl;
151 return kFALSE;
152 }
153
154 char read[128]; // pnaleks: It was 80, but the line length of 98 chars has been seen
155#ifdef GZIP_SUPPORT
156 /*char* ch = */ gzgets(fInputFile, read, 128);
157#else
158 /*char* ch = */ fgets(read, 128, fInputFile);
159#endif
160 Int_t evnr = 0;
161 Float_t b, bimpX, bimpY;
162 Int_t npp = 0; // pnaleks: number of produced particles after cascade and light clusters after coalescence stages
163 // (for fInputFileVersion == 1)
164
165 if (fInputFileVersion == 0) // old version
166 sscanf(read, "%d %f %f %f", &evnr, &b, &bimpX, &bimpY);
167 else // new version
168 sscanf(read, "%d %d %f %f %f", &evnr, &npp, &b, &bimpX, &bimpY);
169
170#ifdef GZIP_SUPPORT
171 if (gzeof(fInputFile)) {
172#else
173 if (feof(fInputFile)) {
174#endif
175 cout << "-I MpdDCMSMMGenerator : End of input file reached." << endl;
176 const Bool_t ZeroSizeEvents = kFALSE;
177 if (ZeroSizeEvents)
178 return kTRUE; // gives zero multiplicity events after EOF
179 else { // gives geant crash after EOF and one empty event in .root file
180#ifdef GZIP_SUPPORT
181 gzclose(fInputFile);
182#else
183 fclose(fInputFile);
184#endif
185 fInputFile = NULL;
186 return kFALSE;
187 }
188 }
189
190 Float_t phi = atan2(bimpY, bimpX);
191
192 // Set event id and impact parameter in MCEvent if not yet done
193 FairMCEventHeader* event = primGen->GetEvent();
194 if (event && (!event->IsSet())) {
195 event->SetEventID(evnr);
196 event->SetB(b);
197 event->MarkSet(kTRUE);
198 event->SetRotZ(phi);
199 }
200
201 TDatabasePDG* dbPDG = TDatabasePDG::Instance();
202
203 Float_t px, py, pz;
204 for (Int_t ibeam = 0; ibeam < 3; ibeam++) { // spectators pz+, spectators pz-, participants.
205 Int_t np = 0;
206#ifdef GZIP_SUPPORT
207 /*ch = */ gzgets(fInputFile, read, 128);
208#else
209 /*ch = */ fgets(read, 128, fInputFile);
210#endif
211 sscanf(read, "%d", &np); // cout<<np<<" "<<endl;
212 // Information for fInputFileVersion == 1:
213 // pnaleks: In this line is also an additional information for ibeam < 2,
214 // which is ignored because next lines contains this datta by fragments
215
216 for (Int_t i = 0; i < np; i++) {
217#ifdef GZIP_SUPPORT
218 /*ch = */ gzgets(fInputFile, read, 128);
219#else
220 /*ch = */ fgets(read, 128, fInputFile);
221#endif
222 Int_t iN, iQ, iS = 0;
223 Float_t xxx = 0., mass;
224 Int_t pdgID_in = 0, iL;
225 Float_t pLabZ /*, pALabZ = 0*/;
226 if (fInputFileVersion == 0) {
227 if (ibeam < 2)
228 sscanf(read, "%d %d %f %f %f %f", &iN, &iQ, &xxx, &px, &py, &pz);
229 else
230 sscanf(read, "%d %d %d %f %f %f %f", &iN, &iQ, &iS, &px, &py, &pz, &mass);
231 } else { // fInputFileVersion == 1
232 Float_t var10, var11;
233 Int_t res = sscanf(read, "%d%d%d%d%d%f%f%f%f%f%f", &iQ, &iL, &iS, &iN, &pdgID_in, &px, &py, &pz, &pLabZ,
234 &var10, &var11);
235 switch (res) {
236 case 10:
237 mass = var10;
238 break;
239 case 11:
240 // pALabZ = var10;
241 mass = var11;
242 break;
243 default:
244 Fatal(__func__, ": data format mismatch, event %d\n", evnr);
245 }
246 } // if (fInputFileVersion == 0)
247
248 // Int_t pdgID=0;
249 Int_t pdgID = pdgID_in;
250 if (fInputFileVersion == 0) {
251 if (ibeam == 2) { // participants
252 Double_t massFactor = 1.;
253 pdgID = FindPDGCodeParticipant(iN, iS, iQ, mass, massFactor);
254 if (massFactor != 1.) {
255 px *= massFactor;
256 py *= massFactor;
257 pz *= massFactor;
258 }
259 // if (pdgID > 1000030000) cout<<pdgID<<" "<<iN<<" "<<iS<<" "<<iQ<<" "<<mass<<endl;
260 } else { // spectators
261 if (fSpectatorsON) {
262 Int_t dN =
263 -999; // difference of number of nucleons iN-NTab between DCM-DCMSMM and registered ions
264 pdgID = FindPDGCodeSpectator(iN, iQ, dN);
265 }
266 }
267 } else { // AZ-250324 - fix error with K0 and K0bar written as K0S and K0L and handle K0 and K0bar
268 if (pdgID == 310 || pdgID == 130 || TMath::Abs(pdgID) == 311) {
269 Double_t kSL = gRandom->Uniform(0., 1.);
270 if (kSL < 0.5)
271 pdgID = 130; // KL
272 else
273 pdgID = 310; // KS
274 }
275 } // AZ
276
277 TParticlePDG* particle = pdgID ? dbPDG->GetParticle(pdgID) : 0;
278 if (particle) {
279 if (fInputFileVersion == 0) {
280 if (fUseLorentzBoost) {
281 // TParticlePDG* particle= dbPDG->GetParticle(pdgID);
282 Double_t m = particle->Mass();
283 Double_t e = pow(px * px + py * py + pz * pz + m * m, 0.5);
284 Double_t pzF = fBoostGamma * (pz + fBoostBeta * e);
285 // Double_t eF = pow( px*px+ py*py+ pzF*pzF+ m*m, 0.5 );
286 // cout<<fBoostGamma<<" "<<fBoostBeta<<" "<<pdgID<<" "<<m<<" "<<pz<<" "<<pzF<<" "<<eF-m<<endl;
287 pz = pzF;
288 }
289 } else // fInputFileVersion == 1
290 pz = pLabZ;
291
292 // Int_t Geant3ID= dbPDG->ConvertPdgToGeant3(pdgID);
293 if (fabs(pz) > 50.)
294 cout << "pz=" << pz << " N=" << iN << " Q=" << iQ << " pdg=" << pdgID << "\n";
295 primGen->AddTrack(pdgID, px, py, pz, 0., 0., 0.);
296 } else {
297 // if( ibeam==2 || (ibeam==2&&fSpectatorsON))
298 if (fInputFileVersion == 1 || ibeam == 2 || fSpectatorsON) {
299 // cout<<"-I MpdDCMSMMGenerator : unknown particle N="<<iN<<" Q="<<iQ<<endl;
300 cout << "-I MpdDCMSMMGenerator : unknown particle N=" << iN << " Q=" << iQ;
301 if (fInputFileVersion == 1)
302 cout << " PDG_in=" << pdgID_in << " ibeam=" << ibeam << endl;
303 else
304 cout << endl;
305 }
306 }
307 }
308 }
309 return kTRUE;
310}
311// ------------------------------------------------------------------------
312
314{
315 if (count <= 0)
316 return kTRUE;
317 for (Int_t ii = 0; ii < count; ii++) {
318 // ---> Check for input file
319 if (!fInputFile) {
320 cout << "-E MpdDCMSMMGenerator: Input file not open! " << endl;
321 return kFALSE;
322 }
323 char read[128];
324#ifdef GZIP_SUPPORT
325 /*char* ch = */ gzgets(fInputFile, read, 128);
326#else
327 /*char* ch = */ fgets(read, 128, fInputFile);
328#endif
329 // Int_t evnr=0; Float_t b, bimpX, bimpY;
330 // sscanf(read, "%d %f %f %f", &evnr, &b, &bimpX, &bimpY);
331 for (Int_t ibeam = 0; ibeam < 3; ibeam++) {
332 Int_t np = 0;
333#ifdef GZIP_SUPPORT
334 /*ch = */ gzgets(fInputFile, read, 128);
335#else
336 /*ch = */ fgets(read, 128, fInputFile);
337#endif
338 sscanf(read, "%d", &np);
339 for (Int_t i = 0; i < np; i++) {
340#ifdef GZIP_SUPPORT
341 /*ch = */ gzgets(fInputFile, read, 128);
342#else
343 /*ch = */ fgets(read, 128, fInputFile);
344#endif
345 // Int_t iN, iQ, iS=0;
346 // Float_t xxx=0., mass, px,py,pz;
347 // if( ibeam < 2) sscanf(read, "%d %d %f %f %f %f", &iN, &iQ, &xxx, &px,&py,&pz); //cout<<np<<" "<<endl;
348 // else sscanf(read, "%d %d %d %f %f %f %f", &iN, &iQ, &iS, &px,&py,&pz, &mass);
349 }
350 }
351 }
352
353 return kTRUE;
354}
355// ------------------------------------------------------------------------
356
357Int_t MpdDCMSMMGenerator::FindPDGCodeParticipant(Int_t A, Int_t S, Int_t Z, Float_t mass, Double_t& massFactor)
358{
359 Int_t k = 0;
360 Int_t aA = abs(A), aS = abs(S), aZ = abs(Z);
361
362 if (aA == 0) { // mesons and gamma =====================================
363 if (S == 0) { // not strange
364 if (aZ == 0) { // neutral
365 const Int_t n000 = 7; // gamma, pi0, eta, rho, omega, etaPri, phi
366 // Float_t M000M[n000]={ 0.0,0.13497,0.54745,0.7690,0.78265,0.95778,1.019455};
367 Float_t M000L[n000] = {0.0, 0.13400, 0.54700, 0.7680, 0.78200, 0.95700, 1.019390};
368 Float_t M000U[n000] = {0.0, 0.13510, 0.54810, 0.7711, 0.78310, 0.95810, 1.019610};
369 Int_t C000[n000] = {22, 111, 211, 113, 223, 331, 333};
370 for (Int_t i = 0; i < n000; i++) {
371 if (mass < M000L[i] || mass > M000U[i])
372 continue;
373 k = C000[i];
374 break;
375 }
376
377 if (k == 0 && mass > 0.4979 && mass < 0.4981)
378 k = 22; // rare particle in pairs in begining of list
379 // invariant mass of these pairs (considering them as gamma) is close to mass of pion
380 } else { // non-stranges charged mesons
381 if (mass > 0.1389 && mass < 0.1401)
382 k = 211 * Z; // pi+/-
383 else if (mass > 0.134 && mass < 0.136)
384 k = 211 * Z; // also pi+/-
385 else if (mass == 0.776)
386 k = 213 * Z; // rho+/-
387 else if (mass > 0.0004 && mass < 0.0011)
388 k = -11 * Z; // e+/e-
389 else if (mass > 0.1055 && mass < 0.1065)
390 k = -13 * Z; // mu+/mu-
391 }
392 } else { // A==0 && S!=0 strange mesons
393 if (aZ == 0) { // neutral K0, K*0
394 if (mass > 0.491 && mass < 0.4981) { // K0 S=1 311, AK0 S=-1 -311
395 Double_t KSL = gRandom->Uniform(0., 1.);
396 if (KSL < 0.5)
397 k = 130; // KL
398 else
399 k = 310; // KS
400 } else if (mass == 0.8922)
401 k = 313 * S; // K*0
402 }
403 if (aZ == 1) {
404 if (mass > 0.491 && mass < 0.4951)
405 k = 321 * Z; // K+
406 else if (mass > 0.8815 && mass > 0.8895)
407 k = 323 * Z; // K*+ (not occured)
408 }
409 }
410 } else if (aA == 1) { // barions ============================================
411 if (S == 0) { // Light baryons: n,p,Delta
412 if (Z == 0) {
413 if ((mass > 0.9389 && mass < 0.940) || mass < -1.)
414 k = 2112 * A; // n
415 else if (mass == 1.232)
416 k = 2114 * A; // Delta0 (not observed)
417 } else if (aZ == 1) {
418 if ((mass > 0.9379 && mass < 0.9401) || mass < -1.)
419 k = 2212 * A; // p
420 else if (mass == 1.232) {
421 if (A * Z == 1)
422 k = 2214 * A; // Delta+ (not observed)
423 else if (A * Z == -1)
424 k = 1114 * A; // Delta- (not observed)
425 }
426 } else if (aZ == 2) {
427 if (mass == 1.232)
428 k = 2224 * A; // Delta++ (not observed)
429 }
430 } else if (aS == 1) { // Lambda, Sigma, Sigma*
431 if (Z == 0) {
432 if (mass > 1.1149 && mass < 1.1161)
433 k = 3122 * A; // Lambda
434 else if (mass > 1.189 && mass < 1.1921)
435 k = 3212 * A; // Sigma0
436 else if (mass == 1.382)
437 k = 3214 * A; // Sigma*0 (not observed)
438 } else if (aZ == 1) { // Sigma-, Sigma+, Sigma*-, Sigma*+
439 if (mass > 1.1889 && mass < 1.1921)
440 k = 3222 * A; // Sigma+
441 else if (mass > 1.196 && mass < 1.1971)
442 k = 3112 * A; // Sigma-
443 else if (mass == 1.3823)
444 k = 3224 * A; // Sigma*+ (not observed)
445 else if (mass == 1.3875)
446 k = 3114 * A; // Sigma*- (not observed)
447 }
448 } else if (aS == 2) { // Xi, Xi*
449 if (Z == 0) { // Xi0, Xi*0
450 if (mass > 1.3139 && mass < 1.3161)
451 k = 3322 * A; // Xi0
452 else if (mass == 1.5318)
453 k = 3324 * A; // Xi*0
454 } else if (aZ == 1) { // Xi-, Xi*-
455 if (mass > 1.3139 && mass < 1.3221)
456 k = 3312 * A; // Xi-
457 else if (mass == 1.535)
458 k = 3314 * A; // Xi*-
459 }
460 } else if (aS == 3) { // Omega-
461 if (aZ == 1) {
462 if (mass > 1.6 && mass < 1.801)
463 k = 3334 * A; // Omega- 1.67245
464 }
465 }
466 } else if (A == 2) { // d is only valid ion (particle). All over are unstable. ==
467
468 if (S == 0) {
469 if (Z == 1 && (mass < -1 || (mass > 1.875 && mass < 1.877)))
470 k = 1000010020; // d
471 } else if (S == -1) { // unstable strange pairs => take the most heavy nucleon
472 if (Z == -1) {
473 k = 3112;
474 massFactor = 1.197436 / mass;
475 } // Sigma- n 2.13699 GeV (not occured)
476 else if (Z == 0 && mass > 2.0545 && mass < 2.0555)
477 {
478 k = 3122;
479 massFactor = 1.11568 / mass;
480 } // Lambda n
481 else if (Z == 1 && mass > 2.0535 && mass < 2.0545)
482 {
483 k = 3122;
484 massFactor = 1.11568 / mass;
485 } // Lambda p
486 else if (Z == 2)
487 {
488 k = 3222;
489 massFactor = 1.18937 / mass;
490 } // Sigma+ p 2.12764 GeV (not occured)
491 } else if (S == -2) {
492 if (Z == -1 && mass > 2.2545 && mass < 2.2615) {
493 k = 3312;
494 massFactor = 1.32132 / mass;
495 } // Xi- n M=2.26088 ( also ghost at 2.255 GeV)
496 else if (Z == 0)
497 {
498 if (mass > 2.2305 && mass < 2.2315) {
499 k = 3122;
500 massFactor = 1.11568 / mass;
501 } // Lambda Lambda M=2.23136
502 else if (mass > 2.2525 && mass < 2.2555)
503 {
504 k = 3322;
505 massFactor = 1.19255 / mass;
506 } // Xi0 n M=2.25446
507 else if (mass > 2.2595 && mass < 2.2605)
508 {
509 k = 3312;
510 massFactor = 1.197436 / mass;
511 } // Xi- p M=2.25959
512 } else if (Z == 1 && mass > 2.2525 && mass < 2.2535) {
513 k = 3322;
514 massFactor = 1.3149 / mass;
515 } // Xi0 p M=2.25317 (132)
516 } else if (S == -3) {
517 if (Z == -1 && mass > 2.4305 && mass < 2.4375) {
518 k = 3312;
519 massFactor = 1.197436 / mass;
520 } // Xi- Lambda M=2.43700 ( also ghost at 2.431 GeV) (37)
521 else if (Z == 0 && mass > 2.4301 && mass < 2.4315)
522 {
523 k = 3322;
524 massFactor = 1.3149 / mass;
525 } // Xi0 Lambda M=2.43058 (36)
526 } else if (S == -4 && Z == -1 && mass > 2.6295 && mass < 2.6365) {
527 k = 3312;
528 massFactor = 1.32132 / mass;
529 } // Xi- Xi0 M=2.63622 (1)
530 } else if (A == 3) {
531 //
532 if (S == 0) {
533 if (Z == 1 && mass > 2.8085 && mass < 2.8095)
534 k = 1000010030; // t
535 else if (Z == 2 && mass > 2.8085 && mass < 2.8095)
536 k = 1000020030; // He3
537 } else if (S == -1) {
538 if (Z == 1 && mass > 2.9915 && mass < 2.9925)
539 k = 1010010030; // H3L
540 }
541 } else if (A == 4) {
542 if (S == 0) {
543 if (Z == 2 && mass > 3.7275 && mass < 3.7285)
544 k = 1000020040; // Alpha
545 // else if( Z==1 && mass>3.7275 && mass<3.7285) k=1000020040; // H4
546 } else if (S == -1) {
547 if (Z == 1 && mass > 3.9245 && mass < 3.9255)
548 k = 1010010040; // H4L
549 if (Z == 2 && mass > 3.9245 && mass < 3.9255)
550 k = 1010020040; // He4L
551 } else if (S == -2) {
552 // if( Z==1 && mass>4.1065 && mass<4.1075 ) k=1020010040; // H4LL (6 per 100k MPD events)
553 }
554 } else if (A == 5) {
555 // if( S==-2 && mass>5.0405 && mass<5.0415) k= 1020020050; // He5LL (0 per 100k MPD events)
556 }
557
558 // if( aA<=2 && k==0) {
559 if (mass > -10 && k == 0) {
560 // if( A==0 && S==0 && Z==0 && fabs( mass)-0.498 < 0.001) {}
561 // else
562 cout << "A=" << A << " S=" << S << " Z=" << Z << " mass= " << mass << " Pdg=" << k << endl;
563 }
564 // k= 1000030060;
565 return k;
566}
567// ------------------------------------------------------------------------
568
569Int_t MpdDCMSMMGenerator::FindPDGCodeSpectator(Int_t N, Int_t Z, Int_t& dN)
570{
571 Int_t k = 0;
572 dN = 0;
573 if (Z == 0 && N == 1)
574 k = 2112; // n
575 else if (Z == 1 && N == 1)
576 k = 2212; // p
577 // else if( Z>=1 && Z<=fZMax ) { // find registeried ion with a same Z
578 // Int_t NTab;
579 // if( N<fN1[Z]) { NTab=fN1[Z]; dN=N-fN1[Z];}
580 // else if( N<=fN2[Z]) { NTab=N; dN=0;} // same B
581 // else { NTab=fN2[Z]; dN=N-fN2[Z];}
582 // k=1000000000+ Z*10000+ NTab*10;
583 // }
584 else if (N >= 2 && N <= fBMax) { // find registeried ion with a same baryon number
585 Int_t ZTab;
586 if (Z < fZ1[N]) {
587 ZTab = fZ1[N];
588 dN = Z - fZ1[N];
589 } else if (Z <= fZ2[N]) {
590 ZTab = Z;
591 dN = 0;
592 } // same Z
593 else
594 {
595 ZTab = fZ2[N];
596 dN = Z - fZ2[N];
597 }
598 k = 1000000000 + ZTab * 10000 + N * 10;
599 }
600 return k;
601}
602// ------------------------------------------------------------------------
603
605{
606 // see ./fairsoft/transport/geant3/TGeant3/TGeant3.cxx::AddParticlesToPdgDataBase()
607 // or ./fairsoft/transport/geant4_vmc/source/physics/src/TG4ParticlesManager::DefineParticles()
608 // Deuteron, Triton, Alpha, HE3, HyperTriton and their ANTI-particle are added there
609 // also Special particles: "Cherenkov", "FeedbackPhoton".
610 // He4L, H4L are added in mpdroot/gconfig/UserDecay.C
611
612 struct gpions
613 {
614 Int_t N, Z;
615 Char_t Name[6];
616 Double_t Mass;
617 };
618 const Int_t NSI = 236;
619 const struct gpions ions[NSI] = {
620 {2, 1, "d", 0.0},
621 {3, 1, "t", 0.0},
622 {3, 2, "He3", 0.0},
623 {4, 2, "He4", 0.0},
624 {5, 2, "He5q", 4.6688534},
625 {6, 2, "He6", 5.6065596},
626 {7, 2, "He7", 6.5465601},
627 {6, 3, "Li6", 5.60305},
628 {7, 3, "Li7", 6.53536},
629 {8, 3, "Li8", 7.4728996},
630 {9, 3, "Li9", 8.40840118},
631 {7, 4, "Be7", 6.53622},
632 {8, 4, "Be8q", 7.4568945},
633 {9, 4, "Be9", 8.39479},
634 {10, 4, "Be10", 9.32754},
635 {11, 4, "Be11", 10.2666},
636 {12, 4, "Be12", 11.203006},
637 {8, 5, "B8", 7.4748744},
638 {9, 5, "B9", 8.3958634},
639 {10, 5, "B10", 9.32699},
640 {11, 5, "B11", 10.25510},
641 {12, 5, "B12", 11.1913},
642 {9, 6, "C9", 8.4123574},
643 {10, 6, "C10", 9.3306397},
644 {11, 6, "C11", 10.257085},
645 {12, 6, "C12", 11.17793},
646 {13, 6, "C13", 12.112548},
647 {14, 6, "C14", 13.043937},
648 {15, 6, "C15", 13.98228},
649 {16, 6, "C16", 14.917600},
650 {11, 7, "N11", 10.257085},
651 {12, 7, "N12", 11.195267},
652 {13, 7, "N13", 12.11477},
653 {14, 7, "N14", 13.04378},
654 {15, 7, "N15", 13.972513},
655 {16, 7, "N16", 14.90959},
656 {17, 7, "N17", 15.843271},
657 // N11, o11-o13 are only unstable isotopes in this list (add for SRC analysis)
658 {11, 8, "O11", 10.257085},
659 {12, 8, "O12", 11.195267},
660 {13, 8, "O13", 12.132536},
661 {14, 8, "O14", 13.048925},
662 {15, 8, "O15", 13.97527},
663 {16, 8, "O16", 14.89917},
664 {17, 8, "O17", 15.834591},
665 {18, 8, "O18", 16.76611},
666 {19, 8, "O19", 17.701723},
667 {19, 9, "F19", 17.69690},
668 {20, 10, "Ne20", 18.62284},
669 {21, 10, "Ne21", 19.555644},
670 {22, 10, "Ne22", 20.484846},
671 {23, 11, "Na23", 21.41483},
672 {24, 12, "Mg24", 22.34193},
673 {25, 12, "Mg25", 23.274160},
674 {26, 12, "Mg26", 24.202632},
675 {27, 13, "Al27", 25.13314},
676 {28, 14, "Si28", 26.06034},
677 {29, 14, "Si29", 26.991434},
678 {30, 14, "Si30", 27.920390},
679 {31, 15, "P31", 28.85188},
680 {32, 16, "S32", 29.78180},
681 {33, 16, "S33", 30.712719},
682 {34, 16, "S34", 31.640868},
683 {35, 17, "Cl35", 32.57328},
684 {36, 18, "Ar36", 33.50356},
685 {37, 18, "Ar37", 34.434334},
686 {38, 18, "Ar38", 35.362061},
687 {39, 19, "K39", 36.29447},
688 {40, 20, "Ca40", 37.22492},
689 {41, 20, "Ca41", 38.156120},
690 {42, 20, "Ca42", 39.084205},
691 {43, 20, "Ca43", 40.015838},
692 {44, 20, "Ca44", 40.944272},
693 {45, 21, "Sc45", 41.87617},
694 {46, 22, "Ti46", 42.804605},
695 {47, 22, "Ti47", 43.735290},
696 /*
697 { 48, 22, "Ti48", 44.66324}, { 49, 22, "Ti49", 45.594652}, { 50, 22, "Ti50", 46.523278},
698 { 51, 23, "V51", 47.45401}, { 52, 24, "Cr52", 48.38228}, { 53, 24, "Cr53", 49.313903},
699 { 54, 24, "Cr54", 50.243749}, { 55, 25, "Mn55", 51.17447}, { 56, 26, "Fe56", 52.10307},
700 { 57, 26, "Fe57", 53.034984}, { 58, 26, "Fe58", 53.964505},
701 { 59, 27, "Co59", 54.89593}, { 60, 28, "Ni60", 55.825174}, { 61, 28, "Ni61", 56.756919},
702 { 62, 28, "Ni62", 57.685888}, { 63, 29, "Cu63", 58.61856}, { 64, 29, "Cu64", 59.550198},
703 { 65, 29, "Cu65", 60.479853}, { 66, 30, "Zn66", 61.409711}, { 67, 30, "Zn67", 62.342224},
704 { 68, 30, "Zn68", 63.271592}, { 69, 31, "Ga69", 64.203765}, { 70, 31, "Ga70", 65.135677},
705 { 71, 31, "Ga71", 66.065941}, { 72, 32, "Ge72", 66.994989}, { 73, 32, "Ge73", 68.857141},
706 { 74, 32, "Ge74", 68.85715}, { 75, 33, "As75", 69.789025}, { 76, 34, "Se76", 70.718300},
707 { 77, 34, "Se77", 71.650446}, { 78, 34, "Se78", 72.579514}, { 79, 34, "Se79", 73.512116},
708 { 80, 34, "Se80", 74.44178}, { 81, 35, "Br81", 75.373047}, { 82, 36, "Kr82", 76.301927},
709 { 83, 36, "Kr83", 77.234029}, { 84, 36, "Kr84", 78.16309}, { 85, 37, "Rb85", 79.09483},
710 { 86, 38, "Sr86", 80.023969}, { 87, 38, "Sr87", 80.955106}, { 88, 38, "Sr88", 81.88358},
711 { 89, 39, "Y89", 82.81527}, { 90, 40, "Zr90", 83.74571}, { 91, 40, "Zr91", 84.678073},
712 { 92, 40, "Zr92", 85.609003}, { 93, 41, "Nb93", 86.541743}, { 94, 41, "Nb94", 87.474081},
713 { 95, 42, "Mo95", 88.404232}, { 96, 42, "Mo96", 89.334643}, { 97, 42, "Mo97", 90.267388},
714 { 98, 42, "Mo98", 91.19832}, { 99, 43, "Tc99", 92.13059}, { 100, 44, "Ru100", 93.060191},
715 { 101, 44, "Ru101", 93.992955}, { 102, 44, "Ru102", 94.923300}, { 103, 45, "Rh103", 95.855870},
716 { 104, 46, "Pd104", 96.785997}, { 105, 46, "Pd105", 97.718468}, { 106, 46, "Pd106", 98.64997},
717 { 107, 47, "Ag107", 99.581467}, { 108, 47, "Ag108", 100.51376}, { 109, 47, "Ag109", 101.444134},
718 { 110, 48, "Cd110", 102.37400}, { 111, 48, "Cd111", 103.30659}, { 112, 48, "Cd112", 104.23676},
719 { 113, 48, "Cd113", 105.16978}, { 114, 48, "Cd114", 106.10997}, { 115, 49, "In115", 107.03228},
720 { 116, 49, "In116", 107.96507}, { 117, 49, "In117", 108.89586}, { 118, 50, "Sn118", 109.82465},
721 { 119, 50, "Sn119", 110.75773}, { 120, 50, "Sn120", 111.68821},
722 { 121, 51, "Sb121", 112.62119}, { 122, 51, "Sb122", 113.55395}, { 123, 51, "Sb123", 114.48455},
723 { 124, 52, "Te124", 115.41474}, { 125, 52, "Te125", 116.3477}, { 126, 52, "Te126", 117.27819},
724 { 127, 53, "I127", 118.21077}, { 128, 53, "I128", 119.14351}, { 129, 53, "I129", 120.07424 },
725 { 130, 54, "Xe130", 121.00435}, { 131, 54, "Xe131", 121.93731}, { 132, 54, "Xe132", 122.86796},
726 { 133, 55, "Cs133", 123.80064}, { 134, 55, "Cs134", 124.73332}, { 135, 55, "Cs135", 125.66412195},
727 { 136, 56, "Ba136", 126.59431}, { 137, 56, "Ba137", 127.52697}, { 138, 56, "Ba138", 128.45793},
728 { 139, 57, "La139", 129.39045}, { 140, 58, "Ce140", 130.32111}, { 141, 59, "Pr141", 131.2546475},
729 { 142, 60, "Nd142", 132.18621}, { 143, 60, "Nd143", 133.11965}, { 144, 60, "Nd144", 134.05140},
730 { 145, 60, "Nd145", 134.9852}, { 146, 60, "Nd146", 135.91721}, { 147, 60, "Nd147", 136.85148},
731 { 148, 61, "Pm148", 137.78426}, { 149, 61, "Pm149", 138.71655},
732 { 150, 62, "Sm150", 139.64706}, { 151, 62, "Sm151", 140.58103}, { 152, 62, "Sm152", 141.51236},
733 { 153, 63, "Eu153", 142.44522}, { 154, 64, "Gd154", 143.37638}, { 155, 64, "Gd155", 144.30951},
734 { 156, 64, "Gd156", 145.24054}, { 157, 64, "Gd157", 146.17374}, { 158, 65, "Tb158", 147.10659},
735 { 159, 65, "Tb159", 148.03802}, { 160, 65, "Tb160", 148.97121}, { 161, 65, "Tb161", 149.90308},
736 { 162, 66, "Dy162", 150.83386}, { 163, 66, "Dy163", 151.76715}, { 164, 66, "Dy164", 152.69909},
737 { 165, 67, "Ho165", 153.63162}, { 166, 68, "Er166", 154.56309}, { 167, 68, "Er167", 155.49622},
738 { 168, 68, "Er168", 156.42801}, { 169, 69, "Tm169", 157.36122}, { 170, 69, "Tm170", 158.29420},
739 { 171, 69, "Tm171", 159.22628}, { 172, 70, "Yb172", 160.15773}, { 173, 70, "Yb173", 161.09092},
740 { 174, 70, "Yb174", 162.02245}, { 175, 71, "Lu175", 162.95630}, { 176, 72, "Hf176", 163.88838},
741 { 177, 72, "Hf177", 164.82157}, { 178, 72, "Hf178", 165.75351}, { 179, 72, "Hf179", 166.68697},
742 { 180, 73, "Ta180", 167.62000}, { 181, 73, "Ta181", 168.55199}, { 182, 74, "W182", 169.48368},
743 { 183, 74, "W183", 170.41705}, { 184, 74, "W184", 171.34924}, { 185, 75, "Re185", 172.28258},
744 { 186, 76, "Os186", 173.21490}, { 187, 76, "Os187", 174.14818}, { 188, 76, "Os188", 175.079754},
745 { 189, 76, "Os189", 176.01340}, { 190, 76, "Os190", 176.94517}, { 191, 77, "Ir191", 177.878667},
746 { 192, 78, "Pt192", 178.81057}, { 193, 78, "Pt193", 179.74388}, { 194, 78, "Pt194", 180.67513},
747 { 195, 78, "Pt195", 181.60855}, { 196, 78, "Pt196", 182.54020}, { 197, 79, "Au197", 183.47324},
748 { 198, 80, "Hg198", 184.40488}, { 199, 80, "Hg199", 185.33778}, { 200, 80, "Hg200", 186.26932},
749 { 201, 80, "Hg201", 187.20265}, { 202, 80, "Hg202", 188.13451}, { 203, 81, "Tl203", 189.06754},
750 { 204, 81, "Tl204", 190.00045}, { 205, 81, "Tl205", 190.93247}, { 206, 82, "Pb206", 191.86400},
751 { 207, 82, "Pb207", 192.79683}, { 208, 82, "Pb208", 193.72907}};
752 */
753 {48, 23, "V48", 44.66324},
754 {49, 23, "V49", 45.594652},
755 {50, 24, "Cr50", 46.523278},
756 {51, 24, "Cr51", 47.45401},
757 {52, 25, "Mn52", 48.38228},
758 {53, 25, "Mn53", 49.313903},
759 {54, 26, "Fe54", 50.243749},
760 {55, 26, "Fe55", 51.17447},
761 {56, 27, "Co56", 52.10307},
762 {57, 27, "Co57", 53.034984},
763 {58, 27, "Co58", 53.964505},
764 {59, 28, "Ni59", 54.89593},
765 {60, 28, "Ni60", 55.825174},
766 {61, 28, "Ni61", 56.756919},
767 {62, 29, "Cu62", 57.685888},
768 {63, 29, "Cu63", 58.61856},
769 {64, 30, "Zn64", 59.550198},
770 {65, 30, "Zn65", 60.479853},
771 {66, 31, "Ga66", 61.409711},
772 {67, 31, "Ga67", 62.342224},
773 {68, 32, "Ge68", 63.271592},
774 {69, 32, "Ge69", 64.203765},
775 {70, 33, "As70", 65.135677},
776 {71, 33, "As71", 66.065941},
777 {72, 34, "Se72", 66.994989},
778 {73, 34, "Se73", 68.857141},
779 {74, 34, "Se74", 68.85715},
780 {75, 35, "Br75", 69.789025},
781 {76, 35, "Br76", 70.718300},
782 {77, 36, "Kr77", 71.650446},
783 {78, 36, "Kr78", 72.579514},
784 {79, 36, "Kr79", 73.512116},
785 {80, 37, "Rb80", 74.44178},
786 {81, 37, "Rb81", 75.373047},
787 {82, 38, "Sr82", 76.301927},
788 {83, 38, "Sr83", 77.234029},
789 {84, 39, "Y84", 78.16309},
790 {85, 39, "Y85", 79.09483},
791 {86, 40, "Zr86", 80.023969},
792 {87, 40, "Zr87", 80.955106},
793 {88, 40, "Zr88", 81.88358},
794 {89, 41, "Nb89", 82.81527},
795 {90, 42, "Mo90", 83.74571},
796 {91, 42, "Mo91", 84.678073},
797 {92, 43, "Tc92", 85.609003},
798 {93, 43, "Tc93", 86.541743},
799 {94, 43, "Tc94", 87.474081},
800 {95, 44, "Ru95", 88.404232},
801 {96, 44, "Ru96", 89.334643},
802 {97, 45, "Rh97", 90.267388},
803 {98, 45, "Rh98", 91.19832},
804 {99, 45, "Rh99", 92.13059},
805 {100, 46, "Pd100", 93.060191},
806 {101, 46, "Pd101", 93.992955},
807 {102, 47, "Ag102", 94.923300},
808 {103, 47, "Ag103", 95.855870},
809 {104, 47, "Ag104", 96.785997},
810 {105, 48, "Cd105", 97.718468},
811 {106, 48, "Cd106", 98.64997},
812 {107, 49, "In107", 99.581467},
813 {108, 49, "In108", 100.51376},
814 {109, 50, "Sn109", 101.444134},
815 {110, 50, "Sn110", 102.37400},
816 {111, 51, "Sb111", 103.30659},
817 {112, 51, "Sb112", 104.23676},
818 {113, 52, "Te113", 105.16978},
819 {114, 52, "Te114", 106.10997},
820 {115, 52, "Te115", 107.03228},
821 {116, 53, "I116", 107.96507},
822 {117, 53, "I117", 108.89586},
823 {118, 54, "Xe118", 109.82465},
824 {119, 54, "Xe119", 110.75773},
825 {120, 55, "Cs120", 111.68821},
826 {121, 55, "Cs121", 112.62119},
827 {122, 56, "Ba122", 113.55395},
828 {123, 56, "Ba123", 114.48455},
829 {124, 56, "Ba124", 115.41474},
830 {125, 57, "La125", 116.3477},
831 {126, 57, "La126", 117.27819},
832 {127, 58, "Ce127", 118.21077},
833 {128, 58, "Ce128", 119.14351},
834 {129, 58, "Ce129", 120.07424},
835 {130, 59, "Pr130", 121.00435},
836 {131, 59, "Pr131", 121.93731},
837 {132, 60, "Nd132", 122.86796},
838 {133, 60, "Nd133", 123.80064},
839 {134, 60, "Nd134", 124.73332},
840 {135, 61, "Pm135", 125.66412195},
841 {136, 62, "Sm136", 126.59431},
842 {137, 62, "Sm137", 127.52697},
843 {138, 62, "Sm138", 128.45793},
844 {139, 63, "Eu139", 129.39045},
845 {140, 63, "Eu140", 130.32111},
846 {141, 64, "Gd141", 131.2546475},
847 {142, 64, "Gd142", 132.18621},
848 {143, 64, "Gd143", 133.11965},
849 {144, 65, "Tb144", 134.05140},
850 {145, 65, "Tb145", 134.9852},
851 {146, 66, "Dy146", 135.91721},
852 {147, 66, "Dy147", 136.85148},
853 {148, 66, "Dy148", 137.78426},
854 {149, 67, "Ho149", 138.71655},
855 {150, 67, "Ho150", 139.64706},
856 {151, 67, "Ho151", 140.58103},
857 {152, 68, "Er152", 141.51236},
858 {153, 68, "Er153", 142.44522},
859 {154, 69, "Tm154", 143.37638},
860 {155, 69, "Tm155", 144.30951},
861 {156, 69, "Tm156", 145.24054},
862 {157, 70, "Yb157", 146.17374},
863 {158, 70, "Yb158", 147.10659},
864 {159, 71, "Lu159", 148.03802},
865 {160, 71, "Lu160", 148.97121},
866 {161, 71, "Lu161", 149.90308},
867 {162, 72, "Hf162", 150.83386},
868 {163, 72, "Hf163", 151.76715},
869 {164, 72, "Hf164", 152.69909},
870 {165, 72, "Hf165", 153.63162},
871 {166, 73, "Ta166", 154.56309},
872 {167, 73, "Ta167", 155.49622},
873 {168, 73, "Ta168", 156.42801},
874 {169, 74, "W169", 157.36122},
875 {170, 74, "W170", 158.29420},
876 {171, 75, "Re171", 159.22628},
877 {172, 75, "Re172", 160.15773},
878 {173, 75, "Re173", 161.09092},
879 {174, 75, "Re174", 162.02245},
880 {175, 75, "Re175", 162.95630},
881 {176, 75, "Re176", 163.88838},
882 {177, 76, "Os177", 164.82157},
883 {178, 76, "Os178", 165.75351},
884 {179, 76, "Os179", 166.68697},
885 {180, 77, "Ir180", 167.62000},
886 {181, 77, "Ir181", 168.55199},
887 {182, 77, "Ir182", 169.48368},
888 {183, 77, "Ir183", 170.41705},
889 {184, 77, "Ir184", 171.34924},
890 {185, 78, "Pt185", 172.28258},
891 {186, 78, "Pt186", 173.21490},
892 {187, 78, "Pt187", 174.14818},
893 {188, 78, "Pt188", 175.079754},
894 {189, 78, "Pt189", 176.01340},
895 {190, 78, "Pt190", 176.94517},
896 {191, 78, "Pt191", 177.878667},
897 {192, 79, "Au192", 178.81057},
898 {193, 79, "Au193", 179.74388},
899 {194, 79, "Au194", 180.67513},
900 {195, 79, "Au195", 181.60855},
901 {196, 79, "Au196", 182.54020},
902 {197, 79, "Au197", 183.47324},
903 {198, 80, "Hg198", 184.40488},
904 {199, 80, "Hg199", 185.33778},
905 {200, 80, "Hg200", 186.26932},
906 {201, 80, "Hg201", 187.20265},
907 {202, 80, "Hg202", 188.13451},
908 {203, 81, "Tl203", 189.06754},
909 {204, 81, "Tl204", 190.00045},
910 {205, 81, "Tl205", 190.93247},
911 {206, 82, "Pb206", 191.86400},
912 {207, 82, "Pb207", 192.79683},
913 {208, 82, "Pb208", 193.72907}};
914 TDatabasePDG* dbPDG = TDatabasePDG::Instance();
915 for (Int_t i = 0; i < NSI; i++) {
916 if (ions[i].Mass < 0.1)
917 continue; // expected that it is already registered
918 Int_t ID = 1000000000 + ions[i].Z * 10000 + ions[i].N * 10;
919 if (dbPDG->GetParticle(ID))
920 continue;
921 FairIon* ion = new FairIon(ions[i].Name, ions[i].Z, ions[i].N, ions[i].Z);
922 FairRunSim::Instance()->AddNewIon(ion);
923 }
924 // Min and max baryon charge for a given Z
925 for (Int_t i = 0; i < NSI; i++) {
926 Int_t z = ions[i].Z;
927 if (z < 0 && z > fZMax)
928 continue;
929 fN2[z] = ions[i].N;
930 }
931 for (Int_t i = NSI - 1; i >= 0; i--) {
932 Int_t z = ions[i].Z;
933 if (z < 0 && z > fZMax)
934 continue;
935 fN1[z] = ions[i].N;
936 }
937
938 // Min and max Z for a given baryon charge
939 for (Int_t i = 0; i < NSI; i++) {
940 Int_t B = ions[i].N;
941 if (B < 0 && B > fBMax)
942 continue;
943 fZ2[B] = ions[i].Z;
944 }
945 for (Int_t i = NSI - 1; i >= 0; i--) {
946 Int_t B = ions[i].N;
947 if (B < 0 && B > fBMax)
948 continue;
949 fZ1[B] = ions[i].Z;
950 }
951
952 for (int i = 2; i <= fBMax; i++)
953 cout << i << ": " << fZ1[i] << "-" << fZ2[i] << " ";
954 cout << "\n";
955
956 return NSI;
957}
friend F32vec4 fabs(const F32vec4 &a)
Definition P4_F32vec4.h:52
int i
Definition P4_F32vec4.h:22
__m128 m
Definition P4_F32vec4.h:27
friend F32vec4 atan2(const F32vec4 &y, const F32vec4 &x)
Definition P4_F32vec4.h:130
Int_t FindPDGCodeSpectator(Int_t N, Int_t Z, Int_t &dN)
Int_t FindPDGCodeParticipant(Int_t A, Int_t S, Int_t Z, Float_t mass, Double_t &massFactor)
Bool_t ReadEvent(FairPrimaryGenerator *primGen)
Bool_t SkipEvents(Int_t count)
STL namespace.