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