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