baseAnalysis.cc

Go to the documentation of this file.

// STL
#include <iostream>
#include <sstream>
#include <cmath>
#include <stdio.h>

#include "HepMC/GenEvent.h"
#include "HepMC/IO_GenEvent.h"
#include "HepMC/GenParticle.h"
#include "HepMC/GenVertex.h"
#include "HepMC/IO_AsciiParticles.h"
#include "HepMC/SimpleVector.h"
#include "HepPDT/ParticleData.hh"
#include "CLHEP/Vector/LorentzVector.h"
#include "fastjet/ClusterSequence.hh"
#include "fastjet/SISConePlugin.hh"

// ROOT 
#include "TH1.h"
#include "TH2.h"
#include "TFile.h"
#include "TMath.h"
#include "TLorentzVector.h"

#include "../include/baseAnalysis.h"

//#define NDEBUG
#include <cassert>

using namespace std;

/**
Constructor
*/
baseAnalysis::baseAnalysis()
{ 
  // initialize pointers!
  m_plugin = 0;
  m_jetDef = 0;
  m_clust_seq = 0;

  //for jet finder
  m_coneRadius          = 0.;
  m_overlapThreshold    = 0.;
  m_jet_ptmin           = 0.;
  m_lepton_ptmin        = 0.;
  m_DeltaR_lepton_track = 0.;
  m_Jetfinder_enabled   = false;
  
  // specify the maximum eta of tracks
  m_max_eta = 0.;
  m_min_pt  = 0.;

  // for calculation of the missing energy
  exMissTruth     = 0.;
  eyMissTruth     = 0.;
  etMissTruth     = 0.;
  etsumMissTruth  = 0.;

  // another pointer to be initialized  
//  m_particleTable = 0;
}

/**
Destructor: delete all the Pointer
*/
baseAnalysis::~baseAnalysis()
{
  // clean-up memory assigned to pointers
  /*if (m_plugin) {
    delete m_plugin;
    m_plugin = 0;
  }
  
  if (m_jetDef) {
    delete m_jetDef;
    m_jetDef = 0;
  }
  
  if ( m_clust_seq ) {
    m_clust_seq = 0;
    delete m_clust_seq;
  }*/
  
/*  for (std::vector<TH1D*>::const_iterator i(m_histVector.begin()); i != m_histVector.end(); ++i) {
    delete (*i);
  }
  m_histVector.resize(0);*/

  //probably not needed anymore, due to averagedHistograms function
/*  for (std::vector<TH1D*>::const_iterator i(m_histVectorVariableBin.begin()); i != m_histVectorVariableBin.end(); ++i) {
    delete (*i);
  }
  m_histVectorVariableBin.resize(0);*/

//   while (!m_histVectorVariableBin.empty()){
//     delete m_histVectorVariableBin.back();
//     m_histVectorVariableBin.pop_back();
//   }
}

// ---------------------------------------------------------------------- 
int baseAnalysis::Init(double maxeta, double minpt) 
{
  // this method is normally implemented in a derived class(es)
  cout << "baseAnalysis: WARNING: You have called the dummy function Init()." << endl; 
  return 0;
}

// ---------------------------------------------------------------------- 
int baseAnalysis::Process(HepMC::GenEvent* hepmcevt) 
{ 
  // this method is normally implemented in a derived class(es)
  cout << "baseAnalysis: WARNING: You have called the dummy function Process()." << endl; 
  return 0;
} 

// ---------------------------------------------------------------------- 
vector< TH1D* > baseAnalysis::averagedHistograms() 
{
  return m_histVector;
}

// ---------------------------------------------------------------------- 
vector<TH1D*>& baseAnalysis::Histograms()
{
  return m_histVector;
}

