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