pythia6.cc

Go to the documentation of this file.
00001 /***********************************************************************
00002 *                                                                      *
00003 * File: pythia6_jet.cc                                                 *
00004 * This is a simple test program.                                       *
00005 * Copyright \ufffd 2005 Torbj\ufffdrn Sj\ufffdstrand                   *
00006 *                                                                      *
00007 * Modified by sj and ca                                                *
00008 *                                                                      *
00009 ***********************************************************************/
00010 //#include "Pythia.h"
00011 
00012 //#include "LHAFortran.h"
00013 
00014 //#include "HepMCInterface.h"
00015 
00016 // this is needed to generate some random numbers
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 // #include "HepMC/IO_Ascii.h"
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 // include our own classes
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   // create analysis Objects and initialize them                                                                                          
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   // init requested analysis modules and jet finder
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   //........................................HEPEVT
00141   //  Standard Pythia >6.1 uses HEPEVT with NMXHEP=4000 entries and 8-byte
00142   //  floating point numbers. We need to explicitly pass this information
00143   //  to the HEPEVT_Wrapper.
00144   pyjets.n=0; // ensure dummyness of the next call
00145   call_pyhepc(1); // mstu(8) is set to NMXHEP in this dummy call (version >=6.404)
00146   //HepMC::HEPEVT_Wrapper::set_max_number_entries(pydat1.mstu[8-1]);
00147   HepMC::HEPEVT_Wrapper::set_max_number_entries(10000);
00148   HepMC::HEPEVT_Wrapper::set_sizeof_real(8);
00149   
00150   //
00151   //........................................PYTHIA INITIALIZATIONS
00152   // (Some platforms may require the initialization of pythia PYDATA block
00153   //  data as external - if you get pythia initialization errors try
00154   //  commenting in/out the below call to initpydata() )
00155   // initpydata();
00156   //
00157   
00158   // set random number seed (mandatory!) ... this may be overwritten by function readConfigFile
00159   pydatr.mrpy[1-1] = currentconfig->rseed ;
00160   
00161   //int tunemode=300;
00162   //pytune( &tunemode );
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   // these three lines are set in the config file Pythia6_Common.config
00170   //   // master switches for FSR, ISR and MI
00171   //   pypars.mstp[61-1]=int(currentconfig->ISR); 
00172   //   pypars.mstp[71-1]=int(currentconfig->FSR); 
00173   //   pypars.mstp[81-1]=int(currentconfig->MI); 
00174 
00175   // Call pythia initialization
00176   call_pyinit( "CMS", "p", "p", 7000.0 );
00177 
00178   call_pystat( 5 );
00179   
00180   //........................................HepMC INITIALIZATIONS
00181   //
00182   // Instantiate an IO strategy for reading from HEPEVT.
00183   HepMC::IO_HEPEVT hepevtio;
00184   //
00185   // Instantiate an IO strategy to write the data to file - it uses the 
00186   //  same ParticleDataTable
00187   //HepMC::IO_Ascii ascii_io("test_pythia_hepmc.dat",ios::out);
00188   
00189   //.......EVENT LOOP:
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();      // generate one event with Pythia
00199     // pythia pyhepc routine converts common PYJETS in common HEPEVT
00200     call_pyhepc(1);
00201     HepMC::GenEvent* evt = hepevtio.read_next_event();
00202     // add some information to the event
00203     evt->set_event_number(iEvent);
00204     evt->set_signal_process_id(20);
00205     // write the event out to the ascii file
00206     //    ascii_io << evt;
00207     
00208     //call the Jetfinder
00209     vector<fastjet::PseudoJet> inclusive_jets = FindJet.GetJet(evt);
00210 
00211     //call analysis
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       //cout<<"pylist"<<endl;
00223       //call_pylist(2);
00224       //cout<<endl;
00225       //cout<<"HepMC Output"<<endl;
00226       //evt->print();
00227       cout<<endl;
00228     }
00229 
00230     //delete all the jet objects, missingEt etc. to be free for the next event
00231     FindJet.ClearEvent(evt);
00232     
00233     // we also need to delete the created event from memory
00234     delete evt;
00235   } // end of event loop
00236   
00237   
00238   // write out some information from Pythia to the screen
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   // clean up analysis modules
00247   for (vector<baseAnalysis*>::const_iterator i(analysis.begin()); i != analysis.end(); ++i) {
00248     delete (*i);
00249   }
00250   analysis.resize(0);
00251 
00252   // clean up config object
00253   delete currentconfig;
00254   
00255   cout<<"pythia6:  program ends"<<endl;
00256   f->Close();
00257 
00258   cout<<"Run successfully!"<<endl;
00259   return 0;
00260 }
00261 

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