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
00015 #include "TH1.h"
00016 #include "TH2.h"
00017 #include "TFile.h"
00018 #include "TMath.h"
00019 #include "TLorentzVector.h"
00020
00021
00022 #include "../include/ttbarAnalysis.h"
00023
00024
00025
00026
00027 ttbarAnalysis::ttbarAnalysis()
00028 {
00029 }
00030
00031
00032 ttbarAnalysis::~ttbarAnalysis()
00033 {
00034 }
00035
00036
00037 int ttbarAnalysis::Init(double tr_max_eta, double tr_min_pt)
00038 {
00039
00040 m_max_eta=tr_max_eta;
00041
00042
00043 m_min_pt=tr_min_pt;
00044
00045
00046 m_outputFileName="ttbarAnalysis.root";
00047 m_outputRootDir="Ttbar";
00048
00049
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
00094 int top_barcode=0;
00095 int tbar_barcode=0;
00096 int W_barcode=0;
00097 int Wbar_barcode=0;
00098
00099
00100 int cmult=0;
00101
00102
00103
00104
00105
00106 m_evtnr->Fill(hepmcevt->event_number());
00107
00108
00109 for ( HepMC::GenEvent::particle_const_iterator p =
00110 hepmcevt->particles_begin(); p != hepmcevt->particles_end(); ++p )
00111 {
00112
00113 int pid = (*p)->pdg_id();
00114
00115
00116
00117 if (pid == 6) top_barcode = (*p)->barcode();
00118 if (pid == -6) tbar_barcode = (*p)->barcode();
00119
00120
00121
00122
00123 if (abs(pid) == 6 && (*p)->end_vertex())
00124 {
00125
00126
00127 if (!(*p)->end_vertex()) continue;
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
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
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
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
00169
00170
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 }
00181
00182
00183 CLHEP::HepLorentzVector lv_leadingJet(0,0,0,0);
00184
00185
00186
00187
00188
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
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
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
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
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
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
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
00280
00281
00282 if (top && top->end_vertex())
00283 {
00284
00285
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
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
00319 if (tbar && tbar->end_vertex())
00320 {
00321
00322
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
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 }