[u/mrichter/AliRoot.git] / THerwig / THerwig6.cxx

// definition of c++ Class THerwig6 to be used in ROOT
// this is a c++ interface to the F77 Herwig6 program
// author: j. g. contreras jgcn@moni.mda.cinvestav.mx
// date: december 22, 2000
//  Class THerwig6 is an interface to the Herwig program
//
// C-----------------------------------------------------------------------
// C                           H E R W I G
// C
// C            a Monte Carlo event generator for simulating
// C        +---------------------------------------------------+
// C        | Hadron Emission Reactions With Interfering Gluons |
// C        +---------------------------------------------------+
// C I.G. Knowles(*), G. Marchesini(+), M.H. Seymour($) and B.R. Webber(#)
// C-----------------------------------------------------------------------
// C with Minimal Supersymmetric Standard Model Matrix Elements by
// C                  S. Moretti($) and K. Odagiri($)
// C-----------------------------------------------------------------------
// C R parity violating Supersymmetric Decays and Matrix Elements by
// C                          P. Richardson(&)
// C-----------------------------------------------------------------------
// C matrix element corrections to top decay and Drell-Yan type processes
// C                         by G. Corcella(+)
// C-----------------------------------------------------------------------
// C Deep Inelastic Scattering and Heavy Flavour Electroproduction by
// C                  G. Abbiendi(@) and L. Stanco(%)
// C-----------------------------------------------------------------------
// C and Jet Photoproduction in Lepton-Hadron Collisions by J. Chyla(~)
// C-----------------------------------------------------------------------
// C(*)  Department of Physics & Astronomy, University of Edinburgh
// C(+)  Dipartimento di Fisica, Universita di Milano
// C($)  Rutherford Appleton Laboratory
// C(#)  Cavendish Laboratory, Cambridge
// C(&)  Department of Physics, University of Oxford
// C(@)  Dipartimento di Fisica, Universita di Bologna
// C(%)  Dipartimento di Fisica, Universita di Padova
// C(~)  Institute of Physics, Prague
// C-----------------------------------------------------------------------
// C                  Version 6.100 - 16th December 1999
// C-----------------------------------------------------------------------
// C Main reference:
// C    G.Marchesini,  B.R.Webber,  G.Abbiendi,  I.G.Knowles,  M.H.Seymour,
// C    and L.Stanco, Computer Physics Communications 67 (1992) 465.
// C-----------------------------------------------------------------------
// C Please send e-mail about  this program  to one of the  authors at the
// C following Internet addresses:
// C    I.Knowles@ed.ac.uk        Giuseppe.Marchesini@mi.infn.it
// C    M.Seymour@rl.ac.uk        webber@hep.phy.cam.ac.uk
// C-----------------------------------------------------------------------


#include "THerwig6.h"
#include "TClonesArray.h"
#include "TParticle.h"
#include "TObjArray.h"


ClassImp(THerwig6)

extern "C" {
  void*  herwig6_common_block_address_(char*, int len);
  void   herwig6_open_fortran_file_ (int* lun, char* name, int);
  void   herwig6_close_fortran_file_(int* lun);
}


