00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017 #include <stdlib.h>
00018 #include <time.h>
00019
00020 #include <stdio.h>
00021 #include <iostream>
00022 #include <fstream>
00023 #include <sstream>
00024 #include "HepMC/PythiaWrapper.h"
00025 #include "HepMC/IO_HEPEVT.h"
00026
00027 #include "HepMC/IO_BaseClass.h"
00028 #include "HepMC/GenEvent.h"
00029 #include "HepMC/SimpleVector.h"
00030 #include "MyPythia6Wrapper.h"
00031
00032 #include "TH1.h"
00033 #include "TFile.h"
00034 #include "TMath.h"
00035 #include <TCanvas.h>
00036
00037 #include "TF1.h"
00038 #include <TH2.h>
00039 #include <TStyle.h>
00040 #include "TLegend.h"
00041 #include <TTree.h>
00042 #include <TString.h>
00043
00044
00045 #include "../../include/Configuration.h"
00046 #include "../../include/JetFinder.h"
00047 #include "../../include/ttbarAnalysis.h"
00048 #include "../../include/ZtautauAnalysis.h"
00049 #include "../../include/WtaunuAnalysis.h"
00050 #include "../../include/bbbarAnalysis.h"
00051 #include "../../include/JetAnalysis.h"
00052 #include "../../include/WplusJetAnalysis.h"
00053 #include "../../include/ZAnalysis.h"
00054 #include "../../include/UEAnalysis.h"
00055 #include "../../include/EtMissAnalysis.h"
00056 #include "../../include/ElasScatAnalysis.h"
00057 #include "../../include/UserAnalysis.h"
00058 #include "readpythiaConfigFile.h"
00059
00060 using namespace std;
00061
00062
00063 int main(int argc, const char* argv[])
00064 {
00065 if(argc <= 1 ) {
00066 cout << "Usage: " << argv[0] << " <Configuration File>" << endl;
00067 return 0;
00068 }
00069
00070
00071 Configuration* currentconfig = new Configuration();
00072 if(currentconfig->read(argv[1]))
00073 {
00074 cout << " ... skip !" << endl;
00075 return 0;
00076 }
00077
00078 JetFinder FindJet;
00079 FindJet.Init(currentconfig->max_eta,currentconfig->min_pt);
00080 FindJet.InitJetFinder( currentconfig->coneRadius, currentconfig->overlapThreshold, currentconfig->jet_ptmin,
00081 currentconfig->lepton_ptmin, currentconfig->DeltaR_lepton_track, currentconfig->jetalgorithm );
00082
00083
00084 vector<baseAnalysis*> analysis;
00085
00086
00087 if (currentconfig->jet) {
00088 cout << "Adding module JetAnalysis" << endl;
00089 analysis.push_back(new JetAnalysis);
00090 }
00091 if (currentconfig->wplusjet) {
00092 cout << "Adding module WplusJetAnalysis" << endl;
00093 analysis.push_back(new WplusJetAnalysis);
00094 }
00095 if(currentconfig->z){
00096 cout << "Adding module ZAnalysis" << endl;
00097 analysis.push_back(new ZAnalysis);
00098 }
00099 if(currentconfig->ztautau){
00100 cout << "Adding module ZtautauAnalysis" << endl;
00101 analysis.push_back(new ZtautauAnalysis);
00102 }
00103 if(currentconfig->wtaunu){
00104 cout << "Adding module WtaunuAnalysis" << endl;
00105 analysis.push_back(new WtaunuAnalysis);
00106 }
00107 if(currentconfig->bbbar){
00108 cout << "Adding module bbbarAnalysis" << endl;
00109 analysis.push_back(new bbbarAnalysis);
00110 }
00111 if(currentconfig->ttbar){
00112 cout << "Adding module ttbarAnalysis" << endl;
00113 analysis.push_back(new ttbarAnalysis);
00114 }
00115 if(currentconfig->ue){
00116 cout << "Adding module UEAnalysis" << endl;
00117 analysis.push_back(new UEAnalysis);
00118 }
00119 if(currentconfig->etmiss){
00120 cout << "Adding module EtMissAnalysis" << endl;
00121 analysis.push_back(new EtMissAnalysis);
00122 }
00123 if(currentconfig->elasScat){
00124 cout << "Adding module ElasScatAnalysis" << endl;
00125 analysis.push_back(new ElasScatAnalysis);
00126 }
00127 if(currentconfig->user){
00128 cout << "Adding module UserAnalysis" << endl;
00129 analysis.push_back(new UserAnalysis);
00130 }
00131
00132
00133
00134 for (vector<baseAnalysis*>::const_iterator i(analysis.begin()); i != analysis.end(); ++i) {
00135 (*i)->Init(currentconfig->max_eta,currentconfig->min_pt);
00136 (*i)->InitJetFinder( currentconfig->coneRadius, currentconfig->overlapThreshold, currentconfig->jet_ptmin,
00137 currentconfig->lepton_ptmin, currentconfig->DeltaR_lepton_track, currentconfig->jetalgorithm );
00138 }
00139
00140
00141
00142
00143
00144 pyjets.n=0;
00145 call_pyhepc(1);
00146
00147 HepMC::HEPEVT_Wrapper::set_max_number_entries(10000);
00148 HepMC::HEPEVT_Wrapper::set_sizeof_real(8);
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159 pydatr.mrpy[1-1] = currentconfig->rseed ;
00160
00161
00162
00163
00164 cout << "Configuration Parameter List:" << endl;
00165
00166 for (unsigned int files=0; files<(currentconfig->ConfigFileNames).size(); files++)
00167 readConfigFile((currentconfig->ConfigFileNames)[files]);
00168
00169
00170
00171
00172
00173
00174
00175
00176 call_pyinit( "CMS", "p", "p", 7000.0 );
00177
00178 call_pystat( 5 );
00179
00180
00181
00182
00183 HepMC::IO_HEPEVT hepevtio;
00184
00185
00186
00187
00188
00189
00190 int nShow = 50;
00191 int nShowPace = max(1,int(currentconfig->nevents)/nShow);
00192
00193 for (unsigned int iEvent = 1 ; iEvent <= currentconfig->nevents; iEvent++) {
00194 if ( iEvent % nShowPace == 0 ) {
00195 cout << "pythia6: event " << iEvent << endl;
00196 }
00197
00198 call_pyevnt();
00199
00200 call_pyhepc(1);
00201 HepMC::GenEvent* evt = hepevtio.read_next_event();
00202
00203 evt->set_event_number(iEvent);
00204 evt->set_signal_process_id(20);
00205
00206
00207
00208
00209 vector<fastjet::PseudoJet> inclusive_jets = FindJet.GetJet(evt);
00210
00211
00212 int ret_all = true;
00213 for (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 cout << "Processing event number " << iEvent << endl;
00221 cout << "throw FAILURE"<<endl;
00222
00223
00224
00225
00226
00227 cout<<endl;
00228 }
00229
00230
00231 FindJet.ClearEvent(evt);
00232
00233
00234 delete evt;
00235 }
00236
00237
00238
00239 call_pystat(1);
00240
00241 TFile* f = new TFile(currentconfig->OutputFileName.c_str(),"RECREATE");
00242 for (vector<baseAnalysis*>::const_iterator i(analysis.begin()); i != analysis.end(); ++i) {
00243 (*i)->Finalize(f);
00244 }
00245
00246
00247 for (vector<baseAnalysis*>::const_iterator i(analysis.begin()); i != analysis.end(); ++i) {
00248 delete (*i);
00249 }
00250 analysis.resize(0);
00251
00252
00253 delete currentconfig;
00254
00255 cout<<"pythia6: program ends"<<endl;
00256 f->Close();
00257
00258 cout<<"Run successfully!"<<endl;
00259 return 0;
00260 }
00261