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