WtaunuAnalysis.cc

Go to the documentation of this file.

#include <iostream>
#include <sstream>
#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 "TMath.h"
#include "TLorentzVector.h"

// fastjet headers
#include "fastjet/PseudoJet.hh"
#include "fastjet/ClusterSequence.hh"
#include "fastjet/JetDefinition.hh"
#include "fastjet/SISConePlugin.hh"

// header for this analysis class
#include "../include/WtaunuAnalysis.h"

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

/// empty default constructor
WtaunuAnalysis::WtaunuAnalysis()
{
}

/// empty default destructor
WtaunuAnalysis::~WtaunuAnalysis()
{
}


int WtaunuAnalysis::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="WtaunuAnalysis.root";
  m_outputRootDir="Wtaunu";

  //declaration of histograms
  m_W_count=initHist("W_count","number of W-bosons in the event","W-bosons multiplicity",6, 0.,3.);
  m_W_charge=initHist("W_charge","charge of W-bosons in the event","W-bosons charge",10, -2., 2.);
  m_W_mass=initHist("W_mass","W mass","m_{W} [GeV]", 80, 60., 100.);
  m_W_mt=initHist(string("W_mt"),string("transverse W mass"),string("m_{T,W} [GeV]"), 120, 0., 120.);
  m_W_trans_mass=initHist(string("W_trans_mass"),string("transverse W mass"),string("m_{T,W} [GeV]"), 120, 0., 120.);
  m_W_pt=initHist("W_pt","transverse momentum of W","p_{T} [GeV]",200, 0., 200.);
  m_W_eta=initHist("W_eta","eta of W and anti W","#eta",60, -6., 6.);
  m_W_phi=initHist("W_phi","phi of W and anti W","#phi",48, -3.15, 3.15);
  m_W_rapidity=initHist("W_rapidity","rapidity y W-Boson","y", 60, -6., 6.);

  m_charged_particle=initHist("charged_particle","number of charged particles in the event","charged particle multiplicity",100, 0., 200.);
  m_charged_particle_pdgID=initHist("charged_particle_pdgID","charged Particles PDG ID","charged particle pid",601, -300., 300.);
  m_ptstable_charged=initHist("ptstable_charged","transveral momentum of charged stable particles","p_{T} [GeV]",200, 0., 20.);
  m_etastable_charged=initHist("etastable_charged","eta of charged stable particles","#eta",100, -6., 6.);
  m_phistable_charged=initHist("phistable_charged","phi of charged stable particles","#phi",100, -3.15, 3.15);

  m_tau_mass=initHist("tau_mass","tau mass","m_{#tau} [GeV]", 50, 0., 3.);
  m_tau_pt=initHist("tau_pt","transverse momentum of tau","p_{T} [GeV]",200, 0., 200.);
  m_tau_eta=initHist(string("tau_eta"),string("eta of tau and anti tau"),string("#eta"),60, -6., 6.);
  m_tau_phi=initHist(string("tau_phi"),string("phi of tau and anti tau"),string("#phi"),48, -3.15, 3.15);
  m_nTrack_tau=initHist(string("nTrack_tau"),string("number of tracks of tau decay"),string("number of tracks"),16, -0.5, 15.5);

  m_Delta_Phi_TauNeutrino=initHist("Delta_Phi_TauNeutrino","phi between tau and neutrino","#Delta #phi",48, -4.15, 4.15);

  m_nu_pt=initHist("nu_pt","transverse momentum of neutrino","p_{T} [GeV]",200, 0., 200.);
  m_nu_px=initHist("nu_px","momentum in x direction of neutrino","p_{x} [GeV]",301, -150., 150.);
  m_nu_py=initHist("nu_py","momentum in y direction of neutrino","p_{y} [GeV]",301, -150., 150.);
  m_nu_pz=initHist("nu_pz","momentum in z direction of neutrino","p_{z} [GeV]",301, -150., 150.);
  m_zoom_nu_px1=initHist("nu_zoom_px", "momentum in x direction of neutrino zoomed", "p_{x} [GeV]",123, -20., 20.);
  m_zoom_nu_py1=initHist("nu_zoom_py", "momentum in y direction of neutrino zoomed", "p_{y} [GeV]",123, -20., 20.);
  m_zoom_nu_pz1=initHist("nu_zoom_pz", "momentum in z direction of neutrino zoomed", "p_{z} [GeV]",123, -20., 20.);

  return true; 
}

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

  //some variables for storage are needed
  int W_barcode=0;
  int Wbar_barcode=0;
  int tau_barcode=0;
  int taubar_barcode=0;
  int nu_tau_barcode=0;
  int nubar_tau_barcode=0;

  //variables to count the number of W bosons and charged particles
  int Wcount=0;
  int charged_particle_counter=0;

  int properStatus = 0;

  bool tau_found = false;
  bool anti_tau_found = false;
  bool neutrino_found = false;
  bool anti_neutrino_found = false;

  //calculating rapidity
  double rapidity= -999;

  // loop over the particles and select pdgid and pt
  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();
      
      //search for decaying W+ bosons
      if (pid == 24 && (*p)->end_vertex())
        { 
          //prove decay of the W+ boson
          if (!(*p)->end_vertex()) continue;
          
          //loop over all decay particles of the W+ boson
          for(HepMC::GenVertex::particle_iterator iter=(*p)->end_vertex()->particles_begin(HepMC::descendants); 
              iter!=(*p)->end_vertex()->particles_end(HepMC::descendants); iter++)
            {
              if( (*iter)->pdg_id() != pid )
                {
                  //search for anti-tau and tau_neutrino as decay products
                  if ((*iter)->pdg_id() == -15) 
                    {
                      taubar_barcode = (*iter)->barcode();
                      anti_tau_found=true;
                    }
                  if ((*iter)->pdg_id() == 16)
                    {
                      nu_tau_barcode = (*iter)->barcode();
                      neutrino_found=true;
                    }
                }
            }
          
          //if W+ decays into anti-tau and tau_neutrino, then fill barcode
          if (anti_tau_found==true && neutrino_found==true )
            {
              W_barcode = (*p)->barcode();
            }
          //calculate rapidity of W boson
          rapidity = getRapidity(p);
        }
         
      //search for decaying W- bosons
      if (pid == -24 && (*p)->end_vertex())
        { 
          //prove decay of the W- boson
          if (!(*p)->end_vertex()) continue;
          
          //loop over all decay particles of the W- boson
          for(HepMC::GenVertex::particle_iterator iter=(*p)->end_vertex()->particles_begin(HepMC::descendants); 
              iter!=(*p)->end_vertex()->particles_end(HepMC::descendants); iter++)
            {
              if( (*iter)->pdg_id() != pid )
                {
                  //search for tau and tau_anti-neutrino as decay products
                  if ((*iter)->pdg_id() == 15) 
                    {
                      tau_barcode = (*iter)->barcode();
                      tau_found=true;
                    }
                  if ((*iter)->pdg_id() == -16)
                    {
                      nubar_tau_barcode = (*iter)->barcode();
                      anti_neutrino_found=true;
                    }
                }
            }

          //if W- decays into tau and anti-tau_neutrino, then fill barcode
          if (tau_found==true && anti_neutrino_found==true )
            {
              Wbar_barcode = (*p)->barcode();
            }
          //calculate rapidity of W boson
          rapidity = getRapidity(p);
        }
      
      //check if final state particle
      if(!IsFinalStateParticle((*p))) continue;
      
      //check all the charge final state particles
      if(!chargedParticle(pid)) continue;
      
      //write out the pt, eta and phi for the charged stable paricles
      m_ptstable_charged->Fill( (*p)->momentum().perp() );
      m_etastable_charged->Fill( (*p)->momentum().eta() );
      m_phistable_charged->Fill( (*p)->momentum().phi() );
      m_charged_particle_pdgID->Fill(pid);
      charged_particle_counter++;
      
    } //end particle loop
  
  //Fill number of charged particles
  m_charged_particle->Fill(charged_particle_counter);
  
  //Define W+ and W- 
  HepMC::GenParticle *W = hepmcevt->barcode_to_particle(W_barcode);
  HepMC::GenParticle *Wbar = hepmcevt->barcode_to_particle(Wbar_barcode);
  
  if (W)
    {
      //calculate transverse mass of the W-boson
      CLHEP::HepLorentzVector W_momentum(0.,0.,0.,0.);
      W_momentum.setPx(W->momentum().px());
      W_momentum.setPy(W->momentum().py());
      W_momentum.setPz(W->momentum().pz());
      W_momentum.setE(W->momentum().e());
      
      m_W_pt->Fill( W->momentum().perp() );
      m_W_eta->Fill( W->momentum().eta() );
      m_W_phi->Fill( W->momentum().phi() );
      m_W_mass->Fill(W->momentum().m());
      m_W_rapidity->Fill(rapidity);
      m_W_mt->Fill(W_momentum.mt());
      m_W_charge->Fill(+1);
      Wcount++;
    }

  if (Wbar)
    {
      //calculate transverse mass of the W-boson
      CLHEP::HepLorentzVector Wbar_momentum(0.,0.,0.,0.);
      Wbar_momentum.setPx(Wbar->momentum().px());
      Wbar_momentum.setPy(Wbar->momentum().py());
      Wbar_momentum.setPz(Wbar->momentum().pz());
      Wbar_momentum.setE(Wbar->momentum().e());
      
      m_W_pt->Fill( Wbar->momentum().perp() );
      m_W_eta->Fill( Wbar->momentum().eta() );
      m_W_phi->Fill( Wbar->momentum().phi() ) ;
      m_W_mass->Fill(Wbar->momentum().m());
      m_W_rapidity->Fill(rapidity);
      m_W_mt->Fill(Wbar_momentum.mt());
      m_W_charge->Fill(-1);
      Wcount++;
    }

  m_W_count->Fill(Wcount);
  
  //decay particles tau and nu
  
  HepMC::GenParticle *tau = hepmcevt->barcode_to_particle(tau_barcode);
  HepMC::GenParticle *taubar = hepmcevt->barcode_to_particle(taubar_barcode);
  HepMC::GenParticle *nu_tau = hepmcevt->barcode_to_particle(nu_tau_barcode);
  HepMC::GenParticle *nubar_tau = hepmcevt->barcode_to_particle(nubar_tau_barcode);

  if (tau && nubar_tau)
    {
      m_tau_pt->Fill( tau->momentum().perp() );
      m_tau_mass->Fill(tau->momentum().m());
      m_tau_eta->Fill(tau->momentum().eta());
      m_tau_phi->Fill(tau->momentum().phi());
      
      m_nu_pt->Fill(nubar_tau->momentum().perp());
      m_nu_px->Fill(nubar_tau->momentum().px());
      m_nu_py->Fill(nubar_tau->momentum().py());
      m_nu_pz->Fill(nubar_tau->momentum().pz());
      m_zoom_nu_px1->Fill(nubar_tau->momentum().px());
      m_zoom_nu_py1->Fill(nubar_tau->momentum().py());
      m_zoom_nu_pz1->Fill(nubar_tau->momentum().pz());
      
      //calculate angle phi between the decay products
      double tau_phi = 0;
      double nubar_tau_phi = 0;
      double DeltaPhi = 0;
      
      tau_phi = tau->momentum().phi();
      nubar_tau_phi = nubar_tau->momentum().phi();
      DeltaPhi = tau_phi - nubar_tau_phi;
      if( DeltaPhi > TMath::Pi() ||  DeltaPhi == TMath::Pi() )
        {
          DeltaPhi = DeltaPhi - 2 * TMath::Pi() ;
        }
      if( DeltaPhi < - TMath::Pi() ||  DeltaPhi == - TMath::Pi() )
        {
          DeltaPhi = DeltaPhi + 2 * TMath::Pi();
        }
      m_Delta_Phi_TauNeutrino->Fill(DeltaPhi);
      
      //calculate transverse mass of W-boson my changing the reference system
      //use Lorentzvectors of the decay particles and set the z-component to zero
      CLHEP::HepLorentzVector lv_transverse_leptonFromW(tau->momentum().px(),tau->momentum().py(),0,tau->momentum().perp());
      CLHEP::HepLorentzVector lv_transverse_neutrinoFromW(nubar_tau->momentum().px(),nubar_tau->momentum().py(),0,nubar_tau->momentum().perp());
      //add the Lorentzvector of the neutrino and lepton
      CLHEP::HepLorentzVector lv_transverse_W(lv_transverse_leptonFromW+lv_transverse_neutrinoFromW);
      //calculate the mass which is according to the reference system the transverse mass
      m_W_trans_mass->Fill(lv_transverse_W.m());
    }

  if (taubar && nu_tau)
    {
      m_tau_pt->Fill( taubar->momentum().perp() );
      m_tau_mass->Fill(taubar->momentum().m());
      m_tau_eta->Fill(taubar->momentum().eta());
      m_tau_phi->Fill(taubar->momentum().phi());
      
      m_nu_pt->Fill(nu_tau->momentum().perp());
      m_nu_px->Fill(nu_tau->momentum().px());
      m_nu_py->Fill(nu_tau->momentum().py());
      m_nu_pz->Fill(nu_tau->momentum().py());
      m_zoom_nu_px1->Fill(nu_tau->momentum().px());
      m_zoom_nu_py1->Fill(nu_tau->momentum().py());
      m_zoom_nu_pz1->Fill(nu_tau->momentum().pz());
      
      double taubar_phi = 0;
      double nu_tau_phi = 0;
      double DeltaPhi2 = 0;
      
      taubar_phi = taubar->momentum().phi();
      nu_tau_phi = nu_tau->momentum().phi();
      DeltaPhi2 = taubar_phi - nu_tau_phi;
      
      if( DeltaPhi2 > TMath::Pi() ||  DeltaPhi2 == TMath::Pi() )
        {
          DeltaPhi2 = DeltaPhi2 - 2 * TMath::Pi() ;
        }
      if( DeltaPhi2 < - TMath::Pi() ||  DeltaPhi2 == - TMath::Pi() )
        {
          DeltaPhi2 = DeltaPhi2 + 2 * TMath::Pi();
        }
      m_Delta_Phi_TauNeutrino->Fill(DeltaPhi2);      

      //calculate transverse mass of W-boson
      CLHEP::HepLorentzVector lv_transverse_leptonFromW(taubar->momentum().px(),taubar->momentum().py(),0,taubar->momentum().perp());
      CLHEP::HepLorentzVector lv_transverse_neutrinoFromW(nu_tau->momentum().px(),nu_tau->momentum().py(),0,nu_tau->momentum().perp());
      CLHEP::HepLorentzVector lv_transverse_W(lv_transverse_leptonFromW+lv_transverse_neutrinoFromW);
      m_W_trans_mass->Fill(lv_transverse_W.m());
    }

  //checking for children of tau
  if (tau && tau->end_vertex())
    {
      // search for decaying taus
      if (!tau->end_vertex()) return false;
      
      //number of tracks in tau decay
      int nTrack = 0; 
      
      //all generations as descendants
      for(HepMC::GenVertex::particle_iterator iChild=tau->end_vertex()->particles_begin(HepMC::descendants); 
          iChild!=tau->end_vertex()->particles_end(HepMC::descendants); iChild++)
        {
          if( (*iChild)->pdg_id() != tau->pdg_id() )
            {
              int iChild_pid = (*iChild)->pdg_id();
              properStatus = 1;
              
              if( chargedParticle(iChild_pid)==true && (*iChild)->status()%1000 == properStatus && !(*iChild)->end_vertex())
                {
                  // search for stable decay products of the taus
                  if ((*iChild)->end_vertex()) continue;
                  nTrack = nTrack + 1;
                }
            }
        }
      m_nTrack_tau->Fill(nTrack);
    }
  
  if (taubar && taubar->end_vertex())
    {
      // search for decaying taubars
      if (!taubar->end_vertex()) return false;
      
      //number of tracks in tau decay
      int nTrack = 0; 
      
      //all generations as descendants
      for(HepMC::GenVertex::particle_iterator iChild=taubar->end_vertex()->particles_begin(HepMC::descendants); 
          iChild!=taubar->end_vertex()->particles_end(HepMC::descendants); iChild++)
        {
          if( (*iChild)->pdg_id() != taubar->pdg_id() )
            {
              int iChild_pid = (*iChild)->pdg_id();
              properStatus = 1;
              
              if( (*iChild)->status()%1000 == properStatus && !(*iChild)->end_vertex() && chargedParticle(iChild_pid)==true)
                {
                  //search for stable decay products of the taubars
                  if ((*iChild)->end_vertex()) continue;
                  nTrack = nTrack + 1;
                }
            }
        }
      m_nTrack_tau->Fill(nTrack);
    }
  
  return true;
} //end of Process

---