FairRoot
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
FairUrqmdGenerator.cxx
Go to the documentation of this file.
1 /********************************************************************************
2  * Copyright (C) 2014 GSI Helmholtzzentrum fuer Schwerionenforschung GmbH *
3  * *
4  * This software is distributed under the terms of the *
5  * GNU Lesser General Public Licence (LGPL) version 3, *
6  * copied verbatim in the file "LICENSE" *
7  ********************************************************************************/
8 // -------------------------------------------------------------------------
9 // ----- FairUrqmdGenerator source file -----
10 // ----- Created 24/06/04 by V. Friese -----
11 // -------------------------------------------------------------------------
12 #include "FairUrqmdGenerator.h"
13 
14 #include "FairLogger.h" // for logging
15 #include "FairMCEventHeader.h" // for FairMCEventHeader
16 #include "FairPrimaryGenerator.h" // for FairPrimaryGenerator
17 
18 #include <TDatabasePDG.h> // for TDatabasePDG
19 #include <TLorentzVector.h> // for TLorentzVector
20 #include <TMath.h> // for Sqrt, sqrt
21 #include <TParticlePDG.h> // for TParticlePDG
22 #include <TString.h> // for TString, operator+
23 #include <TVector3.h> // for TVector3
24 #include <climits> // for INT_MAX
25 #include <cmath> // sqrt
26 #include <cstdlib> // for getenv
27 #include <fstream> // IWYU pragma: keep for ifstream
28 
29 FairUrqmdGenerator::FairUrqmdGenerator()
30  : FairGenerator()
31  , fInputFile(nullptr)
32  , fParticleTable()
33  , fFileName(nullptr)
34 {}
35 
36 FairUrqmdGenerator::FairUrqmdGenerator(const char* fileName)
37  : FairGenerator()
38  , fInputFile(nullptr)
39  , fParticleTable()
40  , fFileName(fileName)
41 {
42  // fFileName = fileName;
43  LOG(info) << "FairUrqmdGenerator: Opening input file " << fileName;
44  fInputFile = fopen(fFileName, "r");
45  if (!fInputFile) {
46  LOG(fatal) << "Cannot open input file.";
47  }
48  ReadConversionTable();
49 }
50 
51 FairUrqmdGenerator::FairUrqmdGenerator(const char* fileName, const char* conversion_table)
52  : FairGenerator()
53  , fInputFile(nullptr)
54  , fParticleTable()
55  , fFileName(fileName)
56 {
57  // fFileName = fileName;
58  LOG(info) << "FairUrqmdGenerator: Opening input file " << fileName;
59  fInputFile = fopen(fFileName, "r");
60  if (!fInputFile) {
61  LOG(fatal) << "Cannot open input file.";
62  }
63  ReadConversionTable(conversion_table);
64 }
65 
66 FairUrqmdGenerator::~FairUrqmdGenerator()
67 {
68  // LOG(debug) << "Enter Destructor of FairUrqmdGenerator";
69  if (fInputFile) {
70  fclose(fInputFile);
71  fInputFile = nullptr;
72  }
73  fParticleTable.clear();
74  // LOG(debug) << "Leave Destructor of FairUrqmdGenerator";
75 }
76 
77 namespace {
78 // silence fgets warnings (with a runtime error message if it actually happens)
79 // for correctness should be handled in the code instead, if this methods are still relevant in the future.
80 void skipBytes(char* bytes, size_t numBytes, FILE* file)
81 {
82  if (fgets(bytes, numBytes, file) == nullptr) {
83  LOG(error) << "Failed reading from file stream";
84  }
85 }
86 
87 } // namespace
88 
89 Bool_t FairUrqmdGenerator::ReadEvent(FairPrimaryGenerator* primGen)
90 {
91  // ---> Check for input file
92  if (!fInputFile) {
93  LOG(error) << "FairUrqmdGenerator: Input file not open! ";
94  return kFALSE;
95  }
96 
97  // ---> Check for primary generator
98  if (!primGen) {
99  LOG(error) << "FairUrqmdGenerator::ReadEvent: "
100  << "No PrimaryGenerator!";
101  return kFALSE;
102  }
103 
104  // ---> Define event variables to be read from file
105  int evnr = 0, ntracks = 0, aProj = 0, zProj = 0, aTarg = 0, zTarg = 0;
106  float b = 0., ekin = 0.;
107 
108  int ityp = 0, i3 = 0, ichg = 0, pid = 0;
109  float ppx = 0., ppy = 0., ppz = 0., m = 0.;
110 
111  // ---> Read and check first event header line from input file
112  char read[200];
113  skipBytes(read, 200, fInputFile);
114  if (feof(fInputFile)) {
115  LOG(info) << "FairUrqmdGenerator : End of input file reached.";
116  fclose(fInputFile);
117  fInputFile = nullptr;
118  return kFALSE;
119  }
120  if (read[0] != 'U') {
121  LOG(error) << "FairUrqmdGenerator: Wrong event header";
122  return kFALSE;
123  }
124 
125  Int_t retval = 0;
126 
127  // ---> Read rest of event header
128  skipBytes(read, 26, fInputFile);
129  retval = fscanf(fInputFile, "%d", &aProj);
130  CheckReturnValue(retval);
131  retval = fscanf(fInputFile, "%d", &zProj);
132  CheckReturnValue(retval);
133  skipBytes(read, 25, fInputFile);
134  retval = fscanf(fInputFile, "%d", &aTarg);
135  CheckReturnValue(retval);
136  retval = fscanf(fInputFile, "%d", &zTarg);
137  CheckReturnValue(retval);
138  skipBytes(read, 200, fInputFile);
139  skipBytes(read, 200, fInputFile);
140  skipBytes(read, 36, fInputFile);
141  retval = fscanf(fInputFile, "%f", &b);
142  CheckReturnValue(retval);
143  skipBytes(read, 200, fInputFile);
144  skipBytes(read, 39, fInputFile);
145  retval = fscanf(fInputFile, "%e", &ekin);
146  CheckReturnValue(retval);
147  skipBytes(read, 200, fInputFile);
148  skipBytes(read, 7, fInputFile);
149  retval = fscanf(fInputFile, "%d", &evnr);
150  CheckReturnValue(retval);
151  skipBytes(read, 200, fInputFile);
152  for (int iline = 0; iline < 8; iline++) {
153  skipBytes(read, 200, fInputFile);
154  }
155  retval = fscanf(fInputFile, "%d", &ntracks);
156  if (ntracks < 0 || ntracks > (INT_MAX - 1))
157  LOG(fatal) << "Error reading the number of events from event header.";
158  CheckReturnValue(retval);
159  skipBytes(read, 200, fInputFile);
160  skipBytes(read, 200, fInputFile);
161 
162  // ---> Calculate beta and gamma for Lorentztransformation
163  TDatabasePDG* pdgDB = TDatabasePDG::Instance();
164  TParticlePDG* kProton = pdgDB->GetParticle(2212);
165  Double_t kProtonMass = kProton->Mass();
166 
167  Double_t eBeam = ekin + kProtonMass;
168  Double_t pBeam = TMath::Sqrt(eBeam * eBeam - kProtonMass * kProtonMass);
169  Double_t betaCM = pBeam / (eBeam + kProtonMass);
170  Double_t gammaCM = TMath::Sqrt(1. / (1. - betaCM * betaCM));
171 
172  LOG(info) << "FairUrqmdGenerator: Event " << evnr << ", b = " << b << " fm, multiplicity " << ntracks
173  << ", ekin: " << ekin;
174 
175  // Set event id and impact parameter in MCEvent if not yet done
176  FairMCEventHeader* event = primGen->GetEvent();
177  if (event && (!event->IsSet())) {
178  event->SetEventID(evnr);
179  event->SetB(b);
180  event->MarkSet(kTRUE);
181  }
182 
183  // ---> Loop over tracks in the current event
184  for (int itrack = 0; itrack < ntracks; itrack++) {
185 
186  // Read momentum and PID from file
187  skipBytes(read, 81, fInputFile);
188  retval = fscanf(fInputFile, "%e", &ppx);
189  CheckReturnValue(retval);
190  retval = fscanf(fInputFile, "%e", &ppy);
191  CheckReturnValue(retval);
192  retval = fscanf(fInputFile, "%e", &ppz);
193  CheckReturnValue(retval);
194  retval = fscanf(fInputFile, "%e", &m);
195  CheckReturnValue(retval);
196  retval = fscanf(fInputFile, "%d", &ityp);
197  CheckReturnValue(retval);
198  retval = fscanf(fInputFile, "%d", &i3);
199  CheckReturnValue(retval);
200  retval = fscanf(fInputFile, "%d", &ichg);
201  CheckReturnValue(retval);
202  skipBytes(read, 200, fInputFile);
203 
204  // Convert UrQMD type and charge to unique pid identifier
205  if (ityp >= 0) {
206  pid = 1000 * (ichg + 2) + ityp;
207  } else {
208  pid = -1000 * (ichg + 2) + ityp;
209  }
210 
211  // Convert Unique PID into PDG particle code
212  if (fParticleTable.find(pid) == fParticleTable.end()) {
213  LOG(warn) << "FairUrqmdGenerator: PID " << ityp << " charge " << ichg << " not found in table (" << pid
214  << ")";
215  continue;
216  }
217  Int_t pdgID = fParticleTable[pid];
218 
219  // Lorentztransformation to lab
220  Double_t mass = Double_t(m);
221  Double_t px = Double_t(ppx);
222  Double_t py = Double_t(ppy);
223  Double_t pz = Double_t(ppz);
224  Double_t e = sqrt(mass * mass + px * px + py * py + pz * pz);
225  pz = gammaCM * (pz + betaCM * e);
226  Double_t ee = sqrt(mass * mass + px * px + py * py + pz * pz);
227 
228  TVector3 aa(px, py, pz);
229  TLorentzVector pp;
230  pp.SetPx(px);
231  pp.SetPy(py);
232  pp.SetPz(pz);
233  pp.SetE(ee);
234 
235  // Give track to PrimaryGenerator
236  primGen->AddTrack(pdgID, px, py, pz, 0., 0., 0.);
237  }
238 
239  return kTRUE;
240 }
241 
242 Bool_t FairUrqmdGenerator::SkipEvents(Int_t count)
243 {
244  if (count <= 0) {
245  return kTRUE;
246  }
247 
248  for (Int_t ii = 0; ii < count; ii++) {
249  // ---> Check for input file
250  if (!fInputFile) {
251  LOG(error) << "FairUrqmdGenerator: Input file not open! ";
252  return kFALSE;
253  }
254 
255  // ---> Define event variables to be read from file
256  int evnr = 0, ntracks = 0, aProj = 0, zProj = 0, aTarg = 0, zTarg = 0;
257  float b = 0., ekin = 0.;
258 
259  // ---> Read and check first event header line from input file
260  char read[200];
261  skipBytes(read, 200, fInputFile);
262  if (feof(fInputFile)) {
263  LOG(info) << "FairUrqmdGenerator : End of input file reached.";
264  fclose(fInputFile);
265  fInputFile = nullptr;
266  return kFALSE;
267  }
268  if (read[0] != 'U') {
269  LOG(error) << "FairUrqmdGenerator: Wrong event header";
270  return kFALSE;
271  }
272 
273  Int_t retval = 0;
274 
275  // ---> Read rest of event header
276  skipBytes(read, 26, fInputFile);
277  retval = fscanf(fInputFile, "%d", &aProj);
278  CheckReturnValue(retval);
279  retval = fscanf(fInputFile, "%d", &zProj);
280  CheckReturnValue(retval);
281  skipBytes(read, 25, fInputFile);
282  retval = fscanf(fInputFile, "%d", &aTarg);
283  CheckReturnValue(retval);
284  retval = fscanf(fInputFile, "%d", &zTarg);
285  CheckReturnValue(retval);
286  skipBytes(read, 200, fInputFile);
287  skipBytes(read, 200, fInputFile);
288  skipBytes(read, 36, fInputFile);
289  retval = fscanf(fInputFile, "%f", &b);
290  CheckReturnValue(retval);
291  skipBytes(read, 200, fInputFile);
292  skipBytes(read, 39, fInputFile);
293  retval = fscanf(fInputFile, "%e", &ekin);
294  CheckReturnValue(retval);
295  skipBytes(read, 200, fInputFile);
296  skipBytes(read, 7, fInputFile);
297  retval = fscanf(fInputFile, "%d", &evnr);
298  CheckReturnValue(retval);
299  skipBytes(read, 200, fInputFile);
300  for (int iline = 0; iline < 8; iline++) {
301  skipBytes(read, 200, fInputFile);
302  }
303  retval = fscanf(fInputFile, "%d", &ntracks);
304  if (ntracks < 0 || ntracks > (INT_MAX - 1))
305  LOG(fatal) << "Error reading the number of events from event header.";
306  CheckReturnValue(retval);
307  skipBytes(read, 200, fInputFile);
308  skipBytes(read, 200, fInputFile);
309 
310  LOG(info) << "FairUrqmdGenerator: Event " << evnr << " skipped!";
311 
312  // ---> Loop over tracks in the current event
313  for (int itrack = 0; itrack < ntracks; itrack++) {
314 
315  // Read momentum and PID from file
316  skipBytes(read, 81, fInputFile);
317  skipBytes(read, 200, fInputFile);
318  }
319  }
320  return kTRUE;
321 }
322 
323 void FairUrqmdGenerator::ReadConversionTable(TString conversion_table)
324 {
325  TString work = getenv("VMCWORKDIR");
326  TString fileName;
327 
328  if (conversion_table.IsNull()) {
329  fileName = work + "/input/urqmd_pdg.dat";
330  } else {
331  fileName = conversion_table.Data();
332  }
333 
334  std::ifstream* pdgconv = new std::ifstream(fileName.Data());
335 
336  if (!pdgconv->good()) {
337  LOG(fatal) << "Could not open Urqmd->PDG input file " << fileName;
338  }
339 
340  Int_t index = 0;
341 
342  Int_t pdgId = 0;
343 
344  while (!pdgconv->eof()) {
345  index = pdgId = 0;
346  *pdgconv >> index >> pdgId;
347  fParticleTable[index] = pdgId;
348  }
349 
350  pdgconv->close();
351  delete pdgconv;
352 
353  LOG(info) << "FairUrqmdGenerator: Particle table for conversion from "
354  << "UrQMD loaded";
355 }
356 
357 void FairUrqmdGenerator::CheckReturnValue(Int_t retval)
358 {
359  if (1 != retval) {
360  LOG(error) << "Error when reading variable from input file";
361  }
362 }
363 
364 ClassImp(FairUrqmdGenerator);
FairUrqmdGenerator::SkipEvents
Bool_t SkipEvents(Int_t count)
Definition: FairUrqmdGenerator.cxx:242
ClassImp
ClassImp(FairEventBuilder)
FairPrimaryGenerator::AddTrack
virtual void AddTrack(Int_t pdgid, Double_t px, Double_t py, Double_t pz, Double_t vx, Double_t vy, Double_t vz, Int_t parent=-1, Bool_t wanttracking=true, Double_t e=-9e9, Double_t tof=0., Double_t weight=0., TMCProcess proc=kPPrimary)
Definition: FairPrimaryGenerator.cxx:279
FairPrimaryGenerator::GetEvent
FairMCEventHeader * GetEvent()
Definition: FairPrimaryGenerator.h:150
FairMCEventHeader::SetEventID
void SetEventID(UInt_t eventId)
Definition: FairMCEventHeader.h:60
FairUrqmdGenerator::ReadEvent
Bool_t ReadEvent(FairPrimaryGenerator *primGen)
Definition: FairUrqmdGenerator.cxx:89