THerwig6::THerwig6() : TGenerator("Herwig6","Herwig6") {

// THerwig6 constructor: creates a TClonesArray in which it will store all
// particles. Note that there may be only one functional THerwig6 object
// at a time, so it's not use to create more than one instance of it.

  delete fParticles; // was allocated as TObjArray in TGenerator

  fParticles = new TClonesArray("TParticle",50);

  // initialize common-blocks
  fHepevt = (Hepevt_t*) herwig6_common_block_address_((char*)"HEPEVT",6);
  fHwbeam = (Hwbeam_t*) herwig6_common_block_address_((char*)"HWBEAM",6);
  fHwbmch = (Hwbmch_t*) herwig6_common_block_address_((char*)"HWBMCH",6);
  fHwproc = (Hwproc_t*) herwig6_common_block_address_((char*)"HWPROC",6);
  fHwpram = (Hwpram_t*) herwig6_common_block_address_((char*)"HWPRAM",6);
  fHwprch = (Hwprch_t*) herwig6_common_block_address_((char*)"HWPRCH",6);
  fHwpart = (Hwpart_t*) herwig6_common_block_address_((char*)"HWPART",6);
  fHwparp = (Hwparp_t*) herwig6_common_block_address_((char*)"HWPARP",6);
  fHwbosc = (Hwbosc_t*) herwig6_common_block_address_((char*)"HWBOSC",6);
  fHwparc = (Hwparc_t*) herwig6_common_block_address_((char*)"HWPARC",6);
  fHwbrch = (Hwbrch_t*) herwig6_common_block_address_((char*)"HWBRCH",6);
  fHwevnt = (Hwevnt_t*) herwig6_common_block_address_((char*)"HWEVNT",6);
  fHwhard = (Hwhard_t*) herwig6_common_block_address_((char*)"HWHARD",6);
  fHwprop = (Hwprop_t*) herwig6_common_block_address_((char*)"HWPROP",6);
  fHwunam = (Hwunam_t*) herwig6_common_block_address_((char*)"HWUNAM",6);
  fHwupdt = (Hwupdt_t*) herwig6_common_block_address_((char*)"HWUPDT",6);
  fHwuwts = (Hwuwts_t*) herwig6_common_block_address_((char*)"HWUWTS",6);
  fHwuclu = (Hwuclu_t*) herwig6_common_block_address_((char*)"HWUCLU",6);
  fHwdist = (Hwdist_t*) herwig6_common_block_address_((char*)"HWDIST",6);
  fHwqdks = (Hwqdks_t*) herwig6_common_block_address_((char*)"HWQDKS",6);
  fHwusud = (Hwusud_t*) herwig6_common_block_address_((char*)"HWUSUD",6);
  fHwsusy = (Hwsusy_t*) herwig6_common_block_address_((char*)"HWSUSY",6);
  fHwrpar = (Hwrpar_t*) herwig6_common_block_address_((char*)"HWRPAR",6);
  fHwminb = (Hwminb_t*) herwig6_common_block_address_((char*)"HWMINB",6);
  fHwclus = (Hwclus_t*) herwig6_common_block_address_((char*)"HWCLUS",6);
  fHwgrav = (Hwgrav_t*) herwig6_common_block_address_((char*)"HWGRAV",6);
  fHw6202 = (Hw6202_t*) herwig6_common_block_address_((char*)"HW6202",6);
  fHw6203 = (Hw6203_t*) herwig6_common_block_address_((char*)"HW6203",6);
  fHw6300 = (Hw6300_t*) herwig6_common_block_address_((char*)"HW6300",6);
  fHwpmrs = (Hwpmrs_t*) herwig6_common_block_address_((char*)"HWPMRS",6);
  fHwcirc = (Hwcirc_t*) herwig6_common_block_address_((char*)"HWCIRC",6);
  fHwdspb = (Hwdspb_t*) herwig6_common_block_address_((char*)"HWDSPB",6);
  fHwdsp2 = (Hwdsp2_t*) herwig6_common_block_address_((char*)"HWDSP2",6);
  fHwdsp3 = (Hwdsp3_t*) herwig6_common_block_address_((char*)"HWDSP3",6);
  fHwdsp4 = (Hwdsp4_t*) herwig6_common_block_address_((char*)"HWDSP4",6);
  fHwdspn = (Hwdspn_t*) herwig6_common_block_address_((char*)"HWDSPN",6);
  fHwspin = (Hwspin_t*) herwig6_common_block_address_((char*)"HWSPIN",6);
  fHwstau = (Hwstau_t*) herwig6_common_block_address_((char*)"HWSTAU",6);
  fHwgupr = (Hwgupr_t*) herwig6_common_block_address_((char*)"HWGUPR",6);
  fHw6500 = (Hw6500_t*) herwig6_common_block_address_((char*)"HW6500",6);
  fHw6504 = (Hw6504_t*) herwig6_common_block_address_((char*)"HW6504",6);
  fHw6506 = (Hw6506_t*) herwig6_common_block_address_((char*)"HW6506",6);
 }

//------------------------------------------------------------------------------
 THerwig6::~THerwig6()
 {
   // Destructor. The data members of TGenerator are delete by itself
 }

//______________________________________________________________________________
void THerwig6::GenerateEvent() 
{
  
  //initialize event
  hwuine_();
  // generate hard subprocess
  hwepro_();
  // generate parton cascades
  hwbgen_();
  // do heavy objects decay
  hwdhob_();
  // do cluster formation
  hwcfor_();
  // do cluster decays
  hwcdec_();
  // do unstable particle decays
  hwdhad_();
  // do heavy flavor hadrons decay
  hwdhvy_();
  // add soft underlying event
  hwmevt_();
  // finish event
  hwufne_();
}

