herwig.cc

Go to the documentation of this file.

// example of generating events with Herwig using HepMC/HerwigWrapper.h 
// Events are read into the HepMC event record from the FORTRAN HEPEVT 
// common block using the IO_HERWIG strategy.

#include <stdlib.h>
#include <time.h>

#include <stdio.h>
#include <iostream>
#include <fstream>
#include <sstream>

#include "MyHerwigWrapper.h"
#include "HepMC/IO_HERWIG.h"
#include "HepMC/IO_Ascii.h"     
#include "HepMC/GenEvent.h"
#include "HepMC/HEPEVT_Wrapper.h"
#include "HepMC/SimpleVector.h"

#include "TH1.h"
#include "TFile.h"
#include "TMath.h"
#include <TCanvas.h>

#include "TF1.h"
#include <TH2.h>
#include <TStyle.h>
#include "TLegend.h"
#include <TTree.h>
#include <TString.h>

using namespace std;

// include our own classes
#include "../../include/TopAnalysis.h"
#include "../../include/TauAnalysis.h"
#include "../../include/DiJetAnalysis.h"
#include "../../include/WplusJetAnalysis.h"
#include "../../include/ZAnalysis.h"
#include "../../include/UEAnalysis.h"
#include "readherwigConfigFile.h"
#include "../../include/Configuration.h"

// a dummy routine called on the termination of the run
extern "C" void hwaend_ (void); 
extern "C" {
  int hwgethepevtsize_(void);
  void hwsetmodbos_(int& i, int& ivalue);
  double hwgetcs(void);
  double hwgetcserr(void);
}


