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