JetAnalysis.cc

Go to the documentation of this file.

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

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

#include "fastjet/PseudoJet.hh"
#include "fastjet/ClusterSequence.hh"
#include "fastjet/JetDefinition.hh"
#include "fastjet/SISConePlugin.hh"

// Analysis header
#include "../include/JetAnalysis.h"

using namespace std;

/**
Constructor
*/
JetAnalysis::JetAnalysis()
{
}

/**
Destructor
*/
JetAnalysis::~JetAnalysis()
{
}

/**
Before using the class initialize
*/
int JetAnalysis::Init(double tr_max_eta, double tr_min_pt )
{
  m_max_eta = tr_max_eta;       // default eta cut
  m_min_pt  = tr_min_pt;        // default pt cut
   
  m_outputFileName  = "JetAnalysis.root";    // output file name
  m_outputRootDir   = "Jet";
  
  //declaration of histograms
  m_jet_count=baseAnalysis::initHist(string("jet_count"), string("Nr of Jets"), string("Nr Jets"),16, -0.5, 15.5);
  m_jet_pt=baseAnalysis::initHist(string("jet_pt"),string("Jet Transverse Momentum"),string("Jet p_{T} [GeV]"), 
                                 200, 0.,1000.);
  m_jet_pt_log=baseAnalysis::initHist(string("jet_pt_logscale"),string("Jet Transverse Momentum -logscale-"),
                                             string("Jet p_{T} [GeV]"), 200, 0.,1000.);
  m_leadingjet_pt=baseAnalysis::initHist(string("leadingjet_pt"),string("Leading Jet Transverse Momentum"),string("Jet p_{T} [GeV]"), 
                                         200, 0.,1000.);
  m_leadingjet_pt_log=baseAnalysis::initHist(string("leadingjet_pt_logscale"),string("Leading Jet Transverse Momentum -logscale-"),
                                             string("Jet p_{T} [GeV]"), 200, 0.,1000.);
  m_secondleadingjet_pt=baseAnalysis::initHist(string("secondleadingjet_pt"),string("Second leading Jet Transverse Momentum"),
                                               string("Jet p_{T} [GeV]"), 200, 0.,1000.);
  m_secondleadingjet_pt_log=baseAnalysis::initHist(string("secondleadingjet_pt_logscale"),
                                                   string("Second leading Jet Transverse Momentum -logscale-"),
                                                   string("Jet p_{T} [GeV]"), 200, 0.,1000.);
  
  m_jet_phi=baseAnalysis::initHist(string("jet_phi"),string("Jet Phi angle"),string("Jet #varphi"), 64, 0., 6.4);
  m_jet_eta=baseAnalysis::initHist(string("jet_eta"),string("Jet Pseudo rapidity"),string("Jet #eta"),60, -6., 6.);
  
  m_charged_particle_multiplicity=baseAnalysis::initHist(string("charged_particle_multiplicity"),string("Number of charged Particles in the Event with #eta <2.5 and p_{T}>0.5 GeV"),string("charged particle multiplicity"),50, 0., 500.);
//   m_charged_particle_temp_pt= new TH1D("Jet_charged_particle_mean_pt","temp histogramm",200,0,200.); 
  m_charged_particle_temp_pt=baseAnalysis::initHist(string("Jet_charged_particle_mean_pt"),string("temp histogramm"),string("Jet charged particle mean p_{T}"),200,0,200.); 
  m_charged_particle_pt=baseAnalysis::initHist("charged_particle_pt","p_{T} of charged particles with #eta <2.5 and p_{T}>0.5 GeV","p_{T}",25,0,25.); 
  m_charged_particle_pt_log=baseAnalysis::initHist("charged_particle_pt_logscale","p_{T} of charged particles with #eta <2.5 and p_{T}>0.5 GeV","p_{T}",25,0,25.); 
  m_charged_particle_mean_pt=baseAnalysis::initHist(string("charged_particle_mean_pt"),string("Average p_{T} charged Particles in the Event with #eta <2.5 and p_{T}>0.5 GeV"),string("charged particle #bar{p}_{T}"),50, 0., 10.);
  m_charged_particle_rms_pt=baseAnalysis::initHist(string("charged_particle_rms_pt"),string("RMS p_{T} charged Particles in the Event with #eta <2.5 and p_{T}>0.5 GeV"),string("charged particle #sigma(p_{T})"),50, 0., 10.);
  //m_charged_particle_pdgID=baseAnalysis::initHist(string("charged_particle_pdgID"),string("charged Particles PDG ID"),string("charged particle pid"),400, -400., 400.);
  
  m_jet_Deta_leadingJet_secondJet=baseAnalysis::initHist(string("jet_Deta_leadingJet_secondJet"),string("#Delta#eta leading and second leading Jet"),string("#Delta#eta"), 60, 0., 6.);
  m_jet_Dphi_leadingJet_secondJet=baseAnalysis::initHist(string("jet_Dphi_leadingJet_secondJet"),string("#Delta#varphi leading and second leading Jet"),string("#Delta#varphi"), 64, 0., 3.2);
  m_jet_DR_leadingJet_secondJet=baseAnalysis::initHist(string("jet_DR_leadingJet_secondJet"),string("#DeltaR leading and second leading Jet"),string("#DeltaR"), 60, 0., 15.);
  
  return true;
}

