More user settings for decays.
[u/mrichter/AliRoot.git] / TDPMjet / TDPMjet.cxx
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 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
73 # define dpmjet_openinp DPMJET_OPENINP
74 # define type_of_call _stdcall
75 #endif
76
77 #ifndef WIN32
78 extern "C" void   type_of_call dt_dtuini(Int_t & , Double_t &, Int_t & , Int_t &, 
79                         Int_t &, Int_t &, Int_t &, Int_t &);
80 extern "C" double type_of_call dt_getemu(Int_t &, Int_t &, Int_t &, Int_t &);
81 extern "C" void   type_of_call dt_kkinc(Int_t &, Int_t &, Int_t &, Int_t &,
82                                     Int_t &, Double_t &, Int_t &, Int_t &);
83 extern "C" void   type_of_call pho_phist(Int_t &, Double_t &);
84 extern "C" void   type_of_call dt_dtuout();
85 extern "C" void   type_of_call dt_rndm(Int_t &);
86 extern "C" void   type_of_call dt_rndmst(Int_t &, Int_t &, Int_t &, Int_t &);
87 extern "C" void   type_of_call dt_rndmin(Int_t &, Int_t &, Int_t &, Int_t &, Int_t &, Int_t &);
88 extern "C" void   type_of_call dt_rndmou(Int_t &, Int_t &, Int_t &, Int_t &, Int_t &, Int_t &);
89 extern "C" void   type_of_call dpmjet_openinp();
90
91 #else
92
93 #endif
94
95 ClassImp(TDPMjet)
96
97
98 //______________________________________________________________________________
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),
114         fDecayAll(0),
115         fProcess(kDpmMb)
116 {
117 // Default Constructor
118 }
119
120 //______________________________________________________________________________
121 TDPMjet::TDPMjet(DpmProcess_t  iproc, Int_t Ip=208, Int_t Ipz=82, Int_t It=208, Int_t Itz=82, 
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),
130       fPpn(0.),
131       fCMEn(CMEn),
132       fIdp(0),
133       fBmin(0.),
134       fBmax(0.),
135       fFCentr(0),
136       fPi0Decay(0),
137       fDecayAll(0),
138       fProcess(iproc)
139 {  
140     printf("TDPMJet Constructor %d %d %d %d \n", Ip, Ipz, It, Itz);
141 }
142
143
144 //______________________________________________________________________________
145 Int_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   }
172   //printf("\n TDPMjet: DPMJET stack contains %d particles", numpart);
173   // printf("\n TDPMjet: Final not decayed particles: %d",    numStabpart);
174   //printf("\n TDPMjet: Total energy: %f GeV          \n",   entot);
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            //
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));
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];
221 /*
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]);
225 */        
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]);
243       } // Particle loop  
244   }
245   return nump;
246 }
247
248
249 //====================== access to dpmjet subroutines =========================
250 //______________________________________________________________________________
251 void TDPMjet::Initialize()
252 {
253 //
254 //  Write standard DPMJET input cards 
255 //
256     FILE* out = fopen("dpmjet.inp","w");
257 //  Projectile and Target definition 
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
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
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.);
279 //  Particle decays
280     if (fPi0Decay) 
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     
284 //
285 //  PHOJET specific
286     fprintf(out, "PHOINPUT\n");
287     fprintf(out, "DEBUG      0 0 0 \n");
288
289     if (fProcess == kDpmMb) {
290         fprintf(out, "PROCESS           1 0 1 1 1 1 1 1\n");
291     } else if (fProcess == kDpmMbNonDiffr) {
292         fprintf(out, "PROCESS           1 0 1 1 0 0 0 1\n");
293     } else if (fProcess == kDpmDiffr) {
294         fprintf(out, "PROCESS           0 0 0 0 1 1 1 0\n");
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");
299     }
300     
301     Int_t iPDG[19] = 
302         {
303 //          K0s   pi0  lam   sig+  sig-  tet0
304         310,  111, 3122, 3222, 3112, 3322,
305 //          tet- om-    D+      D0     Ds+
306         3312, 3334,  411,  421,  431,
307 //          etac lamc+ sigc++ sigc+ sigc0 Ksic+
308         441, 4122, 4222, 4212, 4112, 4232,
309 //         Ksic0 sig0 
310         4132, 3212
311         };
312     
313     
314     Int_t iON = (fDecayAll) ? 1:0;
315     for (Int_t i = 0; i < 19; i++) {
316         fprintf(out, "LUND-DECAY%5d %5d\n",  iPDG[i], iON);    
317     }
318         
319     fprintf(out, "LUND-MSTJ %5d %5d\n",   22, 1);    
320
321     fprintf(out, "ENDINPUT\n");
322 //
323 //  START card
324     fprintf(out, "START            1.0       0.0\n");
325     fprintf(out, "STOP\n");
326     fclose(out);
327     dpmjet_openinp();
328
329 //
330 //  Call DPMJET initialisation
331     Int_t iemu = 0; // No emulsion (default)
332     Dt_Dtuini(1, fEpn, fIp, fIpz, fIt, fItz, fIdp, iemu);
333
334 }
335
336
337 //______________________________________________________________________________
338 void TDPMjet::GenerateEvent()
339 {
340    // Generates one event;
341    fNEvent++;
342    DTEVNO.nevent=fNEvent;
343    Int_t kkmat=-1;
344    Float_t Elab = fEpn;
345    Int_t irej=0;
346    Dt_Kkinc(fIp, fIpz, fIt, fItz, fIdp, Elab, kkmat, irej);
347    if(irej!=0) return;
348 }
349 //______________________________________________________________________________
350 void TDPMjet::Dt_Dtuini(int nevts, double epn, int npmass, int npchar, 
351                         int ntmass, int ntchar, int idp, int iemu)
352 {
353   // Call dmpjet routine DT_DTUINI passing the parameters 
354   // in a way accepted by Fortran routines                                 
355      
356
357    printf("\n-------------------------------------------\n");
358    printf("\n           Dt_Dtuini called with:\n\n");
359    printf(" Projectile  -> A = %d, Z = %d \n",npmass, npchar);
360    printf(" Target      -> A = %d, Z = %d \n",ntmass, ntchar);
361    printf(" Proj. LAB E -> E = %f GeV \n",epn);
362    printf(" nevts = %d, idp = %d, iemu = %d \n",nevts,idp,iemu);
363    printf("\n-------------------------------------------\n");
364
365    dt_dtuini(nevts, epn, npmass, npchar, ntmass, ntchar, idp, iemu);
366     
367 }
368
369 //______________________________________________________________________________
370 void TDPMjet::Dt_Kkinc(int npmass, int npchar, int ntmass, int ntchar, 
371                        int idp, double elab, int kkmat, int irej)
372 {
373   // Call dmpjet routine DT_KKINC passing the parameters 
374   // in a way accepted by Fortran routines                                 
375   dt_kkinc(npmass, npchar, ntmass, ntchar, idp, elab, kkmat, irej);
376 }
377
378 //______________________________________________________________________________
379 void TDPMjet::Pho_Phist(int imode, double weight)
380 {
381   // Call dmpjet routine PHO_PHIST passing the parameters 
382   // in a way accepted by Fortran routines
383   
384   pho_phist(imode,weight);                                 
385
386 }
387
388 //______________________________________________________________________________
389 void TDPMjet::Dt_Dtuout()
390 {
391   // Call dmpjet routine DT_DTUOT passing the parameters 
392   // in a way accepted by Fortran routines                                 
393   
394   dt_dtuout();
395
396 }
397
398 //______________________________________________________________________________
399 Int_t TDPMjet::GetEvNum() const
400 {
401         return DTEVT1.nevhkk;
402 }
403 //______________________________________________________________________________
404 Int_t TDPMjet::GetEntriesNum() const
405 {
406         return DTEVT1.nhkk;
407 }
408 //______________________________________________________________________________
409 Int_t TDPMjet::GetNumStablePc() const
410 {
411         Int_t NumStablePc = 0;
412         for(Int_t i=0; i<DTEVT1.nhkk; i++){
413            if(DTEVT1.isthkk[i] == 1) NumStablePc++;
414         }
415         return NumStablePc;
416 }
417
418 //______________________________________________________________________________
419 Float_t TDPMjet::GetTotEnergy() const
420 {
421         Float_t TotEnergy = 0.;
422         for(Int_t i=0; i<DTEVT1.nhkk; i++){
423           if(DTEVT1.isthkk[i] == 1)
424             TotEnergy += DTEVT1.phkk[i][3]; // PHKK[i][3] <-> PHKK(4,i)
425         }
426         return TotEnergy;
427 }
428
429 //______________________________________________________________________________
430 Int_t TDPMjet::GetStatusCode(Int_t evnum) const 
431 {
432         return DTEVT1.isthkk[evnum];    
433 }
434 //______________________________________________________________________________
435 Int_t TDPMjet::GetPDGCode(Int_t evnum) const   
436 {
437         return DTEVT1.idhkk[evnum];
438 }
439 //______________________________________________________________________________
440 Double_t TDPMjet::Getpx(Int_t evnum) const       
441 {
442         return DTEVT1.phkk[evnum][0];
443 }
444 //______________________________________________________________________________
445 Double_t TDPMjet::Getpy(Int_t evnum) const      
446 {
447         return DTEVT1.phkk[evnum][1];
448 }
449 //______________________________________________________________________________
450 Double_t TDPMjet::Getpz(Int_t evnum) const       
451 {
452         return DTEVT1.phkk[evnum][2];
453 }
454 //______________________________________________________________________________
455 Double_t TDPMjet::GetEnergy(Int_t evnum) const        
456 {
457         return DTEVT1.phkk[evnum][3];
458 }
459 //______________________________________________________________________________
460 Double_t TDPMjet::GetMass(Int_t evnum) const          
461 {
462         return DTEVT1.phkk[evnum][4];
463 }
464 //______________________________________________________________________________
465 Int_t    TDPMjet::GetFragmentA(Int_t evnum) const       
466 {
467         return DTEVT2.idres[evnum];
468 }
469 //______________________________________________________________________________
470 Int_t    TDPMjet::GetFragmentZ(Int_t evnum) const       
471 {
472         return DTEVT2.idxres[evnum];
473 }
474 //______________________________________________________________________________
475 Double_t TDPMjet::GetXSFrac() const           
476 {
477         return DTIMPA.xsfrac;
478 }
479 //______________________________________________________________________________
480 Double_t TDPMjet::GetBImpac() const           
481 {
482         return DTGLCP.bimpac;
483 }
484 //______________________________________________________________________________
485 Double_t TDPMjet::GetProjRadius() const               
486 {
487         return DTGLCP.rproj;
488 }
489 //______________________________________________________________________________
490 Double_t TDPMjet::GetTargRadius() const               
491 {
492         return DTGLCP.rtarg;
493 }
494 //______________________________________________________________________________
495 Int_t TDPMjet::GetProjWounded() const
496 {
497         return DTGLCP.nwasam;
498 }
499 //______________________________________________________________________________
500 Int_t TDPMjet::GetTargWounded() const
501 {
502         return DTGLCP.nwbsam;
503 }
504 //______________________________________________________________________________
505 Int_t TDPMjet::GetProjSpectators() const
506 {
507         return DTGLCP.nwtaac;
508 }
509 //______________________________________________________________________________
510 Int_t TDPMjet::GetTargSpectators() const
511 {
512         return DTGLCP.nwtbac;
513 }
514 //______________________________________________________________________________
515 Int_t TDPMjet::GetProcessCode() const
516 {
517                 return POPRCS.iproce;
518 }
519 //______________________________________________________________________________
520 void TDPMjet::Dt_Rndm(int idummy)
521 {
522         dt_rndm(idummy);
523 }
524
525 //______________________________________________________________________________
526 void TDPMjet::Dt_Rndmst(int na1, int na2, int na3, int nb1)
527 {
528         dt_rndmst(na1, na2, na3, nb1);
529 }
530
531 //______________________________________________________________________________
532 void TDPMjet::Dt_Rndmin(int u, int c, int cd, int cm, int i, int j)
533 {
534         dt_rndmin(u, c, cd, cm, i, j);
535 }
536
537 //______________________________________________________________________________
538 void TDPMjet::Dt_Rndmou(int u, int c, int cd, int cm, int i, int j)
539 {
540         dt_rndmou(u, c, cd, cm, i, j);
541 }
542