cascade.cc.svn-base

Go to the documentation of this file.
00001 /***********************************************************************
00002 *                                                                      *
00003 * File: cascade.cc                                                     *
00004 * This is a simple test program.                                       *
00005 * based on the example cascade program                                 *
00006 * from the MC School at DESY in April 2008                             *
00007 *                                                                      *
00008 * Modified by sj                                                       *
00009 *                                                                      *
00010 ***********************************************************************/
00011 
00012 //////////////////////////////////////////////////////////////////////////
00013 // Events are read into the HepMC event record from the FORTRAN HEPEVT 
00014 // common block using the IO_HEPEVT strategy -- nothing is done with them.
00015 // This program is just used to find the total time required to transfer
00016 // from HEPEVT into the HepMC event record.
00017 //////////////////////////////////////////////////////////////////////////
00018 
00019 #include <iostream>
00020 #include "HepMC/PythiaWrapper.h"
00021 #include "CascadeWrapper.h"
00022 #include "MyCascadeWrapper.h"
00023 #include "HepMC/IO_HEPEVT.h"
00024 //#include "HepMC/IO_Ascii.h"
00025 #include "HepMC/IO_GenEvent.h"
00026 #include "HepMC/GenEvent.h"
00027 #include "PythiaHelper.h"
00028 
00029 #include <stdio.h>
00030 #include "HepMC/IO_AsciiParticles.h"
00031 #include "HepMC/SimpleVector.h"
00032 
00033 #include "TH1.h"
00034 #include "TFile.h"
00035 #include "TMath.h"
00036 #include <TCanvas.h>
00037 
00038 #include "TF1.h"
00039 #include <TH2.h>
00040 #include <TStyle.h>
00041 #include "TLegend.h"
00042 #include <TTree.h>
00043 #include <TString.h>
00044 
00045 // include our own classes
00046 #include "../../include/Configuration.h"
00047 #include "../../include/JetFinder.h"
00048 #include "../../include/ttbarAnalysis.h"
00049 #include "../../include/ZtautauAnalysis.h"
00050 #include "../../include/WtaunuAnalysis.h"
00051 #include "../../include/bbbarAnalysis.h"
00052 #include "../../include/JetAnalysis.h"
00053 #include "../../include/WplusJetAnalysis.h"
00054 #include "../../include/ZAnalysis.h"
00055 #include "../../include/UEAnalysis.h"
00056 #include "../../include/EtMissAnalysis.h"
00057 #include "../../include/ElasScatAnalysis.h"
00058 #include "../../include/UserAnalysis.h"
00059 #include "readCascadeConfigFile.h"
00060 
00061 extern "C" {
00062   extern struct {
00063     int Nevent;
00064   } steer1_;
00065 }
00066 #define steer1 steer1_
00067 
00068 using namespace std;
00069 
00070 // ---------------------------------------------------------------------- 
00071 int main( int argc, const char* argv[] ) 
00072 { 
00073   if ( argc <= 1 ) {
00074     cout << "Usage: " << argv[0] << " <Configuration File>" << endl;
00075     return 0;
00076   }
00077 
00078   Configuration* currentconfig = new Configuration();
00079   if ( currentconfig->read( argv[1] ) ) {
00080     cout << " ... skip !" << endl;
00081     return 0;
00082   }     
00083 
00084   JetFinder FindJet;
00085   FindJet.Init( currentconfig->max_eta, currentconfig->min_pt );
00086   FindJet.InitJetFinder( currentconfig->coneRadius, currentconfig->overlapThreshold, currentconfig->jet_ptmin,
00087       currentconfig->lepton_ptmin, currentconfig->DeltaR_lepton_track, currentconfig->jetalgorithm );
00088 
00089     // create analysis Objects and initialize them  
00090   vector<baseAnalysis*> analysis;
00091 
00092   if (currentconfig->jet) {
00093     cout << "Adding module JetAnalysis" << endl;
00094     analysis.push_back(new JetAnalysis);
00095   }
00096   if (currentconfig->wplusjet) {
00097     cout << "Adding module WplusJetAnalysis" << endl;
00098     analysis.push_back(new WplusJetAnalysis);
00099   }
00100   if(currentconfig->z){
00101     cout << "Adding module ZAnalysis" << endl;
00102     analysis.push_back(new ZAnalysis);
00103   }
00104   if(currentconfig->ztautau){
00105     cout << "Adding module ZtautauAnalysis" << endl;
00106     analysis.push_back(new ZtautauAnalysis);
00107   }
00108   if(currentconfig->wtaunu){
00109     cout << "Adding module WtaunuAnalysis" << endl;
00110     analysis.push_back(new WtaunuAnalysis);
00111   }
00112   if(currentconfig->bbbar){
00113     cout << "Adding module bbbarAnalysis" << endl;
00114     analysis.push_back(new bbbarAnalysis);
00115   }
00116   if(currentconfig->ttbar ){
00117     cout << "Adding module ttbarAnalysis" << endl;
00118     analysis.push_back(new ttbarAnalysis);
00119   }
00120   if(currentconfig->ue){
00121     cout << "Adding module UEAnalysis" << endl;
00122     analysis.push_back(new UEAnalysis);
00123   }  
00124   if(currentconfig->etmiss){
00125     cout << "Adding module EtMissAnalysis" << endl;
00126     analysis.push_back(new EtMissAnalysis);
00127   }
00128   if(currentconfig->elasScat){
00129     cout << "Adding module ElasScatAnalysis" << endl;
00130     analysis.push_back(new ElasScatAnalysis);
00131   }
00132   if(currentconfig->user){
00133     cout << "Adding module UserAnalysis" << endl;
00134     analysis.push_back(new UserAnalysis);
00135   }
00136 
00137     // init requested analysis modules and jet finder
00138   for ( vector<baseAnalysis*>::const_iterator i( analysis.begin() ); i != analysis.end(); ++i ) {
00139     (*i)->Init( currentconfig->max_eta, currentconfig->min_pt );
00140     (*i)->InitJetFinder( currentconfig->coneRadius, currentconfig->overlapThreshold, currentconfig->jet_ptmin,
00141         currentconfig->lepton_ptmin, currentconfig->DeltaR_lepton_track, currentconfig->jetalgorithm );
00142   }
00143 
00144   //
00145   //........................................HEPEVT
00146   // CASCADE (via PYTHIA6) uses HEPEVT with 10000 entries and 
00147   // 8-byte floating point numbers. We need to explicitly pass
00148   // this information to the HEPEVT_Wrapper.
00149   //
00150   //    HepMC::HEPEVT_Wrapper::set_sizeof_int(4);
00151   //    HepMC::HEPEVT_Wrapper::set_sizeof_real(8);
00152   //    HepMC::HEPEVT_Wrapper::set_max_number_entries(10000);
00153   //
00154   pyjets.n=0;     // ensure dummyness of the next call
00155   call_pyhepc(1); // mstu(8) is set to NMXHEP in this dummy call (version >=6.404)
00156   //HepMC::HEPEVT_Wrapper::set_max_number_entries(pydat1.mstu[8-1]);
00157   HepMC::HEPEVT_Wrapper::set_max_number_entries(10000);
00158   HepMC::HEPEVT_Wrapper::set_sizeof_real(8);
00159   
00160   //    
00161   //    initPythia();
00162   //........................................CASCADE INITIALIZATIONS
00163   //  call_rluxgo(lux,iseed,k1,k2)
00164   // lux = luxory level
00165   // iseed = seed for random number generator
00166   //    call_rluxgo(4,87777886,0,0);
00167   //--initialise CASCADE parameters
00168 
00169   call_casini();
00170 
00171   //call steering
00172   //call_steer();
00173 
00174   //-- change standard parameters of CASCADE
00175   call_cascha();
00176   
00177   //change cascade parameters
00178   vector< string > cascade_commands;
00179   // read in commandlines
00180   for (unsigned int files=0; files<(currentconfig->ConfigFileNames).size(); files++)
00181     readConfigFile((currentconfig->ConfigFileNames)[files]);
00182 
00183   //-- change standard parameters of JETSET/PYTHIA
00184   call_pytcha();
00185 
00186   //-- set up for running CASCADE  
00187   call_cascade();
00188 
00189   //-- print result from integration
00190   call_caend(1);   
00191 
00192   //
00193   //........................................HepMC INITIALIZATIONS
00194   //
00195   // Instantiate an IO strategy for reading from HEPEVT.
00196   HepMC::IO_HEPEVT hepevtio;
00197   //
00198   // Instantial an IO strategy to write the data to file - it uses the 
00199   //  same ParticleDataTable
00200   //HepMC::IO_GenEvent ascii_io("example_MyCASCADE.dat",ios::out);
00201   //
00202   //........................................EVENT LOOP
00203   //int Nevent = steer1.Nevent;
00204   int Nevent = currentconfig->nevents;
00205 
00206   for ( int iEvent = 1; iEvent <= Nevent; iEvent++ ) {
00207     if ( iEvent == 1 || iEvent % 100 == 0 ) {
00208       cout << "Processing Event Number "  << iEvent << endl;
00209     }
00210     call_event();      // generate one event with CASCADE
00211 
00212     // pythia pyhepc routine convert common PYJETS in common HEPEVT
00213     call_pyhepc( 1 );
00214     HepMC::GenEvent* hepmcevt = hepevtio.read_next_event();
00215     //
00216     // add some information to the event
00217     hepmcevt->set_event_number(iEvent);
00218     hepmcevt->set_signal_process_id(11);
00219 
00220     //call the Jetfinder
00221     vector<fastjet::PseudoJet> inclusive_jets = FindJet.GetJet(hepmcevt);
00222 
00223     //call analysis
00224     int ret_all = true;
00225     for (vector<baseAnalysis*>::const_iterator i(analysis.begin()); i != analysis.end(); ++i) {
00226       (*i)->SetJet(&inclusive_jets);
00227       int ret = (*i)->Process(hepmcevt);
00228       if (ret == false) ret_all = false;
00229     }
00230     
00231     if (ret_all==false) {
00232       cout << "Processing event number " << iEvent << endl;    
00233       cout << "throw FAILURE"<<endl;
00234       //cout<<endl;
00235       //cout<<"HepMC Output"<<endl;
00236       //hepmcevt->print();
00237       cout<<endl;
00238     }
00239     
00240     //// write the event out to the ascii file
00241     //ascii_io << hepmcevt;
00242     
00243     //delete all the jet objects, missingEt etc. to be free for the next event
00244     FindJet.ClearEvent(hepmcevt);
00245     
00246     // we also need to delete the created event from memory
00247     delete hepmcevt;
00248 
00249   }
00250   //........................................TERMINATION
00251   //  Print out of generated event summary
00252   call_caend(2);
00253   // write out some information from Pythia to the screen
00254   //    call_pystat( 1 );    
00255   
00256   TFile* f = new TFile(currentconfig->OutputFileName.c_str(),"RECREATE");
00257   for (vector<baseAnalysis*>::const_iterator i(analysis.begin()); i != analysis.end(); ++i) {
00258     (*i)->Finalize(f);
00259   }
00260 
00261   // clean up analysis modules
00262   for (vector<baseAnalysis*>::const_iterator i(analysis.begin()); i != analysis.end(); ++i) {
00263     delete (*i);
00264   }
00265   analysis.resize(0);
00266 
00267   // clean up config object
00268   delete currentconfig;
00269 
00270   cout << "cascade: program ends" << endl;
00271   f->Close();
00272 
00273   cout << "Run successfully!" << endl;
00274   return 0;
00275 }
00276 

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