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>

using namespace std;

// include our own classes
#include "../../include/Configuration.h"
#include "../../include/JetFinder.h"
#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 "../../include/EtMissAnalysis.h"
#include "../../include/ElasScatAnalysis.h"
#include "../../include/UserAnalysis.h"
#include "readherwigConfigFile.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;
    }

  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);


  std::vector<baseAnalysis*> analysis;

  // create analysis Objects and initialize them                                                                                          
  if (currentconfig->dijet) {
    std::cout << "Adding module DiJetAnalysis" << std::endl;
    analysis.push_back(new DiJetAnalysis);
  }
  if (currentconfig->wplusjet) {
    std::cout << "Adding module WplusJetAnalysis" << std::endl;
    analysis.push_back(new WplusJetAnalysis);
  }
  if(currentconfig->z){
    std::cout << "Adding module ZAnalysis" << std::endl;
    analysis.push_back(new ZAnalysis);
  }
  if(currentconfig->tau){
    std::cout << "Adding module TauAnalysis" << std::endl;
    analysis.push_back(new TauAnalysis);
  }
  if(currentconfig->top){
    std::cout << "Adding module TopAnalysis" << std::endl;
    analysis.push_back(new TopAnalysis);
  }
  if(currentconfig->ue){
    std::cout << "Adding module UEAnalysis" << std::endl;
    analysis.push_back(new UEAnalysis);
  }  
  if(currentconfig->etmiss){
    std::cout << "Adding module EtMissAnalysis" << std::endl;
    analysis.push_back(new EtMissAnalysis);
  }
  if(currentconfig->elasScat){
    std::cout << "Adding module ElasScatAnalysis" << std::endl;
    analysis.push_back(new ElasScatAnalysis);
  }
  if(currentconfig->user){
    std::cout << "Adding module UserAnalysis" << std::endl;
    analysis.push_back(new UserAnalysis);
  }
  

  // init requested analysis modules and jet finder
  for (std::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);
  }


  //
  //........................................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(); 
  // 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 ) std::cout << "Processing Event Number " 
                                  << iEvent << 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();

    // initialise event
    hwuine();
    // generate hard subproces
    hwepro();
    // generate parton cascades
    hwbgen();
    // Do Jimmy underlying event if required.
    int abort = 0;
    if (abort == 0) {
      hwmsct(&abort);
      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
    }
    
    // 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
    std::vector<fastjet::PseudoJet> inclusive_jets = FindJet.GetJet(evt);

    //call analysis
    int ret_all = true;
    for (std::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) {
      std::cout << "Processing event number " << iEvent << 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 (iEvent<=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, 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 (std::vector<baseAnalysis*>::const_iterator i(analysis.begin()); i != analysis.end(); ++i) {
    (*i)->Finalize(f);
  }

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

  // clean up config object
  delete currentconfig;

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

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

---