ZAnalysis.cc

Go to the documentation of this file.

// ---------------------------------------------------------------------- 

#include <iostream>

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

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

#include "../include/ZAnalysis.h"

using namespace std;

// ---------------------------------------------------------------------- 
ZAnalysis::ZAnalysis() : baseAnalysis()
{}

// ---------------------------------------------------------------------- 
ZAnalysis::~ZAnalysis()
{}

// ---------------------------------------------------------------------- 
/** Z analysis initialization. */
// ---------------------------------------------------------------------- 
int ZAnalysis::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
  
    // default output file name
  m_outputFileName = "ZBosonAnalysis.root";
  m_outputRootDir = "Z";
  
    // book histograms
  m_jet_count   = initHist( "jet_count", "Nr of Jets", "Nr Jets", 16, -0.5, 15.5 );
  m_jet_pt      = initHist( "leadingjet_pt_high", "Leading Jet Transverse Momentum", "Jet p_{T} [GeV]", 200, 0.0, 1000.0 );
  m_jet_pt_log  = initHist( "leadingjet_pt_high_logscale", "Leading Jet Transverse Momentum -logscale-", 
      "Jet p_{T} [GeV]", 200, 0.0, 1000.0 );
 
    // pt histograms 
  m_Z_pt      = initHist( "Z_pt", "Z Transverse Momentum", "Z p_{T} [GeV]", 50, 0.0, 25.0 );
  m_Z_pt_log  = initHist( "Z_pt_logscale", "Z Transverse Momentum -logscale-", "Z p_{T} [GeV]", 50, 0.0, 25.0 );
  m_Z_pt_prop      = initHist( "Z_pt_prop", "Z Transverse Momentum (propagator)", "Z p_{T} [GeV]", 50, 0.0, 25.0 );
  m_Z_pt_prop_log  = initHist( "Z_pt_prop_logscale", "Z Transverse Momentum (propagator) -logscale-", "Z p_{T} [GeV]", 50, 0.0, 25.0 );
  
  m_Z_pt_high     = initHist( "Z_pt_high", "Z Transverse Momentum", "Z p_{T} [GeV]", 80, 0.0, 240.0 );
  m_Z_pt_high_log = initHist( "Z_pt_high_logscale", "Z Transverse Momentum -logscale-", 
      "Z p_{T} [GeV]", 80, 0.0, 240.0 );
  m_Z_pt_high_prop     = initHist( "Z_pt_high_prop", "Z Transverse Momentum (propagator)", "Z p_{T} [GeV]", 80, 0.0, 240.0 );
  m_Z_pt_high_prop_log = initHist( "Z_pt_high_prop_logscale", "Z Transverse Momentum (propagator) -logscale-", 
      "Z p_{T} [GeV]", 80, 0.0, 240.0 );

    // Z eta/phi
  m_Z_eta = initHist( "Z_eta", "#eta Z-Boson", "#eta", 80, -8.0, 8.0 );
  m_Z_phi = initHist( "Z_phi", "#varphi Z-Boson", "#varphi", 64, -3.2, 3.2 );
  m_Z_eta_prop = initHist( "Z_eta_prop", "#eta Z-Boson (propagator)", "#eta", 80, -8.0, 8.0 );
  m_Z_phi_prop = initHist( "Z_phi_prop", "#varphi Z-Boson (propagator)", "#varphi", 64, -3.2, 3.2 );

    // Z mass
  m_Z_mass      = initHist( "Z_mass", "Z Mass", "m_{Z} [GeV]", 80, 70.0, 110.0 );
  m_Z_mass_log  = initHist( "Z_mass_logscale", "Z Mass -logscale-", "m_{Z} [GeV]", 80, 70.0, 110.0 );
  m_Z_mass_prop      = initHist( "Z_mass_prop", "Z Mass (propagator)", "m_{Z} [GeV]", 80, 70.0, 110.0 );
  m_Z_mass_prop_log  = initHist( "Z_mass_prop_logscale", "Z Mass (propagator) -logscale-", "m_{Z} [GeV]", 80, 70.0, 110.0 );
  
  m_charged_particle_multiplicity = initHist( "charged_particle_multiplicity", 
      "Number of charged Particles in the Event", "charged particle multiplicity", 100, 0.0, 300.0 );
  m_charged_particle_temp_pt = initHist( "Z_charged_particle_temp_mean_pt", "temp histogramm", 
      "Z charged particle temp mean p_{T}", 200, 0.0, 100.0 );
  m_charged_particle_pt = initHist( "Z_charged_particle_pt", "p_{T} pf charged particles", 
      "p_{T}", 200, 0.0, 100.0 );

  m_charged_particle_mean_pt = initHist( "charged_particle_mean_pt", "Average p_{T} charged Particles in the Event",
      "charged particle #bar{p}_{T}", 50, 0.0, 10.0 );
  m_charged_particle_rms_pt = initHist( "charged_particle_rms_pt", "RMS p_{T} charged Particles in the Event",
      "charged particle #sigma(p_{T})", 50, 0.0, 10.0 );
  m_charged_particle_pdgID = initHist( "charged_particle_pdgID", "charged Particles PDG ID", 
      "charged particle pid", 601, -300.0, 300.0 );
  
  return true;
}

