herwigpp_2.4.2.cc

Go to the documentation of this file.
00001 //example of Herwig++, generate pp -> ttbar events
00002 //by Zhonghua Qin <zhonghua.qin@desy.de> 
00003 #include "ThePEG/Repository/EventGenerator.h"
00004 #include "ThePEG/Repository/Repository.h"
00005 #include "ThePEG/Persistency/PersistentIStream.h"
00006 #include "ThePEG/Utilities/DynamicLoader.h"
00007 #include "ThePEG/EventRecord/Event.h"
00008 #include "ThePEG/EventRecord/SubProcess.h"
00009 #include "ThePEG/Handlers/XComb.h"
00010 #include "ThePEG/Handlers/EventHandler.h"
00011 #include "ThePEG/PDF/PartonExtractor.h"
00012 #include "ThePEG/PDF/PDF.h"
00013 #include "ThePEG/Vectors/HepMCConverter.h"
00014 #include "ThePEGHepMC_1.5.0.h"
00015 
00016 #include "HepMC/GenEvent.h"
00017 #include "HepMC/IO_GenEvent.h"
00018 #include "HepMC/GenEvent.h"
00019 #include "HepMC/IO_AsciiParticles.h"
00020 #include "HepMC/SimpleVector.h"
00021 
00022 #include <iostream>
00023 #include <stdio.h>
00024 #include <string>
00025 
00026 #include "TH1.h"
00027 #include "TFile.h"
00028 #include "TMath.h"
00029 #include <TCanvas.h>
00030 
00031 #include "TF1.h"
00032 #include <TH2.h>
00033 #include <TStyle.h>
00034 #include "TLegend.h"
00035 #include <TTree.h>
00036 #include <TString.h>
00037 
00038 using namespace ThePEG;
00039 
00040 // include our own classes
00041 #include "../../include/ttbarAnalysis.h"
00042 #include "../../include/ZtautauAnalysis.h"
00043 #include "../../include/WtaunuAnalysis.h"
00044 #include "../../include/bbbarAnalysis.h"
00045 #include "../../include/JetAnalysis.h"
00046 #include "../../include/WplusJetAnalysis.h"
00047 #include "../../include/ZAnalysis.h"
00048 #include "../../include/UEAnalysis.h"
00049 #include "../../include/EtMissAnalysis.h"
00050 #include "../../include/ElasScatAnalysis.h"
00051 #include "../../include/UserAnalysis.h"
00052 #include "readHerwigConfigFile.h"
00053 #include "../../include/Configuration.h"
00054 #include "../../include/JetFinder.h"
00055 
00056 int main(int argc, const char* argv[]) {
00057 
00058   if(argc <= 1 ) {
00059       std::cout << "Usage: " << argv[0] << " <Configuration File>" << std::endl;
00060       return 0;
00061   }
00062   
00063   
00064   Configuration* currentconfig = new Configuration();
00065   if(currentconfig->read(argv[1]))
00066       {
00067           std::cout << " ... skip !" << std::endl;
00068           return 0;
00069       }
00070 
00071   JetFinder FindJet;
00072   FindJet.Init(currentconfig->max_eta,currentconfig->min_pt);
00073   FindJet.InitJetFinder( currentconfig->coneRadius, currentconfig->overlapThreshold, currentconfig->jet_ptmin,
00074     currentconfig->lepton_ptmin,currentconfig->DeltaR_lepton_track, currentconfig->jetalgorithm );
00075 
00076 
00077   std::vector<baseAnalysis*> analysis;
00078 
00079   // create analysis Objects and initialize them                                                                                          
00080   if (currentconfig->jet) {
00081     std::cout << "Adding module JetAnalysis" << std::endl;
00082     analysis.push_back(new JetAnalysis);
00083   }
00084   if (currentconfig->wplusjet) {
00085     std::cout << "Adding module WplusJetAnalysis" << std::endl;
00086     analysis.push_back(new WplusJetAnalysis);
00087   }
00088   if(currentconfig->z){
00089     std::cout << "Adding module ZAnalysis" << std::endl;
00090     analysis.push_back(new ZAnalysis);
00091   }
00092   if(currentconfig->ztautau){ 
00093     std::cout << "Adding module ZtautauAnalysis" << std::endl;
00094     analysis.push_back(new ZtautauAnalysis);
00095   }
00096   if(currentconfig->wtaunu){ 
00097     std::cout << "Adding module WtaunuAnalysis" << std::endl;
00098     analysis.push_back(new WtaunuAnalysis);
00099   }
00100   if(currentconfig->bbbar){ 
00101     std::cout << "Adding module bbbarAnalysis" << std::endl;
00102     analysis.push_back(new bbbarAnalysis);
00103   }
00104   if(currentconfig->ttbar){
00105     std::cout << "Adding module ttbarAnalysis" << std::endl;
00106     analysis.push_back(new ttbarAnalysis);
00107   }
00108   if(currentconfig->ue){
00109     std::cout << "Adding module UEAnalysis" << std::endl;
00110     analysis.push_back(new UEAnalysis);
00111   }  
00112   if(currentconfig->etmiss){
00113     std::cout << "Adding module EtMissAnalysis" << std::endl;
00114     analysis.push_back(new EtMissAnalysis);
00115   }
00116   if(currentconfig->elasScat){
00117     std::cout << "Adding module ElasScatAnalysis" << std::endl;
00118     analysis.push_back(new ElasScatAnalysis);
00119   }
00120   if(currentconfig->user){
00121     std::cout << "Adding module UserAnalysis" << std::endl;
00122     analysis.push_back(new UserAnalysis);
00123   }
00124   
00125 
00126   // init requested analysis modules and jet finder
00127   for (std::vector<baseAnalysis*>::const_iterator i(analysis.begin()); i != analysis.end(); ++i) {
00128     (*i)->Init(currentconfig->max_eta,currentconfig->min_pt);
00129     (*i)->InitJetFinder( currentconfig->coneRadius, currentconfig->overlapThreshold, currentconfig->jet_ptmin,
00130         currentconfig->lepton_ptmin,currentconfig->DeltaR_lepton_track, currentconfig->jetalgorithm );
00131   }
00132 
00133   //Add a path to seach for dynamically linkable libraries
00134   std::cout<<"Add a path to seach for dynamically linkable libraries..."<<std::endl;
00135   ThePEG::DynamicLoader::appendPath(HWMODULEDIR);
00136   ThePEG::DynamicLoader::appendPath(THEPEGMODULEDIR);
00137   
00138   //Load a whole repository from the given file
00139   std::cout<<"Load repository from a given file..."<<std::endl;
00140   const string repopath = string(HWREPODIR) + "/HerwigDefaults.rpo";
00141   ThePEG::Repository::load(repopath);
00142 
00143   //Setup process,interpret the command in cmd and return possible messages
00144   std::cout<<"Setup process,interpret the command in cmd and return possible messages..."<<std::endl;
00145 
00146   // set random number seed (mandatory!) ... this may be overwritten by function readConfigFile
00147   char exec_command[100];
00148   sprintf(exec_command,"set LHCGenerator:RandomNumberGenerator:Seed %ld",currentconfig->rseed); 
00149   ThePEG::Repository::exec("cd /Herwig/Generators", cout);
00150   ThePEG::Repository::exec(exec_command, cout);
00151 
00152   // vector to store herwig commandlines from Configuration file
00153   std::vector< std::string > herwig_commands;
00154   // read in commandlines
00155   for (unsigned int files=0; files<(currentconfig->ConfigFileNames).size(); files++)
00156       readConfigFile(&herwig_commands,(currentconfig->ConfigFileNames)[files]);
00157 
00158   // loop over the commands
00159   for(unsigned int i=0; i<herwig_commands.size(); i++)
00160       ThePEG::Repository::exec( herwig_commands[i].c_str(), cout);
00161 
00162 //switches for FSR, ISR and MI
00163   ThePEG::Repository::exec("cd /Herwig/Shower", cout);
00164 
00165   sprintf(exec_command,"set SplittingGenerator:ISR %c",currentconfig->ISR);
00166   ThePEG::Repository::exec(exec_command, cout);
00167 
00168   sprintf(exec_command,"set SplittingGenerator:FSR %c",currentconfig->FSR);
00169   ThePEG::Repository::exec(exec_command, cout);
00170 
00171   sprintf(exec_command,"set ShowerHandler:MPI %c",currentconfig->MI);
00172   ThePEG::Repository::exec(exec_command, cout);
00173 
00174   //Check sanity of the object during the setup phase
00175   ThePEG::Repository::update();
00176 
00177  // ThePEG generator object.
00178   ThePEG::EGPtr m_hw;
00179 
00180   // ThePEG event object.
00181   ThePEG::EventPtr m_event;
00182 
00183  // Make a "run" object from the config repository.
00184   std::cout<<"make a run object ..."<<std::endl;
00185   ThePEG::EGPtr tmpEG = ThePEG::BaseRepository::GetObject<ThePEG::EGPtr>("/Herwig/Generators/LHCGenerator");
00186   m_hw = ThePEG::Repository::makeRun(tmpEG, "LHC");
00187 
00188  //Herwig++ initializing
00189   std::cout<<"Herwig++ initializing..."<<std::endl;
00190   m_hw->initialize();
00191 
00192  //Loop for events generating and fill histograms
00193   std::cout<<"Generating events and filling histograms..."<<std::endl;
00194   
00195  //Change the value of nEvent to the number of events you want to generate
00196   for (unsigned int iEvent = 0; iEvent < currentconfig->nevents; ++iEvent) {
00197 
00198     if ( iEvent % 100 == 0 ) {
00199       std::cout << "Event number " << iEvent << std::endl;
00200     }
00201  // Run the generator for one event.
00202   assert(m_hw);
00203   m_event = m_hw->shoot();
00204 
00205  // Fill HepMC event from Herwig's internally stored EventPtr.
00206  // Convert the Herwig event into the HepMC GenEvent
00207  // HepMC::GenEvent*  evt = new HepMC::GenEvent();
00208   //  ThePEG::HepMCConverter<HepMC::GenEvent>::convert(*m_event, *evt, false, ThePEG::GeV);
00209   HepMC::GenEvent*  evt= ThePEG::HepMCConverter<HepMC::GenEvent>::convert(*m_event, false, ThePEG::GeV);
00210 
00211   //call the Jetfinder
00212   vector<fastjet::PseudoJet> inclusive_jets = FindJet.GetJet(evt);
00213   
00214   //call analysis
00215   int ret_all = true;
00216   for (std::vector<baseAnalysis*>::const_iterator i(analysis.begin()); i != analysis.end(); ++i) {
00217     (*i)->SetJet(&inclusive_jets);
00218     int ret = (*i)->Process(evt);
00219     if (ret == false) ret_all = false;
00220   }
00221 
00222   if (ret_all==false) {
00223     std::cout << "Processing event number " << iEvent << std::endl;
00224     std::cout << "throw FAILURE"<<std::endl;
00225     //std::cout<<"HepMC Output"<<std::endl;
00226     //evt->print();
00227     std::cout<<std::endl;
00228   }
00229 
00230   //delete all the jet objects, missingEt etc. to be free for the next event
00231   FindJet.ClearEvent(evt);
00232   
00233   //delete event
00234   delete evt;
00235   }
00236 
00237   TFile* f = new TFile(currentconfig->OutputFileName.c_str(),"RECREATE");
00238   for (std::vector<baseAnalysis*>::const_iterator i(analysis.begin()); i != analysis.end(); ++i) {
00239     (*i)->Finalize(f);
00240   }
00241 
00242   // clean up analysis modules
00243   for (std::vector<baseAnalysis*>::const_iterator i(analysis.begin()); i != analysis.end(); ++i) {
00244     delete (*i);
00245   }
00246   analysis.resize(0);
00247 
00248   // clean up config object
00249   delete currentconfig;
00250 
00251   //Tidy up, Herwig++ finalizing
00252   f->Close();
00253   std::cout<<"Herwig++ finalizing..."<<std::endl;
00254   assert(m_hw);
00255   m_hw->finalize();
00256   ThePEG::Repository::cleanup();
00257 
00258   std::cout<<"Run successfully!"<<std::endl;
00259   return 0;
00260 }

Generated on Wed Aug 31 09:44:48 2011 for HepMCAnalysis by  doxygen 1.4.7