WplusJetAnalysis.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"

//WplusJet analysis header
#include "../include/WplusJetAnalysis.h"

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

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

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

/**
Before using the class initialize
*/
int WplusJetAnalysis::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="WplusJetAnalysis.root";
  m_outputRootDir="WplusJet";
  
  //declaration of histograms
  //m_evtnr=baseAnalysis::initHist(string("evtnr"),string("Event number"),string("Eventnumber"),1000, 0., 1000.);
  
  m_W_count=baseAnalysis::initHist(string("W_count"), string("Nr of W-Bosons"), string("Nr W-Boson"),16, -0.5, 15.5);
  m_W_trans_mass=baseAnalysis::initHist(string("W_trans_mass"),string("transvers W Mass"),string("m_{T,W} [GeV]"), 100, 0., 100.);
  m_W_trans_mass_log=baseAnalysis::initHist(string("W_trans_mass_logscale"),string("transvers W Mass -logscale-"),string("m_{T,W} [GeV]"), 100, 0., 100.);
  m_W_pt_high=baseAnalysis::initHist(string("W_pt_high"),string("W Transverse Momentum"),string("W p_{T} [GeV]"), 50, 0., 150.);
  m_W_pt_high_log=baseAnalysis::initHist(string("W_pt_high_logscale"),string("W Transverse Momentum -logscale-"),string("W p_{T} [GeV]"), 50, 0., 150.);
  m_W_pt=baseAnalysis::initHist(string("W_pt"),string("W Transverse Momentum"),string("W p_{T} [GeV]"), 50, 0., 25.);
  m_W_pt_log=baseAnalysis::initHist(string("W_pt_logscale"),string("W Transverse Momentum -logscale-"),string("W p_{T} [GeV]"), 50, 0., 25.);
  m_W_eta=baseAnalysis::initHist(string("W_eta"),string("#eta W-Boson"),string("#eta"), 60, -6., 6.);
  m_W_rapidity=baseAnalysis::initHist(string("W_rapidity"),string("rapidity y W-Boson"),string("y"), 60, -6., 6.);
  m_W_phi=baseAnalysis::initHist(string("W_phi"),string("#varphi W-Boson"),string("#varphi"), 64, -3.2, 3.2);
  
  m_lepton_count=baseAnalysis::initHist(string("lepton_count"), string("Nr of Leptons from W-Bosons"), string("Nr Leptons"),11, -0.5, 10.5);
  m_lepton_pt_high=baseAnalysis::initHist(string("lepton_pt_high"),string("Transverse Momentum Lepton from W-Boson"),string("lepton p_{T} [GeV]"), 50, 0., 150.);
  m_lepton_pt_high_log=baseAnalysis::initHist(string("lepton_pt_high_logscale"),string("Transverse Momentum Lepton from W-Boson -logscale-"),string("lepton p_{T} [GeV]"), 50, 0., 150.);
  m_lepton_pt=baseAnalysis::initHist(string("lepton_pt_high"),string("Transverse Momentum Lepton from W-Boson"),string("lepton p_{T} [GeV]"), 50, 0., 25.);
  m_lepton_pt_log=baseAnalysis::initHist(string("lepton_pt_high_logscale"),string("Transverse Momentum Lepton from W-Boson -logscale-"),string("lepton p_{T} [GeV]"), 50, 0., 25.);
  m_lepton_eta=baseAnalysis::initHist(string("lepton_eta"),string("#eta Lepton from W-Boson"),string("#eta"), 60, -6., 6.);
  m_lepton_phi=baseAnalysis::initHist(string("lepton_phi"),string("#varphi  Lepton from W-Boson"),string("#varphi"), 64, -3.2, 3.2);
  m_lepton_pdgid=baseAnalysis::initHist(string("lepton_pdgid"),string("PDG ID Lepton from W-Boson"),string("PDG ID"), 40, -20., 20.);
  
  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("WplusJet_charged_particle_mean_pt","temp histogramm",200,0,200.);
  m_charged_particle_pt= new TH1D("WplusJet_charged_particle_pt","p_{T} of stable charged particles",200,0,200.);
  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_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_jet_Deta_leadingJet_lepton=baseAnalysis::initHist(string("jet_Deta_leadingJet_lepton"),string("#Delta#eta leading Jet and Lepton from W-Boson"),string("#Delta#eta"), 60, 0., 6.);
