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           printf("%Particle: d %d %d %d \n", i, DTEVT1.idhkk[i],  DTEVT1.isthkk[i], iParent);
 218
 219       } // Particle loop
 220   }
 221   return nump;
 222 }
 223
 224
 225 //====================== access to dpmjet subroutines =========================
 226 //______________________________________________________________________________
 227 void TDPMjet::Initialize()
 228 {
 229 //********************************************************************************
 230 //*Calls DT_DTUINI with the either default parameters or the ones set by the user*
 231 //********************************************************************************
 232
 233 //   printf("\n-------------------------------------------\n");
 234 //   printf("\n         TDPMjet initialized with:\n\n");
 235 //   printf(" Projectile        -> A = %d, Z = %d \n",fIp, fIpz);
 236 //   printf(" Target            -> A = %d, Z = %d \n",fIt, fItz);
 237 //   printf(" Proj. LAB E       -> E = %f GeV \n",fEpn);
 238 //   printf(" CM energy -> Ecm = %f GeV \n",fCMEn);
 239 //   printf("\n-------------------------------------------\n");
 240
 241    Int_t iemu = 0; // No emulsion (default)
 242    Dt_Dtuini(1, fEpn, fIp, fIpz, fIt, fItz, fIdp, iemu);
 243
 244 }
 245
 246
 247 //______________________________________________________________________________
 248 void TDPMjet::GenerateEvent()
 249 {
 250    // Generates one event;
 251    fNEvent++;
 252    DTEVNO.nevent=fNEvent;
 253    //printf("\n fNEvent = %d\n",fNEvent);
 254    Int_t kkmat=-1;
 255    Float_t Elab = fEpn;
 256    Int_t irej=0;
 257    Dt_Kkinc(fIp, fIpz, fIt, fItz, fIdp, Elab, kkmat, irej);
 258    if(irej!=0) return;
 259    //
 260    Int_t imode=2000;
 261    Double_t weight=1.;
 262    Pho_Phist(imode, weight);
 263    Dt_Dtuout();
 264
 265 }
 266 //______________________________________________________________________________
 267 void TDPMjet::Dt_Dtuini(int nevts, double epn, int npmass, int npchar,
 268                         int ntmass, int ntchar, int idp, int iemu)
 269 {
 270   // Call dmpjet routine DT_DTUINI passing the parameters
 271   // in a way accepted by Fortran routines
 272
 273
 274    /*printf("\n-------------------------------------------\n");
 275    printf("\n           Dt_Dtuini called with:\n\n");
 276    printf(" Projectile  -> A = %d, Z = %d \n",npmass, npchar);
 277    printf(" Target      -> A = %d, Z = %d \n",ntmass, ntchar);
 278    printf(" Proj. LAB E -> E = %f GeV \n",epn);
 279    printf(" nevts = %d, idp = %d, iemu = %d \n",nevts,idp,iemu);
 280    printf("\n-------------------------------------------\n");*/
 281
 282    dt_dtuini(nevts, epn, npmass, npchar, ntmass, ntchar, idp, iemu);
 283
 284 }
 285
 286 //______________________________________________________________________________
 287 void TDPMjet::Dt_Kkinc(int npmass, int npchar, int ntmass, int ntchar,
 288                        int idp, double elab, int kkmat, int irej)
 289 {
 290   // Call dmpjet routine DT_KKINC passing the parameters
 291   // in a way accepted by Fortran routines
 292
 293    /*printf("\n-------------------------------------------\n");
 294    printf("\n           Dt_Kkinc called with:\n\n");
 295    printf(" Projectile  -> A = %d, Z = %d \n",npmass, npchar);
 296    printf(" Target      -> A = %d, Z = %d \n",ntmass, ntchar);
 297    printf(" LAB Energy  -> E = %f GeV \n",elab);
 298    printf(" idp = %d,  kkmat = %d, irej = %d \n",idp,kkmat,irej);
 299    printf("\n-------------------------------------------\n");*/
 300
 301   dt_kkinc(npmass, npchar, ntmass, ntchar, idp, elab, kkmat, irej);
 302
 303 }
 304
 305 //______________________________________________________________________________
 306 void TDPMjet::Pho_Phist(int imode, double weight)
 307 {
 308   // Call dmpjet routine PHO_PHIST passing the parameters
 309   // in a way accepted by Fortran routines
 310
 311   pho_phist(imode,weight);
 312
 313 }
 314
 315 //______________________________________________________________________________
 316 void TDPMjet::Dt_Dtuout()
 317 {
 318   // Call dmpjet routine DT_DTUOT passing the parameters
 319   // in a way accepted by Fortran routines
 320
 321   dt_dtuout();
 322
 323 }
 324
 325 //______________________________________________________________________________
 326 Int_t TDPMjet::GetEvNum() const
 327 {
 328         return DTEVT1.nevhkk;
 329 }
 330 //______________________________________________________________________________
 331 Int_t TDPMjet::GetEntriesNum() const
 332 {
 333         return DTEVT1.nhkk;
 334 }
 335 //______________________________________________________________________________
 336 Int_t TDPMjet::GetNumStablePc() const
 337 {
 338         Int_t NumStablePc = 0;
 339         for(Int_t i=0; i<DTEVT1.nhkk; i++){
 340            if(DTEVT1.isthkk[i] == 1) NumStablePc++;
 341         }
 342         return NumStablePc;
 343 }
 344
 345 //______________________________________________________________________________
 346 Float_t TDPMjet::GetTotEnergy() const
 347 {
 348         Float_t TotEnergy = 0.;
 349         for(Int_t i=0; i<DTEVT1.nhkk; i++){
 350           if(DTEVT1.isthkk[i] == 1)
 351             TotEnergy += DTEVT1.phkk[i][3]; // PHKK[i][3] <-> PHKK(4,i)
 352         }
 353         return TotEnergy;
 354 }
 355
 356 //______________________________________________________________________________
 357 Int_t TDPMjet::GetStatusCode(Int_t evnum) const
 358 {
 359         return DTEVT1.isthkk[evnum];
 360 }
 361 //______________________________________________________________________________
 362 Int_t TDPMjet::GetPDGCode(Int_t evnum) const
 363 {
 364         return DTEVT1.idhkk[evnum];
 365 }
 366 //______________________________________________________________________________
 367 Double_t TDPMjet::Getpx(Int_t evnum) const
 368 {
 369         return DTEVT1.phkk[evnum][0];
 370 }
 371 //______________________________________________________________________________
 372 Double_t TDPMjet::Getpy(Int_t evnum) const
 373 {
 374         return DTEVT1.phkk[evnum][1];
 375 }
 376 //______________________________________________________________________________
 377 Double_t TDPMjet::Getpz(Int_t evnum) const
 378 {
 379         return DTEVT1.phkk[evnum][2];
 380 }
 381 //______________________________________________________________________________
 382 Double_t TDPMjet::GetEnergy(Int_t evnum) const
 383 {
 384         return DTEVT1.phkk[evnum][3];
 385 }
 386 //______________________________________________________________________________
 387 Double_t TDPMjet::GetMass(Int_t evnum) const
 388 {
 389         return DTEVT1.phkk[evnum][4];
 390 }
 391 //______________________________________________________________________________
 392 Int_t    TDPMjet::GetFragmentA(Int_t evnum) const
 393 {
 394         return DTEVT2.idres[evnum];
 395 }
 396 //______________________________________________________________________________
 397 Int_t    TDPMjet::GetFragmentZ(Int_t evnum) const
 398 {
 399         return DTEVT2.idxres[evnum];
 400 }
 401 //______________________________________________________________________________
 402 Double_t TDPMjet::GetXSFrac() const
 403 {
 404         return DTIMPA.xsfrac;
 405 }
 406 //______________________________________________________________________________
 407 Double_t TDPMjet::GetBImpac() const
 408 {
 409         return DTGLCP.bimpac;
 410 }
 411 //______________________________________________________________________________
 412 Double_t TDPMjet::GetProjRadius() const
 413 {
 414         return DTGLCP.rproj;
 415 }
 416 //______________________________________________________________________________
 417 Double_t TDPMjet::GetTargRadius() const
 418 {
 419         return DTGLCP.rtarg;
 420 }
 421 //______________________________________________________________________________
 422 Int_t TDPMjet::GetProjWounded() const
 423 {
 424         return DTGLCP.nwasam;
 425 }
 426 //______________________________________________________________________________
 427 Int_t TDPMjet::GetTargWounded() const
 428 {
 429         return DTGLCP.nwbsam;
 430 }
 431 //______________________________________________________________________________
 432 Int_t TDPMjet::GetProjSpectators() const
 433 {
 434         return DTGLCP.nwtaac;
 435 }
 436 //______________________________________________________________________________
 437 Int_t TDPMjet::GetTargSpectators() const
 438 {
 439         return DTGLCP.nwtbac;
 440 }
 441
 442 //______________________________________________________________________________
 443 void TDPMjet::Dt_Rndm(int idummy)
 444 {
 445         dt_rndm(idummy);
 446 }
 447
 448 //______________________________________________________________________________
 449 void TDPMjet::Dt_Rndmst(int na1, int na2, int na3, int nb1)
 450 {
 451         dt_rndmst(na1, na2, na3, nb1);
 452 }
 453
 454 //______________________________________________________________________________
 455 void TDPMjet::Dt_Rndmin(int u, int c, int cd, int cm, int i, int j)
 456 {
 457         dt_rndmin(u, c, cd, cm, i, j);
 458 }
 459
 460 //______________________________________________________________________________
 461 void TDPMjet::Dt_Rndmou(int u, int c, int cd, int cm, int i, int j)
 462 {
 463         dt_rndmou(u, c, cd, cm, i, j);
 464 }
 465