pythia8.cc

Go to the documentation of this file.
00001 /***********************************************************************
00002 *                                                                      *
00003 * File: pythia8_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 #include "HepMC/GenEvent.h"
00017 
00018 // this is needed to generate some random numbers                                                                                                                                                                    
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   //argv[2] = (Pythia8DATA);
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   // create analysis Objects and initialize them                                                                                          
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   // init requested analysis modules and jet finder
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   // Generator. Process selection. LHC initialization.
00141   //Pythia pythia;
00142   Pythia pythia(std::getenv("Pythia8DATA"));
00143 
00144   char processline[64];
00145   // Read in commands from external file.
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   // as done for the random seed, if it is specified
00153   pythia.readString("Random:setSeed = on");
00154   sprintf(processline,"Random:seed = %ld",currentconfig->rseed);
00155   pythia.readString(processline);
00156 
00157   // master switches for FSR, ISR and MI
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   // Settings can also be changed with this command
00168   //pythia.readString("WeakBosonAndParton:all = on");
00169 
00170   pythia.init( 2212, 2212, 14000.);
00171 
00172   // Shorthand for the event and for settings.
00173   Event& event = pythia.event;
00174   Settings& settings = pythia.settings;
00175 
00176   // Extract settings to be used in the main program.
00177   int nEvent = currentconfig->nevents;
00178   //sj//  int nEvent = settings.mode("Main:numberOfEvents");
00179   //int nList = settings.mode("Main:numberToList");
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   // List changed data.
00187   if (showChangedSettings) settings.listChanged();
00188   if (showAllSettings) settings.listAll();
00189 
00190   // List particle data.  
00191   //if (showAllParticleData) ParticleDataTable::listAll();
00192   //starting at least for Pythia8 135 use the following line
00193   if (showAllParticleData) pythia.particleData.listAll();
00194 
00195   // Begin event loop.
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     // Generate events. Quit if too many failures.
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     // List first few events, both hard process and complete events.
00211     //if (iEvent < nList) { 
00212     //if (iEvent < nEvent) { 
00213     //  pythia.process.list();
00214       //event.list();
00215     // }
00216 
00217     HepMC::GenEvent* hepmcevt = new HepMC::GenEvent();
00218     ToHepMC.fill_next_event( event, hepmcevt );
00219 
00220     //hepmcevt->print();
00221 
00222     //call the Jetfinder
00223     vector<fastjet::PseudoJet> inclusive_jets = FindJet.GetJet(hepmcevt);
00224 
00225     //call analysis
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       //std::cout<<std::endl;
00237       //std::cout<<"HepMC Output"<<std::endl;
00238       //hepmcevt->print();
00239       std::cout<<std::endl;
00240     }
00241 
00242     //delete all the jet objects to be free for the next event
00243     FindJet.ClearEvent(hepmcevt);
00244 
00245     // clear memory from this event
00246     delete hepmcevt;
00247     // End of event loop.
00248   }
00249 
00250   // Final statistics.
00251   pythia.statistics();
00252   //  pythia.particleData.listAll();
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   // clean up analysis modules
00261   for (std::vector<baseAnalysis*>::const_iterator i(analysis.begin()); i != analysis.end(); ++i) {
00262     delete (*i);
00263   }
00264   analysis.resize(0);
00265 
00266   // clean up config object
00267   delete currentconfig;
00268 
00269   // Tidy up
00270   std::cout<<"pythia8: program ends"<<std::endl;
00271   f->Close();
00272 
00273   std::cout<<"Run successfully!"<<std::endl;
00274   return 0;
00275 }
00276 

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