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

//Tau analysis header
#include "../include/TauAnalysis.h"

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

//starting functions
TauAnalysis::TauAnalysis()
{

}

TauAnalysis::~TauAnalysis()
{
  while (!m_histVector.empty())
    {
      delete m_histVector.back();
      m_histVector.pop_back();
    }
}

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

  //declaration of histograms
  m_evtnr=initHist(string("evtnr"),string("Event number"),string("Eventnumber"),100, 0., 99.);
  m_taupt=initHist(string("taupt"),string("transversal momentum of tau+ and tau-"),string("p_{T} [GeV]"),50, 0., 100.);
  m_taupt_log=initHist(string("taupt_logscale"),string("transversal momentum of tau+ and tau- - logscale -"),
                       string("p_{T} [GeV]"),50, 0., 100.);
  m_taueta=initHist(string("taueta"),string("eta of tau+ and tau-"),string("#eta"),60, -6., 6.);
  m_tauphi=initHist(string("tauphi"),string("phi of tau+ and tau-"),string("#phi"),48, -3.15, 3.15);
  m_ptstable=initHist(string("ptstable"),string("transveral momentum of stable particles (without neutrinos)"),
                      string("p_{T} [GeV]"),60, 0., 15.);
  m_ptstable_log=initHist(string("ptstable_logscale"),string("transveral momentum of stable particles (without neutrinos) - logscale"),
                          string("p_{T} [GeV]"),60, 0., 15.);
  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_E_stable=initHist(string("E_stable"),string("energy of stable particles (without neutrinos)"),string("E [GeV]"),50, 0., 100.);
  m_E_stable_log=initHist(string("E_stable_logscale"),string("energy of stable particles (without neutrinos) - logscale -"),
                          string("E [GeV]"),50, 0., 100.);
  m_ptstable_charged=initHist(string("ptstable_charged"),string("transveral momentum of charged stable particles"),
                              string("p_{T} [GeV]"),60, 0., 15.);
  m_ptstable_charged_log=initHist(string("ptstable_charged_logscale"),
                                  string("transveral momentum of charged stable particles - logscale -"),
                                  string("p_{T} [GeV]"),60, 0., 15.);
  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_E_stable_charged=initHist(string("E_stable_without_neutrinos"),string("energy of charged stable particles"),
                              string("E [GeV]"),50, 0., 100.);
  m_E_stable_charged_log=initHist(string("E_stable_without_neutrinos_logscale"),
                                  string("energy of charged stable particles - logscale -"),string("E [GeV]"),50, 0., 100.);
  m_cmultpart=initHist(string("cmultpart"),string("number of charged particles in the event"),string("charged particle multiplicity"),
                       120, 0., 600.);
  m_nTrack_tau=initHist(string("nTrack_tau"),string("number of tracks of tau decay"),string("number of tracks"),16, -0.5, 15.5);
  m_pt_hightrack=initHist(string("pt_hightrack"),string("transversal momentum of highest p_{T}-track"),string("p_{T} [GeV]"),
                          100, 0., 50.);
  m_pt_hightrack_log=initHist(string("pt_hightrack_logscale"),string("transversal momentum of highest p_{T}-track - logscale -"),
                              string("p_{T} [GeV]"),100, 0., 50.);
  
  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.);
  
  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;
}

