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"

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

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

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

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"),1000, 0., 1000.);
  m_taupt=initHist(string("taupt"),string("transversal momentum of tau+ and tau-"),string("p_{T} [GeV]"),100, 0., 200.);
  m_taupt_log=initHist(string("taupt_logscale"),string("transversal momentum of tau+ and tau- - logscale -"),
                       string("p_{T} [GeV]"),100, 0., 200.);
  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]"),200, 0., 50.);
  m_ptstable_charged_log=initHist(string("ptstable_charged_logscale"),
                                  string("transveral momentum of charged stable particles - logscale -"),
                                  string("p_{T} [GeV]"),200, 0., 50.);
  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]"),100, 0., 200.);
  m_E_stable_charged_log=initHist(string("E_stable_without_neutrinos_logscale"),
                                  string("energy of charged stable particles - logscale -"),string("E [GeV]"),100, 0., 200.);
  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]"),
                          200, 0., 100.);
  m_pt_hightrack_log=initHist(string("pt_hightrack_logscale"),string("transversal momentum of highest p_{T}-track - logscale -"),
                              string("p_{T} [GeV]"),200, 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("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.);

  return true;
}

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

      //if (pid == 23 && (*p)->status()%1000==2 && (*p)->end_vertex())
      if (pid == 23 && (*p)->end_vertex())
        {
          // search for decaying particles
          // since the particles should decay, they should have a decay vertex
          if(!((*p)->end_vertex())) continue;
//          {
//            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 false;
//          }
      
          // find the last taus in the event and 
          // store their barcodes  to find them back in the event
          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
     if(!IsFinalStateParticle((*p))) 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)) 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);
  
  // check if jets were found
  //if(!m_inclusive_jets.size()) std::cout<<"TauAnalysis::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 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-
  //properStatus = 2;
  
  //if (tau && ( tau->status()%1000 == properStatus && tau->end_vertex()) )
  if (tau && tau->end_vertex())
    {
      // search for decaying taus
      // since the taus should decay, they should have a decay vertex
      if (!tau->end_vertex()) return false;

      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 false;
//          }
          
          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())
                //if(chargedParticle(iChild_pid)==true && !(*iChild)->end_vertex())
                {
                  // search for stable decay products of the taus
                  // since the children of the taus should not decay, they should not have a decay vertex
                  if ((*iChild)->end_vertex()) continue;

                  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) 
      //{
      //std::cout<<"nTrack for tau decay is "<<nTrack<<" for event "<<hepmcevt->event_number()<<std::endl;
      
          m_nTrack_tau->Fill(nTrack);
          m_pt_hightrack->Fill(highpT_Track);
          m_pt_hightrack_log->Fill(highpT_Track);
          //}
    }
  
  //tau+ (tbar)
  //properStatus = 2;
  
  //if (taubar && taubar->status()%1000 == properStatus && taubar->end_vertex() )
  if (taubar && taubar->end_vertex())
    {
      // search for decaying taubars
      // since the antitaus should decay, they should have a decay vertex
      if (!taubar->end_vertex()) return false;

      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 false;
//      }

      //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)
                //if(chargedParticle(iChild_pid)==true && !(*iChild)->end_vertex())
                {
                  //search for stable decay products of the taubars
                  // since the children of the antitaus should not decay, they should not have a decay vertex
                  if ((*iChild)->end_vertex()) continue;

                  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 true;
} //end of Process

---