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
00015
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
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
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
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
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
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
00149 std::cout<<"Setup process,interpret the command in cmd and return possible messages..."<<std::endl;
00150
00151
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
00158 std::vector< std::string > herwig_commands;
00159
00160 for (unsigned int files=0; files<(currentConfig->ConfigFileNames).size(); files++)
00161 readConfigFile(&herwig_commands,(currentConfig->ConfigFileNames)[files]);
00162
00163
00164 for(unsigned int i=0; i<herwig_commands.size(); i++)
00165 ThePEG::Repository::exec( herwig_commands[i].c_str(), cout);
00166
00167
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
00180 ThePEG::Repository::update();
00181
00182
00183 ThePEG::EGPtr m_hw;
00184
00185
00186 ThePEG::EventPtr m_event;
00187
00188
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
00194 std::cout<<"Herwig++ initializing..."<<std::endl;
00195 m_hw->initialize();
00196
00197
00198 std::cout<<"Generating events and filling histograms..."<<std::endl;
00199
00200
00201 for (unsigned int iEvent = 0; iEvent < currentConfig->nevents; ++iEvent) {
00202
00203
00204 assert(m_hw);
00205 m_event = m_hw->shoot();
00206
00207
00208
00209
00210
00211 HepMC::GenEvent* evt= ThePEG::HepMCConverter<HepMC::GenEvent>::convert(*m_event, false, ThePEG::GeV);
00212
00213
00214 vector<fastjet::PseudoJet> inclusive_jets = FindJet.GetJet(evt);
00215
00216
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
00232
00233 std::cout<<std::endl;
00234 }
00235
00236
00237 FindJet.ClearEvent(evt);
00238
00239
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
00249 for (std::vector<baseAnalysis*>::const_iterator i(analysis.begin()); i != analysis.end(); ++i) {
00250 delete (*i);
00251 }
00252 analysis.resize(0);
00253
00254
00255 delete currentConfig;
00256
00257
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