// ---------------------------------------------------------------------- 
/** @return flag indicating that Z boson was indeed found and the FourVector 
of the final state (as the last function's parameter). */
// ---------------------------------------------------------------------- 
bool ZAnalysis::FSVector( HepMC::GenEvent::particle_const_iterator &zBoson, ZAnalysis::EParticleType productPid, 
    HepMC::FourVector &fsVector )
{
  size_t  nNegParts = 0, 
          nPosParts = 0;

  HepMC::GenVertex::particle_iterator posParticle;        
  HepMC::GenVertex::particle_iterator negParticle;        

     // loop over children of the found Z boson
  for ( HepMC::GenVertex::particle_iterator child = (*zBoson)->end_vertex()->particles_begin( HepMC::children );
        child != (*zBoson)->end_vertex()->particles_end( HepMC::children ); ++child ) {

      // look for a pair e+e-, mu+mu-
    int pid = (*child)->pdg_id();
    if ( pid == productPid ) {         // look for specific particles only 
      posParticle = child;
      ++nPosParts;
    }
    if ( pid == -productPid ) {
      negParticle = child;
      ++nNegParts;
    }
  }
  //cout << "ID: " << productPid << ", nPosParticles: " << nPosParticles << ", nNegParticles: " << nNegParticles << endl; 
 
  if ( nPosParts == 1 && nNegParts == 1 )  {
    //cout << "==> Z boson found" << endl;

      // positive particle...  
    for ( HepMC::GenVertex::particle_iterator des = (*posParticle)->end_vertex()->particles_begin( HepMC::descendants );
        des != (*posParticle)->end_vertex()->particles_end( HepMC::descendants ); ++des ) {
      const HepMC::FourVector &desMom = (*des)->momentum();
    
      if ( IsFinalStateParticle( *des) ) {  
        fsVector.setPx( fsVector.px() + desMom.px() );
        fsVector.setPy( fsVector.py() + desMom.py() );
        fsVector.setPz( fsVector.pz() + desMom.pz() );
        fsVector.setE( fsVector.e() + desMom.e() );
      }
    }

      // negative particle...  
    for ( HepMC::GenVertex::particle_iterator des = (*negParticle)->end_vertex()->particles_begin( HepMC::descendants );
        des != (*negParticle)->end_vertex()->particles_end( HepMC::descendants ); ++des ) {
      const HepMC::FourVector &desMom = (*des)->momentum();
    
      if ( IsFinalStateParticle( *des) ) {  
        fsVector.setPx( fsVector.px() + desMom.px() );
        fsVector.setPy( fsVector.py() + desMom.py() );
        fsVector.setPz( fsVector.pz() + desMom.pz() );
        fsVector.setE( fsVector.e() + desMom.e() );
      }
    }

    if ( fsVector.m() == 0 ) {
      return false; 
    }

    return true;
  } 
  
  return false;
} 

