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