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