]> 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
Added debug option
[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),
5fed2e29 119 fProcess(kDpmMb),
120 fFragmentation(kFALSE)
1a52e0ed 121{
f97d4a8e 122// Default Constructor
1a52e0ed 123}
124
125//______________________________________________________________________________
d38751f5 126TDPMjet::TDPMjet(DpmProcess_t iproc, Int_t Ip=208, Int_t Ipz=82, Int_t It=208, Int_t Itz=82,
f97d4a8e 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),
2d22f1fe 135 fPpn(0.),
f97d4a8e 136 fCMEn(CMEn),
137 fIdp(0),
138 fBmin(0.),
139 fBmax(0.),
140 fFCentr(0),
141 fPi0Decay(0),
717765ce 142 fDecayAll(0),
5fed2e29 143 fProcess(iproc),
144 fFragmentation(-1)
1a52e0ed 145{
d38751f5 146 printf("TDPMJet Constructor %d %d %d %d \n", Ip, Ipz, It, Itz);
1a52e0ed 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 }
1a52e0ed 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 //
1a52e0ed 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));
1a52e0ed 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];
d38751f5 224/*
1a52e0ed 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]);
d38751f5 228*/
06f923c8 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*/
1a52e0ed 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]);
1a52e0ed 257 } // Particle loop
258 }
259 return nump;
fb941ba2 260
1a52e0ed 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
e2b2a1dd 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
5fed2e29 296 else if(fIp==1 && fIt>1){ // proton towards A side (directed z>0)
e2b2a1dd 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 }
d38751f5 305// Centrality
5fed2e29 306 if((fIp > 1 || fIt > 1) && fFragmentation)
e2b2a1dd 307 fprintf(out, "CENTRAL %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",-2., fBmin, fBmax, 0., 0., 0.);
5fed2e29 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.);
f97d4a8e 310// Particle decays
311 if (fPi0Decay)
717765ce 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
d38751f5 315//
316// PHOJET specific
317 fprintf(out, "PHOINPUT\n");
318 fprintf(out, "DEBUG 0 0 0 \n");
6a45631b 319
d38751f5 320 if (fProcess == kDpmMb) {
321 fprintf(out, "PROCESS 1 0 1 1 1 1 1 1\n");
322 } else if (fProcess == kDpmMbNonDiffr) {
cab41d67 323 fprintf(out, "PROCESS 1 0 1 0 0 0 0 1\n");
6a45631b 324 } else if (fProcess == kDpmDiffr) {
cab41d67 325 fprintf(out, "PROCESS 0 0 0 1 1 1 1 0\n");
6a45631b 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");
cab41d67 330 } else if (fProcess == kDpmCentralDiffr){
331 fprintf(out, "PROCESS 0 0 0 1 0 0 0 0\n");
d38751f5 332 }
333
06f923c8 334 Int_t iPDG[18] =
717765ce 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
06f923c8 343 4132
717765ce 344 };
345
346
347 Int_t iON = (fDecayAll) ? 1:0;
06f923c8 348 for (Int_t i = 0; i < 8; i++) {
717765ce 349 fprintf(out, "LUND-DECAY%5d %5d\n", iPDG[i], iON);
350 }
351
352 fprintf(out, "LUND-MSTJ %5d %5d\n", 22, 1);
353
d38751f5 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);
db274076 360 dpmjet_openinp();
361
d38751f5 362//
363// Call DPMJET initialisation
364 Int_t iemu = 0; // No emulsion (default)
365 Dt_Dtuini(1, fEpn, fIp, fIpz, fIt, fItz, fIdp, iemu);
1a52e0ed 366
367}
368
369
370//______________________________________________________________________________
371void TDPMjet::GenerateEvent()
372{
373 // Generates one event;
374 fNEvent++;
375 DTEVNO.nevent=fNEvent;
1a52e0ed 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);
1a52e0ed 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
d38751f5 388
389 printf("\n-------------------------------------------\n");
1a52e0ed 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);
d38751f5 395 printf("\n-------------------------------------------\n");
1a52e0ed 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
1a52e0ed 407 dt_kkinc(npmass, npchar, ntmass, ntchar, idp, elab, kkmat, irej);
1a52e0ed 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}
3c349853 536
1a52e0ed 537//______________________________________________________________________________
3c349853 538Int_t TDPMjet::GetProjParticipants() const
1a52e0ed 539{
540 return DTGLCP.nwtaac;
541}
542//______________________________________________________________________________
3c349853 543Int_t TDPMjet::GetTargParticipants() const
1a52e0ed 544{
545 return DTGLCP.nwtbac;
546}
9dcc0beb 547//______________________________________________________________________________
548Int_t TDPMjet::GetProcessCode() const
549{
550 return POPRCS.iproce;
551}
1a52e0ed 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
fb941ba2 576
492d4eba 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];}
fb941ba2 581