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