ttbarAnalysis.cc

Go to the documentation of this file.
00001 #include <iostream>
00002 #include <sstream>
00003 #include <stdio.h>
00004 #include "HepMC/GenEvent.h"
00005 #include "HepMC/IO_GenEvent.h"
00006 #include "HepMC/GenParticle.h"
00007 #include "HepMC/GenVertex.h"
00008 #include "HepMC/IO_AsciiParticles.h"
00009 #include "HepMC/SimpleVector.h"
00010 #include "CLHEP/Vector/LorentzVector.h"
00011 
00012 using namespace std;
00013 
00014 // ROOT headers
00015 #include "TH1.h"
00016 #include "TH2.h"
00017 #include "TFile.h"
00018 #include "TMath.h"
00019 #include "TLorentzVector.h"
00020 
00021 //ttbar analysis header
00022 #include "../include/ttbarAnalysis.h"
00023 
00024 //**********************
00025 
00026 // empty default constructor
00027 ttbarAnalysis::ttbarAnalysis()
00028 {
00029 }
00030 
00031 /// empty default destructor
00032 ttbarAnalysis::~ttbarAnalysis()
00033 {
00034 }
00035 
00036 
00037 int ttbarAnalysis::Init(double tr_max_eta, double tr_min_pt)
00038 {
00039   // specify default eta cut
00040   m_max_eta=tr_max_eta;
00041   
00042   // specify default pt cut
00043   m_min_pt=tr_min_pt;
00044 
00045   // default Output file name
00046   m_outputFileName="ttbarAnalysis.root";
00047   m_outputRootDir="Ttbar";
00048   
00049   //declaration of histograms
00050   m_evtnr=initHist(string("evtnr"),string("Event number"),string("Eventnumber"),1000, 0., 1000.);
00051   m_toppt=initHist(string("toppt"),string("transversal momentum of top and antitop"),string("p_{T} [GeV]"),300, 0., 1500.);
00052   m_toppt_log=initHist(string("toppt_logscale"),string("transversal momentum of top and antitop - logscale -"),string("p_{T} [GeV]"),300, 0., 1500.);
00053   m_topeta=initHist(string("topeta"),string("eta of top and antitop"),string("#eta"),60, -6., 6.);
00054   m_topphi=initHist(string("topphi"),string("phi of top and antitop"),string("#phi"),48, -3.15, 3.15);
00055   m_ptstable=initHist(string("ptstable"),string("transveral momentum of stable particles (without neutrinos)"),
00056                       string("p_{T} [GeV]"),400, 0., 100.);
00057   m_ptstable_log=initHist(string("ptstable_logscale"),
00058                           string("transveral momentum of stable particles (without neutrinos) - logscale -"), 
00059                           string("p_{T} [GeV]"),400, 0., 100.);
00060   m_etastable=initHist(string("etastable"),string("eta of stable particles (without neutrinos)"),string("#eta"),100, -6., 6.);
00061   m_phistable=initHist(string("phistable"),string("phi of stable particles (without neutrinos)"),string("#phi"),100, -3.15, 3.15);
00062   m_ptstable_charged=initHist(string("ptstable_charged"),string("transveral momentum of charged stable particles"),
00063                               string("p_{T} [GeV]"),400, 0., 100.);
00064   m_ptstable_charged_log=initHist(string("ptstable_charged_logscale"),
00065                                   string("transveral momentum of charged stable particles - logscale -"),
00066                                   string("p_{T} [GeV]"),400, 0., 100.);
00067   m_etastable_charged=initHist(string("etastable_charged"),string("eta of charged stable particles"),string("#eta"),100, -6., 6.);
00068   m_phistable_charged=initHist(string("phistable_charged"),string("phi of charged stable particles"),string("#phi"),100, -3.15, 3.15);
00069   m_cmultpart=initHist(string("cmultpart"),string("number of charged particles in the event"),string("charged particle multiplicity"),
00070                        240, 0., 1200.);
00071   m_pttoppair=initHist(string("pttoppair"),string("transveral momentum of top pair"),string("p_{T} [GeV]"),200, 0., 1000.);
00072   m_pttoppair_log=initHist(string("pttoppair_logscale"),string("transveral momentum of top pair - logscale -"),string("p_{T} [GeV]"),200, 0., 1000.);
00073   m_W_W_phi=initHist(string("W_W_phi"),string("phi between W from top and W from tbar"),string("#phi_{W} - #phi_{Wbar}"),
00074                      48, -3.15, 3.15);
00075   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);
00076   m_W_top_phi=initHist(string("W_top_phi"),string("phi between W (from top) and top (also for tbar)"),
00077                        string("#phi_{top} - #phi_{W}"),48, -3.15, 3.15);
00078   m_Wpt=initHist(string("Wpt"),string("transversal momentum of W"),string("p_{T} [GeV]"),200, 0., 1000.);
00079   m_Wpt_log=initHist(string("Wpt_logscale"),string("transversal momentum of W - logscale -"),string("p_{T} [GeV]"),200, 0., 1000.);
00080   m_Weta=initHist(string("Weta"),string("eta of W"),string("#eta"),60, -6., 6.);
00081   m_Wphi=initHist(string("Wphi"),string("phi of W"),string("#phi"),48, -3.15, 3.15);
00082   
00083   m_jet_count=baseAnalysis::initHist(string("jet_count"), string("Nr of Jets"), string("Nr Jets"),16, -0.5, 15.5);
00084   m_jet_pt=baseAnalysis::initHist(string("leadingjet_pt_high"),string("Leading Jet Transverse Momentum"),string("Jet p_{T} [GeV]"), 200, 0.,1000.);
00085   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.);
00086   iSTOP = 0;
00087   return true;
00088 }
00089 
00090 int ttbarAnalysis::Process(HepMC::GenEvent* hepmcevt) 
00091 {
00092 
00093   // Some Variables for storage are needed
00094   int top_barcode=0;
00095   int tbar_barcode=0;
00096   int W_barcode=0;
00097   int Wbar_barcode=0;
00098   
00099   //counter variable
00100   int cmult=0;
00101 
00102   //variable for status of particle
00103   //int properStatus = 2;
00104 
00105   // Fill the eventnumber
00106   m_evtnr->Fill(hepmcevt->event_number());
00107 
00108   // loop over the particles and select pdgid and pt
00109   for ( HepMC::GenEvent::particle_const_iterator p = 
00110           hepmcevt->particles_begin(); p != hepmcevt->particles_end(); ++p )
00111     {
00112       // use the pdg id to identify the particles
00113       int pid = (*p)->pdg_id();
00114 
00115       // find the last top quarks in the event and 
00116       // store their barcodes  to find them back in the event
00117       if (pid == 6) top_barcode = (*p)->barcode();
00118       if (pid == -6)  tbar_barcode = (*p)->barcode();
00119       
00120       //the last W-bosons coming from top 
00121       //properStatus = 2;
00122       //if (abs(pid) == 6 && ( ((*p)->status())%1000 == properStatus) && (*p)->end_vertex())
00123       if (abs(pid) == 6 && (*p)->end_vertex())
00124         {
00125           // search for decaying particles
00126           // since the particles should decay, they should have a decay vertex
00127           if (!(*p)->end_vertex()) continue;
00128 //        if (!((*p)->end_vertex()))
00129 //          {
00130 //            std::cout<<"CRASH!!!!!"<<std::endl;
00131 //            std::cout<<"event number: "<<hepmcevt->event_number()<<std::endl;
00132 //            std::cout<<"pdgId (of top): "<<pid<<std::endl;
00133 //            std::cout<<"particle number (top barcode): "<<(*p)->barcode()<<std::endl;
00134 //            std::cout<<"(*p)->status(): "<<(*p)->status()<<std::endl;
00135 //            std::cout<<"(*p)->end_vertex(): "<<(*p)->end_vertex()<<std::endl;
00136               
00137 //            std::cout<<"skip event"<<std::endl;
00138 //            return false;
00139 //          }
00140           
00141           HepMC::GenVertex::particle_iterator firstChild = (*p)->end_vertex()->particles_begin(HepMC::children);
00142           HepMC::GenVertex::particle_iterator endChild =  (*p)->end_vertex()->particles_end(HepMC::children);
00143           // Now check the decay products and search for a W boson or a Wbar boson
00144           for(HepMC::GenVertex::particle_iterator iter=firstChild ; iter!=endChild; iter++)
00145             {
00146               if( ((*iter)->pdg_id()) == 24 )
00147                 {
00148                   W_barcode = (*iter)->barcode();
00149                 }
00150               
00151               if( ((*iter)->pdg_id()) == -24 )
00152                 {
00153                   Wbar_barcode = (*iter)->barcode();
00154                 }
00155             }
00156         }
00157       // take all stable particles of this event
00158      if(!IsFinalStateParticle((*p))) continue;
00159       
00160       if(! (abs(pid)==12 || abs(pid)==14 || abs(pid)==16) )
00161         {
00162           m_ptstable->Fill( (*p)->momentum().perp() );
00163           m_ptstable_log->Fill( (*p)->momentum().perp() );
00164           m_etastable->Fill( (*p)->momentum().eta() );
00165           m_phistable->Fill( (*p)->momentum().phi() );
00166         }
00167       
00168       // cut out the neutral particles  (charge is not stored in HepMC)
00169       // just remove the neutrinos, gammas and neutrons from 
00170       // the final state
00171       if(!chargedParticle(pid)) continue;
00172       
00173       m_ptstable_charged->Fill( (*p)->momentum().perp() );
00174       m_ptstable_charged_log->Fill( (*p)->momentum().perp() );
00175       m_etastable_charged->Fill( (*p)->momentum().eta() );
00176       m_phistable_charged->Fill( (*p)->momentum().phi() );
00177       
00178       cmult++;
00179       
00180     } // end of loop over particles
00181 
00182   // storage of input particles for jet
00183   CLHEP::HepLorentzVector lv_leadingJet(0,0,0,0);
00184   
00185   // check if jets were found
00186   //if(!m_inclusive_jets.size()) std::cout<<"ttbarAnalysis::no jets found"<<std::endl;
00187 
00188   //Fill the histograms
00189   m_jet_count->Fill(m_inclusive_jets.size());
00190 
00191   if(m_inclusive_jets.size()){
00192     m_jet_pt->Fill(m_inclusive_jets[0].perp());
00193     m_jet_pt_log->Fill(m_inclusive_jets[0].perp());
00194   }
00195     
00196   // now we have the toppair - let's histogramm pt, eta and the pt sum
00197   HepMC::GenParticle *top = hepmcevt->barcode_to_particle(top_barcode);
00198   HepMC::GenParticle *tbar = hepmcevt->barcode_to_particle(tbar_barcode);
00199   if(top && tbar) 
00200     {
00201       m_toppt->Fill( top->momentum().perp() );
00202       m_toppt->Fill( tbar->momentum().perp() );
00203       m_toppt_log->Fill( top->momentum().perp() );
00204       m_toppt_log->Fill( tbar->momentum().perp() );
00205       m_topeta->Fill( top->momentum().eta() );
00206       m_topeta->Fill( tbar->momentum().eta() );
00207       m_topphi->Fill( top->momentum().phi() );
00208       m_topphi->Fill( tbar->momentum().phi() );
00209       
00210       // calculate the pt of the ttbar system 
00211       double px = top->momentum().px();
00212       px += tbar->momentum().px();
00213       double py = top->momentum().py();
00214       py += tbar->momentum().py();
00215       double perp = px*px + py*py;
00216       perp = TMath::Sqrt(perp);
00217       m_pttoppair->Fill(perp);
00218       m_pttoppair_log->Fill(perp);
00219       m_cmultpart->Fill(cmult);
00220 
00221       //calculate phi angle between top and antitop
00222       double top_phi = 0;
00223       double tbar_phi = 0;
00224       double DeltaPhi_top_tbar = 0;
00225 
00226       top_phi = top->momentum().phi();
00227       tbar_phi = tbar->momentum().phi();
00228       DeltaPhi_top_tbar = top_phi - tbar_phi;
00229 
00230       //if( fabs(DeltaPhi_top_tbar) > TMath::Pi() ||  fabs(DeltaPhi_top_tbar) == TMath::Pi() )
00231       if( DeltaPhi_top_tbar > TMath::Pi() ||  DeltaPhi_top_tbar == TMath::Pi() )
00232         {
00233           DeltaPhi_top_tbar = DeltaPhi_top_tbar - 2 * TMath::Pi() ;
00234         }
00235 
00236       if( DeltaPhi_top_tbar < - TMath::Pi() ||  DeltaPhi_top_tbar == - TMath::Pi() )
00237         {
00238           DeltaPhi_top_tbar = DeltaPhi_top_tbar + 2 * TMath::Pi();
00239         }
00240 
00241       m_top_tbar_phi->Fill(DeltaPhi_top_tbar);
00242     }
00243 
00244    //for W
00245    HepMC::GenParticle *W = hepmcevt->barcode_to_particle(W_barcode);
00246    HepMC::GenParticle *Wbar = hepmcevt->barcode_to_particle(Wbar_barcode);
00247 
00248    if(W && Wbar) 
00249      {
00250        m_Wpt->Fill( W->momentum().perp() );
00251        m_Wpt->Fill( Wbar->momentum().perp() );
00252        m_Wpt_log->Fill( W->momentum().perp() );
00253        m_Wpt_log->Fill( Wbar->momentum().perp() );
00254        m_Weta->Fill( W->momentum().eta() );
00255        m_Weta->Fill( Wbar->momentum().eta() );
00256        m_Wphi->Fill( W->momentum().phi() );
00257        m_Wphi->Fill( Wbar->momentum().phi() );
00258        
00259        //calculate phi angle between W and Wbar
00260        double W_phi = 0;
00261        double Wbar_phi = 0;
00262        double DeltaPhi_W_W = 0;
00263        W_phi = W->momentum().phi();
00264        Wbar_phi = Wbar->momentum().phi();
00265        DeltaPhi_W_W = W_phi - Wbar_phi;
00266        
00267        if( DeltaPhi_W_W > TMath::Pi() ||  DeltaPhi_W_W == TMath::Pi() )
00268          {
00269            DeltaPhi_W_W = DeltaPhi_W_W - 2 * TMath::Pi() ;
00270          }
00271        
00272        if( DeltaPhi_W_W < - TMath::Pi() ||  DeltaPhi_W_W == - TMath::Pi() )
00273          {
00274            DeltaPhi_W_W = DeltaPhi_W_W + 2 * TMath::Pi();
00275          }
00276        m_W_W_phi->Fill(DeltaPhi_W_W);
00277      }
00278 
00279    //calculating the angle between top and W
00280    
00281    //if (top && ( top->status()%1000 == 2 && top->end_vertex()) )
00282    if (top && top->end_vertex())
00283      {
00284        //search for decaying tops
00285        // since the tops should decay, they should have a decay vertex
00286        if (!top->end_vertex()) return false;
00287 
00288        double top_phi = 0;   
00289        double W_phi = 0;
00290        double DeltaPhi_W_top = 0;
00291        
00292        HepMC::GenVertex::particle_iterator firstChild = top->end_vertex()->particles_begin(HepMC::children);
00293        HepMC::GenVertex::particle_iterator endChild = top->end_vertex()->particles_end(HepMC::children);
00294 
00295        //loop over all children of top
00296        for(HepMC::GenVertex::particle_iterator iter=firstChild ; iter!=endChild; iter++)
00297          {
00298            if( abs((*iter)->pdg_id()) == 24 )
00299              {
00300                top_phi = top->momentum().phi();
00301                W_phi = (*iter)->momentum().phi();
00302                DeltaPhi_W_top = top_phi - W_phi;
00303                
00304                if( DeltaPhi_W_top > TMath::Pi() ||  DeltaPhi_W_top == TMath::Pi() )
00305                  {
00306                    DeltaPhi_W_top = DeltaPhi_W_top - 2 * TMath::Pi() ;
00307                  }
00308                
00309                if( DeltaPhi_W_top < - TMath::Pi() ||  DeltaPhi_W_top == - TMath::Pi() )
00310                  {
00311                    DeltaPhi_W_top = DeltaPhi_W_top + 2 * TMath::Pi();
00312                  }
00313                m_W_top_phi->Fill(DeltaPhi_W_top);
00314              }
00315          }
00316      }
00317 
00318    //if (tbar && ( tbar->status()%1000 == 2 && tbar->end_vertex()) )
00319    if (tbar && tbar->end_vertex())
00320      {
00321        //search for decaing tbars
00322        // since the tops should decay, they should have a decay vertex
00323        if (!tbar->end_vertex()) return false;
00324 
00325        double top_phi = 0;
00326        double W_phi = 0;
00327        double DeltaPhi_W_top = 0;
00328        
00329        HepMC::GenVertex::particle_iterator firstChild = tbar->end_vertex()->particles_begin(HepMC::children);
00330        HepMC::GenVertex::particle_iterator endChild = tbar->end_vertex()->particles_end(HepMC::children);
00331        //loop over all children of tbar
00332        for(HepMC::GenVertex::particle_iterator iter=firstChild ; iter!=endChild; iter++)
00333          {
00334            if( abs((*iter)->pdg_id()) == 24 )
00335              {
00336                top_phi = tbar->momentum().phi();
00337                W_phi = (*iter)->momentum().phi();
00338                DeltaPhi_W_top = top_phi - W_phi;
00339                
00340                if( DeltaPhi_W_top > TMath::Pi() ||  DeltaPhi_W_top == TMath::Pi() )
00341                  {
00342                    DeltaPhi_W_top = DeltaPhi_W_top - 2 * TMath::Pi() ;
00343                  }
00344                
00345                if( DeltaPhi_W_top < - TMath::Pi() ||  DeltaPhi_W_top == - TMath::Pi() )
00346                  {
00347                    DeltaPhi_W_top = DeltaPhi_W_top + 2 * TMath::Pi();
00348                  }
00349                m_W_top_phi->Fill(DeltaPhi_W_top);
00350              }
00351          }
00352      }
00353    
00354    return true;
00355 } //end Process

Generated on Wed Aug 31 09:44:48 2011 for HepMCAnalysis by  doxygen 1.4.7