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