ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

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

Comparing UserCode/CmsHi/JetAnalysis/src/RecHitComparison.cc (file contents):
Revision 1.1 by yilmaz, Thu Sep 23 13:28:28 2010 UTC vs.
Revision 1.5 by yilmaz, Fri Oct 22 13:34:09 2010 UTC

#	Line 52 | Line 52
52		#include "SimDataFormats/HiGenData/interface/GenHIEvent.h"
53		#include "DataFormats/PatCandidates/interface/Jet.h"
54		#include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
55	+	#include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
56	+	#include "DataFormats/HcalDetId/interface/HcalDetId.h"
57
58		#include "TNtuple.h"
59
#	Line 79 | Line 81 | class RecHitComparison : public edm::EDA
81		edm::Handle<vector<double> > ktRhos;
82		edm::Handle<vector<double> > akRhos;
83
84	<	edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > jets1;
85	<	edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > jets2;
84	>	edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > ebHits1;
85	>	edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > ebHits2;
86	>	edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > eeHits1;
87	>	edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > eeHits2;
88	>
89	>	edm::Handle<HFRecHitCollection> hfHits1;
90	>	edm::Handle<HFRecHitCollection> hfHits2;
91	>	edm::Handle<HBHERecHitCollection> hbheHits1;
92	>	edm::Handle<HBHERecHitCollection> hbheHits2;
93
85	–	//   typedef edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > >::const_iterator EcalIterator;
94		typedef vector<EcalRecHit>::const_iterator EcalIterator;
95	+	typedef vector<HFRecHit>::const_iterator HFIterator;
96	+	typedef vector<HBHERecHit>::const_iterator HBHEIterator;
97
98		edm::Handle<reco::CaloJetCollection> signalJets;
99
100	<	TNtuple* nthit;
100	>	edm::InputTag HcalRecHitHFSrc1_;
101	>	edm::InputTag HcalRecHitHFSrc2_;
102	>	edm::InputTag HcalRecHitHBHESrc1_;
103	>	edm::InputTag HcalRecHitHBHESrc2_;
104	>	edm::InputTag EBSrc1_;
105	>	edm::InputTag EBSrc2_;
106	>	edm::InputTag EESrc1_;
107	>	edm::InputTag EESrc2_;
108	>	edm::InputTag signalTag_;
109	>	edm::InputTag centTag_;
110	>
111	>	TNtuple* ntEB;
112	>	TNtuple* ntEE;
113	>	TNtuple* ntHBHE;
114	>	TNtuple* ntHF;
115		TNtuple* ntjet;
116
117		double cone;
118	+	bool jetsOnly_;
119
120		edm::Service<TFileService> fs;
121		const CentralityBins * cbins_;
#	Line 114 | Line 139 | RecHitComparison::RecHitComparison(const
139		cone(0.5)
140		{
141		//now do what ever initialization is needed
142	+	centTag_ =  iConfig.getUntrackedParameter<edm::InputTag>("centrality",edm::InputTag("hiCentrality","","RECO"));
143	+
144	+	jetsOnly_ = iConfig.getUntrackedParameter<bool>("jetsOnly",false);
145	+	signalTag_ = iConfig.getUntrackedParameter<edm::InputTag>("signalJets",edm::InputTag("iterativeCone5CaloJets","","SIGNAL"));
146	+
147	+	HcalRecHitHFSrc1_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHFRecHitSrc1",edm::InputTag("hfreco"));
148	+	HcalRecHitHFSrc2_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHFRecHitSrc2",edm::InputTag("hfreco"));
149	+	HcalRecHitHBHESrc1_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHBHERecHitSrc1",edm::InputTag("hbhereco"));
150	+	HcalRecHitHBHESrc2_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHBHERecHitSrc2",edm::InputTag("hbhereco"));
151	+	EBSrc1_ = iConfig.getUntrackedParameter<edm::InputTag>("EBRecHitSrc1",edm::InputTag("ecalRecHit","EcalRecHitsEB","RECO"));
152	+	EBSrc2_ = iConfig.getUntrackedParameter<edm::InputTag>("EBRecHitSrc2",edm::InputTag("ecalRecHit","EcalRecHitsEB","SIGNALONLY"));
153	+	EESrc1_ = iConfig.getUntrackedParameter<edm::InputTag>("EERecHitSrc1",edm::InputTag("ecalRecHit","EcalRecHitsEE","RECO"));
154	+	EESrc2_ = iConfig.getUntrackedParameter<edm::InputTag>("EERecHitSrc2",edm::InputTag("ecalRecHit","EcalRecHitsEE","SIGNALONLY"));
155
156		}
157
#	Line 143 | Line 181 | RecHitComparison::analyze(const edm::Eve
181		geo = pGeo.product();
182		}
183
146	–	edm::InputTag jetTag1("ecalRecHit","EcalRecHitsEB","SIGNALONLY");
147	–	edm::InputTag jetTag2("ecalRecHit","EcalRecHitsEB","RECO");
148	–
149	–	edm::InputTag signalTag("iterativeConePu5CaloJets","","SIGNALONLY");
150	–	edm::InputTag centTag("hiCentrality","","RECO");
151	–
184		using namespace edm;
185
186	<	ev.getByLabel(centTag,cent);
187	<	ev.getByLabel(jetTag1,jets1);
188	<	ev.getByLabel(jetTag2,jets2);
189	<	ev.getByLabel(signalTag,signalJets);
186	>	ev.getByLabel(centTag_,cent);
187	>	ev.getByLabel(EBSrc1_,ebHits1);
188	>	ev.getByLabel(EBSrc1_,ebHits2);
189	>	ev.getByLabel(signalTag_,signalJets);
190	>
191	>	ev.getByLabel(HcalRecHitHFSrc1_,hfHits1);
192	>	ev.getByLabel(HcalRecHitHFSrc2_,hfHits2);
193	>	ev.getByLabel(HcalRecHitHBHESrc1_,hbheHits1);
194	>	ev.getByLabel(HcalRecHitHBHESrc2_,hbheHits2);
195	>
196	>	ev.getByLabel(EESrc1_,eeHits1);
197	>	ev.getByLabel(EESrc2_,eeHits2);
198
199		double hf = cent->EtHFhitSum();
200		double sumet = cent->EtMidRapiditySum();
#	Line 190 | Line 230 | RecHitComparison::analyze(const edm::Eve
230		f3.push_back(0);
231		}
232
233	<	for(unsigned int j1 = 0 ; j1 < jets1->size(); ++j1){
233	>	for(unsigned int j1 = 0 ; j1 < ebHits1->size(); ++j1){
234
235	<	const EcalRecHit& jet1 = (*jets1)[j1];
235	>	const EcalRecHit& jet1 = (*ebHits1)[j1];
236		double e1 = jet1.energy();
237
238		const GlobalPoint& pos1=geo->getPosition(jet1.id());
#	Line 238 | Line 278 | RecHitComparison::analyze(const edm::Eve
278
279		GlobalPoint pos2;
280		double e2 = -1;
281	<	int siz = jets2->size();
282	<
243	<	EcalIterator hitit = jets2->find(jet1.id());
244	<	if(hitit != jets2->end()){
281	>	EcalIterator hitit = ebHits2->find(jet1.id());
282	>	if(hitit != ebHits2->end()){
283		e2 = hitit->energy();
284		pos2=geo->getPosition(hitit->id());
285		}else{
#	Line 252 | Line 290 | RecHitComparison::analyze(const edm::Eve
290		double eta2 = pos2.eta();
291		double phi2 = pos2.eta();
292		double et2 = e2*sin(pos2.theta());
293	<
294	<	if(isjet) nthit->Fill(e1,et1,e2,et2,eta2,phi2,sumet,hf,bin,jetpt,drjet);
295	<
293	>	if(!jetsOnly_ ||  isjet) ntEB->Fill(e1,et1,e2,et2,eta2,phi2,sumet,hf,bin,jetpt,drjet);
294	>	}
295	>
296	>	for(unsigned int i = 0; i < eeHits1->size(); ++i){
297	>	const EcalRecHit & jet1= (*eeHits1)[i];
298	>	double e1 = jet1.energy();
299	>	const GlobalPoint& pos1=geo->getPosition(jet1.id());
300	>	double eta1 = pos1.eta();
301	>	double phi1 = pos1.phi();
302	>	double et1 = e1*sin(pos1.theta());
303	>	double drjet = -1;
304	>	double jetpt = -1;
305	>	bool isjet = false;
306	>	int matchedJet = -1;
307	>	for(unsigned int j = 0 ; j < signalJets->size(); ++j){
308	>	const reco::CaloJet & jet = (*signalJets)[j];
309	>	double dr = reco::deltaR(eta1,phi1,jet.eta(),jet.phi());
310	>	if(dr < cone){
311	>	jetpt = jet.pt();
312	>	drjet = dr;
313	>	isjet = true;
314	>	matchedJet = j;
315	>	}
316	>	}
317	>	GlobalPoint pos2;
318	>	double e2 = -1;
319	>	EcalIterator hitit = eeHits2->find(jet1.id());
320	>	if(hitit != eeHits2->end()){
321	>	e2 = hitit->energy();
322	>	pos2=geo->getPosition(hitit->id());
323	>	}else{
324	>	e2 = 0;
325	>	pos2 = pos1;
326	>	}
327	>	double eta2 = pos2.eta();
328	>	double phi2 = pos2.eta();
329	>	double et2 = e2*sin(pos2.theta());
330	>	if(!jetsOnly_ || isjet) ntEE->Fill(e1,et1,e2,et2,eta2,phi2,sumet,hf,bin,jetpt,drjet);
331	>	}
332	>
333	>	for(unsigned int i = 0; i < hbheHits1->size(); ++i){
334	>	const HBHERecHit & jet1= (*hbheHits1)[i];
335	>	double e1 = jet1.energy();
336	>	const GlobalPoint& pos1=geo->getPosition(jet1.id());
337	>	double eta1 = pos1.eta();
338	>	double phi1 = pos1.phi();
339	>	double et1 = e1*sin(pos1.theta());
340	>	double drjet = -1;
341	>	double jetpt = -1;
342	>	bool isjet = false;
343	>	int matchedJet = -1;
344	>	for(unsigned int j = 0 ; j < signalJets->size(); ++j){
345	>	const reco::CaloJet & jet = (*signalJets)[j];
346	>	double dr = reco::deltaR(eta1,phi1,jet.eta(),jet.phi());
347	>	if(dr < cone){
348	>	jetpt = jet.pt();
349	>	drjet = dr;
350	>	isjet = true;
351	>	matchedJet = j;
352	>	}
353	>	}
354	>	GlobalPoint pos2;
355	>	double e2 = -1;
356	>	HBHEIterator hitit = hbheHits2->find(jet1.id());
357	>	if(hitit != hbheHits2->end()){
358	>	e2 = hitit->energy();
359	>	pos2=geo->getPosition(hitit->id());
360	>	}else{
361	>	e2 = 0;
362	>	pos2 = pos1;
363	>	}
364	>	double eta2 = pos2.eta();
365	>	double phi2 = pos2.eta();
366	>	double et2 = e2*sin(pos2.theta());
367	>	if(!jetsOnly_ || isjet) ntHBHE->Fill(e1,et1,e2,et2,eta2,phi2,sumet,hf,bin,jetpt,drjet);
368	>	}
369	>
370	>	for(unsigned int i = 0; i < hfHits1->size(); ++i){
371	>	const HFRecHit & jet1= (*hfHits1)[i];
372	>	double e1 = jet1.energy();
373	>	const GlobalPoint& pos1=geo->getPosition(jet1.id());
374	>	double eta1 = pos1.eta();
375	>	double phi1 = pos1.phi();
376	>	double et1 = e1*sin(pos1.theta());
377	>	double drjet = -1;
378	>	double jetpt = -1;
379	>	bool isjet = false;
380	>	int matchedJet = -1;
381	>	for(unsigned int j = 0 ; j < signalJets->size(); ++j){
382	>	const reco::CaloJet & jet = (*signalJets)[j];
383	>	double dr = reco::deltaR(eta1,phi1,jet.eta(),jet.phi());
384	>	if(dr < cone){
385	>	jetpt = jet.pt();
386	>	drjet = dr;
387	>	isjet = true;
388	>	matchedJet = j;
389	>	}
390	>	}
391	>	GlobalPoint pos2;
392	>	double e2 = -1;
393	>	HFIterator hitit = hfHits2->find(jet1.id());
394	>	if(hitit != hfHits2->end()){
395	>	e2 = hitit->energy();
396	>	pos2=geo->getPosition(hitit->id());
397	>	}else{
398	>	e2 = 0;
399	>	pos2 = pos1;
400	>	}
401	>	double eta2 = pos2.eta();
402	>	double phi2 = pos2.eta();
403	>	double et2 = e2*sin(pos2.theta());
404	>	if(!jetsOnly_ || isjet) ntHF->Fill(e1,et1,e2,et2,eta2,phi2,sumet,hf,bin,jetpt,drjet);
405		}
406
407		for(unsigned int j1 = 0 ; j1 < signalJets->size(); ++j1){
#	Line 273 | Line 420 | RecHitComparison::analyze(const edm::Eve
420		void
421		RecHitComparison::beginJob()
422		{
423	<	nthit = fs->make<TNtuple>("nthit","","e1:et1:e2:et2:eta:phi:sumet:hf:bin:ptjet:drjet");
423	>	ntEB = fs->make<TNtuple>("ntEB","","e1:et1:e2:et2:eta:phi:sumet:hf:bin:ptjet:drjet");
424	>	ntEE = fs->make<TNtuple>("ntEE","","e1:et1:e2:et2:eta:phi:sumet:hf:bin:ptjet:drjet");
425	>	ntHBHE = fs->make<TNtuple>("ntHBHE","","e1:et1:e2:et2:eta:phi:sumet:hf:bin:ptjet:drjet");
426	>	ntHF = fs->make<TNtuple>("ntHF","","e1:et1:e2:et2:eta:phi:sumet:hf:bin:ptjet:drjet");
427		ntjet = fs->make<TNtuple>("ntjet","","bin:pt:eta:ethit:f05:f1:f15:f2:f25:f3:em:emf");
428
429		}

Diff Legend

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