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.22 by yilmaz, Wed Jan 16 17:22:56 2013 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
70		using namespace std;
71
72	<	#define MAXHITS 1000000
72	>	#define MAXHITS 100000
73
74		struct MyRecHit{
75	<
75	>	int depth[MAXHITS];
76		int n;
77
78	+	int ieta[MAXHITS];
79	+	int iphi[MAXHITS];
80	+
81		float e[MAXHITS];
82		float et[MAXHITS];
83		float eta[MAXHITS];
84		float phi[MAXHITS];
85	+	float perp[MAXHITS];
86	+	float emEt[MAXHITS];
87	+	float hadEt[MAXHITS];
88	+
89		bool isjet[MAXHITS];
90	+	float etVtx[MAXHITS];
91	+	float etaVtx[MAXHITS];
92	+	float phiVtx[MAXHITS];
93	+	float perpVtx[MAXHITS];
94	+	float emEtVtx[MAXHITS];
95	+	float hadEtVtx[MAXHITS];
96	+
97	+	float jtpt;
98	+	float jteta;
99	+	float jtphi;
100
101		};
102
103
104	+	struct MyBkg{
105	+	int n;
106	+	float rho[50];
107	+	float sigma[50];
108	+	};
109	+
110
111		//
112		// class declaration
#	Line 88 | Line 116 | class RecHitTreeProducer : public edm::E
116		public:
117		explicit RecHitTreeProducer(const edm::ParameterSet&);
118		~RecHitTreeProducer();
119	<	double       getEt(const DetId &id, double energy);
120	<	double       getEta(const DetId &id);
121	<	double       getPhi(const DetId &id);
119	>	double       getEt(const DetId &id, double energy, reco::Vertex::Point vtx = reco::Vertex::Point(0,0,0));
120	>	double       getEta(const DetId &id, reco::Vertex::Point vtx = reco::Vertex::Point(0,0,0));
121	>	double       getPhi(const DetId &id, reco::Vertex::Point vtx = reco::Vertex::Point(0,0,0));
122	>	double       getPerp(const DetId &id, reco::Vertex::Point vtx = reco::Vertex::Point(0,0,0));
123	>
124	>	math::XYZPoint getPosition(const DetId &id, reco::Vertex::Point vtx = reco::Vertex::Point(0,0,0));
125	>	double getEt(math::XYZPoint pos, double energy);
126	>
127	>	reco::Vertex::Point getVtx(const edm::Event& ev);
128	>
129
130
131		private:
#	Line 110 | Line 145 | class RecHitTreeProducer : public edm::E
145		edm::Handle<HFRecHitCollection> hfHits;
146		edm::Handle<HBHERecHitCollection> hbheHits;
147
148	+	edm::Handle<reco::BasicClusterCollection> bClusters;
149	+	edm::Handle<CaloTowerCollection> towers;
150	+	edm::Handle<reco::VertexCollection> vtxs;
151	+
152		typedef vector<EcalRecHit>::const_iterator EcalIterator;
153
154		edm::Handle<reco::CaloJetCollection> jets;
155	+
156	+	edm::Handle<std::vector<double> > rhos;
157	+	edm::Handle<std::vector<double> > sigmas;
158
159		MyRecHit hbheRecHit;
160		MyRecHit hfRecHit;
161		MyRecHit ebRecHit;
162		MyRecHit eeRecHit;
163	+	MyRecHit myBC;
164	+	MyRecHit myTowers;
165	+	MyRecHit castorRecHit;
166
167	+	MyBkg bkg;
168	+
169	+	TNtuple* nt;
170		TTree* hbheTree;
171		TTree* hfTree;
172		TTree* ebTree;
173		TTree* eeTree;
174	<	double cone;
174	>	TTree* bcTree;
175	>	TTree* towerTree;
176	>	TTree* bkgTree;
177	>	TTree* castorTree;
178
179	+	double cone;
180	+	double hfTowerThreshold_;
181	+	double hfLongThreshold_;
182	+	double hfShortThreshold_;
183	+	double hbheThreshold_;
184	+	double ebThreshold_;
185	+	double eeThreshold_;
186	+
187	+	double hbhePtMin_;
188	+	double hfPtMin_;
189	+	double ebPtMin_;
190	+	double eePtMin_;
191	+	double towerPtMin_;
192	+
193		edm::Service<TFileService> fs;
194		const CentralityBins * cbins_;
195		const CaloGeometry *geo;
#	Line 133 | Line 198 | class RecHitTreeProducer : public edm::E
198		edm::InputTag HcalRecHitHBHESrc_;
199		edm::InputTag EBSrc_;
200		edm::InputTag EESrc_;
201	+	edm::InputTag BCSrc_;
202	+	edm::InputTag TowerSrc_;
203	+	edm::InputTag VtxSrc_;
204	+
205		edm::InputTag JetSrc_;
206	<	bool excludeJets_;
206	>
207	>	edm::InputTag FastJetTag_;
208	>
209	>	bool useJets_;
210	>	bool doBasicClusters_;
211	>	bool doTowers_;
212	>	bool doEcal_;
213	>	bool doHcal_;
214	>	bool doHF_;
215	>	bool doCastor_;
216	>
217	>	bool hasVtx_;
218	>	bool saveBothVtx_;
219	>
220	>	bool doFastJet_;
221	>
222	>	bool doEbyEonly_;
223
224		};
225
#	Line 155 | Line 240 | RecHitTreeProducer::RecHitTreeProducer(c
240		cone(0.5)
241		{
242		//now do what ever initialization is needed
158	–
243		EBSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("EBRecHitSrc",edm::InputTag("ecalRecHit","EcalRecHitsEB"));
244		EESrc_ = iConfig.getUntrackedParameter<edm::InputTag>("EERecHitSrc",edm::InputTag("ecalRecHit","EcalRecHitsEE"));
245		HcalRecHitHFSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHFRecHitSrc",edm::InputTag("hfreco"));
246		HcalRecHitHBHESrc_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHBHERecHitSrc",edm::InputTag("hbhereco"));
247	<	JetSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("JetSrc",edm::InputTag("iterativeCone5CaloJets"));
248	<	excludeJets_ = iConfig.getUntrackedParameter<bool>("excludeJets",false);
247	>	BCSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("BasicClusterSrc1",edm::InputTag("ecalRecHit","EcalRecHitsEB","RECO"));
248	>	TowerSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("towersSrc",edm::InputTag("towerMaker"));
249	>	VtxSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("vtxSrc",edm::InputTag("hiSelectedVertex"));
250	>	JetSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("JetSrc",edm::InputTag("iterativeConePu5CaloJets"));
251	>	useJets_ = iConfig.getUntrackedParameter<bool>("useJets",true);
252	>	doBasicClusters_ = iConfig.getUntrackedParameter<bool>("doBasicClusters",false);
253	>	doTowers_ = iConfig.getUntrackedParameter<bool>("doTowers",true);
254	>	doEcal_ = iConfig.getUntrackedParameter<bool>("doEcal",true);
255	>	doHcal_ = iConfig.getUntrackedParameter<bool>("doHcal",true);
256	>	doHF_ = iConfig.getUntrackedParameter<bool>("doHF",true);
257	>	doCastor_ = iConfig.getUntrackedParameter<bool>("doCASTOR",true);
258	>
259	>	hasVtx_ = iConfig.getUntrackedParameter<bool>("hasVtx",true);
260	>	saveBothVtx_ = iConfig.getUntrackedParameter<bool>("saveBothVtx",false);
261	>
262	>	doFastJet_ = iConfig.getUntrackedParameter<bool>("doFastJet",true);
263	>	FastJetTag_ = iConfig.getUntrackedParameter<edm::InputTag>("FastJetTag",edm::InputTag("kt4CaloJets"));
264	>	doEbyEonly_ = iConfig.getUntrackedParameter<bool>("doEbyEonly",false);
265	>	hfTowerThreshold_ = iConfig.getUntrackedParameter<double>("HFtowerMin",3.);
266	>	hfLongThreshold_ = iConfig.getUntrackedParameter<double>("HFlongMin",0.5);
267	>	hfShortThreshold_ = iConfig.getUntrackedParameter<double>("HFshortMin",0.85);
268	>	hbhePtMin_ = iConfig.getUntrackedParameter<double>("HBHETreePtMin",0);
269	>	hfPtMin_ = iConfig.getUntrackedParameter<double>("HFTreePtMin",0);
270	>	ebPtMin_ = iConfig.getUntrackedParameter<double>("EBTreePtMin",0);
271	>	eePtMin_ = iConfig.getUntrackedParameter<double>("EETreePtMin",0.);
272	>	towerPtMin_ = iConfig.getUntrackedParameter<double>("TowerTreePtMin",0.);
273	>
274	>
275		}
276
277
#	Line 187 | Line 297 | RecHitTreeProducer::analyze(const edm::E
297		hbheRecHit.n = 0;
298		ebRecHit.n = 0;
299		eeRecHit.n = 0;
300	+	myBC.n = 0;
301	+	myTowers.n = 0;
302	+	bkg.n = 0;
303	+
304	+	// get vertex
305	+	reco::Vertex::Point vtx(0,0,0);
306	+	if (hasVtx_) vtx = getVtx(ev);
307	+
308	+	if(doEcal_){
309	+	ev.getByLabel(EBSrc_,ebHits);
310	+	ev.getByLabel(EESrc_,eeHits);
311	+	}
312	+
313	+	if(doHcal_){
314	+	ev.getByLabel(HcalRecHitHBHESrc_,hbheHits);
315	+	}
316	+	if(doHF_){
317	+	ev.getByLabel(HcalRecHitHFSrc_,hfHits);
318	+	}
319
320	<	ev.getByLabel(EBSrc_,ebHits);
321	<	ev.getByLabel(EESrc_,eeHits);
320	>	if(useJets_) {
321	>	ev.getByLabel(JetSrc_,jets);
322	>	}
323
324	<	ev.getByLabel(HcalRecHitHFSrc_,hfHits);
325	<	ev.getByLabel(HcalRecHitHBHESrc_,hbheHits);
324	>	if(doBasicClusters_){
325	>	ev.getByLabel(BCSrc_,bClusters);
326	>	}
327
328	+	if(doTowers_){
329	+	ev.getByLabel(TowerSrc_,towers);
330	+	}
331
332	<	if(!excludeJets_) {
333	<	ev.getByLabel(JetSrc_,jets);
332	>	if(doFastJet_){
333	>	ev.getByLabel(edm::InputTag(FastJetTag_.label(),"rhos",FastJetTag_.process()),rhos);
334	>	ev.getByLabel(edm::InputTag(FastJetTag_.label(),"sigmas",FastJetTag_.process()),sigmas);
335	>	bkg.n = rhos->size();
336	>	for(unsigned int i = 0; i < rhos->size(); ++i){
337	>	bkg.rho[i] = (*rhos)[i];
338	>	bkg.sigma[i] = (*sigmas)[i];
339	>	}
340		}
341
342	<	if(0 && !cbins_) cbins_ = getCentralityBinsFromDB(iSetup);
342	>	if(0 && !cbins_) cbins_ = getCentralityBinsFromDB(iSetup);
343
344		if(!geo){
345		edm::ESHandle<CaloGeometry> pGeo;
346		iSetup.get<CaloGeometryRecord>().get(pGeo);
347		geo = pGeo.product();
348		}
349	<
349	>
350	>	int nHFlongPlus = 0;
351	>	int nHFshortPlus = 0;
352	>	int nHFtowerPlus = 0;
353	>	int nHFlongMinus = 0;
354	>	int nHFshortMinus = 0;
355	>	int nHFtowerMinus = 0;
356	>
357	>	if(doHF_){
358		for(unsigned int i = 0; i < hfHits->size(); ++i){
359		const HFRecHit & hit= (*hfHits)[i];
360		hfRecHit.e[hfRecHit.n] = hit.energy();
361	<	hfRecHit.et[hfRecHit.n] = getEt(hit.id(),hit.energy());
362	<	hfRecHit.eta[hfRecHit.n] = getEta(hit.id());
363	<	hfRecHit.phi[hfRecHit.n] = getPhi(hit.id());
361	>	math::XYZPoint pos = getPosition(hit.id(),vtx);
362	>
363	>	if(!saveBothVtx_){
364	>	hfRecHit.et[hfRecHit.n] = getEt(pos,hit.energy());
365	>	hfRecHit.eta[hfRecHit.n] = pos.eta();
366	>	hfRecHit.phi[hfRecHit.n] = pos.phi();
367	>	hfRecHit.perp[hfRecHit.n] = pos.rho();
368	>	}else{
369	>	hfRecHit.et[hfRecHit.n] = getEt(hit.id(),hit.energy());
370	>	hfRecHit.eta[hfRecHit.n] = getEta(hit.id());
371	>	hfRecHit.phi[hfRecHit.n] = getPhi(hit.id());
372	>	hfRecHit.perp[hfRecHit.n] = getPerp(hit.id());
373	>
374	>	hfRecHit.etVtx[hfRecHit.n] = getEt(pos,hit.energy());
375	>	hfRecHit.etaVtx[hfRecHit.n] = pos.eta();
376	>	hfRecHit.phiVtx[hfRecHit.n] = pos.phi();
377	>	hfRecHit.perpVtx[hfRecHit.n] = pos.rho();
378	>
379	>	}
380	>
381		hfRecHit.isjet[hfRecHit.n] = false;
382	<	if(!excludeJets_){
382	>	hfRecHit.depth[hfRecHit.n] = hit.id().depth();
383	>
384	>	if(hit.id().ieta() > 0){
385	>	if(hit.energy() > hfShortThreshold_ && hit.id().depth() != 1) nHFshortPlus++;
386	>	if(hit.energy() > hfLongThreshold_ && hit.id().depth() == 1) nHFlongPlus++;
387	>	}else{
388	>	if(hit.energy() > hfShortThreshold_ && hit.id().depth() != 1) nHFshortMinus++;
389	>	if(hit.energy() > hfLongThreshold_ && hit.id().depth() == 1) nHFlongMinus++;
390	>	}
391	>
392	>	if(useJets_){
393		for(unsigned int j = 0 ; j < jets->size(); ++j){
394		const reco::Jet& jet = (*jets)[j];
395		double dr = reco::deltaR(hfRecHit.eta[hfRecHit.n],hfRecHit.phi[hfRecHit.n],jet.eta(),jet.phi());
396		if(dr < cone){ hfRecHit.isjet[hfRecHit.n] = true; }
397		}
398		}
399	<	hfRecHit.n++;
399	>	if (hfRecHit.et[hfRecHit.n]>=hfPtMin_) hfRecHit.n++;
400		}
401	<
401	>
402	>	if(doHcal_ && !doEbyEonly_){
403		for(unsigned int i = 0; i < hbheHits->size(); ++i){
404		const HBHERecHit & hit= (*hbheHits)[i];
405	+	if (getEt(hit.id(),hit.energy())<hbhePtMin_) continue;
406	+
407		hbheRecHit.e[hbheRecHit.n] = hit.energy();
408	<	hbheRecHit.et[hbheRecHit.n] = getEt(hit.id(),hit.energy());
409	<	hbheRecHit.eta[hbheRecHit.n] = getEta(hit.id());
410	<	hbheRecHit.phi[hbheRecHit.n] = getPhi(hit.id());
411	<	hbheRecHit.isjet[hbheRecHit.n] = false;
412	<	if(!excludeJets_){
408	>	math::XYZPoint pos = getPosition(hit.id(),vtx);
409	>
410	>	if(!saveBothVtx_){
411	>	hbheRecHit.et[hbheRecHit.n] = getEt(pos,hit.energy());
412	>	hbheRecHit.eta[hbheRecHit.n] = pos.eta();
413	>	hbheRecHit.phi[hbheRecHit.n] = pos.phi();
414	>	hbheRecHit.perp[hbheRecHit.n] = pos.rho();
415	>	}else{
416	>	hbheRecHit.et[hbheRecHit.n] = getEt(hit.id(),hit.energy());
417	>	hbheRecHit.eta[hbheRecHit.n] = getEta(hit.id());
418	>	hbheRecHit.phi[hbheRecHit.n] = getPhi(hit.id());
419	>	hbheRecHit.perp[hbheRecHit.n] = getPerp(hit.id());
420	>
421	>	hbheRecHit.etVtx[hbheRecHit.n] = getEt(pos,hit.energy());
422	>	hbheRecHit.etaVtx[hbheRecHit.n] = pos.eta();
423	>	hbheRecHit.phiVtx[hbheRecHit.n] = pos.phi();
424	>	hbheRecHit.perpVtx[hbheRecHit.n] = pos.rho();
425	>
426	>	}
427	>
428	>	hbheRecHit.isjet[hbheRecHit.n] = false;
429	>	hbheRecHit.depth[hbheRecHit.n] = hit.id().depth();
430	>
431	>	if(useJets_){
432		for(unsigned int j = 0 ; j < jets->size(); ++j){
433		const reco::Jet& jet = (*jets)[j];
434		double dr = reco::deltaR(hbheRecHit.eta[hbheRecHit.n],hbheRecHit.phi[hbheRecHit.n],jet.eta(),jet.phi());
#	Line 240 | Line 437 | RecHitTreeProducer::analyze(const edm::E
437		}
438		hbheRecHit.n++;
439		}
440	<
440	>	}
441	>	}
442	>	if(doEcal_ && !doEbyEonly_){
443		for(unsigned int i = 0; i < ebHits->size(); ++i){
444		const EcalRecHit & hit= (*ebHits)[i];
445	+	if (getEt(hit.id(),hit.energy())<ebPtMin_) continue;
446	+
447		ebRecHit.e[ebRecHit.n] = hit.energy();
448	<	ebRecHit.et[ebRecHit.n] = getEt(hit.id(),hit.energy());
449	<	ebRecHit.eta[ebRecHit.n] = getEta(hit.id());
450	<	ebRecHit.phi[ebRecHit.n] = getPhi(hit.id());
448	>	math::XYZPoint pos = getPosition(hit.id(),vtx);
449	>
450	>	if(!saveBothVtx_){
451	>	ebRecHit.et[ebRecHit.n] = getEt(pos,hit.energy());
452	>	ebRecHit.eta[ebRecHit.n] = pos.eta();
453	>	ebRecHit.phi[ebRecHit.n] = pos.phi();
454	>	ebRecHit.perp[ebRecHit.n] = pos.rho();
455	>	}else{
456	>	ebRecHit.et[ebRecHit.n] = getEt(hit.id(),hit.energy());
457	>	ebRecHit.eta[ebRecHit.n] = getEta(hit.id());
458	>	ebRecHit.phi[ebRecHit.n] = getPhi(hit.id());
459	>	ebRecHit.perp[ebRecHit.n] = getPerp(hit.id());
460	>	ebRecHit.etVtx[ebRecHit.n] = getEt(pos,hit.energy());
461	>	ebRecHit.etaVtx[ebRecHit.n] = pos.eta();
462	>	ebRecHit.phiVtx[ebRecHit.n] = pos.phi();
463	>	ebRecHit.perpVtx[ebRecHit.n] = pos.rho();
464	>
465	>	}
466	>
467		ebRecHit.isjet[ebRecHit.n] = false;
468	<	if(!excludeJets_){
468	>	if(useJets_){
469		for(unsigned int j = 0 ; j < jets->size(); ++j){
470		const reco::Jet& jet = (*jets)[j];
471		double dr = reco::deltaR(ebRecHit.eta[ebRecHit.n],ebRecHit.phi[ebRecHit.n],jet.eta(),jet.phi());
#	Line 260 | Line 477 | RecHitTreeProducer::analyze(const edm::E
477
478		for(unsigned int i = 0; i < eeHits->size(); ++i){
479		const EcalRecHit & hit= (*eeHits)[i];
480	+	if (getEt(hit.id(),hit.energy())<eePtMin_) continue;
481	+
482		eeRecHit.e[eeRecHit.n] = hit.energy();
483	<	eeRecHit.et[eeRecHit.n] = getEt(hit.id(),hit.energy());
484	<	eeRecHit.eta[eeRecHit.n] = getEta(hit.id());
485	<	eeRecHit.phi[eeRecHit.n] = getPhi(hit.id());
483	>	math::XYZPoint pos = getPosition(hit.id(),vtx);
484	>
485	>	if(!saveBothVtx_){
486	>	eeRecHit.et[eeRecHit.n] = getEt(pos,hit.energy());
487	>	eeRecHit.eta[eeRecHit.n] = pos.eta();
488	>	eeRecHit.phi[eeRecHit.n] = pos.phi();
489	>	eeRecHit.perp[eeRecHit.n] = pos.rho();
490	>	}else{
491	>	eeRecHit.et[eeRecHit.n] = getEt(hit.id(),hit.energy());
492	>	eeRecHit.eta[eeRecHit.n] = getEta(hit.id());
493	>	eeRecHit.phi[eeRecHit.n] = getPhi(hit.id());
494	>	eeRecHit.perp[eeRecHit.n] = getPerp(hit.id());
495	>	eeRecHit.etVtx[eeRecHit.n] = getEt(pos,hit.energy());
496	>	eeRecHit.etaVtx[eeRecHit.n] = pos.eta();
497	>	eeRecHit.phiVtx[eeRecHit.n] = pos.phi();
498	>	eeRecHit.perpVtx[eeRecHit.n] = pos.rho();
499	>
500	>	}
501	>
502		eeRecHit.isjet[eeRecHit.n] = false;
503	<	if(!excludeJets_){
503	>
504	>	if(useJets_){
505		for(unsigned int j = 0 ; j < jets->size(); ++j){
506		const reco::Jet& jet = (*jets)[j];
507		double dr = reco::deltaR(eeRecHit.eta[eeRecHit.n],eeRecHit.phi[eeRecHit.n],jet.eta(),jet.phi());
#	Line 274 | Line 510 | RecHitTreeProducer::analyze(const edm::E
510		}
511		eeRecHit.n++;
512		}
513	<
514	<	eeTree->Fill();
515	<	ebTree->Fill();
516	<
517	<	hbheTree->Fill();
518	<	hfTree->Fill();
513	>	}
514	>
515	>	if(doTowers_){
516	>
517	>	for(unsigned int i = 0; i < towers->size(); ++i){
518	>	const CaloTower & hit= (*towers)[i];
519	>	if (getEt(hit.id(),hit.energy())<towerPtMin_) continue;
520	>
521	>	myTowers.et[myTowers.n] = hit.p4(vtx).Et();
522	>	myTowers.eta[myTowers.n] = hit.p4(vtx).Eta();
523	>	myTowers.phi[myTowers.n] = hit.p4(vtx).Phi();
524	>	myTowers.emEt[myTowers.n] = hit.emEt(vtx);
525	>	myTowers.hadEt[myTowers.n] = hit.hadEt(vtx);
526	>
527	>	if (saveBothVtx_) {
528	>	myTowers.e[myTowers.n] = hit.energy();
529	>	myTowers.et[myTowers.n] = getEt(hit.id(),hit.energy());
530	>	myTowers.eta[myTowers.n] = getEta(hit.id());
531	>	myTowers.phi[myTowers.n] = getPhi(hit.id());
532	>	myTowers.isjet[myTowers.n] = false;
533	>	myTowers.etVtx[myTowers.n] = hit.p4(vtx).Et();
534	>	myTowers.etaVtx[myTowers.n] = hit.p4(vtx).Eta();
535	>	myTowers.emEtVtx[myTowers.n] = hit.emEt(vtx);
536	>	myTowers.hadEtVtx[myTowers.n] = hit.hadEt(vtx);
537	>	}
538	>
539	>	myTowers.isjet[myTowers.n] = false;
540	>
541	>	if(hit.ieta() > 29 && hit.energy() > hfTowerThreshold_) nHFtowerPlus++;
542	>	if(hit.ieta() < -29 && hit.energy() > hfTowerThreshold_) nHFtowerMinus++;
543	>
544	>	if(useJets_){
545	>	for(unsigned int j = 0 ; j < jets->size(); ++j){
546	>	const reco::Jet& jet = (*jets)[j];
547	>	double dr = reco::deltaR(myTowers.eta[myTowers.n],myTowers.phi[myTowers.n],jet.eta(),jet.phi());
548	>	if(dr < cone){ myTowers.isjet[myTowers.n] = true; }
549	>	}
550	>	}
551	>	myTowers.n++;
552	>	}
553	>
554	>	}
555	>
556	>	if(doBasicClusters_ && !doEbyEonly_){
557	>	for(unsigned int j = 0 ; j < jets->size(); ++j){
558	>	const reco::Jet& jet = (*jets)[j];
559	>	myBC.n = 0;
560	>	myBC.jtpt = jet.pt();
561	>	myBC.jteta = jet.eta();
562	>	myBC.jtphi = jet.phi();
563	>
564	>	for(unsigned int i = 0; i < bClusters->size(); ++i){
565	>	const reco::BasicCluster & bc= (*bClusters)[i];
566	>	double dr = reco::deltaR(bc.eta(),bc.phi(),jet.eta(),jet.phi());
567	>	if(dr < cone){
568	>	myBC.e[myBC.n] = bc.energy();
569	>	myBC.et[myBC.n] = bc.energy()*sin(bc.position().theta());
570	>	myBC.eta[myBC.n] = bc.eta();
571	>	myBC.phi[myBC.n] = bc.phi();
572	>	myBC.n++;
573	>	}
574	>	}
575	>	bcTree->Fill();
576	>	}
577	>	}
578	>
579	>	if(doCastor_){
580	>
581	>	edm::Handle<CastorRecHitCollection> casrechits;
582	>	try{ ev.getByLabel("castorreco",casrechits); }
583	>	catch(...) { edm::LogWarning(" CASTOR ") << " Cannot get Castor RecHits " << std::endl; }
584	>
585	>	int nhits = 0;
586	>	double energyCastor = 0;
587	>
588	>	if(casrechits.failedToGet()!=0 || !casrechits.isValid()) {
589	>	edm::LogWarning(" CASTOR ") << " Cannot read CastorRecHitCollection" << std::endl;
590	>	} else {
591	>	for(size_t i1 = 0; i1 < casrechits->size(); ++i1) {
592	>	const CastorRecHit & rh = (*casrechits)[i1];
593	>	HcalCastorDetId castorid = rh.id();
594	>	energyCastor += rh.energy();
595	>	if (nhits  < 224) {
596	>	castorRecHit.e[nhits] = rh.energy();
597	>	castorRecHit.iphi[nhits] = castorid.sector();
598	>	castorRecHit.depth[nhits] = castorid.module();
599	>	castorRecHit.phi[nhits] = getPhi(castorid);
600	>	}
601	>
602	>	nhits++;
603	>
604	>	} // end loop castor rechits
605	>	}
606	>
607	>	castorRecHit.n = nhits;
608	>	castorTree->Fill();
609	>	}
610	>
611	>
612	>	if(!doEbyEonly_){
613	>	towerTree->Fill();
614	>
615	>	eeTree->Fill();
616	>	ebTree->Fill();
617	>
618	>	hbheTree->Fill();
619	>	hfTree->Fill();
620	>
621	>	if (doFastJet_) {
622	>	bkgTree->Fill();
623	>	}
624	>	}
625	>
626	>	nt->Fill(nHFtowerPlus,nHFtowerMinus,nHFlongPlus,nHFlongMinus,nHFshortPlus,nHFshortMinus);
627
628		}
629
#	Line 289 | Line 633 | void
633		RecHitTreeProducer::beginJob()
634		{
635
636	<	hbheTree = fs->make<TTree>("hbhe","");
636	>	hbheTree = fs->make<TTree>("hbhe",versionTag);
637		hbheTree->Branch("n",&hbheRecHit.n,"n/I");
638		hbheTree->Branch("e",hbheRecHit.e,"e[n]/F");
639		hbheTree->Branch("et",hbheRecHit.et,"et[n]/F");
640		hbheTree->Branch("eta",hbheRecHit.eta,"eta[n]/F");
641		hbheTree->Branch("phi",hbheRecHit.phi,"phi[n]/F");
642	<	hbheTree->Branch("isjet",hbheRecHit.isjet,"isjet[n]/I");
643	<
644	<	hfTree = fs->make<TTree>("hf","");
642	>	hbheTree->Branch("perp",hbheRecHit.perp,"perp[n]/F");
643	>
644	>	hbheTree->Branch("isjet",hbheRecHit.isjet,"isjet[n]/O");
645	>	hbheTree->Branch("depth",hfRecHit.depth,"depth[n]/I");
646	>
647	>	hfTree = fs->make<TTree>("hf",versionTag);
648		hfTree->Branch("n",&hfRecHit.n,"n/I");
649		hfTree->Branch("e",hfRecHit.e,"e[n]/F");
650		hfTree->Branch("et",hfRecHit.et,"et[n]/F");
651		hfTree->Branch("eta",hfRecHit.eta,"eta[n]/F");
652		hfTree->Branch("phi",hfRecHit.phi,"phi[n]/F");
653	<	hfTree->Branch("isjet",hfRecHit.isjet,"isjet[n]/I");
653	>	hfTree->Branch("perp",hfRecHit.perp,"perp[n]/F");
654	>	hfTree->Branch("depth",hfRecHit.depth,"depth[n]/I");
655	>	hfTree->Branch("isjet",hfRecHit.isjet,"isjet[n]/O");
656
657	<	eeTree = fs->make<TTree>("ee","");
657	>	eeTree = fs->make<TTree>("ee",versionTag);
658		eeTree->Branch("n",&eeRecHit.n,"n/I");
659		eeTree->Branch("e",eeRecHit.e,"e[n]/F");
660		eeTree->Branch("et",eeRecHit.et,"et[n]/F");
661		eeTree->Branch("eta",eeRecHit.eta,"eta[n]/F");
662		eeTree->Branch("phi",eeRecHit.phi,"phi[n]/F");
663	<	eeTree->Branch("isjet",eeRecHit.isjet,"isjet[n]/I");
663	>	eeTree->Branch("perp",eeRecHit.perp,"perp[n]/F");
664	>
665	>	eeTree->Branch("isjet",eeRecHit.isjet,"isjet[n]/O");
666
667	<	ebTree = fs->make<TTree>("eb","");
667	>	ebTree = fs->make<TTree>("eb",versionTag);
668		ebTree->Branch("n",&ebRecHit.n,"n/I");
669		ebTree->Branch("e",ebRecHit.e,"e[n]/F");
670		ebTree->Branch("et",ebRecHit.et,"et[n]/F");
671		ebTree->Branch("eta",ebRecHit.eta,"eta[n]/F");
672		ebTree->Branch("phi",ebRecHit.phi,"phi[n]/F");
673	<	ebTree->Branch("isjet",ebRecHit.isjet,"isjet[n]/I");
673	>	ebTree->Branch("perp",ebRecHit.perp,"perp[n]/F");
674	>
675	>	ebTree->Branch("isjet",ebRecHit.isjet,"isjet[n]/O");
676	>
677	>	towerTree = fs->make<TTree>("tower",versionTag);
678	>	towerTree->Branch("n",&myTowers.n,"n/I");
679	>	towerTree->Branch("e",myTowers.e,"e[n]/F");
680	>	towerTree->Branch("et",myTowers.et,"et[n]/F");
681	>	towerTree->Branch("eta",myTowers.eta,"eta[n]/F");
682	>	towerTree->Branch("phi",myTowers.phi,"phi[n]/F");
683	>	towerTree->Branch("isjet",myTowers.isjet,"isjet[n]/O");
684	>	towerTree->Branch("emEt",myTowers.emEt,"emEt[n]/F");
685	>	towerTree->Branch("hadEt",myTowers.hadEt,"hadEt[n]/F");
686	>
687	>
688	>	if(doCastor_){
689	>	castorTree = fs->make<TTree>("castor",versionTag);
690	>	castorTree->Branch("n",&castorRecHit.n,"n/I");
691	>	castorTree->Branch("e",castorRecHit.e,"e[n]/F");
692	>	castorTree->Branch("iphi",castorRecHit.iphi,"iphi[n]/I");
693	>	castorTree->Branch("phi",castorRecHit.phi,"phi[n]/F");
694	>	castorTree->Branch("depth",castorRecHit.depth,"depth[n]/I");
695	>	}
696	>
697	>	if (saveBothVtx_) {
698	>	towerTree->Branch("etVtx",myTowers.etVtx,"etvtx[n]/F");
699	>	towerTree->Branch("etaVtx",myTowers.etaVtx,"etavtx[n]/F");
700	>	towerTree->Branch("emEtVtx",myTowers.emEtVtx,"emEtVtx[n]/F");
701	>	towerTree->Branch("hadEtVtx",myTowers.hadEtVtx,"hadEtVtx[n]/F");
702	>	}
703	>
704	>
705	>	if(doBasicClusters_){
706	>	bcTree = fs->make<TTree>("bc",versionTag);
707	>	bcTree->Branch("n",&myBC.n,"n/I");
708	>	bcTree->Branch("e",myBC.e,"e[n]/F");
709	>	bcTree->Branch("et",myBC.et,"et[n]/F");
710	>	bcTree->Branch("eta",myBC.eta,"eta[n]/F");
711	>	bcTree->Branch("phi",myBC.phi,"phi[n]/F");
712	>	bcTree->Branch("jtpt",&myBC.jtpt,"jtpt/F");
713	>	bcTree->Branch("jteta",&myBC.jteta,"jteta/F");
714	>	bcTree->Branch("jtphi",&myBC.jtphi,"jtphi/F");
715	>	//     bcTree->Branch("isjet",bcRecHit.isjet,"isjet[n]/O");
716	>	}
717	>
718	>	if(doFastJet_){
719	>	bkgTree = fs->make<TTree>("bkg",versionTag);
720	>	bkgTree->Branch("n",&bkg.n,"n/I");
721	>	bkgTree->Branch("rho",bkg.rho,"rho[n]/F");
722	>	bkgTree->Branch("sigma",bkg.sigma,"sigma[n]/F");
723	>	}
724	>
725	>	nt = fs->make<TNtuple>("ntEvent","","nHFplus:nHFminus:nHFlongPlus:nHFlongMinus:nHFshortPlus:nHFshortMinus");
726
727		}
728
#	Line 328 | Line 731 | void
731		RecHitTreeProducer::endJob() {
732		}
733
734	<	double RecHitTreeProducer::getEt(const DetId &id, double energy){
734	>	math::XYZPoint RecHitTreeProducer::getPosition(const DetId &id, reco::Vertex::Point vtx){
735	>	const GlobalPoint& pos=geo->getPosition(id);
736	>	math::XYZPoint posV(pos.x() - vtx.x(),pos.y() - vtx.y(),pos.z() - vtx.z());
737	>	return posV;
738	>	}
739	>
740	>	double RecHitTreeProducer::getEt(math::XYZPoint pos, double energy){
741	>	double et = energy*sin(pos.theta());
742	>	return et;
743	>	}
744	>
745	>	double RecHitTreeProducer::getEt(const DetId &id, double energy, reco::Vertex::Point vtx){
746		const GlobalPoint& pos=geo->getPosition(id);
747		double et = energy*sin(pos.theta());
748		return et;
749		}
750
751	<	double RecHitTreeProducer::getEta(const DetId &id){
751	>	double RecHitTreeProducer::getEta(const DetId &id, reco::Vertex::Point vtx){
752		const GlobalPoint& pos=geo->getPosition(id);
753		double et = pos.eta();
754		return et;
755		}
756
757	<	double RecHitTreeProducer::getPhi(const DetId &id){
757	>	double RecHitTreeProducer::getPhi(const DetId &id, reco::Vertex::Point vtx){
758		const GlobalPoint& pos=geo->getPosition(id);
759		double et = pos.phi();
760		return et;
761		}
762
763	+	double RecHitTreeProducer::getPerp(const DetId &id, reco::Vertex::Point vtx){
764	+	const GlobalPoint& pos=geo->getPosition(id);
765	+	double et = pos.perp();
766	+	return et;
767	+	}
768
769	+	reco::Vertex::Point RecHitTreeProducer::getVtx(const edm::Event& ev)
770	+	{
771	+	ev.getByLabel(VtxSrc_,vtxs);
772	+	int greatestvtx = 0;
773	+	int nVertex = vtxs->size();
774	+
775	+	for (unsigned int i = 0 ; i< vtxs->size(); ++i){
776	+	unsigned int daughter = (*vtxs)[i].tracksSize();
777	+	if( daughter > (*vtxs)[greatestvtx].tracksSize()) greatestvtx = i;
778	+	//cout <<"Vertex: "<< (*vtxs)[i].position().z()<<" "<<daughter<<endl;
779	+	}
780	+
781	+	if(nVertex<=0){
782	+	return reco::Vertex::Point(0,0,0);
783	+	}
784	+	return (*vtxs)[greatestvtx].position();
785	+	}
786
787		//define this as a plug-in
788		DEFINE_FWK_MODULE(RecHitTreeProducer);

Diff Legend

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