//   m_jet_Deta_leadingJet_lepton_singlejet=baseAnalysis::initHist(string("jet_Deta_leadingJet_lepton_singlejet"),string("#Delta#eta leading Jet (single) and Lepton from W-Boson"),string("#Delta#eta"), 60, 0., 6.);
//   m_jet_Deta_leadingJet_lepton_multijet=baseAnalysis::initHist(string("jet_Deta_leadingJet_lepton_multijet"),string("#Delta#eta leading Jet (multi) and Lepton from W-Boson"),string("#Delta#eta"), 60, 0., 6.);
  
//   m_jet_Dphi_leadingJet_lepton=baseAnalysis::initHist(string("jet_Dphi_leadingJet_lepton"),string("#Delta#varphi leading Jet and Lepton from W-Boson"),string("#Delta#varphi"), 64, 0., 3.2);
//   m_jet_Dphi_leadingJet_lepton_singlejet=baseAnalysis::initHist(string("jet_Dphi_leadingJet_lepton_singlejet"),string("#Delta#varphi leading Jet (single) and Lepton from W-Boson"),string("#Delta#varphi"), 64, 0., 3.2);
//   m_jet_Dphi_leadingJet_lepton_multijet=baseAnalysis::initHist(string("jet_Dphi_leadingJet_lepton_multijet"),string("#Delta#varphi leading Jet (multi) and Lepton from W-Boson"),string("#Delta#varphi"), 64, 0., 3.2);
  
//   m_jet_DR_leadingJet_lepton=baseAnalysis::initHist(string("jet_DR_leadingJet_lepton"),string("#DeltaR leading Jet and Lepton from W-Boson"),string("#DeltaR"), 60, 0., 15.);
//   m_jet_DR_leadingJet_lepton_singlejet=baseAnalysis::initHist(string("jet_DR_leadingJet_lepton_singlejet"),string("#DeltaR leading Jet (single) and Lepton from W-Boson"),string("#DeltaR"), 60, 0., 15.);
//   m_jet_DR_leadingJet_lepton_multijet=baseAnalysis::initHist(string("jet_DR_leadingJet_lepton_multijet"),string("#DeltaR leading Jet (multi) and Lepton from W-Boson"),string("#DeltaR"), 60, 0., 15.);
  
  m_jet_Deta_leadingJet_lepton=baseAnalysis::initHist(string("jet_Deta_leadingJet_lepton"),string("#Delta#eta leading Jet and Lepton from W-Boson"),string("#Delta#eta"), 60, 0., 3.);
  m_jet_Deta_leadingJet_lepton_singlejet=baseAnalysis::initHist(string("jet_Deta_leadingJet_lepton_singlejet"),string("#Delta#eta leading Jet (single) and Lepton from W-Boson"),string("#Delta#eta"), 60, 0., 3.);
  m_jet_Deta_leadingJet_lepton_multijet=baseAnalysis::initHist(string("jet_Deta_leadingJet_lepton_multijet"),string("#Delta#eta leading Jet (multi) and Lepton from W-Boson"),string("#Delta#eta"), 60, 0., 3.);
  
  m_jet_Dphi_leadingJet_lepton=baseAnalysis::initHist(string("jet_Dphi_leadingJet_lepton"),string("#Delta#varphi leading Jet and Lepton from W-Boson"),string("#Delta#varphi"), 64, 0., 1.);
  m_jet_Dphi_leadingJet_lepton_singlejet=baseAnalysis::initHist(string("jet_Dphi_leadingJet_lepton_singlejet"),string("#Delta#varphi leading Jet (single) and Lepton from W-Boson"),string("#Delta#varphi"), 64, 0., 1.);
  m_jet_Dphi_leadingJet_lepton_multijet=baseAnalysis::initHist(string("jet_Dphi_leadingJet_lepton_multijet"),string("#Delta#varphi leading Jet (multi) and Lepton from W-Boson"),string("#Delta#varphi"), 64, 0., 1.);
  
  m_jet_DR_leadingJet_lepton=baseAnalysis::initHist(string("jet_DR_leadingJet_lepton"),string("#DeltaR leading Jet and Lepton from W-Boson"),string("#DeltaR"), 60, 0., 6.);
  m_jet_DR_leadingJet_lepton_singlejet=baseAnalysis::initHist(string("jet_DR_leadingJet_lepton_singlejet"),string("#DeltaR leading Jet (single) and Lepton from W-Boson"),string("#DeltaR"), 60, 0., 6.);
  m_jet_DR_leadingJet_lepton_multijet=baseAnalysis::initHist(string("jet_DR_leadingJet_lepton_multijet"),string("#DeltaR leading Jet (multi) and Lepton from W-Boson"),string("#DeltaR"), 60, 0., 6.);
  
  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 true;
  
}

