CbmRoot
CbmPhsdGenerator.cxx
Go to the documentation of this file.
1 // -------------------------------------------------------------------------
2 // ----- FairUrqmdGenerator source file -----
3 // ----- Created 24/01/14 by V. Vovchenko -----
4 // modified 09/2017 by I.Vassiliev
5 //-------------------------------------------------------------------------
6 #include "CbmPhsdGenerator.h"
7 
8 #include "FairMCEventHeader.h"
9 #include "FairPrimaryGenerator.h"
10 
11 #include "TDatabasePDG.h"
12 #include "TLorentzVector.h"
13 #include "TMCProcess.h"
14 #include "TObjArray.h"
15 #include "TPDGCode.h"
16 #include "TParticle.h"
17 #include "TParticlePDG.h"
18 #include "TRandom.h"
19 #include "TString.h"
20 #include "TVirtualMCStack.h"
21 
22 
23 #include <cstring>
24 #include <iostream>
25 #include <string>
26 
27 using std::cout;
28 using std::endl;
29 using std::string;
30 
31 const Double_t kProtonMass = 0.938271998;
32 
33 
34 // ----- Default constructor ------------------------------------------
36  : FairGenerator()
37  ,
38  //fInputFile(nullptr),
39  fBaryonsFile(nullptr)
40  , fMesonsFile(nullptr)
41  , fParticleTable()
42  , fFileName(nullptr) {}
43 // ------------------------------------------------------------------------
44 
45 
46 // ----- Standard constructor -----------------------------------------
47 CbmPhsdGenerator::CbmPhsdGenerator(const char* fileNameInput,
48  const char* fileNameBaryons,
49  const char* fileNameMesons)
50  : FairGenerator()
51  , fBaryonsFile(nullptr)
52  , fMesonsFile(nullptr)
53  , fDatFile(nullptr)
54  , fParticleTable() {
55  // cout << "-I CbmPhsdGenerator: Opening input file " << fileNameInput << endl;
56  cout << "-I CbmPhsdGenerator: Opening HSD Baryons file " << fileNameBaryons
57  << endl;
58  fBaryonsFile = fopen(fileNameBaryons, "r");
59  if (!fBaryonsFile) {
60  Fatal("FairHSDgenerator", "Cannot open HSD Baryons file.");
61  }
62  cout << "-I CbmPhsdGenerator: Opening HSD Mesons file " << fileNameMesons
63  << endl;
64  fMesonsFile = fopen(fileNameMesons, "r");
65  if (!fMesonsFile) {
66  Fatal("FairHSDgenerator", "Cannot open HSD Mesons file.");
67  }
69  ReadCollisionData(fileNameInput);
70  nextBaryon.init = false;
71  nextMeson.init = false;
72  nextEvent = 1;
73  fReadDat = false;
74 }
75 // ------------------------------------------------------------------------
76 
77 
78 // ----- Standard constructor -----------------------------------------
79 CbmPhsdGenerator::CbmPhsdGenerator(const char* fileNameInput,
80  const char* fileNameDat)
81  : FairGenerator()
82  , fBaryonsFile(nullptr)
83  , fMesonsFile(nullptr)
84  , fDatFile(nullptr)
85  , fParticleTable() {
86  // cout << "-I CbmPhsdGenerator: Opening input file " << fileNameInput << endl;
87  cout << "-I CbmPhsdGenerator: Opening phsd.dat file " << fileNameDat << endl;
88  fDatFile = fopen(fileNameDat, "r");
89  if (!fDatFile) { Fatal("CbmPhsdGenerator", "Cannot open phsd.dat file."); }
90  ReadCollisionData(fileNameInput);
91  nextEvent = 1;
92  fReadDat = true;
93 }
94 // ------------------------------------------------------------------------
95 
96 
97 // ----- Destructor ---------------------------------------------------
99  // cout<<"Enter Destructor of CbmPhsdGenerator"<<endl;
100  /*if ( fInputFile ) {
101  fclose(fInputFile);
102  fInputFile = nullptr;
103  }*/
104  if (fBaryonsFile) {
105  fclose(fBaryonsFile);
106  fBaryonsFile = nullptr;
107  }
108  if (fMesonsFile) {
109  fclose(fMesonsFile);
110  fMesonsFile = nullptr;
111  }
112  if (fDatFile) {
113  fclose(fDatFile);
114  fDatFile = nullptr;
115  }
116  fParticleTable.clear();
117  // cout<<"Leave Destructor of CbmPhsdGenerator"<<endl;
118 }
119 // ------------------------------------------------------------------------
120 
121 
122 // ----- Public method ReadEvent --------------------------------------
123 Bool_t CbmPhsdGenerator::ReadEvent(FairPrimaryGenerator* primGen) {
124  if (fReadDat)
125  return ReadEventDat(primGen);
126  else
127  return ReadEvent300(primGen);
128 }
129 // ------------------------------------------------------------------------
130 
131 
132 // ----- Public method ReadEventDat -----------------------------------
133 Bool_t CbmPhsdGenerator::ReadEventDat(FairPrimaryGenerator* primGen) {
134 
135  // ---> Check for input file
136  if (!fDatFile) {
137  cout << "-E CbmPhsdGenerator: phsd.dat input file is not open! " << endl;
138  return kFALSE;
139  }
140 
141 
142  int nTracks;
143  double b;
144 
145  if (fscanf(fDatFile, "%d %*d %*d %lf%*[^\n]%*c", &nTracks, &b) == EOF) {
146  cout << "-E- CbmPhsdGenerator::ReadEvent: "
147  << "No more events!" << endl;
148  return kFALSE;
149  }
150 
151  fscanf(fDatFile, "%*[^\n]%*c");
152 
153  cout << "Event: " << nextEvent << " nTracks: " << nTracks << "\n";
154 
155 
156  //int pid = 0;
157 
158  // Set event id and impact parameter in MCEvent if not yet done
159  FairMCEventHeader* event = primGen->GetEvent();
160  if (event && (!event->IsSet())) {
161  event->SetEventID(nextEvent);
162  event->SetB(b);
163  event->MarkSet(kTRUE);
164  }
165 
166  for (int i = 0; i < nTracks; ++i) {
167  int ttype, tchr;
168  double tmppz, tmppx, tmppy, tmpp0;
169  fscanf(fDatFile,
170  "%d%d%lf%lf%lf%lf%*[^\n]%*c",
171  &ttype,
172  &tchr,
173  &tmppx,
174  &tmppy,
175  &tmppz,
176  &tmpp0);
177 
178  Int_t pdgID = ttype;
179 
180  // cout << "-I CbmPhsdGenerator: PID " << pdgID << endl;
181 
182  Double_t eBeam = ekin + kProtonMass;
183  Double_t pBeam = TMath::Sqrt(eBeam * eBeam - kProtonMass * kProtonMass);
184  Double_t betaCM = pBeam / (eBeam + kProtonMass);
185  Double_t gammaCM = TMath::Sqrt(1. / (1. - betaCM * betaCM));
186 
187  // Lorentztransformation to lab
188 
189  Double_t px = Double_t(tmppx);
190  Double_t py = Double_t(tmppy);
191  Double_t pz = Double_t(tmppz);
192  Double_t e = Double_t(tmpp0);
193  Double_t mass = sqrt(e * e - px * px - py * py - pz * pz);
194  pz = gammaCM * (pz + betaCM * e);
195  Double_t ee = sqrt(mass * mass + px * px + py * py + pz * pz);
196 
197  TVector3 aa(px, py, pz);
198  TLorentzVector pp;
199  pp.SetPx(px);
200  pp.SetPy(py);
201  pp.SetPz(pz);
202  pp.SetE(ee);
203 
204  // Give track to PrimaryGenerator
205  primGen->AddTrack(pdgID, px, py, pz, 0., 0., 0.);
206  }
207 
208  cout << "-I CbmPhsdGenerator: Event " << nextEvent << ", b = " << b
209  << " fm, multiplicity " << nTracks << ", ekin: " << ekin << endl;
210 
211  nextEvent++;
212 
213  return kTRUE;
214 }
215 // ------------------------------------------------------------------------
216 
217 
218 // ----- Public method ReadEvent300 -----------------------------------
219 Bool_t CbmPhsdGenerator::ReadEvent300(FairPrimaryGenerator* primGen) {
220 
221  // ---> Check for input file
222  if (!fBaryonsFile) {
223  cout << "-E CbmPhsdGenerator: Baryons input file is not open! " << endl;
224  return kFALSE;
225  }
226 
227  if (!fMesonsFile) {
228  cout << "-E CbmPhsdGenerator: Mesons input file is not open! " << endl;
229  return kFALSE;
230  }
231 
232  // ---> Check for primary generator
233  if (!primGen) {
234  cout << "-E- CbmPhsdGenerator::ReadEvent: "
235  << "No PrimaryGenerator!" << endl;
236  return kFALSE;
237  }
238 
239  if (!nextBaryon.init && feof(fBaryonsFile)) {
240  cout << "-E- CbmPhsdGenerator::ReadEvent: "
241  << "No more events!" << endl;
242  return kFALSE;
243  }
244 
245  /*if (!nextMeson.init && feof(fMesonsFile)) {
246  cout << "-E- CbmPhsdGenerator::ReadEvent: "
247  << "No more events!" << endl;
248  return kFALSE;
249  }*/
250 
251  if (!nextBaryon.init) {
252  fscanf(fBaryonsFile,
253  "%d %d %d %d %lf %lf %lf %lf %lf",
254  &nextBaryon.id,
256  &nextBaryon.ISUB,
257  &nextBaryon.IRUN,
258  &nextBaryon.px,
259  &nextBaryon.py,
260  &nextBaryon.pz,
261  &nextBaryon.p0,
262  &nextBaryon.b);
263  fscanf(fBaryonsFile, "%*[^\n]%*c");
265  nextBaryon.init = 1;
266  }
267 
268  if (!nextMeson.init && !feof(fMesonsFile)) {
269  fscanf(fMesonsFile,
270  "%d %d %d %d %lf %lf %lf %lf %lf",
271  &nextMeson.id,
272  &nextMeson.charge,
273  &nextMeson.ISUB,
274  &nextMeson.IRUN,
275  &nextMeson.px,
276  &nextMeson.py,
277  &nextMeson.pz,
278  &nextMeson.p0,
279  &nextMeson.b);
280  fscanf(fMesonsFile, "%*[^\n]%*c");
282  nextMeson.init = 1;
283  }
284 
285 
286  int nTracks = 0, pid = 0;
287  double b = nextBaryon.b;
288 
289  // Set event id and impact parameter in MCEvent if not yet done
290  FairMCEventHeader* event = primGen->GetEvent();
291  if (event && (!event->IsSet())) {
292  event->SetEventID(nextEvent);
293  event->SetB(b);
294  event->MarkSet(kTRUE);
295  }
296 
297  // Read all baryons from current event until we reach the next event in fort.300 or end of file
298  while (nextBaryon.globalEvent == nextEvent) {
299  // Convert HSD type and charge to unique pid identifier which is based on
300  // HSD particle id and charge, calculated separately for baryons, anti-baryons and mesons
301  if (nextBaryon.id >= 0) {
302  pid = nextBaryon.id * 10 + (2 + nextBaryon.charge);
303  }
304  // antibaryons
305  else {
306  pid = -(-nextBaryon.id * 10 + (2 - nextBaryon.charge));
307  }
308 
309  // Convert Unique PID into PDG particle code
310  if (fParticleTable.find(pid) == fParticleTable.end()) {
311  cout << "-W CbmPhsdGenerator: PID " << nextBaryon.id << " charge "
312  << nextBaryon.charge << " not found in table (" << pid << ")"
313  << endl;
314  if (feof(fBaryonsFile)) {
315  nextBaryon.init = 0;
316  break;
317  }
318  fscanf(fBaryonsFile,
319  "%d %d %d %d %lf %lf %lf %lf %lf",
320  &nextBaryon.id,
322  &nextBaryon.ISUB,
323  &nextBaryon.IRUN,
324  &nextBaryon.px,
325  &nextBaryon.py,
326  &nextBaryon.pz,
327  &nextBaryon.p0,
328  &nextBaryon.b);
329  fscanf(fBaryonsFile, "%*[^\n]%*c");
331  continue;
332  }
333  Int_t pdgID = fParticleTable[pid];
334 
335  Double_t eBeam = ekin + kProtonMass;
336  Double_t pBeam = TMath::Sqrt(eBeam * eBeam - kProtonMass * kProtonMass);
337  Double_t betaCM = pBeam / (eBeam + kProtonMass);
338  Double_t gammaCM = TMath::Sqrt(1. / (1. - betaCM * betaCM));
339 
340  // Lorentztransformation to lab
341 
342  Double_t px = Double_t(nextBaryon.px);
343  Double_t py = Double_t(nextBaryon.py);
344  Double_t pz = Double_t(nextBaryon.pz);
345  Double_t e = Double_t(nextBaryon.p0);
346  Double_t mass = sqrt(e * e - px * px - py * py - pz * pz);
347  pz = gammaCM * (pz + betaCM * e);
348  Double_t ee = sqrt(mass * mass + px * px + py * py + pz * pz);
349 
350  TVector3 aa(px, py, pz);
351  TLorentzVector pp;
352  pp.SetPx(px);
353  pp.SetPy(py);
354  pp.SetPz(pz);
355  pp.SetE(ee);
356 
357  // Give track to PrimaryGenerator
358  primGen->AddTrack(pdgID, px, py, pz, 0., 0., 0.);
359 
360  nTracks++;
361 
362  if (feof(fBaryonsFile)) {
363  nextBaryon.init = 0;
364  break;
365  }
366 
367  fscanf(fBaryonsFile,
368  "%d %d %d %d %lf %lf %lf %lf %lf",
369  &nextBaryon.id,
371  &nextBaryon.ISUB,
372  &nextBaryon.IRUN,
373  &nextBaryon.px,
374  &nextBaryon.py,
375  &nextBaryon.pz,
376  &nextBaryon.p0,
377  &nextBaryon.b);
378  fscanf(fBaryonsFile, "%*[^\n]%*c");
380  }
381 
382  // Read all mesons from current event until we reach the next event in fort.301 or end of file
383  while (nextMeson.globalEvent == nextEvent) {
384  // Convert HSD type and charge to unique pid identifier which is based on
385  // HSD particle id and charge, calculated separately for baryons, anti-baryons and mesons
386  { pid = 1000 + nextMeson.id * 10 + (2 + nextMeson.charge); }
387 
388  // Convert Unique PID into PDG particle code
389  if (fParticleTable.find(pid) == fParticleTable.end()) {
390  cout << "-W CbmPhsdGenerator: PID " << nextMeson.id << " charge "
391  << nextMeson.charge << " not found in table (" << pid << ")" << endl;
392  fscanf(fMesonsFile,
393  "%d %d %d %d %lf %lf %lf %lf %lf",
394  &nextMeson.id,
395  &nextMeson.charge,
396  &nextMeson.ISUB,
397  &nextMeson.IRUN,
398  &nextMeson.px,
399  &nextMeson.py,
400  &nextMeson.pz,
401  &nextMeson.p0,
402  &nextMeson.b);
403  fscanf(fMesonsFile, "%*[^\n]%*c");
405  continue;
406  }
407  Int_t pdgID = fParticleTable[pid];
408 
409  Double_t eBeam = ekin + kProtonMass;
410  Double_t pBeam = TMath::Sqrt(eBeam * eBeam - kProtonMass * kProtonMass);
411  Double_t betaCM = pBeam / (eBeam + kProtonMass);
412  Double_t gammaCM = TMath::Sqrt(1. / (1. - betaCM * betaCM));
413 
414  // Lorentztransformation to lab
415 
416  Double_t px = Double_t(nextMeson.px);
417  Double_t py = Double_t(nextMeson.py);
418  Double_t pz = Double_t(nextMeson.pz);
419  Double_t e = Double_t(nextMeson.p0);
420  Double_t mass = sqrt(e * e - px * px - py * py - pz * pz);
421  pz = gammaCM * (pz + betaCM * e);
422  Double_t ee = sqrt(mass * mass + px * px + py * py + pz * pz);
423 
424  TVector3 aa(px, py, pz);
425  TLorentzVector pp;
426  pp.SetPx(px);
427  pp.SetPy(py);
428  pp.SetPz(pz);
429  pp.SetE(ee);
430 
431  // Give track to PrimaryGenerator
432  primGen->AddTrack(pdgID, px, py, pz, 0., 0., 0.);
433 
434  nTracks++;
435 
436  if (feof(fMesonsFile)) {
437  nextMeson.init = 0;
438  break;
439  }
440 
441  fscanf(fMesonsFile,
442  "%d %d %d %d %lf %lf %lf %lf %lf",
443  &nextMeson.id,
444  &nextMeson.charge,
445  &nextMeson.ISUB,
446  &nextMeson.IRUN,
447  &nextMeson.px,
448  &nextMeson.py,
449  &nextMeson.pz,
450  &nextMeson.p0,
451  &nextMeson.b);
452  fscanf(fMesonsFile, "%*[^\n]%*c");
454  }
455 
456  cout << "-I CbmPhsdGenerator: Event " << nextEvent << ", b = " << b
457  << " fm, multiplicity " << nTracks << ", ekin: " << ekin << endl;
458 
459  nextEvent++;
460 
461  return kTRUE;
462 }
463 // ------------------------------------------------------------------------
464 
531 // ------------------------------------------------------------------------
532 
533 // ----- Private method ReadConversionTable ---------------------------
535 
536  TString work = getenv("VMCWORKDIR");
537  TString fileName = work + "/input/hsd_pdg.dat";
538  std::ifstream* pdgconv = new std::ifstream(fileName.Data());
539 
540  if (!(*pdgconv).is_open()) {
541  Fatal("CbmPhsdGenerator", "Particle table for conversion was not found!");
542  //cout << "-W CbmPhsdGenerator: Particle table for conversion was not found!" << endl;
543  }
544 
545  Int_t index = 0;
546  Int_t pdgId = 0;
547 
548  string tmpStr;
549 
550  while (!pdgconv->eof()) {
551  index = pdgId = 0;
552  *pdgconv >> index >> pdgId;
553  std::getline(*pdgconv, tmpStr);
554  fParticleTable[index] = pdgId;
555  }
556 
557  pdgconv->close();
558  delete pdgconv;
559 
560  cout << "-I CbmPhsdGenerator: Particle table for conversion from "
561  << "HSD loaded" << endl;
562 }
563 // ------------------------------------------------------------------------
564 
565 // ----- Private method ReadCollisionData ----------------------------
566 void CbmPhsdGenerator::ReadCollisionData(const char* fileNameInput) {
567 
568  std::ifstream* inputf = new std::ifstream(fileNameInput);
569 
570  // Int_t index = 0;
571  // Int_t pdgId = 0;
572 
573  //TString val, name, tname;
574  string val, name, tname, tmpStr;
575 
576  while (!inputf->eof()) {
577  // index =pdgId =0 ;
578  //*inputf >> index >> pdgId ;
579  *inputf >> val >> name;
580  tname = name.substr(0, 6);
581  //cout << tname << endl;
582  if (tname == "MASSTA")
583  At = atoi(val.substr(0, val.size() - 1).c_str());
584  else if (tname == "MSTAPR")
585  Zt = atoi(val.substr(0, val.size() - 1).c_str());
586  else if (tname == "MASSPR")
587  Ap = atoi(val.substr(0, val.size() - 1).c_str());
588  else if (tname == "MSPRPR")
589  Zp = atoi(val.substr(0, val.size() - 1).c_str());
590  tname = name.substr(0, 4);
591  if (tname == "ELAB") ekin = atof(val.substr(0, val.size() - 1).c_str());
592  tname = name.substr(0, 5);
593  if (tname == "ISUBS") ISUBS = atof(val.substr(0, val.size() - 1).c_str());
594  tname = name.substr(0, 3);
595  if (tname == "NUM") IRUNS = atof(val.substr(0, val.size() - 1).c_str());
596  std::getline(*inputf, tmpStr);
597  }
598 
599  inputf->close();
600  delete inputf;
601 
602  cout << "-I CbmPhsdGenerator: Collision data from "
603  << "HSD loaded" << endl;
604 }
605 // ------------------------------------------------------------------------
606 
607 
CbmPhsdGenerator::ReadConversionTable
void ReadConversionTable()
Definition: CbmPhsdGenerator.cxx:534
CbmPhsdGenerator::fParticleTable
std::map< Int_t, Int_t > fParticleTable
HSD output files.
Definition: CbmPhsdGenerator.h:96
Hadron::p0
double p0
Definition: CbmPhsdGenerator.h:38
CbmPhsdGenerator::Ap
int Ap
Definition: CbmPhsdGenerator.h:106
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:41
Hadron::id
int id
Definition: CbmPhsdGenerator.h:33
CbmPhsdGenerator::nextBaryon
Hadron nextBaryon
Input file name.
Definition: CbmPhsdGenerator.h:101
Hadron::init
bool init
Definition: CbmPhsdGenerator.h:32
CbmPhsdGenerator::Zp
int Zp
Definition: CbmPhsdGenerator.h:106
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
CbmPhsdGenerator::fMesonsFile
FILE * fMesonsFile
Definition: CbmPhsdGenerator.h:92
CbmPhsdGenerator::ISUBS
int ISUBS
Definition: CbmPhsdGenerator.h:107
CbmPhsdGenerator
Definition: CbmPhsdGenerator.h:43
Hadron::pz
double pz
Definition: CbmPhsdGenerator.h:38
ClassImp
ClassImp(CbmPhsdGenerator)
Hadron::IRUN
int IRUN
Definition: CbmPhsdGenerator.h:36
Hadron::px
double px
Definition: CbmPhsdGenerator.h:38
CbmPhsdGenerator::ReadEvent300
Bool_t ReadEvent300(FairPrimaryGenerator *primGen)
Definition: CbmPhsdGenerator.cxx:219
CbmPhsdGenerator::ekin
double ekin
Definition: CbmPhsdGenerator.h:109
Hadron::b
double b
Definition: CbmPhsdGenerator.h:39
CbmPhsdGenerator.h
CbmPhsdGenerator::fDatFile
FILE * fDatFile
HSD output files.
Definition: CbmPhsdGenerator.h:94
Hadron::charge
int charge
Definition: CbmPhsdGenerator.h:34
CbmPhsdGenerator::Zt
int Zt
Definition: CbmPhsdGenerator.h:106
CbmPhsdGenerator::fReadDat
Bool_t fReadDat
Definition: CbmPhsdGenerator.h:90
CbmPhsdGenerator::ReadEvent
Bool_t ReadEvent(FairPrimaryGenerator *primGen)
Definition: CbmPhsdGenerator.cxx:123
CbmPhsdGenerator::nextMeson
Hadron nextMeson
Definition: CbmPhsdGenerator.h:101
CbmPhsdGenerator::CbmPhsdGenerator
CbmPhsdGenerator()
Definition: CbmPhsdGenerator.cxx:35
CbmPhsdGenerator::ReadEventDat
Bool_t ReadEventDat(FairPrimaryGenerator *primGen)
Definition: CbmPhsdGenerator.cxx:133
Hadron::ISUB
int ISUB
Definition: CbmPhsdGenerator.h:35
CbmPhsdGenerator::IRUNS
int IRUNS
Definition: CbmPhsdGenerator.h:108
CbmPhsdGenerator::~CbmPhsdGenerator
~CbmPhsdGenerator()
Definition: CbmPhsdGenerator.cxx:98
CbmPhsdGenerator::nextEvent
int nextEvent
Definition: CbmPhsdGenerator.h:102
CbmPhsdGenerator::fBaryonsFile
FILE * fBaryonsFile
Whether phsd.dat or .300/301 files are used.
Definition: CbmPhsdGenerator.h:92
kProtonMass
const Double_t kProtonMass
Definition: CbmPhsdGenerator.cxx:31
CbmPhsdGenerator::At
int At
Definition: CbmPhsdGenerator.h:106
Hadron::py
double py
Definition: CbmPhsdGenerator.h:38
CbmPhsdGenerator::ReadCollisionData
void ReadCollisionData(const char *fileNameInput)
Definition: CbmPhsdGenerator.cxx:566
Hadron::globalEvent
int globalEvent
Definition: CbmPhsdGenerator.h:37