cascade.cc

Go to the documentation of this file.

/***********************************************************************
*                                                                      *
* File: cascade.cc                                                     *
* This is a simple test program.                                       *
* based on the example cascade program                                 *
* from the MC School at DESY in April 2008                             *
*                                                                      *
* Modified by sj                                                       *
*                                                                      *
***********************************************************************/

//////////////////////////////////////////////////////////////////////////
// Events are read into the HepMC event record from the FORTRAN HEPEVT 
// common block using the IO_HEPEVT strategy -- nothing is done with them.
// This program is just used to find the total time required to transfer
// from HEPEVT into the HepMC event record.
//////////////////////////////////////////////////////////////////////////

#include <iostream>
#include "HepMC/PythiaWrapper.h"
#include "CascadeWrapper.h"
#include "MyCascadeWrapper.h"
#include "HepMC/IO_HEPEVT.h"
//#include "HepMC/IO_Ascii.h"
#include "HepMC/IO_GenEvent.h"
#include "HepMC/GenEvent.h"
#include "PythiaHelper.h"

#include <stdio.h>
#include "HepMC/IO_AsciiParticles.h"
#include "HepMC/SimpleVector.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 "readCascadeConfigFile.h"

extern "C" {
  extern struct {
    int Nevent;
  } steer1_;
}
#define steer1 steer1_

using namespace std;

// ---------------------------------------------------------------------- 
int main( int argc, const char* argv[] ) 
{ 
  if ( argc <= 1 ) {
    cout << "Usage: " << argv[0] << " <Configuration File>" << endl;
    return 0;
  }

  Configuration* currentconfig = new Configuration();
  if ( currentconfig->read( argv[1] ) ) {
    cout << " ... skip !" << 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, 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 and jet finder
  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
  // CASCADE (via PYTHIA6) uses HEPEVT with 10000 entries and 
  // 8-byte floating point numbers. We need to explicitly pass
  // this information to the HEPEVT_Wrapper.
  //
  //    HepMC::HEPEVT_Wrapper::set_sizeof_int(4);
  //    HepMC::HEPEVT_Wrapper::set_sizeof_real(8);
  //    HepMC::HEPEVT_Wrapper::set_max_number_entries(10000);
  //
  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);
  
  //    
  //    initPythia();
  //........................................CASCADE INITIALIZATIONS
  //  call_rluxgo(lux,iseed,k1,k2)
  // lux = luxory level
  // iseed = seed for random number generator
  //    call_rluxgo(4,87777886,0,0);
  //--initialise CASCADE parameters

  call_casini();

  //call steering
  //call_steer();

  //-- change standard parameters of CASCADE
  call_cascha();
  
  //change cascade parameters
  vector< string > cascade_commands;
  // read in commandlines
  for (unsigned int files=0; files<(currentconfig->ConfigFileNames).size(); files++)
    readConfigFile((currentconfig->ConfigFileNames)[files]);

  //-- change standard parameters of JETSET/PYTHIA
  call_pytcha();

  //-- set up for running CASCADE  
  call_cascade();

  //-- print result from integration
  call_caend(1);   

  //
  //........................................HepMC INITIALIZATIONS
  //
  // Instantiate an IO strategy for reading from HEPEVT.
  HepMC::IO_HEPEVT hepevtio;
  //
  // Instantial an IO strategy to write the data to file - it uses the 
  //  same ParticleDataTable
  //HepMC::IO_GenEvent ascii_io("example_MyCASCADE.dat",ios::out);
  //
  //........................................EVENT LOOP
  //int Nevent = steer1.Nevent;
  int Nevent = currentconfig->nevents;

  for ( int iEvent = 1; iEvent <= Nevent; iEvent++ ) {
    if ( iEvent == 1 || iEvent % 100 == 0 ) {
      cout << "Processing Event Number "  << iEvent << endl;
    }
    call_event();      // generate one event with CASCADE

    // pythia pyhepc routine convert common PYJETS in common HEPEVT
    call_pyhepc( 1 );
    HepMC::GenEvent* hepmcevt = hepevtio.read_next_event();
    //
    // add some information to the event
    hepmcevt->set_event_number(iEvent);
    hepmcevt->set_signal_process_id(11);

    //call the Jetfinder
    vector<fastjet::PseudoJet> inclusive_jets = FindJet.GetJet(hepmcevt);

    //call analysis
    int ret_all = true;
    for (vector<baseAnalysis*>::const_iterator i(analysis.begin()); i != analysis.end(); ++i) {
      (*i)->SetJet(&inclusive_jets);
      int ret = (*i)->Process(hepmcevt);
      if (ret == false) ret_all = false;
    }
    
    if (ret_all==false) {
      cout << "Processing event number " << iEvent << endl;    
      cout << "throw FAILURE"<<endl;
      //cout<<endl;
      //cout<<"HepMC Output"<<endl;
      //hepmcevt->print();
      cout<<endl;
    }
    
    //// write the event out to the ascii file
    //ascii_io << hepmcevt;
    
    //delete all the jet objects, missingEt etc. to be free for the next event
    FindJet.ClearEvent(hepmcevt);
    
    // we also need to delete the created event from memory
    delete hepmcevt;

  }
  //........................................TERMINATION
  //  Print out of generated event summary
  call_caend(2);
  // write out some information from Pythia to the screen
  //    call_pystat( 1 );    
  
  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*>::const_iterator i(analysis.begin()); i != analysis.end(); ++i) {
    delete (*i);
  }
  analysis.resize(0);

  // clean up config object
  delete currentconfig;

  cout << "cascade: program ends" << endl;
  f->Close();

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

---