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