]> 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
Rename AliTRDrawStream to AliTRDrawStreamOld
[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_
06f923c8 62# define dt_evtout dt_evtout_
1a52e0ed 63# define type_of_call
64#else
65# define dt_dtuini DT_DTUINI
66# define dt_getemu DT_GETEMU
67# define dt_kkinc DT_KKINC
68# define pho_phist PHO_PHIST
69# define dt_dtuout DT_DTUOUT
70# define dt_rndm DT_RNDM
71# define dt_rndmst DT_RNDMST
72# define dt_rndmin DT_RNDMIN
73# define dt_rndmou DT_RNDMOU
06f923c8 74# define dt_evtout DT_EVTOUT_
db274076 75# define dpmjet_openinp DPMJET_OPENINP
1a52e0ed 76# define type_of_call _stdcall
77#endif
78
06f923c8 79
1a52e0ed 80#ifndef WIN32
81extern "C" void type_of_call dt_dtuini(Int_t & , Double_t &, Int_t & , Int_t &,
82 Int_t &, Int_t &, Int_t &, Int_t &);
83extern "C" double type_of_call dt_getemu(Int_t &, Int_t &, Int_t &, Int_t &);
84extern "C" void type_of_call dt_kkinc(Int_t &, Int_t &, Int_t &, Int_t &,
85 Int_t &, Double_t &, Int_t &, Int_t &);
86extern "C" void type_of_call pho_phist(Int_t &, Double_t &);
87extern "C" void type_of_call dt_dtuout();
06f923c8 88extern "C" int type_of_call dt_evtout(const Int_t &);
1a52e0ed 89extern "C" void type_of_call dt_rndm(Int_t &);
90extern "C" void type_of_call dt_rndmst(Int_t &, Int_t &, Int_t &, Int_t &);
91extern "C" void type_of_call dt_rndmin(Int_t &, Int_t &, Int_t &, Int_t &, Int_t &, Int_t &);
92extern "C" void type_of_call dt_rndmou(Int_t &, Int_t &, Int_t &, Int_t &, Int_t &, Int_t &);
db274076 93extern "C" void type_of_call dpmjet_openinp();
1a52e0ed 94
95#else
96
97#endif
98
99ClassImp(TDPMjet)
100
f97d4a8e 101
1a52e0ed 102//______________________________________________________________________________
f97d4a8e 103 TDPMjet::TDPMjet() :
104 TGenerator("dpmjet","dpmjet"),
105 fNEvent(0),
106 fIp(0),
107 fIpz(0),
108 fIt(0),
109 fItz(0),
110 fEpn(0.),
111 fPpn(0.),
112 fCMEn(0.),
113 fIdp(0),
114 fBmin(0.),
115 fBmax(0.),
116 fFCentr(0),
117 fPi0Decay(0),
717765ce 118 fDecayAll(0),
f97d4a8e 119 fProcess(kDpmMb)
1a52e0ed 120{
f97d4a8e 121// Default Constructor
1a52e0ed 122}
123
124//______________________________________________________________________________
d38751f5 125TDPMjet::TDPMjet(DpmProcess_t iproc, Int_t Ip=208, Int_t Ipz=82, Int_t It=208, Int_t Itz=82,
f97d4a8e 126 Double_t Epn=2700., Double_t CMEn=5400.)
127 : TGenerator("dpmjet","dpmjet"),
128 fNEvent(0),
129 fIp(Ip),
130 fIpz(Ipz),
131 fIt(It),
132 fItz(Itz),
133 fEpn(Epn),
2d22f1fe 134 fPpn(0.),
f97d4a8e 135 fCMEn(CMEn),
136 fIdp(0),
137 fBmin(0.),
138 fBmax(0.),
139 fFCentr(0),
140 fPi0Decay(0),
717765ce 141 fDecayAll(0),
f97d4a8e 142 fProcess(iproc)
1a52e0ed 143{
d38751f5 144 printf("TDPMJet Constructor %d %d %d %d \n", Ip, Ipz, It, Itz);
1a52e0ed 145}
146
147
148//______________________________________________________________________________
149Int_t TDPMjet::ImportParticles(TClonesArray *particles, Option_t *option)
150{
151//
152// Default primary creation method. It reads the /HEPEVT/ common block which
153// has been filled by the GenerateEvent method. If the event generator does
154// not use the HEPEVT common block, This routine has to be overloaded by
155// the subclasses.
156// The function loops on the generated particles and store them in
157// the TClonesArray pointed by the argument particles.
158// The default action is to store only the stable particles
159// This can be demanded explicitly by setting the option = "Final"
160// If the option = "All", all the particles are stored.
161//
162 if(particles==0) return 0;
163 TClonesArray &Particles = *particles;
164 Particles.Clear();
165 Int_t numpart = 0; // Total number of produced particles
166 Int_t numStabpart = 0; // Total number of produced stable particles
167 Double_t entot = 0; // Total energy in final state (from stable particles)
168
169 numpart = DTEVT1.nhkk;
170 for(Int_t i=0; i<numpart; i++){
171 if(DTEVT1.isthkk[i]==1 || DTEVT1.isthkk[i]==-1 || DTEVT1.isthkk[i]==1001){
172 numStabpart++;
173 entot += DTEVT1.phkk[i][3]; // PHKK[i][3] <-> PHKK(4,i)
174 }
175 }
cab41d67 176 printf("\n TDPMjet: DPMJET stack contains %d particles", numpart);
d38751f5 177 // printf("\n TDPMjet: Final not decayed particles: %d", numStabpart);
178 //printf("\n TDPMjet: Total energy: %f GeV \n", entot);
1a52e0ed 179 Int_t nump = 0;
180
181 if(!strcmp(option,"") || !strcmp(option,"Final")){
182 for (Int_t i=0; i < numpart; i++) {
183
184 if (DTEVT1.isthkk[i] == 1) {
185 //
186 // Use the common block values for the TParticle constructor
187 //
1a52e0ed 188 new(Particles[nump]) TParticle(
189 DTEVT1.idhkk[i],
190 DTEVT1.isthkk[i],
191 -1,
192 -1,
193 -1,
194 -1,
195 DTEVT1.phkk[i][0],
196 DTEVT1.phkk[i][1],
197 DTEVT1.phkk[i][2],
198 DTEVT1.phkk[i][3],
199
200 DTEVT1.vhkk[i][0],
201 DTEVT1.vhkk[i][1],
202 DTEVT1.vhkk[i][2],
203 DTEVT1.vhkk[i][3]);
204 nump++;
205 }
206 }
207 }
208 else if(!strcmp(option,"All")){
209 nump = numpart;
210 for (Int_t i=0; i <= numpart; i++){
211
212 // DTEVT1.JMOHKK[i][0] pointer to the entry of the 1st mother of entry i
213 Int_t iParent = DTEVT1.jmohkk[i][0] - 1;
214
215 if(iParent >= 0){
216 TParticle *mother = (TParticle*) (Particles.UncheckedAt(iParent));
1a52e0ed 217 mother->SetLastDaughter(i);
218 if(mother->GetFirstDaughter() == -1) mother->SetFirstDaughter(i);
219 }
220 // --- PDGcode for residual nuclei (idhkk=80000)
221 // --- 10000*Z + 10*A
222 // --- DPMJET -> idres = mass #, idxres = charge
223 if(DTEVT1.idhkk[i] == 80000)
224 DTEVT1.idhkk[i] = 10000*DTEVT2.idxres[i]+10*DTEVT2.idres[i];
d38751f5 225/*
1a52e0ed 226 if(DTEVT2.idxres[i] != 0)
227 printf("\n pc#%d -> A = %d, Z = %d -> PDGcode = %d\n",
228 i,DTEVT2.idres[i],DTEVT2.idxres[i],DTEVT1.idhkk[i]);
d38751f5 229*/
06f923c8 230/*
231 printf("%5d %5d %5d %5d %13.3f %13.3f %13.3f %13.3f \n", i,
232 DTEVT1.idhkk[i],
233 DTEVT1.isthkk[i],
234 iParent,
235 DTEVT1.phkk[i][0],
236 DTEVT1.phkk[i][1],
237 DTEVT1.phkk[i][2],
238 DTEVT1.phkk[i][3]
239 );
240*/
1a52e0ed 241 new(Particles[i]) TParticle(
242 DTEVT1.idhkk[i],
243 DTEVT1.isthkk[i],
244 iParent,
245 -1,
246 -1,
247 -1,
248
249 DTEVT1.phkk[i][0],
250 DTEVT1.phkk[i][1],
251 DTEVT1.phkk[i][2],
252 DTEVT1.phkk[i][3],
253
254 DTEVT1.vhkk[i][0],
255 DTEVT1.vhkk[i][1],
256 DTEVT1.vhkk[i][2],
257 DTEVT1.vhkk[i][3]);
1a52e0ed 258 } // Particle loop
259 }
260 return nump;
261}
262
263
264//====================== access to dpmjet subroutines =========================
265//______________________________________________________________________________
266void TDPMjet::Initialize()
267{
d38751f5 268//
269// Write standard DPMJET input cards
270//
271 FILE* out = fopen("dpmjet.inp","w");
272// Projectile and Target definition
cd19c617 273 if (fIp == 1 && fIpz ==1) {
274 fprintf(out, "PROJPAR PROTON\n");
275 } else if (fIp == 1 && fIpz == -1) {
276 fprintf(out, "PROJPAR APROTON\n");
277 } else {
278 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.);
279 }
280
281 if (fIt == 1 && fItz ==1) {
282 fprintf(out, "TARPAR PROTON\n");
283 } else if (fIt == 1 && fItz == -1) {
284 fprintf(out, "TARPAR APROTON\n");
285 } else {
286 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.);
287 }
288
d38751f5 289// Beam energy and crossing-angle
290 fprintf(out, "BEAM %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",fEpn, fEpn, 0., 0., 0., 0.);
291// Centrality
717765ce 292 if (fIp > 1. && fIt > 1)
293 fprintf(out, "CENTRAL %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",-1., fBmin, fBmax, 0., 0., 0.);
f97d4a8e 294// Particle decays
295 if (fPi0Decay)
717765ce 296 fprintf(out, "PARDECAY %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n", 2., 0., 0., 0., 0., 0.);
297
298
d38751f5 299//
300// PHOJET specific
301 fprintf(out, "PHOINPUT\n");
302 fprintf(out, "DEBUG 0 0 0 \n");
6a45631b 303
d38751f5 304 if (fProcess == kDpmMb) {
305 fprintf(out, "PROCESS 1 0 1 1 1 1 1 1\n");
306 } else if (fProcess == kDpmMbNonDiffr) {
cab41d67 307 fprintf(out, "PROCESS 1 0 1 0 0 0 0 1\n");
6a45631b 308 } else if (fProcess == kDpmDiffr) {
cab41d67 309 fprintf(out, "PROCESS 0 0 0 1 1 1 1 0\n");
6a45631b 310 }else if (fProcess == kDpmSingleDiffr) {
311 fprintf(out, "PROCESS 0 0 0 0 1 1 0 0\n");
312 }else if (fProcess == kDpmDoubleDiffr) {
313 fprintf(out, "PROCESS 0 0 0 0 0 0 1 0\n");
cab41d67 314 } else if (fProcess == kDpmCentralDiffr){
315 fprintf(out, "PROCESS 0 0 0 1 0 0 0 0\n");
d38751f5 316 }
317
06f923c8 318 Int_t iPDG[18] =
717765ce 319 {
320// K0s pi0 lam sig+ sig- tet0
321 310, 111, 3122, 3222, 3112, 3322,
322// tet- om- D+ D0 Ds+
323 3312, 3334, 411, 421, 431,
324// etac lamc+ sigc++ sigc+ sigc0 Ksic+
325 441, 4122, 4222, 4212, 4112, 4232,
326// Ksic0 sig0
06f923c8 327 4132
717765ce 328 };
329
330
331 Int_t iON = (fDecayAll) ? 1:0;
06f923c8 332 for (Int_t i = 0; i < 8; i++) {
717765ce 333 fprintf(out, "LUND-DECAY%5d %5d\n", iPDG[i], iON);
334 }
335
336 fprintf(out, "LUND-MSTJ %5d %5d\n", 22, 1);
337
d38751f5 338 fprintf(out, "ENDINPUT\n");
339//
340// START card
341 fprintf(out, "START 1.0 0.0\n");
342 fprintf(out, "STOP\n");
343 fclose(out);
db274076 344 dpmjet_openinp();
345
d38751f5 346//
347// Call DPMJET initialisation
348 Int_t iemu = 0; // No emulsion (default)
349 Dt_Dtuini(1, fEpn, fIp, fIpz, fIt, fItz, fIdp, iemu);
1a52e0ed 350
351}
352
353
354//______________________________________________________________________________
355void TDPMjet::GenerateEvent()
356{
357 // Generates one event;
358 fNEvent++;
359 DTEVNO.nevent=fNEvent;
1a52e0ed 360 Int_t kkmat=-1;
361 Float_t Elab = fEpn;
362 Int_t irej=0;
363 Dt_Kkinc(fIp, fIpz, fIt, fItz, fIdp, Elab, kkmat, irej);
06f923c8 364// dt_evtout(4);
1a52e0ed 365 if(irej!=0) return;
1a52e0ed 366}
367//______________________________________________________________________________
368void TDPMjet::Dt_Dtuini(int nevts, double epn, int npmass, int npchar,
369 int ntmass, int ntchar, int idp, int iemu)
370{
371 // Call dmpjet routine DT_DTUINI passing the parameters
372 // in a way accepted by Fortran routines
373
d38751f5 374
375 printf("\n-------------------------------------------\n");
1a52e0ed 376 printf("\n Dt_Dtuini called with:\n\n");
377 printf(" Projectile -> A = %d, Z = %d \n",npmass, npchar);
378 printf(" Target -> A = %d, Z = %d \n",ntmass, ntchar);
379 printf(" Proj. LAB E -> E = %f GeV \n",epn);
380 printf(" nevts = %d, idp = %d, iemu = %d \n",nevts,idp,iemu);
d38751f5 381 printf("\n-------------------------------------------\n");
1a52e0ed 382
383 dt_dtuini(nevts, epn, npmass, npchar, ntmass, ntchar, idp, iemu);
384
385}
386
387//______________________________________________________________________________
388void TDPMjet::Dt_Kkinc(int npmass, int npchar, int ntmass, int ntchar,
389 int idp, double elab, int kkmat, int irej)
390{
391 // Call dmpjet routine DT_KKINC passing the parameters
392 // in a way accepted by Fortran routines
1a52e0ed 393 dt_kkinc(npmass, npchar, ntmass, ntchar, idp, elab, kkmat, irej);
1a52e0ed 394}
395
396//______________________________________________________________________________
397void TDPMjet::Pho_Phist(int imode, double weight)
398{
399 // Call dmpjet routine PHO_PHIST passing the parameters
400 // in a way accepted by Fortran routines
401
402 pho_phist(imode,weight);
403
404}
405
406//______________________________________________________________________________
407void TDPMjet::Dt_Dtuout()
408{
409 // Call dmpjet routine DT_DTUOT passing the parameters
410 // in a way accepted by Fortran routines
411
412 dt_dtuout();
413
414}
415
416//______________________________________________________________________________
417Int_t TDPMjet::GetEvNum() const
418{
419 return DTEVT1.nevhkk;
420}
421//______________________________________________________________________________
422Int_t TDPMjet::GetEntriesNum() const
423{
424 return DTEVT1.nhkk;
425}
426//______________________________________________________________________________
427Int_t TDPMjet::GetNumStablePc() const
428{
429 Int_t NumStablePc = 0;
430 for(Int_t i=0; i<DTEVT1.nhkk; i++){
431 if(DTEVT1.isthkk[i] == 1) NumStablePc++;
432 }
433 return NumStablePc;
434}
435
436//______________________________________________________________________________
437Float_t TDPMjet::GetTotEnergy() const
438{
439 Float_t TotEnergy = 0.;
440 for(Int_t i=0; i<DTEVT1.nhkk; i++){
441 if(DTEVT1.isthkk[i] == 1)
442 TotEnergy += DTEVT1.phkk[i][3]; // PHKK[i][3] <-> PHKK(4,i)
443 }
444 return TotEnergy;
445}
446
447//______________________________________________________________________________
448Int_t TDPMjet::GetStatusCode(Int_t evnum) const
449{
450 return DTEVT1.isthkk[evnum];
451}
452//______________________________________________________________________________
453Int_t TDPMjet::GetPDGCode(Int_t evnum) const
454{
455 return DTEVT1.idhkk[evnum];
456}
457//______________________________________________________________________________
458Double_t TDPMjet::Getpx(Int_t evnum) const
459{
460 return DTEVT1.phkk[evnum][0];
461}
462//______________________________________________________________________________
463Double_t TDPMjet::Getpy(Int_t evnum) const
464{
465 return DTEVT1.phkk[evnum][1];
466}
467//______________________________________________________________________________
468Double_t TDPMjet::Getpz(Int_t evnum) const
469{
470 return DTEVT1.phkk[evnum][2];
471}
472//______________________________________________________________________________
473Double_t TDPMjet::GetEnergy(Int_t evnum) const
474{
475 return DTEVT1.phkk[evnum][3];
476}
477//______________________________________________________________________________
478Double_t TDPMjet::GetMass(Int_t evnum) const
479{
480 return DTEVT1.phkk[evnum][4];
481}
482//______________________________________________________________________________
483Int_t TDPMjet::GetFragmentA(Int_t evnum) const
484{
485 return DTEVT2.idres[evnum];
486}
487//______________________________________________________________________________
488Int_t TDPMjet::GetFragmentZ(Int_t evnum) const
489{
490 return DTEVT2.idxres[evnum];
491}
492//______________________________________________________________________________
493Double_t TDPMjet::GetXSFrac() const
494{
495 return DTIMPA.xsfrac;
496}
497//______________________________________________________________________________
498Double_t TDPMjet::GetBImpac() const
499{
500 return DTGLCP.bimpac;
501}
502//______________________________________________________________________________
503Double_t TDPMjet::GetProjRadius() const
504{
505 return DTGLCP.rproj;
506}
507//______________________________________________________________________________
508Double_t TDPMjet::GetTargRadius() const
509{
510 return DTGLCP.rtarg;
511}
512//______________________________________________________________________________
513Int_t TDPMjet::GetProjWounded() const
514{
515 return DTGLCP.nwasam;
516}
517//______________________________________________________________________________
518Int_t TDPMjet::GetTargWounded() const
519{
520 return DTGLCP.nwbsam;
521}
522//______________________________________________________________________________
523Int_t TDPMjet::GetProjSpectators() const
524{
525 return DTGLCP.nwtaac;
526}
527//______________________________________________________________________________
528Int_t TDPMjet::GetTargSpectators() const
529{
530 return DTGLCP.nwtbac;
531}
9dcc0beb 532//______________________________________________________________________________
533Int_t TDPMjet::GetProcessCode() const
534{
535 return POPRCS.iproce;
536}
1a52e0ed 537//______________________________________________________________________________
538void TDPMjet::Dt_Rndm(int idummy)
539{
540 dt_rndm(idummy);
541}
542
543//______________________________________________________________________________
544void TDPMjet::Dt_Rndmst(int na1, int na2, int na3, int nb1)
545{
546 dt_rndmst(na1, na2, na3, nb1);
547}
548
549//______________________________________________________________________________
550void TDPMjet::Dt_Rndmin(int u, int c, int cd, int cm, int i, int j)
551{
552 dt_rndmin(u, c, cd, cm, i, j);
553}
554
555//______________________________________________________________________________
556void TDPMjet::Dt_Rndmou(int u, int c, int cd, int cm, int i, int j)
557{
558 dt_rndmou(u, c, cd, cm, i, j);
559}
560