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.10 by yilmaz, Thu Jan 20 19:41:10 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;
#	Line 68 | Line 77 | struct JRA{
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;
93		float refpt;
#	Line 83 | Line 96 | struct JRAV{
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
#	Line 114 | Line 132 | class HiJetResponseAnalyzer : public edm
132		bool matchNew_;
133		bool sortJets_;
134		bool correctJets_;
135	+	bool getFastJets_;
136
137		double matchR_;
138		double genPtMin_;
#	Line 122 | Line 141 | class HiJetResponseAnalyzer : public edm
141		double n90Min_;
142		double n90hitMin_;
143
144	<	edm::InputTag jetTag_;
144	>	edm::InputTag jetTag_;
145		edm::InputTag matchTag_;
146	+	std::vector<edm::InputTag> matchTags_;
147	+
148	+	JRA jra_;
149	+	std::vector<JRA> jraMatch_;
150
128	–	JRA jra_;
151		TTree* t;
152
153		edm::Handle<edm::GenHIEvent> mc;
#	Line 133 | Line 155 | class HiJetResponseAnalyzer : public edm
155
156		edm::Handle<reco::JetView> jets;
157		edm::Handle<pat::JetCollection> patjets;
158	<	edm::Handle<reco::JetView> matchedJets;
158	>
159	>	FactorizedJetCorrector* jetCorrector_;
160	>	edm::ESWatcher<JetCorrectionsRecord> watchJetCorrectionsRecord_;
161	>
162	>	std::string tags_;
163	>	std::string levels_;
164	>	std::string algo_;
165
166		edm::Service<TFileService> fs;
167
168	+	edm::Handle<vector<double> > ktRhos;
169	+	edm::Handle<vector<double> > akRhos;
170	+	bool doFastJets_;
171	+
172	+	vector<int> doMatchedFastJets_;
173	+	vector<int> correctMatchedJets_;
174	+
175	+	InputTag ktSrc_;
176	+	InputTag akSrc_;
177	+
178	+
179		};
180
181		bool HiJetResponseAnalyzer::selectJet(int i){
#	Line 169 | Line 208 | bool HiJetResponseAnalyzer::selectJet(in
208		HiJetResponseAnalyzer::HiJetResponseAnalyzer(const edm::ParameterSet& iConfig)
209
210		{
211	+
212	+	levels_ = iConfig.getUntrackedParameter<string>("corrLevels","L2Relative:L3Absolute");
213	+
214	+	algo_ = iConfig.getUntrackedParameter<string>("algo","IC5Calo");
215	+	tags_ = "";
216	+
217	+	string l[2] = {"L2Relative","L3Absolute"};
218	+
219	+	for(int i = 0; i <2; ++i){
220	+	edm::FileInPath fip("CondFormats/JetMETObjects/data/Spring10_"+l[i]+"_"+algo_+".txt");
221	+	tags_ += fip.fullPath();
222	+	if(i < 2 - 1)tags_ +=":";
223	+	}
224	+
225	+	jetCorrector_ = new FactorizedJetCorrector(levels_, tags_);
226	+
227		//now do what ever initialization is needed
228		matchR_ = iConfig.getUntrackedParameter<double>("matchR",0.25);
229
#	Line 191 | Line 246 | HiJetResponseAnalyzer::HiJetResponseAnal
246		sortJets_ = iConfig.getUntrackedParameter<bool>("sortJets",true);
247		correctJets_ = iConfig.getUntrackedParameter<bool>("correctJets",false);
248
249	+	correctMatchedJets_ = iConfig.getUntrackedParameter<std::vector<int> >("correctMatchedJets",std::vector<int>(0));
250	+	doMatchedFastJets_ = iConfig.getUntrackedParameter<std::vector<int> >("doMatchedFastJets",std::vector<int>(0));
251	+
252		jetTag_ = iConfig.getUntrackedParameter<edm::InputTag>("src",edm::InputTag("selectedPatJets"));
253		matchTag_ = iConfig.getUntrackedParameter<edm::InputTag>("match",edm::InputTag("selectedPatJets"));
254	+	matchTags_ = iConfig.getUntrackedParameter<std::vector<edm::InputTag> >("matches",std::vector<edm::InputTag>(0));
255	+
256	+	ktSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("ktSrc",edm::InputTag("kt4CaloJets"));
257	+	akSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("akSrc",edm::InputTag("ak5CaloJets"));
258	+	doFastJets_ = iConfig.getUntrackedParameter<bool>("doFastJets",false);
259	+
260	+	getFastJets_ = true;
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;
#	Line 230 | Line 302 | HiJetResponseAnalyzer::analyze(const edm
302		jv.jteta = jet.eta();
303		jv.jtphi = jet.phi();
304		jv.jtcorpt = 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 = getRho(jv.jteta,*ktRhos);
311	+	double akRho = getRho(jv.jteta,*akRhos);
312	+
313	+	jv.rho = ktRho;
314	+	if(correctJets_){
315	+	if(doFastJets_){
316	+	jv.pu = jet.jetArea()*ktRho;
317	+	pt -= jv.pu;
318	+	}
319	+
320	+	jetCorrector_->setJetEta(jet.eta());
321	+	jetCorrector_->setJetPt(pt);
322	+	//      jetCorrector_->setJetE(jet.energy());
323	+	vector<float> corrs = jetCorrector_->getSubCorrections();
324	+	jv.l2 = corrs[0];
325	+	jv.l3 = corrs[1];
326	+	jv.jtcorpt = pt*jv.l2*jv.l3;
327
328	+	}
329		if(usePat_){
330		const pat::Jet& patjet = (*patjets)[j];
331
#	Line 242 | Line 337 | HiJetResponseAnalyzer::analyze(const edm
337		jv.refpt = patjet.genJet()->pt();
338		jv.refeta = patjet.genJet()->eta();
339		jv.refphi = patjet.genJet()->phi();
245	–
340		}else{
341		jv.refpt = -99;
342		jv.refeta = -99;
343		jv.refphi = -99;
344		}
345		}
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	–	}
346		jraV.push_back(jv);
347		}
348
#	Line 272 | Line 353 | HiJetResponseAnalyzer::analyze(const edm
353
354		for(unsigned int i = 0; i < jraV.size(); ++i){
355		JRAV& jv = jraV[i];
356	<	jra_.jtpt[jra_.nref] = jv.jtpt;
357	<	jra_.jteta[jra_.nref] = jv.jteta;
358	<	jra_.jtphi[jra_.nref] = jv.jtphi;
359	<	jra_.jtcorpt[jra_.nref] = jv.jtcorpt;
360	<	jra_.refpt[jra_.nref] = jv.refpt;
361	<	jra_.refeta[jra_.nref] = jv.refeta;
362	<	jra_.refphi[jra_.nref] = jv.refphi;
356	>	const reco::Jet& jet = (*jets)[jv.index];
357	>
358	>	if(matchNew_){
359	>	for(unsigned int im = 0; im < matchTags_.size(); ++im){
360	>	edm::Handle<reco::JetView> matchedJets;
361	>	iEvent.getByLabel(matchTags_[im],matchedJets);
362	>	jraMatch_[im].jtcorpt[i] = -99;
363	>	jraMatch_[im].jtpt[i] = -99;
364	>	jraMatch_[im].jteta[i] = -99;
365	>	jraMatch_[im].jtphi[i] = -99;
366	>	jraMatch_[im].area[i] = -99;
367	>	jraMatch_[im].l2[i] = -99;
368	>	jraMatch_[im].l3[i] = -99;
369	>	jraMatch_[im].pu[i] = -99;
370	>	for(unsigned int m = 0 ; m < matchedJets->size(); ++m){
371	>	const reco::Jet& match = (*matchedJets)[m];
372	>	double dr = reco::deltaR(jet.eta(),jet.phi(),match.eta(),match.phi());
373	>	if(dr < matchR_ && match.pt() > genPtMin_){
374	>	jraMatch_[im].jtcorpt[i] = -99;
375	>	jraMatch_[im].jtpt[i] = match.pt();
376	>	jraMatch_[im].jteta[i] = match.eta();
377	>	jraMatch_[im].jtphi[i] = match.phi();
378	>	jraMatch_[im].area[i] = match.jetArea();
379	>
380	>	double  ktRho = getRho(jraMatch_[im].jteta[i],*ktRhos);
381	>	double akRho  = getRho(jraMatch_[im].jteta[i],*akRhos);
382	>
383	>	jraMatch_[im].rho[i] = ktRho;
384	>
385	>	if((bool)correctMatchedJets_[im]){
386	>
387	>	double ktRho;
388	>	double pt = jraMatch_[im].jtpt[i];
389	>
390	>	if((bool)doMatchedFastJets_[im]){
391	>	jraMatch_[im].pu[i] = jraMatch_[im].area[i]*ktRho;
392	>	pt -= jraMatch_[im].pu[i];
393	>	}
394	>	jetCorrector_->setJetEta(jraMatch_[im].jteta[i]);
395	>	jetCorrector_->setJetPt(pt);
396	>	//              jetCorrector_->setJetE(match.energy());
397	>	vector<float> corrs = jetCorrector_->getSubCorrections();
398	>	jraMatch_[im].l2[i] = corrs[0];
399	>	// Should it be re-set for L3???
400	>	jraMatch_[im].l3[i] = corrs[1];
401	>	jraMatch_[im].jtcorpt[i] = pt*jraMatch_[im].l2[i]*jraMatch_[im].l3[i];
402	>
403	>	}
404	>
405	>
406	>	}
407	>	}
408	>	}
409	>	}
410	>
411	>	jra_.jtpt[i] = jv.jtpt;
412	>	jra_.jteta[i] = jv.jteta;
413	>	jra_.jtphi[i] = jv.jtphi;
414	>	jra_.jtcorpt[i] = jv.jtcorpt;
415	>	jra_.refpt[i] = jv.refpt;
416	>	jra_.refeta[i] = jv.refeta;
417	>	jra_.refphi[i] = jv.refphi;
418	>
419	>	jra_.area[i] = jv.area;
420	>	jra_.pu[i] = jv.pu;
421	>	jra_.rho[i] = jv.rho;
422	>	jra_.l2[i] = jv.l2;
423	>	jra_.l3[i] = jv.l3;
424	>
425		}
426		jra_.nref = jraV.size();
427	<
427	>
428		t->Fill();
429
430		}
431
289	–
432		// ------------ method called once each job just before starting event loop  ------------
433	+
434	+
435		void
436	<	HiJetResponseAnalyzer::beginJob()
437	<	{
436	>	HiJetResponseAnalyzer::beginJob(){
437	>	t= fs->make<TTree>("t","Jet Response Analyzer");
438	>	t->Branch("nref",&jra_.nref,"nref/I");
439	>	t->Branch("jtpt",jra_.jtpt,"jtpt[nref]/F");
440	>	t->Branch("jteta",jra_.jteta,"jteta[nref]/F");
441	>	t->Branch("jtphi",jra_.jtphi,"jtphi[nref]/F");
442	>
443	>	if(correctJets_){
444	>	t->Branch("jtcorpt",jra_.jtcorpt,"jtcorpt[nref]/F");
445	>	t->Branch("l2",jra_.l2,"l2[nref]/F");
446	>	t->Branch("l3",jra_.l3,"l3[nref]/F");
447	>	}
448	>	t->Branch("area",jra_.area,"area[nref]/F");
449	>	t->Branch("pu",jra_.pu,"pu[nref]/F");
450	>	t->Branch("rho",jra_.rho,"rho[nref]/F");
451
452	<	t= fs->make<TTree>("t","Jet Response Analyzer");
453	<	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");
452	>	t->Branch("refpt",jra_.refpt,"refpt[nref]/F");
453	>	t->Branch("refcorpt",jra_.refpt,"refcorpt[nref]/F");
454		t->Branch("refeta",jra_.refeta,"refeta[nref]/F");
305	–	t->Branch("jtphi",jra_.jtphi,"jtphi[nref]/F");
455		t->Branch("refphi",jra_.refphi,"refphi[nref]/F");
456		t->Branch("weight",&jra_.weight,"weight/F");
308	–	t->Branch("bin",&jra_.bin,"bin/I");
457
458	+	jraMatch_.clear();
459	+	for(unsigned int im = 0; im < matchTags_.size(); ++im){
460	+	JRA jrm;
461	+	jraMatch_.push_back(jrm);
462	+	}
463
464	+	for(unsigned int im = 0; im < matchTags_.size(); ++im){
465	+	t->Branch(Form("jtpt%d",im),jraMatch_[im].jtpt,Form("jtpt%d[nref]/F",im));
466	+	t->Branch(Form("jteta%d",im),jraMatch_[im].jteta,Form("jteta%d[nref]/F",im));
467	+	t->Branch(Form("jtphi%d",im),jraMatch_[im].jtphi,Form("jtphi%d[nref]/F",im));
468	+
469	+	if((bool)correctMatchedJets_[im]){
470	+	t->Branch(Form("jtcorpt%d",im),jraMatch_[im].jtcorpt,Form("jtcorpt%d[nref]/F",im));
471	+	t->Branch(Form("l2_%d",im),jraMatch_[im].l2,Form("l2_%d[nref]/F",im));
472	+	t->Branch(Form("l3_%d",im),jraMatch_[im].l3,Form("l3_%d[nref]/F",im));
473	+	}
474	+
475	+	t->Branch(Form("area%d",im),jraMatch_[im].area,Form("area%d[nref]/F",im));
476	+	t->Branch(Form("pu%d",im),jraMatch_[im].pu,Form("pu%d[nref]/F",im));
477	+	t->Branch(Form("rho%d",im),jraMatch_[im].rho,Form("rho%d[nref]/F",im));
478
479	+	}
480	+
481		}
482
483		// ------------ method called once each job just after ending the event loop  ------------

Diff Legend

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