UEAnalysis.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"
#define PI 3.141593
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"

//UE analysis header
#include "../include/UEAnalysis.h"

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

/// empty default constructor
UEAnalysis::UEAnalysis() {

}

/// empty default destructor
UEAnalysis::~UEAnalysis() {

}

int UEAnalysis::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="UEAnalysis.root";
  m_outputRootDir="UE";
  
  char name1[20],title1[20];
  char name2[20],title2[20];
  char name3[20],title3[20];
  //declaration of histograms
  for (int i=0; i<26; i++){
    //just give a different name to the histograms
    //for <Ncharged> 
    sprintf(name1,"NchargedToward_%2d",i);
    sprintf(name2,"NchargedTransverse_%2d",i);
    sprintf(name3,"NchargedAway_%2d",i);
    
    sprintf(title1,"NchargedToward_%2d",i);
    sprintf(title2,"NchargedTransverse_%2d",i);
    sprintf(title2,"NchargedAway_%2d",i);
    
    m_NchargedToward[i] = new TH1D(name1,title1, 20, 0, 100);
    m_NchargedTransverse[i] = new TH1D(name2,title2, 20, 0, 100);
    m_NchargedAway[i] = new TH1D(name3,title3, 20, 0, 100);
//     m_NchargedToward[i] =baseAnalysis::initHist(name1,name1,title1, 20, 0, 100);
//     m_NchargedTransverse[i] =baseAnalysis::initHist(name2,name2,title2, 20, 0, 100);
//     m_NchargedAway[i] =baseAnalysis::initHist(name3,name3,title3, 20, 0, 100);

   //for <Ptsum>
    sprintf(name1,"PtsumToward_%2d",i);
    sprintf(name2,"PtsumTransverse_%2d",i);
    sprintf(name3,"PtsumAway_%2d",i);
    
    sprintf(title1,"PtsumToward_%2d",i);
    sprintf(title2,"PtsumTransverse_%2d",i);
    sprintf(title2,"PtsumAway_%2d",i);
    
    m_PtsumToward[i] = new TH1D(name1,title1, 20, 0, 100);
    m_PtsumTransverse[i] = new TH1D(name2,title2, 20, 0, 100);
    m_PtsumAway[i] = new TH1D(name3,title3, 20, 0, 100);
//     m_PtsumToward[i] =baseAnalysis::initHist(name1,name1,title1, 20, 0, 100);
//     m_PtsumTransverse[i] =baseAnalysis::initHist(name2,name2,title2, 20, 0, 100);
//     m_PtsumAway[i] =baseAnalysis::initHist(name3,name3,title3, 20, 0, 100);

  }


  //setting up the variable bin size
  m_nbin_pT = 25;
  
  for (int i=0; i<11; i++) {
    m_nbinRange_pT[i] = i;
  }
  for (int i=11; i<21; i++) {
    m_nbinRange_pT[i] = 10.0+ 3.0*(i-10);
  }
  m_nbinRange_pT[21] = 50.; 
  m_nbinRange_pT[22] = 60.; 
  m_nbinRange_pT[23] = 70.; 
  m_nbinRange_pT[24] = 80.; 
  m_nbinRange_pT[25] = 120.; 
  
  //booking histograms
  m_NchargedMeanToward = baseAnalysis::initHistVariableBin( string("NchargedMeanToward"), string("NchargedMeanToward"), string("Pt_leadingJet"), m_nbin_pT, m_nbinRange_pT);
  m_NchargedMeanTransverse = baseAnalysis::initHistVariableBin( string("NchargedMeanTransverse"), string("NchargedMeanTransverse"), string("Pt_leadingJet"), m_nbin_pT, m_nbinRange_pT);
  m_NchargedMeanAway = baseAnalysis::initHistVariableBin( string("NchargedMeanAway"), string("NchargedMeanAway"), string("Pt_leadingJet"),m_nbin_pT, m_nbinRange_pT);
 
  m_PtsumMeanToward = baseAnalysis::initHistVariableBin( string("PtsumMeanToward"), string("PtsumMeanToward"), string("Pt_leadingJet"), m_nbin_pT, m_nbinRange_pT);
  m_PtsumMeanTransverse = baseAnalysis::initHistVariableBin( string("PtsumMeanTransverse"), string("PtsumMeanTransverse"), string("Pt_leadingJet"), m_nbin_pT, m_nbinRange_pT);
  m_PtsumMeanAway = baseAnalysis::initHistVariableBin( string("PtsumMeanAway"), string("PtsumMeanAway"), string("Pt_leadingJet"),m_nbin_pT, m_nbinRange_pT);
 
  m_Njet =  baseAnalysis::initHist( string("Njet"), string("Njet"), string("jet Number"), 16, -0.5, 15.5);
  m_Ptjet =  baseAnalysis::initHist( string("Ptjet"), string("jet p_{T}"), string("jet p_{T}"), 50, 0, 200);
  m_Ptjet_log =  baseAnalysis::initHist( string("Ptjet_logscale"), string("jet p_{T}"), string("jet p_{T} - logscale -"), 50, 0, 200);
  m_Ptleadingjet =  baseAnalysis::initHist( string("Ptleadingjet"), string("leading jet p_{T}"), string("leading jet p_{T}"), 50, 0, 200);
  m_Ptleadingjet_log =  baseAnalysis::initHist( string("Ptleadingjet_logscale"), string("leading jet p_{T}"), 
                                                string("leading jet p_{T} - logscale -"), 50, 0, 200);

  m_charged_particle_pt=baseAnalysis::initHist("charged_particle_pt","p_{T} of charged particles","p_{T}",25,0,25.); 

  return true; 
}


