analysis.cc

Go to the documentation of this file.

/***********************************************************************
*                                                                      *
* File:analysis.cc                                                     *  
* This is a simple test program.                                       *
* Copyright \ufffd 2005 Torbj\ufffdrn Sj\ufffdstrand                   *
*                                                                      *
* Modified by sj and ca, ap                                            *
*                                                                      *
***********************************************************************/

#include <iostream>
#include <stdio.h>
#include "HepMC/IO_GenEvent.h"
#include "HepMC/GenEvent.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>

using namespace std;

#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"

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

int main(int argc, const char* argv[]) {
  
  if(argc <= 2 ) { //ap
    std::cout << "Usage: " << argv[0] << " <Configuration File>" << " <name of file to read from>"<< std::endl; //ap
    return 0;
  }
  
  Configuration* currentconfig = new Configuration();
  if(currentconfig->read(argv[1]))
    {
      std::cout << " ... skip !" << std::endl;
      return 0;
    }
  
  //argv[2] = (Pythia8DATA);
  
  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 );

  
  std::vector<baseAnalysis*> analysis;

  // create analysis Objects and initialize them                                                                                          
  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 (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, currentconfig->jetalgorithm );
  }
  
    
  // Extract settings to be used in the main program.
  //int nEvent = currentconfig->nevents;
  int iEvent =0;
  
  //cout << "before ctor" << endl;

    
  HepMC::IO_GenEvent ascii_in(argv[2],std::ios::in);
  HepMC::GenEvent* hepmcevt = ascii_in.read_next_event();
  while ( hepmcevt) {

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

    //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(hepmcevt);
      if (ret == false) ret_all = false;
    }
    iEvent++;
    if (ret_all==false) {
      std::cout << "Processing event number " << iEvent << std::endl;    
      std::cout << "throw FAILURE"<<std::endl;
      std::cout<<std::endl;
      std::cout<<"HepMC Output"<<std::endl;
      //hepmcevt->print();
      std::cout<<std::endl;
    }

    //delete all the jet objects to be free for the next event
    FindJet.ClearEvent(hepmcevt);
    
    // clear memory from this event
    delete hepmcevt;
    // End of event loop
    ascii_in >> hepmcevt;
  }

  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;

  // Tidy up
  std::cout<<" program ends"<<std::endl;
  f->Close();

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

---