/**
InitJetFinder: Initialisation of JetFinder
*/
int baseAnalysis::InitJetFinder(double coneRadius,  double overlapThreshold, double jet_ptmin, double lepton_ptmin, 
  double DeltaR_lepton_track, int jetalgorithm)
{
// jetalgorithm = 1 for siscone, =2 for antikt

    if(m_Jetfinder_enabled){
        // If the Pointer already exist, delete them
        if(m_jetDef) {
            delete m_jetDef;
            m_jetDef=0;
        }
        
        if(m_plugin){
            delete m_plugin;
            m_plugin=0;
        }
    }
    
 // initialise fastjet
  m_coneRadius = coneRadius;
  m_overlapThreshold = overlapThreshold;
  m_jet_ptmin = jet_ptmin;
  m_lepton_ptmin = lepton_ptmin;
  m_DeltaR_lepton_track = DeltaR_lepton_track;

  
  if(jetalgorithm == 1) {
      m_plugin = new fastjet::SISConePlugin(m_coneRadius, m_overlapThreshold);
      m_jetDef = new fastjet::JetDefinition(m_plugin);}
  else if (jetalgorithm==2) {
    m_plugin = 0;
    m_jetDef = new fastjet::JetDefinition(fastjet::antikt_algorithm,m_coneRadius);}
  else
    m_jetDef = 0;

  //   m_plugin = fastjet::antikt_algorithm(m_coneRadius);}


  if(!m_jetDef) return false;

  m_Jetfinder_enabled=true;
  return true;

}



/**
FindJetableParticles
*/
int baseAnalysis::FindJetableParticles(HepMC::GenEvent* hepmcevt)
{
  //delete m_clust_seq;
  //m_clust_seq=0;
  if(!m_Jetfinder_enabled || !m_jetDef){
    cerr << "baseAnalysis: FastJet Algorithm not initialized!" << endl;
    return false;
  }
  CLHEP::HepLorentzVector lv_leadingJet(0,0,0,0);

  for ( HepMC::GenEvent::particle_const_iterator p =
          hepmcevt->particles_begin(); p != hepmcevt->particles_end(); ++p ){

    // use the pdg id to identify the particles
    int pid = (*p)->pdg_id();


    
    // first of all we need stable particle
    if (!(((*p)->status()%1000) == 1)) continue;
    if ((*p)->end_vertex()) continue;
    
    // cut out the neutral particles (charge is not stored in HepMC)
    // just remove the neutrinos and Neutron from the final state
    if(IsNeutrino(pid) || IsNeutron(pid)) continue;

    // cut out all isolated leptons
    // variables for DeltaR calculation
    double deltaPhi = 0;
    double deltaEta = 0;
    double DeltaR = 9999;

    // cut variables
    double DeltaRCut = m_DeltaR_lepton_track;
    double LeptonpTCut = m_lepton_ptmin;

    if(IsElectron(pid) || IsMuon(pid)) {
      int lepton_barcode = (*p)->barcode();
      HepMC::GenParticle *lepton = hepmcevt->barcode_to_particle(lepton_barcode);
      
      // set pt cut for lepton
      if(fabs((*lepton).momentum().perp()) > LeptonpTCut) continue;
      
      // check if it is an isolated lepton
      for ( HepMC::GenEvent::particle_const_iterator q = hepmcevt->particles_begin(); 
            q != hepmcevt->particles_end(); ++q ){
        if(lepton_barcode==(*q)->barcode()) continue;

        // set DeltaR cut
        deltaPhi = (*p)->momentum().phi()-(*q)->momentum().phi();

        if( deltaPhi > TMath::Pi() ||  deltaPhi == TMath::Pi() )
          {
            deltaPhi = deltaPhi - 2 * TMath::Pi() ;
          }
        
        if( deltaPhi < - TMath::Pi() ||  deltaPhi == - TMath::Pi() )
          {
            deltaPhi = deltaPhi + 2 * TMath::Pi();
          }

        deltaEta = (*p)->momentum().eta()-(*q)->momentum().eta();
        double DeltaRnew = sqrt( pow(deltaPhi,2) + pow(deltaEta,2));

        if(DeltaRnew < DeltaR) DeltaR = DeltaRnew;
      }
      if(DeltaR > DeltaRCut) continue;
    }

    // or ... remove all charged particles 
    //if(chargedParticle(pid)==NOTCHARGED) continue;
    
    // They need to have a production Vertex --> we will not select incoming Protons from the Beam
    if (!(*p)->production_vertex()) continue;
    // And since they are stable they should not have a decay vertex.
    if ((*p)->end_vertex()) continue;
    
//     // check wether the track directly comes from a Z-Boson
//     if( trackfromPID(23,(*p)) ) continue;
    
//     // check wether the track directly comes from a W-Boson
//     if( trackfromPID(24,(*p)) ) continue;
//     if( trackfromPID(-24,(*p)) ) continue;
    
    //set the eta cut for charged particles 
    double eta = (*p)->momentum().eta(); 
    if(std::abs(eta) > m_max_eta )  continue;
    
    //set the pt cut for charged particles 
    double pt = (*p)->momentum().perp(); 
    if(pt < m_min_pt) continue;
    
    double px = (*p)->momentum().x();
    double py = (*p)->momentum().y();
    double pz = (*p)->momentum().z();
    double pe = (*p)->momentum().e();
    
    
    // jet finding in stable particles record for charged particles
    fastjet::PseudoJet tempVec;
    tempVec = fastjet::PseudoJet(px,py,pz,pe);
    m_input_particles.push_back(tempVec);
  }
  return true;
}

