Update to HERWIG 6.5-10 and JIMMY. (Y. Delgado)
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 2 Oct 2006 11:05:14 +0000 (11:05 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 2 Oct 2006 11:05:14 +0000 (11:05 +0000)
THerwig/AliGenHerwig.cxx
THerwig/AliGenHerwig.h
THerwig/HCommon.h
THerwig/Herwig6Calls.h
THerwig/THerwig6.cxx
THerwig/THerwig6.h

index 2753e0d..d7d5590 100644 (file)
@@ -37,8 +37,8 @@ ClassImp(AliGenHerwig)
 
   AliGenHerwig::AliGenHerwig() :
     AliGenMC(),
-    fAutPDF("GRV"),
-    fModPDF(5),
+    fAutPDF("LHAPDF"),
+    fModPDF(19070),
     fStrucFunc(kCTEQ5L),
     fKeep(0),
     fDecaysOff(1),
@@ -67,8 +67,8 @@ ClassImp(AliGenHerwig)
 
 AliGenHerwig::AliGenHerwig(Int_t npart)
     :AliGenMC(npart),
-    fAutPDF("GRV"),
-    fModPDF(5),
+    fAutPDF("LHAPDF"),
+    fModPDF(19070),
     fStrucFunc(kCTEQ5L),
     fKeep(0),
     fDecaysOff(1),
@@ -94,14 +94,14 @@ AliGenHerwig::AliGenHerwig(Int_t npart)
 {
     SetTarget();
     SetProjectile();
-    // Set random number generator   
+    // Set random number generator
     AliHerwigRndm::SetHerwigRandom(GetRandom());
 }
 
 AliGenHerwig::AliGenHerwig(const AliGenHerwig & Herwig)
     :AliGenMC(Herwig),
-    fAutPDF("GRV"),
-    fModPDF(5),
+    fAutPDF("LHAPDF"),
+    fModPDF(19070),
     fStrucFunc(kCTEQ5L),
     fKeep(0),
     fDecaysOff(1),
@@ -158,7 +158,7 @@ void AliGenHerwig::Init()
   fHerwig->SetMAXPR(fMaxPr);
   fHerwig->SetMAXER(fMaxErrors);
   fHerwig->SetENSOF(fEnSoft);
-  
+
   fHerwig->SetEV1PR(fEv1Pr);
   fHerwig->SetEV2PR(fEv2Pr);
 
@@ -167,13 +167,13 @@ void AliGenHerwig::Init()
 //       RMASS(2)=0.32
 //       RMASS(3)=0.5
 //       RMASS(4)=1.55
-//       RMASS(5)=4.95
+//       RMASS(5)=4.75
 //       RMASS(6)=174.3
 //       RMASS(13)=0.75
 
   fHerwig->SetRMASS(4,1.2);
   fHerwig->SetRMASS(5,4.75);
-  
+
   if ( fProcess < 0 ) strncpy(VVJIN.QQIN,fFileName.Data(),50);
 
   fHerwig->Hwusta("PI0     ");
@@ -182,55 +182,100 @@ void AliGenHerwig::Init()
   fHerwig->PrepareRun();
 }
 