int UEAnalysis::Process(HepMC::GenEvent* hepmcevt){ 

  // storage of input particles for jet
  double leadingJetPt = 0.0;
  CLHEP::HepLorentzVector lv_leadingJet(0,0,0,0);
  
  // check if jets were found
  // if(!m_inclusive_jets.size()) {
  //   std::cout<<"UEAnalysis::no jets found"<<std::endl;
  //   return 0;
  // }
  if(m_inclusive_jets.size()) {
    // use 4 vector of leading jet to define towards/away/transverse region
    lv_leadingJet.setPx(m_inclusive_jets[0].px());
    lv_leadingJet.setPy(m_inclusive_jets[0].py());
    lv_leadingJet.setPz(m_inclusive_jets[0].pz());
    lv_leadingJet.setE(m_inclusive_jets[0].e());
    
    leadingJetPt = m_inclusive_jets[0].perp(); 
    // fill jet histograms
    for(unsigned int  i=0; i<m_inclusive_jets.size(); i++) {
      m_Ptjet->Fill(m_inclusive_jets[i].perp());
      m_Ptjet_log->Fill(m_inclusive_jets[i].perp());
    }
    m_Njet->Fill(m_inclusive_jets.size()); 
    m_Ptleadingjet->Fill(m_inclusive_jets[0].perp());
    m_Ptleadingjet_log->Fill(m_inclusive_jets[0].perp());
  }
  int NchargedToward[25];
  int NchargedTransverse[25];
  int NchargedAway[25];
  
  double PtsumToward[25];
  double PtsumTransverse[25];
  double PtsumAway[25];
  
  for(int i=0; i<m_nbin_pT; i++){
    NchargedToward[i]=0;
    NchargedTransverse[i]=0;
    NchargedAway[i]=0;
    
    PtsumToward[i]=0;
    PtsumTransverse[i]=0;
    PtsumAway[i]=0;
  }
  
  // find the leadjet bin for this event
  int ileadjet = -1;
  //int ileadjet;
  for(int i=0; i<m_nbin_pT; i++) {
    if( leadingJetPt > m_nbinRange_pT[i] && leadingJetPt < m_nbinRange_pT[i+1])
      ileadjet = i;
  }
  
  // loop over particles in the event
  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();
    
    // cut out the neutral particles (charge is not stored in HepMC)
    // just remove the neutrinos, gammas, neutral pions and
    // neutrons, K0s, ... from the final state
    if(!chargedParticle(pid)) 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;

    // check for final state particles, be careful, maybe there is another eta cut than 2.5 necessary
    if(!IsFinalStateParticle((*p))) continue;

    m_charged_particle_pt->Fill((*p)->momentum().perp());
    
    // to skip all the events with no jet
    if(ileadjet == -1) continue;
    //if(!ileadjet) continue;
    double pt = (*p)->momentum().perp();
    double px = (*p)->momentum().px();
    double py = (*p)->momentum().py();
    double pz = (*p)->momentum().pz();
    double pe = (*p)->momentum().e();
    
    CLHEP::HepLorentzVector lv_chargedParticle(px,py,pz,pe);
    
    //calculate for <Ncharged> vs pt of leading jet at toward region
    if( fabs(lv_chargedParticle.vect().deltaPhi(lv_leadingJet.vect())) < PI/3.0 ) { 
      NchargedToward[ileadjet]++;       
      PtsumToward[ileadjet]+=pt;
    }  
    
    //calculate for <Ncharged> vs pt of leading jet at Transverse region
    if( fabs(lv_chargedParticle.vect().deltaPhi(lv_leadingJet.vect())) >= PI/3.0 &&
        fabs(lv_chargedParticle.vect().deltaPhi(lv_leadingJet.vect())) < 2.0*PI/3.0 ) { 
      NchargedTransverse[ileadjet]++;
      PtsumTransverse[ileadjet]+=pt; 
    }
    //calculate for <Ncharged> vs pt of leading jet at Away region
    if( fabs(lv_chargedParticle.vect().deltaPhi(lv_leadingJet.vect())) >= 2.0*PI/3.0 ) { 
      NchargedAway[ileadjet]++;
      PtsumAway[ileadjet]+=pt;
    }  
      
    m_NchargedToward[ileadjet]->Fill(NchargedToward[ileadjet], 1.);
    m_NchargedTransverse[ileadjet]->Fill(NchargedTransverse[ileadjet], 1.);
    m_NchargedAway[ileadjet]->Fill(NchargedAway[ileadjet], 1.);
    
    m_PtsumToward[ileadjet]->Fill(PtsumToward[ileadjet], 1.);
    m_PtsumTransverse[ileadjet]->Fill(PtsumTransverse[ileadjet], 1.);
    m_PtsumAway[ileadjet]->Fill(PtsumAway[ileadjet], 1.);
    
  }
  
  return true;
} //end of Process

