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.6 by yilmaz, Mon Oct 18 16:13:37 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	+	#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
75	–
89		struct JRAV{
90	<
90	>	int index;
91		float jtpt;
92	<	float jtcorpt;
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		//
92	–	bool compareCorPt(JRAV a, JRAV b) {return a.jtcorpt > b.jtcorpt;}
110		bool comparePt(JRAV a, JRAV b) {return a.jtpt > b.jtpt;}
111
112		class HiJetResponseAnalyzer : public edm::EDAnalyzer {
#	Line 114 | Line 131 | class HiJetResponseAnalyzer : public edm
131		bool matchNew_;
132		bool sortJets_;
133		bool correctJets_;
134	+	bool getFastJets_;
135
136		double matchR_;
137		double genPtMin_;
#	Line 122 | Line 140 | class HiJetResponseAnalyzer : public edm
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
128	–	JRA jra_;
150		TTree* t;
151
152		edm::Handle<edm::GenHIEvent> mc;
#	Line 133 | Line 154 | class HiJetResponseAnalyzer : public edm
154
155		edm::Handle<reco::JetView> jets;
156		edm::Handle<pat::JetCollection> patjets;
157	<	edm::Handle<reco::JetView> matchedJets;
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){
#	Line 169 | 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
#	Line 191 | Line 245 | HiJetResponseAnalyzer::HiJetResponseAnal
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
268
#	Line 217 | Line 287 | HiJetResponseAnalyzer::analyze(const edm
287
288		iEvent.getByLabel(jetTag_,jets);
289		if(usePat_)iEvent.getByLabel(jetTag_,patjets);
220	–	if(matchNew_)iEvent.getByLabel(matchTag_,matchedJets);
221	–
290		std::vector<JRAV> jraV;
291
292	<	for(unsigned int j = 0 ; j < jets->size(); ++j){
292	>	if(getFastJets_){
293	>	iEvent.getByLabel(edm::InputTag(ktSrc_.label(),"rhos"),ktRhos);
294	>	iEvent.getByLabel(edm::InputTag(akSrc_.label(),"rhos"),akRhos);
295	>	}
296
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.jtcorpt = jet.pt();
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.jtpt = patjet.correctedJet("raw").pt();
339	<	jv.jtcorpt = patjet.pt();
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();
245	–
346		}else{
347		jv.refpt = -99;
348		jv.refeta = -99;
349		jv.refphi = -99;
350		}
351		}
252	–
253	–	if(matchNew_){
254	–	for(unsigned int m = 0 ; m < matchedJets->size(); ++m){
255	–	const reco::Jet& match = (*matchedJets)[m];
256	–	double dr = reco::deltaR(jet.eta(),jet.phi(),match.eta(),match.phi());
257	–	if(dr < matchR_){
258	–	jv.refcorpt = -99;
259	–	jv.refpt = match.pt();
260	–	jv.refeta = match.eta();
261	–	jv.refphi = match.phi();
262	–	}
263	–	}
264	–	}
352		jraV.push_back(jv);
353		}
354
355		if(sortJets_){
356	<	if(usePat_ || correctJets_) std::sort(jraV.begin(),jraV.end(),compareCorPt);
270	<	else std::sort(jraV.begin(),jraV.end(),comparePt);
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	<	jra_.jtpt[jra_.nref] = jv.jtpt;
362	<	jra_.jteta[jra_.nref] = jv.jteta;
363	<	jra_.jtphi[jra_.nref] = jv.jtphi;
364	<	jra_.jtcorpt[jra_.nref] = jv.jtcorpt;
365	<	jra_.refpt[jra_.nref] = jv.refpt;
366	<	jra_.refeta[jra_.nref] = jv.refeta;
367	<	jra_.refphi[jra_.nref] = jv.refphi;
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	<
434	>
435		t->Fill();
436
437		}
438
289	–
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");
297	<	t->Branch("hf",&jra_.hf,"hf/F");
298	<	t->Branch("nref",&jra_.nref,"nref/I");
299	<	t->Branch("jtpt",jra_.jtpt,"jtpt[nref]/F");
300	<	t->Branch("jtcorpt",jra_.jtcorpt,"jtcorpt[nref]/F");
301	<	t->Branch("refpt",jra_.refpt,"refpt[nref]/F");
302	<	t->Branch("refcorpt",jra_.refpt,"refcorpt[nref]/F");
303	<	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");
305	–	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");
308	–	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