00001
00002
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.4.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
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
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
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
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
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
00144 std::cout<<"Setup process,interpret the command in cmd and return possible messages..."<<std::endl;
00145
00146
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
00153 std::vector< std::string > herwig_commands;
00154
00155 for (unsigned int files=0; files<(currentconfig->ConfigFileNames).size(); files++)
00156 readConfigFile(&herwig_commands,(currentconfig->ConfigFileNames)[files]);
00157
00158
00159 for(unsigned int i=0; i<herwig_commands.size(); i++)
00160 ThePEG::Repository::exec( herwig_commands[i].c_str(), cout);
00161
00162
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
00175 ThePEG::Repository::update();
00176
00177
00178 ThePEG::EGPtr m_hw;
00179
00180
00181 ThePEG::EventPtr m_event;
00182
00183
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
00189 std::cout<<"Herwig++ initializing..."<<std::endl;
00190 m_hw->initialize();
00191
00192
00193 std::cout<<"Generating events and filling histograms..."<<std::endl;
00194
00195
00196 for (unsigned int iEvent = 0; iEvent < currentconfig->nevents; ++iEvent) {
00197
00198
00199 assert(m_hw);
00200 m_event = m_hw->shoot();
00201
00202
00203
00204
00205
00206 HepMC::GenEvent* evt= ThePEG::HepMCConverter<HepMC::GenEvent>::convert(*m_event, false, ThePEG::GeV);
00207
00208
00209 vector<fastjet::PseudoJet> inclusive_jets = FindJet.GetJet(evt);
00210
00211
00212 int ret_all = true;
00213 for (std::vector<baseAnalysis*>::const_iterator i(analysis.begin()); i != analysis.end(); ++i) {
00214 (*i)->SetJet(&inclusive_jets);
00215 int ret = (*i)->Process(evt);
00216 if (ret == false) ret_all = false;
00217 }
00218
00219 if (ret_all==false) {
00220 std::cout << "Processing event number " << iEvent << std::endl;
00221 std::cout << "throw FAILURE"<<std::endl;
00222
00223
00224 std::cout<<std::endl;
00225 }
00226
00227
00228 FindJet.ClearEvent(evt);
00229
00230
00231 delete evt;
00232 }
00233
00234 TFile* f = new TFile(currentconfig->OutputFileName.c_str(),"RECREATE");
00235 for (std::vector<baseAnalysis*>::const_iterator i(analysis.begin()); i != analysis.end(); ++i) {
00236 (*i)->Finalize(f);
00237 }
00238
00239
00240 for (std::vector<baseAnalysis*>::const_iterator i(analysis.begin()); i != analysis.end(); ++i) {
00241 delete (*i);
00242 }
00243 analysis.resize(0);
00244
00245
00246 delete currentconfig;
00247
00248
00249 std::cout<<"Herwig++ finalizing..."<<std::endl;
00250 assert(m_hw);
00251 m_hw->finalize();
00252 ThePEG::Repository::cleanup();
00253
00254 std::cout<<"Run successfully!"<<std::endl;
00255 return 0;
00256 }