herwigpp_2.3.1.cc

Go to the documentation of this file.

//example of Herwig++, generate pp -> ttbar events
//by Zhonghua Qin <zhonghua.qin@desy.de> 
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Repository/Repository.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Utilities/DynamicLoader.h"
#include "ThePEG/EventRecord/Event.h"
#include "ThePEG/EventRecord/SubProcess.h"
#include "ThePEG/Handlers/XComb.h"
#include "ThePEG/Handlers/EventHandler.h"
#include "ThePEG/PDF/PartonExtractor.h"
#include "ThePEG/PDF/PDF.h"
#include "ThePEG/Vectors/HepMCConverter.h"
#include "ThePEGHepMC_1.4.0.h"

#include "HepMC/GenEvent.h"
#include "HepMC/IO_GenEvent.h"
#include "HepMC/GenEvent.h"
#include "HepMC/IO_AsciiParticles.h"
#include "HepMC/SimpleVector.h"

#include <iostream>
#include <stdio.h>
#include <string>

#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 ThePEG;

// 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"

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


  // Add a path to seach for dynamically linkable libraries
  std::cout<<"Add a path to seach for dynamically linkable libraries..."<<std::endl;
  ThePEG::DynamicLoader::appendPath(HWMODULEDIR);
  ThePEG::DynamicLoader::appendPath(THEPEGMODULEDIR);
  
  //Load a whole repository from the given file
  std::cout<<"Load repository from a given file..."<<std::endl;
  const string repopath = string(HWREPODIR) + "/HerwigDefaults.rpo";
  ThePEG::Repository::load(repopath);

  //Setup process,interpret the command in cmd and return possible messages
  std::cout<<"Setup process,interpret the command in cmd and return possible messages..."<<std::endl;

  // set random number seed (mandatory!) ... this may be overwritten by function readConfigFile
  char exec_command[100];
  sprintf(exec_command,"set LHCGenerator:RandomNumberGenerator:Seed %ld",currentconfig->rseed); 
  ThePEG::Repository::exec("cd /Herwig/Generators", cout);
  ThePEG::Repository::exec(exec_command, cout);

  // vector to store herwig commandlines from Configuration file
  std::vector< std::string > herwig_commands;
  // read in commandlines
  for (unsigned int files=0; files<(currentconfig->ConfigFileNames).size(); files++)
    readConfigFile(&herwig_commands,(currentconfig->ConfigFileNames)[files]);

  // loop over the commands
  for(unsigned int i=0; i<herwig_commands.size(); i++)
    ThePEG::Repository::exec( herwig_commands[i].c_str(), cout);

  //switches for FSR, ISR and MI
  ThePEG::Repository::exec("cd /Herwig/Shower", cout);
  
  sprintf(exec_command,"set SplittingGenerator:ISR %c",currentconfig->ISR);
  ThePEG::Repository::exec(exec_command, cout);
  
  sprintf(exec_command,"set SplittingGenerator:FSR %c",currentconfig->FSR);
  ThePEG::Repository::exec(exec_command, cout);
  
  sprintf(exec_command,"set ShowerHandler:MPI %c",currentconfig->MI);
  ThePEG::Repository::exec(exec_command, cout);


  //Check sanity of the object during the setup phase
  ThePEG::Repository::update();

  // ThePEG generator object.
  ThePEG::EGPtr m_hw;

  // ThePEG event object.
  ThePEG::EventPtr m_event;

  // Make a "run" object from the config repository.
  std::cout<<"make a run object ..."<<std::endl;
  ThePEG::EGPtr tmpEG = ThePEG::BaseRepository::GetObject<ThePEG::EGPtr>("/Herwig/Generators/LHCGenerator");
  m_hw = ThePEG::Repository::makeRun(tmpEG, "LHC");

  // Herwig++ initializing
  std::cout<<"Herwig++ initializing..."<<std::endl;
  m_hw->initialize();

  // Loop for events generating and fill histograms
  std::cout<<"Generating events and filling histograms..."<<std::endl;
  
  //Change the value of nEvent to the number of events you want to generate
  for (unsigned int iEvent = 0; iEvent < currentconfig->nevents; ++iEvent) {

    // Run the generator for one event.
    assert(m_hw);
    m_event = m_hw->shoot();
    
    // Fill HepMC event from Herwig's internally stored EventPtr.
    // Convert the Herwig event into the HepMC GenEvent
    HepMC::GenEvent*  evt= ThePEG::HepMCConverter<HepMC::GenEvent>::convert(*m_event, false, ThePEG::GeV);

    // debug print out
    std::cout << "Processing event number " << iEvent << std::endl;
    if (iEvent < 0) evt->print();

    //call the Jetfinder
    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<<"HepMC Output"<<std::endl;
      //evt->print();
      std::cout<<std::endl;
    }
    
    // delete all the jet objects, missingEt etc. to be free for the next event
    FindJet.ClearEvent(evt);
    
    // delete event
    delete evt;
  }
  
  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, Herwig++ finalizing
  std::cout<<"Herwig++ finalizing..."<<std::endl;
  assert(m_hw);
  m_hw->finalize();
  ThePEG::Repository::cleanup();

  std::cout<<"Run successfully!"<<std::endl;
  return 0;
}

---