alpgen_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"
#include <stdlib.h>       // this is needed to generate some random numbers
#include <time.h>
#include <stdio.h>

#include <iostream>
#include <vector>
#include <fstream>
#include <sstream>

#include "HepMC/PythiaWrapper.h"
#include "HepMC/IO_HEPEVT.h"
// #include "HepMC/IO_Ascii.h"
#include "HepMC/IO_BaseClass.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>

// include our own classes
#include "../../include/Configuration.h"
#include "../../include/JetFinder.h"
#include "../../include/ttbarAnalysis.h"
#include "../../include/ZtautauAnalysis.h"
#include "../../include/WtaunuAnalysis.h"
#include "../../include/bbbarAnalysis.h"
#include "../../include/JetAnalysis.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"
#include "alpgen_spec.h"
#include "test_hepmc.h"

using namespace std;

extern "C" {

  void pyupev_();
};

// ---------------------------------------------------------------------- 
int main( int argc, const char *argv[] ) 
{
  if ( argc <= 1 ) {
    cout << "Usage: " << argv[ 0 ] << " <Configuration File>" << endl;
    return 0;
  }
 
    // read configuration 
  Configuration *currentconfig = new Configuration;
  if ( currentconfig->read( argv[ 1 ] ) ) {
    cout << " ... skip !" << endl;
    return 0;
  }
  
    // initialize Jet Finder
  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, currentconfig->jetalgorithm );

    // create analysis objects and initialize them 
  vector< baseAnalysis* > analysis;

  if ( currentconfig->jet ) {
    cout << "Adding module JetAnalysis" << endl;
    analysis.push_back( new JetAnalysis );
  }
  if ( currentconfig->wplusjet ) {
    cout << "Adding module WplusJetAnalysis" << endl;
    analysis.push_back( new WplusJetAnalysis );
  }
  if ( currentconfig->z ) {
    cout << "Adding module ZAnalysis" << endl;
    analysis.push_back( new ZAnalysis );
  }
  if ( currentconfig->ztautau ) {
    cout << "Adding module ZtautauAnalysis" << endl;
    analysis.push_back( new ZtautauAnalysis );
  }
  if ( currentconfig->wtaunu ) {
    cout << "Adding module WtaunuAnalysis" << endl;
    analysis.push_back( new WtaunuAnalysis );
  }
  if ( currentconfig->bbbar ) {
    cout << "Adding module bbbarAnalysis" << endl;
    analysis.push_back( new bbbarAnalysis );
  }
  if ( currentconfig->ttbar ) {
    cout << "Adding module ttbarAnalysis" << endl;
    analysis.push_back( new ttbarAnalysis );
  }
  if ( currentconfig->ue ) {
    cout << "Adding module UEAnalysis" << endl;
    analysis.push_back( new UEAnalysis );
  }  
  if( currentconfig->etmiss ) {
    cout << "Adding module EtMissAnalysis" << endl;
    analysis.push_back( new EtMissAnalysis );
  }
  if ( currentconfig->elasScat ) {
    cout << "Adding module ElasScatAnalysis" << endl;
    analysis.push_back( new ElasScatAnalysis );
  }
  if ( currentconfig->user ) {
    cout << "Adding module UserAnalysis" << endl;
    analysis.push_back( new UserAnalysis );
  }

    // init requested analysis modules 
  for ( 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, currentconfig->jetalgorithm );
  }

    //........................................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( 4000 );
  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 );
  
  cout << "Configuration Parameter List:" << endl;
  
  for ( size_t iFile = 0; iFile < currentconfig->ConfigFileNames.size(); ++iFile ) {
    readConfigFile( currentconfig->ConfigFileNames.at( iFile ) );
  }

  // 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); 

    // initialize Pythia
  //call_pyinit( "CMS", "p", "p", 14000. );
  call_pyinit( "USER", "", "", 0. );
  //upinit();       // user process
   
  evstat_.nvpass = 0; 
  evstat_.nvcall = 0; 
  evstat_.nvfail = 0;
  
    //........................................HepMC INITIALIZATIONS
    //
    // Instantiate an IO strategy for reading from HEPEVT.
  HepMC::IO_HEPEVT hepevtio;
  HepMC::GenEvent *evt = 0;
  
  //
  // 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);

    // a bit of setup
  //int pyhepcMode = 1;         // PYJETS -> PYHEPC transformation argument
  int pylistMode = 1;         // pylist call argument, see Pythia manual for available options
  int pystatMode = 1;         // pystat call argument, see Pythia manual for available options
    
  bool doFillEvWeight = true;
  bool doFillEvXsec   = true;

    // initialize user process common block
  pyupev_();
 
    //.......EVENT LOOP:
  int nShow = 50;
  int nShowPace = max( 1, static_cast< int >( currentconfig->nevents ) / nShow );
  
  HepMC::IO_GenEvent ascii_io( "myexample_alpgen.dat", std::ios::out );

  for ( size_t iEvent = 1; iEvent <= currentconfig->nevents; ++iEvent ) {
    if ( iEvent % nShowPace == 0) {
      cout << "pythia6:  Now begin event " << iEvent << endl;
    }

      // generate next event with Pythia
    call_pyevnt();      
    //upevnt();
    call_pylist( pylistMode );
    
      // pythia pyhepc routine converts common PYJETS in common HEPEVT
      // get the HEPEVT record of the event
    call_pyhepc( 1 ); //pyhepcMode );
    hepevtio.print();
   
      //ALSFIN: NUP=0 in the event means the end of file
    if ( hepeup_.NUP == 0 ) {
      cout << "end of ER after " << iEvent - 1 << " events." << endl;
      break;
    }
            
      // get the HEPMC record of the event
    evt = hepevtio.read_next_event();
    evt->print();

    if ( !( evt->is_valid() ) ) {
      cerr << "got invalid HepMC event after " << iEvent << " events." << endl;
      cerr << "exit the generation loop now. " << endl;
      break;
    } 
    
    if ( !testHepMC( evt ) ) {
      cerr << "ERROR: Invalid HepMC event! Event #: " << iEvent << endl; 
      cerr << ">> Exit the generation loop now." << endl;
      //break;      !!!
    }


    //hepevtio.fill_next_event();
    
      // add some information to the event ???????????????
    evt->set_event_number( iEvent );
    evt->set_signal_process_id( 20 );

    ascii_io << evt;

      // fill alpgen-specific event weight info
    if ( doFillEvWeight ) {
      fill_eweight( evt );
    }
    
      // fill x-section info
    if ( doFillEvXsec ) {
      fill_exsec( evt );
    }


      // some printout for the first event
/*    if ( iEvent == 1 ) {
      cout << "----------------------------------------------------------------------" << endl;
      cout << " Documentation printout from alpgen_pythia6.cc: " << endl;
      cout << " Listing the Pythia event # " << iEvent << " in Pythia internal PYJETS format: " << endl;
      call_pylist( pylistMode );
        //std::cout<<" Listing the Pythia event # "<<n_events<<" in HepMC format:   "<<std::endl;
      //evt->print();
      cout << "----------------------------------------------------------------------" << endl;
    }*/
                
    
    // write the event out to the ascii file
    //    ascii_io << evt;
    
      // call the Jetfinder
    vector< fastjet::PseudoJet > inclusiveJets = FindJet.GetJet( evt );

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

    if ( !retAll ) {
      cerr << "Processing event number " << iEvent << endl;    
      cerr << "throw FAILURE" << endl;
      cout << ">>> pylist output: " << endl;
      call_pylist( 2 );
      //cout << ">>> HepMC Output" << endl;
      //evt->print();
      cerr << endl;
    }

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

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

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

  cout << "Successfull run!" << endl;
  return 0;
}

---