// ---------------------------------------------------------------------- 
/** This is the main analysis function. Search for Z bosons, get their
properties and fill the histograms. */
// ---------------------------------------------------------------------- 
int ZAnalysis::Process( HepMC::GenEvent* hepmcevt ) 
{
    // reset auxiliary histogramm
  m_charged_particle_temp_pt->Reset();

  size_t cntZ = 0;
  size_t cntChargedParts = 0;

    // loop over particles in the current event 
  for ( HepMC::GenEvent::particle_const_iterator p =  hepmcevt->particles_begin(); p != hepmcevt->particles_end(); ++p ) {
    int pid = (*p)->pdg_id(); 
    //cout << "pid = " << pid << endl;

    HepMC::FourVector eeFs( 0.0, 0.0, 0.0, 0.0 );
    HepMC::FourVector mumuFs( 0.0, 0.0, 0.0, 0.0 );

    if ( pid == 23 && (*p)->end_vertex() ) {   // Z boson PDG id=23 and should decay
        // --------------- propagator
      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 ){
        const HepMC::FourVector &mom = (*p)->momentum();

        m_Z_pt_prop->Fill( mom.perp() );
        m_Z_pt_prop_log->Fill( mom.perp() );
        m_Z_pt_high_prop->Fill( mom.perp() );
        m_Z_pt_high_prop_log->Fill( mom.perp() );
        m_Z_eta_prop->Fill( mom.eta() );
        m_Z_phi_prop->Fill( mom.phi() );
        m_Z_mass_prop->Fill( mom.m() );
        m_Z_mass_prop_log->Fill( mom.m() );
      } // --------------- 

      if ( FSVector( p, ZAnalysis::Electron, eeFs ) ) {  // e+e- 
        FillZHistograms( eeFs );
      } 
      if ( FSVector( p, ZAnalysis::Muon, mumuFs ) ) {    // mu+mu- 
        FillZHistograms( mumuFs );
      } 
      ++cntZ;
    } 
    
    if( !IsFinalStateParticle( *p ) ) {
      continue;
    }

      // 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;
    }
    ++cntChargedParts;
 
      // fill histograms for charged stable particles
    m_charged_particle_pt->Fill( (*p)->momentum().perp() );
    m_charged_particle_temp_pt->Fill( (*p)->momentum().perp() );
    m_charged_particle_pdgID->Fill( pid );
  }

    // fill rest of the histograms
  m_charged_particle_multiplicity->Fill( cntChargedParts );
  
  m_charged_particle_mean_pt->Fill( m_charged_particle_temp_pt->GetMean() );
  m_charged_particle_rms_pt->Fill( m_charged_particle_temp_pt->GetRMS() );
  
  m_jet_count->Fill( m_inclusive_jets.size() ); 

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

  return true;
}

// ---------------------------------------------------------------------- 
/** Fill histograms specific for the Z analysis. */
// ---------------------------------------------------------------------- 
void ZAnalysis::FillZHistograms( const HepMC::FourVector &fsVec )
{
  const double &pt = fsVec.perp();

  //cout << "pt: " << fsVec.perp() << ", m: " << fsVec.m() << ", phi: " << fsVec.phi() << endl;

  m_Z_pt->Fill( pt );
  m_Z_pt_log->Fill( pt );
  m_Z_pt_high->Fill( pt );
  m_Z_pt_high_log->Fill( pt );

  m_Z_phi->Fill( fsVec.phi() );
        m_Z_eta->Fill( fsVec.eta());
        m_Z_mass->Fill( fsVec.m() );
        m_Z_mass_log->Fill( fsVec.m() );
}

---