THerwig/THerwig6.cxx

   1
   2 // definition of c++ Class THerwig6 to be used in ROOT
   3 // this is a c++ interface to the F77 Herwig6 program
   4 // author: j. g. contreras jgcn@moni.mda.cinvestav.mx
   5 // date: december 22, 2000
   6
   7 #include "THerwig6.h"
   8 #include "TClonesArray.h"
   9 #include "TParticle.h"
  10
  11
  12 ClassImp(THerwig6)
  13
  14 extern "C" {
  15   void*  herwig6_common_block_address_(char*, int len);
  16   void   herwig6_open_fortran_file_ (int* lun, char* name, int);
  17   void   herwig6_close_fortran_file_(int* lun);
  18 }
  19
  20
  21 THerwig6::THerwig6() : TGenerator("Herwig6","Herwig6") {
  22
  23 // THerwig6 constructor: creates a TClonesArray in which it will store all
  24 // particles. Note that there may be only one functional THerwig6 object
  25 // at a time, so it's not use to create more than one instance of it.
  26
  27   delete fParticles; // was allocated as TObjArray in TGenerator
  28
  29   fParticles = new TClonesArray("TParticle",50);
  30
  31   // initialize common-blocks
  32   fHepevt = (Hepevt_t*) herwig6_common_block_address_((char*)"HEPEVT",6);
  33   fHwbeam = (Hwbeam_t*) herwig6_common_block_address_((char*)"HWBEAM",6);
  34   fHwbmch = (Hwbmch_t*) herwig6_common_block_address_((char*)"HWBMCH",6);
  35   fHwproc = (Hwproc_t*) herwig6_common_block_address_((char*)"HWPROC",6);
  36   fHwpram = (Hwpram_t*) herwig6_common_block_address_((char*)"HWPRAM",6);
  37   fHwprch = (Hwprch_t*) herwig6_common_block_address_((char*)"HWPRCH",6);
  38   fHwpart = (Hwpart_t*) herwig6_common_block_address_((char*)"HWPART",6);
  39   fHwparp = (Hwparp_t*) herwig6_common_block_address_((char*)"HWPARP",6);
  40   fHwbosc = (Hwbosc_t*) herwig6_common_block_address_((char*)"HWBOSC",6);
  41   fHwparc = (Hwparc_t*) herwig6_common_block_address_((char*)"HWPARC",6);
  42   fHwbrch = (Hwbrch_t*) herwig6_common_block_address_((char*)"HWBRCH",6);
  43   fHwevnt = (Hwevnt_t*) herwig6_common_block_address_((char*)"HWEVNT",6);
  44   fHwhard = (Hwhard_t*) herwig6_common_block_address_((char*)"HWHARD",6);
  45   fHwprop = (Hwprop_t*) herwig6_common_block_address_((char*)"HWPROP",6);
  46   fHwunam = (Hwunam_t*) herwig6_common_block_address_((char*)"HWUNAM",6);
  47   fHwupdt = (Hwupdt_t*) herwig6_common_block_address_((char*)"HWUPDT",6);
  48   fHwuwts = (Hwuwts_t*) herwig6_common_block_address_((char*)"HWUWTS",6);
  49   fHwuclu = (Hwuclu_t*) herwig6_common_block_address_((char*)"HWUCLU",6);
  50   fHwdist = (Hwdist_t*) herwig6_common_block_address_((char*)"HWDIST",6);
  51   fHwqdks = (Hwqdks_t*) herwig6_common_block_address_((char*)"HWQDKS",6);
  52   fHwusud = (Hwusud_t*) herwig6_common_block_address_((char*)"HWUSUD",6);
  53   fHwsusy = (Hwsusy_t*) herwig6_common_block_address_((char*)"HWSUSY",6);
  54   fHwrpar = (Hwrpar_t*) herwig6_common_block_address_((char*)"HWRPAR",6);
  55   fHwminb = (Hwminb_t*) herwig6_common_block_address_((char*)"HWMINB",6);
  56   fHwclus = (Hwclus_t*) herwig6_common_block_address_((char*)"HWCLUS",6);
  57  }
  58
  59 //------------------------------------------------------------------------------
  60  THerwig6::~THerwig6()
  61  {
  62 // Destroys the object, deletes and disposes all TParticles currently on list.
  63
  64     if (fParticles) {
  65       fParticles->Delete();
  66       delete fParticles;
  67       fParticles = 0;
  68    }
  69  }
  70
  71 //______________________________________________________________________________
  72 void THerwig6::GenerateEvent()
  73 {
  74
  75   //initialize event
  76   hwuine_();
  77   // generate hard subprocess
  78   hwepro_();
  79   // generate parton cascades
  80   hwbgen_();
  81   // do heavy objects decay
  82   hwdhob_();
  83   // do cluster formation
  84   hwcfor_();
  85   // do cluster decays
  86   hwcdec_();
  87   // do unstable particle decays
  88   hwdhad_();
  89   // do heavy flavor hadrons decay
  90   hwdhvy_();
  91   // add soft underlying event
  92   hwmevt_();
  93   // finish event
  94   hwufne_();
  95 }
  96
  97 //______________________________________________________________________________
  98 void THerwig6::OpenFortranFile(int lun, char* name) {
  99   herwig6_open_fortran_file_(&lun, name, strlen(name));
 100 }
 101
 102 //______________________________________________________________________________
 103 void THerwig6::CloseFortranFile(int lun) {
 104   herwig6_close_fortran_file_(&lun);
 105 }
 106
 107 void THerwig6::Initialize(const char *beam, const char *target, double pbeam1, double pbeam2, int iproc)
 108
 109 {
 110   // perform the initialization for Herwig6
 111   // sets correct title.
 112   // after calling this method all parameters are set to their default
 113   // values. If you want to modify any parameter you have to set the new
 114   // value after calling Initialize and before PrepareRun.
 115
 116    char  cbeam[8];
 117    strncpy(cbeam,beam,8);
 118    char  ctarget[8];
 119    strncpy(ctarget,target,8);
 120    printf("\n Initializing Herwig !! \n");
 121    if ( (!strncmp(beam, "E+"    ,2)) &&
 122         (!strncmp(beam, "E-"    ,2)) &&
 123         (!strncmp(beam, "MU+"   ,3)) &&
 124         (!strncmp(beam, "MU-"   ,3)) &&
 125         (!strncmp(beam, "NUE"   ,3)) &&
 126         (!strncmp(beam, "NUEB"  ,4)) &&
 127         (!strncmp(beam, "NUMU"  ,4)) &&
 128         (!strncmp(beam, "NMUB"  ,4)) &&
 129         (!strncmp(beam, "NTAU"  ,4)) &&
 130         (!strncmp(beam, "NTAB"  ,4)) &&
 131         (!strncmp(beam, "GAMA"  ,4)) &&
 132         (!strncmp(beam, "P       ",8)) &&
 133         (!strncmp(beam, "PBAR    ",8)) &&
 134         (!strncmp(beam, "N"     ,1)) &&
 135         (!strncmp(beam, "NBAR"  ,4)) &&
 136         (!strncmp(beam, "PI+"   ,3)) &&
 137         (!strncmp(beam, "PI-"   ,3)) ) {
 138       printf("WARNING! In THerwig6:Initialize():\n");
 139       printf(" specified beam=%s is unrecognized .\n",beam);
 140       printf(" resetting to \"P\" .");
 141       sprintf(cbeam,"P");
 142    }
 143
 144    if ( (!strncmp(target, "E+"    ,2)) &&
 145         (!strncmp(target, "E-"    ,2)) &&
 146         (!strncmp(target, "MU+"   ,3)) &&
 147         (!strncmp(target, "MU-"   ,3)) &&
 148         (!strncmp(target, "NUE"   ,3)) &&
 149         (!strncmp(target, "NUEB"  ,4)) &&
 150         (!strncmp(target, "NUMU"  ,4)) &&
 151         (!strncmp(target, "NMUB"  ,4)) &&
 152         (!strncmp(target, "NTAU"  ,4)) &&
 153         (!strncmp(target, "NTAB"  ,4)) &&
 154         (!strncmp(target, "GAMA"  ,4)) &&
 155         (!strncmp(target, "P       ",8)) &&
 156         (!strncmp(target, "PBAR    ",8)) &&
 157         (!strncmp(target, "N"     ,1)) &&
 158         (!strncmp(target, "NBAR"  ,4)) &&
 159         (!strncmp(target, "PI+"   ,3)) &&
 160         (!strncmp(target, "PI-"   ,3)) ) {
 161       printf("WARNING! In THerwig6:Initialize():\n");
 162       printf(" specified target=%s is unrecognized .\n",target);
 163       printf(" resetting to \"P\" .");
 164       sprintf(ctarget,"P");
 165    }
 166
 167    // initialization:
 168    // type of beams
 169    strncpy(fHwbmch->PART1,beam,8);
 170    strncpy(fHwbmch->PART2,target,8);
 171    // momentum of beams
 172    fHwproc->PBEAM1=pbeam1;
 173    fHwproc->PBEAM2=pbeam2;
 174    // process to generate
 175    fHwproc->IPROC=iproc;
 176    // not used in the class definition
 177    fHwproc->MAXEV=1;
 178
 179    // reset all parameters
 180    hwigin_();
 181
 182    // set correct title
 183    char atitle[132];
 184    double win=pbeam1+pbeam2;
 185    printf("\n %s - %s at %g GeV",beam,target,win);
 186    sprintf(atitle,"%s-%s at %g GeV",cbeam,ctarget,win);
 187    SetTitle(atitle);
 188 }
 189
 190 void THerwig6::PrepareRun()
 191 {
 192   // compute parameter dependent constants
 193   hwuinc_();
 194   // initialize elementary processes
 195   hweini_();
 196 }
 197
 198 //______________________________________________________________________________
 199 TObjArray* THerwig6::ImportParticles(Option_t *option)
 200 {
 201 //
 202 //  Default primary creation method. It reads the /HEPEVT/ common block which
 203 //  has been filled by the GenerateEvent method. If the event generator does
 204 //  not use the HEPEVT common block, This routine has to be overloaded by
 205 //  the subclasses.
 206 //  The default action is to store only the stable particles (ISTHEP = 1)
 207 //  This can be demanded explicitly by setting the option = "Final"
 208 //  If the option = "All", all the particles are stored.
 209 //
 210   fParticles->Clear();
 211   Int_t numpart = fHepevt->NHEP;
 212   TClonesArray &a = *((TClonesArray*)fParticles);
 213   if (!strcmp(option,"") || !strcmp(option,"Final")) {
 214     for (Int_t i = 0; i<=numpart; i++) {
 215       if (fHepevt->ISTHEP[i] == 1) {
 216 //
 217 //  Use the common block values for the TParticle constructor
 218 //
 219         new(a[i]) TParticle(
 220                                    fHepevt->IDHEP[i],
 221                                    fHepevt->ISTHEP[i],
 222                                    fHepevt->JMOHEP[i][0]-1,
 223                                    fHepevt->JMOHEP[i][1]-1,
 224                                    fHepevt->JDAHEP[i][0]-1,
 225                                    fHepevt->JDAHEP[i][1]-1,
 226
 227                                    fHepevt->PHEP[i][0],
 228                                    fHepevt->PHEP[i][1],
 229                                    fHepevt->PHEP[i][2],
 230                                    fHepevt->PHEP[i][3],
 231                                    fHepevt->VHEP[i][0],
 232                                    fHepevt->VHEP[i][1],
 233                                    fHepevt->VHEP[i][2],
 234                                    fHepevt->VHEP[i][3]);
 235         }
 236      }
 237   }
 238   else if (!strcmp(option,"All")) {
 239     for (Int_t i = 0; i<=numpart; i++) {
 240       new(a[i]) TParticle(
 241                                    fHepevt->IDHEP[i],
 242                                    fHepevt->ISTHEP[i],
 243                                    fHepevt->JMOHEP[i][0]-1,
 244                                    fHepevt->JMOHEP[i][1]-1,
 245                                    fHepevt->JDAHEP[i][0]-1,
 246                                    fHepevt->JDAHEP[i][1]-1,
 247
 248                                    fHepevt->PHEP[i][0],
 249                                    fHepevt->PHEP[i][1],
 250                                    fHepevt->PHEP[i][2],
 251                                    fHepevt->PHEP[i][3],
 252                                    fHepevt->VHEP[i][0],
 253                                    fHepevt->VHEP[i][1],
 254                                    fHepevt->VHEP[i][2],
 255                                    fHepevt->VHEP[i][3]);
 256     }
 257   }
 258   return fParticles;
 259 }
 260
 261 //______________________________________________________________________________
 262 Int_t THerwig6::ImportParticles(TClonesArray *particles, Option_t *option)
 263 {
 264 //
 265 //  Default primary creation method. It reads the /HEPEVT/ common block which
 266 //  has been filled by the GenerateEvent method. If the event generator does
 267 //  not use the HEPEVT common block, This routine has to be overloaded by
 268 //  the subclasses.
 269 //  The function loops on the generated particles and store them in
 270 //  the TClonesArray pointed by the argument particles.
 271 //  The default action is to store only the stable particles (ISTHEP = 1)
 272 //  This can be demanded explicitly by setting the option = "Final"
 273 //  If the option = "All", all the particles are stored.
 274 //
 275   if (particles == 0) return 0;
 276   TClonesArray &Particles = *particles;
 277   Particles.Clear();
 278   Int_t numpart = fHepevt->NHEP;
 279   if (!strcmp(option,"") || !strcmp(option,"Final")) {
 280     for (Int_t i = 0; i< numpart; i++) {
 281       if (fHepevt->ISTHEP[i] == 1) {
 282 //
 283 //  Use the common block values for the TParticle constructor
 284 //
 285         new(Particles[i]) TParticle(
 286                                    fHepevt->IDHEP[i],
 287                                    fHepevt->ISTHEP[i],
 288                                    fHepevt->JMOHEP[i][0]-1,
 289                                    fHepevt->JMOHEP[i][1]-1,
 290                                    fHepevt->JDAHEP[i][0]-1,
 291                                    fHepevt->JDAHEP[i][1]-1,
 292
 293                                    fHepevt->PHEP[i][0],
 294                                    fHepevt->PHEP[i][1],
 295                                    fHepevt->PHEP[i][2],
 296                                    fHepevt->PHEP[i][3],
 297                                    fHepevt->VHEP[i][0],
 298                                    fHepevt->VHEP[i][1],
 299                                    fHepevt->VHEP[i][2],
 300                                    fHepevt->VHEP[i][3]);
 301         }
 302      }
 303   }
 304   else if (!strcmp(option,"All")) {
 305     for (Int_t i = 0; i< numpart; i++) {
 306       new(Particles[i]) TParticle(
 307                                    fHepevt->IDHEP[i],
 308                                    fHepevt->ISTHEP[i],
 309                                    fHepevt->JMOHEP[i][0]-1,
 310                                    fHepevt->JMOHEP[i][1]-1,
 311                                    fHepevt->JDAHEP[i][0]-1,
 312                                    fHepevt->JDAHEP[i][1]-1,
 313
 314                                    fHepevt->PHEP[i][0],
 315                                    fHepevt->PHEP[i][1],
 316                                    fHepevt->PHEP[i][2],
 317                                    fHepevt->PHEP[i][3],
 318                                    fHepevt->VHEP[i][0],
 319                                    fHepevt->VHEP[i][1],
 320                                    fHepevt->VHEP[i][2],
 321                                    fHepevt->VHEP[i][3]);
 322     }
 323   }
 324   return numpart;
 325 }
 326
 327 void THerwig6::Hwigin()
 328 {
 329   hwigin_();
 330 }
 331
 332 void THerwig6::Hwuinc()
 333 {
 334   hwuinc_();
 335 }
 336
 337 void THerwig6::Hwusta(char* name)
 338
 339 {
 340   hwusta_(name,8);
 341 }
 342
 343 void THerwig6::Hweini()
 344
 345 {
 346   hweini_();
 347 }
 348
 349 void THerwig6::Hwuine()
 350
 351 {
 352   hwuine_();
 353 }
 354
 355 void THerwig6::Hwepro()
 356
 357 {
 358   hwepro_();
 359 }
 360
 361 void THerwig6::Hwbgen()
 362
 363 {
 364   hwbgen_();
 365 }
 366
 367 void THerwig6::Hwdhob()
 368
 369 {
 370   hwdhob_();
 371 }
 372
 373 void THerwig6::Hwcfor()
 374
 375 {
 376   hwcfor_();
 377 }
 378
 379 void THerwig6::Hwcdec()
 380
 381 {
 382   hwcdec_();
 383 }
 384
 385 void THerwig6::Hwdhad()
 386
 387 {
 388   hwdhad_();
 389 }
 390
 391 void THerwig6::Hwdhvy()
 392
 393 {
 394   hwdhvy_();
 395 }
 396
 397 void THerwig6::Hwmevt()
 398
 399 {
 400   hwmevt_();
 401 }
 402
 403 void THerwig6::Hwufne()
 404
 405 {
 406   hwufne_();
 407 }
 408
 409 void THerwig6::Hwefin()
 410
 411 {
 412   hwefin_();
 413 }
 414
 415 void THerwig6::SetupTest()
 416 {
 417   // exampe of running herwig and generating one event
 418   // after changing some options
 419   Initialize("P","PBAR",900.,900.,1500);
 420   // here you can set some parameters
 421   SetPTMIN(15.); // Min pt in hadronic jet production
 422   SetYJMIN(-4.); // Min jet rapidity
 423   SetYJMAX(4.);  // Max jet rapidity
 424   // after you set your wished parameters
 425   // herwig can do its work
 426   PrepareRun();
 427   int n_ev_to_generate=1;
 428   for (int i=0;i<n_ev_to_generate;i++)
 429     {
 430       GenerateEvent();
 431       // do your stuff. For ex:
 432       int n_of_par=GetNumberOfParticles(); // from TGenerator
 433       for (int j=0; j<n_of_par; j++)
 434         {
 435           TParticle* p=GetParticle(j);
 436           // here you do whatever you want with the particle
 437           p->Print();
 438         };
 439     };
 440 }
 441
 442
 443
 444
 445