/**
FindJet: run JetFinder
*/
int baseAnalysis::FindJet(HepMC::GenEvent* hepmcevt) 
{
  if (!m_Jetfinder_enabled || !m_jetDef) {
    cerr << "baseAnalysis: FastJet Algorithm not initialized!" << endl;
    return false;
  }
  
  // run the jet finding
  if (m_input_particles.empty()) {
    if ( !FindJetableParticles(hepmcevt) ) {
      return false;
    }
  }
  
  if (m_input_particles.size()) { 
    m_clust_seq = new fastjet::ClusterSequence(m_input_particles, *m_jetDef);
  } else {
    m_clust_seq = 0;
  }
  
  if (m_input_particles.empty() || m_clust_seq==0) { 
    return false;
  }
  
  // select jets above a threshold, e.g. 0.5 GeV                                                                              
  m_inclusive_jets = sorted_by_pt( m_clust_seq->inclusive_jets( m_jet_ptmin ) );
  
  return true;
}

// ---------------------------------------------------------------------- 
vector< fastjet::PseudoJet > baseAnalysis::GetJet(HepMC::GenEvent* hepmcevt) 
{
  if(m_inclusive_jets.empty()) {
    FindJet(hepmcevt);
    //int inputJet = FindJet(hepmcevt);
    //if(inputJet==FAILURE) std::cout<<"baseAnalysis::no jets found"<<std::endl;
  }
  return m_inclusive_jets;
  //return SUCCESS;
}
  
// ---------------------------------------------------------------------- 
void baseAnalysis::SetJet( vector< fastjet::PseudoJet > *inclusive_jets )
{
  m_inclusive_jets = *inclusive_jets;
}

/**
DeleteJetObject: delete all the jet objects
*/
int baseAnalysis::DeleteJetObject(HepMC::GenEvent* hepmcevt) {
  //if(!m_input_particles.size()&&!m_inclusive_jets.size()&&!m_clust_seq) return FAILURE;
  //if(!(m_input_particles.size()&&m_inclusive_jets.size()&&m_clust_seq)) return FAILURE;

  m_input_particles.clear();
   m_inclusive_jets.clear();
  delete m_clust_seq;
  m_clust_seq=0;

  return true;
}

