CbmRoot
Loading...
Searching...
No Matches
CbmPhsdGenerator.cxx
Go to the documentation of this file.
1/* Copyright (C) 2014-2020 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
2 SPDX-License-Identifier: GPL-3.0-only
3 Authors: Volodymyr Vovchenko, Florian Uhlig [committer] */
4
5// -------------------------------------------------------------------------
6// ----- FairUrqmdGenerator source file -----
7// ----- Created 24/01/14 by V. Vovchenko -----
8// modified 09/2017 by I.Vassiliev
9//-------------------------------------------------------------------------
10#include "CbmPhsdGenerator.h"
11
12#include "FairMCEventHeader.h"
13#include "FairPrimaryGenerator.h"
14
15#include "TDatabasePDG.h"
16#include "TLorentzVector.h"
17#include "TMCProcess.h"
18#include "TObjArray.h"
19#include "TPDGCode.h"
20#include "TParticle.h"
21#include "TParticlePDG.h"
22#include "TRandom.h"
23#include "TString.h"
24#include "TVirtualMCStack.h"
25
26#include <cstring>
27#include <iostream>
28#include <string>
29
30using std::cout;
31using std::endl;
32using std::string;
33
34const Double_t kProtonMass = 0.938271998;
35
36
37// ----- Default constructor ------------------------------------------
38CbmPhsdGenerator::CbmPhsdGenerator()
39 : FairGenerator()
40 ,
41 //fInputFile(nullptr),
42 fBaryonsFile(nullptr)
43 , fMesonsFile(nullptr)
44 , fParticleTable()
45 , fFileName(nullptr)
46{
47}
48// ------------------------------------------------------------------------
49
50
51// ----- Standard constructor -----------------------------------------
52CbmPhsdGenerator::CbmPhsdGenerator(const char* fileNameInput, const char* fileNameBaryons, const char* fileNameMesons)
53 : FairGenerator()
54 , fBaryonsFile(nullptr)
55 , fMesonsFile(nullptr)
56 , fDatFile(nullptr)
57 , fParticleTable()
58{
59 // cout << "-I CbmPhsdGenerator: Opening input file " << fileNameInput << endl;
60 cout << "-I CbmPhsdGenerator: Opening HSD Baryons file " << fileNameBaryons << endl;
61 fBaryonsFile = fopen(fileNameBaryons, "r");
62 if (!fBaryonsFile) { Fatal("FairHSDgenerator", "Cannot open HSD Baryons file."); }
63 cout << "-I CbmPhsdGenerator: Opening HSD Mesons file " << fileNameMesons << endl;
64 fMesonsFile = fopen(fileNameMesons, "r");
65 if (!fMesonsFile) { Fatal("FairHSDgenerator", "Cannot open HSD Mesons file."); }
66 ReadConversionTable();
67 ReadCollisionData(fileNameInput);
68 nextBaryon.init = false;
69 nextMeson.init = false;
70 nextEvent = 1;
71 fReadDat = false;
72}
73// ------------------------------------------------------------------------
74
75
76// ----- Standard constructor -----------------------------------------
77CbmPhsdGenerator::CbmPhsdGenerator(const char* fileNameInput, const char* fileNameDat)
78 : FairGenerator()
79 , fBaryonsFile(nullptr)
80 , fMesonsFile(nullptr)
81 , fDatFile(nullptr)
82 , fParticleTable()
83{
84 // cout << "-I CbmPhsdGenerator: Opening input file " << fileNameInput << endl;
85 cout << "-I CbmPhsdGenerator: Opening phsd.dat file " << fileNameDat << endl;
86 fDatFile = fopen(fileNameDat, "r");
87 if (!fDatFile) { Fatal("CbmPhsdGenerator", "Cannot open phsd.dat file."); }
88 ReadCollisionData(fileNameInput);
89 nextEvent = 1;
90 fReadDat = true;
91}
92// ------------------------------------------------------------------------
93
94
95// ----- Destructor ---------------------------------------------------
96CbmPhsdGenerator::~CbmPhsdGenerator()
97{
98 // cout<<"Enter Destructor of CbmPhsdGenerator"<<endl;
99 /*if ( fInputFile ) {
100 fclose(fInputFile);
101 fInputFile = nullptr;
102 }*/
103 if (fBaryonsFile) {
104 fclose(fBaryonsFile);
105 fBaryonsFile = nullptr;
106 }
107 if (fMesonsFile) {
108 fclose(fMesonsFile);
109 fMesonsFile = nullptr;
110 }
111 if (fDatFile) {
112 fclose(fDatFile);
113 fDatFile = nullptr;
114 }
115 fParticleTable.clear();
116 // cout<<"Leave Destructor of CbmPhsdGenerator"<<endl;
117}
118// ------------------------------------------------------------------------
119
120
121// ----- Public method ReadEvent --------------------------------------
122Bool_t CbmPhsdGenerator::ReadEvent(FairPrimaryGenerator* primGen)
123{
124 if (fReadDat) return ReadEventDat(primGen);
125 else
126 return ReadEvent300(primGen);
127}
128// ------------------------------------------------------------------------
129
130
131// ----- Public method ReadEventDat -----------------------------------
132Bool_t CbmPhsdGenerator::ReadEventDat(FairPrimaryGenerator* primGen)
133{
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, "%d%d%lf%lf%lf%lf%*[^\n]%*c", &ttype, &tchr, &tmppx, &tmppy, &tmppz, &tmpp0);
170
171 Int_t pdgID = ttype;
172
173 // cout << "-I CbmPhsdGenerator: PID " << pdgID << endl;
174
175 Double_t eBeam = ekin + kProtonMass;
176 Double_t pBeam = TMath::Sqrt(eBeam * eBeam - kProtonMass * kProtonMass);
177 Double_t betaCM = pBeam / (eBeam + kProtonMass);
178 Double_t gammaCM = TMath::Sqrt(1. / (1. - betaCM * betaCM));
179
180 // Lorentztransformation to lab
181
182 Double_t px = Double_t(tmppx);
183 Double_t py = Double_t(tmppy);
184 Double_t pz = Double_t(tmppz);
185 Double_t e = Double_t(tmpp0);
186 Double_t mass = sqrt(e * e - px * px - py * py - pz * pz);
187 pz = gammaCM * (pz + betaCM * e);
188 Double_t ee = sqrt(mass * mass + px * px + py * py + pz * pz);
189
190 TVector3 aa(px, py, pz);
191 TLorentzVector pp;
192 pp.SetPx(px);
193 pp.SetPy(py);
194 pp.SetPz(pz);
195 pp.SetE(ee);
196
197 // Give track to PrimaryGenerator
198 primGen->AddTrack(pdgID, px, py, pz, 0., 0., 0.);
199 }
200
201 cout << "-I CbmPhsdGenerator: Event " << nextEvent << ", b = " << b << " fm, multiplicity " << nTracks
202 << ", ekin: " << ekin << endl;
203
204 nextEvent++;
205
206 return kTRUE;
207}
208// ------------------------------------------------------------------------
209
210
211// ----- Public method ReadEvent300 -----------------------------------
212Bool_t CbmPhsdGenerator::ReadEvent300(FairPrimaryGenerator* primGen)
213{
214
215 // ---> Check for input file
216 if (!fBaryonsFile) {
217 cout << "-E CbmPhsdGenerator: Baryons input file is not open! " << endl;
218 return kFALSE;
219 }
220
221 if (!fMesonsFile) {
222 cout << "-E CbmPhsdGenerator: Mesons input file is not open! " << endl;
223 return kFALSE;
224 }
225
226 // ---> Check for primary generator
227 if (!primGen) {
228 cout << "-E- CbmPhsdGenerator::ReadEvent: "
229 << "No PrimaryGenerator!" << endl;
230 return kFALSE;
231 }
232
233 if (!nextBaryon.init && feof(fBaryonsFile)) {
234 cout << "-E- CbmPhsdGenerator::ReadEvent: "
235 << "No more events!" << endl;
236 return kFALSE;
237 }
238
239 /*if (!nextMeson.init && feof(fMesonsFile)) {
240 cout << "-E- CbmPhsdGenerator::ReadEvent: "
241 << "No more events!" << endl;
242 return kFALSE;
243 }*/
244
245 if (!nextBaryon.init) {
246 fscanf(fBaryonsFile, "%d %d %d %d %lf %lf %lf %lf %lf", &nextBaryon.id, &nextBaryon.charge, &nextBaryon.ISUB,
247 &nextBaryon.IRUN, &nextBaryon.px, &nextBaryon.py, &nextBaryon.pz, &nextBaryon.p0, &nextBaryon.b);
248 fscanf(fBaryonsFile, "%*[^\n]%*c");
249 nextBaryon.globalEvent = (nextBaryon.ISUB - 1) * IRUNS + nextBaryon.IRUN;
250 nextBaryon.init = 1;
251 }
252
253 if (!nextMeson.init && !feof(fMesonsFile)) {
254 fscanf(fMesonsFile, "%d %d %d %d %lf %lf %lf %lf %lf", &nextMeson.id, &nextMeson.charge, &nextMeson.ISUB,
255 &nextMeson.IRUN, &nextMeson.px, &nextMeson.py, &nextMeson.pz, &nextMeson.p0, &nextMeson.b);
256 fscanf(fMesonsFile, "%*[^\n]%*c");
257 nextMeson.globalEvent = (nextMeson.ISUB - 1) * IRUNS + nextMeson.IRUN;
258 nextMeson.init = 1;
259 }
260
261
262 int nTracks = 0, pid = 0;
263 double b = nextBaryon.b;
264
265 // Set event id and impact parameter in MCEvent if not yet done
266 FairMCEventHeader* event = primGen->GetEvent();
267 if (event && (!event->IsSet())) {
268 event->SetEventID(nextEvent);
269 event->SetB(b);
270 event->MarkSet(kTRUE);
271 }
272
273 // Read all baryons from current event until we reach the next event in fort.300 or end of file
274 while (nextBaryon.globalEvent == nextEvent) {
275 // Convert HSD type and charge to unique pid identifier which is based on
276 // HSD particle id and charge, calculated separately for baryons, anti-baryons and mesons
277 if (nextBaryon.id >= 0) { pid = nextBaryon.id * 10 + (2 + nextBaryon.charge); }
278 // antibaryons
279 else {
280 pid = -(-nextBaryon.id * 10 + (2 - nextBaryon.charge));
281 }
282
283 // Convert Unique PID into PDG particle code
284 if (fParticleTable.find(pid) == fParticleTable.end()) {
285 cout << "-W CbmPhsdGenerator: PID " << nextBaryon.id << " charge " << nextBaryon.charge << " not found in table ("
286 << pid << ")" << endl;
287 if (feof(fBaryonsFile)) {
288 nextBaryon.init = 0;
289 break;
290 }
291 fscanf(fBaryonsFile, "%d %d %d %d %lf %lf %lf %lf %lf", &nextBaryon.id, &nextBaryon.charge, &nextBaryon.ISUB,
292 &nextBaryon.IRUN, &nextBaryon.px, &nextBaryon.py, &nextBaryon.pz, &nextBaryon.p0, &nextBaryon.b);
293 fscanf(fBaryonsFile, "%*[^\n]%*c");
294 nextBaryon.globalEvent = (nextBaryon.ISUB - 1) * IRUNS + nextBaryon.IRUN;
295 continue;
296 }
297 Int_t pdgID = fParticleTable[pid];
298
299 Double_t eBeam = ekin + kProtonMass;
300 Double_t pBeam = TMath::Sqrt(eBeam * eBeam - kProtonMass * kProtonMass);
301 Double_t betaCM = pBeam / (eBeam + kProtonMass);
302 Double_t gammaCM = TMath::Sqrt(1. / (1. - betaCM * betaCM));
303
304 // Lorentztransformation to lab
305
306 Double_t px = Double_t(nextBaryon.px);
307 Double_t py = Double_t(nextBaryon.py);
308 Double_t pz = Double_t(nextBaryon.pz);
309 Double_t e = Double_t(nextBaryon.p0);
310 Double_t mass = sqrt(e * e - px * px - py * py - pz * pz);
311 pz = gammaCM * (pz + betaCM * e);
312 Double_t ee = sqrt(mass * mass + px * px + py * py + pz * pz);
313
314 TVector3 aa(px, py, pz);
315 TLorentzVector pp;
316 pp.SetPx(px);
317 pp.SetPy(py);
318 pp.SetPz(pz);
319 pp.SetE(ee);
320
321 // Give track to PrimaryGenerator
322 primGen->AddTrack(pdgID, px, py, pz, 0., 0., 0.);
323
324 nTracks++;
325
326 if (feof(fBaryonsFile)) {
327 nextBaryon.init = 0;
328 break;
329 }
330
331 fscanf(fBaryonsFile, "%d %d %d %d %lf %lf %lf %lf %lf", &nextBaryon.id, &nextBaryon.charge, &nextBaryon.ISUB,
332 &nextBaryon.IRUN, &nextBaryon.px, &nextBaryon.py, &nextBaryon.pz, &nextBaryon.p0, &nextBaryon.b);
333 fscanf(fBaryonsFile, "%*[^\n]%*c");
334 nextBaryon.globalEvent = (nextBaryon.ISUB - 1) * IRUNS + nextBaryon.IRUN;
335 }
336
337 // Read all mesons from current event until we reach the next event in fort.301 or end of file
338 while (nextMeson.globalEvent == nextEvent) {
339 // Convert HSD type and charge to unique pid identifier which is based on
340 // HSD particle id and charge, calculated separately for baryons, anti-baryons and mesons
341 {
342 pid = 1000 + nextMeson.id * 10 + (2 + nextMeson.charge);
343 }
344
345 // Convert Unique PID into PDG particle code
346 if (fParticleTable.find(pid) == fParticleTable.end()) {
347 cout << "-W CbmPhsdGenerator: PID " << nextMeson.id << " charge " << nextMeson.charge << " not found in table ("
348 << pid << ")" << endl;
349 fscanf(fMesonsFile, "%d %d %d %d %lf %lf %lf %lf %lf", &nextMeson.id, &nextMeson.charge, &nextMeson.ISUB,
350 &nextMeson.IRUN, &nextMeson.px, &nextMeson.py, &nextMeson.pz, &nextMeson.p0, &nextMeson.b);
351 fscanf(fMesonsFile, "%*[^\n]%*c");
352 nextMeson.globalEvent = (nextMeson.ISUB - 1) * IRUNS + nextMeson.IRUN;
353 continue;
354 }
355 Int_t pdgID = fParticleTable[pid];
356
357 Double_t eBeam = ekin + kProtonMass;
358 Double_t pBeam = TMath::Sqrt(eBeam * eBeam - kProtonMass * kProtonMass);
359 Double_t betaCM = pBeam / (eBeam + kProtonMass);
360 Double_t gammaCM = TMath::Sqrt(1. / (1. - betaCM * betaCM));
361
362 // Lorentztransformation to lab
363
364 Double_t px = Double_t(nextMeson.px);
365 Double_t py = Double_t(nextMeson.py);
366 Double_t pz = Double_t(nextMeson.pz);
367 Double_t e = Double_t(nextMeson.p0);
368 Double_t mass = sqrt(e * e - px * px - py * py - pz * pz);
369 pz = gammaCM * (pz + betaCM * e);
370 Double_t ee = sqrt(mass * mass + px * px + py * py + pz * pz);
371
372 TVector3 aa(px, py, pz);
373 TLorentzVector pp;
374 pp.SetPx(px);
375 pp.SetPy(py);
376 pp.SetPz(pz);
377 pp.SetE(ee);
378
379 // Give track to PrimaryGenerator
380 primGen->AddTrack(pdgID, px, py, pz, 0., 0., 0.);
381
382 nTracks++;
383
384 if (feof(fMesonsFile)) {
385 nextMeson.init = 0;
386 break;
387 }
388
389 fscanf(fMesonsFile, "%d %d %d %d %lf %lf %lf %lf %lf", &nextMeson.id, &nextMeson.charge, &nextMeson.ISUB,
390 &nextMeson.IRUN, &nextMeson.px, &nextMeson.py, &nextMeson.pz, &nextMeson.p0, &nextMeson.b);
391 fscanf(fMesonsFile, "%*[^\n]%*c");
392 nextMeson.globalEvent = (nextMeson.ISUB - 1) * IRUNS + nextMeson.IRUN;
393 }
394
395 cout << "-I CbmPhsdGenerator: Event " << nextEvent << ", b = " << b << " fm, multiplicity " << nTracks
396 << ", ekin: " << ekin << endl;
397
398 nextEvent++;
399
400 return kTRUE;
401}
402// ------------------------------------------------------------------------
403
470// ------------------------------------------------------------------------
471
472// ----- Private method ReadConversionTable ---------------------------
473void CbmPhsdGenerator::ReadConversionTable()
474{
475
476 TString work = getenv("VMCWORKDIR");
477 TString fileName = work + "/input/hsd_pdg.dat";
478 std::ifstream* pdgconv = new std::ifstream(fileName.Data());
479
480 if (!(*pdgconv).is_open()) {
481 Fatal("CbmPhsdGenerator", "Particle table for conversion was not found!");
482 //cout << "-W CbmPhsdGenerator: Particle table for conversion was not found!" << endl;
483 }
484
485 Int_t index = 0;
486 Int_t pdgId = 0;
487
488 string tmpStr;
489
490 while (!pdgconv->eof()) {
491 index = pdgId = 0;
492 *pdgconv >> index >> pdgId;
493 std::getline(*pdgconv, tmpStr);
494 fParticleTable[index] = pdgId;
495 }
496
497 pdgconv->close();
498 delete pdgconv;
499
500 cout << "-I CbmPhsdGenerator: Particle table for conversion from "
501 << "HSD loaded" << endl;
502}
503// ------------------------------------------------------------------------
504
505// ----- Private method ReadCollisionData ----------------------------
506void CbmPhsdGenerator::ReadCollisionData(const char* fileNameInput)
507{
508
509 std::ifstream* inputf = new std::ifstream(fileNameInput);
510
511 // Int_t index = 0;
512 // Int_t pdgId = 0;
513
514 //TString val, name, tname;
515 string val, name, tname, tmpStr;
516
517 while (!inputf->eof()) {
518 // index =pdgId =0 ;
519 //*inputf >> index >> pdgId ;
520 *inputf >> val >> name;
521 tname = name.substr(0, 6);
522 //cout << tname << endl;
523 if (tname == "MASSTA") At = atoi(val.substr(0, val.size() - 1).c_str());
524 else if (tname == "MSTAPR")
525 Zt = atoi(val.substr(0, val.size() - 1).c_str());
526 else if (tname == "MASSPR")
527 Ap = atoi(val.substr(0, val.size() - 1).c_str());
528 else if (tname == "MSPRPR")
529 Zp = atoi(val.substr(0, val.size() - 1).c_str());
530 tname = name.substr(0, 4);
531 if (tname == "ELAB") ekin = atof(val.substr(0, val.size() - 1).c_str());
532 tname = name.substr(0, 5);
533 if (tname == "ISUBS") ISUBS = atof(val.substr(0, val.size() - 1).c_str());
534 tname = name.substr(0, 3);
535 if (tname == "NUM") IRUNS = atof(val.substr(0, val.size() - 1).c_str());
536 std::getline(*inputf, tmpStr);
537 }
538
539 inputf->close();
540 delete inputf;
541
542 cout << "-I CbmPhsdGenerator: Collision data from "
543 << "HSD loaded" << endl;
544}
545// ------------------------------------------------------------------------
546
547
548ClassImp(CbmPhsdGenerator);
kProtonMass
const Double_t kProtonMass
Definition CbmPhsdGenerator.cxx:34
ClassImp
ClassImp(CbmPhsdGenerator)
sqrt
friend fvec sqrt(const fvec &a)
Definition KfSimdPseudo.h:118
CbmPhsdGenerator::fBaryonsFile
FILE * fBaryonsFile
Whether phsd.dat or .300/301 files are used.
Definition CbmPhsdGenerator.h:94
CbmPhsdGenerator::nextBaryon
Hadron nextBaryon
Input file name.
Definition CbmPhsdGenerator.h:103
CbmPhsdGenerator::ReadEvent300
Bool_t ReadEvent300(FairPrimaryGenerator *primGen)
Definition CbmPhsdGenerator.cxx:212
CbmPhsdGenerator::ReadCollisionData
void ReadCollisionData(const char *fileNameInput)
Definition CbmPhsdGenerator.cxx:506
CbmPhsdGenerator::fParticleTable
std::map< Int_t, Int_t > fParticleTable
HSD output files.
Definition CbmPhsdGenerator.h:98
CbmPhsdGenerator::fDatFile
FILE * fDatFile
HSD output files.
Definition CbmPhsdGenerator.h:96
CbmPhsdGenerator::ReadEventDat
Bool_t ReadEventDat(FairPrimaryGenerator *primGen)
Definition CbmPhsdGenerator.cxx:132
CbmPhsdGenerator::ReadEvent
Bool_t ReadEvent(FairPrimaryGenerator *primGen)
Definition CbmPhsdGenerator.cxx:122