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.8 by yilmaz, Tue Nov 2 13:19:06 2010 UTC

#	Line 48 | Line 48
48		#include "FWCore/Utilities/interface/InputTag.h"
49		#include "DataFormats/HeavyIonEvent/interface/CentralityBins.h"
50		#include "DataFormats/CaloTowers/interface/CaloTower.h"
51	<	#include "DataFormats/HeavyIonEvent/interface/Centrality.h"
51	>	#include "DataFormats/HeavyIonEvent/interface/CentralityProvider.h"
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 "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
59
60		#include "TNtuple.h"
61
#	Line 73 | Line 77 | class RecHitComparison : public edm::EDA
77		virtual void endJob() ;
78
79		// ----------member data ---------------------------
76	–
77	–
78	–	edm::Handle<reco::Centrality> cent;
80		edm::Handle<vector<double> > ktRhos;
81		edm::Handle<vector<double> > akRhos;
82
83	<	edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > jets1;
84	<	edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > jets2;
83	>	edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > ebHits1;
84	>	edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > ebHits2;
85	>	edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > eeHits1;
86	>	edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > eeHits2;
87	>
88	>	edm::Handle<HFRecHitCollection> hfHits1;
89	>	edm::Handle<HFRecHitCollection> hfHits2;
90	>	edm::Handle<HBHERecHitCollection> hbheHits1;
91	>	edm::Handle<HBHERecHitCollection> hbheHits2;
92	>
93	>	edm::Handle<reco::BasicClusterCollection> bClusters1;
94	>	edm::Handle<reco::BasicClusterCollection> bClusters2;
95	>
96
85	–	//   typedef edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > >::const_iterator EcalIterator;
97		typedef vector<EcalRecHit>::const_iterator EcalIterator;
98	+	typedef vector<HFRecHit>::const_iterator HFIterator;
99	+	typedef vector<HBHERecHit>::const_iterator HBHEIterator;
100
101		edm::Handle<reco::CaloJetCollection> signalJets;
102
103	<	TNtuple* nthit;
103	>	edm::InputTag HcalRecHitHFSrc1_;
104	>	edm::InputTag HcalRecHitHFSrc2_;
105	>	edm::InputTag HcalRecHitHBHESrc1_;
106	>	edm::InputTag HcalRecHitHBHESrc2_;
107	>	edm::InputTag EBSrc1_;
108	>	edm::InputTag EBSrc2_;
109	>	edm::InputTag EESrc1_;
110	>	edm::InputTag EESrc2_;
111	>
112	>	edm::InputTag BCSrc1_;
113	>	edm::InputTag BCSrc2_;
114	>
115	>	edm::InputTag signalTag_;
116	>
117	>	TNtuple* ntBC;
118	>	TNtuple* ntEB;
119	>	TNtuple* ntEE;
120	>	TNtuple* ntHBHE;
121	>	TNtuple* ntHF;
122		TNtuple* ntjet;
123
124		double cone;
125	+	bool jetsOnly_;
126	+	bool doBasicClusters_;
127
128		edm::Service<TFileService> fs;
129	<	const CentralityBins * cbins_;
129	>	CentralityProvider* centrality_;
130		const CaloGeometry *geo;
131		};
132
#	Line 109 | Line 142 | class RecHitComparison : public edm::EDA
142		// constructors and destructor
143		//
144		RecHitComparison::RecHitComparison(const edm::ParameterSet& iConfig) :
145	<	cbins_(0),
145	>	centrality_(0),
146		geo(0),
147		cone(0.5)
148		{
149		//now do what ever initialization is needed
150	<
150	>	jetsOnly_ = iConfig.getUntrackedParameter<bool>("jetsOnly",false);
151	>	doBasicClusters_ = iConfig.getUntrackedParameter<bool>("doBasicClusters",false);
152	>	signalTag_ = iConfig.getUntrackedParameter<edm::InputTag>("signalJets",edm::InputTag("iterativeCone5CaloJets","","SIGNAL"));
153	>
154	>	HcalRecHitHFSrc1_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHFRecHitSrc1",edm::InputTag("hfreco"));
155	>	HcalRecHitHFSrc2_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHFRecHitSrc2",edm::InputTag("hfreco"));
156	>	HcalRecHitHBHESrc1_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHBHERecHitSrc1",edm::InputTag("hbhereco"));
157	>	HcalRecHitHBHESrc2_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHBHERecHitSrc2",edm::InputTag("hbhereco"));
158	>	EBSrc1_ = iConfig.getUntrackedParameter<edm::InputTag>("EBRecHitSrc1",edm::InputTag("ecalRecHit","EcalRecHitsEB","RECO"));
159	>	EBSrc2_ = iConfig.getUntrackedParameter<edm::InputTag>("EBRecHitSrc2",edm::InputTag("ecalRecHit","EcalRecHitsEB","SIGNALONLY"));
160	>	EESrc1_ = iConfig.getUntrackedParameter<edm::InputTag>("EERecHitSrc1",edm::InputTag("ecalRecHit","EcalRecHitsEE","RECO"));
161	>	EESrc2_ = iConfig.getUntrackedParameter<edm::InputTag>("EERecHitSrc2",edm::InputTag("ecalRecHit","EcalRecHitsEE","SIGNALONLY"));
162	>	BCSrc1_ = iConfig.getUntrackedParameter<edm::InputTag>("BasicClusterSrc1",edm::InputTag("ecalRecHit","EcalRecHitsEB","RECO"));
163	>	BCSrc2_ = iConfig.getUntrackedParameter<edm::InputTag>("BasicClusterSrc2",edm::InputTag("ecalRecHit","EcalRecHitsEB","SIGNALONLY"));
164		}
165
166
#	Line 135 | Line 181 | RecHitComparison::~RecHitComparison()
181		void
182		RecHitComparison::analyze(const edm::Event& ev, const edm::EventSetup& iSetup)
183		{
184	<
139	<	if(!cbins_) cbins_ = getCentralityBinsFromDB(iSetup);
184	>	if(!centrality_) centrality_ = new CentralityProvider(iSetup);
185		if(!geo){
186		edm::ESHandle<CaloGeometry> pGeo;
187		iSetup.get<CaloGeometryRecord>().get(pGeo);
188		geo = pGeo.product();
189		}
190
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	–
191		using namespace edm;
192	+	ev.getByLabel(EBSrc1_,ebHits1);
193	+	ev.getByLabel(EBSrc2_,ebHits2);
194	+	ev.getByLabel(signalTag_,signalJets);
195	+
196	+	ev.getByLabel(HcalRecHitHFSrc1_,hfHits1);
197	+	ev.getByLabel(HcalRecHitHFSrc2_,hfHits2);
198	+	ev.getByLabel(HcalRecHitHBHESrc1_,hbheHits1);
199	+	ev.getByLabel(HcalRecHitHBHESrc2_,hbheHits2);
200	+
201	+	ev.getByLabel(EESrc1_,eeHits1);
202	+	ev.getByLabel(EESrc2_,eeHits2);
203	+
204	+	if(doBasicClusters_){
205	+	ev.getByLabel(BCSrc1_,bClusters1);
206	+	ev.getByLabel(BCSrc2_,bClusters2);
207	+	}
208
209	<	ev.getByLabel(centTag,cent);
210	<	ev.getByLabel(jetTag1,jets1);
211	<	ev.getByLabel(jetTag2,jets2);
157	<	ev.getByLabel(signalTag,signalJets);
158	<
159	<	double hf = cent->EtHFhitSum();
160	<	double sumet = cent->EtMidRapiditySum();
161	<	int run = ev.id().run();
162	<	run = 1;
163	<	int bin = cbins_->getBin(hf);
164	<	int margin = 0;
209	>	centrality_->newEvent(ev,iSetup);
210	>	double hf = centrality_->centralityValue();
211	>	int bin = centrality_->getBin();
212
213		vector<double> fFull;
214		vector<double> f05;
#	Line 190 | Line 237 | RecHitComparison::analyze(const edm::Eve
237		f3.push_back(0);
238		}
239
240	<	for(unsigned int j1 = 0 ; j1 < jets1->size(); ++j1){
240	>	for(unsigned int j1 = 0 ; j1 < ebHits1->size(); ++j1){
241
242	<	const EcalRecHit& jet1 = (*jets1)[j1];
242	>	const EcalRecHit& jet1 = (*ebHits1)[j1];
243		double e1 = jet1.energy();
244
245		const GlobalPoint& pos1=geo->getPosition(jet1.id());
#	Line 238 | Line 285 | RecHitComparison::analyze(const edm::Eve
285
286		GlobalPoint pos2;
287		double e2 = -1;
288	<	int siz = jets2->size();
289	<
243	<	EcalIterator hitit = jets2->find(jet1.id());
244	<	if(hitit != jets2->end()){
288	>	EcalIterator hitit = ebHits2->find(jet1.id());
289	>	if(hitit != ebHits2->end()){
290		e2 = hitit->energy();
291		pos2=geo->getPosition(hitit->id());
292		}else{
#	Line 252 | Line 297 | RecHitComparison::analyze(const edm::Eve
297		double eta2 = pos2.eta();
298		double phi2 = pos2.eta();
299		double et2 = e2*sin(pos2.theta());
300	<
301	<	if(isjet) nthit->Fill(e1,et1,e2,et2,eta2,phi2,sumet,hf,bin,jetpt,drjet);
302	<
300	>	if(!jetsOnly_ ||  isjet) ntEB->Fill(e1,et1,e2,et2,eta2,phi2,hf,bin,jetpt,drjet);
301	>	}
302	>
303	>	for(unsigned int i = 0; i < eeHits1->size(); ++i){
304	>	const EcalRecHit & jet1= (*eeHits1)[i];
305	>	double e1 = jet1.energy();
306	>	const GlobalPoint& pos1=geo->getPosition(jet1.id());
307	>	double eta1 = pos1.eta();
308	>	double phi1 = pos1.phi();
309	>	double et1 = e1*sin(pos1.theta());
310	>	double drjet = -1;
311	>	double jetpt = -1;
312	>	bool isjet = false;
313	>	int matchedJet = -1;
314	>	for(unsigned int j = 0 ; j < signalJets->size(); ++j){
315	>	const reco::CaloJet & jet = (*signalJets)[j];
316	>	double dr = reco::deltaR(eta1,phi1,jet.eta(),jet.phi());
317	>	if(dr < cone){
318	>	jetpt = jet.pt();
319	>	drjet = dr;
320	>	isjet = true;
321	>	matchedJet = j;
322	>	}
323	>	}
324	>	GlobalPoint pos2;
325	>	double e2 = -1;
326	>	EcalIterator hitit = eeHits2->find(jet1.id());
327	>	if(hitit != eeHits2->end()){
328	>	e2 = hitit->energy();
329	>	pos2=geo->getPosition(hitit->id());
330	>	}else{
331	>	e2 = 0;
332	>	pos2 = pos1;
333	>	}
334	>	double eta2 = pos2.eta();
335	>	double phi2 = pos2.eta();
336	>	double et2 = e2*sin(pos2.theta());
337	>	if(!jetsOnly_ || isjet) ntEE->Fill(e1,et1,e2,et2,eta2,phi2,hf,bin,jetpt,drjet);
338	>	}
339	>
340	>	for(unsigned int i = 0; i < hbheHits1->size(); ++i){
341	>	const HBHERecHit & jet1= (*hbheHits1)[i];
342	>	double e1 = jet1.energy();
343	>	const GlobalPoint& pos1=geo->getPosition(jet1.id());
344	>	double eta1 = pos1.eta();
345	>	double phi1 = pos1.phi();
346	>	double et1 = e1*sin(pos1.theta());
347	>	double drjet = -1;
348	>	double jetpt = -1;
349	>	bool isjet = false;
350	>	int matchedJet = -1;
351	>	for(unsigned int j = 0 ; j < signalJets->size(); ++j){
352	>	const reco::CaloJet & jet = (*signalJets)[j];
353	>	double dr = reco::deltaR(eta1,phi1,jet.eta(),jet.phi());
354	>	if(dr < cone){
355	>	jetpt = jet.pt();
356	>	drjet = dr;
357	>	isjet = true;
358	>	matchedJet = j;
359	>	}
360	>	}
361	>	GlobalPoint pos2;
362	>	double e2 = -1;
363	>	HBHEIterator hitit = hbheHits2->find(jet1.id());
364	>	if(hitit != hbheHits2->end()){
365	>	e2 = hitit->energy();
366	>	pos2=geo->getPosition(hitit->id());
367	>	}else{
368	>	e2 = 0;
369	>	pos2 = pos1;
370	>	}
371	>	double eta2 = pos2.eta();
372	>	double phi2 = pos2.eta();
373	>	double et2 = e2*sin(pos2.theta());
374	>	if(!jetsOnly_ || isjet) ntHBHE->Fill(e1,et1,e2,et2,eta2,phi2,hf,bin,jetpt,drjet);
375	>	}
376	>
377	>	for(unsigned int i = 0; i < hfHits1->size(); ++i){
378	>	const HFRecHit & jet1= (*hfHits1)[i];
379	>	double e1 = jet1.energy();
380	>	const GlobalPoint& pos1=geo->getPosition(jet1.id());
381	>	double eta1 = pos1.eta();
382	>	double phi1 = pos1.phi();
383	>	double et1 = e1*sin(pos1.theta());
384	>	double drjet = -1;
385	>	double jetpt = -1;
386	>	bool isjet = false;
387	>	int matchedJet = -1;
388	>	for(unsigned int j = 0 ; j < signalJets->size(); ++j){
389	>	const reco::CaloJet & jet = (*signalJets)[j];
390	>	double dr = reco::deltaR(eta1,phi1,jet.eta(),jet.phi());
391	>	if(dr < cone){
392	>	jetpt = jet.pt();
393	>	drjet = dr;
394	>	isjet = true;
395	>	matchedJet = j;
396	>	}
397	>	}
398	>	GlobalPoint pos2;
399	>	double e2 = -1;
400	>	HFIterator hitit = hfHits2->find(jet1.id());
401	>	if(hitit != hfHits2->end()){
402	>	e2 = hitit->energy();
403	>	pos2=geo->getPosition(hitit->id());
404	>	}else{
405	>	e2 = 0;
406	>	pos2 = pos1;
407	>	}
408	>	double eta2 = pos2.eta();
409	>	double phi2 = pos2.eta();
410	>	double et2 = e2*sin(pos2.theta());
411	>	if(!jetsOnly_ || isjet) ntHF->Fill(e1,et1,e2,et2,eta2,phi2,hf,bin,jetpt,drjet);
412		}
413
414		for(unsigned int j1 = 0 ; j1 < signalJets->size(); ++j1){
#	Line 273 | Line 427 | RecHitComparison::analyze(const edm::Eve
427		void
428		RecHitComparison::beginJob()
429		{
430	<	nthit = fs->make<TNtuple>("nthit","","e1:et1:e2:et2:eta:phi:sumet:hf:bin:ptjet:drjet");
430	>	ntEB = fs->make<TNtuple>("ntEB","","e1:et1:e2:et2:eta:phi:hf:bin:ptjet:drjet");
431	>	ntEE = fs->make<TNtuple>("ntEE","","e1:et1:e2:et2:eta:phi:hf:bin:ptjet:drjet");
432	>	ntHBHE = fs->make<TNtuple>("ntHBHE","","e1:et1:e2:et2:eta:phi:hf:bin:ptjet:drjet");
433	>	ntHF = fs->make<TNtuple>("ntHF","","e1:et1:e2:et2:eta:phi:hf:bin:ptjet:drjet");
434	>
435	>	ntBC = fs->make<TNtuple>("ntBC","","e1:et1:e2:et2:eta:phi:hf:bin:ptjet:drjet");
436	>
437		ntjet = fs->make<TNtuple>("ntjet","","bin:pt:eta:ethit:f05:f1:f15:f2:f25:f3:em:emf");
438
439		}

Diff Legend

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