BmnRoot
Loading...
Searching...
No Matches
MpdUrqmdGenerator.cxx
Go to the documentation of this file.
1// -------------------------------------------------------------------------
2// ----- MpdUrqmdGenerator source file -----
3// ----- Created 24/06/04 by V. Friese -----
4// -------------------------------------------------------------------------
5#include "MpdUrqmdGenerator.h"
6
7#include "FairPrimaryGenerator.h"
8#include "FairMCEventHeader.h"
9#include "constants.h"
10
11#include "TMCProcess.h"
12#include "TObjArray.h"
13#include "TPDGCode.h"
14#include "TParticle.h"
15#include "TRandom.h"
16#include "TString.h"
17#include "TVirtualMCStack.h"
18#include "TLorentzVector.h"
19#include "TDatabasePDG.h"
20#include "TParticlePDG.h"
21
22
23#include <iostream>
24#include <cstring>
25
26#include <stdio.h>
27
28using namespace std;
29using namespace TMath;
30
31//const Double_t kProtonMass = 0.938271998;
32
33
34// ----- Default constructor ------------------------------------------
35MpdUrqmdGenerator::MpdUrqmdGenerator()
36: FairGenerator(),
37 fInputFile(NULL),
38 fParticleTable(),
39 fPhiMin(0.),
40 fPhiMax(0.),
41 fEventPlaneSet(kFALSE),
42 fFileName(NULL),
43 fX(0.), fY(0.), fZ(0.)
44{}
45
46
47// ----- Standard constructor -----------------------------------------
48MpdUrqmdGenerator::MpdUrqmdGenerator(const char* fileName)
49: FairGenerator(),
50 fInputFile(NULL),
51 fParticleTable(),
52 fPhiMin(0.),
53 fPhiMax(0.),
54 fEventPlaneSet(kFALSE),
55 fFileName(fileName),
56 fX(0.), fY(0.), fZ(0.)
57{
58 // fFileName = fileName;
59 cout << "-I MpdUrqmdGenerator: Opening input file " << fFileName << endl;
60 fInputFile = gzopen(fFileName, "rb");
61 if (!fInputFile) {
62 Fatal("MpdUrqmdGenerator", "Cannot open input file.");
63 exit(1);
64 }
65 ReadConversionTable();
66}
67// ------------------------------------------------------------------------
68
69
70
71// ----- Destructor ---------------------------------------------------
72MpdUrqmdGenerator::~MpdUrqmdGenerator() {
73 // cout<<"Enter Destructor of MpdUrqmdGenerator"<<endl;
74 if (fInputFile) {
75#ifdef GZIP_SUPPORT
76 gzclose(fInputFile);
77#else
78 fclose(fInputFile);
79#endif
80 fInputFile = NULL;
81 }
82 fParticleTable.clear();
83 // cout<<"Leave Destructor of MpdUrqmdGenerator"<<endl;
84}
85// ------------------------------------------------------------------------
86
87// ----- Public method ReadEvent --------------------------------------
88
89Bool_t MpdUrqmdGenerator::ReadEvent(FairPrimaryGenerator* primGen) {
90
91 // ---> Check for input file
92 if (!fInputFile) {
93 cout << "-E MpdUrqmdGenerator: Input file not open! " << endl;
94 return kFALSE;
95 }
96
97 // ---> Check for primary generator
98 if (!primGen) {
99 cout << "-E- MpdUrqmdGenerator::ReadEvent: "
100 << "No PrimaryGenerator!" << endl;
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#ifdef GZIP_SUPPORT
114 gzgets(fInputFile, read, 200); // line 1
115#else
116 fgets(read, 200, fInputFile);
117#endif
118 Int_t urqmdVersion = 0;
119 sscanf(read, "UQMD version: %d 1000 %d output_file 14", &urqmdVersion, &urqmdVersion);
120 cout << "URQMD VERSION USED = " << urqmdVersion << endl;
121#ifdef GZIP_SUPPORT
122 if (gzeof(fInputFile)) {
123 cout << "-I MpdUrqmdGenerator : End of input file reached." << endl;
124 gzclose(fInputFile);
125#else
126 if ( feof(fInputFile) ) {
127 cout << "-I MpdUrqmdGenerator : End of input file reached." << endl;
128 fclose(fInputFile);
129#endif
130 fInputFile = NULL;
131 return kFALSE;
132 }
133 if (read[0] != 'U') {
134 cout << "-E MpdUrqmdGenerator: Wrong event header" << endl;
135 return kFALSE;
136 }
137
138 // ---> Read rest of event header
139#ifdef GZIP_SUPPORT
140 gzgets(fInputFile, read, 200); // line 2
141 sscanf(read, "projectile: (mass, char) %d %d target: (mass, char) %d %d",
142 &aProj, &zProj, &aTarg, &zTarg); // line 2
143 gzgets(fInputFile, read, 200); // line 3
144 gzgets(fInputFile, read, 36); // line 4
145 gzgets(fInputFile, read, 200); // line 4
146 sscanf(read, "%f", &b); // line 4
147 gzgets(fInputFile, read, 39); // line 5
148 gzgets(fInputFile, read, 200); // line 5
149 sscanf(read, "%e", &ekin); // line 5
150 gzgets(fInputFile, read, 7); // line 6
151 gzgets(fInputFile, read, 200); // line 6
152 sscanf(read, "%d", &evnr); // line 6
153
154 for (Int_t iline = 0; iline < ((urqmdVersion == 30400) ? 11 : 8); iline++)
155 gzgets(fInputFile, read, 200);
156
157 gzgets(fInputFile, read, 200); // line 18
158 sscanf(read, "%d", &ntracks); // line 18
159 gzgets(fInputFile, read, 200); // line 19
160#else
161 fgets(read, 26, fInputFile);
162 fscanf(fInputFile, "%d", &aProj);
163 fscanf(fInputFile, "%d", &zProj);
164 fgets(read, 25, fInputFile);
165 fscanf(fInputFile, "%d", &aTarg);
166 fscanf(fInputFile, "%d", &zTarg);
167 fgets(read, 200, fInputFile);
168 fgets(read, 200, fInputFile);
169 fgets(read, 36, fInputFile);
170 fscanf(fInputFile, "%f", &b);
171 fgets(read, 200, fInputFile);
172 fgets(read, 39, fInputFile);
173 fscanf(fInputFile, "%e", &ekin);
174 fgets(read, 200, fInputFile);
175 fgets(read, 7, fInputFile);
176 fscanf(fInputFile, "%d", &evnr);
177 fgets(read, 200, fInputFile);
178 for (int iline=0; iline<8; iline++) { fgets(read, 200,fInputFile); }
179 fscanf(fInputFile, "%d", &ntracks);
180 fgets(read, 200, fInputFile);
181 fgets(read, 200, fInputFile);
182#endif
183
184 // ---> Calculate beta and gamma for Lorentztransformation
185 TDatabasePDG* pdgDB = TDatabasePDG::Instance();
186 TParticlePDG* kProton = pdgDB->GetParticle(2212);
187 Double_t kProtonMass = kProton->Mass();
188
189 Double_t eBeam = ekin + kProtonMass;
190 Double_t pBeam = TMath::Sqrt(eBeam * eBeam - kProtonMass * kProtonMass);
191 Double_t betaCM = pBeam / (eBeam + kProtonMass);
192 Double_t gammaCM = TMath::Sqrt(1. / (1. - betaCM * betaCM));
193
194 cout << "-I MpdUrqmdGenerator: Event " << evnr << ", b = " << b
195 << " fm, multiplicity " << ntracks << ", ekin: " << ekin << endl;
196
197
198 Double_t phi = 0.;
199 // ---> Generate rotation angle
200 if (fEventPlaneSet) {
201 phi = gRandom->Uniform(fPhiMin, fPhiMax);
202 }
203
204 // Set event id and impact parameter in MCEvent if not yet done
205 FairMCEventHeader* event = primGen->GetEvent();
206 if (event && (!event->IsSet())) {
207 event->SetEventID(evnr);
208 event->SetB(b);
209 event->MarkSet(kTRUE);
210 event->SetRotZ(phi);
211 }
212
213
214 // ---> Loop over tracks in the current event
215 for (int itrack = 0; itrack < ntracks; itrack++) {
216
217 // Read momentum and PID from file
218#ifdef GZIP_SUPPORT
219 gzgets(fInputFile, read, 81);
220 gzgets(fInputFile, read, 200);
221 sscanf(read, "%e %e %e %e %d %d %d", &ppx, &ppy, &ppz, &m, &ityp, &i3, &ichg);
222#else
223 fgets(read, 81, fInputFile);
224 fscanf(fInputFile, "%e", &ppx);
225 fscanf(fInputFile, "%e", &ppy);
226 fscanf(fInputFile, "%e", &ppz);
227 fscanf(fInputFile, "%e", &m);
228 fscanf(fInputFile, "%d", &ityp);
229 fscanf(fInputFile, "%d", &i3);
230 fscanf(fInputFile, "%d", &ichg);
231 fgets(read, 200, fInputFile);
232#endif
233
234 // Convert UrQMD type and charge to unique pid identifier
235 if (ityp >= 0) {
236 pid = 1000 * (ichg + 2) + ityp;
237 } else {
238 pid = -1000 * (ichg + 2) + ityp;
239 }
240
241 // Convert Unique PID into PDG particle code
242 if (fParticleTable.find(pid) == fParticleTable.end()) {
243 cout << "-W MpdUrqmdGenerator: PID " << ityp << " charge "
244 << ichg << " not found in table (" << pid << ")" << endl;
245 continue;
246 }
247 Int_t pdgID = fParticleTable[pid];
248
249 // Lorentztransformation to lab
250 Double_t mass = Double_t(m);
251 Double_t px = Double_t(ppx);
252 Double_t py = Double_t(ppy);
253 Double_t pz = Double_t(ppz);
254 Double_t e = sqrt(mass * mass + px * px + py * py + pz * pz);
255 if (gCoordinateSystem == sysLaboratory)
256 pz = gammaCM * (pz + betaCM * e);
257 Double_t ee = sqrt(mass * mass + px * px + py * py + pz * pz);
258
259 if (fEventPlaneSet) {
260 Double_t pt = Sqrt(px * px + py * py);
261 Double_t azim = ATan2(py, px);
262 azim += phi;
263 px = pt * Cos(azim);
264 py = pt * Sin(azim);
265 }
266
267 TLorentzVector pp;
268 pp.SetPx(px);
269 pp.SetPy(py);
270 pp.SetPz(pz);
271 pp.SetE(ee);
272
273 // Give track to PrimaryGenerator
274 //primGen->AddTrack(pdgID, px, py, pz, 0., 0., 0.);
275 primGen->AddTrack(pdgID, px, py, pz, fX, fY, fZ);
276
277 }
278
279 return kTRUE;
280}
281// ------------------------------------------------------------------------
282
283
284// ----- Public method ReadEvent --------------------------------------
285
286Bool_t MpdUrqmdGenerator::SkipEvents(Int_t count) {
287 if (count <= 0) {
288 return kTRUE;
289 }
290
291 for (Int_t ii = 0; ii < count; ii++) {
292 // ---> Check for input file
293 if (!fInputFile) {
294 cout << "-E MpdUrqmdGenerator: Input file not open! " << endl;
295 return kFALSE;
296 }
297
298 // ---> Define event variables to be read from file
299 int evnr = 0, ntracks = 0, aProj = 0, zProj = 0, aTarg = 0, zTarg = 0;
300 float b = 0., ekin = 0.;
301
302 // ---> Read and check first event header line from input file
303 char read[200];
304#ifdef GZIP_SUPPORT
305 gzgets(fInputFile, read, 200); // line 1
306#else
307 fgets(read, 200, fInputFile);
308#endif
309 Int_t urqmdVersion = 0;
310 sscanf(read, "UQMD version: %d 1000 %d output_file 14", &urqmdVersion, &urqmdVersion);
311 cout << "URQMD VERSION USED = " << urqmdVersion << endl;
312
313#ifdef GZIP_SUPPORT
314 if (gzeof(fInputFile)) {
315 cout << "-I MpdUrqmdGenerator : End of input file reached." << endl;
316 gzclose(fInputFile);
317#else
318 if ( feof(fInputFile) ) {
319 cout << "-I MpdUrqmdGenerator : End of input file reached." << endl;
320 fclose(fInputFile);
321#endif
322 fInputFile = NULL;
323 return kFALSE;
324 }
325 if (read[0] != 'U') {
326 cout << "-E MpdUrqmdGenerator: Wrong event header" << endl;
327 return kFALSE;
328 }
329
330 // ---> Read rest of event header
331#ifdef GZIP_SUPPORT
332 gzgets(fInputFile, read, 200); // line 2
333 sscanf(read, "projectile: (mass, char) %d %d target: (mass, char) %d %d",
334 &aProj, &zProj, &aTarg, &zTarg); // line 2
335 gzgets(fInputFile, read, 200); // line 3
336 gzgets(fInputFile, read, 36); // line 4
337 gzgets(fInputFile, read, 200); // line 4
338 sscanf(read, "%f", &b); // line 4
339 gzgets(fInputFile, read, 39); // line 5
340 gzgets(fInputFile, read, 200); // line 5
341 sscanf(read, "%e", &ekin); // line 5
342 gzgets(fInputFile, read, 7); // line 6
343 gzgets(fInputFile, read, 200); // line 6
344 sscanf(read, "%d", &evnr); // line 6
345
346 for (Int_t iline = 0; iline < ((urqmdVersion == 30400) ? 11 : 8); iline++)
347 gzgets(fInputFile, read, 200);
348
349 gzgets(fInputFile, read, 200); // line 15
350 sscanf(read, "%d", &ntracks); // line 15
351 gzgets(fInputFile, read, 200); // line 16
352#else
353 fgets(read, 26, fInputFile);
354 fscanf(fInputFile, "%d", &aProj);
355 fscanf(fInputFile, "%d", &zProj);
356 fgets(read, 25, fInputFile);
357 fscanf(fInputFile, "%d", &aTarg);
358 fscanf(fInputFile, "%d", &zTarg);
359 fgets(read, 200, fInputFile);
360 fgets(read, 200, fInputFile);
361 fgets(read, 36, fInputFile);
362 fscanf(fInputFile, "%f", &b);
363 fgets(read, 200, fInputFile);
364 fgets(read, 39, fInputFile);
365 fscanf(fInputFile, "%e", &ekin);
366 fgets(read, 200, fInputFile);
367 fgets(read, 7, fInputFile);
368 fscanf(fInputFile, "%d", &evnr);
369 fgets(read, 200, fInputFile);
370 for (int iline=0; iline<8; iline++) { fgets(read, 200,fInputFile); }
371 fscanf(fInputFile, "%d", &ntracks);
372 fgets(read, 200, fInputFile);
373 fgets(read, 200, fInputFile);
374#endif
375
376 cout << "-I MpdUrqmdGenerator: Event " << evnr << " skipped!" << endl;
377
378 // ---> Loop over tracks in the current event
379 for (int itrack = 0; itrack < ntracks; itrack++) {
380
381 // Read momentum and PID from file
382#ifdef GZIP_SUPPORT
383 gzgets(fInputFile, read, 200);
384#else
385 fgets(read, 81, fInputFile);
386 fgets(read, 200, fInputFile);
387#endif
388 }
389 }
390 return kTRUE;
391}
392// ------------------------------------------------------------------------
393
394// ----- Private method ReadConverisonTable ---------------------------
395
396void MpdUrqmdGenerator::ReadConversionTable() {
397
398 TString work = getenv("VMCWORKDIR");
399 TString fileName = work + "/input/urqmd_pdg.dat";
400 ifstream* pdgconv = new ifstream(fileName.Data());
401
402 Int_t index = 0;
403 Int_t pdgId = 0;
404
405 while (!pdgconv->eof()) {
406 index = pdgId = 0;
407 *pdgconv >> index >> pdgId;
408 fParticleTable[index] = pdgId;
409 }
410
411 pdgconv->close();
412 delete pdgconv;
413
414 cout << "-I MpdUrqmdGenerator: Particle table for conversion from "
415 << "UrQMD loaded" << endl;
416
417}
418// ------------------------------------------------------------------------
419
420void MpdUrqmdGenerator::SetEventPlane(Double_t phiMin, Double_t phiMax) {
421 fPhiMin = phiMin;
422 fPhiMax = phiMax;
423 fEventPlaneSet = kTRUE;
424}
fPhiMin
Double_t fPhiMin
Definition BmnParticleEqualizer.cxx:47
fPhiMax
Double_t fPhiMax
Definition BmnParticleEqualizer.cxx:47
kProtonMass
const Double_t kProtonMass
Definition MpdLAQGSMGenerator.cxx:34
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition P4_F32vec4.h:34
m
__m128 m
Definition P4_F32vec4.h:27
MpdUrqmdGenerator::ReadEvent
Bool_t ReadEvent(FairPrimaryGenerator *primGen)
Definition MpdUrqmdGenerator.cxx:89
MpdUrqmdGenerator::SkipEvents
Bool_t SkipEvents(Int_t count)
Definition MpdUrqmdGenerator.cxx:286
MpdUrqmdGenerator::SetEventPlane
void SetEventPlane(Double_t phiMin, Double_t phiMax)
Definition MpdUrqmdGenerator.cxx:420
gCoordinateSystem
const CoordinateSystem gCoordinateSystem
Definition constants.h:2
sysLaboratory
@ sysLaboratory
Definition constants.h:1
std
STL namespace.