THerwig/THerwig6.cxx

   1 // definition of c++ Class THerwig6 to be used in ROOT
   2 // this is a c++ interface to the F77 Herwig6 program
   3 // author: j. g. contreras jgcn@moni.mda.cinvestav.mx
   4 // date: december 22, 2000
   5 //  Class THerwig6 is an interface to the Herwig program
   6 //
   7 // C-----------------------------------------------------------------------
   8 // C                           H E R W I G
   9 // C
  10 // C            a Monte Carlo event generator for simulating
  11 // C        +---------------------------------------------------+
  12 // C        | Hadron Emission Reactions With Interfering Gluons |
  13 // C        +---------------------------------------------------+
  14 // C I.G. Knowles(*), G. Marchesini(+), M.H. Seymour($) and B.R. Webber(#)
  15 // C-----------------------------------------------------------------------
  16 // C with Minimal Supersymmetric Standard Model Matrix Elements by
  17 // C                  S. Moretti($) and K. Odagiri($)
  18 // C-----------------------------------------------------------------------
  19 // C R parity violating Supersymmetric Decays and Matrix Elements by
  20 // C                          P. Richardson(&)
  21 // C-----------------------------------------------------------------------
  22 // C matrix element corrections to top decay and Drell-Yan type processes
  23 // C                         by G. Corcella(+)
  24 // C-----------------------------------------------------------------------
  25 // C Deep Inelastic Scattering and Heavy Flavour Electroproduction by
  26 // C                  G. Abbiendi(@) and L. Stanco(%)
  27 // C-----------------------------------------------------------------------
  28 // C and Jet Photoproduction in Lepton-Hadron Collisions by J. Chyla(~)
  29 // C-----------------------------------------------------------------------
  30 // C(*)  Department of Physics & Astronomy, University of Edinburgh
  31 // C(+)  Dipartimento di Fisica, Universita di Milano
  32 // C($)  Rutherford Appleton Laboratory
  33 // C(#)  Cavendish Laboratory, Cambridge
  34 // C(&)  Department of Physics, University of Oxford
  35 // C(@)  Dipartimento di Fisica, Universita di Bologna
  36 // C(%)  Dipartimento di Fisica, Universita di Padova
  37 // C(~)  Institute of Physics, Prague
  38 // C-----------------------------------------------------------------------
  39 // C                  Version 6.100 - 16th December 1999
  40 // C-----------------------------------------------------------------------
  41 // C Main reference:
  42 // C    G.Marchesini,  B.R.Webber,  G.Abbiendi,  I.G.Knowles,  M.H.Seymour,
  43 // C    and L.Stanco, Computer Physics Communications 67 (1992) 465.
  44 // C-----------------------------------------------------------------------
  45 // C Please send e-mail about  this program  to one of the  authors at the
  46 // C following Internet addresses:
  47 // C    I.Knowles@ed.ac.uk        Giuseppe.Marchesini@mi.infn.it
  48 // C    M.Seymour@rl.ac.uk        webber@hep.phy.cam.ac.uk
  49 // C-----------------------------------------------------------------------
  50
  51
  52 #include "THerwig6.h"
  53 #include "TClonesArray.h"
  54 #include "TParticle.h"
  55 #include "TObjArray.h"
  56
  57
  58 ClassImp(THerwig6)
  59
  60 extern "C" {
  61   void*  herwig6_common_block_address_(char*, int len);
  62   void   herwig6_open_fortran_file_ (int* lun, char* name, int);
  63   void   herwig6_close_fortran_file_(int* lun);
  64 }
  65
  66
  67 THerwig6::THerwig6() : TGenerator("Herwig6","Herwig6") {
  68
  69 // THerwig6 constructor: creates a TClonesArray in which it will store all
  70 // particles. Note that there may be only one functional THerwig6 object
  71 // at a time, so it's not use to create more than one instance of it.
  72
  73   delete fParticles; // was allocated as TObjArray in TGenerator
  74
  75   fParticles = new TClonesArray("TParticle",50);
  76
  77   // initialize common-blocks
  78   fHepevt = (Hepevt_t*) herwig6_common_block_address_((char*)"HEPEVT",6);
  79   fHwbeam = (Hwbeam_t*) herwig6_common_block_address_((char*)"HWBEAM",6);
  80   fHwbmch = (Hwbmch_t*) herwig6_common_block_address_((char*)"HWBMCH",6);
  81   fHwproc = (Hwproc_t*) herwig6_common_block_address_((char*)"HWPROC",6);
  82   fHwpram = (Hwpram_t*) herwig6_common_block_address_((char*)"HWPRAM",6);
  83   fHwprch = (Hwprch_t*) herwig6_common_block_address_((char*)"HWPRCH",6);
  84   fHwpart = (Hwpart_t*) herwig6_common_block_address_((char*)"HWPART",6);
  85   fHwparp = (Hwparp_t*) herwig6_common_block_address_((char*)"HWPARP",6);
  86   fHwbosc = (Hwbosc_t*) herwig6_common_block_address_((char*)"HWBOSC",6);
  87   fHwparc = (Hwparc_t*) herwig6_common_block_address_((char*)"HWPARC",6);
  88   fHwbrch = (Hwbrch_t*) herwig6_common_block_address_((char*)"HWBRCH",6);
  89   fHwevnt = (Hwevnt_t*) herwig6_common_block_address_((char*)"HWEVNT",6);
  90   fHwhard = (Hwhard_t*) herwig6_common_block_address_((char*)"HWHARD",6);
  91   fHwprop = (Hwprop_t*) herwig6_common_block_address_((char*)"HWPROP",6);
  92   fHwunam = (Hwunam_t*) herwig6_common_block_address_((char*)"HWUNAM",6);
  93   fHwupdt = (Hwupdt_t*) herwig6_common_block_address_((char*)"HWUPDT",6);
  94   fHwuwts = (Hwuwts_t*) herwig6_common_block_address_((char*)"HWUWTS",6);
  95   fHwuclu = (Hwuclu_t*) herwig6_common_block_address_((char*)"HWUCLU",6);
  96   fHwdist = (Hwdist_t*) herwig6_common_block_address_((char*)"HWDIST",6);
  97   fHwqdks = (Hwqdks_t*) herwig6_common_block_address_((char*)"HWQDKS",6);
  98   fHwusud = (Hwusud_t*) herwig6_common_block_address_((char*)"HWUSUD",6);
  99   fHwsusy = (Hwsusy_t*) herwig6_common_block_address_((char*)"HWSUSY",6);
 100   fHwrpar = (Hwrpar_t*) herwig6_common_block_address_((char*)"HWRPAR",6);
 101   fHwminb = (Hwminb_t*) herwig6_common_block_address_((char*)"HWMINB",6);
 102   fHwclus = (Hwclus_t*) herwig6_common_block_address_((char*)"HWCLUS",6);
 103   fHwgrav = (Hwgrav_t*) herwig6_common_block_address_((char*)"HWGRAV",6);
 104   fHw6202 = (Hw6202_t*) herwig6_common_block_address_((char*)"HW6202",6);
 105   fHw6203 = (Hw6203_t*) herwig6_common_block_address_((char*)"HW6203",6);
 106   fHw6300 = (Hw6300_t*) herwig6_common_block_address_((char*)"HW6300",6);
 107   fHwpmrs = (Hwpmrs_t*) herwig6_common_block_address_((char*)"HWPMRS",6);
 108   fHwcirc = (Hwcirc_t*) herwig6_common_block_address_((char*)"HWCIRC",6);
 109   fHwdspb = (Hwdspb_t*) herwig6_common_block_address_((char*)"HWDSPB",6);
 110   fHwdsp2 = (Hwdsp2_t*) herwig6_common_block_address_((char*)"HWDSP2",6);
 111   fHwdsp3 = (Hwdsp3_t*) herwig6_common_block_address_((char*)"HWDSP3",6);
 112   fHwdsp4 = (Hwdsp4_t*) herwig6_common_block_address_((char*)"HWDSP4",6);
 113   fHwdspn = (Hwdspn_t*) herwig6_common_block_address_((char*)"HWDSPN",6);
 114   fHwspin = (Hwspin_t*) herwig6_common_block_address_((char*)"HWSPIN",6);
 115   fHwstau = (Hwstau_t*) herwig6_common_block_address_((char*)"HWSTAU",6);
 116   fHwgupr = (Hwgupr_t*) herwig6_common_block_address_((char*)"HWGUPR",6);
 117   fHw6500 = (Hw6500_t*) herwig6_common_block_address_((char*)"HW6500",6);
 118   fHw6504 = (Hw6504_t*) herwig6_common_block_address_((char*)"HW6504",6);
 119   fHw6506 = (Hw6506_t*) herwig6_common_block_address_((char*)"HW6506",6);
 120  }
 121
 122 //------------------------------------------------------------------------------
 123  THerwig6::~THerwig6()
 124  {
 125    // Destructor. The data members of TGenerator are delete by itself
 126  }
 127
 128 //______________________________________________________________________________
 129 void THerwig6::GenerateEvent()
 130 {
 131
 132   //initialize event
 133   hwuine_();
 134   // generate hard subprocess
 135   hwepro_();
 136   // generate parton cascades
 137   hwbgen_();
 138   // do heavy objects decay
 139   hwdhob_();
 140   // do cluster formation
 141   hwcfor_();
 142   // do cluster decays
 143   hwcdec_();
 144   // do unstable particle decays
 145   hwdhad_();
 146   // do heavy flavor hadrons decay
 147   hwdhvy_();
 148   // add soft underlying event
 149   hwmevt_();
 150   // finish event
 151   hwufne_();
 152 }
 153
 154 //______________________________________________________________________________
 155 void THerwig6::OpenFortranFile(int lun, char* name) {
 156   herwig6_open_fortran_file_(&lun, name, strlen(name));
 157 }
 158
 159 //______________________________________________________________________________
 160 void THerwig6::CloseFortranFile(int lun) {
 161   herwig6_close_fortran_file_(&lun);
 162 }
 163
 164 void THerwig6::Initialize(const char *beam, const char *target, double pbeam1, double pbeam2, int iproc)
 165
 166 {
 167   // perform the initialization for Herwig6
 168   // sets correct title.
 169   // after calling this method all parameters are set to their default
 170   // values. If you want to modify any parameter you have to set the new
 171   // value after calling Initialize and before PrepareRun.
 172
 173    char  cbeam[8];
 174    strncpy(cbeam,beam,8);
 175    char  ctarget[8];
 176    strncpy(ctarget,target,8);
 177    printf("\n Initializing Herwig !! \n");
 178    if ( (!strncmp(beam, "E+"    ,2)) &&
 179         (!strncmp(beam, "E-"    ,2)) &&
 180         (!strncmp(beam, "MU+"   ,3)) &&
 181         (!strncmp(beam, "MU-"   ,3)) &&
 182         (!strncmp(beam, "NUE"   ,3)) &&
 183         (!strncmp(beam, "NUEB"  ,4)) &&
 184         (!strncmp(beam, "NUMU"  ,4)) &&
 185         (!strncmp(beam, "NMUB"  ,4)) &&
 186         (!strncmp(beam, "NTAU"  ,4)) &&
 187         (!strncmp(beam, "NTAB"  ,4)) &&
 188         (!strncmp(beam, "GAMA"  ,4)) &&
 189         (!strncmp(beam, "P       ",8)) &&
 190         (!strncmp(beam, "PBAR    ",8)) &&
 191         (!strncmp(beam, "N"     ,1)) &&
 192         (!strncmp(beam, "NBAR"  ,4)) &&
 193         (!strncmp(beam, "PI+"   ,3)) &&
 194         (!strncmp(beam, "PI-"   ,3)) ) {
 195       printf("WARNING! In THerwig6:Initialize():\n");
 196       printf(" specified beam=%s is unrecognized .\n",beam);
 197       printf(" resetting to \"P\" .");
 198       sprintf(cbeam,"P");
 199    }
 200
 201    if ( (!strncmp(target, "E+"    ,2)) &&
 202         (!strncmp(target, "E-"    ,2)) &&
 203         (!strncmp(target, "MU+"   ,3)) &&
 204         (!strncmp(target, "MU-"   ,3)) &&
 205         (!strncmp(target, "NUE"   ,3)) &&
 206         (!strncmp(target, "NUEB"  ,4)) &&
 207         (!strncmp(target, "NUMU"  ,4)) &&
 208         (!strncmp(target, "NMUB"  ,4)) &&
 209         (!strncmp(target, "NTAU"  ,4)) &&
 210         (!strncmp(target, "NTAB"  ,4)) &&
 211         (!strncmp(target, "GAMA"  ,4)) &&
 212         (!strncmp(target, "P       ",8)) &&
 213         (!strncmp(target, "PBAR    ",8)) &&
 214         (!strncmp(target, "N"     ,1)) &&
 215         (!strncmp(target, "NBAR"  ,4)) &&
 216         (!strncmp(target, "PI+"   ,3)) &&
 217         (!strncmp(target, "PI-"   ,3)) ) {
 218       printf("WARNING! In THerwig6:Initialize():\n");
 219       printf(" specified target=%s is unrecognized .\n",target);
 220       printf(" resetting to \"P\" .");
 221       sprintf(ctarget,"P");
 222    }
 223
 224    // initialization:
 225    // type of beams
 226    strncpy(fHwbmch->PART1,beam,8);
 227    strncpy(fHwbmch->PART2,target,8);
 228    // momentum of beams
 229    fHwproc->PBEAM1=pbeam1;
 230    fHwproc->PBEAM2=pbeam2;
 231    // process to generate
 232    fHwproc->IPROC=iproc;
 233    // not used in the class definition
 234    fHwproc->MAXEV=1;
 235
 236    // reset all parameters
 237    hwigin_();
 238
 239    // set correct title
 240    char atitle[132];
 241    double win=pbeam1+pbeam2;
 242    printf("\n %s - %s at %g GeV",beam,target,win);
 243    sprintf(atitle,"%s-%s at %g GeV",cbeam,ctarget,win);
 244    SetTitle(atitle);
 245 }
 246
 247 void THerwig6::PrepareRun()
 248 {
 249   // compute parameter dependent constants
 250   hwuinc_();
 251   // initialize elementary processes
 252   hweini_();
 253 }
 254
 255 //______________________________________________________________________________
 256 TObjArray* THerwig6::ImportParticles(Option_t *option)
 257 {
 258 //
 259 //  Default primary creation method. It reads the /HEPEVT/ common block which
 260 //  has been filled by the GenerateEvent method. If the event generator does
 261 //  not use the HEPEVT common block, This routine has to be overloaded by
 262 //  the subclasses.
 263 //  The default action is to store only the stable particles (ISTHEP = 1)
 264 //  This can be demanded explicitly by setting the option = "Final"
 265 //  If the option = "All", all the particles are stored.
 266 //
 267   fParticles->Clear();
 268   Int_t numpart = fHepevt->NHEP;
 269   TClonesArray &a = *((TClonesArray*)fParticles);
 270   if (!strcmp(option,"") || !strcmp(option,"Final")) {
 271     for (Int_t i = 0; i<=numpart; i++) {
 272       if (fHepevt->ISTHEP[i] == 1) {
 273 //
 274 //  Use the common block values for the TParticle constructor
 275 //
 276         new(a[i]) TParticle(
 277                                    fHepevt->IDHEP[i],
 278                                    fHepevt->ISTHEP[i],
 279                                    fHepevt->JMOHEP[i][0]-1,
 280                                    fHepevt->JMOHEP[i][1]-1,
 281                                    fHepevt->JDAHEP[i][0]-1,
 282                                    fHepevt->JDAHEP[i][1]-1,
 283
 284                                    fHepevt->PHEP[i][0],
 285                                    fHepevt->PHEP[i][1],
 286                                    fHepevt->PHEP[i][2],
 287                                    fHepevt->PHEP[i][3],
 288                                    fHepevt->VHEP[i][0],
 289                                    fHepevt->VHEP[i][1],
 290                                    fHepevt->VHEP[i][2],
 291                                    fHepevt->VHEP[i][3]);
 292         }
 293      }
 294   }
 295   else if (!strcmp(option,"All")) {
 296     for (Int_t i = 0; i<=numpart; i++) {
 297       new(a[i]) TParticle(
 298                                    fHepevt->IDHEP[i],
 299                                    fHepevt->ISTHEP[i],
 300                                    fHepevt->JMOHEP[i][0]-1,
 301                                    fHepevt->JMOHEP[i][1]-1,
 302                                    fHepevt->JDAHEP[i][0]-1,
 303                                    fHepevt->JDAHEP[i][1]-1,
 304
 305                                    fHepevt->PHEP[i][0],
 306                                    fHepevt->PHEP[i][1],
 307                                    fHepevt->PHEP[i][2],
 308                                    fHepevt->PHEP[i][3],
 309                                    fHepevt->VHEP[i][0],
 310                                    fHepevt->VHEP[i][1],
 311                                    fHepevt->VHEP[i][2],
 312                                    fHepevt->VHEP[i][3]);
 313     }
 314   }
 315   return fParticles;
 316 }
 317
 318 //______________________________________________________________________________
 319 Int_t THerwig6::ImportParticles(TClonesArray *particles, Option_t *option)
 320 {
 321 //
 322 //  Default primary creation method. It reads the /HEPEVT/ common block which
 323 //  has been filled by the GenerateEvent method. If the event generator does
 324 //  not use the HEPEVT common block, This routine has to be overloaded by
 325 //  the subclasses.
 326 //  The function loops on the generated particles and store them in
 327 //  the TClonesArray pointed by the argument particles.
 328 //  The default action is to store only the stable particles (ISTHEP = 1)
 329 //  This can be demanded explicitly by setting the option = "Final"
 330 //  If the option = "All", all the particles are stored.
 331 //
 332   if (particles == 0) return 0;
 333   TClonesArray &refParticles = *particles;
 334   refParticles.Clear();
 335   Int_t numpart = fHepevt->NHEP;
 336   if (!strcmp(option,"") || !strcmp(option,"Final")) {
 337     for (Int_t i = 0; i< numpart; i++) {
 338       if (fHepevt->ISTHEP[i] == 1) {
 339 //
 340 //  Use the common block values for the TParticle constructor
 341 //
 342         new(refParticles[i]) TParticle(
 343                                    fHepevt->IDHEP[i],
 344                                    fHepevt->ISTHEP[i],
 345                                    fHepevt->JMOHEP[i][0]-1,
 346                                    fHepevt->JMOHEP[i][1]-1,
 347                                    fHepevt->JDAHEP[i][0]-1,
 348                                    fHepevt->JDAHEP[i][1]-1,
 349
 350                                    fHepevt->PHEP[i][0],
 351                                    fHepevt->PHEP[i][1],
 352                                    fHepevt->PHEP[i][2],
 353                                    fHepevt->PHEP[i][3],
 354                                    fHepevt->VHEP[i][0],
 355                                    fHepevt->VHEP[i][1],
 356                                    fHepevt->VHEP[i][2],
 357                                    fHepevt->VHEP[i][3]);
 358         }
 359      }
 360   }
 361   else if (!strcmp(option,"All")) {
 362     for (Int_t i = 0; i< numpart; i++) {
 363       new(refParticles[i]) TParticle(
 364                                    fHepevt->IDHEP[i],
 365                                    fHepevt->ISTHEP[i],
 366                                    fHepevt->JMOHEP[i][0]-1,
 367                                    fHepevt->JMOHEP[i][1]-1,
 368                                    fHepevt->JDAHEP[i][0]-1,
 369                                    fHepevt->JDAHEP[i][1]-1,
 370
 371                                    fHepevt->PHEP[i][0],
 372                                    fHepevt->PHEP[i][1],
 373                                    fHepevt->PHEP[i][2],
 374                                    fHepevt->PHEP[i][3],
 375                                    fHepevt->VHEP[i][0],
 376                                    fHepevt->VHEP[i][1],
 377                                    fHepevt->VHEP[i][2],
 378                                    fHepevt->VHEP[i][3]);
 379     }
 380   }
 381   return numpart;
 382 }
 383
 384 void THerwig6::Hwigin()
 385 {
 386   hwigin_();
 387 }
 388
 389 void THerwig6::Hwuinc()
 390 {
 391   hwuinc_();
 392 }
 393
 394 void THerwig6::Hwusta(char* name)
 395
 396 {
 397   hwusta_(name,8);
 398 }
 399
 400 void THerwig6::Hweini()
 401
 402 {
 403   hweini_();
 404 }
 405
 406 void THerwig6::Hwuine()
 407
 408 {
 409   hwuine_();
 410 }
 411
 412 void THerwig6::Hwepro()
 413
 414 {
 415   hwepro_();
 416 }
 417
 418 void THerwig6::Hwbgen()
 419
 420 {
 421   hwbgen_();
 422 }
 423
 424 void THerwig6::Hwdhob()
 425
 426 {
 427   hwdhob_();
 428 }
 429
 430 void THerwig6::Hwcfor()
 431
 432 {
 433   hwcfor_();
 434 }
 435
 436 void THerwig6::Hwcdec()
 437
 438 {
 439   hwcdec_();
 440 }
 441
 442 void THerwig6::Hwdhad()
 443
 444 {
 445   hwdhad_();
 446 }
 447
 448 void THerwig6::Hwdhvy()
 449
 450 {
 451   hwdhvy_();
 452 }
 453
 454 void THerwig6::Hwmevt()
 455
 456 {
 457   hwmevt_();
 458 }
 459
 460 void THerwig6::Hwufne()
 461
 462 {
 463   hwufne_();
 464 }
 465
 466 void THerwig6::Hwefin()
 467
 468 {
 469   hwefin_();
 470 }
 471
 472 void THerwig6::Hwiodk(int iopt)
 473
 474 {
 475   hwiodk_(iopt);
 476 }
 477
 478 void THerwig6::SetupTest()
 479 {
 480   // exampe of running herwig and generating one event
 481   // after changing some options
 482   Initialize("P","PBAR",900.,900.,1500);
 483   // here you can set some parameters
 484   SetPTMIN(15.); // Min pt in hadronic jet production
 485   SetYJMIN(-4.); // Min jet rapidity
 486   SetYJMAX(4.);  // Max jet rapidity
 487   // after you set your wished parameters
 488   // herwig can do its work
 489   PrepareRun();
 490   int nEvToGenerate=1;
 491   for (int i=0;i<nEvToGenerate;i++)
 492     {
 493       GenerateEvent();
 494       // do your stuff. For ex:
 495       int nOfPar=GetNumberOfParticles(); // from TGenerator
 496       for (int j=0; j<nOfPar; j++)
 497         {
 498           TParticle* p=GetParticle(j);
 499           // here you do whatever you want with the particle
 500           p->Print();
 501         };
 502     };
 503 }
 504
 505
 506
 507
 508