TopAnalysis.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 "CLHEP/Vector/LorentzVector.h"

using namespace std;

// ROOT headers
#include "TH1.h"
#include "TH2.h"
#include "TFile.h"
#include "TMath.h"
#include "TLorentzVector.h"

//top analysis header
#include "../include/TopAnalysis.h"

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

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

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


int TopAnalysis::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="TopAnalysis.root";
  m_outputRootDir="Top";
  
  //declaration of histograms
  m_evtnr=initHist(string("evtnr"),string("Event number"),string("Eventnumber"),1000, 0., 1000.);
  m_toppt=initHist(string("toppt"),string("transversal momentum of top and antitop"),string("p_{T} [GeV]"),300, 0., 1500.);
  m_toppt_log=initHist(string("toppt_logscale"),string("transversal momentum of top and antitop - logscale -"),string("p_{T} [GeV]"),300, 0., 1500.);
  m_topeta=initHist(string("topeta"),string("eta of top and antitop"),string("#eta"),60, -6., 6.);
  m_topphi=initHist(string("topphi"),string("phi of top and antitop"),string("#phi"),48, -3.15, 3.15);
  m_ptstable=initHist(string("ptstable"),string("transveral momentum of stable particles (without neutrinos)"),
                      string("p_{T} [GeV]"),400, 0., 100.);
  m_ptstable_log=initHist(string("ptstable_logscale"),
                          string("transveral momentum of stable particles (without neutrinos) - logscale -"), 
                          string("p_{T} [GeV]"),400, 0., 100.);
  m_etastable=initHist(string("etastable"),string("eta of stable particles (without neutrinos)"),string("#eta"),100, -6., 6.);
  m_phistable=initHist(string("phistable"),string("phi of stable particles (without neutrinos)"),string("#phi"),100, -3.15, 3.15);
  m_ptstable_charged=initHist(string("ptstable_charged"),string("transveral momentum of charged stable particles"),
                              string("p_{T} [GeV]"),400, 0., 100.);
  m_ptstable_charged_log=initHist(string("ptstable_charged_logscale"),
                                  string("transveral momentum of charged stable particles - logscale -"),
                                  string("p_{T} [GeV]"),400, 0., 100.);
  m_etastable_charged=initHist(string("etastable_charged"),string("eta of charged stable particles"),string("#eta"),100, -6., 6.);
  m_phistable_charged=initHist(string("phistable_charged"),string("phi of charged stable particles"),string("#phi"),100, -3.15, 3.15);
  m_cmultpart=initHist(string("cmultpart"),string("number of charged particles in the event"),string("charged particle multiplicity"),
                       240, 0., 1200.);
  m_pttoppair=initHist(string("pttoppair"),string("transveral momentum of top pair"),string("p_{T} [GeV]"),200, 0., 1000.);
  m_pttoppair_log=initHist(string("pttoppair_logscale"),string("transveral momentum of top pair - logscale -"),string("p_{T} [GeV]"),200, 0., 1000.);
  m_W_W_phi=initHist(string("W_W_phi"),string("phi between W from top and W from tbar"),string("#phi_{W} - #phi_{Wbar}"),
                     48, -3.15, 3.15);
  m_top_tbar_phi=initHist(string("top_tbar_phi"),string("phi between top and tbar"),string("#phi_{top} - #phi_{tbar}"),48, -3.15, 3.15);
  m_W_top_phi=initHist(string("W_top_phi"),string("phi between W (from top) and top (also for tbar)"),
                       string("#phi_{top} - #phi_{W}"),48, -3.15, 3.15);
  m_Wpt=initHist(string("Wpt"),string("transversal momentum of W"),string("p_{T} [GeV]"),200, 0., 1000.);
  m_Wpt_log=initHist(string("Wpt_logscale"),string("transversal momentum of W - logscale -"),string("p_{T} [GeV]"),200, 0., 1000.);
  m_Weta=initHist(string("Weta"),string("eta of W"),string("#eta"),60, -6., 6.);
  m_Wphi=initHist(string("Wphi"),string("phi of W"),string("#phi"),48, -3.15, 3.15);
  
  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.);
  iSTOP = 0;
  return true;
}

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

  // Some Variables for storage are needed
  int top_barcode=0;
  int tbar_barcode=0;
  int W_barcode=0;
  int Wbar_barcode=0;
  
  //counter variable
  int cmult=0;

  //variable for status of particle
  //int properStatus = 2;

  // Fill the eventnumber
  m_evtnr->Fill(hepmcevt->event_number());

  // 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();

      // find the last top quarks in the event and 
      // store their barcodes  to find them back in the event
      if (pid == 6) top_barcode = (*p)->barcode();
      if (pid == -6)  tbar_barcode = (*p)->barcode();
      
      //the last W-bosons coming from top 
      //properStatus = 2;
      //if (abs(pid) == 6 && ( ((*p)->status())%1000 == properStatus) && (*p)->end_vertex())
      if (abs(pid) == 6 && (*p)->end_vertex())
        {
          // search for decaying particles
          // since the particles should decay, they should have a decay vertex
          if (!(*p)->end_vertex()) continue;
//        if (!((*p)->end_vertex()))
//          {
//            std::cout<<"CRASH!!!!!"<<std::endl;
//            std::cout<<"event number: "<<hepmcevt->event_number()<<std::endl;
//            std::cout<<"pdgId (of top): "<<pid<<std::endl;
//            std::cout<<"particle number (top barcode): "<<(*p)->barcode()<<std::endl;
//            std::cout<<"(*p)->status(): "<<(*p)->status()<<std::endl;
//            std::cout<<"(*p)->end_vertex(): "<<(*p)->end_vertex()<<std::endl;
              
//            std::cout<<"skip event"<<std::endl;
//            return false;
//          }
          
          HepMC::GenVertex::particle_iterator firstChild = (*p)->end_vertex()->particles_begin(HepMC::children);
          HepMC::GenVertex::particle_iterator endChild =  (*p)->end_vertex()->particles_end(HepMC::children);
          // Now check the decay products and search for a W boson or a Wbar boson
          for(HepMC::GenVertex::particle_iterator iter=firstChild ; iter!=endChild; iter++)
            {
              if( ((*iter)->pdg_id()) == 24 )
                {
                  W_barcode = (*iter)->barcode();
                }
              
              if( ((*iter)->pdg_id()) == -24 )
                {
                  Wbar_barcode = (*iter)->barcode();
                }
            }
        }
      // take all stable particles of this event
     if(!IsFinalStateParticle((*p))) continue;
      
      if(! (abs(pid)==12 || abs(pid)==14 || abs(pid)==16) )
        {
          m_ptstable->Fill( (*p)->momentum().perp() );
          m_ptstable_log->Fill( (*p)->momentum().perp() );
          m_etastable->Fill( (*p)->momentum().eta() );
          m_phistable->Fill( (*p)->momentum().phi() );
        }
      
      // cut out the neutral particles  (charge is not stored in HepMC)
      // just remove the neutrinos, gammas and neutrons from 
      // the final state
      if(!chargedParticle(pid)) continue;
      
      m_ptstable_charged->Fill( (*p)->momentum().perp() );
      m_ptstable_charged_log->Fill( (*p)->momentum().perp() );
      m_etastable_charged->Fill( (*p)->momentum().eta() );
      m_phistable_charged->Fill( (*p)->momentum().phi() );
      
      cmult++;
      
    } // end of loop over particles

  // 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<<"TopAnalysis::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());
  }
    
  // now we have the toppair - let's histogramm pt, eta and the pt sum
  HepMC::GenParticle *top = hepmcevt->barcode_to_particle(top_barcode);
  HepMC::GenParticle *tbar = hepmcevt->barcode_to_particle(tbar_barcode);
  if(top && tbar) 
    {
      m_toppt->Fill( top->momentum().perp() );
      m_toppt->Fill( tbar->momentum().perp() );
      m_toppt_log->Fill( top->momentum().perp() );
      m_toppt_log->Fill( tbar->momentum().perp() );
      m_topeta->Fill( top->momentum().eta() );
      m_topeta->Fill( tbar->momentum().eta() );
      m_topphi->Fill( top->momentum().phi() );
      m_topphi->Fill( tbar->momentum().phi() );
      
      // calculate the pt of the ttbar system 
      double px = top->momentum().px();
      px += tbar->momentum().px();
      double py = top->momentum().py();
      py += tbar->momentum().py();
      double perp = px*px + py*py;
      perp = TMath::Sqrt(perp);
      m_pttoppair->Fill(perp);
      m_pttoppair_log->Fill(perp);
      m_cmultpart->Fill(cmult);

      //calculate phi angle between top and antitop
      double top_phi = 0;
      double tbar_phi = 0;
      double DeltaPhi_top_tbar = 0;

      top_phi = top->momentum().phi();
      tbar_phi = tbar->momentum().phi();
      DeltaPhi_top_tbar = top_phi - tbar_phi;

      //if( fabs(DeltaPhi_top_tbar) > TMath::Pi() ||  fabs(DeltaPhi_top_tbar) == TMath::Pi() )
      if( DeltaPhi_top_tbar > TMath::Pi() ||  DeltaPhi_top_tbar == TMath::Pi() )
        {
          DeltaPhi_top_tbar = DeltaPhi_top_tbar - 2 * TMath::Pi() ;
        }

      if( DeltaPhi_top_tbar < - TMath::Pi() ||  DeltaPhi_top_tbar == - TMath::Pi() )
        {
          DeltaPhi_top_tbar = DeltaPhi_top_tbar + 2 * TMath::Pi();
        }

      m_top_tbar_phi->Fill(DeltaPhi_top_tbar);
    }

   //for W
   HepMC::GenParticle *W = hepmcevt->barcode_to_particle(W_barcode);
   HepMC::GenParticle *Wbar = hepmcevt->barcode_to_particle(Wbar_barcode);

   if(W && Wbar) 
     {
       m_Wpt->Fill( W->momentum().perp() );
       m_Wpt->Fill( Wbar->momentum().perp() );
       m_Wpt_log->Fill( W->momentum().perp() );
       m_Wpt_log->Fill( Wbar->momentum().perp() );
       m_Weta->Fill( W->momentum().eta() );
       m_Weta->Fill( Wbar->momentum().eta() );
       m_Wphi->Fill( W->momentum().phi() );
       m_Wphi->Fill( Wbar->momentum().phi() );
       
       //calculate phi angle between W and Wbar
       double W_phi = 0;
       double Wbar_phi = 0;
       double DeltaPhi_W_W = 0;
       W_phi = W->momentum().phi();
       Wbar_phi = Wbar->momentum().phi();
       DeltaPhi_W_W = W_phi - Wbar_phi;
       
       if( DeltaPhi_W_W > TMath::Pi() ||  DeltaPhi_W_W == TMath::Pi() )
         {
           DeltaPhi_W_W = DeltaPhi_W_W - 2 * TMath::Pi() ;
         }
       
       if( DeltaPhi_W_W < - TMath::Pi() ||  DeltaPhi_W_W == - TMath::Pi() )
         {
           DeltaPhi_W_W = DeltaPhi_W_W + 2 * TMath::Pi();
         }
       m_W_W_phi->Fill(DeltaPhi_W_W);
     }

   //calculating the angle between top and W
   
   //if (top && ( top->status()%1000 == 2 && top->end_vertex()) )
   if (top && top->end_vertex())
     {
       //search for decaying tops
       // since the tops should decay, they should have a decay vertex
       if (!top->end_vertex()) return false;

       double top_phi = 0;   
       double W_phi = 0;
       double DeltaPhi_W_top = 0;
       
       HepMC::GenVertex::particle_iterator firstChild = top->end_vertex()->particles_begin(HepMC::children);
       HepMC::GenVertex::particle_iterator endChild = top->end_vertex()->particles_end(HepMC::children);

       //loop over all children of top
       for(HepMC::GenVertex::particle_iterator iter=firstChild ; iter!=endChild; iter++)
         {
           if( abs((*iter)->pdg_id()) == 24 )
             {
               top_phi = top->momentum().phi();
               W_phi = (*iter)->momentum().phi();
               DeltaPhi_W_top = top_phi - W_phi;
               
               if( DeltaPhi_W_top > TMath::Pi() ||  DeltaPhi_W_top == TMath::Pi() )
                 {
                   DeltaPhi_W_top = DeltaPhi_W_top - 2 * TMath::Pi() ;
                 }
               
               if( DeltaPhi_W_top < - TMath::Pi() ||  DeltaPhi_W_top == - TMath::Pi() )
                 {
                   DeltaPhi_W_top = DeltaPhi_W_top + 2 * TMath::Pi();
                 }
               m_W_top_phi->Fill(DeltaPhi_W_top);
             }
         }
     }

   //if (tbar && ( tbar->status()%1000 == 2 && tbar->end_vertex()) )
   if (tbar && tbar->end_vertex())
     {
       //search for decaing tbars
       // since the tops should decay, they should have a decay vertex
       if (!tbar->end_vertex()) return false;

       double top_phi = 0;
       double W_phi = 0;
       double DeltaPhi_W_top = 0;
       
       HepMC::GenVertex::particle_iterator firstChild = tbar->end_vertex()->particles_begin(HepMC::children);
       HepMC::GenVertex::particle_iterator endChild = tbar->end_vertex()->particles_end(HepMC::children);
       //loop over all children of tbar
       for(HepMC::GenVertex::particle_iterator iter=firstChild ; iter!=endChild; iter++)
         {
           if( abs((*iter)->pdg_id()) == 24 )
             {
               top_phi = tbar->momentum().phi();
               W_phi = (*iter)->momentum().phi();
               DeltaPhi_W_top = top_phi - W_phi;
               
               if( DeltaPhi_W_top > TMath::Pi() ||  DeltaPhi_W_top == TMath::Pi() )
                 {
                   DeltaPhi_W_top = DeltaPhi_W_top - 2 * TMath::Pi() ;
                 }
               
               if( DeltaPhi_W_top < - TMath::Pi() ||  DeltaPhi_W_top == - TMath::Pi() )
                 {
                   DeltaPhi_W_top = DeltaPhi_W_top + 2 * TMath::Pi();
                 }
               m_W_top_phi->Fill(DeltaPhi_W_top);
             }
         }
     }
   
   return true;
} //end Process

---