//______________________________________________________________________________
void THerwig6::OpenFortranFile(int lun, char* name) {
  herwig6_open_fortran_file_(&lun, name, strlen(name));
}

//______________________________________________________________________________
void THerwig6::CloseFortranFile(int lun) {
  herwig6_close_fortran_file_(&lun);
}

void THerwig6::Initialize(const char *beam, const char *target, double pbeam1, double pbeam2, int iproc)

{
  // perform the initialization for Herwig6
  // sets correct title. 
  // after calling this method all parameters are set to their default
  // values. If you want to modify any parameter you have to set the new
  // value after calling Initialize and before PrepareRun.

   char  cbeam[8];
   strncpy(cbeam,beam,8);
   char  ctarget[8];
   strncpy(ctarget,target,8);
   printf("\n Initializing Herwig !! \n");
   if ( (!strncmp(beam, "E+"    ,2)) &&
        (!strncmp(beam, "E-"    ,2)) &&
        (!strncmp(beam, "MU+"   ,3)) &&
        (!strncmp(beam, "MU-"   ,3)) &&
        (!strncmp(beam, "NUE"   ,3)) &&
        (!strncmp(beam, "NUEB"  ,4)) &&
        (!strncmp(beam, "NUMU"  ,4)) &&
        (!strncmp(beam, "NMUB"  ,4)) &&
        (!strncmp(beam, "NTAU"  ,4)) &&
        (!strncmp(beam, "NTAB"  ,4)) &&
        (!strncmp(beam, "GAMA"  ,4)) &&
        (!strncmp(beam, "P       ",8)) &&
        (!strncmp(beam, "PBAR    ",8)) &&
        (!strncmp(beam, "N"     ,1)) &&
        (!strncmp(beam, "NBAR"  ,4)) &&
        (!strncmp(beam, "PI+"   ,3)) &&
        (!strncmp(beam, "PI-"   ,3)) ) {
      printf("WARNING! In THerwig6:Initialize():\n");
      printf(" specified beam=%s is unrecognized .\n",beam);
      printf(" resetting to \"P\" .");
      sprintf(cbeam,"P");
   }

   if ( (!strncmp(target, "E+"    ,2)) &&
        (!strncmp(target, "E-"    ,2)) &&
        (!strncmp(target, "MU+"   ,3)) &&
        (!strncmp(target, "MU-"   ,3)) &&
        (!strncmp(target, "NUE"   ,3)) &&
        (!strncmp(target, "NUEB"  ,4)) &&
        (!strncmp(target, "NUMU"  ,4)) &&
        (!strncmp(target, "NMUB"  ,4)) &&
        (!strncmp(target, "NTAU"  ,4)) &&
        (!strncmp(target, "NTAB"  ,4)) &&
        (!strncmp(target, "GAMA"  ,4)) &&
        (!strncmp(target, "P       ",8)) &&
        (!strncmp(target, "PBAR    ",8)) &&
        (!strncmp(target, "N"     ,1)) &&
        (!strncmp(target, "NBAR"  ,4)) &&
        (!strncmp(target, "PI+"   ,3)) &&
        (!strncmp(target, "PI-"   ,3)) ) {
      printf("WARNING! In THerwig6:Initialize():\n");
      printf(" specified target=%s is unrecognized .\n",target);
      printf(" resetting to \"P\" .");
      sprintf(ctarget,"P");
   }

   // initialization:
   // type of beams
   strncpy(fHwbmch->PART1,beam,8);
   strncpy(fHwbmch->PART2,target,8);
   // momentum of beams
   fHwproc->PBEAM1=pbeam1;
   fHwproc->PBEAM2=pbeam2;
   // process to generate
   fHwproc->IPROC=iproc;
   // not used in the class definition 
   fHwproc->MAXEV=1;

   // reset all parameters
   hwigin_();
   
   // set correct title
   char atitle[132];
   double win=pbeam1+pbeam2;
   printf("\n %s - %s at %g GeV",beam,target,win);
   sprintf(atitle,"%s-%s at %g GeV",cbeam,ctarget,win);
   SetTitle(atitle); 
}

void THerwig6::PrepareRun()
{
  // compute parameter dependent constants
  hwuinc_();
  // initialize elementary processes
  hweini_();
}