std::vector<TH1D*> UEAnalysis::averagedHistograms() {

  //Fill Histograms at Toward, Transverse and Away region
  for(int i=0; i<m_nbin_pT; i++){
    m_NchargedMeanToward->SetBinContent( i+1, m_NchargedToward[i]->GetMean() );  
    m_NchargedMeanToward->SetBinError( i+1, m_NchargedToward[i]->GetRMS() );  
    
    m_NchargedMeanTransverse->SetBinContent( i+1, m_NchargedTransverse[i]->GetMean() );  
    m_NchargedMeanTransverse->SetBinError( i+1, m_NchargedTransverse[i]->GetRMS() );  
    
    m_NchargedMeanAway->SetBinContent( i+1, m_NchargedAway[i]->GetMean() );  
    m_NchargedMeanAway->SetBinError( i+1, m_NchargedAway[i]->GetRMS() );  
    
    m_PtsumMeanToward->SetBinContent( i+1, m_PtsumToward[i]->GetMean() );  
    m_PtsumMeanToward->SetBinError( i+1, m_PtsumToward[i]->GetRMS() );  
    
    m_PtsumMeanTransverse->SetBinContent( i+1, m_PtsumTransverse[i]->GetMean() );  
    m_PtsumMeanTransverse->SetBinError( i+1, m_PtsumTransverse[i]->GetRMS() );  
    
    m_PtsumMeanAway->SetBinContent( i+1, m_PtsumAway[i]->GetMean() );  
    m_PtsumMeanAway->SetBinError( i+1, m_PtsumAway[i]->GetRMS() ); 
  }
  //m_histVector.insert(m_histVector.end(), m_histVectorVariableBin.begin(), m_histVectorVariableBin.end());
  //m_histVector.insert(m_histVector.begin()+2, m_histVectorVariableBin.begin(), m_histVectorVariableBin.end());

  return m_histVector;
} //end of averagedHistograms

// ---------------------------------------------------------------------- 
int UEAnalysis::Finalize()
{
  return baseAnalysis::Finalize();
}

// ---------------------------------------------------------------------- 
int UEAnalysis::Finalize(TFile* output){
  
  averagedHistograms();  

  baseAnalysis::Finalize(output);
  
  return true;
}

---