ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/CmsHi/JetAnalysis/src/HiJetResponseAnalyzer.cc

Comparing UserCode/CmsHi/JetAnalysis/src/HiJetResponseAnalyzer.cc (file contents):
Revision 1.5 by yilmaz, Wed Sep 29 15:59:10 2010 UTC vs.
Revision 1.11 by yilmaz, Fri Apr 1 12:09:17 2011 UTC

#	Line 40 | Line 40
40		#include "DataFormats/JetReco/interface/CaloJetCollection.h"
41		#include "DataFormats/JetReco/interface/GenJetCollection.h"
42
43	+	#include "CondFormats/JetMETObjects/interface/FactorizedJetCorrector.h"
44	+	#include "JetMETCorrections/Objects/interface/JetCorrectionsRecord.h"
45	+	#include "CondFormats/JetMETObjects/interface/JetCorrectionUncertainty.h"
46	+
47		#include "CommonTools/UtilAlgos/interface/TFileService.h"
48		#include "FWCore/ServiceRegistry/interface/Service.h"
49	+	#include "DataFormats/Math/interface/deltaR.h"
50	+	#include "FWCore/ParameterSet/interface/FileInPath.h"
51	+	#include "FWCore/Framework/interface/ESWatcher.h"
52
53	<
53	>	#include "CmsHi/JetAnalysis/interface/RhoGetter.h"
54		#include "TTree.h"
55
56		using namespace std;
57	+	using namespace edm;
58
59		static const int MAXJETS = 500;
60
#	Line 55 | Line 63 | struct etdr{
63		double dr;
64		};
65
66	<	struct JRA{
66	>	class JRA{
67
68	+	public:
69		int nref;
70		int bin;
71		float b;
72		float hf;
73		float jtpt[MAXJETS];
74	<	float jtcorpt[MAXJETS];
74	>	float jtrawpt[MAXJETS];
75		float refpt[MAXJETS];
76		float jteta[MAXJETS];
77		float refeta[MAXJETS];
78		float jtphi[MAXJETS];
79		float refphi[MAXJETS];
80	+	float l2[MAXJETS];
81	+	float l3[MAXJETS];
82	+	float area[MAXJETS];
83	+	float pu[MAXJETS];
84	+	float rho[MAXJETS];
85
86		float weight;
87		};
88
89	+	struct JRAV{
90	+	int index;
91	+	float jtpt;
92	+	float jtrawpt;
93	+	float refpt;
94	+	float refcorpt;
95	+	float jteta;
96	+	float refeta;
97	+	float jtphi;
98	+	float refphi;
99	+	float l2;
100	+	float l3;
101	+	float area;
102	+	float pu;
103	+	float rho;
104	+
105	+	};
106	+
107		//
108		// class declaration
109		//
110	+	bool comparePt(JRAV a, JRAV b) {return a.jtpt > b.jtpt;}
111
112		class HiJetResponseAnalyzer : public edm::EDAnalyzer {
113		public:
#	Line 87 | Line 120 | class HiJetResponseAnalyzer : public edm
120		virtual void analyze(const edm::Event&, const edm::EventSetup&);
121		virtual void endJob() ;
122		bool selectJet(int i);
90	–
123		// ----------member data ---------------------------
124
125	<	bool usePat_;
126	<	bool doMC_;
127	<	bool filterJets_;
128	<	bool diJetsOnly_;
129	<	bool matchDiJets_;
125	>	bool usePat_;
126	>	bool doMC_;
127	>	bool filterJets_;
128	>	bool diJetsOnly_;
129	>	bool matchDiJets_;
130	>	bool matchPatGen_;
131	>	bool matchNew_;
132	>	bool sortJets_;
133	>	bool correctJets_;
134	>	bool getFastJets_;
135
136	+	double matchR_;
137		double genPtMin_;
138		double ptMin_;
139		double emfMin_;
140		double n90Min_;
141		double n90hitMin_;
142
143	<	edm::InputTag jetTag_;
143	>	edm::InputTag jetTag_;
144	>	edm::InputTag matchTag_;
145	>	std::vector<edm::InputTag> matchTags_;
146	>
147	>	JRA jra_;
148	>	std::vector<JRA> jraMatch_;
149
107	–	JRA jra_;
150		TTree* t;
151
152		edm::Handle<edm::GenHIEvent> mc;
#	Line 113 | Line 155 | class HiJetResponseAnalyzer : public edm
155		edm::Handle<reco::JetView> jets;
156		edm::Handle<pat::JetCollection> patjets;
157
158	+	FactorizedJetCorrector* jetCorrector_;
159	+	edm::ESWatcher<JetCorrectionsRecord> watchJetCorrectionsRecord_;
160	+
161	+	std::string tags_;
162	+	std::string levels_;
163	+	std::string algo_;
164	+
165		edm::Service<TFileService> fs;
166
167	+	edm::Handle<vector<double> > ktRhos;
168	+	edm::Handle<vector<double> > akRhos;
169	+	bool doFastJets_;
170	+
171	+	vector<int> doMatchedFastJets_;
172	+	vector<int> correctMatchedJets_;
173	+
174	+	InputTag ktSrc_;
175	+	InputTag akSrc_;
176	+
177	+
178		};
179
180		bool HiJetResponseAnalyzer::selectJet(int i){
121	–
181		const reco::Jet& jet = (*jets)[i];
123	–
182		if(usePat_){
125	–	//      cout<<"a"<<endl;
183		const pat::Jet& patjet = (*patjets)[i];
127	–	//      cout<<"b"<<endl;
128	–
184		if(patjet.emEnergyFraction() <= emfMin_) return false;
130	–	//      cout<<"c"<<endl;
131	–
185		if(patjet.jetID().n90Hits <= n90hitMin_) return false;
133	–	//      cout<<"d"<<endl;
134	–
186		if(doMC_){
187
188		}
#	Line 139 | Line 190 | bool HiJetResponseAnalyzer::selectJet(in
190		}
191
192		return true;
142	–
193		}
194
195
146	–
196		//
197		// constants, enums and typedefs
198		//
#	Line 158 | Line 207 | bool HiJetResponseAnalyzer::selectJet(in
207		HiJetResponseAnalyzer::HiJetResponseAnalyzer(const edm::ParameterSet& iConfig)
208
209		{
210	+
211	+	levels_ = iConfig.getUntrackedParameter<string>("corrLevels","L2Relative:L3Absolute");
212	+
213	+	algo_ = iConfig.getUntrackedParameter<string>("algo","IC5Calo");
214	+	tags_ = "";
215	+
216	+	string l[2] = {"L2Relative","L3Absolute"};
217	+
218	+	for(int i = 0; i <2; ++i){
219	+	edm::FileInPath fip("CondFormats/JetMETObjects/data/Spring10_"+l[i]+"_"+algo_+".txt");
220	+	tags_ += fip.fullPath();
221	+	if(i < 2 - 1)tags_ +=":";
222	+	}
223	+
224	+	jetCorrector_ = new FactorizedJetCorrector(levels_, tags_);
225	+
226		//now do what ever initialization is needed
227	+	matchR_ = iConfig.getUntrackedParameter<double>("matchR",0.25);
228
229		ptMin_ = iConfig.getUntrackedParameter<double>("ptMin",0);
230	<	genPtMin_ = iConfig.getUntrackedParameter<double>("genPtMin",0);
231	<	emfMin_ = iConfig.getUntrackedParameter<double>("emfMin",0);
232	<	n90Min_ = iConfig.getUntrackedParameter<double>("n90Min",0);
233	<	n90hitMin_ = iConfig.getUntrackedParameter<double>("n90hitMin",0);
230	>	genPtMin_ = iConfig.getUntrackedParameter<double>("genPtMin",20);
231	>	emfMin_ = iConfig.getUntrackedParameter<double>("emfMin",0.01);
232	>	n90Min_ = iConfig.getUntrackedParameter<double>("n90Min",1);
233	>	n90hitMin_ = iConfig.getUntrackedParameter<double>("n90hitMin",1);
234
235		filterJets_ = iConfig.getUntrackedParameter<bool>("filterJets",true);
236	<	diJetsOnly_ = iConfig.getUntrackedParameter<bool>("diJetsOnly",true);
237	<	matchDiJets_ = iConfig.getUntrackedParameter<bool>("usePat",true);
236	>	diJetsOnly_ = iConfig.getUntrackedParameter<bool>("diJetsOnly",false);
237	>	matchDiJets_ = iConfig.getUntrackedParameter<bool>("matchDiJets",false);
238	>	matchPatGen_ = iConfig.getUntrackedParameter<bool>("matchPatGen",false);
239	>
240	>	matchNew_ = iConfig.getUntrackedParameter<bool>("matchNew",false);
241	>
242		usePat_ = iConfig.getUntrackedParameter<bool>("usePat",true);
243		doMC_ = iConfig.getUntrackedParameter<bool>("doMC",true);
244	+
245	+	sortJets_ = iConfig.getUntrackedParameter<bool>("sortJets",true);
246	+	correctJets_ = iConfig.getUntrackedParameter<bool>("correctJets",false);
247	+
248	+	correctMatchedJets_ = iConfig.getUntrackedParameter<std::vector<int> >("correctMatchedJets",std::vector<int>(0));
249	+	doMatchedFastJets_ = iConfig.getUntrackedParameter<std::vector<int> >("doMatchedFastJets",std::vector<int>(0));
250	+
251		jetTag_ = iConfig.getUntrackedParameter<edm::InputTag>("src",edm::InputTag("selectedPatJets"));
252	+	matchTag_ = iConfig.getUntrackedParameter<edm::InputTag>("match",edm::InputTag("selectedPatJets"));
253	+	matchTags_ = iConfig.getUntrackedParameter<std::vector<edm::InputTag> >("matches",std::vector<edm::InputTag>(0));
254	+
255	+	ktSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("ktSrc",edm::InputTag("kt4CaloJets"));
256	+	akSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("akSrc",edm::InputTag("ak5CaloJets"));
257	+	doFastJets_ = iConfig.getUntrackedParameter<bool>("doFastJets",false);
258	+
259	+	getFastJets_ = iConfig.getUntrackedParameter<bool>("getFastJets",false);
260	+
261	+	for(unsigned int i = 0; i < doMatchedFastJets_.size(); ++i){
262	+	getFastJets_ = getFastJets_ || (bool)doMatchedFastJets_[i];
263	+	}
264	+
265
266		}
267
#	Line 197 | Line 287 | HiJetResponseAnalyzer::analyze(const edm
287
288		iEvent.getByLabel(jetTag_,jets);
289		if(usePat_)iEvent.getByLabel(jetTag_,patjets);
290	+	std::vector<JRAV> jraV;
291
292	<	jra_.nref = 0;
293	<	for(unsigned int j = 0 ; j < jets->size(); ++j){
294	<
295	<	//      const pat::Jet& jet = (*jets)[j];
205	<	if(filterJets_ && !selectJet(j)) continue;
206	<
207	<	const reco::Jet& jet = (*jets)[j];
208	<
209	<	jra_.jtpt[jra_.nref] = jet.pt();
210	<	jra_.jteta[jra_.nref] = jet.eta();
211	<	jra_.jtphi[jra_.nref] = jet.phi();
212	<	jra_.jtcorpt[jra_.nref] = jet.pt();
292	>	if(getFastJets_){
293	>	iEvent.getByLabel(edm::InputTag(ktSrc_.label(),"rhos"),ktRhos);
294	>	iEvent.getByLabel(edm::InputTag(akSrc_.label(),"rhos"),akRhos);
295	>	}
296
297	<	if(usePat_){
298	<	const pat::Jet& patjet = (*patjets)[j];
297	>	for(unsigned int j = 0 ; j < jets->size(); ++j){
298	>	if(filterJets_ && !selectJet(j)) continue;
299	>	const reco::Jet& jet = (*jets)[j];
300	>	JRAV jv;
301	>	jv.jtpt = jet.pt();
302	>	jv.jteta = jet.eta();
303	>	jv.jtphi = jet.phi();
304	>	jv.jtrawpt = jet.pt();
305	>	jv.area = jet.jetArea();
306	>	jv.pu  = jet.pileup();
307	>	jv.index = j;
308	>
309	>	double pt = jet.pt();
310	>	double ktRho = -1, akRho = -1;
311	>	if(getFastJets_){
312	>	ktRho = getRho(jv.jteta,*ktRhos);
313	>	akRho = getRho(jv.jteta,*akRhos);
314	>	}
315	>
316	>	jv.rho = akRho;
317	>
318	>	if(doFastJets_){
319	>	jv.pu = jet.jetArea()*akRho;
320	>	pt -= jv.pu;
321	>	}
322	>	jv.jtpt = pt;
323	>
324	>	if(correctJets_){
325	>	jetCorrector_->setJetEta(jet.eta());
326	>	jetCorrector_->setJetPt(pt);
327	>	//      jetCorrector_->setJetE(jet.energy());
328	>	vector<float> corrs = jetCorrector_->getSubCorrections();
329	>	jv.l2 = corrs[0];
330	>	jv.l3 = corrs[1];
331	>	jv.jtpt = pt*jv.l2*jv.l3;
332	>	}
333	>
334	>
335	>	if(usePat_){
336	>	const pat::Jet& patjet = (*patjets)[j];
337	>
338	>	jv.jtrawpt = patjet.correctedJet("raw").pt();
339	>	jv.jtpt = patjet.pt();
340	>
341	>	if(doMC_ && matchPatGen_ && patjet.genJet() != 0){
342	>	if(patjet.genJet()->pt() < genPtMin_) continue;
343	>	jv.refpt = patjet.genJet()->pt();
344	>	jv.refeta = patjet.genJet()->eta();
345	>	jv.refphi = patjet.genJet()->phi();
346	>	}else{
347	>	jv.refpt = -99;
348	>	jv.refeta = -99;
349	>	jv.refphi = -99;
350	>	}
351	>	}
352	>	jraV.push_back(jv);
353	>	}
354
355	<	jra_.jtpt[jra_.nref] = patjet.correctedJet("raw").pt();
356	<	jra_.jtcorpt[jra_.nref] = patjet.pt();
219	<
220	<	if(doMC_ && patjet.genJet() != 0){
221	<	if(patjet.genJet()->pt() < genPtMin_) continue;
222	<	jra_.refpt[jra_.nref] = patjet.genJet()->pt();
223	<	jra_.refeta[jra_.nref] = patjet.genJet()->eta();
224	<	jra_.refphi[jra_.nref] = patjet.genJet()->phi();
225	<
226	<	}else{
227	<	jra_.refpt[jra_.nref] = -99;
228	<	jra_.refeta[jra_.nref] = -99;
229	<	jra_.refphi[jra_.nref] = -99;
230	<	}
231	<	}
232	<
233	<	jra_.nref++;
234	<
355	>	if(sortJets_){
356	>	std::sort(jraV.begin(),jraV.end(),comparePt);
357		}
358
359	+	for(unsigned int i = 0; i < jraV.size(); ++i){
360	+	JRAV& jv = jraV[i];
361	+	const reco::Jet& jet = (*jets)[jv.index];
362	+
363	+	if(matchNew_){
364	+	for(unsigned int im = 0; im < matchTags_.size(); ++im){
365	+	edm::Handle<reco::JetView> matchedJets;
366	+	iEvent.getByLabel(matchTags_[im],matchedJets);
367	+	jraMatch_[im].jtrawpt[i] = -99;
368	+	jraMatch_[im].jtpt[i] = -99;
369	+	jraMatch_[im].jteta[i] = -99;
370	+	jraMatch_[im].jtphi[i] = -99;
371	+	jraMatch_[im].area[i] = -99;
372	+	jraMatch_[im].l2[i] = -99;
373	+	jraMatch_[im].l3[i] = -99;
374	+	jraMatch_[im].pu[i] = -99;
375	+	for(unsigned int m = 0 ; m < matchedJets->size(); ++m){
376	+	const reco::Jet& match = (*matchedJets)[m];
377	+	double dr = reco::deltaR(jet.eta(),jet.phi(),match.eta(),match.phi());
378	+	if(dr < matchR_ && match.pt() > genPtMin_){
379	+	jraMatch_[im].jtpt[i] = match.pt();
380	+	jraMatch_[im].jtrawpt[i] = match.pt();
381	+	jraMatch_[im].jteta[i] = match.eta();
382	+	jraMatch_[im].jtphi[i] = match.phi();
383	+	jraMatch_[im].area[i] = match.jetArea();
384	+
385	+	double ktRhoM = -1, akRhoM = -1;
386	+	if(getFastJets_){
387	+	ktRhoM = getRho(jraMatch_[im].jteta[i],*ktRhos);
388	+	akRhoM = getRho(jraMatch_[im].jteta[i],*akRhos);
389	+	}
390	+
391	+	jraMatch_[im].rho[i] = akRhoM;
392	+	double pt = jraMatch_[im].jtpt[i];
393	+
394	+	if((bool)doMatchedFastJets_[im]){
395	+	jraMatch_[im].pu[i] = jraMatch_[im].area[i]*akRhoM;
396	+	pt -= jraMatch_[im].pu[i];
397	+	}
398	+
399	+	jraMatch_[im].jtpt[i] = pt;
400	+
401	+	if((bool)correctMatchedJets_[im]){
402	+
403	+	jetCorrector_->setJetEta(jraMatch_[im].jteta[i]);
404	+	jetCorrector_->setJetPt(pt);
405	+	//              jetCorrector_->setJetE(match.energy());
406	+	vector<float> corrs = jetCorrector_->getSubCorrections();
407	+	jraMatch_[im].l2[i] = corrs[0];
408	+	// Should it be re-set for L3???
409	+	jraMatch_[im].l3[i] = corrs[1];
410	+	jraMatch_[im].jtpt[i] = pt*jraMatch_[im].l2[i]*jraMatch_[im].l3[i];
411	+	}
412	+
413	+	}
414	+	}
415	+	}
416	+	}
417	+
418	+	jra_.jtpt[i] = jv.jtpt;
419	+	jra_.jteta[i] = jv.jteta;
420	+	jra_.jtphi[i] = jv.jtphi;
421	+	jra_.jtrawpt[i] = jv.jtrawpt;
422	+	jra_.refpt[i] = jv.refpt;
423	+	jra_.refeta[i] = jv.refeta;
424	+	jra_.refphi[i] = jv.refphi;
425	+
426	+	jra_.area[i] = jv.area;
427	+	jra_.pu[i] = jv.pu;
428	+	jra_.rho[i] = jv.rho;
429	+	jra_.l2[i] = jv.l2;
430	+	jra_.l3[i] = jv.l3;
431	+
432	+	}
433	+	jra_.nref = jraV.size();
434	+
435		t->Fill();
436
437		}
438
241	–
439		// ------------ method called once each job just before starting event loop  ------------
440	+
441	+
442		void
443	<	HiJetResponseAnalyzer::beginJob()
444	<	{
443	>	HiJetResponseAnalyzer::beginJob(){
444	>	t= fs->make<TTree>("t","Jet Response Analyzer");
445	>	t->Branch("nref",&jra_.nref,"nref/I");
446	>	t->Branch("jtpt",jra_.jtpt,"jtpt[nref]/F");
447	>	t->Branch("jteta",jra_.jteta,"jteta[nref]/F");
448	>	t->Branch("jtphi",jra_.jtphi,"jtphi[nref]/F");
449	>	t->Branch("jtrawpt",jra_.jtrawpt,"jtrawpt[nref]/F");
450	>
451	>	if(correctJets_){
452	>	t->Branch("l2",jra_.l2,"l2[nref]/F");
453	>	t->Branch("l3",jra_.l3,"l3[nref]/F");
454	>	}
455	>	t->Branch("area",jra_.area,"area[nref]/F");
456	>	t->Branch("pu",jra_.pu,"pu[nref]/F");
457	>	t->Branch("rho",jra_.rho,"rho[nref]/F");
458
459	<	t= fs->make<TTree>("t","Jet Response Analyzer");
460	<	t->Branch("b",&jra_.b,"b/F");
249	<	t->Branch("hf",&jra_.hf,"hf/F");
250	<	t->Branch("nref",&jra_.nref,"nref/I");
251	<	t->Branch("jtpt",jra_.jtpt,"jtpt[nref]/F");
252	<	t->Branch("jtcorpt",jra_.jtcorpt,"jtcorpt[nref]/F");
253	<	t->Branch("refpt",jra_.refpt,"refpt[nref]/F");
254	<	t->Branch("jteta",jra_.jteta,"jteta[nref]/F");
459	>	t->Branch("refpt",jra_.refpt,"refpt[nref]/F");
460	>	t->Branch("refcorpt",jra_.refpt,"refcorpt[nref]/F");
461		t->Branch("refeta",jra_.refeta,"refeta[nref]/F");
256	–	t->Branch("jtphi",jra_.jtphi,"jtphi[nref]/F");
462		t->Branch("refphi",jra_.refphi,"refphi[nref]/F");
463		t->Branch("weight",&jra_.weight,"weight/F");
259	–	t->Branch("bin",&jra_.bin,"bin/I");
464
465	+	jraMatch_.clear();
466	+	for(unsigned int im = 0; im < matchTags_.size(); ++im){
467	+	JRA jrm;
468	+	jraMatch_.push_back(jrm);
469	+	}
470
471	+	for(unsigned int im = 0; im < matchTags_.size(); ++im){
472	+	t->Branch(Form("jtpt%d",im),jraMatch_[im].jtpt,Form("jtpt%d[nref]/F",im));
473	+	t->Branch(Form("jteta%d",im),jraMatch_[im].jteta,Form("jteta%d[nref]/F",im));
474	+	t->Branch(Form("jtphi%d",im),jraMatch_[im].jtphi,Form("jtphi%d[nref]/F",im));
475	+	t->Branch(Form("jtrawpt%d",im),jraMatch_[im].jtrawpt,Form("jtrawpt%d[nref]/F",im));
476	+
477	+	if((bool)correctMatchedJets_[im]){
478	+	t->Branch(Form("l2_%d",im),jraMatch_[im].l2,Form("l2_%d[nref]/F",im));
479	+	t->Branch(Form("l3_%d",im),jraMatch_[im].l3,Form("l3_%d[nref]/F",im));
480	+	}
481	+
482	+	t->Branch(Form("area%d",im),jraMatch_[im].area,Form("area%d[nref]/F",im));
483	+	t->Branch(Form("pu%d",im),jraMatch_[im].pu,Form("pu%d[nref]/F",im));
484	+	t->Branch(Form("rho%d",im),jraMatch_[im].rho,Form("rho%d[nref]/F",im));
485
486	+	}
487	+
488		}
489
490		// ------------ method called once each job just after ending the event loop  ------------

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines