]> 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
Corrected storing of number of participants in header.
[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 }
1a52e0ed 176 Int_t nump = 0;
177
178 if(!strcmp(option,"") || !strcmp(option,"Final")){
179 for (Int_t i=0; i < numpart; i++) {
180
181 if (DTEVT1.isthkk[i] == 1) {
182 //
183 // Use the common block values for the TParticle constructor
184 //
1a52e0ed 185 new(Particles[nump]) TParticle(
186 DTEVT1.idhkk[i],
187 DTEVT1.isthkk[i],
188 -1,
189 -1,
190 -1,
191 -1,
192 DTEVT1.phkk[i][0],
193 DTEVT1.phkk[i][1],
194 DTEVT1.phkk[i][2],
195 DTEVT1.phkk[i][3],
196
197 DTEVT1.vhkk[i][0],
198 DTEVT1.vhkk[i][1],
199 DTEVT1.vhkk[i][2],
200 DTEVT1.vhkk[i][3]);
201 nump++;
202 }
203 }
204 }
205 else if(!strcmp(option,"All")){
206 nump = numpart;
207 for (Int_t i=0; i <= numpart; i++){
208
209 // DTEVT1.JMOHKK[i][0] pointer to the entry of the 1st mother of entry i
210 Int_t iParent = DTEVT1.jmohkk[i][0] - 1;
211
212 if(iParent >= 0){
213 TParticle *mother = (TParticle*) (Particles.UncheckedAt(iParent));
1a52e0ed 214 mother->SetLastDaughter(i);
215 if(mother->GetFirstDaughter() == -1) mother->SetFirstDaughter(i);
216 }
217 // --- PDGcode for residual nuclei (idhkk=80000)
218 // --- 10000*Z + 10*A
219 // --- DPMJET -> idres = mass #, idxres = charge
220 if(DTEVT1.idhkk[i] == 80000)
221 DTEVT1.idhkk[i] = 10000*DTEVT2.idxres[i]+10*DTEVT2.idres[i];
d38751f5 222/*
1a52e0ed 223 if(DTEVT2.idxres[i] != 0)
224 printf("\n pc#%d -> A = %d, Z = %d -> PDGcode = %d\n",
225 i,DTEVT2.idres[i],DTEVT2.idxres[i],DTEVT1.idhkk[i]);
d38751f5 226*/
06f923c8 227/*
228 printf("%5d %5d %5d %5d %13.3f %13.3f %13.3f %13.3f \n", i,
229 DTEVT1.idhkk[i],
230 DTEVT1.isthkk[i],
231 iParent,
232 DTEVT1.phkk[i][0],
233 DTEVT1.phkk[i][1],
234 DTEVT1.phkk[i][2],
235 DTEVT1.phkk[i][3]
236 );
237*/
1a52e0ed 238 new(Particles[i]) TParticle(
239 DTEVT1.idhkk[i],
240 DTEVT1.isthkk[i],
241 iParent,
242 -1,
243 -1,
244 -1,
245
246 DTEVT1.phkk[i][0],
247 DTEVT1.phkk[i][1],
248 DTEVT1.phkk[i][2],
249 DTEVT1.phkk[i][3],
250
251 DTEVT1.vhkk[i][0],
252 DTEVT1.vhkk[i][1],
253 DTEVT1.vhkk[i][2],
254 DTEVT1.vhkk[i][3]);
1a52e0ed 255 } // Particle loop
256 }
257 return nump;
fb941ba2 258
1a52e0ed 259}
260
261
262//====================== access to dpmjet subroutines =========================
263//______________________________________________________________________________
264void TDPMjet::Initialize()
265{
d38751f5 266//
267// Write standard DPMJET input cards
268//
269 FILE* out = fopen("dpmjet.inp","w");
270// Projectile and Target definition
cd19c617 271 if (fIp == 1 && fIpz ==1) {
272 fprintf(out, "PROJPAR PROTON\n");
273 } else if (fIp == 1 && fIpz == -1) {
274 fprintf(out, "PROJPAR APROTON\n");
275 } else {
276 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.);
277 }
278
279 if (fIt == 1 && fItz ==1) {
280 fprintf(out, "TARPAR PROTON\n");
281 } else if (fIt == 1 && fItz == -1) {
282 fprintf(out, "TARPAR APROTON\n");
283 } else {
284 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.);
285 }
286
d38751f5 287// Beam energy and crossing-angle
288 fprintf(out, "BEAM %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",fEpn, fEpn, 0., 0., 0., 0.);
289// Centrality
717765ce 290 if (fIp > 1. && fIt > 1)
291 fprintf(out, "CENTRAL %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",-1., fBmin, fBmax, 0., 0., 0.);
f97d4a8e 292// Particle decays
293 if (fPi0Decay)
717765ce 294 fprintf(out, "PARDECAY %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n", 2., 0., 0., 0., 0., 0.);
295
296
d38751f5 297//
298// PHOJET specific
299 fprintf(out, "PHOINPUT\n");
300 fprintf(out, "DEBUG 0 0 0 \n");
6a45631b 301
d38751f5 302 if (fProcess == kDpmMb) {
303 fprintf(out, "PROCESS 1 0 1 1 1 1 1 1\n");
304 } else if (fProcess == kDpmMbNonDiffr) {
cab41d67 305 fprintf(out, "PROCESS 1 0 1 0 0 0 0 1\n");
6a45631b 306 } else if (fProcess == kDpmDiffr) {
cab41d67 307 fprintf(out, "PROCESS 0 0 0 1 1 1 1 0\n");
6a45631b 308 }else if (fProcess == kDpmSingleDiffr) {
309 fprintf(out, "PROCESS 0 0 0 0 1 1 0 0\n");
310 }else if (fProcess == kDpmDoubleDiffr) {
311 fprintf(out, "PROCESS 0 0 0 0 0 0 1 0\n");
cab41d67 312 } else if (fProcess == kDpmCentralDiffr){
313 fprintf(out, "PROCESS 0 0 0 1 0 0 0 0\n");
d38751f5 314 }
315
06f923c8 316 Int_t iPDG[18] =
717765ce 317 {
318// K0s pi0 lam sig+ sig- tet0
319 310, 111, 3122, 3222, 3112, 3322,
320// tet- om- D+ D0 Ds+
321 3312, 3334, 411, 421, 431,
322// etac lamc+ sigc++ sigc+ sigc0 Ksic+
323 441, 4122, 4222, 4212, 4112, 4232,
324// Ksic0 sig0
06f923c8 325 4132
717765ce 326 };
327
328
329 Int_t iON = (fDecayAll) ? 1:0;
06f923c8 330 for (Int_t i = 0; i < 8; i++) {
717765ce 331 fprintf(out, "LUND-DECAY%5d %5d\n", iPDG[i], iON);
332 }
333
334 fprintf(out, "LUND-MSTJ %5d %5d\n", 22, 1);
335
d38751f5 336 fprintf(out, "ENDINPUT\n");
337//
338// START card
339 fprintf(out, "START 1.0 0.0\n");
340 fprintf(out, "STOP\n");
341 fclose(out);
db274076 342 dpmjet_openinp();
343
d38751f5 344//
345// Call DPMJET initialisation
346 Int_t iemu = 0; // No emulsion (default)
347 Dt_Dtuini(1, fEpn, fIp, fIpz, fIt, fItz, fIdp, iemu);
1a52e0ed 348
349}
350
351
352//______________________________________________________________________________
353void TDPMjet::GenerateEvent()
354{
355 // Generates one event;
356 fNEvent++;
357 DTEVNO.nevent=fNEvent;
1a52e0ed 358 Int_t kkmat=-1;
359 Float_t Elab = fEpn;
360 Int_t irej=0;
361 Dt_Kkinc(fIp, fIpz, fIt, fItz, fIdp, Elab, kkmat, irej);
06f923c8 362// dt_evtout(4);
1a52e0ed 363 if(irej!=0) return;
1a52e0ed 364}
365//______________________________________________________________________________
366void TDPMjet::Dt_Dtuini(int nevts, double epn, int npmass, int npchar,
367 int ntmass, int ntchar, int idp, int iemu)
368{
369 // Call dmpjet routine DT_DTUINI passing the parameters
370 // in a way accepted by Fortran routines
371
d38751f5 372
373 printf("\n-------------------------------------------\n");
1a52e0ed 374 printf("\n Dt_Dtuini called with:\n\n");
375 printf(" Projectile -> A = %d, Z = %d \n",npmass, npchar);
376 printf(" Target -> A = %d, Z = %d \n",ntmass, ntchar);
377 printf(" Proj. LAB E -> E = %f GeV \n",epn);
378 printf(" nevts = %d, idp = %d, iemu = %d \n",nevts,idp,iemu);
d38751f5 379 printf("\n-------------------------------------------\n");
1a52e0ed 380
381 dt_dtuini(nevts, epn, npmass, npchar, ntmass, ntchar, idp, iemu);
382
383}
384
385//______________________________________________________________________________
386void TDPMjet::Dt_Kkinc(int npmass, int npchar, int ntmass, int ntchar,
387 int idp, double elab, int kkmat, int irej)
388{
389 // Call dmpjet routine DT_KKINC passing the parameters
390 // in a way accepted by Fortran routines
1a52e0ed 391 dt_kkinc(npmass, npchar, ntmass, ntchar, idp, elab, kkmat, irej);
1a52e0ed 392}
393
394//______________________________________________________________________________
395void TDPMjet::Pho_Phist(int imode, double weight)
396{
397 // Call dmpjet routine PHO_PHIST passing the parameters
398 // in a way accepted by Fortran routines
399
400 pho_phist(imode,weight);
401
402}
403
404//______________________________________________________________________________
405void TDPMjet::Dt_Dtuout()
406{
407 // Call dmpjet routine DT_DTUOT passing the parameters
408 // in a way accepted by Fortran routines
409
410 dt_dtuout();
411
412}
413
414//______________________________________________________________________________
415Int_t TDPMjet::GetEvNum() const
416{
417 return DTEVT1.nevhkk;
418}
419//______________________________________________________________________________
420Int_t TDPMjet::GetEntriesNum() const
421{
422 return DTEVT1.nhkk;
423}
424//______________________________________________________________________________
425Int_t TDPMjet::GetNumStablePc() const
426{
427 Int_t NumStablePc = 0;
428 for(Int_t i=0; i<DTEVT1.nhkk; i++){
429 if(DTEVT1.isthkk[i] == 1) NumStablePc++;
430 }
431 return NumStablePc;
432}
433
434//______________________________________________________________________________
435Float_t TDPMjet::GetTotEnergy() const
436{
437 Float_t TotEnergy = 0.;
438 for(Int_t i=0; i<DTEVT1.nhkk; i++){
439 if(DTEVT1.isthkk[i] == 1)
440 TotEnergy += DTEVT1.phkk[i][3]; // PHKK[i][3] <-> PHKK(4,i)
441 }
442 return TotEnergy;
443}
444
445//______________________________________________________________________________
446Int_t TDPMjet::GetStatusCode(Int_t evnum) const
447{
448 return DTEVT1.isthkk[evnum];
449}
450//______________________________________________________________________________
451Int_t TDPMjet::GetPDGCode(Int_t evnum) const
452{
453 return DTEVT1.idhkk[evnum];
454}
455//______________________________________________________________________________
456Double_t TDPMjet::Getpx(Int_t evnum) const
457{
458 return DTEVT1.phkk[evnum][0];
459}
460//______________________________________________________________________________
461Double_t TDPMjet::Getpy(Int_t evnum) const
462{
463 return DTEVT1.phkk[evnum][1];
464}
465//______________________________________________________________________________
466Double_t TDPMjet::Getpz(Int_t evnum) const
467{
468 return DTEVT1.phkk[evnum][2];
469}
470//______________________________________________________________________________
471Double_t TDPMjet::GetEnergy(Int_t evnum) const
472{
473 return DTEVT1.phkk[evnum][3];
474}
475//______________________________________________________________________________
476Double_t TDPMjet::GetMass(Int_t evnum) const
477{
478 return DTEVT1.phkk[evnum][4];
479}
480//______________________________________________________________________________
481Int_t TDPMjet::GetFragmentA(Int_t evnum) const
482{
483 return DTEVT2.idres[evnum];
484}
485//______________________________________________________________________________
486Int_t TDPMjet::GetFragmentZ(Int_t evnum) const
487{
488 return DTEVT2.idxres[evnum];
489}
490//______________________________________________________________________________
491Double_t TDPMjet::GetXSFrac() const
492{
493 return DTIMPA.xsfrac;
494}
495//______________________________________________________________________________
496Double_t TDPMjet::GetBImpac() const
497{
498 return DTGLCP.bimpac;
499}
500//______________________________________________________________________________
501Double_t TDPMjet::GetProjRadius() const
502{
503 return DTGLCP.rproj;
504}
505//______________________________________________________________________________
506Double_t TDPMjet::GetTargRadius() const
507{
508 return DTGLCP.rtarg;
509}
510//______________________________________________________________________________
511Int_t TDPMjet::GetProjWounded() const
512{
513 return DTGLCP.nwasam;
514}
515//______________________________________________________________________________
516Int_t TDPMjet::GetTargWounded() const
517{
518 return DTGLCP.nwbsam;
519}
3c349853 520
1a52e0ed 521//______________________________________________________________________________
3c349853 522Int_t TDPMjet::GetProjParticipants() const
1a52e0ed 523{
524 return DTGLCP.nwtaac;
525}
526//______________________________________________________________________________
3c349853 527Int_t TDPMjet::GetTargParticipants() const
1a52e0ed 528{
529 return DTGLCP.nwtbac;
530}
9dcc0beb 531//______________________________________________________________________________
532Int_t TDPMjet::GetProcessCode() const
533{
534 return POPRCS.iproce;
535}
1a52e0ed 536//______________________________________________________________________________
537void TDPMjet::Dt_Rndm(int idummy)
538{
539 dt_rndm(idummy);
540}
541
542//______________________________________________________________________________
543void TDPMjet::Dt_Rndmst(int na1, int na2, int na3, int nb1)
544{
545 dt_rndmst(na1, na2, na3, nb1);
546}
547
548//______________________________________________________________________________
549void TDPMjet::Dt_Rndmin(int u, int c, int cd, int cm, int i, int j)
550{
551 dt_rndmin(u, c, cd, cm, i, j);
552}
553
554//______________________________________________________________________________
555void TDPMjet::Dt_Rndmou(int u, int c, int cd, int cm, int i, int j)
556{
557 dt_rndmou(u, c, cd, cm, i, j);
558}
559
fb941ba2 560
561Int_t TDPMjet::NHEP() const {return POEVT1.nhep;}
562Int_t TDPMjet::ISTHEP(Int_t i) const {return POEVT1.isthep[i];}
563Int_t TDPMjet::IDHEP(Int_t i) const {return POEVT1.idhep[i];}
564Int_t TDPMjet::PHEP(Int_t i, Int_t j) const {return POEVT1.phep[i][j];}
565