ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

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

Comparing UserCode/CmsHi/JetAnalysis/src/RecHitTreeProducer.cc (file contents):
Revision 1.1 by yilmaz, Sat Oct 16 13:46:13 2010 UTC vs.
Revision 1.9 by yilmaz, Tue Nov 9 22:54:45 2010 UTC

#	Line 57 | Line 57
57		#include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
58		#include "DataFormats/HcalDetId/interface/HcalDetId.h"
59		#include "DataFormats/EcalDetId/interface/EcalDetIdCollections.h"
60	+	#include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
61		#include "DataFormats/DetId/interface/DetId.h"
62
63
#	Line 64 | Line 65
65
66		using namespace std;
67
68	<	#define MAXHITS 100000
68	>	#define MAXHITS 1000000
69
70		struct MyRecHit{
71	<
71	>	int depth[MAXHITS];
72		int n;
73
74		float e[MAXHITS];
75		float et[MAXHITS];
76		float eta[MAXHITS];
77		float phi[MAXHITS];
78	+	bool isjet[MAXHITS];
79	+
80	+	float jtpt;
81	+	float jteta;
82	+	float jtphi;
83
84		};
85
86
87	+	struct MyBkg{
88	+	int n;
89	+	float rho[50];
90	+	float sigma[50];
91	+	};
92	+
93
94		//
95		// class declaration
#	Line 99 | Line 111 | class RecHitTreeProducer : public edm::E
111
112		// ----------member data ---------------------------
113
102	–
114		edm::Handle<reco::Centrality> cent;
115		edm::Handle<vector<double> > ktRhos;
116		edm::Handle<vector<double> > akRhos;
117
118	<	edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > jets1;
119	<	edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > jets2;
118	>	edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > ebHits;
119	>	edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > eeHits;
120	>
121	>	edm::Handle<HFRecHitCollection> hfHits;
122	>	edm::Handle<HBHERecHitCollection> hbheHits;
123	>
124	>	edm::Handle<reco::BasicClusterCollection> bClusters;
125	>	edm::Handle<CaloTowerCollection> towers;
126
110	–	//   typedef edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > >::const_iterator EcalIterator;
111	–	typedef vector<EcalRecHit>::const_iterator EcalIterator;
127
128	<	edm::Handle<reco::CaloJetCollection> signalJets;
128	>	typedef vector<EcalRecHit>::const_iterator EcalIterator;
129	>
130	>	edm::Handle<reco::CaloJetCollection> jets;
131
132	+	edm::Handle<std::vector<double> > rhos;
133	+	edm::Handle<std::vector<double> > sigmas;
134	+
135		MyRecHit hbheRecHit;
136		MyRecHit hfRecHit;
137	+	MyRecHit ebRecHit;
138	+	MyRecHit eeRecHit;
139	+	MyRecHit myBC;
140	+	MyRecHit myTowers;
141	+	MyBkg bkg;
142
143	+	TNtuple* nt;
144		TTree* hbheTree;
145		TTree* hfTree;
146	<
147	<	double cone;
146	>	TTree* ebTree;
147	>	TTree* eeTree;
148	>	TTree* bcTree;
149	>	TTree* towerTree;
150	>	TTree* bkgTree;
151	>
152	>	double cone;
153	>	double hfTowerThreshold_;
154	>	double hfLongThreshold_;
155	>	double hfShortThreshold_;
156	>	double hbheThreshold_;
157	>	double ebThreshold_;
158	>	double eeThreshold_;
159
160		edm::Service<TFileService> fs;
161		const CentralityBins * cbins_;
#	Line 126 | Line 163 | class RecHitTreeProducer : public edm::E
163
164		edm::InputTag HcalRecHitHFSrc_;
165		edm::InputTag HcalRecHitHBHESrc_;
166	+	edm::InputTag EBSrc_;
167	+	edm::InputTag EESrc_;
168	+	edm::InputTag BCSrc_;
169	+	edm::InputTag TowerSrc_;
170	+
171	+	edm::InputTag JetSrc_;
172
173	+	edm::InputTag FastJetTag_;
174
175	+	bool useJets_;
176	+	bool doBasicClusters_;
177	+	bool doTowers_;
178	+	bool doEcal_;
179	+	bool doHcal_;
180	+
181	+	bool doFastJet_;
182	+
183	+	bool doEbyEonly_;
184		};
185
186		//
#	Line 147 | Line 200 | RecHitTreeProducer::RecHitTreeProducer(c
200		cone(0.5)
201		{
202		//now do what ever initialization is needed
203	+	EBSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("EBRecHitSrc",edm::InputTag("ecalRecHit","EcalRecHitsEB"));
204	+	EESrc_ = iConfig.getUntrackedParameter<edm::InputTag>("EERecHitSrc",edm::InputTag("ecalRecHit","EcalRecHitsEE"));
205		HcalRecHitHFSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHFRecHitSrc",edm::InputTag("hfreco"));
206		HcalRecHitHBHESrc_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHBHERecHitSrc",edm::InputTag("hbhereco"));
207	+	BCSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("BasicClusterSrc1",edm::InputTag("ecalRecHit","EcalRecHitsEB","RECO"));
208	+	TowerSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("towersSrc",edm::InputTag("towerMaker"));
209	+	JetSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("JetSrc",edm::InputTag("iterativeCone5CaloJets"));
210	+	useJets_ = iConfig.getUntrackedParameter<bool>("useJets",true);
211	+	doBasicClusters_ = iConfig.getUntrackedParameter<bool>("doBasicClusters",false);
212	+	doTowers_ = iConfig.getUntrackedParameter<bool>("doTowers",true);
213	+	doEcal_ = iConfig.getUntrackedParameter<bool>("doEcal",true);
214	+	doHcal_ = iConfig.getUntrackedParameter<bool>("doHcal",true);
215	+	doFastJet_ = iConfig.getUntrackedParameter<bool>("doFastJet",true);
216	+	FastJetTag_ = iConfig.getUntrackedParameter<edm::InputTag>("FastJetTag",edm::InputTag("kt4CaloJets"));
217	+	doEbyEonly_ = iConfig.getUntrackedParameter<bool>("doEbyEonly",false);
218	+	hfTowerThreshold_ = iConfig.getUntrackedParameter<double>("HFtowerMin",3.);
219	+	hfLongThreshold_ = iConfig.getUntrackedParameter<double>("HFlongMin",3.);
220	+	hfShortThreshold_ = iConfig.getUntrackedParameter<double>("HFshortMin",3.);
221		}
222
223
#	Line 172 | Line 241 | RecHitTreeProducer::analyze(const edm::E
241
242		hfRecHit.n = 0;
243		hbheRecHit.n = 0;
244	<
245	<	edm::Handle<HFRecHitCollection> hfHits;
246	<	edm::Handle<HBHERecHitCollection> hbheHits;
247	<
244	>	ebRecHit.n = 0;
245	>	eeRecHit.n = 0;
246	>	myBC.n = 0;
247	>	myTowers.n = 0;
248	>	bkg.n = 0;
249	>
250	>	if(doEcal_){
251	>	ev.getByLabel(EBSrc_,ebHits);
252	>	ev.getByLabel(EESrc_,eeHits);
253	>	}
254	>	if(doHcal_){
255		ev.getByLabel(HcalRecHitHFSrc_,hfHits);
256		ev.getByLabel(HcalRecHitHBHESrc_,hbheHits);
257	<
258	<
259	<	if(0 && !cbins_) cbins_ = getCentralityBinsFromDB(iSetup);
257	>	}
258	>	if(useJets_) {
259	>	ev.getByLabel(JetSrc_,jets);
260	>	}
261	>
262	>	if(doBasicClusters_){
263	>	ev.getByLabel(BCSrc_,bClusters);
264	>	}
265	>
266	>	if(doTowers_){
267	>	ev.getByLabel(TowerSrc_,towers);
268	>	}
269	>
270	>	if(doFastJet_){
271	>	ev.getByLabel(edm::InputTag(FastJetTag_.label(),"rhos",FastJetTag_.process()),rhos);
272	>	ev.getByLabel(edm::InputTag(FastJetTag_.label(),"sigmas",FastJetTag_.process()),sigmas);
273	>	bkg.n = rhos->size();
274	>	for(unsigned int i = 0; i < rhos->size(); ++i){
275	>	bkg.rho[i] = (*rhos)[i];
276	>	bkg.sigma[i] = (*sigmas)[i];
277	>	}
278	>	}
279	>
280	>	if(0 && !cbins_) cbins_ = getCentralityBinsFromDB(iSetup);
281
282		if(!geo){
283		edm::ESHandle<CaloGeometry> pGeo;
#	Line 188 | Line 285 | RecHitTreeProducer::analyze(const edm::E
285		geo = pGeo.product();
286		}
287
288	+	int nHFlong = 0;
289	+	int nHFshort = 0;
290	+	int nHFtower = 0;
291	+
292	+	if(doHcal_){
293		for(unsigned int i = 0; i < hfHits->size(); ++i){
294		const HFRecHit & hit= (*hfHits)[i];
295		hfRecHit.e[hfRecHit.n] = hit.energy();
296		hfRecHit.et[hfRecHit.n] = getEt(hit.id(),hit.energy());
297		hfRecHit.eta[hfRecHit.n] = getEta(hit.id());
298		hfRecHit.phi[hfRecHit.n] = getPhi(hit.id());
299	<	//     hfRecHit.ieta[hfRecHit.n] = hit.id().ieta();
300	<	//     hfRecHit.iphi[hfRecHit.n] = hit.id().iphi();
301	<
299	>	hfRecHit.isjet[hfRecHit.n] = false;
300	>	hfRecHit.depth[hfRecHit.n] = hit.id().depth();
301	>	if(hit.energy() > hfShortThreshold_ && hit.id().depth() != 1) nHFshort++;
302	>	if(hit.energy() > hfLongThreshold_ && hit.id().depth() == 1) nHFlong++;
303	>
304	>	if(useJets_){
305	>	for(unsigned int j = 0 ; j < jets->size(); ++j){
306	>	const reco::Jet& jet = (*jets)[j];
307	>	double dr = reco::deltaR(hfRecHit.eta[hfRecHit.n],hfRecHit.phi[hfRecHit.n],jet.eta(),jet.phi());
308	>	if(dr < cone){ hfRecHit.isjet[hfRecHit.n] = true; }
309	>	}
310	>	}
311		hfRecHit.n++;
312		}
313	<
313	>	if(!doEbyEonly_){
314		for(unsigned int i = 0; i < hbheHits->size(); ++i){
315		const HBHERecHit & hit= (*hbheHits)[i];
316	<	hbheRecHit.e[hfRecHit.n] = hit.energy();
317	<	hbheRecHit.et[hfRecHit.n] = getEt(hit.id(),hit.energy());
318	<	hbheRecHit.eta[hfRecHit.n] = getEta(hit.id());
319	<	hbheRecHit.phi[hfRecHit.n] = getPhi(hit.id());
320	<	//     hbheRecHit.ieta[hfRecHit.n] = hit.id().ieta();
321	<	//     hbheRecHit.iphi[hfRecHit.n] = hit.id().iphi();
322	<
316	>	hbheRecHit.e[hbheRecHit.n] = hit.energy();
317	>	hbheRecHit.et[hbheRecHit.n] = getEt(hit.id(),hit.energy());
318	>	hbheRecHit.eta[hbheRecHit.n] = getEta(hit.id());
319	>	hbheRecHit.phi[hbheRecHit.n] = getPhi(hit.id());
320	>	hbheRecHit.isjet[hbheRecHit.n] = false;
321	>	if(useJets_){
322	>	for(unsigned int j = 0 ; j < jets->size(); ++j){
323	>	const reco::Jet& jet = (*jets)[j];
324	>	double dr = reco::deltaR(hbheRecHit.eta[hbheRecHit.n],hbheRecHit.phi[hbheRecHit.n],jet.eta(),jet.phi());
325	>	if(dr < cone){ hbheRecHit.isjet[hbheRecHit.n] = true; }
326	>	}
327	>	}
328		hbheRecHit.n++;
329		}
330	+	}
331	+	}
332	+	if(doEcal_ && !doEbyEonly_){
333	+	for(unsigned int i = 0; i < ebHits->size(); ++i){
334	+	const EcalRecHit & hit= (*ebHits)[i];
335	+	ebRecHit.e[ebRecHit.n] = hit.energy();
336	+	ebRecHit.et[ebRecHit.n] = getEt(hit.id(),hit.energy());
337	+	ebRecHit.eta[ebRecHit.n] = getEta(hit.id());
338	+	ebRecHit.phi[ebRecHit.n] = getPhi(hit.id());
339	+	ebRecHit.isjet[ebRecHit.n] = false;
340	+	if(useJets_){
341	+	for(unsigned int j = 0 ; j < jets->size(); ++j){
342	+	const reco::Jet& jet = (*jets)[j];
343	+	double dr = reco::deltaR(ebRecHit.eta[ebRecHit.n],ebRecHit.phi[ebRecHit.n],jet.eta(),jet.phi());
344	+	if(dr < cone){ ebRecHit.isjet[ebRecHit.n] = true; }
345	+	}
346	+	}
347	+	ebRecHit.n++;
348	+	}
349	+
350	+	for(unsigned int i = 0; i < eeHits->size(); ++i){
351	+	const EcalRecHit & hit= (*eeHits)[i];
352	+	eeRecHit.e[eeRecHit.n] = hit.energy();
353	+	eeRecHit.et[eeRecHit.n] = getEt(hit.id(),hit.energy());
354	+	eeRecHit.eta[eeRecHit.n] = getEta(hit.id());
355	+	eeRecHit.phi[eeRecHit.n] = getPhi(hit.id());
356	+	eeRecHit.isjet[eeRecHit.n] = false;
357	+	if(useJets_){
358	+	for(unsigned int j = 0 ; j < jets->size(); ++j){
359	+	const reco::Jet& jet = (*jets)[j];
360	+	double dr = reco::deltaR(eeRecHit.eta[eeRecHit.n],eeRecHit.phi[eeRecHit.n],jet.eta(),jet.phi());
361	+	if(dr < cone){ eeRecHit.isjet[eeRecHit.n] = true; }
362	+	}
363	+	}
364	+	eeRecHit.n++;
365	+	}
366	+	}
367
368	<	hbheTree->Fill();
369	<	hfTree->Fill();
368	>	if(doTowers_){
369	>
370	>	for(unsigned int i = 0; i < towers->size(); ++i){
371	>	const CaloTower & hit= (*towers)[i];
372	>	myTowers.e[myTowers.n] = hit.energy();
373	>	myTowers.et[myTowers.n] = getEt(hit.id(),hit.energy());
374	>	myTowers.eta[myTowers.n] = getEta(hit.id());
375	>	myTowers.phi[myTowers.n] = getPhi(hit.id());
376	>	myTowers.isjet[myTowers.n] = false;
377	>
378	>	if(hit.ieta() > 29 && hit.energy() > hfTowerThreshold_) nHFtower++;
379	>
380	>	if(useJets_){
381	>	for(unsigned int j = 0 ; j < jets->size(); ++j){
382	>	const reco::Jet& jet = (*jets)[j];
383	>	double dr = reco::deltaR(myTowers.eta[myTowers.n],myTowers.phi[myTowers.n],jet.eta(),jet.phi());
384	>	if(dr < cone){ myTowers.isjet[myTowers.n] = true; }
385	>	}
386	>	}
387	>	myTowers.n++;
388	>	}
389	>
390	>	}
391	>
392	>	if(doBasicClusters_ && !doEbyEonly_){
393	>	for(unsigned int j = 0 ; j < jets->size(); ++j){
394	>	const reco::Jet& jet = (*jets)[j];
395	>	myBC.n = 0;
396	>	myBC.jtpt = jet.pt();
397	>	myBC.jteta = jet.eta();
398	>	myBC.jtphi = jet.phi();
399	>
400	>	for(unsigned int i = 0; i < bClusters->size(); ++i){
401	>	const reco::BasicCluster & bc= (*bClusters)[i];
402	>	double dr = reco::deltaR(bc.eta(),bc.phi(),jet.eta(),jet.phi());
403	>	if(dr < cone){
404	>	myBC.e[myBC.n] = bc.energy();
405	>	myBC.et[myBC.n] = bc.energy()*sin(bc.position().theta());
406	>	myBC.eta[myBC.n] = bc.eta();
407	>	myBC.phi[myBC.n] = bc.phi();
408	>	myBC.n++;
409	>	}
410	>	}
411	>	bcTree->Fill();
412	>	}
413	>	}
414
415	+	if(!doEbyEonly_){
416	+	towerTree->Fill();
417	+
418	+	eeTree->Fill();
419	+	ebTree->Fill();
420	+
421	+	hbheTree->Fill();
422	+	hfTree->Fill();
423	+
424	+	if (doFastJet_) {
425	+	bkgTree->Fill();
426	+	}
427	+	}
428	+
429	+	nt->Fill(nHFtower,nHFlong,nHFshort);
430	+
431		}
432
433
#	Line 222 | Line 435 | RecHitTreeProducer::analyze(const edm::E
435		void
436		RecHitTreeProducer::beginJob()
437		{
438	<
438	>
439		hbheTree = fs->make<TTree>("hbhe","");
440		hbheTree->Branch("n",&hbheRecHit.n,"n/I");
441		hbheTree->Branch("e",hbheRecHit.e,"e[n]/F");
442		hbheTree->Branch("et",hbheRecHit.et,"et[n]/F");
443		hbheTree->Branch("eta",hbheRecHit.eta,"eta[n]/F");
444		hbheTree->Branch("phi",hbheRecHit.phi,"phi[n]/F");
445	<
445	>	hbheTree->Branch("isjet",hbheRecHit.isjet,"isjet[n]/O");
446	>
447		hfTree = fs->make<TTree>("hf","");
448		hfTree->Branch("n",&hfRecHit.n,"n/I");
449		hfTree->Branch("e",hfRecHit.e,"e[n]/F");
450		hfTree->Branch("et",hfRecHit.et,"et[n]/F");
451		hfTree->Branch("eta",hfRecHit.eta,"eta[n]/F");
452		hfTree->Branch("phi",hfRecHit.phi,"phi[n]/F");
453	+	hfTree->Branch("depth",hfRecHit.depth,"depth[n]/I");
454	+	hfTree->Branch("isjet",hfRecHit.isjet,"isjet[n]/O");
455	+
456	+	eeTree = fs->make<TTree>("ee","");
457	+	eeTree->Branch("n",&eeRecHit.n,"n/I");
458	+	eeTree->Branch("e",eeRecHit.e,"e[n]/F");
459	+	eeTree->Branch("et",eeRecHit.et,"et[n]/F");
460	+	eeTree->Branch("eta",eeRecHit.eta,"eta[n]/F");
461	+	eeTree->Branch("phi",eeRecHit.phi,"phi[n]/F");
462	+	eeTree->Branch("isjet",eeRecHit.isjet,"isjet[n]/O");
463	+
464	+	ebTree = fs->make<TTree>("eb","");
465	+	ebTree->Branch("n",&ebRecHit.n,"n/I");
466	+	ebTree->Branch("e",ebRecHit.e,"e[n]/F");
467	+	ebTree->Branch("et",ebRecHit.et,"et[n]/F");
468	+	ebTree->Branch("eta",ebRecHit.eta,"eta[n]/F");
469	+	ebTree->Branch("phi",ebRecHit.phi,"phi[n]/F");
470	+	ebTree->Branch("isjet",ebRecHit.isjet,"isjet[n]/O");
471	+
472	+	towerTree = fs->make<TTree>("tower","");
473	+	towerTree->Branch("n",&myTowers.n,"n/I");
474	+	towerTree->Branch("e",myTowers.e,"e[n]/F");
475	+	towerTree->Branch("et",myTowers.et,"et[n]/F");
476	+	towerTree->Branch("eta",myTowers.eta,"eta[n]/F");
477	+	towerTree->Branch("phi",myTowers.phi,"phi[n]/F");
478	+	towerTree->Branch("isjet",myTowers.isjet,"isjet[n]/O");
479	+
480	+
481	+	if(doBasicClusters_){
482	+	bcTree = fs->make<TTree>("bc","");
483	+	bcTree->Branch("n",&myBC.n,"n/I");
484	+	bcTree->Branch("e",myBC.e,"e[n]/F");
485	+	bcTree->Branch("et",myBC.et,"et[n]/F");
486	+	bcTree->Branch("eta",myBC.eta,"eta[n]/F");
487	+	bcTree->Branch("phi",myBC.phi,"phi[n]/F");
488	+	bcTree->Branch("jtpt",&myBC.jtpt,"jtpt/F");
489	+	bcTree->Branch("jteta",&myBC.jteta,"jteta/F");
490	+	bcTree->Branch("jtphi",&myBC.jtphi,"jtphi/F");
491	+	//     bcTree->Branch("isjet",bcRecHit.isjet,"isjet[n]/O");
492	+	}
493	+
494	+	if(doFastJet_){
495	+	bkgTree = fs->make<TTree>("bkg","");
496	+	bkgTree->Branch("n",&bkg.n,"n/I");
497	+	bkgTree->Branch("rho",bkg.rho,"rho[n]/F");
498	+	bkgTree->Branch("sigma",bkg.sigma,"sigma[n]/F");
499	+	}
500	+
501	+	nt = fs->make<TNtuple>("ntEvent","","nHF:nHFlong:nHFshort");
502
503		}
504
#	Line 249 | Line 512 | double RecHitTreeProducer::getEt(const D
512		double et = energy*sin(pos.theta());
513		return et;
514		}
515	<
515	>
516		double RecHitTreeProducer::getEta(const DetId &id){
517		const GlobalPoint& pos=geo->getPosition(id);
518		double et = pos.eta();

Diff Legend

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