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"),1000, 0., 1000.);
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_temp_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_temp_pt=baseAnalysis::initHist(string("WplusJet_charged_particle_temp_mean_pt"),string("temp histogramm"),string("WplusJet charged particle temp mean p_{T}"),200,0,200.);
00087   m_charged_particle_pt=baseAnalysis::initHist(string("WplusJet_charged_particle_pt"),string("p_{T} of stable charged particles"),string("p_{T}"),200,0,200.);
00088   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.);
00089   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.);
00090   //m_charged_particle_pdgID=baseAnalysis::initHist(string("charged_particle_pdgID"),string("charged Particles PDG ID"),string("charged particle pid"),400, -400., 400.);
00091   
00092   m_jet_count=baseAnalysis::initHist(string("jet_count"), string("Nr of Jets"), string("Nr Jets"),16, -0.5, 15.5);
00093   m_jet_pt=baseAnalysis::initHist(string("leadingjet_pt_high"),string("Leading Jet Transverse Momentum"),string("Jet p_{T} [GeV]"), 200, 0.,1000.);
00094   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.);
00095 
00096 //   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.);
00097 //   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.);
00098 //   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.);
00099   
00100 //   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);
00101 //   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);
00102 //   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);
00103   
00104 //   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.);
00105 //   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.);
00106 //   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.);
00107   
00108   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.);
00109   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.);
00110   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.);
00111   
00112   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.);
00113   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.);
00114   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.);
00115   
00116   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.);
00117   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.);
00118   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.);
00119   
00120   for(unsigned int hist=0; hist< m_histVector.size(); hist++)
00121     {
00122       std::ostringstream counter_s;
00123       counter_s << hist;
00124       std::string histogram_name=m_outputRootDir;
00125       histogram_name+="_";
00126       histogram_name+=counter_s.str();
00127       histogram_name+="_";
00128       histogram_name+=m_histVector[hist]->GetName();
00129       m_histVector[hist]->SetName(histogram_name.c_str());
00130     }
00131   
00132   return true;
00133   
00134 }
00135 
00136 /**
00137 This is the main analysis function. 
00138 Search the particle, get their properties and the histograms.
00139 */
00140 int WplusJetAnalysis::Process(HepMC::GenEvent* hepmcevt) 
00141 {
00142   
00143   // Some Variables or storage classes are needed
00144   CLHEP::HepLorentzVector lv_leptonFromW(0,0,0,0);
00145   CLHEP::HepLorentzVector lv_neutrinoFromW(0,0,0,0);
00146   HepMC::GenParticle* leptonFromW=0;
00147   HepMC::GenParticle* neutrinoFromW=0;
00148   
00149   // storage of input particles for jet
00150   //int properStatus = 2;
00151   
00152   // Reset the histogramm
00153   m_charged_particle_temp_pt->Reset();
00154   
00155   // Fill the eventnumber
00156   //m_evtnr->Fill((hepmcevt->event_number())%100);
00157   
00158   // initialize the counter variables
00159   int Wcount=0;
00160   int leptoncount=0;
00161   int chargeParticlecount=0;
00162   
00163   // loop over the particles and select pdgid and pt
00164   for ( HepMC::GenEvent::particle_const_iterator p =  hepmcevt->particles_begin(); p != hepmcevt->particles_end(); ++p )
00165     {
00166       
00167       HepMC::GenParticle* WBoson=0;
00168       
00169       // use the pdg id to identify the particles
00170       int pid = (*p)->pdg_id();
00171       
00172       // W decays
00173       // look for W-Bosons (pdg = 24) 
00174       //if( abs(pid) == 24 && ((*p)->status())%1000 == properStatus && (*p)->end_vertex()){
00175       if( abs(pid) == 24 && (*p)->end_vertex()){
00176 
00177         // search for decaying particles
00178         // since the particles decay, they should have an end vertex
00179         if (!(*p)->end_vertex()) continue;
00180 
00181         HepMC::GenVertex::particle_iterator firstChild = (*p)->end_vertex()->particles_begin(HepMC::children);
00182         HepMC::GenVertex::particle_iterator endChild =  (*p)->end_vertex()->particles_end(HepMC::children);
00183 
00184         bool LastWBoson = true;
00185         for(HepMC::GenVertex::particle_iterator iter=firstChild ; iter!=endChild; iter++)
00186           {
00187             if( (*iter)->pdg_id() == pid ){
00188               LastWBoson = false;
00189             }
00190           }
00191 
00192         // plot Pt, eta and phi of the W (generator values)
00193 //      if(LastWBoson == true) {
00194 //        WBoson = (*p);
00195 //        m_W_pt->Fill((*p)->momentum().perp());
00196 //        m_W_pt_log->Fill((*p)->momentum().perp());
00197 //        m_W_pt_high->Fill((*p)->momentum().perp());
00198 //        m_W_pt_high_log->Fill((*p)->momentum().perp());
00199 //        m_W_eta->Fill((*p)->momentum().eta());
00200 //        m_W_rapidity->Fill(baseAnalysis::getRapidity(p));
00201 //        m_W_phi->Fill((*p)->momentum().phi());
00202 //        Wcount++;
00203 //      }
00204 
00205         if(LastWBoson == false) continue;
00206 
00207         WBoson = (*p);
00208         m_W_pt->Fill((*p)->momentum().perp());
00209         m_W_pt_log->Fill((*p)->momentum().perp());
00210         m_W_pt_high->Fill((*p)->momentum().perp());
00211         m_W_pt_high_log->Fill((*p)->momentum().perp());
00212         m_W_eta->Fill((*p)->momentum().eta());
00213         m_W_rapidity->Fill(baseAnalysis::getRapidity(p));
00214         m_W_phi->Fill((*p)->momentum().phi());
00215         Wcount++;
00216         
00217         bool WtoLepton=false;
00218         // Now check the decay products and search for a lepton and the neutrino
00219         for(HepMC::GenVertex::particle_iterator iter=firstChild ; iter!=endChild; iter++)
00220           {
00221             // set the properties of the neutrino Lorentzvector 
00222             if(abs((*iter)->pdg_id())==12 || abs((*iter)->pdg_id())==14 || abs((*iter)->pdg_id())==16){
00223               neutrinoFromW=(*iter);
00224               lv_neutrinoFromW.setPx((*iter)->momentum().px());
00225               lv_neutrinoFromW.setPy((*iter)->momentum().py());
00226               lv_neutrinoFromW.setPz((*iter)->momentum().pz());
00227               lv_neutrinoFromW.setE((*iter)->momentum().e());
00228             }
00229             // set the properties of the lepton Lorentzvector 
00230             if(abs((*iter)->pdg_id())==11 || abs((*iter)->pdg_id())==13 || abs((*iter)->pdg_id())==15) {
00231               leptonFromW=*iter;
00232               WtoLepton=true;
00233               m_lepton_pdgid->Fill((*iter)->pdg_id());
00234               lv_leptonFromW.setPx((*iter)->momentum().px());
00235               lv_leptonFromW.setPy((*iter)->momentum().py());
00236               lv_leptonFromW.setPz((*iter)->momentum().pz());
00237               lv_leptonFromW.setE((*iter)->momentum().e());
00238             }
00239           }
00240         
00241         if(!WtoLepton) continue;
00242         
00243         // fill the lepton histograms
00244         m_lepton_pt->Fill(lv_leptonFromW.perp());
00245         m_lepton_pt_log->Fill(lv_leptonFromW.perp());
00246         m_lepton_pt_high->Fill(lv_leptonFromW.perp());
00247         m_lepton_pt_high_log->Fill(lv_leptonFromW.perp());
00248         m_lepton_eta->Fill(lv_leptonFromW.eta());
00249         m_lepton_phi->Fill(lv_leptonFromW.phi());
00250         leptoncount++;
00251         
00252         // Use the standard Lorentzvector function to get the transverse Mass of the W
00253         // First set the z-component of the lepton and neutrino to zero and adjust the transverse energy (=transverse momentum)
00254         // then add the two lorentzvectors to get the transverse mass of the W-Boson
00255         CLHEP::HepLorentzVector lv_transverse_leptonFromW(lv_leptonFromW.px(),lv_leptonFromW.py(),0,lv_leptonFromW.perp());
00256         CLHEP::HepLorentzVector lv_transverse_neutrinoFromW(lv_neutrinoFromW.px(),lv_neutrinoFromW.py(),0,lv_neutrinoFromW.perp());
00257         CLHEP::HepLorentzVector lv_transverse_W(lv_transverse_leptonFromW);
00258         lv_transverse_W+=lv_transverse_neutrinoFromW;
00259         m_W_trans_mass->Fill(lv_transverse_W.m());
00260         m_W_trans_mass_log->Fill(lv_transverse_W.m());
00261         
00262       }
00263       
00264       // cut out the neutral particles (charge is not stored in HepMC)
00265       // just remove the neutrinos, gammas, neutral pions and
00266       // neutrons, K0s, ... from the final state
00267       if(!chargedParticle(pid)) continue;
00268       
00269       // take all stable particles
00270       if(!IsFinalStateParticle((*p))) continue;
00271       
00272       // remove pi0 from charged particle list
00273       //if(pid==111) continue;
00274       
00275       m_charged_particle_pt->Fill((*p)->momentum().perp());
00276       m_charged_particle_temp_pt->Fill((*p)->momentum().perp());
00277       chargeParticlecount++;
00278       //m_charged_particle_pdgID->Fill(pid);
00279       
00280     } 
00281   
00282   // storage of input particles for jet
00283   CLHEP::HepLorentzVector lv_leadingJet(0,0,0,0);
00284   
00285   // check if jets were found
00286   //if(!m_inclusive_jets.size()) std::cout<<"WplusJetAnalysis::no jets found"<<std::endl;
00287 
00288   //Fill the histograms
00289   m_W_count->Fill(Wcount);
00290   m_lepton_count->Fill(leptoncount);
00291   m_charged_particle_multiplicity->Fill(chargeParticlecount);
00292   m_charged_particle_mean_pt->Fill(m_charged_particle_temp_pt->GetMean());
00293   m_charged_particle_rms_pt->Fill(m_charged_particle_temp_pt->GetRMS());
00294   
00295   m_jet_count->Fill(m_inclusive_jets.size());
00296 
00297   // not needed, but to be on the safe side
00298   m_charged_particle_temp_pt->Reset();
00299   
00300   if(m_inclusive_jets.size()){
00301 
00302     m_jet_pt->Fill(m_inclusive_jets[0].perp());
00303     m_jet_pt_log->Fill(m_inclusive_jets[0].perp());
00304     
00305     CLHEP::HepLorentzVector lv_firstJet(m_inclusive_jets[0].px(), m_inclusive_jets[0].py(), m_inclusive_jets[0].pz(), m_inclusive_jets[0].e());
00306     
00307     if(leptonFromW)  {
00308       if(leptoncount>1) 
00309         std::cout << "WplusJetAnalysis: WARNING: You have more than one Lepton from W-Boson in the event!" << std::endl;
00310       double perp1 = lv_leptonFromW.perp();
00311       double perp2 = lv_firstJet.perp();
00312       double eta1 = lv_leptonFromW.eta();
00313       double eta2 = lv_firstJet.eta();
00314       double phi1 = lv_leptonFromW.phi();
00315       double phi2 = lv_firstJet.phi();
00316       // skip leptons which are not isolated, because they are in the jet and in the W decay
00317       if(fabs(lv_firstJet.deltaR(lv_leptonFromW)) < 0.0005 && abs(lv_firstJet.deltaR(lv_leptonFromW)) > -0.0005
00318          && perp1 == perp2 && eta1==eta2 && phi1==phi2) return false;
00319 
00320       m_jet_DR_leadingJet_lepton->Fill(fabs(lv_firstJet.deltaR(lv_leptonFromW)));
00321       m_jet_Deta_leadingJet_lepton->Fill(fabs(lv_firstJet.eta()-lv_leptonFromW.eta()));
00322       m_jet_Dphi_leadingJet_lepton->Fill(fabs(lv_firstJet.vect().deltaPhi(lv_leptonFromW.vect())));
00323 
00324       // Fill the angle plots for a lepton according to the leading Jet
00325       // ... if we have more than one Jet
00326       if(m_inclusive_jets.size() > 1 ) {
00327         m_jet_DR_leadingJet_lepton_multijet->Fill(fabs(lv_firstJet.deltaR(lv_leptonFromW)));
00328         m_jet_Deta_leadingJet_lepton_multijet->Fill(fabs(lv_firstJet.eta()-lv_leptonFromW.eta()));
00329         m_jet_Dphi_leadingJet_lepton_multijet->Fill(fabs(lv_firstJet.vect().deltaPhi(lv_leptonFromW.vect())));
00330       }
00331       else{       // ... or only one Jet
00332         m_jet_DR_leadingJet_lepton_singlejet->Fill(fabs(lv_firstJet.deltaR(lv_leptonFromW)));
00333         m_jet_Deta_leadingJet_lepton_singlejet->Fill(fabs(lv_firstJet.eta()-lv_leptonFromW.eta()));
00334         m_jet_Dphi_leadingJet_lepton_singlejet->Fill(fabs(lv_firstJet.vect().deltaPhi(lv_leptonFromW.vect())));
00335       }
00336     }
00337   }
00338 
00339   return true;
00340 }