int TauAnalysis::Process(HepMC::GenEvent* hepmcevt) 
{
  
  int cmult = 0; //charged particle multiplicity
  int tau_barcode = 0;
  int taubar_barcode = 0;
  
  //variable for status of particle
  int properStatus = 0;
  
  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 taus in the event and 
      // store their barcodes  to find them back in the event
      
      if (pid == 23)
        {
          if (!((*p)->end_vertex()))
            {
              std::cout<<"CRASH!!!!!"<<std::endl;
              std::cout<<"event number: "<<hepmcevt->event_number()<<std::endl;
              std::cout<<"pdgId (of Z): "<<pid<<std::endl;
              std::cout<<"particle number (Z 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 FAILURE;
            }
          for(HepMC::GenVertex::particle_iterator iChild=(*p)->end_vertex()->particles_begin(HepMC::descendants); 
              iChild!=(*p)->end_vertex()->particles_end(HepMC::descendants); iChild++)
            {
              if( (*iChild)->pdg_id() != pid )
                {
                  if ((*iChild)->pdg_id() == 15) tau_barcode = (*iChild)->barcode();
                  if ((*iChild)->pdg_id() == -15)  taubar_barcode = (*iChild)->barcode();
                }
            }
        }
      
      // take all stable particles of this event
      properStatus = 1;
      
      if( !(*p)->status() == properStatus) continue;
      
      if(! ( abs(pid)==12 || abs(pid)==14 || abs(pid)==16 ) ) //energy stable particles without neutrinos
        {
          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() );
          m_E_stable->Fill( (*p)->momentum().e() );
          m_E_stable_log->Fill( (*p)->momentum().e() );
        }
      
      // 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)==NOTCHARGED) 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() );
      m_E_stable_charged->Fill( (*p)->momentum().e() );
      m_E_stable_charged_log->Fill( (*p)->momentum().e() );
      
      cmult++;
    }  // end of loop over particles
  
  // 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());
  }
  
  // now we have the taupair - let's histogramm pt, eta and the pt sum
  HepMC::GenParticle *tau = hepmcevt->barcode_to_particle(tau_barcode);
  HepMC::GenParticle *taubar = hepmcevt->barcode_to_particle(taubar_barcode);
  if(tau && taubar) 
    {
      m_taupt->Fill( tau->momentum().perp() );
      m_taupt->Fill( taubar->momentum().perp() );
      m_taupt_log->Fill( tau->momentum().perp() );
      m_taupt_log->Fill( taubar->momentum().perp() );
      m_taueta->Fill( tau->momentum().eta() );
      m_taueta->Fill( taubar->momentum().eta() );
      m_tauphi->Fill( tau->momentum().phi() );
      m_tauphi->Fill( taubar->momentum().phi() );
      
      m_cmultpart->Fill(cmult);
    }
  
  //checking for children of tau, needed for nTrack and the highest pT track
  
  //tau-
  
  if (tau && tau->end_vertex())
    {
      int nTrack = 0; //number of tracks in tau decay
      double highpT_Track = 0.; //track with hightest pT coming from tau decay
      
      //important for the tracks (in the following the commented lines apply to the first generation of tau decay:
      //if you look only for the first generation as descendant (child), than the child can be stable or not stable,
      //only if you look for more than one generation as descendants, you need the proper status on stable;
      //this means that the final decay product of the tau has to be stable
      
      //only first generation
      //for(HepMC::GenVertex::particle_iterator iChild=tau->end_vertex()->particles_begin(HepMC::children); 
      //  iChild!=tau->end_vertex()->particles_end(HepMC::children); iChild++)
      
      //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 (!tau->end_vertex())
            {
              std::cout<<"CRASH!!!!!"<<std::endl;
              std::cout<<"event number: "<<hepmcevt->event_number()<<std::endl;
              std::cout<<"particle number (tau barcode): "<<tau->barcode()<<std::endl;
              std::cout<<"tau->status(): "<<tau->status()<<std::endl;
              std::cout<<"tau->end_vertex(): "<<tau->end_vertex()<<std::endl;
              
              std::cout<<"skip event"<<std::endl;
              return FAILURE;
            }
          
          if( (*iChild)->pdg_id() != tau->pdg_id() )
            {
              int iChild_pid = (*iChild)->pdg_id();
              properStatus = 1;
              
              if( (*iChild)->status() == properStatus && chargedParticle(iChild_pid)==CHARGED)
                {
                  nTrack = nTrack + 1;
                  if ( (*iChild)->momentum().perp() > highpT_Track )
                    {
                      highpT_Track = (*iChild)->momentum().perp();
                    }
                }
            }
        }

      //only fill nTrack and pt of highest track if a tau decays
      if (nTrack > 0) 
        {
          m_nTrack_tau->Fill(nTrack);
          m_pt_hightrack->Fill(highpT_Track);
          m_pt_hightrack_log->Fill(highpT_Track);
        }
    }
  
  //tau+ (tbar)
  
  if (taubar && taubar->end_vertex())
    {
      int nTrack = 0; //number of tracks in tau decay
      double highpT_Track = 0.; //track with hightest pT coming from tau decay
      
      //important for the tracks (in the following the commented lines apply to the first generation of tau decay:
      //if you look only for the first generation as descendant (child), than the child can be stable or not stable,
      //only if you look for more than one generation as descendants, you need the proper status on stable;
      //this means that the final decay product of the tau has to be stable
      
      //only first generation
      //for(HepMC::GenVertex::particle_iterator iChild=taubar->end_vertex()->particles_begin(HepMC::children); 
      //  iChild!=taubar->end_vertex()->particles_end(HepMC::children); iChild++)
      
      if (!taubar->end_vertex())
        {
          std::cout<<"CRASH!!!!!"<<std::endl;
          std::cout<<"event number: "<<hepmcevt->event_number()<<std::endl;
          std::cout<<"particle number (tau barcode): "<<taubar->barcode()<<std::endl;
          std::cout<<"taubar->status(): "<<taubar->status()<<std::endl;
          std::cout<<"taubar->end_vertex(): "<<taubar->end_vertex()<<std::endl;
          
          std::cout<<"skip event"<<std::endl;
          return FAILURE;
        }

      //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() == properStatus && chargedParticle(iChild_pid)==CHARGED)
                {
                  nTrack = nTrack + 1;
                  if ( (*iChild)->momentum().perp() > highpT_Track )
                    {
                      highpT_Track = (*iChild)->momentum().perp();
                    }
                }
            }
          
        }
      //only fill nTrack and pt of highest track if a tau decays
      if (nTrack > 0) 
        {
          m_nTrack_tau->Fill(nTrack);
          m_pt_hightrack->Fill(highpT_Track);
          m_pt_hightrack_log->Fill(highpT_Track);
        }
    }
  return SUCCESS;
} //end of Process

---