HepMCAnalysis: src/WplusJetAnalysis.cc Source File

00001 #include <iostream>
00002 #include <stdio.h>
00003 #include "HepMC/GenEvent.h"
00004 #include "HepMC/IO_GenEvent.h"
00005 #include "HepMC/GenParticle.h"
00006 #include "HepMC/GenVertex.h"
00007 #include "HepMC/IO_AsciiParticles.h"
00008 #include "HepMC/SimpleVector.h"
00009 #include "HepPDT/ParticleData.hh"
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 #include "fastjet/PseudoJet.hh"
00022 #include "fastjet/ClusterSequence.hh"
00023 #include "fastjet/JetDefinition.hh"
00024 #include "fastjet/SISConePlugin.hh"
00025 
00026 //WplusJet analysis header
00027 #include "../include/WplusJetAnalysis.h"
00028 
00029 //**********************
00030 
00031 /**
00032 Constructor
00033 */
00034 WplusJetAnalysis::WplusJetAnalysis()
00035 {
00036 
00037 }
00038 
00039 /**
00040 Destructor: delete all the Pointer
00041 */
00042 WplusJetAnalysis::~WplusJetAnalysis()
00043 {
00044   while (!m_histVector.empty()){
00045     delete m_histVector.back();
00046     m_histVector.pop_back();
00047   }
00048 }
00049 
00050 /**
00051 Before using the class initialize
00052 */
00053 int WplusJetAnalysis::Init(double tr_max_eta, double tr_min_pt)
00054 {
00055   // specify default eta cut
00056   m_max_eta=tr_max_eta;
00057   
00058   //specify default pt cut
00059   m_min_pt=tr_min_pt;
00060   
00061   // default Output file name
00062   m_outputFileName="WplusJetAnalysis.root";
00063   m_outputRootDir="WplusJet";
00064   
00065   //declaration of histograms
00066   //m_evtnr=baseAnalysis::initHist(string("evtnr"),string("Event number"),string("Eventnumber"),100, 0., 100.);
00067   
00068   m_W_count=baseAnalysis::initHist(string("W_count"), string("Nr of W-Bosons"), string("Nr W-Boson"),16, -0.5, 15.5);
00069   m_W_trans_mass=baseAnalysis::initHist(string("W_trans_mass"),string("transvers W Mass"),string("m_{T,W} [GeV]"), 100, 0., 100.);
00070   m_W_trans_mass_log=baseAnalysis::initHist(string("W_trans_mass_logscale"),string("transvers W Mass -logscale-"),string("m_{T,W} [GeV]"), 100, 0., 100.);
00071   m_W_pt=baseAnalysis::initHist(string("W_pt"),string("W Transverse Momentum"),string("W p_{T} [GeV]"), 50, 0., 25.);
00072   m_W_pt_log=baseAnalysis::initHist(string("W_pt_logscale"),string("W Transverse Momentum -logscale-"),string("W p_{T} [GeV]"), 50, 0., 25.);
00073   m_W_pt_high=baseAnalysis::initHist(string("W_pt_high"),string("W Transverse Momentum"),string("W p_{T} [GeV]"), 50, 0., 150.);
00074   m_W_pt_high_log=baseAnalysis::initHist(string("W_pt_high_logscale"),string("W Transverse Momentum -logscale-"),string("W p_{T} [GeV]"), 50, 0., 150.);
00075   m_W_eta=baseAnalysis::initHist(string("W_eta"),string("#eta W-Boson"),string("#eta"), 60, -6., 6.);
00076   m_W_rapidity=baseAnalysis::initHist(string("W_rapidity"),string("rapidity y W-Boson"),string("y"), 60, -6., 6.);
00077   m_W_phi=baseAnalysis::initHist(string("W_phi"),string("#varphi W-Boson"),string("#varphi"), 64, -3.2, 3.2);
00078   
00079   m_lepton_count=baseAnalysis::initHist(string("lepton_count"), string("Nr of Leptons from W-Bosons"), string("Nr Leptons"),11, -0.5, 10.5);
00080   m_lepton_pt=baseAnalysis::initHist(string("lepton_pt_high"),string("Transverse Momentum Lepton from W-Boson"),string("lepton p_{T} [GeV]"), 50, 0., 25.);
00081   m_lepton_pt_log=baseAnalysis::initHist(string("lepton_pt_high_logscale"),string("Transverse Momentum Lepton from W-Boson -logscale-"),string("lepton p_{T} [GeV]"), 50, 0., 25.);
00082   m_lepton_pt_high=baseAnalysis::initHist(string("lepton_pt_high"),string("Transverse Momentum Lepton from W-Boson"),string("lepton p_{T} [GeV]"), 50, 0., 150.);
00083   m_lepton_pt_high_log=baseAnalysis::initHist(string("lepton_pt_high_logscale"),string("Transverse Momentum Lepton from W-Boson -logscale-"),string("lepton p_{T} [GeV]"), 50, 0., 150.);
00084   m_lepton_eta=baseAnalysis::initHist(string("lepton_eta"),string("#eta Lepton from W-Boson"),string("#eta"), 60, -6., 6.);
00085   m_lepton_phi=baseAnalysis::initHist(string("lepton_phi"),string("#varphi  Lepton from W-Boson"),string("#varphi"), 64, -3.2, 3.2);
00086   m_lepton_pdgid=baseAnalysis::initHist(string("lepton_pdgid"),string("PDG ID Lepton from W-Boson"),string("PDG ID"), 40, -20., 20.);
00087   
00088   m_charged_particle_multiplicity=baseAnalysis::initHist(string("charged_particle_multiplicity"),string("Number of charged Particles in the Event"),string("charged particle multiplicity"),100, 0., 300.);
00089   m_charged_particle_temp_pt= new TH1D("WplusJet_charged_particle_mean_pt","temp histogramm",200,0,200.);
00090   m_charged_particle_pt= new TH1D("WplusJet_charged_particle_pt","p_{T} of stable charged particles",200,0,200.);
00091   m_charged_particle_mean_pt=baseAnalysis::initHist(string("charged_particle_mean_pt"),string("Avarage p_{T} charged Particles in the Event"),string("charged particle #bar{p}_{T}"),50, 0., 10.);
00092   m_charged_particle_rms_pt=baseAnalysis::initHist(string("charged_particle_rms_pt"),string("RMS p_{T} charged Particles in the Event"),string("charged particle #sigma(p_{T})"),50, 0., 10.);
00093   //m_charged_particle_pdgID=baseAnalysis::initHist(string("charged_particle_pdgID"),string("charged Particles PDG ID"),string("charged particle pid"),400, -400., 400.);
00094   
00095   m_jet_count=baseAnalysis::initHist(string("jet_count"), string("Nr of Jets"), string("Nr Jets"),16, -0.5, 15.5);
00096   m_jet_pt=baseAnalysis::initHist(string("jet_pt_high"),string("Jet Transverse Momentum"),string("Jet p_{T} [GeV]"), 50, 0.,150.);
00097   m_jet_pt_log=baseAnalysis::initHist(string("jet_pt_high_logscale"),string("Jet Transverse Momentum -logscale-"),string("Jet p_{T} [GeV]"), 50, 0.,150.);
00098 
00099   m_jet_Deta_leadingJet_lepton=baseAnalysis::initHist(string("jet_Deta_leadingJet_lepton"),string("#Delta#eta leading Jet and Lepton from W-Boson"),string("#Delta#eta"), 60, 0., 6.);
00100   m_jet_Deta_leadingJet_lepton_singlejet=baseAnalysis::initHist(string("jet_Deta_leadingJet_lepton_singlejet"),string("#Delta#eta leading Jet (single) and Lepton from W-Boson"),string("#Delta#eta"), 60, 0., 6.);
00101   m_jet_Deta_leadingJet_lepton_multijet=baseAnalysis::initHist(string("jet_Deta_leadingJet_lepton_multijet"),string("#Delta#eta leading Jet (multi) and Lepton from W-Boson"),string("#Delta#eta"), 60, 0., 6.);
00102   
00103   m_jet_Dphi_leadingJet_lepton=baseAnalysis::initHist(string("jet_Dphi_leadingJet_lepton"),string("#Delta#varphi leading Jet and Lepton from W-Boson"),string("#Delta#varphi"), 64, 0., 3.2);
00104   m_jet_Dphi_leadingJet_lepton_singlejet=baseAnalysis::initHist(string("jet_Dphi_leadingJet_lepton_singlejet"),string("#Delta#varphi leading Jet (single) and Lepton from W-Boson"),string("#Delta#varphi"), 64, 0., 3.2);
00105   m_jet_Dphi_leadingJet_lepton_multijet=baseAnalysis::initHist(string("jet_Dphi_leadingJet_lepton_multijet"),string("#Delta#varphi leading Jet (multi) and Lepton from W-Boson"),string("#Delta#varphi"), 64, 0., 3.2);
00106   
00107   m_jet_DR_leadingJet_lepton=baseAnalysis::initHist(string("jet_DR_leadingJet_lepton"),string("#DeltaR leading Jet and Lepton from W-Boson"),string("#DeltaR"), 60, 0., 15.);
00108   m_jet_DR_leadingJet_lepton_singlejet=baseAnalysis::initHist(string("jet_DR_leadingJet_lepton_singlejet"),string("#DeltaR leading Jet (single) and Lepton from W-Boson"),string("#DeltaR"), 60, 0., 15.);
00109   m_jet_DR_leadingJet_lepton_multijet=baseAnalysis::initHist(string("jet_DR_leadingJet_lepton_multijet"),string("#DeltaR leading Jet (multi) and Lepton from W-Boson"),string("#DeltaR"), 60, 0., 15.);
00110   
00111   for(unsigned int hist=0; hist< m_histVector.size(); hist++)
00112     {
00113       std::ostringstream counter_s;
00114       counter_s << hist;
00115       std::string histogram_name=m_outputRootDir;
00116       histogram_name+="_";
00117       histogram_name+=counter_s.str();
00118       histogram_name+="_";
00119       histogram_name+=m_histVector[hist]->GetName();
00120       m_histVector[hist]->SetName(histogram_name.c_str());
00121     }
00122   
00123   return SUCCESS;
00124   
00125 }
00126 
00127 /**
00128 This is the main analysis function. 
00129 Search the particle, get their properties and the histograms.
00130 */
00131 int WplusJetAnalysis::Process(HepMC::GenEvent* hepmcevt) 
00132 {
00133   
00134   // Some Variables or storage classes are needed
00135   CLHEP::HepLorentzVector lv_leptonFromW(0,0,0,0);
00136   CLHEP::HepLorentzVector lv_neutrinoFromW(0,0,0,0);
00137   HepMC::GenParticle* leptonFromW=0;
00138   HepMC::GenParticle* neutrinoFromW=0;
00139   
00140   
00141   // Reset the histogramm
00142   m_charged_particle_temp_pt->Reset();
00143   
00144   // Fill the eventnumber
00145   //m_evtnr->Fill((hepmcevt->event_number())%100);
00146   
00147   // initialize the counter variables
00148   int Wcount=0;
00149   int leptoncount=0;
00150   int chargeParticlecount=0;
00151   
00152   // loop over the particles and select pdgid and pt
00153   for ( HepMC::GenEvent::particle_const_iterator p =  hepmcevt->particles_begin(); p != hepmcevt->particles_end(); ++p )
00154     {
00155       
00156       HepMC::GenParticle* WBoson=0;
00157       
00158       // use the pdg id to identify the particles
00159       int pid = (*p)->pdg_id();
00160       
00161       // W decays
00162       // look for W-Bosons (pdg = 24) 
00163       if( abs(pid) == 24 && (*p)->end_vertex() ){
00164 
00165         if (!(*p)->end_vertex()) continue;
00166 
00167         HepMC::GenVertex::particle_iterator firstChild = (*p)->end_vertex()->particles_begin(HepMC::children);
00168         HepMC::GenVertex::particle_iterator endChild =  (*p)->end_vertex()->particles_end(HepMC::children);
00169         // plot Pt, eta and phi of the W (generator values)
00170         if( (*firstChild)->pdg_id() != pid ){
00171           WBoson = (*p);
00172           m_W_pt->Fill((*p)->momentum().perp());
00173           m_W_pt_log->Fill((*p)->momentum().perp());
00174           m_W_pt_high->Fill((*p)->momentum().perp());
00175           m_W_pt_high_log->Fill((*p)->momentum().perp());
00176           m_W_eta->Fill((*p)->momentum().eta());
00177           m_W_rapidity->Fill(baseAnalysis::getRapidity(p));
00178           m_W_phi->Fill((*p)->momentum().phi());
00179           Wcount++;
00180         }
00181         
00182         bool WtoLepton=false;
00183         // Now check the decay products and search for a lepton and the neutrino
00184         for(HepMC::GenVertex::particle_iterator iter=firstChild ; iter!=endChild; iter++)
00185           {
00186             // set the properties of the neutrino Lorentzvector 
00187             if(abs((*iter)->pdg_id())==12 || abs((*iter)->pdg_id())==14 || abs((*iter)->pdg_id())==16){
00188               neutrinoFromW=(*iter);
00189               lv_neutrinoFromW.setPx((*iter)->momentum().px());
00190               lv_neutrinoFromW.setPy((*iter)->momentum().py());
00191               lv_neutrinoFromW.setPz((*iter)->momentum().pz());
00192               lv_neutrinoFromW.setE((*iter)->momentum().e());
00193             }
00194             // set the properties of the lepton Lorentzvector 
00195             if(abs((*iter)->pdg_id())==11 || abs((*iter)->pdg_id())==13 || abs((*iter)->pdg_id())==15) {
00196               leptonFromW=*iter;
00197               WtoLepton=true;
00198               m_lepton_pdgid->Fill((*iter)->pdg_id());
00199               lv_leptonFromW.setPx((*iter)->momentum().px());
00200               lv_leptonFromW.setPy((*iter)->momentum().py());
00201               lv_leptonFromW.setPz((*iter)->momentum().pz());
00202               lv_leptonFromW.setE((*iter)->momentum().e());
00203             }
00204           }
00205         
00206         if(!WtoLepton) continue;
00207 
00208         // fill the lepton histograms
00209         m_lepton_pt->Fill(lv_leptonFromW.perp());
00210         m_lepton_pt_log->Fill(lv_leptonFromW.perp());
00211         m_lepton_pt_high->Fill(lv_leptonFromW.perp());
00212         m_lepton_pt_high_log->Fill(lv_leptonFromW.perp());
00213         m_lepton_eta->Fill(lv_leptonFromW.eta());
00214         m_lepton_phi->Fill(lv_leptonFromW.phi());
00215         leptoncount++;
00216         
00217         // Use the standard Lorentzvector function to get the transverse Mass of the W
00218         // First set the z-component of the lepton and neutrino to zero and adjust the transverse energy (=transverse momentum)
00219         // then add the two lorentzvectors to get the transverse mass of the W-Boson
00220         CLHEP::HepLorentzVector lv_transverse_leptonFromW(lv_leptonFromW.px(),lv_leptonFromW.py(),0,lv_leptonFromW.perp());
00221         CLHEP::HepLorentzVector lv_transverse_neutrinoFromW(lv_neutrinoFromW.px(),lv_neutrinoFromW.py(),0,lv_neutrinoFromW.perp());
00222         CLHEP::HepLorentzVector lv_transverse_W(lv_transverse_leptonFromW);
00223         lv_transverse_W+=lv_transverse_neutrinoFromW;
00224         m_W_trans_mass->Fill(lv_transverse_W.m());
00225         m_W_trans_mass_log->Fill(lv_transverse_W.m());
00226         
00227       }
00228 
00229       // cut out the neutral particles (charge is not stored in HepMC)
00230       // just remove the neutrinos, gammas, neutral pions and
00231       // neutrons, K0s, ... from the final state
00232       if(chargedParticle(pid)==NOTCHARGED) continue;
00233       
00234       // now collect the tracks for the Jet algorithm
00235       // first of all we need stable particle
00236       if (!((*p)->status() == 1)) continue;
00237       // They need to have a production Vertex --> we will not select incoming Protons from the Beam 
00238       if (!(*p)->production_vertex()) continue;
00239       // And since they are stable they should not have a decay vertex. 
00240       if ((*p)->end_vertex()) continue;
00241       
00242       // fiducial range eta cut
00243       if(abs((*p)->momentum().eta()) > m_max_eta) continue;
00244       
00245       // minimal pt for tracks                                                                                                
00246       if((*p)->momentum().perp() < m_min_pt) continue;
00247       
00248       // check wether the track directly comes from a Z-Boson      
00249       if(((abs(pid)==11) || (abs(pid)==13)) && trackfromPID(23,(*p))) continue;
00250       
00251       // check wether the track directly comes from a W-Boson
00252       if(((abs(pid)==11) || (abs(pid)==13)) && trackfromPID(24,(*p))) continue;
00253       if(((abs(pid)==11) || (abs(pid)==13)) && trackfromPID(-24,(*p))) continue;
00254       
00255       // remove pi0 from charged particle list
00256       //if(pid==111) continue;
00257       
00258       m_charged_particle_pt->Fill((*p)->momentum().perp());
00259       m_charged_particle_temp_pt->Fill((*p)->momentum().perp());
00260       chargeParticlecount++;
00261       //m_charged_particle_pdgID->Fill(pid);
00262       
00263     } 
00264   
00265   // storage of input particles for jet
00266   CLHEP::HepLorentzVector lv_leadingJet(0,0,0,0);
00267   
00268   // if jet finder wasn't run before do it now:
00269   if(!m_inclusive_jets.size()) {
00270     int test = FindJet(hepmcevt);
00271     if(!test) {
00272       cout <<"no jets found - return failure" << endl;
00273       //return FAILURE;
00274     }
00275   }
00276   
00277   //Fill the histograms
00278   m_W_count->Fill(Wcount);
00279   m_lepton_count->Fill(leptoncount);
00280   m_charged_particle_multiplicity->Fill(chargeParticlecount);
00281   m_charged_particle_mean_pt->Fill(m_charged_particle_temp_pt->GetMean());
00282   m_charged_particle_rms_pt->Fill(m_charged_particle_temp_pt->GetRMS());
00283   
00284   m_jet_count->Fill(m_inclusive_jets.size());
00285 
00286   // not needed, but to be on the safe side
00287   m_charged_particle_temp_pt->Reset();
00288   
00289   if(m_inclusive_jets.size()){
00290 
00291     m_jet_pt->Fill(m_inclusive_jets[0].perp());
00292     m_jet_pt_log->Fill(m_inclusive_jets[0].perp());
00293     
00294     CLHEP::HepLorentzVector lv_firstJet(m_inclusive_jets[0].px(), m_inclusive_jets[0].py(), m_inclusive_jets[0].pz(), m_inclusive_jets[0].e());
00295     
00296     if(leptonFromW)  {
00297       if(leptoncount>1) 
00298         std::cout << "WplusJetAnalysis: WARNING: You have more than one Lepton from W-Boson in the event!" << std::endl;
00299       m_jet_DR_leadingJet_lepton->Fill(abs(lv_firstJet.deltaR(lv_leptonFromW)));
00300       m_jet_Deta_leadingJet_lepton->Fill(abs(lv_firstJet.eta()-lv_leptonFromW.eta()));
00301       m_jet_Dphi_leadingJet_lepton->Fill(abs(lv_firstJet.vect().deltaPhi(lv_leptonFromW.vect())));
00302       
00303       // Fill the angle plots for a lepton according to the leading Jet
00304       // ... if we have more than one Jet
00305       if(m_inclusive_jets.size() > 1 ) {
00306         m_jet_DR_leadingJet_lepton_multijet->Fill(abs(lv_firstJet.deltaR(lv_leptonFromW)));
00307         m_jet_Deta_leadingJet_lepton_multijet->Fill(abs(lv_firstJet.eta()-lv_leptonFromW.eta()));
00308         m_jet_Dphi_leadingJet_lepton_multijet->Fill(abs(lv_firstJet.vect().deltaPhi(lv_leptonFromW.vect())));
00309       }
00310       else{       // ... or only one Jet
00311         m_jet_DR_leadingJet_lepton_singlejet->Fill(abs(lv_firstJet.deltaR(lv_leptonFromW)));
00312         m_jet_Deta_leadingJet_lepton_singlejet->Fill(abs(lv_firstJet.eta()-lv_leptonFromW.eta()));
00313         m_jet_Dphi_leadingJet_lepton_singlejet->Fill(abs(lv_firstJet.vect().deltaPhi(lv_leptonFromW.vect())));
00314       }
00315     }
00316     
00317   }
00318 
00319   return SUCCESS;
00320 }