00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include "Pythia.h"
00011
00012 #include "LHAFortran.h"
00013
00014 #include "HepMCInterface.h"
00015
00016 #include "HepMC/GenEvent.h"
00017
00018
00019 #include <stdlib.h>
00020 #include <time.h>
00021
00022 #include <iostream>
00023 #include <stdio.h>
00024 #include "HepMC/IO_GenEvent.h"
00025 #include "HepMC/GenEvent.h"
00026 #include "HepMC/IO_AsciiParticles.h"
00027 #include "HepMC/SimpleVector.h"
00028
00029 #include "TH1.h"
00030 #include "TFile.h"
00031 #include "TMath.h"
00032 #include <TCanvas.h>
00033
00034 #include "TF1.h"
00035 #include <TH2.h>
00036 #include <TStyle.h>
00037 #include "TLegend.h"
00038 #include <TTree.h>
00039 #include <TString.h>
00040
00041 using namespace Pythia8;
00042 using namespace std;
00043
00044 #include "../../include/Configuration.h"
00045 #include "../../include/JetFinder.h"
00046 #include "../../include/JetAnalysis.h"
00047 #include "../../include/WplusJetAnalysis.h"
00048 #include "../../include/ZAnalysis.h"
00049 #include "../../include/ZtautauAnalysis.h"
00050 #include "../../include/WtaunuAnalysis.h"
00051 #include "../../include/bbbarAnalysis.h"
00052 #include "../../include/ttbarAnalysis.h"
00053 #include "../../include/UEAnalysis.h"
00054 #include "../../include/EtMissAnalysis.h"
00055 #include "../../include/ElasScatAnalysis.h"
00056 #include "../../include/UserAnalysis.h"
00057
00058
00059
00060 int main(int argc, const char* argv[]) {
00061
00062 if(argc <= 1 ) {
00063 std::cout << "Usage: " << argv[0] << " <Configuration File>" << std::endl;
00064 return 0;
00065 }
00066
00067
00068 Configuration* currentconfig = new Configuration();
00069 if(currentconfig->read(argv[1]))
00070 {
00071 std::cout << " ... skip !" << std::endl;
00072 return 0;
00073 }
00074
00075
00076
00077 JetFinder FindJet;
00078 FindJet.Init(currentconfig->max_eta,currentconfig->min_pt);
00079 FindJet.InitJetFinder( currentconfig->coneRadius, currentconfig->overlapThreshold, currentconfig->jet_ptmin,
00080 currentconfig->lepton_ptmin, currentconfig->DeltaR_lepton_track, currentconfig->jetalgorithm );
00081
00082 std::vector<baseAnalysis*> analysis;
00083
00084
00085 if (currentconfig->jet) {
00086 std::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 for (std::vector<baseAnalysis*>::const_iterator i(analysis.begin()); i != analysis.end(); ++i) {
00132 (*i)->Init(currentconfig->max_eta,currentconfig->min_pt);
00133 (*i)->InitJetFinder( currentconfig->coneRadius, currentconfig->overlapThreshold, currentconfig->jet_ptmin,
00134 currentconfig->lepton_ptmin, currentconfig->DeltaR_lepton_track, currentconfig->jetalgorithm );
00135 }
00136
00137
00138 HepMC::I_Pythia8 ToHepMC;
00139
00140
00141
00142 Pythia pythia(std::getenv("Pythia8DATA"));
00143
00144 char processline[64];
00145
00146 for (unsigned int files=0; files<currentconfig->ConfigFileNames.size(); files++)
00147 {
00148 std::cout << "pythia8: Will use settings from file " << currentconfig->ConfigFileNames[files].c_str() << std::endl;
00149 pythia.readFile(currentconfig->ConfigFileNames[files].c_str());
00150 }
00151
00152
00153 pythia.readString("Random:setSeed = on");
00154 sprintf(processline,"Random:seed = %ld",currentconfig->rseed);
00155 pythia.readString(processline);
00156
00157
00158 sprintf(processline,"PartonLevel:FSR = %d",int(currentconfig->FSR));
00159 pythia.readString(processline);
00160
00161 sprintf(processline,"PartonLevel:ISR = %d",int(currentconfig->ISR));
00162 pythia.readString(processline);
00163
00164 sprintf(processline,"PartonLevel:MI = %d",int(currentconfig->MI));
00165 pythia.readString(processline);
00166
00167
00168
00169
00170 pythia.init( 2212, 2212, 14000.);
00171
00172
00173 Event& event = pythia.event;
00174 Settings& settings = pythia.settings;
00175
00176
00177 int nEvent = currentconfig->nevents;
00178
00179
00180 int nShow = settings.mode("Main:timesToShow");
00181 int nAbort = settings.mode("Main:timesAllowErrors");
00182 bool showChangedSettings = settings.flag("Main:showChangedSettings");
00183 bool showAllSettings = settings.flag("Main:showAllSettings");
00184 bool showAllParticleData = settings.flag("Main:showAllParticleData");
00185
00186
00187 if (showChangedSettings) settings.listChanged();
00188 if (showAllSettings) settings.listAll();
00189
00190
00191
00192
00193 if (showAllParticleData) pythia.particleData.listAll();
00194
00195
00196 int nShowPace = max(1,nEvent/nShow);
00197 int iAbort = 0;
00198 for (int iEvent = 0; iEvent < nEvent; ++iEvent) {
00199
00200 if (iEvent%nShowPace == 0)
00201 cout << "pythia8: Now begin event " << iEvent << endl;
00202
00203
00204 if (!pythia.next()) {
00205 if (++iAbort < nAbort) continue;
00206 cout << "pythia8: Event generation aborted prematurely, owing to error!\n";
00207 break;
00208 }
00209
00210
00211
00212
00213
00214
00215
00216
00217 HepMC::GenEvent* hepmcevt = new HepMC::GenEvent();
00218 ToHepMC.fill_next_event( event, hepmcevt );
00219
00220
00221
00222
00223 vector<fastjet::PseudoJet> inclusive_jets = FindJet.GetJet(hepmcevt);
00224
00225
00226 int ret_all = true;
00227 for (std::vector<baseAnalysis*>::const_iterator i(analysis.begin()); i != analysis.end(); ++i) {
00228 (*i)->SetJet(&inclusive_jets);
00229 int ret = (*i)->Process(hepmcevt);
00230 if (ret == false) ret_all = false;
00231 }
00232
00233 if (ret_all==false) {
00234 std::cout << "Processing event number " << iEvent << std::endl;
00235 std::cout << "throw FAILURE"<<std::endl;
00236
00237
00238
00239 std::cout<<std::endl;
00240 }
00241
00242
00243 FindJet.ClearEvent(hepmcevt);
00244
00245
00246 delete hepmcevt;
00247
00248 }
00249
00250
00251 pythia.statistics();
00252
00253 pythia.settings.writeFile("allSettings.cmnd",true);
00254
00255 TFile* f = new TFile(currentconfig->OutputFileName.c_str(),"RECREATE");
00256 for (std::vector<baseAnalysis*>::const_iterator i(analysis.begin()); i != analysis.end(); ++i) {
00257 (*i)->Finalize(f);
00258 }
00259
00260
00261 for (std::vector<baseAnalysis*>::const_iterator i(analysis.begin()); i != analysis.end(); ++i) {
00262 delete (*i);
00263 }
00264 analysis.resize(0);
00265
00266
00267 delete currentconfig;
00268
00269
00270 std::cout<<"pythia8: program ends"<<std::endl;
00271 f->Close();
00272
00273 std::cout<<"Run successfully!"<<std::endl;
00274 return 0;
00275 }
00276