/**
This is the main analysis function. 
Search the particle, get their properties and the histograms.
*/
int JetAnalysis::Process(HepMC::GenEvent* hepmcevt) 
{
  unsigned int chargeParticlecount=0;  
  
  // Reset the histogramm                        
  m_charged_particle_temp_pt->Reset();
  // Fill the eventnumber
  //m_evtnr->Fill((hepmcevt->event_number())%100);
  
  // loop over the particles and select pdgid and pt
  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();
      
      // we want stable charged particles

      // just remove the neutrinos, gammas and neutrons, ... from the
      // final state
      if(!chargedParticle(pid)) continue;
      
      // take all stable particles of this event
     if(!IsFinalStateParticle((*p))) continue;

      m_charged_particle_temp_pt->Fill((*p)->momentum().perp());
      m_charged_particle_pt->Fill((*p)->momentum().perp());
      m_charged_particle_pt_log->Fill((*p)->momentum().perp());
      chargeParticlecount++;
      //m_charged_particle_pdgID->Fill(pid);
  }
  
  // storage of input particles for jet
  CLHEP::HepLorentzVector lv_leadingJet(0,0,0,0);
  
  // check if jets were found
  //if(!m_inclusive_jets.size()) std::cout<<"JetAnalysis::no jets found"<<std::endl;
  
  //Fill the histograms
  m_jet_count->Fill(m_inclusive_jets.size());
  m_charged_particle_multiplicity->Fill(chargeParticlecount);
  m_charged_particle_mean_pt->Fill(m_charged_particle_temp_pt->GetMean());
  m_charged_particle_rms_pt->Fill(m_charged_particle_temp_pt->GetRMS());
  
  // not needed, but to be on the safe side
  m_charged_particle_temp_pt->Reset();
  
  if(m_inclusive_jets.size()){

    for(unsigned int  i=0; i<m_inclusive_jets.size(); i++) {
      m_jet_pt->Fill(m_inclusive_jets[i].perp());
      m_jet_pt_log->Fill(m_inclusive_jets[i].perp());
    }

    m_leadingjet_pt->Fill(m_inclusive_jets[0].perp());
    m_leadingjet_pt_log->Fill(m_inclusive_jets[0].perp());
    m_secondleadingjet_pt->Fill(m_inclusive_jets[1].perp());
    m_secondleadingjet_pt_log->Fill(m_inclusive_jets[1].perp());
    m_jet_phi->Fill(m_inclusive_jets[0].phi());
    m_jet_eta->Fill(m_inclusive_jets[0].eta());
    
    CLHEP::HepLorentzVector lv_firstJet(m_inclusive_jets[0].px(), m_inclusive_jets[0].py(), m_inclusive_jets[0].pz(), m_inclusive_jets[0].e());
    if(m_inclusive_jets.size() > 1 ) {
      CLHEP::HepLorentzVector lv_secondJet(m_inclusive_jets[1].px(), m_inclusive_jets[1].py(), m_inclusive_jets[1].pz(), m_inclusive_jets[1].e());   
      m_jet_DR_leadingJet_secondJet->Fill(fabs(lv_firstJet.deltaR(lv_secondJet)));
      m_jet_Deta_leadingJet_secondJet->Fill(fabs(lv_firstJet.eta()-lv_secondJet.eta()));
      m_jet_Dphi_leadingJet_secondJet->Fill(fabs(lv_firstJet.vect().deltaPhi(lv_secondJet.vect())));
    }
  }
  return true;
}

---