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.8 by yilmaz, Thu Oct 21 16:24:04 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		};
#	Line 75 | Line 89 | struct JRA{
89		struct JRAV{
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		//
91	–	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 113 | 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 135 | 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 171 | 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 193 | 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 222 | Line 289 | HiJetResponseAnalyzer::analyze(const edm
289		if(usePat_)iEvent.getByLabel(jetTag_,patjets);
290		std::vector<JRAV> jraV;
291
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];
#	Line 229 | Line 301 | HiJetResponseAnalyzer::analyze(const edm
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;
#	Line 248 | Line 349 | HiJetResponseAnalyzer::analyze(const edm
349		jv.refphi = -99;
350		}
351		}
251	–
352		jraV.push_back(jv);
353		}
354
355		if(sortJets_){
356	<	if(usePat_ || correctJets_) std::sort(jraV.begin(),jraV.end(),compareCorPt);
257	<	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){
#	Line 263 | Line 362 | HiJetResponseAnalyzer::analyze(const edm
362
363		if(matchNew_){
364		for(unsigned int im = 0; im < matchTags_.size(); ++im){
365	<	iEvent.getByLabel(matchTags_[im],matchedJets);
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_){
379	<	jraMatch_[im].jtcorpt[i] = -99;
380	<	jraMatch_[im].jtpt[i] = match.pt();
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		}
#	Line 280 | Line 418 | HiJetResponseAnalyzer::analyze(const edm
418		jra_.jtpt[i] = jv.jtpt;
419		jra_.jteta[i] = jv.jteta;
420		jra_.jtphi[i] = jv.jtphi;
421	<	jra_.jtcorpt[i] = jv.jtcorpt;
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
#	Line 297 | Line 442 | HiJetResponseAnalyzer::analyze(const edm
442		void
443		HiJetResponseAnalyzer::beginJob(){
444		t= fs->make<TTree>("t","Jet Response Analyzer");
300	–	t->Branch("b",&jra_.b,"b/F");
301	–	t->Branch("hf",&jra_.hf,"hf/F");
445		t->Branch("nref",&jra_.nref,"nref/I");
446		t->Branch("jtpt",jra_.jtpt,"jtpt[nref]/F");
304	–	t->Branch("jtcorpt",jra_.jtcorpt,"jtcorpt[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->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");
462		t->Branch("refphi",jra_.refphi,"refphi[nref]/F");
463		t->Branch("weight",&jra_.weight,"weight/F");
464	<	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));
317	–	t->Branch(Form("jtcorpt%d",im),jraMatch_[im].jtcorpt,Form("jtcorpt%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		}

Diff Legend

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