/**
FindMissingEt: calculate the missing transversal energy of each event
*/
int baseAnalysis::FindMissingEt(HepMC::GenEvent* hepmcevt) 
{
  int properStatus = 1; //stable particles
  
  exMissTruth     = 0.0;
  eyMissTruth     = 0.0;
  etMissTruth     = 0.0;
  etsumMissTruth  = 0.0;
  
  // loop over the particles in GenEvent and select pdgid and pt
  for ( HepMC::GenEvent::particle_const_iterator p =  hepmcevt->particles_begin(); p != hepmcevt->particles_end(); ++p ) {
    int pid = (*p)->pdg_id();

    // skip all not stable particles
    if ( (*p)->status()%1000 != properStatus ) continue;
    if ( (*p)->end_vertex() ) continue;

    // fiducial range eta cut on 2.5
    if(fabs((*p)->momentum().eta()) > 2.5) continue;
    // minimum pt for tracks on 0.5 GeV
    if((*p)->momentum().perp() < 0.5) continue;

    double px  = (*p)->momentum().px();
    double py  = (*p)->momentum().py();
    double pt  = (*p)->momentum().perp();

    if ( IsNeutrino(pid) || abs(pid)==1000039 || abs(pid) ==1000022 || abs(pid) ==39 || abs(pid) ==5100022 ) {
      exMissTruth     += px;    
      eyMissTruth     += py;
      etsumMissTruth  += pt;
    }
  }
  
  etMissTruth = sqrt(exMissTruth*exMissTruth + eyMissTruth*eyMissTruth); 
  
  return true;
}

/**
ClearMissingEt: setting the variables to calculate missing energy to zero
*/
int baseAnalysis::ClearMissingEt(HepMC::GenEvent* hepmcevt) 
{
  exMissTruth     = 0.0;
  eyMissTruth     = 0.0;
  etMissTruth     = 0.0;
  etsumMissTruth  = 0.0;

  return true;
}

/** 
Check if final state particle: returns 0 if it is not a final state particle, 1 otherwise
*/
int baseAnalysis::IsFinalStateParticle(HepMC::GenParticle *p) 
{
  // first of all we need stable particle
  if ( !( (*p).status()%1000 == 1 ) ) return false;

  // And since they are stable they should not have a decay vertex. 
  if ( (*p).end_vertex() ) return false;
    
  // fiducial range eta cut on 2.5
  if( fabs((*p).momentum().eta()) > 2.5 ) return false;
  
  // minimum pt for tracks on 0.5
  if( (*p).momentum().perp() < 0.5 ) return false;

  return true;
}

/**
ClearEvent: 
  delete and clear some variables from the event
*/
int baseAnalysis::ClearEvent(HepMC::GenEvent* hepmcevt) 
{
  baseAnalysis::DeleteJetObject(hepmcevt);
  baseAnalysis::ClearMissingEt(hepmcevt);

  return true;
}

// ---------------------------------------------------------------------- 
int baseAnalysis::setOutputFileName(const string& filename)
{ 
  m_outputFileName = filename; 
  return 0;
}
  
// ---------------------------------------------------------------------- 
int baseAnalysis::setOutputRootDir(const string& dirname)
{ 
  m_outputRootDir = dirname; 
  return 0;
}

// ---------------------------------------------------------------------- 
bool baseAnalysis::IsNeutrino(int pid) 
{ 
  if( abs(pid)==12 || abs(pid)==14 || abs(pid)==16) return true; 
  else return false; 
}

// ---------------------------------------------------------------------- 
bool baseAnalysis::IsGamma(int pid) 
{ 
  if( abs(pid)==22 ) return true; 
  else return false;
}

// ---------------------------------------------------------------------- 
bool baseAnalysis::IsNeutron(int pid) 
{ 
  if( abs(pid)==2112 ) return true; 
  else return false; 
}

// ---------------------------------------------------------------------- 
bool baseAnalysis::IsK0(int pid) 
{ 
  if( abs(pid)==130 || abs(pid)==310 || abs(pid)==311 ) return true; 
  else return false; 
}

// ---------------------------------------------------------------------- 
bool baseAnalysis::IsPi0(int pid) 
{ 
  if( abs(pid)==111 || abs(pid)==113 || abs(pid)==221 ) return true; 
  else return false; 
}

// ---------------------------------------------------------------------- 
bool baseAnalysis::IsElectron(int pid) 
{ 
  if( abs(pid)==11 ) return true; 
  else return false; 
}

// ---------------------------------------------------------------------- 
bool baseAnalysis::IsMuon(int pid) 
{ 
  if( abs(pid)==13 ) return true; 
  else return false; 
}

// ---------------------------------------------------------------------- 
bool baseAnalysis::IsGluon(int pid) 
{ 
  if( abs(pid)==21 ) return true; 
  else return false; 
}

