baseAnalysis.cc

Go to the documentation of this file.

#include <iostream>
#include <sstream>
#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"

using namespace std;

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

//top analysis header
#include "../include/baseAnalysis.h"

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

/**
Constructor
*/
baseAnalysis::baseAnalysis()
{ 
m_Jetfinder_enabled=false;
}

/**
Destructor: delete all the Pointer
*/
baseAnalysis::~baseAnalysis()
{
  while (!m_histVector.empty()){
    delete m_histVector.back();
    m_histVector.pop_back();
  }

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

}

/**
InitJetFinder: Initialisation of JetFinder
*/
int baseAnalysis::InitJetFinder(double coneRadius, double overlapThreshold, double jet_ptmin)
{

    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_plugin = new fastjet::SISConePlugin(m_coneRadius, m_overlapThreshold);
  m_jetDef = new fastjet::JetDefinition(m_plugin);

  if(!m_jetDef) return FAILURE;

  m_Jetfinder_enabled=true;
  return SUCCESS;

}

/**
FindJetableParticles
*/
int baseAnalysis::FindJetableParticles(HepMC::GenEvent* hepmcevt)
{

  //delete m_clust_seq;
  //m_clust_seq=0;
  if(!m_Jetfinder_enabled || !m_jetDef){
    std::cout << "ERROR: FastJet Algorithm not initialized!" << std::endl;
    return FAILURE;
  }

  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();
    
    // now collect the tracks for the Jet algorithm
    // first of all we need stable particle
    if (!((*p)->status() == 1)) continue;
    
    // cut out the neutral particles (charge is not stored in HepMC)
    // just remove the neutrinos, gammas and Neutron ... from the final state
    if(IsNeutrino(pid) || IsGamma(pid) || IsNeutron(pid)) 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(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 SUCCESS;
}

/**
FindJet: run JetFinder
*/
int baseAnalysis::FindJet(HepMC::GenEvent* hepmcevt) {
  
  if(!m_Jetfinder_enabled || !m_jetDef){
    std::cout << "ERROR: FastJet Algorithm not initialized!" << std::endl;
    return FAILURE;
  }
  
  m_input_particles.clear();
  m_inclusive_jets.clear();
  m_clust_seq=0;
  
  // run the jet finding
  if(!m_input_particles.size()) {
    int input = baseAnalysis::FindJetableParticles(hepmcevt); 
    if(!input) 
      return FAILURE;
  }
  
  if(m_input_particles.size())
    m_clust_seq = new fastjet::ClusterSequence(m_input_particles,*m_jetDef);
  
  if(!m_input_particles.size() || !m_clust_seq) 
    return FAILURE;
  
  // 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 SUCCESS;
}

/**
DeleteJetObject: delete all the jet objects
*/
int baseAnalysis::DeleteJetObject(HepMC::GenEvent* hepmcevt) {
  m_input_particles.clear();
  m_inclusive_jets.clear();
  delete m_clust_seq;
  m_clust_seq=0;

  return SUCCESS;
}

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

/**
In the final step all the histogramms are stored in a rootfile.
The name of the rootfile can be set with the function setOutpuFileName(const char* filename). 
 */
int baseAnalysis::Finalize(TFile* output) 
{
  if (! output->cd(m_outputRootDir.c_str())) {
    std::cout << "The Directory " << m_outputRootDir.c_str() << " does not exist." << std::endl;
    std::cout << "The Directory will be created" << std::endl;
    TDirectory* 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(); 
  
  for ( vector< TH1D * >::const_iterator iter = m_histVectorVariableBin.begin(); iter!=m_histVectorVariableBin.end(); iter++ ) 
    (*iter)->Write(); 
  
  return SUCCESS;
}

/**
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(string name, string title, string xlabel, int nrBins, double xmin, double xmax)
{
  TH1D * histo=new TH1D(name.c_str(),title.c_str(),nrBins,xmin,xmax);
  histo->GetXaxis()->SetTitle(xlabel.c_str());
  histo->GetYaxis()->SetTitle("Count");
  baseAnalysis::m_histVector.push_back(histo);
  return histo;
}

/**
histograms with variable bin size
*/
TH1D * baseAnalysis::initHistVariableBin(string name, string title, string xlabel, int nbin, Double_t nbinRange[])
{
  TH1D * histoVariableBin=new TH1D(name.c_str(),title.c_str(),nbin,nbinRange);
  histoVariableBin->GetXaxis()->SetTitle(xlabel.c_str());
  histoVariableBin->GetYaxis()->SetTitle("Count");
  baseAnalysis::m_histVectorVariableBin.push_back(histoVariableBin);
  return histoVariableBin;
}

/**
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)
{
  bool found=false;
  int maxLoops=maxGenerations;
  int loop=1;
  
  // initialize the maximal number of mothers (Levels) to be checked
  if(maxGenerations<0) maxLoops=1000;
  
  // 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;
  
  // 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) && (loop<maxLoops))
    {
      current_mother = current_vertex->particles_begin(HepMC::parents);
      
      if((*current_mother)==mother) 
        found=true;
      else
        current_vertex=(*current_mother)->production_vertex();
      
      loop++;
    }
  return found;
}


/** 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) 
{
  bool found=false;
  int maxLoops=maxGenerations;
  int loop=1;
  
  // initialize the maximal number of mothers (Levels) to be checked
  if(maxGenerations<0) maxLoops=1000;
  
  // 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;
  
  // 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) && (loop<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();
      
      loop++;
    }
  return found;
}

---