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