herwigpp_2.5.0.cc

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

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