00001
00002
00003
00004
00005
00006 #include "RecoJets/JetAnalyzers/interface/JetValidation.h"
00007 #include "RecoJets/JetAlgorithms/interface/JetAlgoHelper.h"
00008 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
00009 #include "DataFormats/JetReco/interface/GenJetCollection.h"
00010 #include "DataFormats/JetReco/interface/CaloJet.h"
00011 #include "DataFormats/JetReco/interface/GenJet.h"
00012 #include "DataFormats/JetReco/interface/JetTracksAssociation.h"
00013 #include "DataFormats/Math/interface/deltaR.h"
00014 #include "FWCore/Framework/interface/Event.h"
00015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00016 #include <TROOT.h>
00017 #include <TSystem.h>
00018 #include <TFile.h>
00019 #include <cmath>
00020 using namespace edm;
00021 using namespace reco;
00022 using namespace std;
00024 JetValidation::JetValidation(edm::ParameterSet const& cfg)
00025 {
00026 dRmatch = cfg.getParameter<double> ("dRmatch");
00027 PtMin = cfg.getParameter<double> ("PtMin");
00028 Njets = cfg.getParameter<int> ("Njets");
00029 MCarlo = cfg.getParameter<bool> ("MCarlo");
00030 genAlgo = cfg.getParameter<string> ("genAlgo");
00031 calAlgo = cfg.getParameter<string> ("calAlgo");
00032 jetTracksAssociator = cfg.getParameter<string> ("jetTracksAssociator");
00033 histoFileName = cfg.getParameter<string> ("histoFileName");
00034 }
00036 void JetValidation::beginJob()
00037 {
00038 m_file = new TFile(histoFileName.c_str(),"RECREATE");
00039
00040 m_HistNames1D["CaloJetMulti"] = new TH1F("CaloJetMulti","Multiplicity of CaloJets",100,0,100);
00041 m_HistNames1D["ptCalo"] = new TH1F("ptCalo","p_{T} of CaloJets",7000,0,7000);
00042 m_HistNames1D["etaCalo"] = new TH1F("etaCalo","#eta of CaloJets",100,-5.0,5.0);
00043 m_HistNames1D["phiCalo"] = new TH1F("phiCalo","#phi of CaloJets",72,-M_PI, M_PI);
00044 m_HistNames1D["m2jCalo"] = new TH1F("m2jCalo","Dijet Mass of leading CaloJets",7000,0,14000);
00045 m_HistNames1D["nTracks"] = new TH1F("nTracks","Number of tracks associated with a jet",100,0,100);
00046 m_HistNames1D["chargeFraction"] = new TH1F("chargeFraction","Fraction of charged tracks pt",500,0,5);
00047 m_HistNames1D["emEnergyFraction"] = new TH1F("emEnergyFraction","Jets EM Fraction",110,0,1.1);
00048 m_HistNames1D["emEnergyInEB"] = new TH1F("emEnergyInEB","Jets emEnergyInEB",7000,0,14000);
00049 m_HistNames1D["emEnergyInEE"] = new TH1F("emEnergyInEE","Jets emEnergyInEE",7000,0,14000);
00050 m_HistNames1D["emEnergyInHF"] = new TH1F("emEnergyInHF","Jets emEnergyInHF",7000,0,14000);
00051 m_HistNames1D["hadEnergyInHB"] = new TH1F("hadEnergyInHB","Jets hadEnergyInHB",7000,0,14000);
00052 m_HistNames1D["hadEnergyInHE"] = new TH1F("hadEnergyInHE","Jets hadEnergyInHE",7000,0,14000);
00053 m_HistNames1D["hadEnergyInHF"] = new TH1F("hadEnergyInHF","Jets hadEnergyInHF",7000,0,14000);
00054 m_HistNames1D["hadEnergyInHO"] = new TH1F("hadEnergyInHO","Jets hadEnergyInHO",7000,0,14000);
00055 m_HistNamesProfile["EBfractionVsEta"] = new TProfile("EBfractionVsEta","Jets EBfraction vs #eta",100,-5.0,5.0);
00056 m_HistNamesProfile["EEfractionVsEta"] = new TProfile("EEfractionVsEta","Jets EEfraction vs #eta",100,-5.0,5.0);
00057 m_HistNamesProfile["HBfractionVsEta"] = new TProfile("HBfractionVsEta","Jets HBfraction vs #eta",100,-5.0,5.0);
00058 m_HistNamesProfile["HOfractionVsEta"] = new TProfile("HOfractionVsEta","Jets HOfraction vs #eta",100,-5.0,5.0);
00059 m_HistNamesProfile["HEfractionVsEta"] = new TProfile("HEfractionVsEta","Jets HEfraction vs #eta",100,-5.0,5.0);
00060 m_HistNamesProfile["HFfractionVsEta"] = new TProfile("HFfractionVsEta","Jets HFfraction vs #eta",100,-5.0,5.0);
00061 m_HistNamesProfile["CaloEnergyVsEta"] = new TProfile("CaloEnergyVsEta","CaloJets Energy Vs. Eta",100,-5.0,5.0);
00062 m_HistNamesProfile["emEnergyVsEta"] = new TProfile("emEnergyVsEta","Jets EM Energy Vs. Eta",100,-5.0,5.0);
00063 m_HistNamesProfile["hadEnergyVsEta"] = new TProfile("hadEnergyVsEta","Jets HAD Energy Vs. Eta",100,-5.0,5.0);
00064 if (MCarlo)
00065 {
00066 m_HistNames1D["GenJetMulti"] = new TH1F("GenJetMulti","Multiplicity of GenJets",100,0,100);
00067 m_HistNames1D["ptHat"] = new TH1F("ptHat","p_{T}hat",7000,0,7000);
00068 m_HistNames1D["ptGen"] = new TH1F("ptGen","p_{T} of GenJets",7000,0,7000);
00069 m_HistNames1D["etaGen"] = new TH1F("etaGen","#eta of GenJets",100,-5.0,5.0);
00070 m_HistNames1D["phiGen"] = new TH1F("phiGen","#phi of GenJets",72,-M_PI, M_PI);
00071 m_HistNames1D["m2jGen"] = new TH1F("m2jGen","Dijet Mass of leading GenJets",7000,0,14000);
00072 m_HistNames1D["dR"] = new TH1F("dR","GenJets dR with matched CaloJet",200,0,1);
00073 m_HistNamesProfile["GenEnergyVsEta"] = new TProfile("GenEnergyVsEta","GenJets Energy Vs. Eta",100,-5.0,5.0);
00074 m_HistNamesProfile["respVsPtBarrel"] = new TProfile("respVsPtBarrel","CaloJet Response of GenJets in Barrel",7000,0,7000);
00075 m_HistNamesProfile["CaloErespVsEta"] = new TProfile("CaloErespVsEta","Jets Energy Response Vs. Eta",100,-5.0,5.0);
00076 m_HistNamesProfile["emErespVsEta"] = new TProfile("emErespVsEta","Jets EM Energy Response Vs. Eta",100,-5.0,5.0);
00077 m_HistNamesProfile["hadErespVsEta"] = new TProfile("hadErespVsEta","Jets HAD Energy Response Vs. Eta",100,-5.0,5.0);
00078 }
00079 }
00081 void JetValidation::analyze(edm::Event const& evt, edm::EventSetup const& iSetup)
00082 {
00083 math::XYZTLorentzVector p4jet[2];
00084 int jetInd,jetCounter,nTracks;
00085 double dRmin,dR,e,eta,emEB,emEE,emHF,hadHB,hadHE,hadHO,hadHF,pt,phi,pthat,chf;
00086 Handle<CaloJetCollection> caljets;
00087 Handle<GenJetCollection> genjets;
00088 Handle<double> genEventScale;
00089 Handle<JetTracksAssociation::Container> jetTracks;
00090 CaloJetCollection::const_iterator i_caljet;
00091 GenJetCollection::const_iterator i_genjet;
00092 evt.getByLabel(calAlgo,caljets);
00093 evt.getByLabel(jetTracksAssociator,jetTracks);
00094 jetInd = 0;
00095 jetCounter = 0;
00096 if (caljets->size()==0)
00097 cout<<"WARNING: NO calo jets in event "<<evt.id().event()<<", Run "<<evt.id().run()<<" !!!!"<<endl;
00098 for(i_caljet = caljets->begin(); i_caljet != caljets->end() && jetInd<Njets; ++i_caljet)
00099 {
00100 e = i_caljet->energy();
00101 pt = i_caljet->pt();
00102 phi = i_caljet->phi();
00103 eta = i_caljet->eta();
00104 emEB = i_caljet->emEnergyInEB();
00105 emEE = i_caljet->emEnergyInEE();
00106 emHF = i_caljet->emEnergyInHF();
00107 hadHB = i_caljet->hadEnergyInHB();
00108 hadHE = i_caljet->hadEnergyInHE();
00109 hadHO = i_caljet->hadEnergyInHO();
00110 hadHF = i_caljet->hadEnergyInHF();
00111 nTracks = JetTracksAssociation::tracksNumber(*jetTracks,*i_caljet);
00112 chf = (JetTracksAssociation::tracksP4(*jetTracks,*i_caljet)).pt()/pt;
00113 if (jetInd<2)
00114 p4jet[jetInd] = i_caljet->p4();
00115 if (pt>PtMin)
00116 {
00117 FillHist1D("ptCalo",pt);
00118 FillHist1D("etaCalo",eta);
00119 FillHist1D("phiCalo",phi);
00120 FillHist1D("emEnergyFraction",i_caljet->emEnergyFraction());
00121 FillHist1D("nTracks",nTracks);
00122 FillHist1D("chargeFraction",chf);
00123 FillHist1D("emEnergyInEB",emEB);
00124 FillHist1D("emEnergyInEE",emEE);
00125 FillHist1D("emEnergyInHF",emHF);
00126 FillHist1D("hadEnergyInHB",hadHB);
00127 FillHist1D("hadEnergyInHE",hadHE);
00128 FillHist1D("hadEnergyInHF",hadHF);
00129 FillHist1D("hadEnergyInHO",hadHO);
00130 FillHistProfile("EBfractionVsEta",eta,emEB/e);
00131 FillHistProfile("EEfractionVsEta",eta,emEE/e);
00132 FillHistProfile("HBfractionVsEta",eta,hadHB/e);
00133 FillHistProfile("HOfractionVsEta",eta,hadHO/e);
00134 FillHistProfile("HEfractionVsEta",eta,hadHE/e);
00135 FillHistProfile("HFfractionVsEta",eta,(hadHF+emHF)/e);
00136 FillHistProfile("CaloEnergyVsEta",eta,e);
00137 FillHistProfile("emEnergyVsEta",eta,emEB+emEE+emHF);
00138 FillHistProfile("hadEnergyVsEta",eta,hadHB+hadHO+hadHE+hadHF);
00139 jetCounter++;
00140 }
00141 jetInd++;
00142 }
00143 FillHist1D("CaloJetMulti",jetCounter);
00144 if (jetInd>1)
00145 FillHist1D("m2jCalo",(p4jet[0]+p4jet[1]).mass());
00146 if (MCarlo)
00147 {
00148 evt.getByLabel(genAlgo,genjets);
00149 evt.getByLabel("genEventScale",genEventScale);
00150 pthat = *genEventScale;
00151 FillHist1D("ptHat",pthat);
00152 CaloJet MatchedJet;
00153 jetInd = 0;
00154 if (genjets->size()==0)
00155 cout<<"WARNING: NO gen jets in event "<<evt.id().event()<<", Run "<<evt.id().run()<<" !!!!"<<endl;
00156 for(i_genjet = genjets->begin(); i_genjet != genjets->end() && jetInd<Njets; ++i_genjet)
00157 {
00158 if (jetInd<2)
00159 p4jet[jetInd] = i_genjet->p4();
00160 FillHist1D("ptGen",i_genjet->pt());
00161 FillHist1D("etaGen",i_genjet->eta());
00162 FillHist1D("phiGen",i_genjet->phi());
00163 FillHistProfile("GenEnergyVsEta",i_genjet->eta(),i_genjet->energy());
00164 dRmin=1000.0;
00165 for(i_caljet = caljets->begin(); i_caljet != caljets->end(); ++i_caljet)
00166 {
00167 dR = deltaR(i_caljet->eta(),i_caljet->phi(),i_genjet->eta(),i_genjet->phi());
00168 if (dR<dRmin)
00169 {
00170 dRmin = dR;
00171 MatchedJet = *i_caljet;
00172 }
00173 }
00174 FillHist1D("dR",dRmin);
00175 e = MatchedJet.energy();
00176 pt = MatchedJet.pt();
00177 eta = MatchedJet.eta();
00178 phi = MatchedJet.phi();
00179 emEB = MatchedJet.emEnergyInEB();
00180 emEE = MatchedJet.emEnergyInEE();
00181 emHF = MatchedJet.emEnergyInHF();
00182 hadHB = MatchedJet.hadEnergyInHB();
00183 hadHE = MatchedJet.hadEnergyInHE();
00184 hadHO = MatchedJet.hadEnergyInHO();
00185 hadHF = MatchedJet.hadEnergyInHF();
00186 if (dRmin<dRmatch && pt>PtMin)
00187 {
00188 FillHistProfile("CaloErespVsEta",eta,e/i_genjet->energy());
00189 FillHistProfile("emErespVsEta",eta,(emEB+emEE+emHF)/i_genjet->energy());
00190 FillHistProfile("hadErespVsEta",eta,(hadHB+hadHO+hadHE+hadHF)/i_genjet->energy());
00191 if (fabs(i_genjet->eta())<1.)
00192 FillHistProfile("respVsPtBarrel",i_genjet->pt(),pt/i_genjet->pt());
00193 }
00194 jetInd++;
00195 }
00196 FillHist1D("GenJetMulti",jetInd);
00197 if (jetInd>1)
00198 FillHist1D("m2jGen",(p4jet[0]+p4jet[1]).mass());
00199 }
00200 }
00202 void JetValidation::endJob()
00203 {
00205 if (m_file !=0)
00206 {
00207 m_file->cd();
00208 for (std::map<TString, TH1*>::iterator hid = m_HistNames1D.begin(); hid != m_HistNames1D.end(); hid++)
00209 hid->second->Write();
00210 for (std::map<TString, TH2*>::iterator hid = m_HistNames2D.begin(); hid != m_HistNames2D.end(); hid++)
00211 hid->second->Write();
00212 for (std::map<TString, TProfile*>::iterator hid = m_HistNamesProfile.begin(); hid != m_HistNamesProfile.end(); hid++)
00213 hid->second->Write();
00214 delete m_file;
00215 m_file = 0;
00216 }
00217 }
00219 void JetValidation::FillHist1D(const TString& histName,const Double_t& value)
00220 {
00221 std::map<TString, TH1*>::iterator hid=m_HistNames1D.find(histName);
00222 if (hid==m_HistNames1D.end())
00223 std::cout << "%fillHist -- Could not find histogram with name: " << histName << std::endl;
00224 else
00225 hid->second->Fill(value);
00226 }
00228 void JetValidation::FillHist2D(const TString& histName,const Double_t& valuex,const Double_t& valuey)
00229 {
00230 std::map<TString, TH2*>::iterator hid=m_HistNames2D.find(histName);
00231 if (hid==m_HistNames2D.end())
00232 std::cout << "%fillHist -- Could not find histogram with name: " << histName << std::endl;
00233 else
00234 hid->second->Fill(valuex,valuey);
00235 }
00237 void JetValidation::FillHistProfile(const TString& histName,const Double_t& valuex,const Double_t& valuey)
00238 {
00239 std::map<TString, TProfile*>::iterator hid=m_HistNamesProfile.find(histName);
00240 if (hid==m_HistNamesProfile.end())
00241 std::cout << "%fillHist -- Could not find histogram with name: " << histName << std::endl;
00242 else
00243 hid->second->Fill(valuex,valuey);
00244 }
00245 #include "FWCore/Framework/interface/MakerMacros.h"
00246 DEFINE_FWK_MODULE(JetValidation);
1.7.3