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