//______________________________________________________________________________
TObjArray* THerwig6::ImportParticles(Option_t *option)
{
//
//  Default primary creation method. It reads the /HEPEVT/ common block which
//  has been filled by the GenerateEvent method. If the event generator does
//  not use the HEPEVT common block, This routine has to be overloaded by
//  the subclasses.
//  The default action is to store only the stable particles (ISTHEP = 1)
//  This can be demanded explicitly by setting the option = "Final"
//  If the option = "All", all the particles are stored.
//
  fParticles->Clear();
  Int_t numpart = fHepevt->NHEP;
  TClonesArray &a = *((TClonesArray*)fParticles);
  if (!strcmp(option,"") || !strcmp(option,"Final")) {
    for (Int_t i = 0; i<=numpart; i++) {
      if (fHepevt->ISTHEP[i] == 1) {
//
//  Use the common block values for the TParticle constructor
//
        new(a[i]) TParticle(
                                   fHepevt->IDHEP[i],
                                   fHepevt->ISTHEP[i],
                                   fHepevt->JMOHEP[i][0]-1,
                                   fHepevt->JMOHEP[i][1]-1,
                                   fHepevt->JDAHEP[i][0]-1,
                                   fHepevt->JDAHEP[i][1]-1,

                                   fHepevt->PHEP[i][0],
                                   fHepevt->PHEP[i][1],
                                   fHepevt->PHEP[i][2],
                                   fHepevt->PHEP[i][3],
                                   fHepevt->VHEP[i][0],
                                   fHepevt->VHEP[i][1],
                                   fHepevt->VHEP[i][2],
                                   fHepevt->VHEP[i][3]);
        }
     }
  }
  else if (!strcmp(option,"All")) {
    for (Int_t i = 0; i<=numpart; i++) {
      new(a[i]) TParticle(
                                   fHepevt->IDHEP[i],
                                   fHepevt->ISTHEP[i],
                                   fHepevt->JMOHEP[i][0]-1,
                                   fHepevt->JMOHEP[i][1]-1,
                                   fHepevt->JDAHEP[i][0]-1,
                                   fHepevt->JDAHEP[i][1]-1,

                                   fHepevt->PHEP[i][0],
                                   fHepevt->PHEP[i][1],
                                   fHepevt->PHEP[i][2],
                                   fHepevt->PHEP[i][3],
                                   fHepevt->VHEP[i][0],
                                   fHepevt->VHEP[i][1],
                                   fHepevt->VHEP[i][2],
                                   fHepevt->VHEP[i][3]);
    }
  }
  return fParticles;
}

//______________________________________________________________________________
Int_t THerwig6::ImportParticles(TClonesArray *particles, Option_t *option)
{
//
//  Default primary creation method. It reads the /HEPEVT/ common block which
//  has been filled by the GenerateEvent method. If the event generator does
//  not use the HEPEVT common block, This routine has to be overloaded by
//  the subclasses.
//  The function loops on the generated particles and store them in
//  the TClonesArray pointed by the argument particles.
//  The default action is to store only the stable particles (ISTHEP = 1)
//  This can be demanded explicitly by setting the option = "Final"
//  If the option = "All", all the particles are stored.
//
  if (particles == 0) return 0;
  TClonesArray &refParticles = *particles;
  refParticles.Clear();
  Int_t numpart = fHepevt->NHEP;
  if (!strcmp(option,"") || !strcmp(option,"Final")) {
    for (Int_t i = 0; i< numpart; i++) {
      if (fHepevt->ISTHEP[i] == 1) {
//
//  Use the common block values for the TParticle constructor
//
        new(refParticles[i]) TParticle(
                                   fHepevt->IDHEP[i],
                                   fHepevt->ISTHEP[i],
                                   fHepevt->JMOHEP[i][0]-1,
                                   fHepevt->JMOHEP[i][1]-1,
                                   fHepevt->JDAHEP[i][0]-1,
                                   fHepevt->JDAHEP[i][1]-1,

                                   fHepevt->PHEP[i][0],
                                   fHepevt->PHEP[i][1],
                                   fHepevt->PHEP[i][2],
                                   fHepevt->PHEP[i][3],
                                   fHepevt->VHEP[i][0],
                                   fHepevt->VHEP[i][1],
                                   fHepevt->VHEP[i][2],
                                   fHepevt->VHEP[i][3]);
        }
     }
  }
  else if (!strcmp(option,"All")) {
    for (Int_t i = 0; i< numpart; i++) {
      new(refParticles[i]) TParticle(
                                   fHepevt->IDHEP[i],
                                   fHepevt->ISTHEP[i],
                                   fHepevt->JMOHEP[i][0]-1,
                                   fHepevt->JMOHEP[i][1]-1,
                                   fHepevt->JDAHEP[i][0]-1,
                                   fHepevt->JDAHEP[i][1]-1,

                                   fHepevt->PHEP[i][0],
                                   fHepevt->PHEP[i][1],
                                   fHepevt->PHEP[i][2],
                                   fHepevt->PHEP[i][3],
                                   fHepevt->VHEP[i][0],
                                   fHepevt->VHEP[i][1],
                                   fHepevt->VHEP[i][2],
                                   fHepevt->VHEP[i][3]);
    }
  }
  return numpart;
}