/**
This is the main analysis function. 
Search the particle, get their properties and the histograms.
*/
int WplusJetAnalysis::Process(HepMC::GenEvent* hepmcevt) 
{
  
  // Some Variables or storage classes are needed
  CLHEP::HepLorentzVector lv_leptonFromW(0,0,0,0);
  CLHEP::HepLorentzVector lv_neutrinoFromW(0,0,0,0);
  HepMC::GenParticle* leptonFromW=0;
  HepMC::GenParticle* neutrinoFromW=0;
  
  // storage of input particles for jet
  //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 Wcount=0;
  int leptoncount=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* WBoson=0;
      
      // use the pdg id to identify the particles
      int pid = (*p)->pdg_id();
      
      // W decays
      // look for W-Bosons (pdg = 24) 
      //if( abs(pid) == 24 && ((*p)->status())%1000 == properStatus && (*p)->end_vertex()){
      if( abs(pid) == 24 && (*p)->end_vertex()){

        // search for decaying particles
        // since the particles decay, they should have an end 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);

        bool LastWBoson = true;
        for(HepMC::GenVertex::particle_iterator iter=firstChild ; iter!=endChild; iter++)
          {
            if( (*iter)->pdg_id() == pid ){
              LastWBoson = false;
            }
          }

        // plot Pt, eta and phi of the W (generator values)
        if(LastWBoson == true) {
          WBoson = (*p);
          m_W_pt->Fill((*p)->momentum().perp());
          m_W_pt_log->Fill((*p)->momentum().perp());
          m_W_pt_high->Fill((*p)->momentum().perp());
          m_W_pt_high_log->Fill((*p)->momentum().perp());
          m_W_eta->Fill((*p)->momentum().eta());
          m_W_rapidity->Fill(baseAnalysis::getRapidity(p));
          m_W_phi->Fill((*p)->momentum().phi());
          Wcount++;
        }
        
        bool WtoLepton=false;
        // Now check the decay products and search for a lepton and the neutrino
        for(HepMC::GenVertex::particle_iterator iter=firstChild ; iter!=endChild; iter++)
          {
            // set the properties of the neutrino Lorentzvector 
            if(abs((*iter)->pdg_id())==12 || abs((*iter)->pdg_id())==14 || abs((*iter)->pdg_id())==16){
              neutrinoFromW=(*iter);
              lv_neutrinoFromW.setPx((*iter)->momentum().px());
              lv_neutrinoFromW.setPy((*iter)->momentum().py());
              lv_neutrinoFromW.setPz((*iter)->momentum().pz());
              lv_neutrinoFromW.setE((*iter)->momentum().e());
            }
            // set the properties of the lepton Lorentzvector 
            if(abs((*iter)->pdg_id())==11 || abs((*iter)->pdg_id())==13 || abs((*iter)->pdg_id())==15) {
              leptonFromW=*iter;
              WtoLepton=true;
              m_lepton_pdgid->Fill((*iter)->pdg_id());
              lv_leptonFromW.setPx((*iter)->momentum().px());
              lv_leptonFromW.setPy((*iter)->momentum().py());
              lv_leptonFromW.setPz((*iter)->momentum().pz());
              lv_leptonFromW.setE((*iter)->momentum().e());
            }
          }
        
        if(!WtoLepton) continue;
        
        // fill the lepton histograms
        m_lepton_pt->Fill(lv_leptonFromW.perp());
        m_lepton_pt_log->Fill(lv_leptonFromW.perp());
        m_lepton_pt_high->Fill(lv_leptonFromW.perp());
        m_lepton_pt_high_log->Fill(lv_leptonFromW.perp());
        m_lepton_eta->Fill(lv_leptonFromW.eta());
        m_lepton_phi->Fill(lv_leptonFromW.phi());
        leptoncount++;
        
        // Use the standard Lorentzvector function to get the transverse Mass of the W
        // First set the z-component of the lepton and neutrino to zero and adjust the transverse energy (=transverse momentum)
        // then add the two lorentzvectors to get the transverse mass of the W-Boson
        CLHEP::HepLorentzVector lv_transverse_leptonFromW(lv_leptonFromW.px(),lv_leptonFromW.py(),0,lv_leptonFromW.perp());
        CLHEP::HepLorentzVector lv_transverse_neutrinoFromW(lv_neutrinoFromW.px(),lv_neutrinoFromW.py(),0,lv_neutrinoFromW.perp());
        CLHEP::HepLorentzVector lv_transverse_W(lv_transverse_leptonFromW);
        lv_transverse_W+=lv_transverse_neutrinoFromW;
        m_W_trans_mass->Fill(lv_transverse_W.m());
        m_W_trans_mass_log->Fill(lv_transverse_W.m());
        
      }
      
      // 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;
      
      // take all stable particles
      if(!IsFinalStateParticle((*p))) continue;
      
      // remove pi0 from charged particle list
      //if(pid==111) continue;
      
      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<<"WplusJetAnalysis::no jets found"<<std::endl;

  //Fill the histograms
  m_W_count->Fill(Wcount);
  m_lepton_count->Fill(leptoncount);
  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());
  
  m_jet_count->Fill(m_inclusive_jets.size());

  // not needed, but to be on the safe side
  m_charged_particle_temp_pt->Reset();
  
  if(m_inclusive_jets.size()){

    m_jet_pt->Fill(m_inclusive_jets[0].perp());
    m_jet_pt_log->Fill(m_inclusive_jets[0].perp());
    
    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(leptonFromW)  {
      if(leptoncount>1) 
        std::cout << "WplusJetAnalysis: WARNING: You have more than one Lepton from W-Boson in the event!" << std::endl;
      double perp1 = lv_leptonFromW.perp();
      double perp2 = lv_firstJet.perp();
      double eta1 = lv_leptonFromW.eta();
      double eta2 = lv_firstJet.eta();
      double phi1 = lv_leptonFromW.phi();
      double phi2 = lv_firstJet.phi();
      // skip leptons which are not isolated, because they are in the jet and in the W decay
      if(fabs(lv_firstJet.deltaR(lv_leptonFromW)) < 0.0005 && abs(lv_firstJet.deltaR(lv_leptonFromW)) > -0.0005
         && perp1 == perp2 && eta1==eta2 && phi1==phi2) return false;

      m_jet_DR_leadingJet_lepton->Fill(fabs(lv_firstJet.deltaR(lv_leptonFromW)));
      m_jet_Deta_leadingJet_lepton->Fill(fabs(lv_firstJet.eta()-lv_leptonFromW.eta()));
      m_jet_Dphi_leadingJet_lepton->Fill(fabs(lv_firstJet.vect().deltaPhi(lv_leptonFromW.vect())));

      // Fill the angle plots for a lepton according to the leading Jet
      // ... if we have more than one Jet
      if(m_inclusive_jets.size() > 1 ) {
        m_jet_DR_leadingJet_lepton_multijet->Fill(fabs(lv_firstJet.deltaR(lv_leptonFromW)));
        m_jet_Deta_leadingJet_lepton_multijet->Fill(fabs(lv_firstJet.eta()-lv_leptonFromW.eta()));
        m_jet_Dphi_leadingJet_lepton_multijet->Fill(fabs(lv_firstJet.vect().deltaPhi(lv_leptonFromW.vect())));
      }
      else{       // ... or only one Jet
        m_jet_DR_leadingJet_lepton_singlejet->Fill(fabs(lv_firstJet.deltaR(lv_leptonFromW)));
        m_jet_Deta_leadingJet_lepton_singlejet->Fill(fabs(lv_firstJet.eta()-lv_leptonFromW.eta()));
        m_jet_Dphi_leadingJet_lepton_singlejet->Fill(fabs(lv_firstJet.vect().deltaPhi(lv_leptonFromW.vect())));
      }
    }
  }

  return true;
}

---