]> 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
Moved to STEERBase
[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),
2d22f1fe 129 fPpn(0.),
f97d4a8e 130 fCMEn(CMEn),
131 fIdp(0),
132 fBmin(0.),
133 fBmax(0.),
134 fFCentr(0),
135 fPi0Decay(0),
136 fProcess(iproc)
1a52e0ed 137{
d38751f5 138 printf("TDPMJet Constructor %d %d %d %d \n", Ip, Ipz, It, Itz);
1a52e0ed 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 }
d38751f5 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);
1a52e0ed 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 //
1a52e0ed 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));
1a52e0ed 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];
d38751f5 219/*
1a52e0ed 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]);
d38751f5 223*/
1a52e0ed 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]);
1a52e0ed 241 } // Particle loop
242 }
243 return nump;
244}
245
246
247//====================== access to dpmjet subroutines =========================
248//______________________________________________________________________________
249void TDPMjet::Initialize()
250{
d38751f5 251//
252// Write standard DPMJET input cards
253//
254 FILE* out = fopen("dpmjet.inp","w");
255// Projectile and Target definition
256 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.);
257 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.);
258// Beam energy and crossing-angle
259 fprintf(out, "BEAM %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",fEpn, fEpn, 0., 0., 0., 0.);
260// Centrality
261 fprintf(out, "CENTRAL %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",-1., fBmin, fBmax, 0., 0., 0.);
f97d4a8e 262// Particle decays
263 if (fPi0Decay)
264 fprintf(out, "PARDECAY %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n", 2., 0., 0., 0., 0., 0.);
d38751f5 265//
266// PHOJET specific
267 fprintf(out, "PHOINPUT\n");
268 fprintf(out, "DEBUG 0 0 0 \n");
6a45631b 269
d38751f5 270 if (fProcess == kDpmMb) {
271 fprintf(out, "PROCESS 1 0 1 1 1 1 1 1\n");
272 } else if (fProcess == kDpmMbNonDiffr) {
273 fprintf(out, "PROCESS 1 0 1 1 0 0 0 1\n");
6a45631b 274 } else if (fProcess == kDpmDiffr) {
275 fprintf(out, "PROCESS 0 0 0 0 1 1 1 0\n");
276 }else if (fProcess == kDpmSingleDiffr) {
277 fprintf(out, "PROCESS 0 0 0 0 1 1 0 0\n");
278 }else if (fProcess == kDpmDoubleDiffr) {
279 fprintf(out, "PROCESS 0 0 0 0 0 0 1 0\n");
d38751f5 280 }
281
282 fprintf(out, "ENDINPUT\n");
283//
284// START card
285 fprintf(out, "START 1.0 0.0\n");
286 fprintf(out, "STOP\n");
287 fclose(out);
db274076 288 dpmjet_openinp();
289
d38751f5 290//
291// Call DPMJET initialisation
292 Int_t iemu = 0; // No emulsion (default)
293 Dt_Dtuini(1, fEpn, fIp, fIpz, fIt, fItz, fIdp, iemu);
1a52e0ed 294
295}
296
297
298//______________________________________________________________________________
299void TDPMjet::GenerateEvent()
300{
301 // Generates one event;
302 fNEvent++;
303 DTEVNO.nevent=fNEvent;
1a52e0ed 304 Int_t kkmat=-1;
305 Float_t Elab = fEpn;
306 Int_t irej=0;
307 Dt_Kkinc(fIp, fIpz, fIt, fItz, fIdp, Elab, kkmat, irej);
308 if(irej!=0) return;
1a52e0ed 309}
310//______________________________________________________________________________
311void TDPMjet::Dt_Dtuini(int nevts, double epn, int npmass, int npchar,
312 int ntmass, int ntchar, int idp, int iemu)
313{
314 // Call dmpjet routine DT_DTUINI passing the parameters
315 // in a way accepted by Fortran routines
316
d38751f5 317
318 printf("\n-------------------------------------------\n");
1a52e0ed 319 printf("\n Dt_Dtuini called with:\n\n");
320 printf(" Projectile -> A = %d, Z = %d \n",npmass, npchar);
321 printf(" Target -> A = %d, Z = %d \n",ntmass, ntchar);
322 printf(" Proj. LAB E -> E = %f GeV \n",epn);
323 printf(" nevts = %d, idp = %d, iemu = %d \n",nevts,idp,iemu);
d38751f5 324 printf("\n-------------------------------------------\n");
1a52e0ed 325
326 dt_dtuini(nevts, epn, npmass, npchar, ntmass, ntchar, idp, iemu);
327
328}
329
330//______________________________________________________________________________
331void TDPMjet::Dt_Kkinc(int npmass, int npchar, int ntmass, int ntchar,
332 int idp, double elab, int kkmat, int irej)
333{
334 // Call dmpjet routine DT_KKINC passing the parameters
335 // in a way accepted by Fortran routines
1a52e0ed 336 dt_kkinc(npmass, npchar, ntmass, ntchar, idp, elab, kkmat, irej);
1a52e0ed 337}
338
339//______________________________________________________________________________
340void TDPMjet::Pho_Phist(int imode, double weight)
341{
342 // Call dmpjet routine PHO_PHIST passing the parameters
343 // in a way accepted by Fortran routines
344
345 pho_phist(imode,weight);
346
347}
348
349//______________________________________________________________________________
350void TDPMjet::Dt_Dtuout()
351{
352 // Call dmpjet routine DT_DTUOT passing the parameters
353 // in a way accepted by Fortran routines
354
355 dt_dtuout();
356
357}
358
359//______________________________________________________________________________
360Int_t TDPMjet::GetEvNum() const
361{
362 return DTEVT1.nevhkk;
363}
364//______________________________________________________________________________
365Int_t TDPMjet::GetEntriesNum() const
366{
367 return DTEVT1.nhkk;
368}
369//______________________________________________________________________________
370Int_t TDPMjet::GetNumStablePc() const
371{
372 Int_t NumStablePc = 0;
373 for(Int_t i=0; i<DTEVT1.nhkk; i++){
374 if(DTEVT1.isthkk[i] == 1) NumStablePc++;
375 }
376 return NumStablePc;
377}
378
379//______________________________________________________________________________
380Float_t TDPMjet::GetTotEnergy() const
381{
382 Float_t TotEnergy = 0.;
383 for(Int_t i=0; i<DTEVT1.nhkk; i++){
384 if(DTEVT1.isthkk[i] == 1)
385 TotEnergy += DTEVT1.phkk[i][3]; // PHKK[i][3] <-> PHKK(4,i)
386 }
387 return TotEnergy;
388}
389
390//______________________________________________________________________________
391Int_t TDPMjet::GetStatusCode(Int_t evnum) const
392{
393 return DTEVT1.isthkk[evnum];
394}
395//______________________________________________________________________________
396Int_t TDPMjet::GetPDGCode(Int_t evnum) const
397{
398 return DTEVT1.idhkk[evnum];
399}
400//______________________________________________________________________________
401Double_t TDPMjet::Getpx(Int_t evnum) const
402{
403 return DTEVT1.phkk[evnum][0];
404}
405//______________________________________________________________________________
406Double_t TDPMjet::Getpy(Int_t evnum) const
407{
408 return DTEVT1.phkk[evnum][1];
409}
410//______________________________________________________________________________
411Double_t TDPMjet::Getpz(Int_t evnum) const
412{
413 return DTEVT1.phkk[evnum][2];
414}
415//______________________________________________________________________________
416Double_t TDPMjet::GetEnergy(Int_t evnum) const
417{
418 return DTEVT1.phkk[evnum][3];
419}
420//______________________________________________________________________________
421Double_t TDPMjet::GetMass(Int_t evnum) const
422{
423 return DTEVT1.phkk[evnum][4];
424}
425//______________________________________________________________________________
426Int_t TDPMjet::GetFragmentA(Int_t evnum) const
427{
428 return DTEVT2.idres[evnum];
429}
430//______________________________________________________________________________
431Int_t TDPMjet::GetFragmentZ(Int_t evnum) const
432{
433 return DTEVT2.idxres[evnum];
434}
435//______________________________________________________________________________
436Double_t TDPMjet::GetXSFrac() const
437{
438 return DTIMPA.xsfrac;
439}
440//______________________________________________________________________________
441Double_t TDPMjet::GetBImpac() const
442{
443 return DTGLCP.bimpac;
444}
445//______________________________________________________________________________
446Double_t TDPMjet::GetProjRadius() const
447{
448 return DTGLCP.rproj;
449}
450//______________________________________________________________________________
451Double_t TDPMjet::GetTargRadius() const
452{
453 return DTGLCP.rtarg;
454}
455//______________________________________________________________________________
456Int_t TDPMjet::GetProjWounded() const
457{
458 return DTGLCP.nwasam;
459}
460//______________________________________________________________________________
461Int_t TDPMjet::GetTargWounded() const
462{
463 return DTGLCP.nwbsam;
464}
465//______________________________________________________________________________
466Int_t TDPMjet::GetProjSpectators() const
467{
468 return DTGLCP.nwtaac;
469}
470//______________________________________________________________________________
471Int_t TDPMjet::GetTargSpectators() const
472{
473 return DTGLCP.nwtbac;
474}
9dcc0beb 475//______________________________________________________________________________
476Int_t TDPMjet::GetProcessCode() const
477{
478 return POPRCS.iproce;
479}
1a52e0ed 480//______________________________________________________________________________
481void TDPMjet::Dt_Rndm(int idummy)
482{
483 dt_rndm(idummy);
484}
485
486//______________________________________________________________________________
487void TDPMjet::Dt_Rndmst(int na1, int na2, int na3, int nb1)
488{
489 dt_rndmst(na1, na2, na3, nb1);
490}
491
492//______________________________________________________________________________
493void TDPMjet::Dt_Rndmin(int u, int c, int cd, int cm, int i, int j)
494{
495 dt_rndmin(u, c, cd, cm, i, j);
496}
497
498//______________________________________________________________________________
499void TDPMjet::Dt_Rndmou(int u, int c, int cd, int cm, int i, int j)
500{
501 dt_rndmou(u, c, cd, cm, i, j);
502}
503