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