+void AliGenHerwig::InitJimmy()
+{
+// Initialisation
+  fTarget.Resize(8);
+  fProjectile.Resize(8);
+  SetMC(new THerwig6());
+  fHerwig=(THerwig6*) fMCEvGen;
+  // initialize common blocks
+  fHerwig->InitializeJimmy(fProjectile, fTarget, fMomentum1, fMomentum2, fProcess);
+  // reset parameters according to user needs
+  InitPDF();
+  fHerwig->SetPTMIN(fPtHardMin);
+  fHerwig->SetPTRMS(fPtRMS);
+  fHerwig->SetMAXPR(fMaxPr);
+  fHerwig->SetMAXER(fMaxErrors);
+  fHerwig->SetENSOF(fEnSoft);
+
+  fHerwig->SetEV1PR(fEv1Pr);
+  fHerwig->SetEV2PR(fEv2Pr);
+
+// C---D,U,S,C,B,T QUARK AND GLUON MASSES (IN THAT ORDER)
+//       RMASS(1)=0.32
+//       RMASS(2)=0.32
+//       RMASS(3)=0.5
+//       RMASS(4)=1.55
+//       RMASS(5)=4.75
+//       RMASS(6)=174.3
+//       RMASS(13)=0.75
+
+  fHerwig->SetRMASS(4,1.2);
+  fHerwig->SetRMASS(5,4.75);
+
+  if ( fProcess < 0 ) strncpy(VVJIN.QQIN,fFileName.Data(),50);
+
+  fHerwig->Hwusta("PI0     ");
+
+  // compute parameter dependent constants
+  fHerwig->PrepareRunJimmy();
+}
+
 void AliGenHerwig::InitPDF()
 {
   switch(fStrucFunc)
     {
-//    case kGRVLO:
-//      fModPDF=5;
-//      fAutPDF="GRV";
-//      break;
-//    case kGRVHO:
-//      fModPDF=6;
-//      fAutPDF="GRV";
-//      break;
+// ONLY USES LHAPDF STRUCTURE FUNCTIONS
     case kGRVLO98:
-      fModPDF=12;
-      fAutPDF="GRV";
+      fModPDF=80060;
+      fAutPDF="HWLHAPDF";
       break;
-//    case kMRSDminus:
-//      fModPDF=31;
-//      fAutPDF="MRS";
-//      break;
-//    case kMRSD0:
-//      fModPDF=30;
-//      fAutPDF="MRS";
-//      break;
-//    case kMRSG:
-//      fModPDF=41;
-//      fAutPDF="MRS";
-//      break;
-//    case kMRSTcgLO:
-//      fModPDF=72;
-//      fAutPDF="MRS";
-//      break;
-    case kCTEQ4M:
-      fModPDF=34;
-      fAutPDF="CTEQ";
+    case kCTEQ6:
+      fModPDF=10040;
+      fAutPDF="HWLHAPDF";
+      break;
+    case kCTEQ61:
+      fModPDF=10100;
+      fAutPDF="HWLHAPDF";
+      break;
+    case kCTEQ6m:
+      fModPDF=10050;
+      fAutPDF="HWLHAPDF";
+      break;
+    case kCTEQ6l:
+      fModPDF=10041;
+      fAutPDF="HWLHAPDF";
+      break;
+    case kCTEQ6ll:
+      fModPDF=10042;
+      fAutPDF="HWLHAPDF";
+      break;
+    case kCTEQ5M:
+      fModPDF=19050;
+      fAutPDF="HWLHAPDF";
       break;
     case kCTEQ5L:
-      fModPDF=46;
-      fAutPDF="CTEQ";
+      fModPDF=19070;
+      fAutPDF="HWLHAPDF";
+      break;
+    case kCTEQ4M:
+      fModPDF=19150;
+      fAutPDF="HWLHAPDF";
+      break;
+    case kCTEQ4L:
+      fModPDF=19170;
+      fAutPDF="HWLHAPDF";
+      break;
+    case kMRST2004nlo:
+      fModPDF=20400;
+      fAutPDF="HWLHAPDF";
       break;
-//    case kCTEQ5M:
-//      fModPDF=48;
-//      fAutPDF="CTEQ";
-//      break;
     default:
       cerr << "This structure function is not inplemented " << fStrucFunc << endl;
       break;
     }
-  fAutPDF.Resize(20);      
+  fAutPDF.Resize(20);
   fHerwig->SetMODPDF(1,fModPDF);
   fHerwig->SetMODPDF(2,fModPDF);
   fHerwig->SetAUTPDF(1,fAutPDF);
@@ -245,7 +290,7 @@ void AliGenHerwig::Generate()
   Float_t origin[3]=   {0,0,0};
   Float_t origin0[3]=  {0,0,0};
   Float_t p[4], random[6];
-  
+
   static TClonesArray *particles;
   //  converts from mm/c to s
   const Float_t kconv=0.001/2.999792458e8;
@@ -254,9 +299,9 @@ void AliGenHerwig::Generate()
   Int_t jev=0;
   Int_t j, kf, ks, imo;
   kf=0;
-  
+
   if(!particles) particles=new TClonesArray("TParticle",10000);
-  
+
   fTrials=0;
   for (j=0;j<3;j++) origin0[j]=fOrigin[j];
   if(fVertexSmear==kPerEvent) {
@@ -266,7 +311,7 @@ void AliGenHerwig::Generate()
        TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
     }
   }
