DiJetAnalysis.cc

Go to the documentation of this file.

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

using namespace std;

// ROOT headers
#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"

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

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

/**
Constructor
*/
DiJetAnalysis::DiJetAnalysis()
{

}

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

/**
Before using the class initialize
*/
int DiJetAnalysis::Init(double tr_max_eta, double tr_min_pt )
{
  
  // specify default eta cut
  m_max_eta=tr_max_eta;
  
  // specify default pt cut
  m_min_pt=tr_min_pt;
  
  // default Output file name
  m_outputFileName="DiJetAnalysis.root";
  m_outputRootDir="DiJet";
  
  //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]"), 50, 0.,25.);
  //m_jet_pt_log=baseAnalysis::initHist(string("jet_pt_logscale"),string("Jet Transverse Momentum -logscale-"),string("Jet p_{T} [GeV]"), 50, 0.,50.);
  //m_jet_pt_high=baseAnalysis::initHist(string("jet_pt_high"),string("Jet Transverse Momentum"),string("Jet p_{T} [GeV]"), 50, 0.,150.);
  //m_jet_pt_high_log=baseAnalysis::initHist(string("jet_pt_high_logscale"),string("Jet Transverse Momentum -logscale-"),string("Jet p_{T} [GeV]"), 50, 0.,150.);
  m_jet_pt=baseAnalysis::initHist(string("jet_pt"),string("Jet Transverse Momentum"),string("Jet p_{T} [GeV]"), 
                                 50, 0.,150.);
  m_jet_pt_log=baseAnalysis::initHist(string("jet_pt_logscale"),string("Jet Transverse Momentum -logscale-"),
                                             string("Jet p_{T} [GeV]"), 50, 0.,150.);
  m_leadingjet_pt=baseAnalysis::initHist(string("leadingjet_pt"),string("Leading Jet Transverse Momentum"),string("Jet p_{T} [GeV]"), 
                                         50, 0.,150.);
  m_leadingjet_pt_log=baseAnalysis::initHist(string("leadingjet_pt_logscale"),string("Leading Jet Transverse Momentum -logscale-"),
                                             string("Jet p_{T} [GeV]"), 50, 0.,150.);
  m_secondleadingjet_pt=baseAnalysis::initHist(string("secondleadingjet_pt"),string("Second leading Jet Transverse Momentum"),
                                               string("Jet p_{T} [GeV]"), 50, 0.,150.);
  m_secondleadingjet_pt_log=baseAnalysis::initHist(string("secondleadingjet_pt_logscale"),
                                                   string("Second leading Jet Transverse Momentum -logscale-"),
                                                   string("Jet p_{T} [GeV]"), 50, 0.,150.);
  
  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"),string("charged particle multiplicity"),100, 0., 1000.);
  m_charged_particle_temp_pt= new TH1D("DiJet_charged_particle_mean_pt","temp histogramm",200,0,200.); 
  m_charged_particle_pt=baseAnalysis::initHist("charged_particle_pt","p_{T} of charged particles","p_{T}",25,0,25.); 
  m_charged_particle_mean_pt=baseAnalysis::initHist(string("charged_particle_mean_pt"),string("Avarage p_{T} charged Particles in the Event"),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"),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.);
  
  for(unsigned int hist=0; hist< m_histVector.size(); hist++)
    {
      std::ostringstream counter_s;
      counter_s << hist;
      std::string histogram_name=m_outputRootDir;
      histogram_name+="_";
      histogram_name+=counter_s.str();
      histogram_name+="_";
      histogram_name+=m_histVector[hist]->GetName();
      m_histVector[hist]->SetName(histogram_name.c_str());
    }
  return SUCCESS;
}

/**
This is the main analysis function. 
Search the particle, get their properties and the histograms.
*/
int DiJetAnalysis::Process(HepMC::GenEvent* hepmcevt) 
{
  // define some variables
  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();
      
      // just remove the neutrinos, gammas and neutrons, ... from the
      // final state
      if(chargedParticle(pid)==NOTCHARGED) continue;
      
      // we want stable charged particles
      if (!((*p)->status() == 1)) continue;
      
      m_charged_particle_temp_pt->Fill((*p)->momentum().perp());
      m_charged_particle_pt->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);
  
  // run the jetfinder
  if(!m_inclusive_jets.size()) {
    int test = FindJet(hepmcevt);
    if(test<=0){
      cout <<"no jets found - return failure" << endl;
      //return FAILURE;
    }
  }
  
  //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_pt_high->Fill(m_inclusive_jets[0].perp());
    //m_jet_pt_high_log->Fill(m_inclusive_jets[0].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(abs(lv_firstJet.deltaR(lv_secondJet)));
      m_jet_Deta_leadingJet_secondJet->Fill(abs(lv_firstJet.eta()-lv_secondJet.eta()));
      m_jet_Dphi_leadingJet_secondJet->Fill(abs(lv_firstJet.vect().deltaPhi(lv_secondJet.vect())));
    }
  }
  return SUCCESS;
}

---