herwig.cc

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////
// Matt.Dobbs@Cern.CH, October 2002
// 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>

// include our own classes
#include "../../include/Configuration.h"
#include "../../include/JetFinder.h"
#include "../../include/ttbarAnalysis.h"
#include "../../include/ZtautauAnalysis.h"
#include "../../include/WtaunuAnalysis.h"
#include "../../include/bbbarAnalysis.h"
#include "../../include/JetAnalysis.h"
#include "../../include/WplusJetAnalysis.h"
#include "../../include/ZAnalysis.h"
#include "../../include/UEAnalysis.h"
#include "../../include/EtMissAnalysis.h"
#include "../../include/ElasScatAnalysis.h"
#include "../../include/UserAnalysis.h"
#include "readherwigConfigFile.h"

using namespace std;

// 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 ) {
    cout << "Usage: " << argv[0] << " <Configuration File>" << endl;
    return 0;
  }

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

  JetFinder FindJet;
  FindJet.Init( currentconfig->max_eta, currentconfig->min_pt );
  FindJet.InitJetFinder( currentconfig->coneRadius, currentconfig->overlapThreshold, currentconfig->jet_ptmin,
    currentconfig->lepton_ptmin, currentconfig->DeltaR_lepton_track, currentconfig->jetalgorithm );

    // create analysis Objects and initialize them                                                                                          
  vector<baseAnalysis*> analysis;

  if (currentconfig->jet) {
    cout << "Adding module JetAnalysis" << endl;
    analysis.push_back(new JetAnalysis);
  }
  if (currentconfig->wplusjet) {
    cout << "Adding module WplusJetAnalysis" << endl;
    analysis.push_back(new WplusJetAnalysis);
  }
  if(currentconfig->z){
    cout << "Adding module ZAnalysis" << endl;
    analysis.push_back(new ZAnalysis);
  }
  if(currentconfig->ztautau){
    cout << "Adding module ZtautauAnalysis" << endl;
    analysis.push_back(new ZtautauAnalysis);
  }
  if(currentconfig->wtaunu){
    cout << "Adding module WtaunuAnalysis" << endl;
    analysis.push_back(new WtaunuAnalysis);
  }
  if(currentconfig->bbbar){
    cout << "Adding module bbbarAnalysis" << endl;
    analysis.push_back(new bbbarAnalysis);
  }
  if(currentconfig->ttbar){
    cout << "Adding module ttbarAnalysis" << endl;
    analysis.push_back(new ttbarAnalysis);
  }
  if(currentconfig->ue){
    cout << "Adding module UEAnalysis" << endl;
    analysis.push_back(new UEAnalysis);
  }  
  if(currentconfig->etmiss){
    cout << "Adding module EtMissAnalysis" << endl;
    analysis.push_back(new EtMissAnalysis);
  }
  if(currentconfig->elasScat){
    cout << "Adding module ElasScatAnalysis" << endl;
    analysis.push_back(new ElasScatAnalysis);
  }
  if(currentconfig->user){
    cout << "Adding module UserAnalysis" << endl;
    analysis.push_back(new UserAnalysis);
  }

  // init requested analysis modules and jet finder
  for (vector<baseAnalysis*>::const_iterator i(analysis.begin()); i != analysis.end(); ++i) {
    (*i)->Init(currentconfig->max_eta,currentconfig->min_pt);
    (*i)->InitJetFinder( currentconfig->coneRadius, currentconfig->overlapThreshold, currentconfig->jet_ptmin,
      currentconfig->lepton_ptmin, currentconfig->DeltaR_lepton_track, currentconfig->jetalgorithm );
  }


  //
  //........................................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' : ' ';
  }

  // INITIALISE OTHER COMMON BLOCKS
  hwigin();
  jimmin();

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

  // master switches for ISR
  // true switches ISR off!
  //  hwpram.nospac=int(currentconfig->ISR); //switch on/off ISR
    
  // number of events to print
  hwevnt.maxpr = 1; 
  // compute parameter-dependent constants
  hwuinc(); 
  // initialise elementary process
  hweini(); 
  // initialise jimmy event
  jminit();
  //........................................HepMC INITIALIZATIONS
  //
  // Instantiate an IO strategy for reading from HEPEVT.
  HepMC::IO_HERWIG hepevtio;
  //
  //........................................EVENT LOOP
  for ( int iEvent = 1; iEvent <= hwproc.maxev; iEvent++ ) {
    if ( iEvent%50==1 ) cout << "Processing Event Number " 
                                  << iEvent << 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();

    // initialise event
    hwuine();
    // generate hard subproces
    hwepro();
    // generate parton cascades
    hwbgen();

    // Do Jimmy underlying event if required.
    int abort = 0;
    hwmsct(&abort);
    if (abort != 0) {
      cout << "Event aborted by hwmsct." << endl;
      hwufne();  // finish event
      continue;
    }

    //    if (abort == 0) {
      hwdhob();  // Do heavy quark decays
      hwcfor();  // Do cluster formation
      hwcdec();  // Do cluster decays
      hwdhad();  // Do unstable particle decays
      hwdhvy();  // Do heavy flavor decays
      hwmevt();  // Add soft underlying event
      hwufne();  // finish event
    //}
      if (hwevnt.ierror==-1 || hwevnt.ierror>99) {
        cout << "Event Error, Number: " << " Code: " << hwevnt.ierror << endl;
        continue;
      }
    
    // Print HEPEVT format event here 
    // HepMC::HEPEVT_Wrapper::print_hepevt();

    HepMC::GenEvent* evt = hepevtio.read_next_event();
    // add some information to the event
    evt->set_event_number(iEvent);
    evt->set_signal_process_id(20);
   
    // the status codes should not be changed, otherwise you change
    // the output of the MC generators (status code variabl) and this
    // should not be done
    //     //change the status code in HepMC format event.  
    //     for ( HepMC::GenEvent::particle_const_iterator p =
    //      evt->particles_begin(); p != evt->particles_end(); ++p ){
    //       if( abs((*p)->status())>1 )
    //  (*p)->set_status(2);
    //     }
       
    //call the Jetfinder
    vector<fastjet::PseudoJet> inclusive_jets = FindJet.GetJet(evt);

    //call analysis
    int ret_all = true;
    for (vector<baseAnalysis*>::const_iterator i(analysis.begin()); i != analysis.end(); ++i) {
      (*i)->SetJet(&inclusive_jets);
      int ret = (*i)->Process(evt);
      if (ret == false) ret_all = false;
    }

    if (ret_all==false) {
      cout << "Processing event number " << iEvent << endl;
      cout << "throw FAILURE"<<endl;
      //cout<<"pylist"<<endl;
      // call_pylist(2);
      cout<<endl;
      //cout<<"HepMC Output"<<endl;
      // evt->print();
      cout<<endl;
    }


    if (iEvent<=hwevnt.maxpr) {
      cout << "\n\n This is the FIXED version of HEPEVT as "
                << "coded in IO_HERWIG " << endl;
      //    HepMC::HEPEVT_Wrapper::print_hepevt();
      //    evt->print();
    }
    
    //delete all the jet objects, missingEt etc. to be free for the next event
    FindJet.ClearEvent(evt);
      
    // we also need to delete the created event from memory
    delete evt;

  }
  //........................................TERMINATION
  hwefin();

 
  TFile* f = new TFile(currentconfig->OutputFileName.c_str(),"RECREATE");
  for (vector<baseAnalysis*>::const_iterator i(analysis.begin()); i != analysis.end(); ++i) {
    (*i)->Finalize(f);
  }

  // clean up analysis modules
  for (vector<baseAnalysis*>::const_iterator i(analysis.begin()); i != analysis.end(); ++i) {
    delete (*i);
  }
  analysis.resize(0);

  // clean up config object
  delete currentconfig;

  // Tidy up
  cout<<"Herwig:  program ends"<<endl;
  f->Close();
 
  cout<<"Run successfully!"<<endl;
  return 0;
}

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

---