]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TDPMjet/TDPMjet.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Open input files before initialisation of dpmjet.
[u/mrichter/AliRoot.git] / TDPMjet / TDPMjet.cxx
CommitLineData
1a52e0ed 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"
07c71c2e 42//*KEEP,TParticle,T=C++
1a52e0ed 43#include "TParticle.h"
07c71c2e 44#include "TClonesArray.h"
1a52e0ed 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_
db274076 61# define dpmjet_openinp dpmjet_openinp_
1a52e0ed 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
db274076 73# define dpmjet_openinp DPMJET_OPENINP
1a52e0ed 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 &);
db274076 89extern "C" void type_of_call dpmjet_openinp();
1a52e0ed 90
91#else
92
93#endif
94
95ClassImp(TDPMjet)
96
f97d4a8e 97
1a52e0ed 98//______________________________________________________________________________
f97d4a8e 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)
1a52e0ed 115{
f97d4a8e 116// Default Constructor
1a52e0ed 117}
118
119//______________________________________________________________________________
d38751f5 120TDPMjet::TDPMjet(DpmProcess_t iproc, Int_t Ip=208, Int_t Ipz=82, Int_t It=208, Int_t Itz=82,
f97d4a8e 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 fCMEn(CMEn),
130 fIdp(0),
131 fBmin(0.),
132 fBmax(0.),
133 fFCentr(0),
134 fPi0Decay(0),
135 fProcess(iproc)
1a52e0ed 136{
d38751f5 137 printf("TDPMJet Constructor %d %d %d %d \n", Ip, Ipz, It, Itz);
1a52e0ed 138}
139
140
141//______________________________________________________________________________
142Int_t TDPMjet::ImportParticles(TClonesArray *particles, Option_t *option)
143{
144//
145// Default primary creation method. It reads the /HEPEVT/ common block which
146// has been filled by the GenerateEvent method. If the event generator does
147// not use the HEPEVT common block, This routine has to be overloaded by
148// the subclasses.
149// The function loops on the generated particles and store them in
150// the TClonesArray pointed by the argument particles.
151// The default action is to store only the stable particles
152// This can be demanded explicitly by setting the option = "Final"
153// If the option = "All", all the particles are stored.
154//
155 if(particles==0) return 0;
156 TClonesArray &Particles = *particles;
157 Particles.Clear();
158 Int_t numpart = 0; // Total number of produced particles
159 Int_t numStabpart = 0; // Total number of produced stable particles
160 Double_t entot = 0; // Total energy in final state (from stable particles)
161
162 numpart = DTEVT1.nhkk;
163 for(Int_t i=0; i<numpart; i++){
164 if(DTEVT1.isthkk[i]==1 || DTEVT1.isthkk[i]==-1 || DTEVT1.isthkk[i]==1001){
165 numStabpart++;
166 entot += DTEVT1.phkk[i][3]; // PHKK[i][3] <-> PHKK(4,i)
167 }
168 }
d38751f5 169 //printf("\n TDPMjet: DPMJET stack contains %d particles", numpart);
170 // printf("\n TDPMjet: Final not decayed particles: %d", numStabpart);
171 //printf("\n TDPMjet: Total energy: %f GeV \n", entot);
1a52e0ed 172 Int_t nump = 0;
173
174 if(!strcmp(option,"") || !strcmp(option,"Final")){
175 for (Int_t i=0; i < numpart; i++) {
176
177 if (DTEVT1.isthkk[i] == 1) {
178 //
179 // Use the common block values for the TParticle constructor
180 //
1a52e0ed 181 new(Particles[nump]) TParticle(
182 DTEVT1.idhkk[i],
183 DTEVT1.isthkk[i],
184 -1,
185 -1,
186 -1,
187 -1,
188 DTEVT1.phkk[i][0],
189 DTEVT1.phkk[i][1],
190 DTEVT1.phkk[i][2],
191 DTEVT1.phkk[i][3],
192
193 DTEVT1.vhkk[i][0],
194 DTEVT1.vhkk[i][1],
195 DTEVT1.vhkk[i][2],
196 DTEVT1.vhkk[i][3]);
197 nump++;
198 }
199 }
200 }
201 else if(!strcmp(option,"All")){
202 nump = numpart;
203 for (Int_t i=0; i <= numpart; i++){
204
205 // DTEVT1.JMOHKK[i][0] pointer to the entry of the 1st mother of entry i
206 Int_t iParent = DTEVT1.jmohkk[i][0] - 1;
207
208 if(iParent >= 0){
209 TParticle *mother = (TParticle*) (Particles.UncheckedAt(iParent));
1a52e0ed 210 mother->SetLastDaughter(i);
211 if(mother->GetFirstDaughter() == -1) mother->SetFirstDaughter(i);
212 }
213 // --- PDGcode for residual nuclei (idhkk=80000)
214 // --- 10000*Z + 10*A
215 // --- DPMJET -> idres = mass #, idxres = charge
216 if(DTEVT1.idhkk[i] == 80000)
217 DTEVT1.idhkk[i] = 10000*DTEVT2.idxres[i]+10*DTEVT2.idres[i];
d38751f5 218/*
1a52e0ed 219 if(DTEVT2.idxres[i] != 0)
220 printf("\n pc#%d -> A = %d, Z = %d -> PDGcode = %d\n",
221 i,DTEVT2.idres[i],DTEVT2.idxres[i],DTEVT1.idhkk[i]);
d38751f5 222*/
1a52e0ed 223 new(Particles[i]) TParticle(
224 DTEVT1.idhkk[i],
225 DTEVT1.isthkk[i],
226 iParent,
227 -1,
228 -1,
229 -1,
230
231 DTEVT1.phkk[i][0],
232 DTEVT1.phkk[i][1],
233 DTEVT1.phkk[i][2],
234 DTEVT1.phkk[i][3],
235
236 DTEVT1.vhkk[i][0],
237 DTEVT1.vhkk[i][1],
238 DTEVT1.vhkk[i][2],
239 DTEVT1.vhkk[i][3]);
1a52e0ed 240 } // Particle loop
241 }
242 return nump;
243}
244
245
246//====================== access to dpmjet subroutines =========================
247//______________________________________________________________________________
248void TDPMjet::Initialize()
249{
d38751f5 250//
251// Write standard DPMJET input cards
252//
253 FILE* out = fopen("dpmjet.inp","w");
254// Projectile and Target definition
255 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.);
256 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.);
257// Beam energy and crossing-angle
258 fprintf(out, "BEAM %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",fEpn, fEpn, 0., 0., 0., 0.);
259// Centrality
260 fprintf(out, "CENTRAL %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",-1., fBmin, fBmax, 0., 0., 0.);
f97d4a8e 261// Particle decays
262 if (fPi0Decay)
263 fprintf(out, "PARDECAY %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n", 2., 0., 0., 0., 0., 0.);
d38751f5 264//
265// PHOJET specific
266 fprintf(out, "PHOINPUT\n");
267 fprintf(out, "DEBUG 0 0 0 \n");
6a45631b 268
d38751f5 269 if (fProcess == kDpmMb) {
270 fprintf(out, "PROCESS 1 0 1 1 1 1 1 1\n");
271 } else if (fProcess == kDpmMbNonDiffr) {
272 fprintf(out, "PROCESS 1 0 1 1 0 0 0 1\n");
6a45631b 273 } else if (fProcess == kDpmDiffr) {
274 fprintf(out, "PROCESS 0 0 0 0 1 1 1 0\n");
275 }else if (fProcess == kDpmSingleDiffr) {
276 fprintf(out, "PROCESS 0 0 0 0 1 1 0 0\n");
277 }else if (fProcess == kDpmDoubleDiffr) {
278 fprintf(out, "PROCESS 0 0 0 0 0 0 1 0\n");
d38751f5 279 }
280
281 fprintf(out, "ENDINPUT\n");
282//
283// START card
284 fprintf(out, "START 1.0 0.0\n");
285 fprintf(out, "STOP\n");
286 fclose(out);
db274076 287 dpmjet_openinp();
288
d38751f5 289//
290// Call DPMJET initialisation
291 Int_t iemu = 0; // No emulsion (default)
292 Dt_Dtuini(1, fEpn, fIp, fIpz, fIt, fItz, fIdp, iemu);
1a52e0ed 293
294}
295
296
297//______________________________________________________________________________
298void TDPMjet::GenerateEvent()
299{
300 // Generates one event;
301 fNEvent++;
302 DTEVNO.nevent=fNEvent;
1a52e0ed 303 Int_t kkmat=-1;
304 Float_t Elab = fEpn;
305 Int_t irej=0;
306 Dt_Kkinc(fIp, fIpz, fIt, fItz, fIdp, Elab, kkmat, irej);
307 if(irej!=0) return;
1a52e0ed 308}
309//______________________________________________________________________________
310void TDPMjet::Dt_Dtuini(int nevts, double epn, int npmass, int npchar,
311 int ntmass, int ntchar, int idp, int iemu)
312{
313 // Call dmpjet routine DT_DTUINI passing the parameters
314 // in a way accepted by Fortran routines
315
d38751f5 316
317 printf("\n-------------------------------------------\n");
1a52e0ed 318 printf("\n Dt_Dtuini called with:\n\n");
319 printf(" Projectile -> A = %d, Z = %d \n",npmass, npchar);
320 printf(" Target -> A = %d, Z = %d \n",ntmass, ntchar);
321 printf(" Proj. LAB E -> E = %f GeV \n",epn);
322 printf(" nevts = %d, idp = %d, iemu = %d \n",nevts,idp,iemu);
d38751f5 323 printf("\n-------------------------------------------\n");
1a52e0ed 324
325 dt_dtuini(nevts, epn, npmass, npchar, ntmass, ntchar, idp, iemu);
326
327}
328
329//______________________________________________________________________________
330void TDPMjet::Dt_Kkinc(int npmass, int npchar, int ntmass, int ntchar,
331 int idp, double elab, int kkmat, int irej)
332{
333 // Call dmpjet routine DT_KKINC passing the parameters
334 // in a way accepted by Fortran routines
1a52e0ed 335 dt_kkinc(npmass, npchar, ntmass, ntchar, idp, elab, kkmat, irej);
1a52e0ed 336}
337
338//______________________________________________________________________________
339void TDPMjet::Pho_Phist(int imode, double weight)
340{
341 // Call dmpjet routine PHO_PHIST passing the parameters
342 // in a way accepted by Fortran routines
343
344 pho_phist(imode,weight);
345
346}
347
348//______________________________________________________________________________
349void TDPMjet::Dt_Dtuout()
350{
351 // Call dmpjet routine DT_DTUOT passing the parameters
352 // in a way accepted by Fortran routines
353
354 dt_dtuout();
355
356}
357
358//______________________________________________________________________________
359Int_t TDPMjet::GetEvNum() const
360{
361 return DTEVT1.nevhkk;
362}
363//______________________________________________________________________________
364Int_t TDPMjet::GetEntriesNum() const
365{
366 return DTEVT1.nhkk;
367}
368//______________________________________________________________________________
369Int_t TDPMjet::GetNumStablePc() const
370{
371 Int_t NumStablePc = 0;
372 for(Int_t i=0; i<DTEVT1.nhkk; i++){
373 if(DTEVT1.isthkk[i] == 1) NumStablePc++;
374 }
375 return NumStablePc;
376}
377
378//______________________________________________________________________________
379Float_t TDPMjet::GetTotEnergy() const
380{
381 Float_t TotEnergy = 0.;
382 for(Int_t i=0; i<DTEVT1.nhkk; i++){
383 if(DTEVT1.isthkk[i] == 1)
384 TotEnergy += DTEVT1.phkk[i][3]; // PHKK[i][3] <-> PHKK(4,i)
385 }
386 return TotEnergy;
387}
388
389//______________________________________________________________________________
390Int_t TDPMjet::GetStatusCode(Int_t evnum) const
391{
392 return DTEVT1.isthkk[evnum];
393}
394//______________________________________________________________________________
395Int_t TDPMjet::GetPDGCode(Int_t evnum) const
396{
397 return DTEVT1.idhkk[evnum];
398}
399//______________________________________________________________________________
400Double_t TDPMjet::Getpx(Int_t evnum) const
401{
402 return DTEVT1.phkk[evnum][0];
403}
404//______________________________________________________________________________
405Double_t TDPMjet::Getpy(Int_t evnum) const
406{
407 return DTEVT1.phkk[evnum][1];
408}
409//______________________________________________________________________________
410Double_t TDPMjet::Getpz(Int_t evnum) const
411{
412 return DTEVT1.phkk[evnum][2];
413}
414//______________________________________________________________________________
415Double_t TDPMjet::GetEnergy(Int_t evnum) const
416{
417 return DTEVT1.phkk[evnum][3];
418}
419//______________________________________________________________________________
420Double_t TDPMjet::GetMass(Int_t evnum) const
421{
422 return DTEVT1.phkk[evnum][4];
423}
424//______________________________________________________________________________
425Int_t TDPMjet::GetFragmentA(Int_t evnum) const
426{
427 return DTEVT2.idres[evnum];
428}
429//______________________________________________________________________________
430Int_t TDPMjet::GetFragmentZ(Int_t evnum) const
431{
432 return DTEVT2.idxres[evnum];
433}
434//______________________________________________________________________________
435Double_t TDPMjet::GetXSFrac() const
436{
437 return DTIMPA.xsfrac;
438}
439//______________________________________________________________________________
440Double_t TDPMjet::GetBImpac() const
441{
442 return DTGLCP.bimpac;
443}
444//______________________________________________________________________________
445Double_t TDPMjet::GetProjRadius() const
446{
447 return DTGLCP.rproj;
448}
449//______________________________________________________________________________
450Double_t TDPMjet::GetTargRadius() const
451{
452 return DTGLCP.rtarg;
453}
454//______________________________________________________________________________
455Int_t TDPMjet::GetProjWounded() const
456{
457 return DTGLCP.nwasam;
458}
459//______________________________________________________________________________
460Int_t TDPMjet::GetTargWounded() const
461{
462 return DTGLCP.nwbsam;
463}
464//______________________________________________________________________________
465Int_t TDPMjet::GetProjSpectators() const
466{
467 return DTGLCP.nwtaac;
468}
469//______________________________________________________________________________
470Int_t TDPMjet::GetTargSpectators() const
471{
472 return DTGLCP.nwtbac;
473}
9dcc0beb 474//______________________________________________________________________________
475Int_t TDPMjet::GetProcessCode() const
476{
477 return POPRCS.iproce;
478}
1a52e0ed 479//______________________________________________________________________________
480void TDPMjet::Dt_Rndm(int idummy)
481{
482 dt_rndm(idummy);
483}
484
485//______________________________________________________________________________
486void TDPMjet::Dt_Rndmst(int na1, int na2, int na3, int nb1)
487{
488 dt_rndmst(na1, na2, na3, nb1);
489}
490
491//______________________________________________________________________________
492void TDPMjet::Dt_Rndmin(int u, int c, int cd, int cm, int i, int j)
493{
494 dt_rndmin(u, c, cd, cm, i, j);
495}
496
497//______________________________________________________________________________
498void TDPMjet::Dt_Rndmou(int u, int c, int cd, int cm, int i, int j)
499{
500 dt_rndmou(u, c, cd, cm, i, j);
501}
502