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