reader.cc

Go to the documentation of this file.

/**
 *
 * @short Program to read HepMC format from file and run analysis modules
 *
 * @author Wolfgang Ehrenfeld, DESY, Germany
 *
 * @version \$Id: reader.cc,v 1.3 2009-05-27 11:06:59 sjohnert Exp $
 *
 */

#include "HepMC/IO_GenEvent.h"
#include "HepMC/GenEvent.h"

#include <iostream>
#include <fstream>

// include our analysis modules
#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 "TFile.h"

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


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


  // Loop over events from files and fill histograms
  cout<<"Reading events and filling histograms..."<<endl;
  
  size_t iEvent(0);
  for (size_t i(0); iEvent < currentconfig->nevents && i<currentconfig->InputFileNames.size(); ++i) {

    // open input stream
    ifstream istr(currentconfig->InputFileNames[i].c_str());
    if (!istr) {
      cerr << "reader: can not open file "  << currentconfig->InputFileNames[i] << endl;
      exit(-1);
    }

    HepMC::IO_GenEvent ascii_in(istr);
    
    // now read the file
    HepMC::GenEvent* evt = ascii_in.read_next_event();

    while (evt && (iEvent < currentconfig->nevents)) {
      ++iEvent;

      /* --- analysis of event --- */

      // debug print out
      if (iEvent%50==0) {
        cout << "Processing event number " << iEvent << endl;
        // evt->print();
      }

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

      // call active analysis modules
      int ret_all = true;
      for (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;
      }
    
      // catch errors
      if (ret_all==false) {
        cout << "Processing event number " << iEvent << endl;
        cout << "throw FAILURE"<<endl;
        //cout<<"HepMC Output"<<endl;
        //evt->print();
        cout<<endl;
      }
    
      // delete all the jet objects, missingEt etc. to be free for the next event
      FindJet.ClearEvent(evt);
          
      //delete event and read in next
      delete evt;
      ascii_in >> evt;

    } // while loop over events

  } // loop over input fules
  cout << "reader: " << iEvent << " events found." << endl;

  // finalise step
  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 and configuration
  for (vector<baseAnalysis*>::const_iterator i(analysis.begin()); i != analysis.end(); ++i) {
    delete (*i);
  }
  analysis.resize(0);

  delete currentconfig;

  cout << "reader: finished successfully!" << endl;
  return 0;
}

---