]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - TDPMjet/TDPMjet.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
providing the RawReader already to the configuration step in order to allow modules...
[u/mrichter/AliRoot.git] / TDPMjet / TDPMjet.cxx
... / ...
CommitLineData
1//*KEEP,CopyRight,T=C.
2/*************************************************************************
3 * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
4 * All rights reserved. *
5 * *
6 * For the licensing terms see $ROOTSYS/LICENSE. *
7 * For the list of contributors see $ROOTSYS/README/CREDITS. *
8 *************************************************************************/
9//*KEND.
10
11//* +-------------------------------------------------------------+
12//* | |
13//* | |
14//* | DPMJET 3.0 |
15//* | |
16//* | |
17//* | S. Roesler+), R. Engel#), J. Ranft*) |
18//* | |
19//* | +) CERN, TIS-RP |
20//* | CH-1211 Geneva 23, Switzerland |
21//* | Email: Stefan.Roesler@cern.ch |
22//* | |
23//* | #) University of Delaware, BRI |
24//* | Newark, DE 19716, USA |
25//* | |
26//* | *) University of Siegen, Dept. of Physics |
27//* | D-57068 Siegen, Germany |
28//* | |
29//* | |
30//* | http://home.cern.ch/sroesler/dpmjet3.html |
31//* | |
32//* | |
33//* | Monte Carlo models used for event generation: |
34//* | PHOJET 1.12, JETSET 7.4 and LEPTO 6.5.1 |
35//* | |
36//* +-------------------------------------------------------------+
37
38//*KEEP,TDPMjet.
39#include "TDPMjet.h"
40//*KEEP,DPMCOMMON.
41#include "DPMcommon.h"
42//*KEEP,TParticle,T=C++
43#include "TParticle.h"
44#include "TClonesArray.h"
45//*KEND.
46
47//*KEEP,TROOT.
48#include "TROOT.h"
49//*KEND.
50
51#ifndef WIN32
52# define dt_dtuini dt_dtuini_
53# define dt_getemu de_getemu_
54# define dt_kkinc dt_kkinc_
55# define pho_phist pho_phist_
56# define dt_dtuout dt_dtuout_
57# define dt_rndm dt_rndm_
58# define dt_rndmst dt_rndmst_
59# define dt_rndmin dt_rndmin_
60# define dt_rndmou dt_rndmou_
61# define dpmjet_openinp dpmjet_openinp_
62# define type_of_call
63#else
64# define dt_dtuini DT_DTUINI
65# define dt_getemu DT_GETEMU
66# define dt_kkinc DT_KKINC
67# define pho_phist PHO_PHIST
68# define dt_dtuout DT_DTUOUT
69# define dt_rndm DT_RNDM
70# define dt_rndmst DT_RNDMST
71# define dt_rndmin DT_RNDMIN
72# define dt_rndmou DT_RNDMOU
73# define dpmjet_openinp DPMJET_OPENINP
74# define type_of_call _stdcall
75#endif
76
77#ifndef WIN32
78extern "C" void type_of_call dt_dtuini(Int_t & , Double_t &, Int_t & , Int_t &,
79 Int_t &, Int_t &, Int_t &, Int_t &);
80extern "C" double type_of_call dt_getemu(Int_t &, Int_t &, Int_t &, Int_t &);
81extern "C" void type_of_call dt_kkinc(Int_t &, Int_t &, Int_t &, Int_t &,
82 Int_t &, Double_t &, Int_t &, Int_t &);
83extern "C" void type_of_call pho_phist(Int_t &, Double_t &);
84extern "C" void type_of_call dt_dtuout();
85extern "C" void type_of_call dt_rndm(Int_t &);
86extern "C" void type_of_call dt_rndmst(Int_t &, Int_t &, Int_t &, Int_t &);
87extern "C" void type_of_call dt_rndmin(Int_t &, Int_t &, Int_t &, Int_t &, Int_t &, Int_t &);
88extern "C" void type_of_call dt_rndmou(Int_t &, Int_t &, Int_t &, Int_t &, Int_t &, Int_t &);
89extern "C" void type_of_call dpmjet_openinp();
90
91#else
92
93#endif
94
95ClassImp(TDPMjet)
96
97
98//______________________________________________________________________________
99 TDPMjet::TDPMjet() :
100 TGenerator("dpmjet","dpmjet"),
101 fNEvent(0),
102 fIp(0),
103 fIpz(0),
104 fIt(0),
105 fItz(0),
106 fEpn(0.),
107 fPpn(0.),
108 fCMEn(0.),
109 fIdp(0),
110 fBmin(0.),
111 fBmax(0.),
112 fFCentr(0),
113 fPi0Decay(0),
114 fProcess(kDpmMb)
115{
116// Default Constructor
117}
118
119//______________________________________________________________________________
120TDPMjet::TDPMjet(DpmProcess_t iproc, Int_t Ip=208, Int_t Ipz=82, Int_t It=208, Int_t Itz=82,
121 Double_t Epn=2700., Double_t CMEn=5400.)
122 : TGenerator("dpmjet","dpmjet"),
123 fNEvent(0),
124 fIp(Ip),
125 fIpz(Ipz),
126 fIt(It),
127 fItz(Itz),
128 fEpn(Epn),
129 fPpn(0.),
130 fCMEn(CMEn),
131 fIdp(0),
132 fBmin(0.),
133 fBmax(0.),
134 fFCentr(0),
135 fPi0Decay(0),
136 fProcess(iproc)
137{
138 printf("TDPMJet Constructor %d %d %d %d \n", Ip, Ipz, It, Itz);
139}
140
141
142//______________________________________________________________________________
143Int_t TDPMjet::ImportParticles(TClonesArray *particles, Option_t *option)
144{
145//
146// Default primary creation method. It reads the /HEPEVT/ common block which
147// has been filled by the GenerateEvent method. If the event generator does
148// not use the HEPEVT common block, This routine has to be overloaded by
149// the subclasses.
150// The function loops on the generated particles and store them in
151// the TClonesArray pointed by the argument particles.
152// The default action is to store only the stable particles
153// This can be demanded explicitly by setting the option = "Final"
154// If the option = "All", all the particles are stored.
155//
156 if(particles==0) return 0;
157 TClonesArray &Particles = *particles;
158 Particles.Clear();
159 Int_t numpart = 0; // Total number of produced particles
160 Int_t numStabpart = 0; // Total number of produced stable particles
161 Double_t entot = 0; // Total energy in final state (from stable particles)
162
163 numpart = DTEVT1.nhkk;
164 for(Int_t i=0; i<numpart; i++){
165 if(DTEVT1.isthkk[i]==1 || DTEVT1.isthkk[i]==-1 || DTEVT1.isthkk[i]==1001){
166 numStabpart++;
167 entot += DTEVT1.phkk[i][3]; // PHKK[i][3] <-> PHKK(4,i)
168 }
169 }
170 //printf("\n TDPMjet: DPMJET stack contains %d particles", numpart);
171 // printf("\n TDPMjet: Final not decayed particles: %d", numStabpart);
172 //printf("\n TDPMjet: Total energy: %f GeV \n", entot);
173 Int_t nump = 0;
174
175 if(!strcmp(option,"") || !strcmp(option,"Final")){
176 for (Int_t i=0; i < numpart; i++) {
177
178 if (DTEVT1.isthkk[i] == 1) {
179 //
180 // Use the common block values for the TParticle constructor
181 //
182 new(Particles[nump]) TParticle(
183 DTEVT1.idhkk[i],
184 DTEVT1.isthkk[i],
185 -1,
186 -1,
187 -1,
188 -1,
189 DTEVT1.phkk[i][0],
190 DTEVT1.phkk[i][1],
191 DTEVT1.phkk[i][2],
192 DTEVT1.phkk[i][3],
193
194 DTEVT1.vhkk[i][0],
195 DTEVT1.vhkk[i][1],
196 DTEVT1.vhkk[i][2],
197 DTEVT1.vhkk[i][3]);
198 nump++;
199 }
200 }
201 }
202 else if(!strcmp(option,"All")){
203 nump = numpart;
204 for (Int_t i=0; i <= numpart; i++){
205
206 // DTEVT1.JMOHKK[i][0] pointer to the entry of the 1st mother of entry i
207 Int_t iParent = DTEVT1.jmohkk[i][0] - 1;
208
209 if(iParent >= 0){
210 TParticle *mother = (TParticle*) (Particles.UncheckedAt(iParent));
211 mother->SetLastDaughter(i);
212 if(mother->GetFirstDaughter() == -1) mother->SetFirstDaughter(i);
213 }
214 // --- PDGcode for residual nuclei (idhkk=80000)
215 // --- 10000*Z + 10*A
216 // --- DPMJET -> idres = mass #, idxres = charge
217 if(DTEVT1.idhkk[i] == 80000)
218 DTEVT1.idhkk[i] = 10000*DTEVT2.idxres[i]+10*DTEVT2.idres[i];
219/*
220 if(DTEVT2.idxres[i] != 0)
221 printf("\n pc#%d -> A = %d, Z = %d -> PDGcode = %d\n",
222 i,DTEVT2.idres[i],DTEVT2.idxres[i],DTEVT1.idhkk[i]);
223*/
224 new(Particles[i]) TParticle(
225 DTEVT1.idhkk[i],
226 DTEVT1.isthkk[i],
227 iParent,
228 -1,
229 -1,
230 -1,
231
232 DTEVT1.phkk[i][0],
233 DTEVT1.phkk[i][1],
234 DTEVT1.phkk[i][2],
235 DTEVT1.phkk[i][3],
236
237 DTEVT1.vhkk[i][0],
238 DTEVT1.vhkk[i][1],
239 DTEVT1.vhkk[i][2],
240 DTEVT1.vhkk[i][3]);
241 } // Particle loop
242 }
243 return nump;
244}
245
246
247//====================== access to dpmjet subroutines =========================
248//______________________________________________________________________________
249void TDPMjet::Initialize()
250{
251//
252// Write standard DPMJET input cards
253//
254 FILE* out = fopen("dpmjet.inp","w");
255// Projectile and Target definition
256 if (fIp == 1 && fIpz ==1) {
257 fprintf(out, "PROJPAR PROTON\n");
258 } else if (fIp == 1 && fIpz == -1) {
259 fprintf(out, "PROJPAR APROTON\n");
260 } else {
261 fprintf(out, "PROJPAR %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n", (Float_t) fIp, (Float_t) fIpz, 0., 0., 0., 0.);
262 }
263
264 if (fIt == 1 && fItz ==1) {
265 fprintf(out, "TARPAR PROTON\n");
266 } else if (fIt == 1 && fItz == -1) {
267 fprintf(out, "TARPAR APROTON\n");
268 } else {
269 fprintf(out, "TARPAR %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n", (Float_t) fIt, (Float_t) fItz, 0., 0., 0., 0.);
270 }
271
272// Beam energy and crossing-angle
273 fprintf(out, "BEAM %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",fEpn, fEpn, 0., 0., 0., 0.);
274// Centrality
275 fprintf(out, "CENTRAL %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",-1., fBmin, fBmax, 0., 0., 0.);
276// Particle decays
277 if (fPi0Decay)
278 fprintf(out, "PARDECAY %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n", 2., 0., 0., 0., 0., 0.);
279//
280// PHOJET specific
281 fprintf(out, "PHOINPUT\n");
282 fprintf(out, "DEBUG 0 0 0 \n");
283
284 if (fProcess == kDpmMb) {
285 fprintf(out, "PROCESS 1 0 1 1 1 1 1 1\n");
286 } else if (fProcess == kDpmMbNonDiffr) {
287 fprintf(out, "PROCESS 1 0 1 1 0 0 0 1\n");
288 } else if (fProcess == kDpmDiffr) {
289 fprintf(out, "PROCESS 0 0 0 0 1 1 1 0\n");
290 }else if (fProcess == kDpmSingleDiffr) {
291 fprintf(out, "PROCESS 0 0 0 0 1 1 0 0\n");
292 }else if (fProcess == kDpmDoubleDiffr) {
293 fprintf(out, "PROCESS 0 0 0 0 0 0 1 0\n");
294 }
295
296 fprintf(out, "ENDINPUT\n");
297//
298// START card
299 fprintf(out, "START 1.0 0.0\n");
300 fprintf(out, "STOP\n");
301 fclose(out);
302 dpmjet_openinp();
303
304//
305// Call DPMJET initialisation
306 Int_t iemu = 0; // No emulsion (default)
307 Dt_Dtuini(1, fEpn, fIp, fIpz, fIt, fItz, fIdp, iemu);
308
309}
310
311
312//______________________________________________________________________________
313void TDPMjet::GenerateEvent()
314{
315 // Generates one event;
316 fNEvent++;
317 DTEVNO.nevent=fNEvent;
318 Int_t kkmat=-1;
319 Float_t Elab = fEpn;
320 Int_t irej=0;
321 Dt_Kkinc(fIp, fIpz, fIt, fItz, fIdp, Elab, kkmat, irej);
322 if(irej!=0) return;
323}
324//______________________________________________________________________________
325void TDPMjet::Dt_Dtuini(int nevts, double epn, int npmass, int npchar,
326 int ntmass, int ntchar, int idp, int iemu)
327{
328 // Call dmpjet routine DT_DTUINI passing the parameters
329 // in a way accepted by Fortran routines
330
331
332 printf("\n-------------------------------------------\n");
333 printf("\n Dt_Dtuini called with:\n\n");
334 printf(" Projectile -> A = %d, Z = %d \n",npmass, npchar);
335 printf(" Target -> A = %d, Z = %d \n",ntmass, ntchar);
336 printf(" Proj. LAB E -> E = %f GeV \n",epn);
337 printf(" nevts = %d, idp = %d, iemu = %d \n",nevts,idp,iemu);
338 printf("\n-------------------------------------------\n");
339
340 dt_dtuini(nevts, epn, npmass, npchar, ntmass, ntchar, idp, iemu);
341
342}
343
344//______________________________________________________________________________
345void TDPMjet::Dt_Kkinc(int npmass, int npchar, int ntmass, int ntchar,
346 int idp, double elab, int kkmat, int irej)
347{
348 // Call dmpjet routine DT_KKINC passing the parameters
349 // in a way accepted by Fortran routines
350 dt_kkinc(npmass, npchar, ntmass, ntchar, idp, elab, kkmat, irej);
351}
352
353//______________________________________________________________________________
354void TDPMjet::Pho_Phist(int imode, double weight)
355{
356 // Call dmpjet routine PHO_PHIST passing the parameters
357 // in a way accepted by Fortran routines
358
359 pho_phist(imode,weight);
360
361}
362
363//______________________________________________________________________________
364void TDPMjet::Dt_Dtuout()
365{
366 // Call dmpjet routine DT_DTUOT passing the parameters
367 // in a way accepted by Fortran routines
368
369 dt_dtuout();
370
371}
372
373//______________________________________________________________________________
374Int_t TDPMjet::GetEvNum() const
375{
376 return DTEVT1.nevhkk;
377}
378//______________________________________________________________________________
379Int_t TDPMjet::GetEntriesNum() const
380{
381 return DTEVT1.nhkk;
382}
383//______________________________________________________________________________
384Int_t TDPMjet::GetNumStablePc() const
385{
386 Int_t NumStablePc = 0;
387 for(Int_t i=0; i<DTEVT1.nhkk; i++){
388 if(DTEVT1.isthkk[i] == 1) NumStablePc++;
389 }
390 return NumStablePc;
391}
392
393//______________________________________________________________________________
394Float_t TDPMjet::GetTotEnergy() const
395{
396 Float_t TotEnergy = 0.;
397 for(Int_t i=0; i<DTEVT1.nhkk; i++){
398 if(DTEVT1.isthkk[i] == 1)
399 TotEnergy += DTEVT1.phkk[i][3]; // PHKK[i][3] <-> PHKK(4,i)
400 }
401 return TotEnergy;
402}
403
404//______________________________________________________________________________
405Int_t TDPMjet::GetStatusCode(Int_t evnum) const
406{
407 return DTEVT1.isthkk[evnum];
408}
409//______________________________________________________________________________
410Int_t TDPMjet::GetPDGCode(Int_t evnum) const
411{
412 return DTEVT1.idhkk[evnum];
413}
414//______________________________________________________________________________
415Double_t TDPMjet::Getpx(Int_t evnum) const
416{
417 return DTEVT1.phkk[evnum][0];
418}
419//______________________________________________________________________________
420Double_t TDPMjet::Getpy(Int_t evnum) const
421{
422 return DTEVT1.phkk[evnum][1];
423}
424//______________________________________________________________________________
425Double_t TDPMjet::Getpz(Int_t evnum) const
426{
427 return DTEVT1.phkk[evnum][2];
428}
429//______________________________________________________________________________
430Double_t TDPMjet::GetEnergy(Int_t evnum) const
431{
432 return DTEVT1.phkk[evnum][3];
433}
434//______________________________________________________________________________
435Double_t TDPMjet::GetMass(Int_t evnum) const
436{
437 return DTEVT1.phkk[evnum][4];
438}
439//______________________________________________________________________________
440Int_t TDPMjet::GetFragmentA(Int_t evnum) const
441{
442 return DTEVT2.idres[evnum];
443}
444//______________________________________________________________________________
445Int_t TDPMjet::GetFragmentZ(Int_t evnum) const
446{
447 return DTEVT2.idxres[evnum];
448}
449//______________________________________________________________________________
450Double_t TDPMjet::GetXSFrac() const
451{
452 return DTIMPA.xsfrac;
453}
454//______________________________________________________________________________
455Double_t TDPMjet::GetBImpac() const
456{
457 return DTGLCP.bimpac;
458}
459//______________________________________________________________________________
460Double_t TDPMjet::GetProjRadius() const
461{
462 return DTGLCP.rproj;
463}
464//______________________________________________________________________________
465Double_t TDPMjet::GetTargRadius() const
466{
467 return DTGLCP.rtarg;
468}
469//______________________________________________________________________________
470Int_t TDPMjet::GetProjWounded() const
471{
472 return DTGLCP.nwasam;
473}
474//______________________________________________________________________________
475Int_t TDPMjet::GetTargWounded() const
476{
477 return DTGLCP.nwbsam;
478}
479//______________________________________________________________________________
480Int_t TDPMjet::GetProjSpectators() const
481{
482 return DTGLCP.nwtaac;
483}
484//______________________________________________________________________________
485Int_t TDPMjet::GetTargSpectators() const
486{
487 return DTGLCP.nwtbac;
488}
489//______________________________________________________________________________
490Int_t TDPMjet::GetProcessCode() const
491{
492 return POPRCS.iproce;
493}
494//______________________________________________________________________________
495void TDPMjet::Dt_Rndm(int idummy)
496{
497 dt_rndm(idummy);
498}
499
500//______________________________________________________________________________
501void TDPMjet::Dt_Rndmst(int na1, int na2, int na3, int nb1)
502{
503 dt_rndmst(na1, na2, na3, nb1);
504}
505
506//______________________________________________________________________________
507void TDPMjet::Dt_Rndmin(int u, int c, int cd, int cm, int i, int j)
508{
509 dt_rndmin(u, c, cd, cm, i, j);
510}
511
512//______________________________________________________________________________
513void TDPMjet::Dt_Rndmou(int u, int c, int cd, int cm, int i, int j)
514{
515 dt_rndmou(u, c, cd, cm, i, j);
516}
517