// ---------------------------------------------------------------------- 
bool baseAnalysis::chargedParticle(int pid) 
{ 
  if( IsNeutrino(pid) || IsGamma(pid) || IsNeutron(pid) || IsK0(pid) || IsPi0(pid) || IsGluon(pid) || 
    abs(pid)== 94 || abs(pid)==1000039 || abs(pid) ==1000022 || abs(pid) ==39 || abs(pid) ==5100022 ) {
      //above neutral particles and some "special" neutral particles
    return false;
  }  
  
  return true;
}
  
// ---------------------------------------------------------------------- 
int baseAnalysis::setMaxEtaCut(double maxeta)
{ 
  m_max_eta = maxeta; 
  return 0;
}

// ---------------------------------------------------------------------- 
int baseAnalysis::setMinPtCut(double minpt)
{ 
  m_min_pt = minpt; 
  return 0;
}

/**
Finalize:
In the final step all the histogramms are stored in a rootfile.
The name of the rootfile can be set with the function setOutputFileName(const string& fileName). 
 */
int baseAnalysis::Finalize() 
{
  //write the output in a root file
  TFile f(m_outputFileName.c_str(), "RECREATE");
  Finalize(&f);
  f.Close();
  
  return true;
}

/**
In the final step all the histogramms are stored in a rootfile.
The name of the rootfile can be set with the function setOutputFileName(const string &fileName). 
 */
int baseAnalysis::Finalize(TFile* output) 
{
  TDirectory* dir = output->GetDirectory(m_outputRootDir.c_str());
  if (dir) {
    dir->cd();
  } else {
    cout << "The directory " << m_outputRootDir << " will be created" << endl;
    dir = output->mkdir(m_outputRootDir.c_str());
    dir->cd();
  }

  // loop over the standard vector
  for ( vector< TH1D * >::const_iterator iter = m_histVector.begin(); iter!=m_histVector.end(); iter++ ) 
    (*iter)->Write(); 
  
  // probably not needed anymore, due to averagedHistograms function
  //   for ( vector< TH1D * >::const_iterator iter = m_histVectorVariableBin.begin(); iter!=m_histVectorVariableBin.end(); iter++ ) 
  //     (*iter)->Write(); 
  
  return true;
}

/**
popHisto:
return the last entry  of the histogram vector, needed for ATHENA (ATLAS software) implementation
*/
//TH1D* baseAnalysis::popHisto(std::vector<TH1D*> m_histVector) {
TH1D* baseAnalysis::popHisto() 
{
  // return NULL pointer if no histogram
  if (m_histVector.empty()) return 0;

  // get last element from vector (don't change vector)
  TH1D* temp = m_histVector.back();
  std::cout<<"m_histVector.back()"<< temp;
  // delete last entry from vector
  m_histVector.pop_back();
  return temp;
}

/**
InitHist:
  To get a better handling of the several Histogramms. This class initializes the histograms (name, title, binning, x-axis label) 
  and collects them in the std::vector m_histVector. By doing this looping over all the 
  histograms becomes much more convenient.
*/
TH1D* baseAnalysis::initHist(const string &name, const string &title, const string &xlabel, 
    int nrBins, double xmin, double xmax)
{
  TDirectory* dir = gDirectory->GetDirectory(("/"+m_outputRootDir).c_str());
  if (dir) {
    dir->cd();
  } else {
    gDirectory->cd("/");
    dir = gDirectory->mkdir(m_outputRootDir.c_str());
    dir->cd();
  }

    // convert integer into string
  int histCnt = m_histCounter[ m_outputRootDir ];
  ostringstream ssHistCnt;
  ssHistCnt << histCnt; 

  TH1D* histo = new TH1D((m_outputRootDir + "_" + ssHistCnt.str() + "_" + name).c_str(), title.c_str(), 
    nrBins, xmin, xmax);
  //TH1D* histo=new TH1D((m_outputRootDir+"_"+intToString(m_histCounter[m_outputRootDir]) + "_" + name).c_str(),
  //title.c_str(),nrBins,xmin,xmax);

  histo->GetXaxis()->SetTitle(xlabel.c_str());
  histo->GetYaxis()->SetTitle("Count");
  m_histVector.push_back(histo);

  ++m_histCounter[m_outputRootDir];

  return histo;
}

