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.4 by yilmaz, Wed Oct 20 15:51:18 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;
#	Line 114 | Line 138 | RecHitComparison::RecHitComparison(const
138		cone(0.5)
139		{
140		//now do what ever initialization is needed
141	+	centTag_ =  iConfig.getUntrackedParameter<edm::InputTag>("centrality",edm::InputTag("hiCentrality","","RECO"));
142	+
143	+	signalTag_ = iConfig.getUntrackedParameter<edm::InputTag>("signalJets",edm::InputTag("iterativeCone5CaloJets","","SIGNAL"));
144	+
145	+	HcalRecHitHFSrc1_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHFRecHitSrc1",edm::InputTag("hfreco"));
146	+	HcalRecHitHFSrc2_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHFRecHitSrc2",edm::InputTag("hfreco"));
147	+	HcalRecHitHBHESrc1_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHBHERecHitSrc1",edm::InputTag("hbhereco"));
148	+	HcalRecHitHBHESrc2_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHBHERecHitSrc2",edm::InputTag("hbhereco"));
149	+	EBSrc1_ = iConfig.getUntrackedParameter<edm::InputTag>("EBRecHitSrc1",edm::InputTag("ecalRecHit","EcalRecHitsEB","RECO"));
150	+	EBSrc2_ = iConfig.getUntrackedParameter<edm::InputTag>("EBRecHitSrc2",edm::InputTag("ecalRecHit","EcalRecHitsEB","SIGNALONLY"));
151	+	EESrc1_ = iConfig.getUntrackedParameter<edm::InputTag>("EERecHitSrc1",edm::InputTag("ecalRecHit","EcalRecHitsEE","RECO"));
152	+	EESrc2_ = iConfig.getUntrackedParameter<edm::InputTag>("EERecHitSrc2",edm::InputTag("ecalRecHit","EcalRecHitsEE","SIGNALONLY"));
153
154		}
155
#	Line 143 | Line 179 | RecHitComparison::analyze(const edm::Eve
179		geo = pGeo.product();
180		}
181
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	–
182		using namespace edm;
183
184	<	ev.getByLabel(centTag,cent);
185	<	ev.getByLabel(jetTag1,jets1);
186	<	ev.getByLabel(jetTag2,jets2);
187	<	ev.getByLabel(signalTag,signalJets);
184	>	ev.getByLabel(centTag_,cent);
185	>	ev.getByLabel(EBSrc1_,ebHits1);
186	>	ev.getByLabel(EBSrc1_,ebHits2);
187	>	ev.getByLabel(signalTag_,signalJets);
188	>
189	>	ev.getByLabel(HcalRecHitHFSrc1_,hfHits1);
190	>	ev.getByLabel(HcalRecHitHFSrc2_,hfHits2);
191	>	ev.getByLabel(HcalRecHitHBHESrc1_,hbheHits1);
192	>	ev.getByLabel(HcalRecHitHBHESrc2_,hbheHits2);
193	>
194	>	ev.getByLabel(EESrc1_,eeHits1);
195	>	ev.getByLabel(EESrc2_,eeHits2);
196
197		double hf = cent->EtHFhitSum();
198		double sumet = cent->EtMidRapiditySum();
#	Line 190 | Line 228 | RecHitComparison::analyze(const edm::Eve
228		f3.push_back(0);
229		}
230
231	<	for(unsigned int j1 = 0 ; j1 < jets1->size(); ++j1){
231	>	for(unsigned int j1 = 0 ; j1 < ebHits1->size(); ++j1){
232
233	<	const EcalRecHit& jet1 = (*jets1)[j1];
233	>	const EcalRecHit& jet1 = (*ebHits1)[j1];
234		double e1 = jet1.energy();
235
236		const GlobalPoint& pos1=geo->getPosition(jet1.id());
#	Line 238 | Line 276 | RecHitComparison::analyze(const edm::Eve
276
277		GlobalPoint pos2;
278		double e2 = -1;
279	<	int siz = jets2->size();
280	<
243	<	EcalIterator hitit = jets2->find(jet1.id());
244	<	if(hitit != jets2->end()){
279	>	EcalIterator hitit = ebHits2->find(jet1.id());
280	>	if(hitit != ebHits2->end()){
281		e2 = hitit->energy();
282		pos2=geo->getPosition(hitit->id());
283		}else{
#	Line 252 | Line 288 | RecHitComparison::analyze(const edm::Eve
288		double eta2 = pos2.eta();
289		double phi2 = pos2.eta();
290		double et2 = e2*sin(pos2.theta());
291	<
292	<	if(isjet) nthit->Fill(e1,et1,e2,et2,eta2,phi2,sumet,hf,bin,jetpt,drjet);
293	<
291	>	if(isjet) ntEB->Fill(e1,et1,e2,et2,eta2,phi2,sumet,hf,bin,jetpt,drjet);
292	>	}
293	>
294	>	for(unsigned int i = 0; i < eeHits1->size(); ++i){
295	>	const EcalRecHit & jet1= (*eeHits1)[i];
296	>	double e1 = jet1.energy();
297	>	const GlobalPoint& pos1=geo->getPosition(jet1.id());
298	>	double eta1 = pos1.eta();
299	>	double phi1 = pos1.phi();
300	>	double et1 = e1*sin(pos1.theta());
301	>	double drjet = -1;
302	>	double jetpt = -1;
303	>	bool isjet = false;
304	>	int matchedJet = -1;
305	>	for(unsigned int j = 0 ; j < signalJets->size(); ++j){
306	>	const reco::CaloJet & jet = (*signalJets)[j];
307	>	double dr = reco::deltaR(eta1,phi1,jet.eta(),jet.phi());
308	>	if(dr < cone){
309	>	jetpt = jet.pt();
310	>	drjet = dr;
311	>	isjet = true;
312	>	matchedJet = j;
313	>	}
314	>	}
315	>	GlobalPoint pos2;
316	>	double e2 = -1;
317	>	EcalIterator hitit = eeHits2->find(jet1.id());
318	>	if(hitit != eeHits2->end()){
319	>	e2 = hitit->energy();
320	>	pos2=geo->getPosition(hitit->id());
321	>	}else{
322	>	e2 = 0;
323	>	pos2 = pos1;
324	>	}
325	>	double eta2 = pos2.eta();
326	>	double phi2 = pos2.eta();
327	>	double et2 = e2*sin(pos2.theta());
328	>	if(isjet) ntEE->Fill(e1,et1,e2,et2,eta2,phi2,sumet,hf,bin,jetpt,drjet);
329	>	}
330	>
331	>	for(unsigned int i = 0; i < hbheHits1->size(); ++i){
332	>	const HBHERecHit & jet1= (*hbheHits1)[i];
333	>	double e1 = jet1.energy();
334	>	const GlobalPoint& pos1=geo->getPosition(jet1.id());
335	>	double eta1 = pos1.eta();
336	>	double phi1 = pos1.phi();
337	>	double et1 = e1*sin(pos1.theta());
338	>	double drjet = -1;
339	>	double jetpt = -1;
340	>	bool isjet = false;
341	>	int matchedJet = -1;
342	>	for(unsigned int j = 0 ; j < signalJets->size(); ++j){
343	>	const reco::CaloJet & jet = (*signalJets)[j];
344	>	double dr = reco::deltaR(eta1,phi1,jet.eta(),jet.phi());
345	>	if(dr < cone){
346	>	jetpt = jet.pt();
347	>	drjet = dr;
348	>	isjet = true;
349	>	matchedJet = j;
350	>	}
351	>	}
352	>	GlobalPoint pos2;
353	>	double e2 = -1;
354	>	HBHEIterator hitit = hbheHits2->find(jet1.id());
355	>	if(hitit != hbheHits2->end()){
356	>	e2 = hitit->energy();
357	>	pos2=geo->getPosition(hitit->id());
358	>	}else{
359	>	e2 = 0;
360	>	pos2 = pos1;
361	>	}
362	>	double eta2 = pos2.eta();
363	>	double phi2 = pos2.eta();
364	>	double et2 = e2*sin(pos2.theta());
365	>	if(isjet) ntHBHE->Fill(e1,et1,e2,et2,eta2,phi2,sumet,hf,bin,jetpt,drjet);
366	>	}
367	>
368	>	for(unsigned int i = 0; i < hfHits1->size(); ++i){
369	>	const HFRecHit & jet1= (*hfHits1)[i];
370	>	double e1 = jet1.energy();
371	>	const GlobalPoint& pos1=geo->getPosition(jet1.id());
372	>	double eta1 = pos1.eta();
373	>	double phi1 = pos1.phi();
374	>	double et1 = e1*sin(pos1.theta());
375	>	double drjet = -1;
376	>	double jetpt = -1;
377	>	bool isjet = false;
378	>	int matchedJet = -1;
379	>	for(unsigned int j = 0 ; j < signalJets->size(); ++j){
380	>	const reco::CaloJet & jet = (*signalJets)[j];
381	>	double dr = reco::deltaR(eta1,phi1,jet.eta(),jet.phi());
382	>	if(dr < cone){
383	>	jetpt = jet.pt();
384	>	drjet = dr;
385	>	isjet = true;
386	>	matchedJet = j;
387	>	}
388	>	}
389	>	GlobalPoint pos2;
390	>	double e2 = -1;
391	>	HFIterator hitit = hfHits2->find(jet1.id());
392	>	if(hitit != hfHits2->end()){
393	>	e2 = hitit->energy();
394	>	pos2=geo->getPosition(hitit->id());
395	>	}else{
396	>	e2 = 0;
397	>	pos2 = pos1;
398	>	}
399	>	double eta2 = pos2.eta();
400	>	double phi2 = pos2.eta();
401	>	double et2 = e2*sin(pos2.theta());
402	>	if(isjet) ntHF->Fill(e1,et1,e2,et2,eta2,phi2,sumet,hf,bin,jetpt,drjet);
403		}
404
405		for(unsigned int j1 = 0 ; j1 < signalJets->size(); ++j1){
#	Line 273 | Line 418 | RecHitComparison::analyze(const edm::Eve
418		void
419		RecHitComparison::beginJob()
420		{
421	<	nthit = fs->make<TNtuple>("nthit","","e1:et1:e2:et2:eta:phi:sumet:hf:bin:ptjet:drjet");
421	>	ntEB = fs->make<TNtuple>("ntEB","","e1:et1:e2:et2:eta:phi:sumet:hf:bin:ptjet:drjet");
422	>	ntEE = fs->make<TNtuple>("ntEE","","e1:et1:e2:et2:eta:phi:sumet:hf:bin:ptjet:drjet");
423	>	ntHBHE = fs->make<TNtuple>("ntHBHE","","e1:et1:e2:et2:eta:phi:sumet:hf:bin:ptjet:drjet");
424	>	ntHF = fs->make<TNtuple>("ntHF","","e1:et1:e2:et2:eta:phi:sumet:hf:bin:ptjet:drjet");
425		ntjet = fs->make<TNtuple>("ntjet","","bin:pt:eta:ethit:f05:f1:f15:f2:f25:f3:em:emf");
426
427		}

Diff Legend

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