pythia6.cc

Go to the documentation of this file.

/***********************************************************************
*                                                                      *
* File: pythia6_jet.cc                                                 *
* This is a simple test program.                                       *
* Copyright \ufffd 2005 Torbj\ufffdrn Sj\ufffdstrand                   *
*                                                                      *
* Modified by sj and ca                                                *
*                                                                      *
***********************************************************************/
//#include "Pythia.h"

//#include "LHAFortran.h"

//#include "HepMCInterface.h"

// this is needed to generate some random numbers
#include <stdlib.h>
#include <time.h>

#include <stdio.h>
#include <iostream>
#include <fstream>
#include <sstream>
#include "HepMC/PythiaWrapper.h"
#include "HepMC/IO_HEPEVT.h"
#include "HepMC/IO_Ascii.h"
#include "HepMC/GenEvent.h"
#include "HepMC/SimpleVector.h"
#include "MyPythia6Wrapper.h"

#include "TH1.h"
#include "TFile.h"
#include "TMath.h"
#include <TCanvas.h>

#include "TF1.h"
#include <TH2.h>
#include <TStyle.h>
#include "TLegend.h"
#include <TTree.h>
#include <TString.h>

using namespace std;

// include our own classes
#include "../../include/Configuration.h"
#include "../../include/JetFinder.h"
#include "../../include/TopAnalysis.h"
#include "../../include/TauAnalysis.h"
#include "../../include/DiJetAnalysis.h"
#include "../../include/WplusJetAnalysis.h"
#include "../../include/ZAnalysis.h"
#include "../../include/UEAnalysis.h"
#include "../../include/EtMissAnalysis.h"
#include "../../include/ElasScatAnalysis.h"
#include "../../include/UserAnalysis.h"
#include "readpythiaConfigFile.h"

/**************************************************************************/

