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: delete all the Pointer
*/
ZAnalysis::~ZAnalysis()
{
  while (!m_histVector.empty()){
    delete m_histVector.back();
    m_histVector.pop_back();
  }
}

/**
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"),100, 0., 100.);
  
  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_high"),string("Jet Transverse Momentum"),string("Jet p_{T} [GeV]"), 50, 0.,150.);
  m_jet_pt_log=baseAnalysis::initHist(string("jet_pt_high_logscale"),string("Jet Transverse Momentum -logscale-"),string("Jet p_{T} [GeV]"), 50, 0.,150.);
  
  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=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.);
  m_Z_pt_high=baseAnalysis::initHist(string("Z_pt_high"),string("Z Transverse Momentum"),string("Z p_{T} [GeV]"), 50, 0., 150.);
  m_Z_pt_high_log=baseAnalysis::initHist(string("Z_pt_high_logscale"),string("Z Transverse Momentum -logscale-"),string("Z p_{T} [GeV]"), 50, 0., 150.);
  
  //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_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"),100, 200., 300.);
  
  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 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)->end_vertex()){

        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++;
        }
      }
      // 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)==NOTCHARGED) continue;

      // now collect the tracks for the Jet algorithm
      // 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
      if(abs((*p)->momentum().eta()) > m_max_eta) continue;
      // minimum pt for tracks
      if((*p)->momentum().perp() < m_min_pt) continue;
      // check wether the track directly comes from a Z-Boson
      if(((abs(pid)==11) || (abs(pid)==13)) && trackfromPID(23,(*p))) continue;
      // check wether the track directly comes from a W-Boson      
      if(((abs(pid)==11) || (abs(pid)==13)) && trackfromPID(24,(*p))) continue;
      if(((abs(pid)==11) || (abs(pid)==13)) && trackfromPID(-24,(*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);
  
  // 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()); 

  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 SUCCESS;
}

---