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.3 by nart, Wed Oct 20 15:01:11 2010 UTC vs.
Revision 1.19 by yilmaz, Wed Nov 30 21:26:08 2011 UTC

#	Line 12 | Line 12
12		*/
13		//
14		// Original Author:  Yetkin Yilmaz
15	+	// Modified: Frank Ma, Yen-Jie Lee
16		//         Created:  Tue Sep  7 11:38:19 EDT 2010
17		// $Id$
18		//
19		//
20
21	<
21	>	#define versionTag "v1"
22		// system include files
23		#include <memory>
24		#include <vector>
#	Line 57 | Line 58
58		#include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
59		#include "DataFormats/HcalDetId/interface/HcalDetId.h"
60		#include "DataFormats/EcalDetId/interface/EcalDetIdCollections.h"
61	+	#include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
62		#include "DataFormats/DetId/interface/DetId.h"
63
64	+	#include "DataFormats/VertexReco/interface/Vertex.h"
65	+	#include "DataFormats/VertexReco/interface/VertexFwd.h"
66	+
67
68		#include "TNtuple.h"
69
#	Line 67 | Line 72 | using namespace std;
72		#define MAXHITS 1000000
73
74		struct MyRecHit{
75	<
75	>	int depth[MAXHITS];
76		int n;
77
78		float e[MAXHITS];
79		float et[MAXHITS];
80		float eta[MAXHITS];
81		float phi[MAXHITS];
82	+	float perp[MAXHITS];
83	+	float emEt[MAXHITS];
84	+	float hadEt[MAXHITS];
85	+
86		bool isjet[MAXHITS];
87	+	float etVtx[MAXHITS];
88	+	float etaVtx[MAXHITS];
89	+	float phiVtx[MAXHITS];
90	+	float perpVtx[MAXHITS];
91	+	float emEtVtx[MAXHITS];
92	+	float hadEtVtx[MAXHITS];
93	+
94	+	float jtpt;
95	+	float jteta;
96	+	float jtphi;
97
98		};
99
100
101	+	struct MyBkg{
102	+	int n;
103	+	float rho[50];
104	+	float sigma[50];
105	+	};
106	+
107
108		//
109		// class declaration
#	Line 88 | Line 113 | class RecHitTreeProducer : public edm::E
113		public:
114		explicit RecHitTreeProducer(const edm::ParameterSet&);
115		~RecHitTreeProducer();
116	<	double       getEt(const DetId &id, double energy);
117	<	double       getEta(const DetId &id);
118	<	double       getPhi(const DetId &id);
116	>	double       getEt(const DetId &id, double energy, reco::Vertex::Point vtx = reco::Vertex::Point(0,0,0));
117	>	double       getEta(const DetId &id, reco::Vertex::Point vtx = reco::Vertex::Point(0,0,0));
118	>	double       getPhi(const DetId &id, reco::Vertex::Point vtx = reco::Vertex::Point(0,0,0));
119	>	double       getPerp(const DetId &id, reco::Vertex::Point vtx = reco::Vertex::Point(0,0,0));
120	>
121	>	math::XYZPoint getPosition(const DetId &id, reco::Vertex::Point vtx = reco::Vertex::Point(0,0,0));
122	>	double getEt(math::XYZPoint pos, double energy);
123	>
124	>	reco::Vertex::Point getVtx(const edm::Event& ev);
125	>
126
127
128		private:
#	Line 110 | Line 142 | class RecHitTreeProducer : public edm::E
142		edm::Handle<HFRecHitCollection> hfHits;
143		edm::Handle<HBHERecHitCollection> hbheHits;
144
145	+	edm::Handle<reco::BasicClusterCollection> bClusters;
146	+	edm::Handle<CaloTowerCollection> towers;
147	+	edm::Handle<reco::VertexCollection> vtxs;
148	+
149		typedef vector<EcalRecHit>::const_iterator EcalIterator;
150
151		edm::Handle<reco::CaloJetCollection> jets;
152	+
153	+	edm::Handle<std::vector<double> > rhos;
154	+	edm::Handle<std::vector<double> > sigmas;
155
156		MyRecHit hbheRecHit;
157		MyRecHit hfRecHit;
158		MyRecHit ebRecHit;
159		MyRecHit eeRecHit;
160	+	MyRecHit myBC;
161	+	MyRecHit myTowers;
162	+	MyBkg bkg;
163
164	+	TNtuple* nt;
165		TTree* hbheTree;
166		TTree* hfTree;
167		TTree* ebTree;
168		TTree* eeTree;
169	<	double cone;
169	>	TTree* bcTree;
170	>	TTree* towerTree;
171	>	TTree* bkgTree;
172
173	+	double cone;
174	+	double hfTowerThreshold_;
175	+	double hfLongThreshold_;
176	+	double hfShortThreshold_;
177	+	double hbheThreshold_;
178	+	double ebThreshold_;
179	+	double eeThreshold_;
180	+
181	+	double hbhePtMin_;
182	+	double hfPtMin_;
183	+	double ebPtMin_;
184	+	double eePtMin_;
185	+	double towerPtMin_;
186	+
187		edm::Service<TFileService> fs;
188		const CentralityBins * cbins_;
189		const CaloGeometry *geo;
#	Line 133 | Line 192 | class RecHitTreeProducer : public edm::E
192		edm::InputTag HcalRecHitHBHESrc_;
193		edm::InputTag EBSrc_;
194		edm::InputTag EESrc_;
195	+	edm::InputTag BCSrc_;
196	+	edm::InputTag TowerSrc_;
197	+	edm::InputTag VtxSrc_;
198	+
199		edm::InputTag JetSrc_;
200	<	bool excludeJets_;
200	>
201	>	edm::InputTag FastJetTag_;
202	>
203	>	bool useJets_;
204	>	bool doBasicClusters_;
205	>	bool doTowers_;
206	>	bool doEcal_;
207	>	bool doHcal_;
208	>	bool doHF_;
209	>	bool hasVtx_;
210	>	bool saveBothVtx_;
211	>
212	>	bool doFastJet_;
213	>
214	>	bool doEbyEonly_;
215
216		};
217
#	Line 155 | Line 232 | RecHitTreeProducer::RecHitTreeProducer(c
232		cone(0.5)
233		{
234		//now do what ever initialization is needed
158	–
235		EBSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("EBRecHitSrc",edm::InputTag("ecalRecHit","EcalRecHitsEB"));
236		EESrc_ = iConfig.getUntrackedParameter<edm::InputTag>("EERecHitSrc",edm::InputTag("ecalRecHit","EcalRecHitsEE"));
237		HcalRecHitHFSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHFRecHitSrc",edm::InputTag("hfreco"));
238		HcalRecHitHBHESrc_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHBHERecHitSrc",edm::InputTag("hbhereco"));
239	<	JetSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("JetSrc",edm::InputTag("iterativeCone5CaloJets"));
240	<	excludeJets_ = iConfig.getUntrackedParameter<bool>("excludeJets",false);
239	>	BCSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("BasicClusterSrc1",edm::InputTag("ecalRecHit","EcalRecHitsEB","RECO"));
240	>	TowerSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("towersSrc",edm::InputTag("towerMaker"));
241	>	VtxSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("vtxSrc",edm::InputTag("hiSelectedVertex"));
242	>	JetSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("JetSrc",edm::InputTag("iterativeConePu5CaloJets"));
243	>	useJets_ = iConfig.getUntrackedParameter<bool>("useJets",true);
244	>	doBasicClusters_ = iConfig.getUntrackedParameter<bool>("doBasicClusters",false);
245	>	doTowers_ = iConfig.getUntrackedParameter<bool>("doTowers",true);
246	>	doEcal_ = iConfig.getUntrackedParameter<bool>("doEcal",true);
247	>	doHcal_ = iConfig.getUntrackedParameter<bool>("doHcal",true);
248	>	doHF_ = iConfig.getUntrackedParameter<bool>("doHF",true);
249	>	hasVtx_ = iConfig.getUntrackedParameter<bool>("hasVtx",true);
250	>	saveBothVtx_ = iConfig.getUntrackedParameter<bool>("saveBothVtx",false);
251	>
252	>	doFastJet_ = iConfig.getUntrackedParameter<bool>("doFastJet",true);
253	>	FastJetTag_ = iConfig.getUntrackedParameter<edm::InputTag>("FastJetTag",edm::InputTag("kt4CaloJets"));
254	>	doEbyEonly_ = iConfig.getUntrackedParameter<bool>("doEbyEonly",false);
255	>	hfTowerThreshold_ = iConfig.getUntrackedParameter<double>("HFtowerMin",3.);
256	>	hfLongThreshold_ = iConfig.getUntrackedParameter<double>("HFlongMin",0.5);
257	>	hfShortThreshold_ = iConfig.getUntrackedParameter<double>("HFshortMin",0.85);
258	>	hbhePtMin_ = iConfig.getUntrackedParameter<double>("HBHETreePtMin",0);
259	>	hfPtMin_ = iConfig.getUntrackedParameter<double>("HFTreePtMin",0);
260	>	ebPtMin_ = iConfig.getUntrackedParameter<double>("EBTreePtMin",0);
261	>	eePtMin_ = iConfig.getUntrackedParameter<double>("EETreePtMin",0.);
262	>	towerPtMin_ = iConfig.getUntrackedParameter<double>("TowerTreePtMin",0.);
263	>
264	>
265		}
266
267
#	Line 187 | Line 287 | RecHitTreeProducer::analyze(const edm::E
287		hbheRecHit.n = 0;
288		ebRecHit.n = 0;
289		eeRecHit.n = 0;
290	+	myBC.n = 0;
291	+	myTowers.n = 0;
292	+	bkg.n = 0;
293	+
294	+	// get vertex
295	+	reco::Vertex::Point vtx(0,0,0);
296	+	if (hasVtx_) vtx = getVtx(ev);
297	+
298	+	if(doEcal_){
299	+	ev.getByLabel(EBSrc_,ebHits);
300	+	ev.getByLabel(EESrc_,eeHits);
301	+	}
302	+
303	+	if(doHcal_){
304	+	ev.getByLabel(HcalRecHitHBHESrc_,hbheHits);
305	+	}
306	+	if(doHF_){
307	+	ev.getByLabel(HcalRecHitHFSrc_,hfHits);
308	+	}
309
310	<	ev.getByLabel(EBSrc_,ebHits);
311	<	ev.getByLabel(EESrc_,eeHits);
310	>	if(useJets_) {
311	>	ev.getByLabel(JetSrc_,jets);
312	>	}
313
314	<	ev.getByLabel(HcalRecHitHFSrc_,hfHits);
315	<	ev.getByLabel(HcalRecHitHBHESrc_,hbheHits);
314	>	if(doBasicClusters_){
315	>	ev.getByLabel(BCSrc_,bClusters);
316	>	}
317
318	+	if(doTowers_){
319	+	ev.getByLabel(TowerSrc_,towers);
320	+	}
321
322	<	if(!excludeJets_) {
323	<	ev.getByLabel(JetSrc_,jets);
322	>	if(doFastJet_){
323	>	ev.getByLabel(edm::InputTag(FastJetTag_.label(),"rhos",FastJetTag_.process()),rhos);
324	>	ev.getByLabel(edm::InputTag(FastJetTag_.label(),"sigmas",FastJetTag_.process()),sigmas);
325	>	bkg.n = rhos->size();
326	>	for(unsigned int i = 0; i < rhos->size(); ++i){
327	>	bkg.rho[i] = (*rhos)[i];
328	>	bkg.sigma[i] = (*sigmas)[i];
329	>	}
330		}
331
332	<	if(0 && !cbins_) cbins_ = getCentralityBinsFromDB(iSetup);
332	>	if(0 && !cbins_) cbins_ = getCentralityBinsFromDB(iSetup);
333
334		if(!geo){
335		edm::ESHandle<CaloGeometry> pGeo;
336		iSetup.get<CaloGeometryRecord>().get(pGeo);
337		geo = pGeo.product();
338		}
339	<
339	>
340	>	int nHFlongPlus = 0;
341	>	int nHFshortPlus = 0;
342	>	int nHFtowerPlus = 0;
343	>	int nHFlongMinus = 0;
344	>	int nHFshortMinus = 0;
345	>	int nHFtowerMinus = 0;
346	>
347	>	if(doHF_){
348		for(unsigned int i = 0; i < hfHits->size(); ++i){
349		const HFRecHit & hit= (*hfHits)[i];
350		hfRecHit.e[hfRecHit.n] = hit.energy();
351	<	hfRecHit.et[hfRecHit.n] = getEt(hit.id(),hit.energy());
352	<	hfRecHit.eta[hfRecHit.n] = getEta(hit.id());
353	<	hfRecHit.phi[hfRecHit.n] = getPhi(hit.id());
351	>	math::XYZPoint pos = getPosition(hit.id(),vtx);
352	>
353	>	if(!saveBothVtx_){
354	>	hfRecHit.et[hfRecHit.n] = getEt(pos,hit.energy());
355	>	hfRecHit.eta[hfRecHit.n] = pos.eta();
356	>	hfRecHit.phi[hfRecHit.n] = pos.phi();
357	>	hfRecHit.perp[hfRecHit.n] = pos.rho();
358	>	}else{
359	>	hfRecHit.et[hfRecHit.n] = getEt(hit.id(),hit.energy());
360	>	hfRecHit.eta[hfRecHit.n] = getEta(hit.id());
361	>	hfRecHit.phi[hfRecHit.n] = getPhi(hit.id());
362	>	hfRecHit.perp[hfRecHit.n] = getPerp(hit.id());
363	>
364	>	hfRecHit.etVtx[hfRecHit.n] = getEt(pos,hit.energy());
365	>	hfRecHit.etaVtx[hfRecHit.n] = pos.eta();
366	>	hfRecHit.phiVtx[hfRecHit.n] = pos.phi();
367	>	hfRecHit.perpVtx[hfRecHit.n] = pos.rho();
368	>
369	>	}
370	>
371		hfRecHit.isjet[hfRecHit.n] = false;
372	<	if(!excludeJets_){
372	>	hfRecHit.depth[hfRecHit.n] = hit.id().depth();
373	>
374	>	if(hit.id().ieta() > 0){
375	>	if(hit.energy() > hfShortThreshold_ && hit.id().depth() != 1) nHFshortPlus++;
376	>	if(hit.energy() > hfLongThreshold_ && hit.id().depth() == 1) nHFlongPlus++;
377	>	}else{
378	>	if(hit.energy() > hfShortThreshold_ && hit.id().depth() != 1) nHFshortMinus++;
379	>	if(hit.energy() > hfLongThreshold_ && hit.id().depth() == 1) nHFlongMinus++;
380	>	}
381	>
382	>	if(useJets_){
383		for(unsigned int j = 0 ; j < jets->size(); ++j){
384		const reco::Jet& jet = (*jets)[j];
385		double dr = reco::deltaR(hfRecHit.eta[hfRecHit.n],hfRecHit.phi[hfRecHit.n],jet.eta(),jet.phi());
386		if(dr < cone){ hfRecHit.isjet[hfRecHit.n] = true; }
387		}
388		}
389	<	hfRecHit.n++;
389	>	if (hfRecHit.et[hfRecHit.n]>=hfPtMin_) hfRecHit.n++;
390		}
391	<
391	>
392	>	if(doHcal_ && !doEbyEonly_){
393		for(unsigned int i = 0; i < hbheHits->size(); ++i){
394		const HBHERecHit & hit= (*hbheHits)[i];
395	+	if (getEt(hit.id(),hit.energy())<hbhePtMin_) continue;
396	+
397		hbheRecHit.e[hbheRecHit.n] = hit.energy();
398	<	hbheRecHit.et[hbheRecHit.n] = getEt(hit.id(),hit.energy());
399	<	hbheRecHit.eta[hbheRecHit.n] = getEta(hit.id());
400	<	hbheRecHit.phi[hbheRecHit.n] = getPhi(hit.id());
401	<	hbheRecHit.isjet[hbheRecHit.n] = false;
402	<	if(!excludeJets_){
398	>	math::XYZPoint pos = getPosition(hit.id(),vtx);
399	>
400	>	if(!saveBothVtx_){
401	>	hbheRecHit.et[hbheRecHit.n] = getEt(pos,hit.energy());
402	>	hbheRecHit.eta[hbheRecHit.n] = pos.eta();
403	>	hbheRecHit.phi[hbheRecHit.n] = pos.phi();
404	>	hbheRecHit.perp[hbheRecHit.n] = pos.rho();
405	>	}else{
406	>	hbheRecHit.et[hbheRecHit.n] = getEt(hit.id(),hit.energy());
407	>	hbheRecHit.eta[hbheRecHit.n] = getEta(hit.id());
408	>	hbheRecHit.phi[hbheRecHit.n] = getPhi(hit.id());
409	>	hbheRecHit.perp[hbheRecHit.n] = getPerp(hit.id());
410	>
411	>	hbheRecHit.etVtx[hbheRecHit.n] = getEt(pos,hit.energy());
412	>	hbheRecHit.etaVtx[hbheRecHit.n] = pos.eta();
413	>	hbheRecHit.phiVtx[hbheRecHit.n] = pos.phi();
414	>	hbheRecHit.perpVtx[hbheRecHit.n] = pos.rho();
415	>
416	>	}
417	>
418	>	hbheRecHit.isjet[hbheRecHit.n] = false;
419	>	hbheRecHit.depth[hbheRecHit.n] = hit.id().depth();
420	>
421	>	if(useJets_){
422		for(unsigned int j = 0 ; j < jets->size(); ++j){
423		const reco::Jet& jet = (*jets)[j];
424		double dr = reco::deltaR(hbheRecHit.eta[hbheRecHit.n],hbheRecHit.phi[hbheRecHit.n],jet.eta(),jet.phi());
#	Line 240 | Line 427 | RecHitTreeProducer::analyze(const edm::E
427		}
428		hbheRecHit.n++;
429		}
430	<
430	>	}
431	>	}
432	>	if(doEcal_ && !doEbyEonly_){
433		for(unsigned int i = 0; i < ebHits->size(); ++i){
434		const EcalRecHit & hit= (*ebHits)[i];
435	+	if (getEt(hit.id(),hit.energy())<ebPtMin_) continue;
436	+
437		ebRecHit.e[ebRecHit.n] = hit.energy();
438	<	ebRecHit.et[ebRecHit.n] = getEt(hit.id(),hit.energy());
439	<	ebRecHit.eta[ebRecHit.n] = getEta(hit.id());
440	<	ebRecHit.phi[ebRecHit.n] = getPhi(hit.id());
438	>	math::XYZPoint pos = getPosition(hit.id(),vtx);
439	>
440	>	if(!saveBothVtx_){
441	>	ebRecHit.et[ebRecHit.n] = getEt(pos,hit.energy());
442	>	ebRecHit.eta[ebRecHit.n] = pos.eta();
443	>	ebRecHit.phi[ebRecHit.n] = pos.phi();
444	>	ebRecHit.perp[ebRecHit.n] = pos.rho();
445	>	}else{
446	>	ebRecHit.et[ebRecHit.n] = getEt(hit.id(),hit.energy());
447	>	ebRecHit.eta[ebRecHit.n] = getEta(hit.id());
448	>	ebRecHit.phi[ebRecHit.n] = getPhi(hit.id());
449	>	ebRecHit.perp[ebRecHit.n] = getPerp(hit.id());
450	>	ebRecHit.etVtx[ebRecHit.n] = getEt(pos,hit.energy());
451	>	ebRecHit.etaVtx[ebRecHit.n] = pos.eta();
452	>	ebRecHit.phiVtx[ebRecHit.n] = pos.phi();
453	>	ebRecHit.perpVtx[ebRecHit.n] = pos.rho();
454	>
455	>	}
456	>
457		ebRecHit.isjet[ebRecHit.n] = false;
458	<	if(!excludeJets_){
458	>	if(useJets_){
459		for(unsigned int j = 0 ; j < jets->size(); ++j){
460		const reco::Jet& jet = (*jets)[j];
461		double dr = reco::deltaR(ebRecHit.eta[ebRecHit.n],ebRecHit.phi[ebRecHit.n],jet.eta(),jet.phi());
#	Line 260 | Line 467 | RecHitTreeProducer::analyze(const edm::E
467
468		for(unsigned int i = 0; i < eeHits->size(); ++i){
469		const EcalRecHit & hit= (*eeHits)[i];
470	+	if (getEt(hit.id(),hit.energy())<eePtMin_) continue;
471	+
472		eeRecHit.e[eeRecHit.n] = hit.energy();
473	<	eeRecHit.et[eeRecHit.n] = getEt(hit.id(),hit.energy());
474	<	eeRecHit.eta[eeRecHit.n] = getEta(hit.id());
475	<	eeRecHit.phi[eeRecHit.n] = getPhi(hit.id());
473	>	math::XYZPoint pos = getPosition(hit.id(),vtx);
474	>
475	>	if(!saveBothVtx_){
476	>	eeRecHit.et[eeRecHit.n] = getEt(pos,hit.energy());
477	>	eeRecHit.eta[eeRecHit.n] = pos.eta();
478	>	eeRecHit.phi[eeRecHit.n] = pos.phi();
479	>	eeRecHit.perp[eeRecHit.n] = pos.rho();
480	>	}else{
481	>	eeRecHit.et[eeRecHit.n] = getEt(hit.id(),hit.energy());
482	>	eeRecHit.eta[eeRecHit.n] = getEta(hit.id());
483	>	eeRecHit.phi[eeRecHit.n] = getPhi(hit.id());
484	>	eeRecHit.perp[eeRecHit.n] = getPerp(hit.id());
485	>	eeRecHit.etVtx[eeRecHit.n] = getEt(pos,hit.energy());
486	>	eeRecHit.etaVtx[eeRecHit.n] = pos.eta();
487	>	eeRecHit.phiVtx[eeRecHit.n] = pos.phi();
488	>	eeRecHit.perpVtx[eeRecHit.n] = pos.rho();
489	>
490	>	}
491	>
492		eeRecHit.isjet[eeRecHit.n] = false;
493	<	if(!excludeJets_){
493	>
494	>	if(useJets_){
495		for(unsigned int j = 0 ; j < jets->size(); ++j){
496		const reco::Jet& jet = (*jets)[j];
497		double dr = reco::deltaR(eeRecHit.eta[eeRecHit.n],eeRecHit.phi[eeRecHit.n],jet.eta(),jet.phi());
#	Line 274 | Line 500 | RecHitTreeProducer::analyze(const edm::E
500		}
501		eeRecHit.n++;
502		}
503	<
504	<	eeTree->Fill();
505	<	ebTree->Fill();
506	<
507	<	hbheTree->Fill();
508	<	hfTree->Fill();
503	>	}
504	>
505	>	if(doTowers_){
506	>
507	>	for(unsigned int i = 0; i < towers->size(); ++i){
508	>	const CaloTower & hit= (*towers)[i];
509	>	if (getEt(hit.id(),hit.energy())<towerPtMin_) continue;
510	>
511	>	myTowers.et[myTowers.n] = hit.p4(vtx).Et();
512	>	myTowers.eta[myTowers.n] = hit.p4(vtx).Eta();
513	>	myTowers.phi[myTowers.n] = hit.p4(vtx).Phi();
514	>	myTowers.emEt[myTowers.n] = hit.emEt(vtx);
515	>	myTowers.hadEt[myTowers.n] = hit.hadEt(vtx);
516	>
517	>	if (saveBothVtx_) {
518	>	myTowers.e[myTowers.n] = hit.energy();
519	>	myTowers.et[myTowers.n] = getEt(hit.id(),hit.energy());
520	>	myTowers.eta[myTowers.n] = getEta(hit.id());
521	>	myTowers.phi[myTowers.n] = getPhi(hit.id());
522	>	myTowers.isjet[myTowers.n] = false;
523	>	myTowers.etVtx[myTowers.n] = hit.p4(vtx).Et();
524	>	myTowers.etaVtx[myTowers.n] = hit.p4(vtx).Eta();
525	>	myTowers.emEtVtx[myTowers.n] = hit.emEt(vtx);
526	>	myTowers.hadEtVtx[myTowers.n] = hit.hadEt(vtx);
527	>	}
528	>
529	>	myTowers.isjet[myTowers.n] = false;
530	>
531	>	if(hit.ieta() > 29 && hit.energy() > hfTowerThreshold_) nHFtowerPlus++;
532	>	if(hit.ieta() < -29 && hit.energy() > hfTowerThreshold_) nHFtowerMinus++;
533	>
534	>	if(useJets_){
535	>	for(unsigned int j = 0 ; j < jets->size(); ++j){
536	>	const reco::Jet& jet = (*jets)[j];
537	>	double dr = reco::deltaR(myTowers.eta[myTowers.n],myTowers.phi[myTowers.n],jet.eta(),jet.phi());
538	>	if(dr < cone){ myTowers.isjet[myTowers.n] = true; }
539	>	}
540	>	}
541	>	myTowers.n++;
542	>	}
543	>
544	>	}
545	>
546	>	if(doBasicClusters_ && !doEbyEonly_){
547	>	for(unsigned int j = 0 ; j < jets->size(); ++j){
548	>	const reco::Jet& jet = (*jets)[j];
549	>	myBC.n = 0;
550	>	myBC.jtpt = jet.pt();
551	>	myBC.jteta = jet.eta();
552	>	myBC.jtphi = jet.phi();
553	>
554	>	for(unsigned int i = 0; i < bClusters->size(); ++i){
555	>	const reco::BasicCluster & bc= (*bClusters)[i];
556	>	double dr = reco::deltaR(bc.eta(),bc.phi(),jet.eta(),jet.phi());
557	>	if(dr < cone){
558	>	myBC.e[myBC.n] = bc.energy();
559	>	myBC.et[myBC.n] = bc.energy()*sin(bc.position().theta());
560	>	myBC.eta[myBC.n] = bc.eta();
561	>	myBC.phi[myBC.n] = bc.phi();
562	>	myBC.n++;
563	>	}
564	>	}
565	>	bcTree->Fill();
566	>	}
567	>	}
568	>
569	>	if(!doEbyEonly_){
570	>	towerTree->Fill();
571	>
572	>	eeTree->Fill();
573	>	ebTree->Fill();
574	>
575	>	hbheTree->Fill();
576	>	hfTree->Fill();
577	>
578	>	if (doFastJet_) {
579	>	bkgTree->Fill();
580	>	}
581	>	}
582	>
583	>	nt->Fill(nHFtowerPlus,nHFtowerMinus,nHFlongPlus,nHFlongMinus,nHFshortPlus,nHFshortMinus);
584
585		}
586
#	Line 289 | Line 590 | void
590		RecHitTreeProducer::beginJob()
591		{
592
593	<	hbheTree = fs->make<TTree>("hbhe","");
593	>	hbheTree = fs->make<TTree>("hbhe",versionTag);
594		hbheTree->Branch("n",&hbheRecHit.n,"n/I");
595		hbheTree->Branch("e",hbheRecHit.e,"e[n]/F");
596		hbheTree->Branch("et",hbheRecHit.et,"et[n]/F");
597		hbheTree->Branch("eta",hbheRecHit.eta,"eta[n]/F");
598		hbheTree->Branch("phi",hbheRecHit.phi,"phi[n]/F");
599	<	hbheTree->Branch("isjet",hbheRecHit.isjet,"isjet[n]/I");
600	<
601	<	hfTree = fs->make<TTree>("hf","");
599	>	hbheTree->Branch("perp",hbheRecHit.perp,"perp[n]/F");
600	>
601	>	hbheTree->Branch("isjet",hbheRecHit.isjet,"isjet[n]/O");
602	>	hbheTree->Branch("depth",hfRecHit.depth,"depth[n]/I");
603	>
604	>	hfTree = fs->make<TTree>("hf",versionTag);
605		hfTree->Branch("n",&hfRecHit.n,"n/I");
606		hfTree->Branch("e",hfRecHit.e,"e[n]/F");
607		hfTree->Branch("et",hfRecHit.et,"et[n]/F");
608		hfTree->Branch("eta",hfRecHit.eta,"eta[n]/F");
609		hfTree->Branch("phi",hfRecHit.phi,"phi[n]/F");
610	<	hfTree->Branch("isjet",hfRecHit.isjet,"isjet[n]/I");
610	>	hfTree->Branch("perp",hfRecHit.perp,"perp[n]/F");
611	>	hfTree->Branch("depth",hfRecHit.depth,"depth[n]/I");
612	>	hfTree->Branch("isjet",hfRecHit.isjet,"isjet[n]/O");
613
614	<	eeTree = fs->make<TTree>("ee","");
614	>	eeTree = fs->make<TTree>("ee",versionTag);
615		eeTree->Branch("n",&eeRecHit.n,"n/I");
616		eeTree->Branch("e",eeRecHit.e,"e[n]/F");
617		eeTree->Branch("et",eeRecHit.et,"et[n]/F");
618		eeTree->Branch("eta",eeRecHit.eta,"eta[n]/F");
619		eeTree->Branch("phi",eeRecHit.phi,"phi[n]/F");
620	<	eeTree->Branch("isjet",eeRecHit.isjet,"isjet[n]/I");
620	>	eeTree->Branch("perp",eeRecHit.perp,"perp[n]/F");
621	>
622	>	eeTree->Branch("isjet",eeRecHit.isjet,"isjet[n]/O");
623
624	<	ebTree = fs->make<TTree>("eb","");
624	>	ebTree = fs->make<TTree>("eb",versionTag);
625		ebTree->Branch("n",&ebRecHit.n,"n/I");
626		ebTree->Branch("e",ebRecHit.e,"e[n]/F");
627		ebTree->Branch("et",ebRecHit.et,"et[n]/F");
628		ebTree->Branch("eta",ebRecHit.eta,"eta[n]/F");
629		ebTree->Branch("phi",ebRecHit.phi,"phi[n]/F");
630	<	ebTree->Branch("isjet",ebRecHit.isjet,"isjet[n]/I");
630	>	ebTree->Branch("perp",ebRecHit.perp,"perp[n]/F");
631	>
632	>	ebTree->Branch("isjet",ebRecHit.isjet,"isjet[n]/O");
633	>
634	>	towerTree = fs->make<TTree>("tower",versionTag);
635	>	towerTree->Branch("n",&myTowers.n,"n/I");
636	>	towerTree->Branch("e",myTowers.e,"e[n]/F");
637	>	towerTree->Branch("et",myTowers.et,"et[n]/F");
638	>	towerTree->Branch("eta",myTowers.eta,"eta[n]/F");
639	>	towerTree->Branch("phi",myTowers.phi,"phi[n]/F");
640	>	towerTree->Branch("isjet",myTowers.isjet,"isjet[n]/O");
641	>	towerTree->Branch("emEt",myTowers.emEt,"emEt[n]/F");
642	>	towerTree->Branch("hadEt",myTowers.hadEt,"hadEt[n]/F");
643	>
644	>	if (saveBothVtx_) {
645	>	towerTree->Branch("etVtx",myTowers.etVtx,"etvtx[n]/F");
646	>	towerTree->Branch("etaVtx",myTowers.etaVtx,"etavtx[n]/F");
647	>	towerTree->Branch("emEtVtx",myTowers.emEtVtx,"emEtVtx[n]/F");
648	>	towerTree->Branch("hadEtVtx",myTowers.hadEtVtx,"hadEtVtx[n]/F");
649	>	}
650	>
651	>
652	>	if(doBasicClusters_){
653	>	bcTree = fs->make<TTree>("bc",versionTag);
654	>	bcTree->Branch("n",&myBC.n,"n/I");
655	>	bcTree->Branch("e",myBC.e,"e[n]/F");
656	>	bcTree->Branch("et",myBC.et,"et[n]/F");
657	>	bcTree->Branch("eta",myBC.eta,"eta[n]/F");
658	>	bcTree->Branch("phi",myBC.phi,"phi[n]/F");
659	>	bcTree->Branch("jtpt",&myBC.jtpt,"jtpt/F");
660	>	bcTree->Branch("jteta",&myBC.jteta,"jteta/F");
661	>	bcTree->Branch("jtphi",&myBC.jtphi,"jtphi/F");
662	>	//     bcTree->Branch("isjet",bcRecHit.isjet,"isjet[n]/O");
663	>	}
664	>
665	>	if(doFastJet_){
666	>	bkgTree = fs->make<TTree>("bkg",versionTag);
667	>	bkgTree->Branch("n",&bkg.n,"n/I");
668	>	bkgTree->Branch("rho",bkg.rho,"rho[n]/F");
669	>	bkgTree->Branch("sigma",bkg.sigma,"sigma[n]/F");
670	>	}
671	>
672	>	nt = fs->make<TNtuple>("ntEvent","","nHFplus:nHFminus:nHFlongPlus:nHFlongMinus:nHFshortPlus:nHFshortMinus");
673
674		}
675
#	Line 328 | Line 678 | void
678		RecHitTreeProducer::endJob() {
679		}
680
681	<	double RecHitTreeProducer::getEt(const DetId &id, double energy){
681	>	math::XYZPoint RecHitTreeProducer::getPosition(const DetId &id, reco::Vertex::Point vtx){
682	>	const GlobalPoint& pos=geo->getPosition(id);
683	>	math::XYZPoint posV(pos.x() - vtx.x(),pos.y() - vtx.y(),pos.z() - vtx.z());
684	>	return posV;
685	>	}
686	>
687	>	double RecHitTreeProducer::getEt(math::XYZPoint pos, double energy){
688	>	double et = energy*sin(pos.theta());
689	>	return et;
690	>	}
691	>
692	>	double RecHitTreeProducer::getEt(const DetId &id, double energy, reco::Vertex::Point vtx){
693		const GlobalPoint& pos=geo->getPosition(id);
694		double et = energy*sin(pos.theta());
695		return et;
696		}
697
698	<	double RecHitTreeProducer::getEta(const DetId &id){
698	>	double RecHitTreeProducer::getEta(const DetId &id, reco::Vertex::Point vtx){
699		const GlobalPoint& pos=geo->getPosition(id);
700		double et = pos.eta();
701		return et;
702		}
703
704	<	double RecHitTreeProducer::getPhi(const DetId &id){
704	>	double RecHitTreeProducer::getPhi(const DetId &id, reco::Vertex::Point vtx){
705		const GlobalPoint& pos=geo->getPosition(id);
706		double et = pos.phi();
707		return et;
708		}
709
710	+	double RecHitTreeProducer::getPerp(const DetId &id, reco::Vertex::Point vtx){
711	+	const GlobalPoint& pos=geo->getPosition(id);
712	+	double et = pos.perp();
713	+	return et;
714	+	}
715
716	+	reco::Vertex::Point RecHitTreeProducer::getVtx(const edm::Event& ev)
717	+	{
718	+	ev.getByLabel(VtxSrc_,vtxs);
719	+	int greatestvtx = 0;
720	+	int nVertex = vtxs->size();
721	+
722	+	for (unsigned int i = 0 ; i< vtxs->size(); ++i){
723	+	unsigned int daughter = (*vtxs)[i].tracksSize();
724	+	if( daughter > (*vtxs)[greatestvtx].tracksSize()) greatestvtx = i;
725	+	//cout <<"Vertex: "<< (*vtxs)[i].position().z()<<" "<<daughter<<endl;
726	+	}
727	+
728	+	if(nVertex<=0){
729	+	return reco::Vertex::Point(0,0,0);
730	+	}
731	+	return (*vtxs)[greatestvtx].position();
732	+	}
733
734		//define this as a plug-in
735		DEFINE_FWK_MODULE(RecHitTreeProducer);

Diff Legend

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