-  
+
   while(1)
     {
        fHerwig->GenerateEvent();
@@ -274,19 +319,19 @@ void AliGenHerwig::Generate()
        fHerwig->ImportParticles(particles,"All");
        Int_t np = particles->GetEntriesFast()-1;
        if (np == 0 ) continue;
-       
+
        Int_t nc=0;
-       
+
        Int_t * newPos = new Int_t[np];
        for (Int_t i = 0; i<np; i++) *(newPos+i)=-1;
-       
+
        for (Int_t i = 0; i<np; i++) {
            TParticle *  iparticle       = (TParticle *) particles->At(i);
            imo = iparticle->GetFirstMother();
            kf        = iparticle->GetPdgCode();
            ks        = iparticle->GetStatusCode();
-           if (ks != 3 && 
-               KinematicSelection(iparticle,0)) 
+           if (ks != 3 &&
+               KinematicSelection(iparticle,0))
            {
                nc++;
                p[0]=iparticle->Px();
@@ -301,7 +346,7 @@ void AliGenHerwig::Generate()
                Int_t   trackIt = (ks == 1) && fTrackIt;
                PushTrack(trackIt, iparent, kf,
                          p[0], p[1], p[2], p[3],
-                         origin[0], origin[1], origin[2], 
+                         origin[0], origin[1], origin[2],
                          tof,
                          polar[0], polar[1], polar[2],
                          kPPrimary, nt, fHerwig->GetEVWGT(), ks);
@@ -336,7 +381,7 @@ void AliGenHerwig::AdjustWeights()
         part->SetWeight(part->GetWeight()*fKineBias);
     }
 }
+
 
 void AliGenHerwig::KeepFullEvent()
 {
@@ -356,9 +401,9 @@ Bool_t AliGenHerwig::DaughtersSelection(TParticle* iparticle, TClonesArray* part
     Bool_t selected=kFALSE;
     if (hasDaughters) {
        imin=iparticle->GetFirstDaughter();
-       imax=iparticle->GetLastDaughter();       
+       imax=iparticle->GetLastDaughter();
        for (i=imin; i<= imax; i++){
-           TParticle *  jparticle       = (TParticle *) particles->At(i);      
+           TParticle *  jparticle       = (TParticle *) particles->At(i);
            Int_t ip=jparticle->GetPdgCode();
            if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
                selected=kTRUE; break;
@@ -380,7 +425,7 @@ Bool_t AliGenHerwig::SelectFlavor(Int_t pid)
 // 4: charm and beauty
 // 5: beauty
     if (fFlavor == 0) return kTRUE;
-    
+
     Int_t ifl=TMath::Abs(pid/100);
     if (ifl > 10) ifl/=10;
     return (fFlavor == ifl);
@@ -391,9 +436,9 @@ Bool_t AliGenHerwig::Stable(TParticle*  particle)
 // Return true for a stable particle
 //
     Int_t kf = TMath::Abs(particle->GetPdgCode());
-    
+
     if ( (particle->GetFirstDaughter() < 0 ) || (kf == 1000*fFlavor+122))
-        
+
     {
        return kTRUE;
     } else {
@@ -414,4 +459,10 @@ AliGenHerwig& AliGenHerwig::operator=(const  AliGenHerwig& rhs)
     return (*this);
 }
 
+void AliGenHerwig::FinishRunJimmy()
+{
+  fHerwig->Hwefin();
+  fHerwig->Jmefin();
+
+}
 
index 3d67815..e15a9e7 100644 (file)
@@ -34,25 +34,27 @@ class AliGenHerwig : public AliGenMC
     virtual ~AliGenHerwig();
     virtual void    Generate();
     virtual void    Init();
+    virtual void    InitJimmy();
     // set centre of mass energy
     virtual void    SetBeamMomenta(Float_t p1=7000., Float_t p2=7000.)
        {fMomentum1 = p1; fMomentum2 = p2;}
-    virtual void    SetProcess(Int_t proc)            {fProcess = proc;}    
+    virtual void    SetProcess(Int_t proc)            {fProcess = proc;}
     virtual void    KeepFullEvent();
     virtual void    SetDecaysOff(Int_t flag=1)        {fDecaysOff = flag;}
     virtual void    SetTrigger(Int_t flag=kNoTrigger) {fTrigger   = flag;}
-    virtual void    SetFlavor(Int_t flag=0)           {fFlavor    = flag;}    
-    virtual void    SetSelectAll(Int_t flag=0)        {fSelectAll = flag;}    
+    virtual void    SetFlavor(Int_t flag=0)           {fFlavor    = flag;}
+    virtual void    SetSelectAll(Int_t flag=0)        {fSelectAll = flag;}
     AliGenHerwig &  operator=(const AliGenHerwig & rhs);
-    virtual void    SetStrucFunc(StrucFunc_t func = kCTEQ5L) 
+    virtual void    SetStrucFunc(StrucFunc_t func = kCTEQ5L)
       {fStrucFunc = func;}
     virtual void    SetPtHardMin(Double_t pt) {fPtHardMin=pt;}
     virtual void    SetPtRMS(Double_t pt) {fPtRMS=pt;}
     virtual void    SetMaxPr(Int_t i) {fMaxPr=i;}
     virtual void    SetMaxErrors(Int_t i) {fMaxErrors=i;}
     virtual void    FinishRun();
+    virtual void    FinishRunJimmy();
     virtual void    SetEnSoft(Double_t e) {fEnSoft=e;}
-    
+
     virtual void    SetHardProcessFile(TString filename) {fFileName=TString(filename);};
     virtual void    SetEventListRange(Int_t eventFirst=-1, Int_t eventLast=-1);
 
@@ -76,7 +78,7 @@ class AliGenHerwig : public AliGenMC
     Float_t     fXsection;       // Cross-section
     THerwig6    *fHerwig;        // Herwig
     Int_t       fProcess;        // Process number
-    Double_t    fPtHardMin;      // lower pT-hard cut 
+    Double_t    fPtHardMin;      // lower pT-hard cut
     Double_t    fPtRMS;          // intrinsic pt of incoming hadrons
     Int_t       fMaxPr;          // maximum number of events to print out
     Int_t       fMaxErrors;      // maximum number of errors allowed
@@ -84,7 +86,7 @@ class AliGenHerwig : public AliGenMC
     Int_t       fEv1Pr;          // first event to be printed
     Int_t       fEv2Pr;          // last event to be printed
     TString     fFileName;       //!Name of file to read from hard scattering
-      
+
  private:
     // check if particle is selected as parent particle
     Bool_t ParentSelected(Int_t ip);
@@ -96,7 +98,7 @@ class AliGenHerwig : public AliGenMC
     Bool_t DaughtersSelection(TParticle* iparticle, TClonesArray* particles);
     // check if stable
     Bool_t Stable(TParticle*  particle);
-    
+
     void InitPDF();
 
     ClassDef(AliGenHerwig,1) // AliGenerator interface to Herwig
index 63d2c3b..f92fe56 100644 (file)
@@ -640,6 +640,15 @@ extern "C" {
   void  hwefin_();
 }
 
+// subroutines to be call by JIMMY
+
+extern "C" {
+  void  jminit_();
+  void  jimmin_();
+  void  jmefin_();
+}
+
+
 
 
 
index e3b4d8f..c1b9890 100644 (file)
@@ -643,6 +643,14 @@ extern "C" {
   void  hwefin_();
 }
 
+// subroutines to be call by JIMMY
+
+extern "C" {
+  void  jminit_();
+  void  jimmin_();
+  void  jmefin_();
+}
+
 
 #endif
 
index 7a2dcc4..47ac6a3 100644 (file)
@@ -119,7 +119,7 @@ THerwig6::THerwig6() : TGenerator("Herwig6","Herwig6"),
   // initialize common-blocks
  }
 
-THerwig6::THerwig6(const THerwig6 & source): TGenerator(source), 
+THerwig6::THerwig6(const THerwig6 & source): TGenerator(source),
   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)),
@@ -172,9 +172,9 @@ THerwig6::THerwig6(const THerwig6 & source): TGenerator(source),
  }
 
 //______________________________________________________________________________
-void THerwig6::GenerateEvent() 
+void THerwig6::GenerateEvent()
 {
-  
+
   //initialize event
   hwuine_();
   // generate hard subprocess
@@ -196,7 +196,6 @@ void THerwig6::GenerateEvent()
   // finish event
   hwufne_();
 }
-
 //______________________________________________________________________________
 void THerwig6::OpenFortranFile(int lun, char* name) {
   herwig6_open_fortran_file_(&lun, name, strlen(name));
@@ -211,7 +210,7 @@ void THerwig6::Initialize(const char *beam, const char *target, double pbeam1, d
 
 {
   // perform the initialization for Herwig6
-  // sets correct title. 
+  // 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.
@@ -276,18 +275,103 @@ void THerwig6::Initialize(const char *beam, const char *target, double pbeam1, d
    fHwproc->PBEAM2=pbeam2;
    // process to generate
    fHwproc->IPROC=iproc;
-   // not used in the class definition 
+   // 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); 
+   SetTitle(atitle);
+}
+
+void THerwig6::InitializeJimmy(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_();
+   // JIMMY initialization
+   jimmin_();
+
+   // 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()
@@ -298,6 +382,15 @@ void THerwig6::PrepareRun()
   hweini_();
 }
 
+void THerwig6::PrepareRunJimmy()
+{
+  // compute parameter dependent constants
+  hwuinc_();
+  // initialize elementary processes
+  hweini_();
+  // more initializations for JIMMY
+  jminit_();
+}
 //______________________________________________________________________________
 TObjArray* THerwig6::ImportParticles(Option_t *option)
 {
@@ -548,6 +641,24 @@ void THerwig6::SetupTest()
     };
 }
 
+// Jimmy subroutines
+
+void THerwig6::Jminit()
+{
+  jminit_();
+}
+
+void THerwig6::Jimmin()
+{
+  jimmin_();
+}
+
+void THerwig6::Jmefin()
+{
+  jmefin_();
+}
+
+
 
 
 
index 11528d9..f55dcd1 100644 (file)
@@ -60,10 +60,11 @@ C-----------------------------------------------------------------------
 
 typedef enum
 {
-   kHwCharm         =  1704, 
+   kHwCharm         =  1704,
    kHwBeauty        =  1705,
    kHwCharmMCATNLO  = -1704,
-   kHwBeautyMCATNLO = -1705
+   kHwBeautyMCATNLO = -1705,
+   kHwJetsMCATNLO   = -1396
 } Process_t;
 
 class TObjArray;
@@ -155,7 +156,7 @@ typedef struct {
   double VFCH[2][16];
   double VCKM[3][3];
   double VGCUT;
-  double VQCUT;   
+  double VQCUT;
   double VPCUT;
   double ZBINM;
   double EFFMIN;
@@ -456,13 +457,13 @@ typedef struct {
   double COTB;
   double ZMIXSS[4][4];
   double ZMXNSS[4][4];
-  double ZSGNSS[4]; 
+  double ZSGNSS[4];
   double LFCH[16];
   double RFCH[16];
   double SLFCH[4][16];
-  double SRFCH[4][16]; 
+  double SRFCH[4][16];
   double WMXUSS[2][2];
-  double WMXVSS[2][2]; 
+  double WMXVSS[2][2];
   double WSGNSS[2];
   double QMIXSS[2][2][6];
   double LMIXSS[2][2][6];
@@ -676,7 +677,7 @@ typedef struct {
   int     LHNEVT[MAXHRP];
   int     ITYPLH;
   int     LHSOFT;
-  int     LHGLSF;    
+  int     LHGLSF;
 } Hwgupr_t;
 
 typedef struct {
@@ -714,6 +715,12 @@ extern "C" {
   void  hwiodk_(int);
 }
 
+// JIMMY4.2
+extern "C" {
+  void  jimmin_();
+  void  jminit_();
+  void  jmefin_();
+}
 
 
 
@@ -740,9 +747,9 @@ public:
   int         GetNhep          () const     { return fHepevt->NHEP; }
   int         GetISTHEP    (int i)const     { return fHepevt->ISTHEP[i-1]; }
   int         GetIDHEP     (int i)const     { return fHepevt->IDHEP[i-1]; }
-  int         GetJMOHEP (int i, int j) const 
+  int         GetJMOHEP (int i, int j) const
     { return fHepevt->JMOHEP[i-1][j-1]; }
-  int         GetJDAHEP (int i, int j) const 
+  int         GetJDAHEP (int i, int j) const
     { return fHepevt->JDAHEP[i-1][j-1]; }
   double      GetPHEP   (int i, int j) const
     { return fHepevt->PHEP[i-1][j-1]; }
@@ -762,8 +769,8 @@ public:
   Hwbmch_t*   GetHwbmch        ()           { return fHwbmch; }
   char*       GetPART1         () const     { return fHwbmch->PART1; }
   char*       GetPART2         () const     { return fHwbmch->PART2; }
-  
-  
+
+
   // /HWPROC/
   Hwproc_t*   GetHwproc        ()           { return fHwproc; }
   double      GetEBEAM1        () const     { return fHwproc->EBEAM1; }
@@ -798,7 +805,7 @@ public:
   double      GetPTRMS         () const     { return fHwpram->PTRMS; }
   void        SetPTRMS    (double p)        { fHwpram->PTRMS = p; }
   double      GetENSOF         () const     { return fHwpram->ENSOF; }
-  void        SetENSOF    (double e)        { fHwpram->ENSOF = e; } 
+  void        SetENSOF    (double e)        { fHwpram->ENSOF = e; }
   int         GetIPRINT        () const     { return fHwpram->IPRINT; }
   void        SetIPRINT   (int i)           { fHwpram->IPRINT = i; }
   int         GetMODPDF   (int i) const     { return fHwpram->MODPDF[i-1];}
@@ -814,7 +821,7 @@ public:
 
   // /HWPART/
   Hwpart_t*   GetHwpart        ()           { return fHwpart; }
-  
+
   // /HWPARP/
   Hwparp_t*   GetHwparp        ()           { return fHwparp; }
 
@@ -832,10 +839,10 @@ public:
   double      GetAVWGT         () const     { return fHwevnt->AVWGT; }
   int         GetMAXPR         () const     { return fHwevnt->MAXPR; }
   void        SetMAXPR    (int i)           { fHwevnt->MAXPR = i; }
-  
+
   void        SetEV1PR    (int i)           { fHwevnt->EV1PR = i; }
   void        SetEV2PR    (int i)           { fHwevnt->EV2PR = i; }
-  
+
   int         GetMAXER         () const     { return fHwevnt->MAXER; }
   void        SetMAXER    (int i)           { fHwevnt->MAXER = i; }
   int         GetNRN      (int i) const     { return fHwevnt->NRN[i-1]; }
@@ -843,9 +850,9 @@ public:
   double      GetEVWGT         () const     { return fHwevnt->EVWGT; }
 
   int         GetIDHW     (int i) const     { return fHwevnt->IDHW[i]; }
-  
+
   int         GetIERROR        () const     { return fHwevnt->IERROR; }
-  
+
   // /HWHARD/
   Hwhard_t*   GetHwhard        ()           { return fHwhard; }
   double      GetPTMIN         () const     { return fHwhard->PTMIN; }
@@ -877,7 +884,7 @@ public:
 
   void        GetRNAME (int i, char a[9])   { for (int j=0;j<8;j++) a[j] = fHwunam->RNAME[i][j]; a[8] = '\0';}
 /*  char*       GetRNAME(int i) { return fHwunam->RNAME[i]; }*/
-  
+
   // /HWUPDT/
   Hwupdt_t*   GetHwupdt        ()           { return fHwupdt; }
 
@@ -908,8 +915,8 @@ public:
   // /HWCLUS/
   Hwclus_t*   GetHwclus        ()           { return fHwclus; }
 
-  // Herwig6 routines  
-  // the user would call 
+  // Herwig6 routines
+  // the user would call
   //   Initialize
   //   change by himself the parameters s/he wants
   //   Hwusta to make stable the particles s/he wants
@@ -919,7 +926,9 @@ public:
 
   void             GenerateEvent();
   void             Initialize(const char *beam, const char *target, double pbeam1, double pbeam2, int iproc);
+  void             InitializeJimmy(const char *beam, const char *target, double pbeam1, double pbeam2, int iproc);
   void             PrepareRun();
+  void             PrepareRunJimmy();
   void             OpenFortranFile(int lun, char* name);
   void             CloseFortranFile(int lun);
   Int_t            ImportParticles(TClonesArray *particles, Option_t *option="");
@@ -942,6 +951,10 @@ public:
   void             Hwefin();
   void             Hwiodk(int iopt);
   void             SetupTest();
+ // Jimmy subroutines:
+  void             Jminit();
+  void             Jimmin();
+  void             Jmefin();
 protected:
 
   Hepevt_t* fHepevt; // Standard hep common block
@@ -955,14 +968,14 @@ protected:
   Hwparp_t* fHwparp; // Parton polarization common
   Hwbosc_t* fHwbosc; // Electroweak boson common
   Hwparc_t* fHwparc; // Parton colour common
-  Hwbrch_t* fHwbrch; // Branching common 
-  Hwevnt_t* fHwevnt; // Event common 
-  Hwhard_t* fHwhard; // Hard subprocess common 
+  Hwbrch_t* fHwbrch; // Branching common
+  Hwevnt_t* fHwevnt; // Event common
+  Hwhard_t* fHwhard; // Hard subprocess common
   Hwprop_t* fHwprop; // Particle properties
   Hwunam_t* fHwunam; // Particle properties
   Hwupdt_t* fHwupdt; // Particle decays
   Hwuwts_t* fHwuwts; // Weights used in cluster decays
-  Hwuclu_t* fHwuclu; // Parameters for cluster decays 
+  Hwuclu_t* fHwuclu; // Parameters for cluster decays
   Hwdist_t* fHwdist; // Variables controling mixing and vertex information
   Hwqdks_t* fHwqdks; // Arrays for temporarily storing heavy-b,c-hadrons decaying partonicaly
   Hwusud_t* fHwusud; // Parameters for Sudakov form factors