int main(int argc, const char* argv[]) {
    
  if(argc <= 1 ) {
    std::cout << "Usage: " << argv[0] << " <Configuration File>" << std::endl;
    return 0;
  }
  
  
  Configuration* currentconfig = new Configuration();
  if(currentconfig->read(argv[1]))
    {
      std::cout << " ... skip !" << std::endl;
      return 0;
    }
  
  JetFinder FindJet;
  FindJet.Init(currentconfig->max_eta,currentconfig->min_pt);
  FindJet.InitJetFinder(currentconfig->coneRadius,currentconfig->overlapThreshold,currentconfig->jet_ptmin,currentconfig->lepton_ptmin,currentconfig->DeltaR_lepton_track);


  std::vector<baseAnalysis*> analysis;

  // create analysis Objects and initialize them                                                                                          
  if (currentconfig->dijet) {
    std::cout << "Adding module DiJetAnalysis" << std::endl;
    analysis.push_back(new DiJetAnalysis);
  }
  if (currentconfig->wplusjet) {
    std::cout << "Adding module WplusJetAnalysis" << std::endl;
    analysis.push_back(new WplusJetAnalysis);
  }
  if(currentconfig->z){
    std::cout << "Adding module ZAnalysis" << std::endl;
    analysis.push_back(new ZAnalysis);
  }
  if(currentconfig->tau){
    std::cout << "Adding module TauAnalysis" << std::endl;
    analysis.push_back(new TauAnalysis);
  }
  if(currentconfig->top){
    std::cout << "Adding module TopAnalysis" << std::endl;
    analysis.push_back(new TopAnalysis);
  }
  if(currentconfig->ue){
    std::cout << "Adding module UEAnalysis" << std::endl;
    analysis.push_back(new UEAnalysis);
  }  
  if(currentconfig->etmiss){
    std::cout << "Adding module EtMissAnalysis" << std::endl;
    analysis.push_back(new EtMissAnalysis);
  }
  if(currentconfig->elasScat){
    std::cout << "Adding module ElasScatAnalysis" << std::endl;
    analysis.push_back(new ElasScatAnalysis);
  }
  if(currentconfig->user){
    std::cout << "Adding module UserAnalysis" << std::endl;
    analysis.push_back(new UserAnalysis);
  }
  

  // init requested analysis modules and jet finder
  for (std::vector<baseAnalysis*>::const_iterator i(analysis.begin()); i != analysis.end(); ++i) {
    (*i)->Init(currentconfig->max_eta,currentconfig->min_pt);
    (*i)->InitJetFinder(currentconfig->coneRadius,currentconfig->overlapThreshold,currentconfig->jet_ptmin,currentconfig->lepton_ptmin,currentconfig->DeltaR_lepton_track);
  }

  //........................................HEPEVT
  //  Standard Pythia >6.1 uses HEPEVT with NMXHEP=4000 entries and 8-byte
  //  floating point numbers. We need to explicitly pass this information
  //  to the HEPEVT_Wrapper.
  pyjets.n=0; // ensure dummyness of the next call
  call_pyhepc(1); // mstu(8) is set to NMXHEP in this dummy call (version >=6.404)
  //HepMC::HEPEVT_Wrapper::set_max_number_entries(pydat1.mstu[8-1]);
  HepMC::HEPEVT_Wrapper::set_max_number_entries(10000);
  HepMC::HEPEVT_Wrapper::set_sizeof_real(8);
  
  //
  //........................................PYTHIA INITIALIZATIONS
  // (Some platforms may require the initialization of pythia PYDATA block
  //  data as external - if you get pythia initialization errors try
  //  commenting in/out the below call to initpydata() )
  // initpydata();
  //
  
  // set random number seed (mandatory!) ... this may be overwritten by function readConfigFile
  pydatr.mrpy[1-1] = currentconfig->rseed ;
  
  //int tunemode=300;
  //pytune( &tunemode );
  
  std::cout << "Configuration Parameter List:" << std::endl;
  
  for (unsigned int files=0; files<(currentconfig->ConfigFileNames).size(); files++)
    readConfigFile((currentconfig->ConfigFileNames)[files]);

  // these three lines are set in the config file Pythia6_Common.config
  //   // master switches for FSR, ISR and MI
  //   pypars.mstp[61-1]=int(currentconfig->ISR); 
  //   pypars.mstp[71-1]=int(currentconfig->FSR); 
  //   pypars.mstp[81-1]=int(currentconfig->MI); 

  // Call pythia initialization
  call_pyinit("CMS", "p", "p", 14000.);
  
  //........................................HepMC INITIALIZATIONS
  //
  // Instantiate an IO strategy for reading from HEPEVT.
  HepMC::IO_HEPEVT hepevtio;
  //
  // Instantiate an IO strategy to write the data to file - it uses the 
  //  same ParticleDataTable
  //HepMC::IO_Ascii ascii_io("test_pythia_hepmc.dat",std::ios::out);
  
  //.......EVENT LOOP:
  
  int nShow = 50;
  int nShowPace = max(1,int(currentconfig->nevents)/nShow);
  
  for (unsigned int iEvent = 1 ; iEvent <= currentconfig->nevents; iEvent++) {
    
    // Begin event loop. 
    if (iEvent%nShowPace == 0) cout << "pythia6:  Now begin event "
                                    << iEvent << endl;
    
    call_pyevnt();      // generate one event with Pythia
    // pythia pyhepc routine converts common PYJETS in common HEPEVT
    call_pyhepc(1);
    HepMC::GenEvent* evt = hepevtio.read_next_event();
    // add some information to the event
    evt->set_event_number(iEvent);
    evt->set_signal_process_id(20);
    // write the event out to the ascii file
    //    ascii_io << evt;
    
    //call the Jetfinder
    vector<fastjet::PseudoJet> inclusive_jets = FindJet.GetJet(evt);

    //call analysis
    int ret_all = true;
    for (std::vector<baseAnalysis*>::const_iterator i(analysis.begin()); i != analysis.end(); ++i) {
      (*i)->SetJet(&inclusive_jets);
      int ret = (*i)->Process(evt);
      if (ret == false) ret_all = false;
    }

    if (ret_all==false) {
      std::cout << "Processing event number " << iEvent << std::endl;    
      std::cout << "throw FAILURE"<<std::endl;
      //std::cout<<"pylist"<<std::endl;
      //call_pylist(2);
      //std::cout<<std::endl;
      //std::cout<<"HepMC Output"<<std::endl;
      //evt->print();
      std::cout<<std::endl;
    }

    //delete all the jet objects, missingEt etc. to be free for the next event
    FindJet.ClearEvent(evt);
    
    // we also need to delete the created event from memory
    delete evt;
  } // end of event loop
  
  
  // write out some information from Pythia to the screen
  call_pystat(1);    
  
  TFile* f = new TFile(currentconfig->OutputFileName.c_str(),"RECREATE");
  for (std::vector<baseAnalysis*>::const_iterator i(analysis.begin()); i != analysis.end(); ++i) {
    (*i)->Finalize(f);
  }

  // clean up analysis modules
  for (std::vector<baseAnalysis*>::const_iterator i(analysis.begin()); i != analysis.end(); ++i) {
    delete (*i);
  }
  analysis.resize(0);

  // clean up config object
  delete currentconfig;
  
  std::cout<<"pythia6:  program ends"<<std::endl;
  f->Close();

  std::cout<<"Run successfully!"<<std::endl;
  return 0;
  
}

---