8 #include "FairMCEventHeader.h"
9 #include "FairPrimaryGenerator.h"
11 #include "TDatabasePDG.h"
12 #include "TLorentzVector.h"
13 #include "TMCProcess.h"
14 #include "TObjArray.h"
16 #include "TParticle.h"
17 #include "TParticlePDG.h"
20 #include "TVirtualMCStack.h"
40 , fMesonsFile(nullptr)
42 , fFileName(nullptr) {}
48 const char* fileNameBaryons,
49 const char* fileNameMesons)
51 , fBaryonsFile(nullptr)
52 , fMesonsFile(nullptr)
56 cout <<
"-I CbmPhsdGenerator: Opening HSD Baryons file " << fileNameBaryons
60 Fatal(
"FairHSDgenerator",
"Cannot open HSD Baryons file.");
62 cout <<
"-I CbmPhsdGenerator: Opening HSD Mesons file " << fileNameMesons
66 Fatal(
"FairHSDgenerator",
"Cannot open HSD Mesons file.");
80 const char* fileNameDat)
82 , fBaryonsFile(nullptr)
83 , fMesonsFile(nullptr)
87 cout <<
"-I CbmPhsdGenerator: Opening phsd.dat file " << fileNameDat << endl;
89 if (!
fDatFile) { Fatal(
"CbmPhsdGenerator",
"Cannot open phsd.dat file."); }
137 cout <<
"-E CbmPhsdGenerator: phsd.dat input file is not open! " << endl;
145 if (fscanf(
fDatFile,
"%d %*d %*d %lf%*[^\n]%*c", &nTracks, &b) == EOF) {
146 cout <<
"-E- CbmPhsdGenerator::ReadEvent: "
147 <<
"No more events!" << endl;
153 cout <<
"Event: " <<
nextEvent <<
" nTracks: " << nTracks <<
"\n";
159 FairMCEventHeader*
event = primGen->GetEvent();
160 if (event && (!event->IsSet())) {
163 event->MarkSet(kTRUE);
166 for (
int i = 0;
i < nTracks; ++
i) {
168 double tmppz, tmppx, tmppy, tmpp0;
170 "%d%d%lf%lf%lf%lf%*[^\n]%*c",
185 Double_t gammaCM = TMath::Sqrt(1. / (1. - betaCM * betaCM));
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);
197 TVector3 aa(px, py, pz);
205 primGen->AddTrack(pdgID, px, py, pz, 0., 0., 0.);
208 cout <<
"-I CbmPhsdGenerator: Event " <<
nextEvent <<
", b = " << b
209 <<
" fm, multiplicity " << nTracks <<
", ekin: " <<
ekin << endl;
223 cout <<
"-E CbmPhsdGenerator: Baryons input file is not open! " << endl;
228 cout <<
"-E CbmPhsdGenerator: Mesons input file is not open! " << endl;
234 cout <<
"-E- CbmPhsdGenerator::ReadEvent: "
235 <<
"No PrimaryGenerator!" << endl;
240 cout <<
"-E- CbmPhsdGenerator::ReadEvent: "
241 <<
"No more events!" << endl;
253 "%d %d %d %d %lf %lf %lf %lf %lf",
270 "%d %d %d %d %lf %lf %lf %lf %lf",
286 int nTracks = 0, pid = 0;
290 FairMCEventHeader*
event = primGen->GetEvent();
291 if (event && (!event->IsSet())) {
294 event->MarkSet(kTRUE);
311 cout <<
"-W CbmPhsdGenerator: PID " <<
nextBaryon.
id <<
" charge "
319 "%d %d %d %d %lf %lf %lf %lf %lf",
338 Double_t gammaCM = TMath::Sqrt(1. / (1. - betaCM * betaCM));
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);
350 TVector3 aa(px, py, pz);
358 primGen->AddTrack(pdgID, px, py, pz, 0., 0., 0.);
368 "%d %d %d %d %lf %lf %lf %lf %lf",
390 cout <<
"-W CbmPhsdGenerator: PID " <<
nextMeson.
id <<
" charge "
393 "%d %d %d %d %lf %lf %lf %lf %lf",
412 Double_t gammaCM = TMath::Sqrt(1. / (1. - betaCM * betaCM));
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);
424 TVector3 aa(px, py, pz);
432 primGen->AddTrack(pdgID, px, py, pz, 0., 0., 0.);
442 "%d %d %d %d %lf %lf %lf %lf %lf",
456 cout <<
"-I CbmPhsdGenerator: Event " <<
nextEvent <<
", b = " << b
457 <<
" fm, multiplicity " << nTracks <<
", ekin: " <<
ekin << endl;
536 TString work = getenv(
"VMCWORKDIR");
537 TString fileName = work +
"/input/hsd_pdg.dat";
538 std::ifstream* pdgconv =
new std::ifstream(fileName.Data());
540 if (!(*pdgconv).is_open()) {
541 Fatal(
"CbmPhsdGenerator",
"Particle table for conversion was not found!");
550 while (!pdgconv->eof()) {
552 *pdgconv >> index >> pdgId;
553 std::getline(*pdgconv, tmpStr);
560 cout <<
"-I CbmPhsdGenerator: Particle table for conversion from "
561 <<
"HSD loaded" << endl;
568 std::ifstream* inputf =
new std::ifstream(fileNameInput);
574 string val, name, tname, tmpStr;
576 while (!inputf->eof()) {
579 *inputf >> val >> name;
580 tname = name.substr(0, 6);
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);
602 cout <<
"-I CbmPhsdGenerator: Collision data from "
603 <<
"HSD loaded" << endl;