CMSSW Reference Manual

00001 // -*- C++ -*-
00002 //
00003 // Package:    RecoEgamma/Examples
00004 // Class:      GsfElectronMCAnalyzer
00005 // 
00013 //
00014 // Original Author:  Ursula Berthon
00015 //         Created:  Mon Mar 27 13:22:06 CEST 2006
00016 // $Id: GsfElectronMCAnalyzer.cc,v 1.4 2008/09/13 08:30:26 charlot Exp $
00017 //
00018 //
00019 
00020 // user include files
00021 #include "RecoEgamma/Examples/plugins/GsfElectronMCAnalyzer.h"
00022 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00023 #include "FWCore/Framework/interface/EDAnalyzer.h"
00024 #include "FWCore/Framework/interface/Event.h"
00025 #include "FWCore/Framework/interface/MakerMacros.h"
00026 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00027 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
00028 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
00029 #include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
00030 #include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
00031 #include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
00032 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
00033 #include "SimDataFormats/HepMCProduct/interface/HepMCProduct.h"
00034 
00035 #include <iostream>
00036 #include "TMath.h"
00037 #include "TFile.h"
00038 #include "TH1F.h"
00039 #include "TH1I.h"
00040 #include "TH2F.h"
00041 #include "TProfile.h"
00042 #include "TTree.h"
00043 #include <iostream>
00044 
00045 using namespace reco;
00046  
00047 GsfElectronMCAnalyzer::GsfElectronMCAnalyzer(const edm::ParameterSet& conf)
00048 {
00049 
00050   outputFile_ = conf.getParameter<std::string>("outputFile");
00051   histfile_ = new TFile(outputFile_.c_str(),"RECREATE");
00052   electronCollection_=conf.getParameter<edm::InputTag>("electronCollection");
00053   mcTruthCollection_ = conf.getParameter<edm::InputTag>("mcTruthCollection");
00054   maxPt_ = conf.getParameter<double>("MaxPt");
00055   maxAbsEta_ = conf.getParameter<double>("MaxAbsEta");
00056   deltaR_ = conf.getParameter<double>("DeltaR");
00057   etamin=conf.getParameter<double>("Etamin");
00058   etamax=conf.getParameter<double>("Etamax");
00059   phimin=conf.getParameter<double>("Phimin");
00060   phimax=conf.getParameter<double>("Phimax");
00061   ptmax=conf.getParameter<double>("Ptmax");
00062   pmax=conf.getParameter<double>("Pmax");
00063   eopmax=conf.getParameter<double>("Eopmax");
00064   eopmaxsht=conf.getParameter<double>("Eopmaxsht");
00065   detamin=conf.getParameter<double>("Detamin");
00066   detamax=conf.getParameter<double>("Detamax");
00067   dphimin=conf.getParameter<double>("Dphimin");
00068   dphimax=conf.getParameter<double>("Dphimax");
00069   detamatchmin=conf.getParameter<double>("Detamatchmin");
00070   detamatchmax=conf.getParameter<double>("Detamatchmax");
00071   dphimatchmin=conf.getParameter<double>("Dphimatchmin");
00072   dphimatchmax=conf.getParameter<double>("Dphimatchmax");
00073   fhitsmax=conf.getParameter<double>("Fhitsmax");
00074   lhitsmax=conf.getParameter<double>("Lhitsmax");
00075   nbineta=conf.getParameter<int>("Nbineta");
00076   nbineta2D=conf.getParameter<int>("Nbineta2D");
00077   nbinp=conf.getParameter<int>("Nbinp");
00078   nbinpt=conf.getParameter<int>("Nbinpt");
00079   nbinp2D=conf.getParameter<int>("Nbinp2D");
00080   nbinpt2D=conf.getParameter<int>("Nbinpt2D");
00081   nbinpteff=conf.getParameter<int>("Nbinpteff");
00082   nbinphi=conf.getParameter<int>("Nbinphi");
00083   nbinphi2D=conf.getParameter<int>("Nbinphi2D");
00084   nbineop=conf.getParameter<int>("Nbineop");
00085   nbineop2D=conf.getParameter<int>("Nbineop2D");
00086   nbinfhits=conf.getParameter<int>("Nbinfhits");
00087   nbinlhits=conf.getParameter<int>("Nbinlhits");
00088   nbinxyz=conf.getParameter<int>("Nbinxyz");
00089   nbindeta=conf.getParameter<int>("Nbindeta");
00090   nbindphi=conf.getParameter<int>("Nbindphi");
00091   nbindetamatch=conf.getParameter<int>("Nbindetamatch");
00092   nbindphimatch=conf.getParameter<int>("Nbindphimatch");
00093   nbindetamatch2D=conf.getParameter<int>("Nbindetamatch2D");
00094   nbindphimatch2D=conf.getParameter<int>("Nbindphimatch2D");
00095 }  
00096   
00097 GsfElectronMCAnalyzer::~GsfElectronMCAnalyzer()
00098 {
00099  
00100   // do anything here that needs to be done at desctruction time
00101   // (e.g. close files, deallocate resources etc.)
00102   histfile_->Write();
00103   histfile_->Close();
00104 }
00105 
00106 void GsfElectronMCAnalyzer::beginJob(edm::EventSetup const&iSetup){
00107 
00108   histfile_->cd();
00109   
00110   
00111   // mc truth  
00112 
00113   h_mcNum              = new TH1F( "h_mcNum",              "# mc particles",    nbinfhits,0.,fhitsmax );
00114   h_eleNum             = new TH1F( "h_mcNum_ele",             "# mc electrons",             nbinfhits,0.,fhitsmax);
00115   h_gamNum             = new TH1F( "h_mcNum_gam",             "# mc gammas",             nbinfhits,0.,fhitsmax);
00116     
00117   // rec event
00118   
00119   histNum_= new TH1F("h_recEleNum","# rec electrons",20, 0.,20.);
00120   
00121   // mc  
00122   h_simEta             = new TH1F( "h_mc_eta",             "mc #eta",           nbineta,etamin,etamax); 
00123   h_simAbsEta             = new TH1F( "h_mc_abseta",             "mc |#eta|",           nbineta/2,0.,etamax); 
00124   h_simP               = new TH1F( "h_mc_P",               "mc p",              nbinp,0.,pmax); 
00125   h_simPt               = new TH1F( "h_mc_Pt",               "mc pt",            nbinpteff,5.,ptmax); 
00126   h_simPhi               = new TH1F( "h_mc_phi",               "mc phi",        nbinphi,phimin,phimax); 
00127   h_simZ      = new TH1F( "h_mc_z",      "mc z ",    50, -25, 25 );
00128 
00129   // ctf tracks
00130   h_ctf_foundHitsVsEta      = new TH2F( "h_ctf_foundHitsVsEta",      "ctf track # found hits vs eta",  nbineta2D,etamin,etamax,nbinfhits,0.,fhitsmax);
00131   h_ctf_lostHitsVsEta       = new TH2F( "h_ctf_lostHitsVsEta",       "ctf track # lost hits vs eta",   nbineta2D,etamin,etamax,nbinlhits,0.,lhitsmax);
00132   
00133   // all electrons  
00134   h_ele_EoverP_all       = new TH1F( "h_ele_EoverP_all",       "all ele E/p at vertex",  nbineop,0.,eopmax);
00135   h_ele_TIP_all       = new TH1F( "h_ele_TIP_all",       "all ele tip at vertex",  100,0.,0.2);
00136   h_ele_vertexPt_all       = new TH1F( "h_ele_vertexPt_all",       "all ele p_{T} at vertex",  nbinpteff,5.,ptmax);
00137   h_ele_vertexEta_all      = new TH1F( "h_ele_vertexEta_all",      "all ele #eta at vertex",    nbineta,etamin,etamax);
00138   h_ele_mee_all      = new TH1F( "h_ele_mee_all", "all ele pairs invariant mass", 100, 0., 150. );
00139 
00140   // matched electrons
00141   h_ele_charge         = new TH1F( "h_ele_charge",         "ele charge",             5,-2.,2.);   
00142   h_ele_chargeVsEta    = new TH2F( "h_ele_chargeVsEta",         "ele charge vs eta", nbineta2D,etamin,etamax,5,-2.,2.);   
00143   h_ele_chargeVsPhi    = new TH2F( "h_ele_chargeVsPhi",         "ele charge vs phi", nbinphi2D,phimin,phimax,5,-2.,2.);   
00144   h_ele_chargeVsPt    = new TH2F( "h_ele_chargeVsPt",         "ele charge vs pt", nbinpt,0.,100.,5,-2.,2.);   
00145   h_ele_vertexP        = new TH1F( "h_ele_vertexP",        "ele p at vertex",       nbinp,0.,pmax);
00146   h_ele_vertexPt       = new TH1F( "h_ele_vertexPt",       "ele p_{T} at vertex",  nbinpt,0.,ptmax);
00147   h_ele_vertexPtVsEta   = new TH2F( "h_ele_vertexPtVsEta",       "ele p_{T} at vertex vs eta",nbinpt2D,etamin,etamax,nbinpt2D,0.,ptmax);
00148   h_ele_vertexPtVsPhi   = new TH2F( "h_ele_vertexPtVsPhi",       "ele p_{T} at vertex vs phi",nbinphi2D,phimin,phimax,nbinpt2D,0.,ptmax);
00149   h_ele_simPt_matched       = new TH1F( "h_ele_simPt_matched",       "sim p_{T}, matched ",  nbinpteff,5.,ptmax);
00150   h_ele_vertexEta      = new TH1F( "h_ele_vertexEta",      "ele #eta at vertex",    nbineta,etamin,etamax);
00151   h_ele_vertexEtaVsPhi  = new TH2F( "h_ele_vertexEtaVsPhi",      "ele #eta at vertex vs phi",nbineta2D,etamin,etamax,nbinphi2D,phimin,phimax );
00152   h_ele_simAbsEta_matched      = new TH1F( "h_ele_simAbsEta_matched",      "sim |#eta|, matched ",    nbineta/2,0.,2.5);
00153   h_ele_simEta_matched      = new TH1F( "h_ele_simEta_matched",      "sim #eta,matched ",    nbineta,etamin,etamax);
00154   h_ele_simPhi_matched               = new TH1F( "h_ele_simPhi_matched",               "sim phi, matched ",        nbinphi,phimin,phimax); 
00155   h_ele_vertexPhi      = new TH1F( "h_ele_vertexPhi",      "ele #phi at vertex",    nbinphi,phimin,phimax);
00156   h_ele_vertexX      = new TH1F( "h_ele_vertexX",      "ele x at vertex",    nbinxyz,-0.1,0.1 );
00157   h_ele_vertexY      = new TH1F( "h_ele_vertexY",      "ele y at vertex",    nbinxyz,-0.1,0.1 );
00158   h_ele_vertexZ      = new TH1F( "h_ele_vertexZ",      "ele z at vertex",    nbinxyz,-25, 25 );
00159   h_ele_simZ_matched      = new TH1F( "h_ele_simZ_matched",      "sim z, matched ",    nbinxyz,-25,25);
00160   h_ele_vertexTIP      = new TH1F( "h_ele_vertexTIP",      "ele TIP",    90,0.,0.15);
00161   h_ele_vertexTIPVsEta      = new TH2F( "h_ele_vertexTIPVsEta",      "ele TIP vs eta", nbineta2D,etamin,etamax,45,0.,0.15);
00162   h_ele_vertexTIPVsPhi      = new TH2F( "h_ele_vertexTIPVsPhi",      "ele TIP vs phi", nbinphi2D,phimin,phimax,45,0.,0.15);
00163   h_ele_vertexTIPVsPt      = new TH2F( "h_ele_vertexTIPVsPt",      "ele TIP vs Pt", nbinpt2D,0.,ptmax,45,0.,0.15);
00164   h_ele_PoPtrue        = new TH1F( "h_ele_PoPtrue",        "ele P/Ptrue @ vertex", 75,0.,1.5);
00165   h_ele_PoPtrueVsEta   = new TH2F( "h_ele_PoPtrueVsEta",        "ele P/Ptrue @ vertex vs eta", nbineta2D,etamin,etamax,50,0.,1.5);
00166   h_ele_PoPtrueVsPhi   = new TH2F( "h_ele_PoPtrueVsPhi",        "ele P/Ptrue @ vertex vs phi", nbinphi2D,phimin,phimax,50,0.,1.5);
00167   h_ele_PoPtrueVsPt   = new TH2F( "h_ele_PoPtrueVsPt",        "ele P/Ptrue @ vertex vs eta", nbinpt2D,0.,ptmax,50,0.,1.5);
00168   h_ele_PoPtrue_barrel         = new TH1F( "h_ele_PoPtrue_barrel",        "ele P/Ptrue @ vertex, barrel",75,0.,1.5);
00169   h_ele_PoPtrue_endcaps        = new TH1F( "h_ele_PoPtrue_endcaps",        "ele P/Ptrue @ vertex, endcaps",75,0.,1.5);
00170   h_ele_EtaMnEtaTrue   = new TH1F( "h_ele_EtaMnEtaTrue",   "ele #eta_{rec} - #eta_{sim} @ vertex",nbindeta,detamin,detamax);
00171   h_ele_EtaMnEtaTrueVsEta   = new TH2F( "h_ele_EtaMnEtaTrueVsEta",   "ele #eta_{rec} - #eta_{sim} @ vertex vs eta",nbineta2D,etamin,etamax,nbindeta/2,detamin,detamax);
00172   h_ele_EtaMnEtaTrueVsPhi   = new TH2F( "h_ele_EtaMnEtaTrueVsPhi",   "ele #eta_{rec} - #eta_{sim} @ vertex vs phi",nbinphi2D,phimin,phimax,nbindeta/2,detamin,detamax);
00173   h_ele_EtaMnEtaTrueVsPt   = new TH2F( "h_ele_EtaMnEtaTrueVsPt",   "ele #eta_{rec} - #eta_{sim} @ vertex vs pt",nbinpt,0.,ptmax,nbindeta/2,detamin,detamax);
00174   h_ele_PhiMnPhiTrue   = new TH1F( "h_ele_PhiMnPhiTrue",   "ele #phi_{rec} - #phi_{sim} @ vertex",nbindphi,dphimin,dphimax);
00175   h_ele_PhiMnPhiTrue2   = new TH1F( "h_ele_PhiMnPhiTrue2",   "ele #phi_{rec} - #phi_{sim} @ vertex",nbindphimatch2D,dphimatchmin,dphimatchmax);
00176   h_ele_PhiMnPhiTrueVsEta   = new TH2F( "h_ele_PhiMnPhiTrueVsEta",   "ele #phi_{rec} - #phi_{sim} @ vertex vs eta",nbineta2D,etamin,etamax,nbindphi/2,dphimin,dphimax);
00177   h_ele_PhiMnPhiTrueVsPhi   = new TH2F( "h_ele_PhiMnPhiTrueVsPhi",   "ele #phi_{rec} - #phi_{sim} @ vertex vs phi",nbinphi2D,phimin,phimax,nbindphi/2,dphimin,dphimax);
00178   h_ele_PhiMnPhiTrueVsPt   = new TH2F( "h_ele_PhiMnPhiTrueVsPt",   "ele #phi_{rec} - #phi_{sim} @ vertex vs pt",nbinpt2D,0.,ptmax,nbindphi/2,dphimin,dphimax);
00179 
00180   // matched electron, superclusters
00181   histSclEn_ = new TH1F("h_scl_energy","ele supercluster energy",nbinp,0.,pmax);
00182   histSclEoEtrue_barrel = new TH1F("h_scl_EoEtrue_barrel","ele supercluster energy over true energy, barrel",50,0.2,1.2);
00183   histSclEoEtrue_endcaps = new TH1F("h_scl_EoEtrue_endcaps","ele supercluster energy over true energy, endcaps",50,0.2,1.2);
00184   histSclEt_ = new TH1F("h_scl_et","ele supercluster transverse energy",nbinpt,0.,ptmax);
00185   histSclEtVsEta_ = new TH2F("h_scl_etVsEta","ele supercluster transverse energy vs eta",nbineta2D,etamin,etamax,nbinpt,0.,ptmax);
00186   histSclEtVsPhi_ = new TH2F("h_scl_etVsPhi","ele supercluster transverse energy vs phi",nbinphi2D,phimin,phimax,nbinpt,0.,ptmax);
00187   histSclEtaVsPhi_ = new TH2F("h_scl_etaVsPhi","ele supercluster eta vs phi",nbinphi2D,phimin,phimax,nbineta2D,etamin,etamax);
00188   histSclEta_ = new TH1F("h_scl_eta","ele supercluster eta",nbineta,etamin,etamax);
00189   histSclPhi_ = new TH1F("h_scl_phi","ele supercluster phi",nbinphi,phimin,phimax);
00190 
00191   // matched electron, gsf tracks
00192   h_ele_foundHits      = new TH1F( "h_ele_foundHits",      "ele track # found hits",      nbinfhits,0.,fhitsmax);
00193   h_ele_foundHitsVsEta      = new TH2F( "h_ele_foundHitsVsEta",      "ele track # found hits vs eta",  nbineta2D,etamin,etamax,nbinfhits,0.,fhitsmax);
00194   h_ele_foundHitsVsPhi      = new TH2F( "h_ele_foundHitsVsPhi",      "ele track # found hits vs phi",  nbinphi2D,phimin,phimax,nbinfhits,0.,fhitsmax);
00195   h_ele_foundHitsVsPt      = new TH2F( "h_ele_foundHitsVsPt",      "ele track # found hits vs pt",  nbinpt2D,0.,ptmax,nbinfhits,0.,fhitsmax);
00196   h_ctf_foundHits      = new TH1F( "h_ctf_foundHits",      "ctf track # found hits",      nbinfhits,0.,fhitsmax);
00197   h_ele_lostHits       = new TH1F( "h_ele_lostHits",       "ele track # lost hits",       5,0.,5.);
00198   h_ele_lostHitsVsEta       = new TH2F( "h_ele_lostHitsVsEta",       "ele track # lost hits vs eta",   nbineta2D,etamin,etamax,nbinlhits,0.,lhitsmax);
00199   h_ele_lostHitsVsPhi       = new TH2F( "h_ele_lostHitsVsPhi",       "ele track # lost hits vs eta",   nbinphi2D,phimin,phimax,nbinlhits,0.,lhitsmax);
00200   h_ele_lostHitsVsPt       = new TH2F( "h_ele_lostHitsVsPt",       "ele track # lost hits vs eta",   nbinpt2D,0.,ptmax,nbinlhits,0.,lhitsmax);
00201   h_ele_chi2           = new TH1F( "h_ele_chi2",           "ele track #chi^{2}",         100,0.,15.);   
00202   h_ele_chi2VsEta           = new TH2F( "h_ele_chi2VsEta",           "ele track #chi^{2} vs eta",  nbineta2D,etamin,etamax,50,0.,15.);   
00203   h_ele_chi2VsPhi           = new TH2F( "h_ele_chi2VsPhi",           "ele track #chi^{2} vs phi",  nbinphi2D,phimin,phimax,50,0.,15.);   
00204   h_ele_chi2VsPt           = new TH2F( "h_ele_chi2VsPt",           "ele track #chi^{2} vs pt",  nbinpt2D,0.,ptmax,50,0.,15.);   
00205   h_ele_PinMnPout      = new TH1F( "h_ele_PinMnPout",      "ele track inner p - outer p, mean"   ,nbinp,0.,200.);
00206   h_ele_PinMnPout_mode      = new TH1F( "h_ele_PinMnPout_mode",      "ele track inner p - outer p, mode"   ,nbinp,0.,100.);
00207   h_ele_PinMnPoutVsEta_mode = new TH2F( "h_ele_PinMnPoutVsEta_mode",      "ele track inner p - outer p vs eta, mode" ,nbineta2D, etamin,etamax,nbinp2D,0.,100.);
00208   h_ele_PinMnPoutVsPhi_mode = new TH2F( "h_ele_PinMnPoutVsPhi_mode",      "ele track inner p - outer p vs phi, mode" ,nbinphi2D, phimin,phimax,nbinp2D,0.,100.);
00209   h_ele_PinMnPoutVsPt_mode = new TH2F( "h_ele_PinMnPoutVsPt_mode",      "ele track inner p - outer p vs pt, mode" ,nbinpt2D, 0.,ptmax,nbinp2D,0.,100.);
00210   h_ele_PinMnPoutVsE_mode = new TH2F( "h_ele_PinMnPoutVsE_mode",      "ele track inner p - outer p vs E, mode" ,nbinp2D, 0.,200.,nbinp2D,0.,100.);
00211   h_ele_PinMnPoutVsChi2_mode = new TH2F( "h_ele_PinMnPoutVsChi2_mode",      "ele track inner p - outer p vs track chi2, mode" ,50, 0.,20.,nbinp2D,0.,100.);
00212   h_ele_outerP         = new TH1F( "h_ele_outerP",         "ele track outer p, mean",          nbinp,0.,pmax);
00213   h_ele_outerP_mode         = new TH1F( "h_ele_outerP_mode",         "ele track outer p, mode",          nbinp,0.,pmax);
00214   h_ele_outerPVsEta_mode         = new TH2F( "h_ele_outerPVsEta_mode",         "ele track outer p vs eta mode", nbineta2D,etamin,etamax,50,0.,pmax);
00215   h_ele_outerPt        = new TH1F( "h_ele_outerPt",        "ele track outer p_{T}, mean",      nbinpt,0.,ptmax);
00216   h_ele_outerPt_mode        = new TH1F( "h_ele_outerPt_mode",        "ele track outer p_{T}, mode",      nbinpt,0.,ptmax);
00217   h_ele_outerPtVsEta_mode        = new TH2F( "h_ele_outerPtVsEta_mode", "ele track outer p_{T} vs eta, mode", nbineta2D,etamin,etamax,nbinpt2D,0.,ptmax);
00218   h_ele_outerPtVsPhi_mode        = new TH2F( "h_ele_outerPtVsPhi_mode", "ele track outer p_{T} vs phi, mode", nbinphi2D,phimin,phimax,nbinpt2D,0.,ptmax);
00219   h_ele_outerPtVsPt_mode        = new TH2F( "h_ele_outerPtVsPt_mode", "ele track outer p_{T} vs pt, mode", nbinpt2D,0.,100.,nbinpt2D,0.,ptmax);
00220   
00221   // matched electrons, matching 
00222   h_ele_EoP            = new TH1F( "h_ele_EoP",            "ele E/P_{vertex}",        nbineop,0.,eopmax);
00223   h_ele_EoPVsEta            = new TH2F( "h_ele_EoPVsEta",            "ele E/P_{vertex} vs eta",  nbineta2D,etamin,etamax,nbineop2D,0.,eopmaxsht);
00224   h_ele_EoPVsPhi            = new TH2F( "h_ele_EoPVsPhi",            "ele E/P_{vertex} vs phi",  nbinphi2D,phimin,phimax,nbineop2D,0.,eopmaxsht);
00225   h_ele_EoPVsE            = new TH2F( "h_ele_EoPVsE",            "ele E/P_{vertex} vs E",  50,0.,pmax ,50,0.,5.);
00226   h_ele_EoPout         = new TH1F( "h_ele_EoPout",         "ele E/P_{out}",           nbineop,0.,eopmax);
00227   h_ele_EoPoutVsEta         = new TH2F( "h_ele_EoPoutVsEta",         "ele E/P_{out} vs eta",    nbineta2D,etamin,etamax,nbineop2D,0.,eopmaxsht);
00228   h_ele_EoPoutVsPhi         = new TH2F( "h_ele_EoPoutVsPhi",         "ele E/P_{out} vs phi",    nbinphi2D,phimin,phimax,nbineop2D,0.,eopmaxsht);
00229   h_ele_EoPoutVsE         = new TH2F( "h_ele_EoPoutVsE",         "ele E/P_{out} vs E",    nbinp2D,0.,pmax,nbineop2D,0.,eopmaxsht);
00230   h_ele_dEtaSc_propVtx = new TH1F( "h_ele_dEtaSc_propVtx", "ele #eta_{sc} - #eta_{tr} - prop from vertex",      nbindetamatch,detamatchmin,detamatchmax);
00231   h_ele_dEtaScVsEta_propVtx = new TH2F( "h_ele_dEtaScVsEta_propVtx", "ele #eta_{sc} - #eta_{tr} vs eta, prop from vertex", nbineta2D,etamin,etamax,nbindetamatch2D,detamatchmin,detamatchmax);
00232   h_ele_dEtaScVsPhi_propVtx = new TH2F( "h_ele_dEtaScVsPhi_propVtx", "ele #eta_{sc} - #eta_{tr} vs phi, prop from vertex", nbinphi2D,phimin,phimax,nbindetamatch2D,detamatchmin,detamatchmax);
00233   h_ele_dEtaScVsPt_propVtx = new TH2F( "h_ele_dEtaScVsPt_propVtx", "ele #eta_{sc} - #eta_{tr} vs pt, prop from vertex", nbinpt2D,0.,ptmax,nbindetamatch2D,detamatchmin,detamatchmax);
00234   h_ele_dPhiSc_propVtx = new TH1F( "h_ele_dPhiSc_propVtx", "ele #phi_{sc} - #phi_{tr} - prop from vertex",      nbindphimatch,dphimatchmin,dphimatchmax);
00235   h_ele_dPhiScVsEta_propVtx = new TH2F( "h_ele_dPhiScVsEta_propVtx", "ele #phi_{sc} - #phi_{tr} vs eta, prop from vertex", nbineta2D,etamin,etamax,nbindphimatch2D,dphimatchmin,dphimatchmax);
00236   h_ele_dPhiScVsPhi_propVtx = new TH2F( "h_ele_dPhiScVsPhi_propVtx", "ele #phi_{sc} - #phi_{tr} vs phi, prop from vertex", nbinphi2D,phimin,phimax,nbindphimatch2D,dphimatchmin,dphimatchmax);
00237   h_ele_dPhiScVsPt_propVtx = new TH2F( "h_ele_dPhiScVsPt_propVtx", "ele #phi_{sc} - #phi_{tr} vs pt, prop from vertex", nbinpt2D,0.,ptmax,nbindphimatch2D,dphimatchmin,dphimatchmax);
00238   h_ele_dEtaCl_propOut = new TH1F( "h_ele_dEtaCl_propOut", "ele #eta_{cl} - #eta_{tr} - prop from outermost",   nbindetamatch,detamatchmin,detamatchmax);
00239   h_ele_dEtaClVsEta_propOut = new TH2F( "h_ele_dEtaClVsEta_propOut", "ele #eta_{cl} - #eta_{tr} vs eta, prop from out", nbineta2D,etamin,etamax,nbindetamatch2D,detamatchmin,detamatchmax);
00240   h_ele_dEtaClVsPhi_propOut = new TH2F( "h_ele_dEtaClVsPhi_propOut", "ele #eta_{cl} - #eta_{tr} vs phi, prop from out", nbinphi2D,phimin,phimax,nbindetamatch2D,detamatchmin,detamatchmax);
00241   h_ele_dEtaClVsPt_propOut = new TH2F( "h_ele_dEtaScVsPt_propOut", "ele #eta_{cl} - #eta_{tr} vs pt, prop from out", nbinpt2D,0.,ptmax,nbindetamatch2D,detamatchmin,detamatchmax);
00242   h_ele_dPhiCl_propOut = new TH1F( "h_ele_dPhiCl_propOut", "ele #phi_{cl} - #phi_{tr} - prop from outermost",   nbindphimatch,dphimatchmin,dphimatchmax);
00243   h_ele_dPhiClVsEta_propOut = new TH2F( "h_ele_dPhiClVsEta_propOut", "ele #phi_{cl} - #phi_{tr} vs eta, prop from out", nbineta2D,etamin,etamax,nbindphimatch2D,dphimatchmin,dphimatchmax);
00244   h_ele_dPhiClVsPhi_propOut = new TH2F( "h_ele_dPhiClVsPhi_propOut", "ele #phi_{cl} - #phi_{tr} vs phi, prop from out", nbinphi2D,phimin,phimax,nbindphimatch2D,dphimatchmin,dphimatchmax);
00245   h_ele_dPhiClVsPt_propOut = new TH2F( "h_ele_dPhiSClsPt_propOut", "ele #phi_{cl} - #phi_{tr} vs pt, prop from out", nbinpt2D,0.,ptmax,nbindphimatch2D,dphimatchmin,dphimatchmax);
00246   
00247   h_ele_HoE = new TH1F("h_ele_HoE", "ele H/E", 55,-0.05,0.5) ;
00248   h_ele_HoEVsEta = new TH2F("h_ele_HoEVsEta", "ele H/E vs eta", nbineta,etamin,etamax,55,-0.05,0.5) ;
00249   h_ele_HoEVsPhi = new TH2F("h_ele_HoEVsPhi", "ele H/E vs phi", nbinphi2D,phimin,phimax,55,-0.05,0.5) ;
00250   h_ele_HoEVsE = new TH2F("h_ele_HoEVsE", "ele H/E vs E", nbinp, 0.,300.,55,-0.05,0.5) ;
00251  
00252   // classes  
00253   h_ele_classes = new TH1F( "h_ele_classes", "electron classes",      150,0.0,150.);
00254   h_ele_eta = new TH1F( "h_ele_eta", "electron eta",  nbineta/2,0.0,etamax);
00255   h_ele_eta_golden = new TH1F( "h_ele_eta_golden", "electron eta golden",  nbineta/2,0.0,etamax);
00256   h_ele_eta_bbrem = new TH1F( "h_ele_eta_bbrem", "electron eta bbrem",  nbineta/2,0.0,etamax);
00257   h_ele_eta_narrow = new TH1F( "h_ele_eta_narrow", "electron eta narrow",  nbineta/2,0.0,etamax);
00258   h_ele_eta_shower = new TH1F( "h_ele_eta_show", "electron eta showering",  nbineta/2,0.0,etamax);
00259   h_ele_PinVsPoutGolden_mode = new TH2F( "h_ele_PinVsPoutGolden_mode",      "ele track inner p vs outer p vs eta, golden, mode" ,nbinp2D,0.,pmax,50,0.,pmax);
00260   h_ele_PinVsPoutShowering0_mode = new TH2F( "h_ele_PinVsPoutShowering0_mode",      "ele track inner p vs outer p vs eta, showering0, mode" ,nbinp2D,0.,pmax,50,0.,pmax);
00261   h_ele_PinVsPoutShowering1234_mode = new TH2F( "h_ele_PinVsPoutShowering1234_mode",      "ele track inner p vs outer p vs eta, showering1234, mode" ,nbinp2D,0.,pmax,50,0.,pmax);
00262   h_ele_PinVsPoutGolden_mean = new TH2F( "h_ele_PinVsPoutGolden_mean",      "ele track inner p vs outer p vs eta, golden, mean" ,nbinp2D,0.,pmax,50,0.,pmax);
00263   h_ele_PinVsPoutShowering0_mean = new TH2F( "h_ele_PinVsPoutShowering0_mean",      "ele track inner p vs outer p vs eta, showering0, mean" ,nbinp2D,0.,pmax,50,0.,pmax);
00264   h_ele_PinVsPoutShowering1234_mean = new TH2F( "h_ele_PinVsPoutShowering1234_mean",      "ele track inner p vs outer p vs eta, showering1234, mean" ,nbinp2D,0.,pmax,50,0.,pmax);
00265   h_ele_PtinVsPtoutGolden_mode = new TH2F( "h_ele_PtinVsPtoutGolden_mode",      "ele track inner pt vs outer pt vs eta, golden, mode" ,nbinpt2D,0.,ptmax,50,0.,ptmax);
00266   h_ele_PtinVsPtoutShowering0_mode = new TH2F( "h_ele_PtinVsPtoutShowering0_mode",      "ele track inner pt vs outer pt vs eta, showering0, mode" ,nbinpt2D,0.,ptmax,50,0.,ptmax);
00267   h_ele_PtinVsPtoutShowering1234_mode = new TH2F( "h_ele_PtinVsPtoutShowering1234_mode",      "ele track inner pt vs outer pt vs eta, showering1234, mode" ,nbinpt2D,0.,ptmax,50,0.,ptmax);
00268   h_ele_PtinVsPtoutGolden_mean = new TH2F( "h_ele_PtinVsPtoutGolden_mean",      "ele track inner pt vs outer pt vs eta, golden, mean" ,nbinpt2D,0.,ptmax,50,0.,ptmax);
00269   h_ele_PtinVsPtoutShowering0_mean = new TH2F( "h_ele_PtinVsPtoutShowering0_mean",      "ele track inner pt vs outer pt vs eta, showering0, mean" ,nbinpt2D,0.,ptmax,50,0.,ptmax);
00270   h_ele_PtinVsPtoutShowering1234_mean = new TH2F( "h_ele_PtinVsPtoutShowering1234_mean",      "ele track inner pt vs outer pt vs eta, showering1234, mean" ,nbinpt2D,0.,ptmax,50,0.,ptmax);
00271   histSclEoEtrueGolden_barrel = new TH1F("h_scl_EoEtrue golden, barrel","ele supercluster energy over true energy, golden, barrel",100,0.2,1.2);
00272   histSclEoEtrueGolden_endcaps = new TH1F("h_scl_EoEtrue golden, endcaps","ele supercluster energy over true energy, golden, endcaps",100,0.2,1.2);
00273   histSclEoEtrueShowering0_barrel = new TH1F("h_scl_EoEtrue showering0, barrel","ele supercluster energy over true energy, showering0, barrel",100,0.2,1.2);
00274   histSclEoEtrueShowering0_endcaps = new TH1F("h_scl_EoEtrue showering0, endcaps","ele supercluster energy over true energy, showering0, endcaps",100,0.2,1.2);
00275   histSclEoEtrueShowering1234_barrel = new TH1F("h_scl_EoEtrue showering1234, barrel","ele supercluster energy over true energy, showering1234, barrel",100,0.2,1.2);
00276   histSclEoEtrueShowering1234_endcaps = new TH1F("h_scl_EoEtrue showering1234, endcaps","ele supercluster energy over true energy, showering1234, endcaps",100,0.2,1.2);
00277 
00278   // fbrem
00279   h_ele_fbremVsEta_mode = new TProfile( "h_ele_fbremvsEtamode","mean pout/pin vs eta, mode",nbineta2D,etamin,etamax,0.,1.);
00280   h_ele_fbremVsEta_mean = new TProfile( "h_ele_fbremvsEtamean","mean pout/pin vs eta, mean",nbineta2D,etamin,etamax,0.,1.);
00281   
00282   // histos titles
00283   h_mcNum              -> GetXaxis()-> SetTitle("# MC particles");
00284   h_eleNum             -> GetXaxis()-> SetTitle("# MC ele");
00285   h_gamNum             -> GetXaxis()-> SetTitle("# MC gammas");
00286   h_simEta             -> GetXaxis()-> SetTitle("true #eta");
00287   h_simP               -> GetXaxis()-> SetTitle("true p (GeV/c)");
00288   h_ele_foundHits      -> GetXaxis()-> SetTitle("# hits");   
00289   h_ele_lostHits       -> GetXaxis()-> SetTitle("# lost hits");   
00290   h_ele_chi2           -> GetXaxis()-> SetTitle("#Chi^{2}");   
00291   h_ele_charge         -> GetXaxis()-> SetTitle("charge");   
00292   h_ele_vertexP        -> GetXaxis()-> SetTitle("p_{vertex} (GeV/c)");
00293   h_ele_vertexPt       -> GetXaxis()-> SetTitle("p_{T vertex} (GeV/c)");
00294   h_ele_vertexEta      -> GetXaxis()-> SetTitle("#eta");  
00295   h_ele_vertexPhi      -> GetXaxis()-> SetTitle("#phi");   
00296   h_ele_PoPtrue        -> GetXaxis()-> SetTitle("P/P_{true}");
00297   h_ele_EtaMnEtaTrue   -> GetXaxis()-> SetTitle("#eta_{rec} - #eta_{true}");
00298   h_ele_PhiMnPhiTrue   -> GetXaxis()-> SetTitle("#phi_{rec} - #phi_{true}");
00299   h_ele_PinMnPout      -> GetXaxis()-> SetTitle("p_{vertex} - p_{out} (GeV)");
00300   h_ele_PinMnPout_mode      -> GetXaxis()-> SetTitle("p_{vertex} - p_{out}, mode (GeV)");
00301   h_ele_outerP         -> GetXaxis()-> SetTitle("p_{out} (GeV/c)");
00302   h_ele_outerP_mode         -> GetXaxis()-> SetTitle("p_{out} (GeV/c)");
00303   h_ele_outerPt        -> GetXaxis()-> SetTitle("p_{T out} (GeV/c)");
00304   h_ele_outerPt_mode        -> GetXaxis()-> SetTitle("p_{T out} (GeV/c)");
00305   h_ele_EoP            -> GetXaxis()-> SetTitle("E/p_{vertex}");
00306   h_ele_EoPout         -> GetXaxis()-> SetTitle("E/p_{out}");
00307 
00308 }     
00309 
00310 void
00311 GsfElectronMCAnalyzer::endJob(){
00312   
00313   histfile_->cd();
00314   std::cout << "efficiency calculation " << std::endl; 
00315   // efficiency vs eta
00316   TH1F *h_ele_etaEff = (TH1F*)h_ele_simEta_matched->Clone("h_ele_etaEff");
00317   h_ele_etaEff->Reset();
00318   h_ele_etaEff->Divide(h_ele_simEta_matched,h_simEta,1,1);
00319   h_ele_etaEff->Print();
00320   h_ele_etaEff->GetXaxis()->SetTitle("#eta");
00321   h_ele_etaEff->GetYaxis()->SetTitle("eff");
00322   
00323   // efficiency vs z
00324   TH1F *h_ele_zEff = (TH1F*)h_ele_simZ_matched->Clone("h_ele_zEff");
00325   h_ele_zEff->Reset();
00326   h_ele_zEff->Divide(h_ele_simZ_matched,h_simZ,1,1);
00327   h_ele_zEff->Print();
00328   h_ele_zEff->GetXaxis()->SetTitle("z");
00329   h_ele_zEff->GetYaxis()->SetTitle("eff");
00330 
00331   // efficiency vs |eta|
00332   TH1F *h_ele_absetaEff = (TH1F*)h_ele_simAbsEta_matched->Clone("h_ele_absetaEff");
00333   h_ele_absetaEff->Reset();
00334   h_ele_absetaEff->Divide(h_ele_simAbsEta_matched,h_simAbsEta,1,1);
00335   h_ele_absetaEff->GetXaxis()->SetTitle("|#eta|");
00336   h_ele_absetaEff->GetYaxis()->SetTitle("eff");
00337 
00338   // efficiency vs pt
00339   TH1F *h_ele_ptEff = (TH1F*)h_ele_simPt_matched->Clone("h_ele_ptEff");
00340   h_ele_ptEff->Reset();
00341   h_ele_ptEff->Divide(h_ele_simPt_matched,h_simPt,1,1);
00342   h_ele_ptEff->GetXaxis()->SetTitle("p_T");
00343   h_ele_ptEff->GetYaxis()->SetTitle("eff");
00344 
00345   // efficiency vs phi
00346   TH1F *h_ele_phiEff = (TH1F*)h_ele_simPhi_matched->Clone("h_ele_phiEff");
00347   h_ele_phiEff->Reset();
00348   h_ele_phiEff->Divide(h_ele_simPhi_matched,h_simPhi,1,1);
00349   h_ele_phiEff->GetXaxis()->SetTitle("phi");
00350   h_ele_phiEff->GetYaxis()->SetTitle("eff");
00351 
00352   // rec/gen all electrons
00353   TH1F *h_ele_etaEff_all = (TH1F*)h_ele_vertexEta_all->Clone("h_ele_etaEff_all");
00354   h_ele_etaEff_all->Reset();
00355   h_ele_etaEff_all->Divide(h_ele_vertexEta_all,h_simEta,1,1);
00356   h_ele_etaEff_all->Print();
00357   h_ele_etaEff_all->GetXaxis()->SetTitle("#eta");
00358   h_ele_etaEff_all->GetYaxis()->SetTitle("rec/gen");
00359   TH1F *h_ele_ptEff_all = (TH1F*)h_ele_vertexPt_all->Clone("h_ele_ptEff_all");
00360   h_ele_ptEff_all->Reset();
00361   h_ele_ptEff_all->Divide(h_ele_vertexPt_all,h_simPt,1,1);
00362   h_ele_ptEff_all->Print();
00363   h_ele_ptEff_all->GetXaxis()->SetTitle("p_{T}");
00364   h_ele_ptEff_all->GetYaxis()->SetTitle("rec/gen");
00365 
00366   // classes
00367   TH1F *h_ele_eta_goldenFrac = (TH1F*)h_ele_eta_golden->Clone("h_ele_eta_goldenFrac");
00368   h_ele_eta_goldenFrac->Reset();
00369   h_ele_eta_goldenFrac->Divide(h_ele_eta_golden,h_ele_eta,1,1);
00370   TH1F *h_ele_eta_bbremFrac = (TH1F*)h_ele_eta_bbrem->Clone("h_ele_eta_bbremFrac");
00371   h_ele_eta_bbremFrac->Reset();
00372   h_ele_eta_bbremFrac->Divide(h_ele_eta_bbrem,h_ele_eta,1,1);
00373   TH1F *h_ele_eta_narrowFrac = (TH1F*)h_ele_eta_narrow->Clone("h_ele_eta_narrowFrac");
00374   h_ele_eta_narrowFrac->Reset();
00375   h_ele_eta_narrowFrac->Divide(h_ele_eta_narrow,h_ele_eta,1,1);
00376   TH1F *h_ele_eta_showerFrac = (TH1F*)h_ele_eta_shower->Clone("h_ele_eta_showerFrac");
00377   h_ele_eta_showerFrac->Reset();
00378   h_ele_eta_showerFrac->Divide(h_ele_eta_shower,h_ele_eta,1,1);
00379   
00380   // fbrem
00381   TH1F *h_ele_xOverX0VsEta = new TH1F( "h_ele_xOverx0VsEta","mean X/X_0 vs eta",nbineta/2,0.0,2.5);
00382   for (int ibin=1;ibin<h_ele_fbremVsEta_mean->GetNbinsX()+1;ibin++) {
00383     double xOverX0 = 0.;
00384     if (h_ele_fbremVsEta_mean->GetBinContent(ibin)>0.) xOverX0 = -log(h_ele_fbremVsEta_mean->GetBinContent(ibin));
00385     h_ele_xOverX0VsEta->SetBinContent(ibin,xOverX0);
00386   }
00387    
00388   //profiles from 2D histos
00389   TProfile *p_ele_PoPtrueVsEta = h_ele_PoPtrueVsEta->ProfileX();
00390   p_ele_PoPtrueVsEta->Write();
00391   TProfile *p_ele_PoPtrueVsPhi = h_ele_PoPtrueVsPhi->ProfileX();
00392   p_ele_PoPtrueVsPhi->Write();
00393   TProfile *p_ele_EtaMnEtaTrueVsEta = h_ele_EtaMnEtaTrueVsEta->ProfileX();
00394   p_ele_EtaMnEtaTrueVsEta->Write();
00395   TProfile *p_ele_EtaMnEtaTrueVsPhi = h_ele_EtaMnEtaTrueVsPhi->ProfileX();
00396   p_ele_EtaMnEtaTrueVsPhi->Write();
00397   TProfile *p_ele_PhiMnPhiTrueVsEta = h_ele_PhiMnPhiTrueVsEta->ProfileX();
00398   p_ele_PhiMnPhiTrueVsEta->Write();
00399   TProfile *p_ele_PhiMnPhiTrueVsPhi = h_ele_PhiMnPhiTrueVsPhi->ProfileX();
00400   p_ele_PhiMnPhiTrueVsPhi->Write();
00401   TProfile *p_ele_vertexPtVsEta = h_ele_vertexPtVsEta->ProfileX();
00402   p_ele_vertexPtVsEta->Write();
00403   TProfile *p_ele_vertexPtVsPhi = h_ele_vertexPtVsPhi->ProfileX();
00404   p_ele_vertexPtVsPhi->Write();
00405   TProfile *p_ele_EoPVsEta = h_ele_EoPVsEta->ProfileX();
00406   p_ele_EoPVsEta->Write();
00407   TProfile *p_ele_EoPVsPhi = h_ele_EoPVsPhi->ProfileX();
00408   p_ele_EoPVsPhi->Write();
00409   TProfile *p_ele_EoPoutVsEta = h_ele_EoPoutVsEta->ProfileX();
00410   p_ele_EoPoutVsEta->Write();
00411   TProfile *p_ele_EoPoutVsPhi = h_ele_EoPoutVsPhi->ProfileX();
00412   p_ele_EoPoutVsPhi->Write();
00413   TProfile *p_ele_HoEVsEta = h_ele_HoEVsEta->ProfileX();
00414   p_ele_HoEVsEta->Write();
00415   TProfile *p_ele_HoEVsPhi = h_ele_HoEVsPhi->ProfileX();
00416   p_ele_HoEVsPhi->Write();
00417   TProfile *p_ele_chi2VsEta = h_ele_chi2VsEta->ProfileX();
00418   p_ele_chi2VsEta->Write();
00419   TProfile *p_ele_chi2VsPhi = h_ele_chi2VsPhi->ProfileX();
00420   p_ele_chi2VsPhi->Write();
00421   TProfile *p_ele_foundHitsVsEta = h_ele_foundHitsVsEta->ProfileX();
00422   p_ele_foundHitsVsEta->Write();
00423   TProfile *p_ele_foundHitsVsPhi = h_ele_foundHitsVsPhi->ProfileX();
00424   p_ele_foundHitsVsPhi->Write();
00425   TProfile *p_ele_lostHitsVsEta = h_ele_lostHitsVsEta->ProfileX();
00426   p_ele_lostHitsVsEta->Write();
00427   TProfile *p_ele_lostHitsVsPhi = h_ele_lostHitsVsPhi->ProfileX();
00428   p_ele_lostHitsVsPhi->Write();
00429   
00430   // mc truth  
00431 
00432   h_mcNum->Write();
00433   h_eleNum->Write();
00434   h_gamNum->Write();
00435     
00436   // rec event
00437   
00438   histNum_->Write();
00439   
00440   // mc  
00441   h_simEta->Write();
00442   h_simAbsEta->Write();
00443   h_simP->Write();
00444   h_simPt->Write();
00445 
00446   // ctf tracks
00447   h_ctf_foundHitsVsEta->Write();
00448   h_ctf_lostHitsVsEta->Write();
00449   
00450   // all electrons  
00451   h_ele_EoverP_all->Write();
00452   h_ele_TIP_all->Write();
00453   h_ele_vertexPt_all->Write();
00454   h_ele_vertexEta_all->Write();
00455   h_ele_mee_all->Write();
00456 
00457   // matched electrons
00458   h_ele_charge->Write();
00459   h_ele_chargeVsEta->Write();
00460   h_ele_chargeVsPhi->Write();
00461   h_ele_chargeVsPt->Write();
00462   h_ele_vertexP->Write();
00463   h_ele_vertexPt->Write();
00464   h_ele_vertexPtVsEta->Write();
00465   h_ele_vertexPtVsPhi->Write();
00466   h_ele_simPt_matched->Write();
00467   h_ele_vertexEta->Write();
00468   h_ele_vertexEtaVsPhi->Write();
00469   h_ele_simAbsEta_matched->Write();
00470   h_ele_simEta_matched->Write();
00471   h_ele_vertexPhi->Write();
00472   h_ele_vertexX->Write();
00473   h_ele_vertexY ->Write();
00474   h_ele_vertexZ->Write();
00475   h_ele_vertexTIP->Write();
00476   h_ele_simZ_matched->Write();
00477   h_ele_vertexTIPVsEta->Write();
00478   h_ele_vertexTIPVsPhi->Write();
00479   h_ele_vertexTIPVsPt->Write();
00480   h_ele_PoPtrue->Write();
00481   h_ele_PoPtrueVsEta ->Write();
00482   h_ele_PoPtrueVsPhi->Write();
00483   h_ele_PoPtrueVsPt->Write();
00484   h_ele_PoPtrue_barrel ->Write();
00485   h_ele_PoPtrue_endcaps->Write();
00486   h_ele_EtaMnEtaTrue->Write();
00487   h_ele_EtaMnEtaTrueVsEta ->Write();
00488   h_ele_EtaMnEtaTrueVsPhi->Write();
00489   h_ele_EtaMnEtaTrueVsPt->Write();
00490   h_ele_PhiMnPhiTrue ->Write();
00491   h_ele_PhiMnPhiTrue2 ->Write();
00492   h_ele_PhiMnPhiTrueVsEta->Write();
00493   h_ele_PhiMnPhiTrueVsPhi->Write();
00494   h_ele_PhiMnPhiTrueVsPt->Write();
00495 
00496   // matched electron, superclusters
00497   histSclEn_->Write();
00498   histSclEoEtrue_barrel->Write();
00499   histSclEoEtrue_endcaps->Write();
00500   histSclEt_->Write();
00501   histSclEtVsEta_->Write();
00502   histSclEtVsPhi_->Write();
00503   histSclEtaVsPhi_ ->Write();
00504   histSclEta_->Write();
00505   histSclPhi_->Write();
00506 
00507   // matched electron, gsf tracks
00508   h_ele_foundHits->Write();
00509   h_ele_foundHitsVsEta->Write();
00510   h_ele_foundHitsVsPhi->Write();
00511   h_ele_foundHitsVsPt->Write();
00512   h_ctf_foundHits->Write();
00513   h_ele_lostHits->Write();
00514   h_ele_lostHitsVsEta->Write();
00515   h_ele_lostHitsVsPhi->Write();
00516   h_ele_lostHitsVsPt->Write();
00517   h_ele_chi2 ->Write();
00518   h_ele_chi2VsEta ->Write();
00519   h_ele_chi2VsPhi ->Write();
00520   h_ele_chi2VsPt->Write();
00521   h_ele_PinMnPout->Write();
00522   h_ele_PinMnPout_mode->Write();
00523   h_ele_PinMnPoutVsEta_mode->Write();
00524   h_ele_PinMnPoutVsPhi_mode->Write();
00525   h_ele_PinMnPoutVsPt_mode->Write();
00526   h_ele_PinMnPoutVsE_mode->Write();
00527   h_ele_PinMnPoutVsChi2_mode->Write();
00528   h_ele_outerP ->Write();
00529   h_ele_outerP_mode->Write();
00530   h_ele_outerPVsEta_mode->Write();
00531   h_ele_outerPt->Write();
00532   h_ele_outerPt_mode ->Write();
00533   h_ele_outerPtVsEta_mode->Write();
00534   h_ele_outerPtVsPhi_mode->Write();
00535   h_ele_outerPtVsPt_mode->Write();
00536   
00537   // matched electrons, matching 
00538   h_ele_EoP ->Write();
00539   h_ele_EoPVsEta ->Write();
00540   h_ele_EoPVsPhi->Write();
00541   h_ele_EoPVsE->Write();
00542   h_ele_EoPout->Write();
00543   h_ele_EoPoutVsEta->Write();
00544   h_ele_EoPoutVsPhi->Write();
00545   h_ele_EoPoutVsE ->Write();
00546   h_ele_dEtaSc_propVtx->Write();
00547   h_ele_dEtaScVsEta_propVtx->Write();
00548   h_ele_dEtaScVsPhi_propVtx->Write();
00549   h_ele_dEtaScVsPt_propVtx ->Write();
00550   h_ele_dPhiSc_propVtx->Write();
00551   h_ele_dPhiScVsEta_propVtx ->Write();
00552   h_ele_dPhiScVsPhi_propVtx->Write();
00553   h_ele_dPhiScVsPt_propVtx->Write();
00554   h_ele_dEtaCl_propOut->Write();
00555   h_ele_dEtaClVsEta_propOut->Write();
00556   h_ele_dEtaClVsPhi_propOut->Write();
00557   h_ele_dEtaClVsPt_propOut->Write();
00558   h_ele_dPhiCl_propOut->Write();
00559   h_ele_dPhiClVsEta_propOut->Write();
00560   h_ele_dPhiClVsPhi_propOut->Write();
00561   h_ele_dPhiClVsPt_propOut->Write();
00562   
00563   h_ele_HoE->Write();
00564   h_ele_HoEVsEta->Write();
00565   h_ele_HoEVsPhi->Write();
00566   h_ele_HoEVsE->Write();
00567  
00568   // classes  
00569   h_ele_classes->Write();
00570   h_ele_eta->Write();
00571   h_ele_eta_golden->Write();
00572   h_ele_eta_bbrem->Write();
00573   h_ele_eta_narrow->Write();
00574   h_ele_eta_shower->Write();
00575   h_ele_PinVsPoutGolden_mode->Write();
00576   h_ele_PinVsPoutShowering0_mode->Write();
00577   h_ele_PinVsPoutShowering1234_mode->Write();
00578   h_ele_PinVsPoutGolden_mean->Write();
00579   h_ele_PinVsPoutShowering0_mean->Write();
00580   h_ele_PinVsPoutShowering1234_mean->Write();
00581   h_ele_PtinVsPtoutGolden_mode->Write();
00582   h_ele_PtinVsPtoutShowering0_mode->Write();
00583   h_ele_PtinVsPtoutShowering1234_mode->Write();
00584   h_ele_PtinVsPtoutGolden_mean->Write();
00585   h_ele_PtinVsPtoutShowering0_mean->Write();
00586   h_ele_PtinVsPtoutShowering1234_mean->Write();
00587   histSclEoEtrueGolden_barrel->Write();
00588   histSclEoEtrueGolden_endcaps->Write();
00589   histSclEoEtrueShowering0_barrel->Write();
00590   histSclEoEtrueShowering0_endcaps->Write();
00591   histSclEoEtrueShowering1234_barrel->Write();
00592   histSclEoEtrueShowering1234_endcaps->Write();
00593 
00594   // fbrem
00595   h_ele_fbremVsEta_mode->Write();
00596   h_ele_fbremVsEta_mean->Write();
00597   h_ele_etaEff->Write();
00598   h_ele_zEff->Write();
00599   h_ele_phiEff->Write();
00600   h_ele_absetaEff->Write();
00601   h_ele_ptEff->Write();
00602   h_ele_etaEff_all->Write();
00603   h_ele_ptEff_all->Write();
00604   h_ele_eta_goldenFrac->Write();
00605   h_ele_eta_bbremFrac->Write();
00606   h_ele_eta_narrowFrac->Write();
00607   h_ele_eta_showerFrac->Write();
00608   h_ele_xOverX0VsEta->Write();
00609   
00610 }
00611 
00612 void
00613 GsfElectronMCAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
00614 {
00615   std::cout << "analyzing new event " << std::endl;
00616   // get electrons
00617   
00618   edm::Handle<GsfElectronCollection> gsfElectrons;
00619   iEvent.getByLabel(electronCollection_,gsfElectrons); 
00620   edm::LogInfo("")<<"\n\n =================> Treating event "<<iEvent.id()<<" Number of electrons "<<gsfElectrons.product()->size();
00621 
00622   edm::Handle<edm::HepMCProduct> hepMC;
00623   iEvent.getByLabel(mcTruthCollection_,hepMC);
00624 
00625   histNum_->Fill((*gsfElectrons).size());
00626   
00627   // all rec electrons
00628   for (reco::GsfElectronCollection::const_iterator gsfIter=gsfElectrons->begin();
00629    gsfIter!=gsfElectrons->end(); gsfIter++){
00630     // preselect electrons
00631     if (gsfIter->pt()>maxPt_ || fabs(gsfIter->eta())>maxAbsEta_) continue;
00632     double d = gsfIter->gsfTrack()->vertex().x()*gsfIter->gsfTrack()->vertex().x()+
00633      gsfIter->gsfTrack()->vertex().y()*gsfIter->gsfTrack()->vertex().y();
00634     h_ele_TIP_all     -> Fill( sqrt(d) );
00635     h_ele_EoverP_all     -> Fill( gsfIter->eSuperClusterOverP() );
00636     h_ele_vertexEta_all     -> Fill( gsfIter->eta() );
00637     h_ele_vertexPt_all      -> Fill( gsfIter->pt() );
00638     // mee
00639     for (reco::GsfElectronCollection::const_iterator gsfIter2=gsfIter+1;
00640      gsfIter2!=gsfElectrons->end(); gsfIter2++){
00641         math::XYZTLorentzVector p12 = (*gsfIter).p4()+(*gsfIter2).p4();
00642         float mee2 = p12.Dot(p12);
00643         h_ele_mee_all -> Fill(sqrt(mee2));       
00644     }
00645   }
00646    
00647   // association mc-reco
00648   HepMC::GenParticle* genPc=0;
00649   const HepMC::GenEvent *myGenEvent = hepMC->GetEvent();
00650   int mcNum=0, gamNum=0, eleNum=0;
00651   HepMC::FourVector pAssSim;
00652       
00653   for ( HepMC::GenEvent::particle_const_iterator mcIter=myGenEvent->particles_begin(); mcIter != myGenEvent->particles_end(); mcIter++ ) {
00654     
00655     // number of mc particles
00656     mcNum++;
00657 
00658     // counts photons
00659     if ((*mcIter)->pdg_id() == 22 ){ gamNum++; }       
00660 
00661     // select electrons
00662     if ( (*mcIter)->pdg_id() == 11 || (*mcIter)->pdg_id() == -11 ){       
00663 
00664       // single primary electrons or electrons from Zs or Ws
00665       HepMC::GenParticle* mother = 0;
00666       if ( (*mcIter)->production_vertex() )  {
00667        if ( (*mcIter)->production_vertex()->particles_begin(HepMC::parents) != 
00668            (*mcIter)->production_vertex()->particles_end(HepMC::parents))  
00669             mother = *((*mcIter)->production_vertex()->particles_begin(HepMC::parents));
00670       } 
00671       if ( ((mother == 0) || ((mother != 0) && (mother->pdg_id() == 23))
00672                           || ((mother != 0) && (mother->pdg_id() == 32))
00673                           || ((mother != 0) && (fabs(mother->pdg_id()) == 24)))) {       
00674    
00675       genPc=(*mcIter);
00676       pAssSim = genPc->momentum();
00677 
00678       if (pAssSim.perp()> maxPt_ || fabs(pAssSim.eta())> maxAbsEta_) continue;
00679       
00680       // suppress the endcaps
00681       //if (fabs(pAssSim.eta()) > 1.5) continue;
00682       // select central z
00683       //if ( fabs((*mcIter)->production_vertex()->position().z())>50.) continue;
00684  
00685       eleNum++;
00686       h_simEta -> Fill( pAssSim.eta() );
00687       h_simAbsEta -> Fill( fabs(pAssSim.eta()) );
00688       h_simP   -> Fill( pAssSim.t() );
00689       h_simPt   -> Fill( pAssSim.perp() );
00690       h_simPhi   -> Fill( pAssSim.phi() );
00691       h_simZ   -> Fill( (*mcIter)->production_vertex()->position().z()/10. );
00692         
00693       // looking for the best matching gsf electron
00694       bool okGsfFound = false;
00695       double gsfOkRatio = 999999.;
00696 
00697       // find best matched electron
00698       reco::GsfElectron bestGsfElectron;
00699       for (reco::GsfElectronCollection::const_iterator gsfIter=gsfElectrons->begin();
00700        gsfIter!=gsfElectrons->end(); gsfIter++){
00701         
00702         double deltaR = sqrt(pow((gsfIter->eta()-pAssSim.eta()),2) + pow((gsfIter->phi()-pAssSim.phi()),2));
00703         if ( deltaR < deltaR_ ){
00704         if ( (genPc->pdg_id() == 11) && (gsfIter->charge() < 0.) || (genPc->pdg_id() == -11) &&
00705         (gsfIter->charge() > 0.) ){
00706           double tmpGsfRatio = gsfIter->p()/pAssSim.t();
00707           if ( fabs(tmpGsfRatio-1) < fabs(gsfOkRatio-1) ) {
00708             gsfOkRatio = tmpGsfRatio;
00709             bestGsfElectron=*gsfIter;
00710             okGsfFound = true;
00711           } 
00712         } 
00713         } 
00714       } // loop over rec ele to look for the best one   
00715 
00716       // analysis when the mc track is found
00717      if (okGsfFound){
00718 
00719         // electron related distributions
00720         h_ele_charge        -> Fill( bestGsfElectron.charge() );
00721         h_ele_chargeVsEta        -> Fill( bestGsfElectron.eta(),bestGsfElectron.charge() );
00722         h_ele_chargeVsPhi        -> Fill( bestGsfElectron.phi(),bestGsfElectron.charge() );
00723         h_ele_chargeVsPt        -> Fill( bestGsfElectron.pt(),bestGsfElectron.charge() );
00724         h_ele_vertexP       -> Fill( bestGsfElectron.p() );
00725         h_ele_vertexPt      -> Fill( bestGsfElectron.pt() );
00726         h_ele_vertexPtVsEta      -> Fill(  bestGsfElectron.eta(),bestGsfElectron.pt() );
00727         h_ele_vertexPtVsPhi      -> Fill(  bestGsfElectron.phi(),bestGsfElectron.pt() );
00728         h_ele_vertexEta     -> Fill( bestGsfElectron.eta() );
00729         // generated distributions for matched electrons
00730         h_ele_simPt_matched      -> Fill( pAssSim.perp() );
00731         h_ele_simPhi_matched   -> Fill( pAssSim.phi() );
00732         h_ele_simAbsEta_matched     -> Fill( fabs(pAssSim.eta()) );
00733         h_ele_simEta_matched     -> Fill( pAssSim.eta() );
00734         h_ele_vertexEtaVsPhi     -> Fill(  bestGsfElectron.phi(),bestGsfElectron.eta() );
00735         h_ele_vertexPhi     -> Fill( bestGsfElectron.phi() );
00736         h_ele_vertexX     -> Fill( bestGsfElectron.vertex().x() );
00737         h_ele_vertexY     -> Fill( bestGsfElectron.vertex().y() );
00738         h_ele_vertexZ     -> Fill( bestGsfElectron.vertex().z() );
00739         h_ele_simZ_matched   -> Fill( (*mcIter)->production_vertex()->position().z()/10. );
00740         double d = bestGsfElectron.gsfTrack()->vertex().x()*bestGsfElectron.gsfTrack()->vertex().x()+
00741          bestGsfElectron.gsfTrack()->vertex().y()*bestGsfElectron.gsfTrack()->vertex().y();
00742         d = sqrt(d); 
00743         h_ele_vertexTIP     -> Fill( d );
00744         h_ele_vertexTIPVsEta     -> Fill(  bestGsfElectron.eta(), d );
00745         h_ele_vertexTIPVsPhi     -> Fill(  bestGsfElectron.phi(), d );
00746         h_ele_vertexTIPVsPt     -> Fill(  bestGsfElectron.pt(), d );    
00747         h_ele_EtaMnEtaTrue  -> Fill( bestGsfElectron.eta()-pAssSim.eta());
00748         h_ele_EtaMnEtaTrueVsEta  -> Fill( bestGsfElectron.eta(), bestGsfElectron.eta()-pAssSim.eta());
00749         h_ele_EtaMnEtaTrueVsPhi  -> Fill( bestGsfElectron.phi(), bestGsfElectron.eta()-pAssSim.eta());
00750         h_ele_EtaMnEtaTrueVsPt  -> Fill( bestGsfElectron.pt(), bestGsfElectron.eta()-pAssSim.eta());
00751         h_ele_PhiMnPhiTrue  -> Fill( bestGsfElectron.phi()-pAssSim.phi());
00752         h_ele_PhiMnPhiTrue2  -> Fill( bestGsfElectron.phi()-pAssSim.phi());
00753         h_ele_PhiMnPhiTrueVsEta  -> Fill( bestGsfElectron.eta(), bestGsfElectron.phi()-pAssSim.phi());
00754         h_ele_PhiMnPhiTrueVsPhi  -> Fill( bestGsfElectron.phi(), bestGsfElectron.phi()-pAssSim.phi());
00755         h_ele_PhiMnPhiTrueVsPt  -> Fill( bestGsfElectron.pt(), bestGsfElectron.phi()-pAssSim.phi());
00756         h_ele_PoPtrue       -> Fill( bestGsfElectron.p()/pAssSim.t());
00757         h_ele_PoPtrueVsEta       -> Fill( bestGsfElectron.eta(), bestGsfElectron.p()/pAssSim.t());
00758         h_ele_PoPtrueVsPhi       -> Fill( bestGsfElectron.phi(), bestGsfElectron.p()/pAssSim.t());
00759         h_ele_PoPtrueVsPt       -> Fill( bestGsfElectron.py(), bestGsfElectron.p()/pAssSim.t());
00760         if (bestGsfElectron.classification() < 100) h_ele_PoPtrue_barrel       -> Fill( bestGsfElectron.p()/pAssSim.t());
00761         if (bestGsfElectron.classification() >= 100) h_ele_PoPtrue_endcaps       -> Fill( bestGsfElectron.p()/pAssSim.t());
00762 
00763         // supercluster related distributions
00764         reco::SuperClusterRef sclRef = bestGsfElectron.superCluster();
00765         histSclEn_->Fill(sclRef->energy());
00766         double R=TMath::Sqrt(sclRef->x()*sclRef->x() + sclRef->y()*sclRef->y() +sclRef->z()*sclRef->z());
00767         double Rt=TMath::Sqrt(sclRef->x()*sclRef->x() + sclRef->y()*sclRef->y());
00768         histSclEt_->Fill(sclRef->energy()*(Rt/R));
00769         histSclEtVsEta_->Fill(sclRef->eta(),sclRef->energy()*(Rt/R));
00770         histSclEtVsPhi_->Fill(sclRef->phi(),sclRef->energy()*(Rt/R));
00771         if (bestGsfElectron.classification() < 100)  histSclEoEtrue_barrel->Fill(sclRef->energy()/pAssSim.t());
00772         if (bestGsfElectron.classification() >= 100)  histSclEoEtrue_endcaps->Fill(sclRef->energy()/pAssSim.t());
00773         histSclEta_->Fill(sclRef->eta());
00774         histSclEtaVsPhi_->Fill(sclRef->phi(),sclRef->eta());
00775         histSclPhi_->Fill(sclRef->phi());
00776 
00777         // track related distributions
00778         h_ele_foundHits     -> Fill( bestGsfElectron.gsfTrack()->numberOfValidHits() );
00779         h_ele_foundHitsVsEta     -> Fill( bestGsfElectron.eta(), bestGsfElectron.gsfTrack()->numberOfValidHits() );
00780         h_ele_foundHitsVsPhi     -> Fill( bestGsfElectron.phi(), bestGsfElectron.gsfTrack()->numberOfValidHits() );
00781         h_ele_foundHitsVsPt     -> Fill( bestGsfElectron.pt(), bestGsfElectron.gsfTrack()->numberOfValidHits() );
00782         h_ele_lostHits      -> Fill( bestGsfElectron.gsfTrack()->numberOfLostHits() );
00783         h_ele_lostHitsVsEta      -> Fill( bestGsfElectron.eta(), bestGsfElectron.gsfTrack()->numberOfLostHits() );
00784         h_ele_lostHitsVsPhi      -> Fill( bestGsfElectron.phi(), bestGsfElectron.gsfTrack()->numberOfLostHits() );
00785         h_ele_lostHitsVsPt      -> Fill( bestGsfElectron.pt(), bestGsfElectron.gsfTrack()->numberOfLostHits() );
00786         h_ele_chi2          -> Fill( bestGsfElectron.gsfTrack()->normalizedChi2() );  
00787         h_ele_chi2VsEta          -> Fill( bestGsfElectron.eta(), bestGsfElectron.gsfTrack()->normalizedChi2() );  
00788         h_ele_chi2VsPhi          -> Fill( bestGsfElectron.phi(), bestGsfElectron.gsfTrack()->normalizedChi2() );  
00789         h_ele_chi2VsPt          -> Fill( bestGsfElectron.pt(), bestGsfElectron.gsfTrack()->normalizedChi2() );  
00790         // from gsf track interface, hence using mean
00791         h_ele_PinMnPout     -> Fill( bestGsfElectron.gsfTrack()->innerMomentum().R() - bestGsfElectron.gsfTrack()->outerMomentum().R() );
00792         h_ele_outerP        -> Fill( bestGsfElectron.gsfTrack()->outerMomentum().R() );
00793         h_ele_outerPt       -> Fill( bestGsfElectron.gsfTrack()->outerMomentum().Rho() );
00794         // from electron interface, hence using mode
00795         h_ele_PinMnPout_mode     -> Fill( bestGsfElectron.trackMomentumAtVtx().R() - bestGsfElectron.trackMomentumOut().R() );
00796         h_ele_PinMnPoutVsEta_mode     -> Fill(  bestGsfElectron.eta(), bestGsfElectron.trackMomentumAtVtx().R() - bestGsfElectron.trackMomentumOut().R() );
00797         h_ele_PinMnPoutVsPhi_mode     -> Fill(  bestGsfElectron.phi(), bestGsfElectron.trackMomentumAtVtx().R() - bestGsfElectron.trackMomentumOut().R() );
00798         h_ele_PinMnPoutVsPt_mode     -> Fill(  bestGsfElectron.pt(), bestGsfElectron.trackMomentumAtVtx().R() - bestGsfElectron.trackMomentumOut().R() );
00799         h_ele_PinMnPoutVsE_mode     -> Fill(  bestGsfElectron.caloEnergy(), bestGsfElectron.trackMomentumAtVtx().R() - bestGsfElectron.trackMomentumOut().R() );
00800         h_ele_PinMnPoutVsChi2_mode     -> Fill(  bestGsfElectron.gsfTrack()->normalizedChi2(), bestGsfElectron.trackMomentumAtVtx().R() - bestGsfElectron.trackMomentumOut().R() );
00801         h_ele_outerP_mode        -> Fill( bestGsfElectron.trackMomentumOut().R() );
00802         h_ele_outerPVsEta_mode        -> Fill(bestGsfElectron.eta(),  bestGsfElectron.trackMomentumOut().R() );
00803         h_ele_outerPt_mode       -> Fill( bestGsfElectron.trackMomentumOut().Rho() );
00804         h_ele_outerPtVsEta_mode       -> Fill(bestGsfElectron.eta(),  bestGsfElectron.trackMomentumOut().Rho() );
00805         h_ele_outerPtVsPhi_mode       -> Fill(bestGsfElectron.phi(),  bestGsfElectron.trackMomentumOut().Rho() );
00806         h_ele_outerPtVsPt_mode       -> Fill(bestGsfElectron.pt(),  bestGsfElectron.trackMomentumOut().Rho() );
00807         
00808         // match distributions 
00809         h_ele_EoP    -> Fill( bestGsfElectron.eSuperClusterOverP() );
00810         h_ele_EoPVsEta    -> Fill(bestGsfElectron.eta(),  bestGsfElectron.eSuperClusterOverP() );
00811         h_ele_EoPVsPhi    -> Fill(bestGsfElectron.phi(),  bestGsfElectron.eSuperClusterOverP() );
00812         h_ele_EoPVsE    -> Fill(bestGsfElectron.caloEnergy(),  bestGsfElectron.eSuperClusterOverP() );
00813         h_ele_EoPout -> Fill( bestGsfElectron.eSeedClusterOverPout() );
00814         h_ele_EoPoutVsEta -> Fill( bestGsfElectron.eta(), bestGsfElectron.eSeedClusterOverPout() );
00815         h_ele_EoPoutVsPhi -> Fill( bestGsfElectron.phi(), bestGsfElectron.eSeedClusterOverPout() );
00816         h_ele_EoPoutVsE -> Fill( bestGsfElectron.caloEnergy(), bestGsfElectron.eSeedClusterOverPout() );
00817         h_ele_dEtaSc_propVtx -> Fill(bestGsfElectron.deltaEtaSuperClusterTrackAtVtx());
00818         h_ele_dEtaScVsEta_propVtx -> Fill( bestGsfElectron.eta(),bestGsfElectron.deltaEtaSuperClusterTrackAtVtx());
00819         h_ele_dEtaScVsPhi_propVtx -> Fill(bestGsfElectron.phi(),bestGsfElectron.deltaEtaSuperClusterTrackAtVtx());
00820         h_ele_dEtaScVsPt_propVtx -> Fill(bestGsfElectron.pt(),bestGsfElectron.deltaEtaSuperClusterTrackAtVtx());
00821         h_ele_dPhiSc_propVtx -> Fill(bestGsfElectron.deltaPhiSuperClusterTrackAtVtx()); 
00822         h_ele_dPhiScVsEta_propVtx -> Fill( bestGsfElectron.eta(),bestGsfElectron.deltaPhiSuperClusterTrackAtVtx()); 
00823         h_ele_dPhiScVsPhi_propVtx -> Fill(bestGsfElectron.phi(),bestGsfElectron.deltaPhiSuperClusterTrackAtVtx()); 
00824         h_ele_dPhiScVsPt_propVtx -> Fill(bestGsfElectron.pt(),bestGsfElectron.deltaPhiSuperClusterTrackAtVtx()); 
00825         h_ele_dEtaCl_propOut -> Fill(bestGsfElectron.deltaEtaSeedClusterTrackAtCalo()); 
00826         h_ele_dEtaClVsEta_propOut -> Fill( bestGsfElectron.eta(),bestGsfElectron.deltaEtaSeedClusterTrackAtCalo()); 
00827         h_ele_dEtaClVsPhi_propOut -> Fill(bestGsfElectron.phi(),bestGsfElectron.deltaEtaSeedClusterTrackAtCalo()); 
00828         h_ele_dEtaClVsPt_propOut -> Fill(bestGsfElectron.pt(),bestGsfElectron.deltaEtaSeedClusterTrackAtCalo()); 
00829         h_ele_dPhiCl_propOut -> Fill(bestGsfElectron.deltaPhiSeedClusterTrackAtCalo()); 
00830         h_ele_dPhiClVsEta_propOut -> Fill( bestGsfElectron.eta(),bestGsfElectron.deltaPhiSeedClusterTrackAtCalo()); 
00831         h_ele_dPhiClVsPhi_propOut -> Fill(bestGsfElectron.phi(),bestGsfElectron.deltaPhiSeedClusterTrackAtCalo()); 
00832         h_ele_dPhiClVsPt_propOut -> Fill(bestGsfElectron.pt(),bestGsfElectron.deltaPhiSeedClusterTrackAtCalo()); 
00833         h_ele_HoE -> Fill(bestGsfElectron.hadronicOverEm());
00834         h_ele_HoEVsEta -> Fill( bestGsfElectron.eta(),bestGsfElectron.hadronicOverEm());
00835         h_ele_HoEVsPhi -> Fill(bestGsfElectron.phi(),bestGsfElectron.hadronicOverEm());
00836         h_ele_HoEVsE -> Fill(bestGsfElectron.caloEnergy(),bestGsfElectron.hadronicOverEm());
00837          
00838         //classes
00839         int eleClass = bestGsfElectron.classification();
00840         h_ele_classes ->Fill(eleClass); 
00841 /*
00842         if (bestGsfElectron.classification() == 0)  histSclEoEtrueGolden_barrel->Fill(sclRef->energy()/pAssSim.t());
00843         if (bestGsfElectron.classification() == 100)  histSclEoEtrueGolden_endcaps->Fill(sclRef->energy()/pAssSim.t());
00844         if (bestGsfElectron.classification() == 30)  histSclEoEtrueShowering0_barrel->Fill(sclRef->energy()/pAssSim.t());
00845         if (bestGsfElectron.classification() == 130)  histSclEoEtrueShowering0_endcaps->Fill(sclRef->energy()/pAssSim.t());
00846         if (bestGsfElectron.classification() == 31 || bestGsfElectron.classification() == 32  || bestGsfElectron.classification() == 33 || eleClass == 34 )  histSclEoEtrueShowering1234_barrel->Fill(sclRef->energy()/pAssSim.t());
00847         if (bestGsfElectron.classification() == 131 || bestGsfElectron.classification() == 132 || bestGsfElectron.classification() == 133 || eleClass == 134 )  histSclEoEtrueShowering1234_endcaps->Fill(sclRef->energy()/pAssSim.t());
00848 */
00849         eleClass = eleClass%100; // get rid of barrel/endcap distinction
00850         h_ele_eta->Fill(fabs(bestGsfElectron.eta()));
00851         if (eleClass == 0) h_ele_eta_golden ->Fill(fabs(bestGsfElectron.eta()));
00852         if (eleClass == 10) h_ele_eta_bbrem ->Fill(fabs(bestGsfElectron.eta()));
00853         if (eleClass == 20) h_ele_eta_narrow ->Fill(fabs(bestGsfElectron.eta()));
00854         if (eleClass == 30 || eleClass == 31 || eleClass == 32  || eleClass == 33 || eleClass == 34 ) h_ele_eta_shower ->Fill(fabs(bestGsfElectron.eta()));
00855 
00856         //fbrem 
00857         double fbrem_mean =  bestGsfElectron.gsfTrack()->outerMomentum().R()/bestGsfElectron.gsfTrack()->innerMomentum().R();
00858         double fbrem_mode =  bestGsfElectron.trackMomentumOut().R()/bestGsfElectron.trackMomentumAtVtx().R();
00859         h_ele_fbremVsEta_mode->Fill(bestGsfElectron.eta(),fbrem_mode);
00860         h_ele_fbremVsEta_mean->Fill(bestGsfElectron.eta(),fbrem_mean);
00861  
00862         if (eleClass == 0) h_ele_PinVsPoutGolden_mode -> Fill(bestGsfElectron.trackMomentumOut().R(), bestGsfElectron.trackMomentumAtVtx().R());
00863         if (eleClass == 30)
00864          h_ele_PinVsPoutShowering0_mode -> Fill(bestGsfElectron.trackMomentumOut().R(), bestGsfElectron.trackMomentumAtVtx().R());
00865         if (eleClass == 31 || eleClass == 32  || eleClass == 33 || eleClass == 34 )
00866          h_ele_PinVsPoutShowering1234_mode -> Fill(bestGsfElectron.trackMomentumOut().R(), bestGsfElectron.trackMomentumAtVtx().R());
00867         if (eleClass == 0) h_ele_PinVsPoutGolden_mean -> Fill(bestGsfElectron.gsfTrack()->outerMomentum().R(), bestGsfElectron.gsfTrack()->innerMomentum().R());
00868         if (eleClass == 30)
00869          h_ele_PinVsPoutShowering0_mean ->  Fill(bestGsfElectron.gsfTrack()->outerMomentum().R(), bestGsfElectron.gsfTrack()->innerMomentum().R());
00870         if (eleClass == 31 || eleClass == 32  || eleClass == 33 || eleClass == 34 )
00871          h_ele_PinVsPoutShowering1234_mean ->  Fill(bestGsfElectron.gsfTrack()->outerMomentum().R(), bestGsfElectron.gsfTrack()->innerMomentum().R());
00872         if (eleClass == 0) h_ele_PtinVsPtoutGolden_mode -> Fill(bestGsfElectron.trackMomentumOut().Rho(), bestGsfElectron.trackMomentumAtVtx().Rho());
00873         if (eleClass == 30 )
00874          h_ele_PtinVsPtoutShowering0_mode -> Fill(bestGsfElectron.trackMomentumOut().Rho(), bestGsfElectron.trackMomentumAtVtx().Rho());
00875         if (eleClass == 31 || eleClass == 32  || eleClass == 33 || eleClass == 34 )
00876          h_ele_PtinVsPtoutShowering1234_mode -> Fill(bestGsfElectron.trackMomentumOut().Rho(), bestGsfElectron.trackMomentumAtVtx().Rho());
00877         if (eleClass == 0) h_ele_PtinVsPtoutGolden_mean -> Fill(bestGsfElectron.gsfTrack()->outerMomentum().Rho(), bestGsfElectron.gsfTrack()->innerMomentum().Rho());
00878         if (eleClass == 30 )
00879          h_ele_PtinVsPtoutShowering0_mean ->  Fill(bestGsfElectron.gsfTrack()->outerMomentum().Rho(), bestGsfElectron.gsfTrack()->innerMomentum().Rho());
00880         if (eleClass == 31 || eleClass == 32  || eleClass == 33 || eleClass == 34 )
00881          h_ele_PtinVsPtoutShowering1234_mean ->  Fill(bestGsfElectron.gsfTrack()->outerMomentum().Rho(), bestGsfElectron.gsfTrack()->innerMomentum().Rho());
00882 
00883       } // gsf electron found
00884 
00885     } // mc particle found
00886 
00887     }
00888 
00889   } // loop over mc particle 
00890   
00891   h_mcNum->Fill(mcNum);
00892   h_eleNum->Fill(eleNum);
00893 
00894 }
00895 
00896