/**
InitHistVariableBin: 
  histograms with variable bin size
*/
TH1D* baseAnalysis::initHistVariableBin(const string &name, const string &title, const string &xlabel, 
    int nbin, Double_t nbinRange[])
{
  TDirectory* dir = gDirectory->GetDirectory( ("/"+m_outputRootDir).c_str() );
  if (dir) {
    dir->cd();
  } else {
    gDirectory->cd("/");
    dir = gDirectory->mkdir(m_outputRootDir.c_str());
    dir->cd();
  }

  TH1D* histoVariableBin = new TH1D(name.c_str(), title.c_str(), nbin, nbinRange);
  histoVariableBin->GetXaxis()->SetTitle(xlabel.c_str());
  histoVariableBin->GetYaxis()->SetTitle("Count");
  m_histVector.push_back(histoVariableBin);

  return histoVariableBin;
}

/**
CheckDaughter:
  This function returns true if the particle daughter is the decay product of the particle mother. 
  The function scans the mother particle up to a certain level,
  which can be specified with maxGenerations (negativ number --> all levels).
*/
bool baseAnalysis::checkDaughter(HepMC::GenParticle* mother, HepMC::GenParticle* daughter, int maxGenerations)
{
  // initialize the maximal number of mothers (Levels) to be checked
  const size_t maxLoops = ((maxGenerations < 0) ? 1000 : maxGenerations);
  
  // initialize current Production Vertex and mother particles
  // There should always be only one mother particle for decaying particles
  HepMC::GenVertex* current_vertex = daughter->production_vertex();
  HepMC::GenVertex::particle_iterator current_mother;
  
  bool found = false;
  size_t loopCnt = 1;

  // iterate through the mother particles and compare with mother
  // If found match, than stop the loop, otherwise try until there is no production vertex
  // or until you have reached the max number of loops specified
  while ((current_vertex) && !(found) && (loopCnt < maxLoops)) {
    current_mother = current_vertex->particles_begin(HepMC::parents);
      
    if ((*current_mother) == mother) {
        found = true;
    } else {
        current_vertex = (*current_mother)->production_vertex();
    }
    // check next level  
    ++loopCnt;
  }
  return found;
}

/** 
TrackFromPID: 
  check if the track comes from a specific particle e.g. pid(W-Boson)=23 
*/
bool baseAnalysis::trackfromPID(int pid, HepMC::GenParticle* track, int maxGenerations) 
{
  // initialize the maximal number of mothers (Levels) to be checked
  const size_t maxLoops = ((maxGenerations < 0) ? 1000 : maxGenerations);
   
  // initialize current Production Vertex and mother particles
  // There should always be only one mother particle for decaying particles
  HepMC::GenVertex* current_vertex=track->production_vertex();
  HepMC::GenVertex::particle_iterator current_mother;
  
  bool found = false;
  size_t loopCnt = 1;

  // iterate through the mother particles and compare with mother
  // If found match, than stop the loop, otherwise try until there is no production vertex
  // or until you have reached the max number of loops specified
  while ((current_vertex) && !(found) && (loopCnt < maxLoops)) {
    current_mother = current_vertex->particles_begin(HepMC::parents);
      
    //if no mother anymore, jump out the loop
    if (*current_mother == 0) break;     

    if ((*current_mother)->pdg_id() == pid) {
            found = true;
    } else {
        current_vertex = (*current_mother)->production_vertex();
    }
    // check next level  
    ++loopCnt;
  }

  return found;
}

/**
GetRapidity: 
  calculate rapidity of a given particle 
*/
double baseAnalysis::getRapidity(const HepMC::GenEvent::particle_const_iterator &p) 
{
  double e  = (*p)->momentum().e();
  double pz = (*p)->momentum().pz();
  double rapidity = 0.5 * log((e + pz) / (e - pz));
  return rapidity;
}

---