ZAnalysis.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/ZAnalysis.h"

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

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

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

/**
Before using the class initialize
*/
int ZAnalysis::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="ZBosonAnalysis.root";
  m_outputRootDir="Z";
  
  //declaration of histograms
  //m_evtnr=baseAnalysis::initHist(string("evtnr"),string("Event number"),string("Eventnumber"),1000, 0., 1000.);
  
  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("leadingjet_pt_high"),string("Leading Jet Transverse Momentum"),string("Jet p_{T} [GeV]"), 200, 0.,1000.);
  m_jet_pt_log=baseAnalysis::initHist(string("leadingjet_pt_high_logscale"),string("Leading Jet Transverse Momentum -logscale-"),string("Jet p_{T} [GeV]"), 200, 0.,1000.);
  
  m_Z_count=baseAnalysis::initHist(string("Z_count"), string("Nr of Z-Bosons"), string("Nr Z-Boson"),16, -0.5, 15.5);
  m_Z_mass=baseAnalysis::initHist(string("Z_mass"),string("Z Mass"),string("m_{Z} [GeV]"), 80, 70., 110.);
  m_Z_mass_log=baseAnalysis::initHist(string("Z_mass_logscale"),string("Z Mass -logscale-"),string("m_{Z} [GeV]"), 80, 70., 110.);
  
  m_Z_pt_high=baseAnalysis::initHist(string("Z_pt_high"),string("Z Transverse Momentum"),string("Z p_{T} [GeV]"), 80, 0., 240.);
  m_Z_pt_high_log=baseAnalysis::initHist(string("Z_pt_high_logscale"),string("Z Transverse Momentum -logscale-"),string("Z p_{T} [GeV]"), 80, 0., 240.);
  m_Z_pt=baseAnalysis::initHist(string("Z_pt_high"),string("Z Transverse Momentum"),string("Z p_{T} [GeV]"), 50, 0., 25.);
  m_Z_pt_log=baseAnalysis::initHist(string("Z_pt_high_logscale"),string("Z Transverse Momentum -logscale-"),string("Z p_{T} [GeV]"), 50, 0., 25.);
  
  //change the range to (-8,8)
  m_Z_eta=baseAnalysis::initHist(string("Z_eta"),string("#eta Z-Boson"),string("#eta"), 80, -8., 8.);
  m_Z_phi=baseAnalysis::initHist(string("Z_phi"),string("#varphi Z-Boson"),string("#varphi"), 64, -3.2, 3.2);

  m_charged_particle_multiplicity=baseAnalysis::initHist(string("charged_particle_multiplicity"),string("Number of charged Particles in the Event"),string("charged particle multiplicity"),100, 0., 300.);
//   m_charged_particle_temp_pt= new TH1D("Z_charged_particle_mean_pt","temp histogramm",200,0,100.);
//   m_charged_particle_pt= new TH1D("Z_charged_particle_pt","p_{T} pf charged particles",200,0,100.);
  m_charged_particle_temp_pt=baseAnalysis::initHist(string("Z_charged_particle_temp_mean_pt"),string("temp histogramm"),string("Z charged particle temp mean p_{T}"),200,0,100.);
  m_charged_particle_pt=baseAnalysis::initHist(string("Z_charged_particle_pt"),string("p_{T} pf charged particles"),string("p_{T}"),200,0,100.);
  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_charged_particle_pdgID=baseAnalysis::initHist(string("charged_particle_pdgID"),string("charged Particles PDG ID"),string("charged particle pid"),601, -300., 300.);
  
  return true;
}

/**
This is the main analysis function. 
Search the particle, get their properties and the histograms.
*/

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

  // Some Variables or storage classes are needed
  //int properStatus = 2;

  // Reset the histogramm
  m_charged_particle_temp_pt->Reset();
  
  // Fill the eventnumber
  //m_evtnr->Fill((hepmcevt->event_number())%100);
  
  // initialize the counter variables
  int Zcount=0;
  int chargeParticlecount=0;

  // loop over the particles and select pdgid and pt
  for ( HepMC::GenEvent::particle_const_iterator p =  hepmcevt->particles_begin(); p != hepmcevt->particles_end(); ++p )
    {
      HepMC::GenParticle* ZBoson=0;

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


      // Z decays
      // look for Z-Bosons (pdg = 23) 
      //if( pid == 23 && ((*p)->status())%1000 == properStatus && (*p)->end_vertex()){
      if( pid == 23 && (*p)->end_vertex()){

        //search for decaying particles
        // since the particles should decay they should have a decay vertex
        if (!(*p)->end_vertex()) continue;

        HepMC::GenVertex::particle_iterator firstChild = (*p)->end_vertex()->particles_begin(HepMC::children);
        HepMC::GenVertex::particle_iterator endChild =  (*p)->end_vertex()->particles_end(HepMC::children);

        // plot Pt, eta, phi and mass of the Z (generator values)
        if( ((*firstChild)->pdg_id()) != pid ){
          ZBoson = (*p);

          m_Z_pt->Fill((*p)->momentum().perp());
          m_Z_pt_log->Fill((*p)->momentum().perp());
          m_Z_pt_high->Fill((*p)->momentum().perp());
          m_Z_pt_high_log->Fill((*p)->momentum().perp());
          m_Z_eta->Fill((*p)->momentum().eta());
          m_Z_phi->Fill((*p)->momentum().phi());
          m_Z_mass->Fill((*p)->momentum().m());
          m_Z_mass_log->Fill((*p)->momentum().m());
          Zcount++;
        }
      }

      // select final state particles

      // 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;
      //       // these lines are comment out and implemented in the function IsFinalStateParticle in the baseAnalysis:
      //       // first of all we need stable particle
      //       if (!((*p)->status() == 1)) 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;
      
      //       // fiducial range eta cut on 2.5
      //       if(fabs((*p)->momentum().eta()) > 2.5) continue;
      //       // minimum pt for tracks on 0.5
      //       if((*p)->momentum().perp() < 0.5) continue;
      
      // take all stable particles
      if(!IsFinalStateParticle((*p))) {
        continue;
      }
      // do not count pi0, K0_L, K0_S and K0
      //if((abs(pid)==111) || (abs(pid)==130) || (abs(pid)==310) || (abs(pid)==311)) continue;
      
      //fill values for charged stable particles
      m_charged_particle_pt->Fill((*p)->momentum().perp());
      m_charged_particle_temp_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);
  
  // check if jets were found
  //if(!m_inclusive_jets.size()) std::cout<<"ZAnalysis::no jets found"<<std::endl;

  
  //Fill the histograms
  m_jet_count->Fill(m_inclusive_jets.size()); 

  if(m_inclusive_jets.size()){
    m_jet_pt->Fill(m_inclusive_jets[0].perp());
    m_jet_pt_log->Fill(m_inclusive_jets[0].perp());
  }

  //Fill the histograms
  m_Z_count->Fill(Zcount);
  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());
  
  chargeParticlecount=0;
  
  // not needed, but to be on the safe side
  m_charged_particle_temp_pt->Reset();
  
  return true;
}

---