void THerwig6::Hwigin()
{
  hwigin_();
}

void THerwig6::Hwuinc()
{
  hwuinc_();
}

void THerwig6::Hwusta(char* name)

{
  hwusta_(name,8);
}

void THerwig6::Hweini()

{
  hweini_();
}

void THerwig6::Hwuine()

{
  hwuine_();
}

void THerwig6::Hwepro()

{
  hwepro_();
}

void THerwig6::Hwbgen()

{
  hwbgen_();
}

void THerwig6::Hwdhob()

{
  hwdhob_();
}

void THerwig6::Hwcfor()

{
  hwcfor_();
}

void THerwig6::Hwcdec()

{
  hwcdec_();
}

void THerwig6::Hwdhad()

{
  hwdhad_();
}

void THerwig6::Hwdhvy()

{
  hwdhvy_();
}

void THerwig6::Hwmevt()

{
  hwmevt_();
}

void THerwig6::Hwufne()

{
  hwufne_();
}

void THerwig6::Hwefin()

{
  hwefin_();
}

void THerwig6::Hwiodk(int iopt)

{
  hwiodk_(iopt);
}

void THerwig6::SetupTest()
{
  // exampe of running herwig and generating one event
  // after changing some options
  Initialize("P","PBAR",900.,900.,1500);
  // here you can set some parameters
  SetPTMIN(15.); // Min pt in hadronic jet production
  SetYJMIN(-4.); // Min jet rapidity
  SetYJMAX(4.);  // Max jet rapidity
  // after you set your wished parameters
  // herwig can do its work
  PrepareRun();
  int nEvToGenerate=1;
  for (int i=0;i<nEvToGenerate;i++)
    {
      GenerateEvent();
      // do your stuff. For ex:
      int nOfPar=GetNumberOfParticles(); // from TGenerator
      for (int j=0; j<nOfPar; j++)
	{
	  TParticle* p=GetParticle(j);
	  // here you do whatever you want with the particle
	  p->Print();
	};
    };
}
Commit	Line	Data
2d8d4543	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
884ccf6a	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
2d8d4543	51
2d8d4543	52	#include "THerwig6.h"
	53	#include "TClonesArray.h"
	54	#include "TParticle.h"
884ccf6a	55	#include "TObjArray.h"
2d8d4543	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);
e2054d85	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);
2d8d4543	120	}
	121
	122	//------------------------------------------------------------------------------
	123	THerwig6::~THerwig6()
	124	{
bbc19b09	125	// Destructor. The data members of TGenerator are delete by itself
2d8d4543	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;
884ccf6a	333	TClonesArray &refParticles = *particles;
884ccf6a	334	refParticles.Clear();
2d8d4543	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	//
884ccf6a	342	new(refParticles[i]) TParticle(
2d8d4543	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++) {
884ccf6a	363	new(refParticles[i]) TParticle(
2d8d4543	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
e2054d85	472	void THerwig6::Hwiodk(int iopt)
	473
	474	{
	475	hwiodk_(iopt);
	476	}
	477
2d8d4543	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();
884ccf6a	490	int nEvToGenerate=1;
884ccf6a	491	for (int i=0;i<nEvToGenerate;i++)
2d8d4543	492	{
	493	GenerateEvent();
	494	// do your stuff. For ex:
884ccf6a	495	int nOfPar=GetNumberOfParticles(); // from TGenerator
884ccf6a	496	for (int j=0; j<nOfPar; j++)
2d8d4543	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