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