int main(int argc, const char* argv[]) {

  if(argc <= 1 ) {
    std::cout << "Usage: " << argv[0] << " <Configuration File>" << std::endl;
    return 0;
  }


  Configuration* currentconfig = new Configuration();
  if(currentconfig->read(argv[1]))
    {
      std::cout << " ... skip !" << std::endl;
      return 0;
    }

  DiJetAnalysis dijet;
  WplusJetAnalysis wplusjet;
  ZAnalysis z0;
  TauAnalysis tau;
  TopAnalysis top;
  UEAnalysis ue;

  // create analysis Objects and initialize them
  if(currentconfig->dijet){
    dijet.Init(currentconfig->max_eta,currentconfig->min_pt);
    dijet.InitJetFinder(currentconfig->coneRadius,currentconfig->overlapThreshold,currentconfig->jet_ptmin);
  }

  if(currentconfig->wplusjet){
    wplusjet.Init(currentconfig->max_eta,currentconfig->min_pt);
    wplusjet.InitJetFinder(currentconfig->coneRadius,currentconfig->overlapThreshold,currentconfig->jet_ptmin);
  }

  if(currentconfig->z){
    z0.Init(currentconfig->max_eta,currentconfig->min_pt);
    z0.InitJetFinder(currentconfig->coneRadius,currentconfig->overlapThreshold,currentconfig->jet_ptmin);
  }
  if(currentconfig->tau){
    tau.Init(currentconfig->max_eta,currentconfig->min_pt);
    tau.InitJetFinder(currentconfig->coneRadius,currentconfig->overlapThreshold,currentconfig->jet_ptmin);
  }
  if(currentconfig->top){
    top.Init(currentconfig->max_eta,currentconfig->min_pt);
    top.InitJetFinder(currentconfig->coneRadius,currentconfig->overlapThreshold,currentconfig->jet_ptmin);
  }
  if(currentconfig->ue){
    ue.Init(currentconfig->max_eta,currentconfig->min_pt);
    ue.InitJetFinder(currentconfig->coneRadius,currentconfig->overlapThreshold,currentconfig->jet_ptmin);
  }

    //
    //........................................HEPEVT
    // herwig-6510-2 uses HEPEVT with 10000 entries and 8-byte floating point
    //  numbers. We need to explicitly pass this information to the 
    //  HEPEVT_Wrapper.
    //
    HepMC::HEPEVT_Wrapper::set_max_number_entries(hwgethepevtsize_());
    HepMC::HEPEVT_Wrapper::set_sizeof_real(8);

    //random seed
     hwevnt.nrn[0]= currentconfig->rseed;

    // number of events
    hwproc.maxev = int(currentconfig->nevents); 

    // tell it what the beam particles are:
    for ( unsigned int i = 0; i < 8; ++i ) {
        hwbmch.part1[i] = (i < 1) ? 'P' : ' ';
        hwbmch.part2[i] = (i < 1) ? 'P' : ' ';
    }

    std::cout << "Configuration Parameter List:" << std::endl;
    for (unsigned int files=0; files<(currentconfig->ConfigFileNames).size(); files++)
      readConfigFile((currentconfig->ConfigFileNames)[files]);

    // master switches for ISR 
    hwpram.nospac=int(currentconfig->ISR); //switch on/off ISR

   // INITIALISE OTHER COMMON BLOCKS
    hwigin();
   // number of events to print
    hwevnt.maxpr = 1; 
   // compute parameter-dependent constants
    hwuinc(); 
   // initialise elementary process
    hweini(); 

    //........................................HepMC INITIALIZATIONS
    //
    // Instantiate an IO strategy for reading from HEPEVT.
    HepMC::IO_HERWIG hepevtio;
    //
    //........................................EVENT LOOP
    for ( int i = 1; i <= hwproc.maxev; i++ ) {
        if ( i%50==1 ) std::cout << "Processing Event Number " 
                                 << i << std::endl;
        // initialise event
        hwuine();
        // generate hard subproces
        hwepro();
        // generate parton cascades
        hwbgen();
        // do heavy object decays
        hwdhob();
        // do cluster formation
        hwcfor();
        // do cluster decays
        hwcdec();
        // do unstable particle decays
        hwdhad();
        // do heavy flavour hadron decays
        hwdhvy();
        // add soft underlying event if needed
        hwmevt();
        // finish event
        hwufne();

        HepMC::GenEvent* evt = hepevtio.read_next_event();
        // add some information to the event
        evt->set_event_number(i);
        evt->set_signal_process_id(20);

    int ret_dijet = SUCCESS;
    int ret_wplusjet = SUCCESS;
    int ret_z0 = SUCCESS;
    int ret_tau = SUCCESS;
    int ret_top = SUCCESS;
    int ret_ue = SUCCESS;

    //call analysis
    if(currentconfig->dijet)       ret_dijet = dijet.Process(evt);
    if(currentconfig->wplusjet)    ret_wplusjet = wplusjet.Process(evt);
    if(currentconfig->z)           ret_z0 = z0.Process(evt);
    if(currentconfig->tau)  ret_tau = tau.Process(evt);
    if(currentconfig->top)  ret_top = top.Process(evt);
    if(currentconfig->ue)   ret_ue = ue.Process(evt);


    if (ret_dijet==FAILURE || ret_wplusjet==FAILURE || ret_z0==FAILURE
        || ret_tau==FAILURE || ret_top==FAILURE || ret_ue==FAILURE)
      {
        std::cout << "Processing event number " << i << std::endl;
        std::cout << "throw FAILURE"<<std::endl;
        //std::cout<<"pylist"<<std::endl;
       // call_pylist(2);
        std::cout<<std::endl;
        //std::cout<<"HepMC Output"<<std::endl;
       // evt->print();
        std::cout<<std::endl;
      }


        if (i<=hwevnt.maxpr) {
            std::cout << "\n\n This is the FIXED version of HEPEVT as "
                      << "coded in IO_HERWIG " << std::endl;
        //    HepMC::HEPEVT_Wrapper::print_hepevt();
        //    evt->print();
        }

        //delete all the jet objects to be free for the next event
       //if(currentconfig->dijet)       dijet.DeleteJetObject(evt);
        //if(currentconfig->wplusjet)    wplusjet.DeleteJetObject(evt);
        //if(currentconfig->z)         z0.DeleteJetObject(evt);
        //if(currentconfig->tau)       tau.DeleteJetObject(evt);
        //if(currentconfig->top)       top.DeleteJetObject(evt);
        //if(currentconfig->ue)        ue.DeleteJetObject(evt);

        // we also need to delete the created event from memory
        delete evt;
    }
    //........................................TERMINATION
    hwefin();

  TFile* f = new TFile(currentconfig->OutputFileName.c_str(),"RECREATE");
  if(currentconfig->dijet)       dijet.Finalize(f);
  if(currentconfig->wplusjet)    wplusjet.Finalize(f);
  if(currentconfig->z)         z0.Finalize(f);
  if(currentconfig->tau)       tau.Finalize(f);
  if(currentconfig->top)       top.Finalize(f);
  if(currentconfig->ue)        ue.finalize(f);

  std::cout<<"Herwig:  program ends"<<std::endl;
  f->Close();
 
 // Done.
    return 0;
}

// a dummy user-supplied subroutine called when the run terminates:
void hwaend_ (void) { return; }

---