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