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   printf("\n TDPMjet: DPMJET stack contains %d particles", numpart);
 177   // printf("\n TDPMjet: Final not decayed particles: %d",    numStabpart);
 178   //printf("\n TDPMjet: Total energy: %f GeV          \n",   entot);
 179   Int_t nump = 0;
 180
 181   if(!strcmp(option,"") || !strcmp(option,"Final")){
 182       for (Int_t i=0; i < numpart; i++) {
 183
 184           if (DTEVT1.isthkk[i] == 1) {
 185            //
 186            //  Use the common block values for the TParticle constructor
 187            //
 188               new(Particles[nump]) TParticle(
 189                   DTEVT1.idhkk[i],
 190                   DTEVT1.isthkk[i],
 191                   -1,
 192                   -1,
 193                   -1,
 194                   -1,
 195                   DTEVT1.phkk[i][0],
 196                   DTEVT1.phkk[i][1],
 197                   DTEVT1.phkk[i][2],
 198                   DTEVT1.phkk[i][3],
 199
 200                   DTEVT1.vhkk[i][0],
 201                   DTEVT1.vhkk[i][1],
 202                   DTEVT1.vhkk[i][2],
 203                   DTEVT1.vhkk[i][3]);
 204               nump++;
 205           }
 206       }
 207   }
 208   else if(!strcmp(option,"All")){
 209       nump = numpart;
 210       for (Int_t i=0; i <= numpart; i++){
 211
 212           // DTEVT1.JMOHKK[i][0] pointer to the entry of the 1st mother of entry i
 213           Int_t iParent = DTEVT1.jmohkk[i][0] - 1;
 214
 215           if(iParent >= 0){
 216               TParticle *mother = (TParticle*) (Particles.UncheckedAt(iParent));
 217               mother->SetLastDaughter(i);
 218               if(mother->GetFirstDaughter() == -1) mother->SetFirstDaughter(i);
 219           }
 220           // --- PDGcode for residual nuclei (idhkk=80000)
 221           // ---        10000*Z + 10*A
 222           // --- DPMJET -> idres = mass #, idxres = charge
 223           if(DTEVT1.idhkk[i] == 80000)
 224               DTEVT1.idhkk[i] = 10000*DTEVT2.idxres[i]+10*DTEVT2.idres[i];
 225 /*
 226           if(DTEVT2.idxres[i] != 0)
 227               printf("\n        pc#%d -> A = %d, Z = %d -> PDGcode = %d\n",
 228                      i,DTEVT2.idres[i],DTEVT2.idxres[i],DTEVT1.idhkk[i]);
 229 */
 230 /*
 231           printf("%5d %5d %5d %5d %13.3f %13.3f %13.3f %13.3f \n", i,
 232                  DTEVT1.idhkk[i],
 233                  DTEVT1.isthkk[i],
 234                  iParent,
 235                  DTEVT1.phkk[i][0],
 236                  DTEVT1.phkk[i][1],
 237                  DTEVT1.phkk[i][2],
 238                  DTEVT1.phkk[i][3]
 239                  );
 240 */
 241           new(Particles[i]) TParticle(
 242               DTEVT1.idhkk[i],
 243               DTEVT1.isthkk[i],
 244               iParent,
 245               -1,
 246               -1,
 247               -1,
 248
 249               DTEVT1.phkk[i][0],
 250               DTEVT1.phkk[i][1],
 251               DTEVT1.phkk[i][2],
 252               DTEVT1.phkk[i][3],
 253
 254               DTEVT1.vhkk[i][0],
 255               DTEVT1.vhkk[i][1],
 256               DTEVT1.vhkk[i][2],
 257               DTEVT1.vhkk[i][3]);
 258       } // Particle loop
 259   }
 260   return nump;
 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, "BEAM      %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",fEpn, fEpn, 0., 0., 0., 0.);
 291 //  Centrality
 292     if (fIp > 1. && fIt > 1)
 293         fprintf(out, "CENTRAL   %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",-1., fBmin, fBmax, 0., 0., 0.);
 294 //  Particle decays
 295     if (fPi0Decay)
 296         fprintf(out, "PARDECAY  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n", 2., 0., 0., 0., 0., 0.);
 297
 298
 299 //
 300 //  PHOJET specific
 301     fprintf(out, "PHOINPUT\n");
 302     fprintf(out, "DEBUG      0 0 0 \n");
 303
 304     if (fProcess == kDpmMb) {
 305         fprintf(out, "PROCESS           1 0 1 1 1 1 1 1\n");
 306     } else if (fProcess == kDpmMbNonDiffr) {
 307         fprintf(out, "PROCESS           1 0 1 0 0 0 0 1\n");
 308     } else if (fProcess == kDpmDiffr) {
 309         fprintf(out, "PROCESS           0 0 0 1 1 1 1 0\n");
 310     }else if (fProcess == kDpmSingleDiffr) {
 311         fprintf(out, "PROCESS           0 0 0 0 1 1 0 0\n");
 312     }else if (fProcess == kDpmDoubleDiffr) {
 313         fprintf(out, "PROCESS           0 0 0 0 0 0 1 0\n");
 314     } else if (fProcess == kDpmCentralDiffr){
 315         fprintf(out, "PROCESS           0 0 0 1 0 0 0 0\n");
 316     }
 317
 318     Int_t iPDG[18] =
 319         {
 320 //          K0s   pi0  lam   sig+  sig-  tet0
 321         310,  111, 3122, 3222, 3112, 3322,
 322 //          tet- om-    D+      D0     Ds+
 323         3312, 3334,  411,  421,  431,
 324 //          etac lamc+ sigc++ sigc+ sigc0 Ksic+
 325         441, 4122, 4222, 4212, 4112, 4232,
 326 //         Ksic0 sig0
 327         4132
 328         };
 329
 330
 331     Int_t iON = (fDecayAll) ? 1:0;
 332     for (Int_t i = 0; i < 8; i++) {
 333         fprintf(out, "LUND-DECAY%5d %5d\n",  iPDG[i], iON);
 334     }
 335
 336     fprintf(out, "LUND-MSTJ %5d %5d\n",   22, 1);
 337
 338     fprintf(out, "ENDINPUT\n");
 339 //
 340 //  START card
 341     fprintf(out, "START            1.0       0.0\n");
 342     fprintf(out, "STOP\n");
 343     fclose(out);
 344     dpmjet_openinp();
 345
 346 //
 347 //  Call DPMJET initialisation
 348     Int_t iemu = 0; // No emulsion (default)
 349     Dt_Dtuini(1, fEpn, fIp, fIpz, fIt, fItz, fIdp, iemu);
 350
 351 }
 352
 353
 354 //______________________________________________________________________________
 355 void TDPMjet::GenerateEvent()
 356 {
 357    // Generates one event;
 358    fNEvent++;
 359    DTEVNO.nevent=fNEvent;
 360    Int_t kkmat=-1;
 361    Float_t Elab = fEpn;
 362    Int_t irej=0;
 363    Dt_Kkinc(fIp, fIpz, fIt, fItz, fIdp, Elab, kkmat, irej);
 364 //   dt_evtout(4);
 365    if(irej!=0) return;
 366 }
 367 //______________________________________________________________________________
 368 void TDPMjet::Dt_Dtuini(int nevts, double epn, int npmass, int npchar,
 369                         int ntmass, int ntchar, int idp, int iemu)
 370 {
 371   // Call dmpjet routine DT_DTUINI passing the parameters
 372   // in a way accepted by Fortran routines
 373
 374
 375    printf("\n-------------------------------------------\n");
 376    printf("\n           Dt_Dtuini called with:\n\n");
 377    printf(" Projectile  -> A = %d, Z = %d \n",npmass, npchar);
 378    printf(" Target      -> A = %d, Z = %d \n",ntmass, ntchar);
 379    printf(" Proj. LAB E -> E = %f GeV \n",epn);
 380    printf(" nevts = %d, idp = %d, iemu = %d \n",nevts,idp,iemu);
 381    printf("\n-------------------------------------------\n");
 382
 383    dt_dtuini(nevts, epn, npmass, npchar, ntmass, ntchar, idp, iemu);
 384
 385 }
 386
 387 //______________________________________________________________________________
 388 void TDPMjet::Dt_Kkinc(int npmass, int npchar, int ntmass, int ntchar,
 389                        int idp, double elab, int kkmat, int irej)
 390 {
 391   // Call dmpjet routine DT_KKINC passing the parameters
 392   // in a way accepted by Fortran routines
 393   dt_kkinc(npmass, npchar, ntmass, ntchar, idp, elab, kkmat, irej);
 394 }
 395
 396 //______________________________________________________________________________
 397 void TDPMjet::Pho_Phist(int imode, double weight)
 398 {
 399   // Call dmpjet routine PHO_PHIST passing the parameters
 400   // in a way accepted by Fortran routines
 401
 402   pho_phist(imode,weight);
 403
 404 }
 405
 406 //______________________________________________________________________________
 407 void TDPMjet::Dt_Dtuout()
 408 {
 409   // Call dmpjet routine DT_DTUOT passing the parameters
 410   // in a way accepted by Fortran routines
 411
 412   dt_dtuout();
 413
 414 }
 415
 416 //______________________________________________________________________________
 417 Int_t TDPMjet::GetEvNum() const
 418 {
 419         return DTEVT1.nevhkk;
 420 }
 421 //______________________________________________________________________________
 422 Int_t TDPMjet::GetEntriesNum() const
 423 {
 424         return DTEVT1.nhkk;
 425 }
 426 //______________________________________________________________________________
 427 Int_t TDPMjet::GetNumStablePc() const
 428 {
 429         Int_t NumStablePc = 0;
 430         for(Int_t i=0; i<DTEVT1.nhkk; i++){
 431            if(DTEVT1.isthkk[i] == 1) NumStablePc++;
 432         }
 433         return NumStablePc;
 434 }
 435
 436 //______________________________________________________________________________
 437 Float_t TDPMjet::GetTotEnergy() const
 438 {
 439         Float_t TotEnergy = 0.;
 440         for(Int_t i=0; i<DTEVT1.nhkk; i++){
 441           if(DTEVT1.isthkk[i] == 1)
 442             TotEnergy += DTEVT1.phkk[i][3]; // PHKK[i][3] <-> PHKK(4,i)
 443         }
 444         return TotEnergy;
 445 }
 446
 447 //______________________________________________________________________________
 448 Int_t TDPMjet::GetStatusCode(Int_t evnum) const
 449 {
 450         return DTEVT1.isthkk[evnum];
 451 }
 452 //______________________________________________________________________________
 453 Int_t TDPMjet::GetPDGCode(Int_t evnum) const
 454 {
 455         return DTEVT1.idhkk[evnum];
 456 }
 457 //______________________________________________________________________________
 458 Double_t TDPMjet::Getpx(Int_t evnum) const
 459 {
 460         return DTEVT1.phkk[evnum][0];
 461 }
 462 //______________________________________________________________________________
 463 Double_t TDPMjet::Getpy(Int_t evnum) const
 464 {
 465         return DTEVT1.phkk[evnum][1];
 466 }
 467 //______________________________________________________________________________
 468 Double_t TDPMjet::Getpz(Int_t evnum) const
 469 {
 470         return DTEVT1.phkk[evnum][2];
 471 }
 472 //______________________________________________________________________________
 473 Double_t TDPMjet::GetEnergy(Int_t evnum) const
 474 {
 475         return DTEVT1.phkk[evnum][3];
 476 }
 477 //______________________________________________________________________________
 478 Double_t TDPMjet::GetMass(Int_t evnum) const
 479 {
 480         return DTEVT1.phkk[evnum][4];
 481 }
 482 //______________________________________________________________________________
 483 Int_t    TDPMjet::GetFragmentA(Int_t evnum) const
 484 {
 485         return DTEVT2.idres[evnum];
 486 }
 487 //______________________________________________________________________________
 488 Int_t    TDPMjet::GetFragmentZ(Int_t evnum) const
 489 {
 490         return DTEVT2.idxres[evnum];
 491 }
 492 //______________________________________________________________________________
 493 Double_t TDPMjet::GetXSFrac() const
 494 {
 495         return DTIMPA.xsfrac;
 496 }
 497 //______________________________________________________________________________
 498 Double_t TDPMjet::GetBImpac() const
 499 {
 500         return DTGLCP.bimpac;
 501 }
 502 //______________________________________________________________________________
 503 Double_t TDPMjet::GetProjRadius() const
 504 {
 505         return DTGLCP.rproj;
 506 }
 507 //______________________________________________________________________________
 508 Double_t TDPMjet::GetTargRadius() const
 509 {
 510         return DTGLCP.rtarg;
 511 }
 512 //______________________________________________________________________________
 513 Int_t TDPMjet::GetProjWounded() const
 514 {
 515         return DTGLCP.nwasam;
 516 }
 517 //______________________________________________________________________________
 518 Int_t TDPMjet::GetTargWounded() const
 519 {
 520         return DTGLCP.nwbsam;
 521 }
 522 //______________________________________________________________________________
 523 Int_t TDPMjet::GetProjSpectators() const
 524 {
 525         return DTGLCP.nwtaac;
 526 }
 527 //______________________________________________________________________________
 528 Int_t TDPMjet::GetTargSpectators() const
 529 {
 530         return DTGLCP.nwtbac;
 531 }
 532 //______________________________________________________________________________
 533 Int_t TDPMjet::GetProcessCode() const
 534 {
 535                 return POPRCS.iproce;
 536 }
 537 //______________________________________________________________________________
 538 void TDPMjet::Dt_Rndm(int idummy)
 539 {
 540         dt_rndm(idummy);
 541 }
 542
 543 //______________________________________________________________________________
 544 void TDPMjet::Dt_Rndmst(int na1, int na2, int na3, int nb1)
 545 {
 546         dt_rndmst(na1, na2, na3, nb1);
 547 }
 548
 549 //______________________________________________________________________________
 550 void TDPMjet::Dt_Rndmin(int u, int c, int cd, int cm, int i, int j)
 551 {
 552         dt_rndmin(u, c, cd, cm, i, j);
 553 }
 554
 555 //______________________________________________________________________________
 556 void TDPMjet::Dt_Rndmou(int u, int c, int cd, int cm, int i, int j)
 557 {
 558         dt_rndmou(u, c, cd, cm, i, j);
 559 }
 560