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.4 by yilmaz, Fri Oct 22 14:06:15 2010 UTC vs.
Revision 1.23 by yilmaz, Sat Jan 19 16:25:13 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 55 | Line 56
56
57		#include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
58		#include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
59	+	#include "DataFormats/HcalDigi/interface/HcalDigiCollections.h"
60		#include "DataFormats/HcalDetId/interface/HcalDetId.h"
61		#include "DataFormats/EcalDetId/interface/EcalDetIdCollections.h"
62	+	#include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
63		#include "DataFormats/DetId/interface/DetId.h"
64
65	+	#include "DataFormats/VertexReco/interface/Vertex.h"
66	+	#include "DataFormats/VertexReco/interface/VertexFwd.h"
67	+
68
69		#include "TNtuple.h"
70
71		using namespace std;
72
73	<	#define MAXHITS 1000000
73	>	#define MAXHITS 100000
74
69	–	struct MyRecHit{
75
76	+	struct MyRecHit{
77	+	int depth[MAXHITS];
78		int n;
79
80	+	int ieta[MAXHITS];
81	+	int iphi[MAXHITS];
82	+
83		float e[MAXHITS];
84		float et[MAXHITS];
85		float eta[MAXHITS];
86		float phi[MAXHITS];
87	+	float perp[MAXHITS];
88	+	float emEt[MAXHITS];
89	+	float hadEt[MAXHITS];
90	+
91		bool isjet[MAXHITS];
92	+	float etVtx[MAXHITS];
93	+	float etaVtx[MAXHITS];
94	+	float phiVtx[MAXHITS];
95	+	float perpVtx[MAXHITS];
96	+	float emEtVtx[MAXHITS];
97	+	float hadEtVtx[MAXHITS];
98	+
99	+	float jtpt;
100	+	float jteta;
101	+	float jtphi;
102
103		};
104
105	+	struct MyZDCRecHit{
106	+	int n;
107	+	float  e[18];
108	+	int    zside[18];
109	+	int    section [18];
110	+	int    channel[18];
111	+	int    saturation[18];
112	+	};
113	+
114	+	struct MyZDCDigi{
115	+	int    n;
116	+	int    nts;
117	+	float  chargefC[18][10];
118	+	int    adc[18][10];
119	+	int    zside[18][10];
120	+	int    section [18][10];
121	+	int    channel[18][10];
122	+	int    ts[18][10];
123	+	};
124	+
125	+	struct MyBkg{
126	+	int n;
127	+	float rho[50];
128	+	float sigma[50];
129	+	};
130
131
132		//
#	Line 88 | Line 137 | class RecHitTreeProducer : public edm::E
137		public:
138		explicit RecHitTreeProducer(const edm::ParameterSet&);
139		~RecHitTreeProducer();
140	<	double       getEt(const DetId &id, double energy);
141	<	double       getEta(const DetId &id);
142	<	double       getPhi(const DetId &id);
140	>	double       getEt(const DetId &id, double energy, reco::Vertex::Point vtx = reco::Vertex::Point(0,0,0));
141	>	double       getEta(const DetId &id, reco::Vertex::Point vtx = reco::Vertex::Point(0,0,0));
142	>	double       getPhi(const DetId &id, reco::Vertex::Point vtx = reco::Vertex::Point(0,0,0));
143	>	double       getPerp(const DetId &id, reco::Vertex::Point vtx = reco::Vertex::Point(0,0,0));
144	>
145	>	math::XYZPoint getPosition(const DetId &id, reco::Vertex::Point vtx = reco::Vertex::Point(0,0,0));
146	>	double getEt(math::XYZPoint pos, double energy);
147	>
148	>	reco::Vertex::Point getVtx(const edm::Event& ev);
149	>
150
151
152		private:
#	Line 110 | Line 166 | class RecHitTreeProducer : public edm::E
166		edm::Handle<HFRecHitCollection> hfHits;
167		edm::Handle<HBHERecHitCollection> hbheHits;
168
169	+	edm::Handle<reco::BasicClusterCollection> bClusters;
170	+	edm::Handle<CaloTowerCollection> towers;
171	+	edm::Handle<reco::VertexCollection> vtxs;
172	+
173		typedef vector<EcalRecHit>::const_iterator EcalIterator;
174
175		edm::Handle<reco::CaloJetCollection> jets;
176	+
177	+	edm::Handle<std::vector<double> > rhos;
178	+	edm::Handle<std::vector<double> > sigmas;
179
180		MyRecHit hbheRecHit;
181		MyRecHit hfRecHit;
182		MyRecHit ebRecHit;
183		MyRecHit eeRecHit;
184	+	MyRecHit myBC;
185	+	MyRecHit myTowers;
186	+	MyRecHit castorRecHit;
187	+
188	+	MyZDCRecHit zdcRecHit;
189	+	MyZDCDigi zdcDigi;
190	+
191	+	MyBkg bkg;
192
193	+	TNtuple* nt;
194		TTree* hbheTree;
195		TTree* hfTree;
196		TTree* ebTree;
197		TTree* eeTree;
198	<	double cone;
198	>	TTree* bcTree;
199	>	TTree* towerTree;
200	>	TTree* bkgTree;
201	>	TTree* castorTree;
202	>	TTree* zdcRecHitTree;
203	>	TTree* zdcDigiTree;
204
205	+	double cone;
206	+	double hfTowerThreshold_;
207	+	double hfLongThreshold_;
208	+	double hfShortThreshold_;
209	+	double hbheThreshold_;
210	+	double ebThreshold_;
211	+	double eeThreshold_;
212	+
213	+	double hbhePtMin_;
214	+	double hfPtMin_;
215	+	double ebPtMin_;
216	+	double eePtMin_;
217	+	double towerPtMin_;
218	+
219		edm::Service<TFileService> fs;
220		const CentralityBins * cbins_;
221		const CaloGeometry *geo;
#	Line 133 | Line 224 | class RecHitTreeProducer : public edm::E
224		edm::InputTag HcalRecHitHBHESrc_;
225		edm::InputTag EBSrc_;
226		edm::InputTag EESrc_;
227	+	edm::InputTag BCSrc_;
228	+	edm::InputTag TowerSrc_;
229	+	edm::InputTag VtxSrc_;
230	+
231		edm::InputTag JetSrc_;
232	<	bool excludeJets_;
232	>
233	>	edm::InputTag FastJetTag_;
234	>
235	>	bool useJets_;
236	>	bool doBasicClusters_;
237	>	bool doTowers_;
238	>	bool doEcal_;
239	>	bool doHcal_;
240	>	bool doHF_;
241	>	bool doCastor_;
242	>	bool doZDCRecHit_;
243	>	bool doZDCDigi_;
244	>
245	>	bool hasVtx_;
246	>	bool saveBothVtx_;
247	>
248	>	bool doFastJet_;
249	>
250	>	bool doEbyEonly_;
251
252		};
253
#	Line 155 | Line 268 | RecHitTreeProducer::RecHitTreeProducer(c
268		cone(0.5)
269		{
270		//now do what ever initialization is needed
158	–
271		EBSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("EBRecHitSrc",edm::InputTag("ecalRecHit","EcalRecHitsEB"));
272		EESrc_ = iConfig.getUntrackedParameter<edm::InputTag>("EERecHitSrc",edm::InputTag("ecalRecHit","EcalRecHitsEE"));
273		HcalRecHitHFSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHFRecHitSrc",edm::InputTag("hfreco"));
274		HcalRecHitHBHESrc_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHBHERecHitSrc",edm::InputTag("hbhereco"));
275	<	JetSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("JetSrc",edm::InputTag("iterativeCone5CaloJets"));
276	<	excludeJets_ = iConfig.getUntrackedParameter<bool>("excludeJets",false);
275	>	BCSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("BasicClusterSrc1",edm::InputTag("ecalRecHit","EcalRecHitsEB","RECO"));
276	>	TowerSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("towersSrc",edm::InputTag("towerMaker"));
277	>	VtxSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("vtxSrc",edm::InputTag("hiSelectedVertex"));
278	>	JetSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("JetSrc",edm::InputTag("iterativeConePu5CaloJets"));
279	>	useJets_ = iConfig.getUntrackedParameter<bool>("useJets",true);
280	>	doBasicClusters_ = iConfig.getUntrackedParameter<bool>("doBasicClusters",false);
281	>	doTowers_ = iConfig.getUntrackedParameter<bool>("doTowers",true);
282	>	doEcal_ = iConfig.getUntrackedParameter<bool>("doEcal",true);
283	>	doHcal_ = iConfig.getUntrackedParameter<bool>("doHcal",true);
284	>	doHF_ = iConfig.getUntrackedParameter<bool>("doHF",true);
285	>	doCastor_ = iConfig.getUntrackedParameter<bool>("doCASTOR",true);
286	>	doZDCRecHit_ = iConfig.getUntrackedParameter<bool>("doZDCRecHit",true);
287	>	doZDCDigi_ = iConfig.getUntrackedParameter<bool>("doZDCDigi",true);
288	>
289	>	hasVtx_ = iConfig.getUntrackedParameter<bool>("hasVtx",true);
290	>	saveBothVtx_ = iConfig.getUntrackedParameter<bool>("saveBothVtx",false);
291	>
292	>	doFastJet_ = iConfig.getUntrackedParameter<bool>("doFastJet",true);
293	>	FastJetTag_ = iConfig.getUntrackedParameter<edm::InputTag>("FastJetTag",edm::InputTag("kt4CaloJets"));
294	>	doEbyEonly_ = iConfig.getUntrackedParameter<bool>("doEbyEonly",false);
295	>	hfTowerThreshold_ = iConfig.getUntrackedParameter<double>("HFtowerMin",3.);
296	>	hfLongThreshold_ = iConfig.getUntrackedParameter<double>("HFlongMin",0.5);
297	>	hfShortThreshold_ = iConfig.getUntrackedParameter<double>("HFshortMin",0.85);
298	>	hbhePtMin_ = iConfig.getUntrackedParameter<double>("HBHETreePtMin",0);
299	>	hfPtMin_ = iConfig.getUntrackedParameter<double>("HFTreePtMin",0);
300	>	ebPtMin_ = iConfig.getUntrackedParameter<double>("EBTreePtMin",0);
301	>	eePtMin_ = iConfig.getUntrackedParameter<double>("EETreePtMin",0.);
302	>	towerPtMin_ = iConfig.getUntrackedParameter<double>("TowerTreePtMin",0.);
303	>
304	>
305		}
306
307
#	Line 187 | Line 327 | RecHitTreeProducer::analyze(const edm::E
327		hbheRecHit.n = 0;
328		ebRecHit.n = 0;
329		eeRecHit.n = 0;
330	+	myBC.n = 0;
331	+	myTowers.n = 0;
332	+	bkg.n = 0;
333	+
334	+	// get vertex
335	+	reco::Vertex::Point vtx(0,0,0);
336	+	if (hasVtx_) vtx = getVtx(ev);
337	+
338	+	if(doEcal_){
339	+	ev.getByLabel(EBSrc_,ebHits);
340	+	ev.getByLabel(EESrc_,eeHits);
341	+	}
342
343	<	ev.getByLabel(EBSrc_,ebHits);
344	<	ev.getByLabel(EESrc_,eeHits);
343	>	if(doHcal_){
344	>	ev.getByLabel(HcalRecHitHBHESrc_,hbheHits);
345	>	}
346	>	if(doHF_){
347	>	ev.getByLabel(HcalRecHitHFSrc_,hfHits);
348	>	}
349	>
350	>	if(useJets_) {
351	>	ev.getByLabel(JetSrc_,jets);
352	>	}
353
354	<	ev.getByLabel(HcalRecHitHFSrc_,hfHits);
355	<	ev.getByLabel(HcalRecHitHBHESrc_,hbheHits);
354	>	if(doBasicClusters_){
355	>	ev.getByLabel(BCSrc_,bClusters);
356	>	}
357
358	+	if(doTowers_){
359	+	ev.getByLabel(TowerSrc_,towers);
360	+	}
361
362	<	if(!excludeJets_) {
363	<	ev.getByLabel(JetSrc_,jets);
362	>	if(doFastJet_){
363	>	ev.getByLabel(edm::InputTag(FastJetTag_.label(),"rhos",FastJetTag_.process()),rhos);
364	>	ev.getByLabel(edm::InputTag(FastJetTag_.label(),"sigmas",FastJetTag_.process()),sigmas);
365	>	bkg.n = rhos->size();
366	>	for(unsigned int i = 0; i < rhos->size(); ++i){
367	>	bkg.rho[i] = (*rhos)[i];
368	>	bkg.sigma[i] = (*sigmas)[i];
369	>	}
370		}
371
372	<	if(0 && !cbins_) cbins_ = getCentralityBinsFromDB(iSetup);
372	>	if(0 && !cbins_) cbins_ = getCentralityBinsFromDB(iSetup);
373
374		if(!geo){
375		edm::ESHandle<CaloGeometry> pGeo;
376		iSetup.get<CaloGeometryRecord>().get(pGeo);
377		geo = pGeo.product();
378		}
379	<
379	>
380	>	int nHFlongPlus = 0;
381	>	int nHFshortPlus = 0;
382	>	int nHFtowerPlus = 0;
383	>	int nHFlongMinus = 0;
384	>	int nHFshortMinus = 0;
385	>	int nHFtowerMinus = 0;
386	>
387	>	if(doHF_){
388		for(unsigned int i = 0; i < hfHits->size(); ++i){
389		const HFRecHit & hit= (*hfHits)[i];
390		hfRecHit.e[hfRecHit.n] = hit.energy();
391	<	hfRecHit.et[hfRecHit.n] = getEt(hit.id(),hit.energy());
392	<	hfRecHit.eta[hfRecHit.n] = getEta(hit.id());
393	<	hfRecHit.phi[hfRecHit.n] = getPhi(hit.id());
391	>	math::XYZPoint pos = getPosition(hit.id(),vtx);
392	>
393	>	if(!saveBothVtx_){
394	>	hfRecHit.et[hfRecHit.n] = getEt(pos,hit.energy());
395	>	hfRecHit.eta[hfRecHit.n] = pos.eta();
396	>	hfRecHit.phi[hfRecHit.n] = pos.phi();
397	>	hfRecHit.perp[hfRecHit.n] = pos.rho();
398	>	}else{
399	>	hfRecHit.et[hfRecHit.n] = getEt(hit.id(),hit.energy());
400	>	hfRecHit.eta[hfRecHit.n] = getEta(hit.id());
401	>	hfRecHit.phi[hfRecHit.n] = getPhi(hit.id());
402	>	hfRecHit.perp[hfRecHit.n] = getPerp(hit.id());
403	>
404	>	hfRecHit.etVtx[hfRecHit.n] = getEt(pos,hit.energy());
405	>	hfRecHit.etaVtx[hfRecHit.n] = pos.eta();
406	>	hfRecHit.phiVtx[hfRecHit.n] = pos.phi();
407	>	hfRecHit.perpVtx[hfRecHit.n] = pos.rho();
408	>
409	>	}
410	>
411		hfRecHit.isjet[hfRecHit.n] = false;
412	<	if(!excludeJets_){
412	>	hfRecHit.depth[hfRecHit.n] = hit.id().depth();
413	>
414	>	if(hit.id().ieta() > 0){
415	>	if(hit.energy() > hfShortThreshold_ && hit.id().depth() != 1) nHFshortPlus++;
416	>	if(hit.energy() > hfLongThreshold_ && hit.id().depth() == 1) nHFlongPlus++;
417	>	}else{
418	>	if(hit.energy() > hfShortThreshold_ && hit.id().depth() != 1) nHFshortMinus++;
419	>	if(hit.energy() > hfLongThreshold_ && hit.id().depth() == 1) nHFlongMinus++;
420	>	}
421	>
422	>	if(useJets_){
423		for(unsigned int j = 0 ; j < jets->size(); ++j){
424		const reco::Jet& jet = (*jets)[j];
425		double dr = reco::deltaR(hfRecHit.eta[hfRecHit.n],hfRecHit.phi[hfRecHit.n],jet.eta(),jet.phi());
426		if(dr < cone){ hfRecHit.isjet[hfRecHit.n] = true; }
427		}
428		}
429	<	hfRecHit.n++;
429	>	if (hfRecHit.et[hfRecHit.n]>=hfPtMin_) hfRecHit.n++;
430		}
431	<
431	>
432	>	if(doHcal_ && !doEbyEonly_){
433		for(unsigned int i = 0; i < hbheHits->size(); ++i){
434		const HBHERecHit & hit= (*hbheHits)[i];
435	+	if (getEt(hit.id(),hit.energy())<hbhePtMin_) continue;
436	+
437		hbheRecHit.e[hbheRecHit.n] = hit.energy();
438	<	hbheRecHit.et[hbheRecHit.n] = getEt(hit.id(),hit.energy());
439	<	hbheRecHit.eta[hbheRecHit.n] = getEta(hit.id());
440	<	hbheRecHit.phi[hbheRecHit.n] = getPhi(hit.id());
441	<	hbheRecHit.isjet[hbheRecHit.n] = false;
442	<	if(!excludeJets_){
438	>	math::XYZPoint pos = getPosition(hit.id(),vtx);
439	>
440	>	if(!saveBothVtx_){
441	>	hbheRecHit.et[hbheRecHit.n] = getEt(pos,hit.energy());
442	>	hbheRecHit.eta[hbheRecHit.n] = pos.eta();
443	>	hbheRecHit.phi[hbheRecHit.n] = pos.phi();
444	>	hbheRecHit.perp[hbheRecHit.n] = pos.rho();
445	>	}else{
446	>	hbheRecHit.et[hbheRecHit.n] = getEt(hit.id(),hit.energy());
447	>	hbheRecHit.eta[hbheRecHit.n] = getEta(hit.id());
448	>	hbheRecHit.phi[hbheRecHit.n] = getPhi(hit.id());
449	>	hbheRecHit.perp[hbheRecHit.n] = getPerp(hit.id());
450	>
451	>	hbheRecHit.etVtx[hbheRecHit.n] = getEt(pos,hit.energy());
452	>	hbheRecHit.etaVtx[hbheRecHit.n] = pos.eta();
453	>	hbheRecHit.phiVtx[hbheRecHit.n] = pos.phi();
454	>	hbheRecHit.perpVtx[hbheRecHit.n] = pos.rho();
455	>
456	>	}
457	>
458	>	hbheRecHit.isjet[hbheRecHit.n] = false;
459	>	hbheRecHit.depth[hbheRecHit.n] = hit.id().depth();
460	>
461	>	if(useJets_){
462		for(unsigned int j = 0 ; j < jets->size(); ++j){
463		const reco::Jet& jet = (*jets)[j];
464		double dr = reco::deltaR(hbheRecHit.eta[hbheRecHit.n],hbheRecHit.phi[hbheRecHit.n],jet.eta(),jet.phi());
#	Line 240 | Line 467 | RecHitTreeProducer::analyze(const edm::E
467		}
468		hbheRecHit.n++;
469		}
470	<
470	>	}
471	>	}
472	>	if(doEcal_ && !doEbyEonly_){
473		for(unsigned int i = 0; i < ebHits->size(); ++i){
474		const EcalRecHit & hit= (*ebHits)[i];
475	+	if (getEt(hit.id(),hit.energy())<ebPtMin_) continue;
476	+
477		ebRecHit.e[ebRecHit.n] = hit.energy();
478	<	ebRecHit.et[ebRecHit.n] = getEt(hit.id(),hit.energy());
479	<	ebRecHit.eta[ebRecHit.n] = getEta(hit.id());
480	<	ebRecHit.phi[ebRecHit.n] = getPhi(hit.id());
478	>	math::XYZPoint pos = getPosition(hit.id(),vtx);
479	>
480	>	if(!saveBothVtx_){
481	>	ebRecHit.et[ebRecHit.n] = getEt(pos,hit.energy());
482	>	ebRecHit.eta[ebRecHit.n] = pos.eta();
483	>	ebRecHit.phi[ebRecHit.n] = pos.phi();
484	>	ebRecHit.perp[ebRecHit.n] = pos.rho();
485	>	}else{
486	>	ebRecHit.et[ebRecHit.n] = getEt(hit.id(),hit.energy());
487	>	ebRecHit.eta[ebRecHit.n] = getEta(hit.id());
488	>	ebRecHit.phi[ebRecHit.n] = getPhi(hit.id());
489	>	ebRecHit.perp[ebRecHit.n] = getPerp(hit.id());
490	>	ebRecHit.etVtx[ebRecHit.n] = getEt(pos,hit.energy());
491	>	ebRecHit.etaVtx[ebRecHit.n] = pos.eta();
492	>	ebRecHit.phiVtx[ebRecHit.n] = pos.phi();
493	>	ebRecHit.perpVtx[ebRecHit.n] = pos.rho();
494	>
495	>	}
496	>
497		ebRecHit.isjet[ebRecHit.n] = false;
498	<	if(!excludeJets_){
498	>	if(useJets_){
499		for(unsigned int j = 0 ; j < jets->size(); ++j){
500		const reco::Jet& jet = (*jets)[j];
501		double dr = reco::deltaR(ebRecHit.eta[ebRecHit.n],ebRecHit.phi[ebRecHit.n],jet.eta(),jet.phi());
#	Line 260 | Line 507 | RecHitTreeProducer::analyze(const edm::E
507
508		for(unsigned int i = 0; i < eeHits->size(); ++i){
509		const EcalRecHit & hit= (*eeHits)[i];
510	+	if (getEt(hit.id(),hit.energy())<eePtMin_) continue;
511	+
512		eeRecHit.e[eeRecHit.n] = hit.energy();
513	<	eeRecHit.et[eeRecHit.n] = getEt(hit.id(),hit.energy());
514	<	eeRecHit.eta[eeRecHit.n] = getEta(hit.id());
515	<	eeRecHit.phi[eeRecHit.n] = getPhi(hit.id());
513	>	math::XYZPoint pos = getPosition(hit.id(),vtx);
514	>
515	>	if(!saveBothVtx_){
516	>	eeRecHit.et[eeRecHit.n] = getEt(pos,hit.energy());
517	>	eeRecHit.eta[eeRecHit.n] = pos.eta();
518	>	eeRecHit.phi[eeRecHit.n] = pos.phi();
519	>	eeRecHit.perp[eeRecHit.n] = pos.rho();
520	>	}else{
521	>	eeRecHit.et[eeRecHit.n] = getEt(hit.id(),hit.energy());
522	>	eeRecHit.eta[eeRecHit.n] = getEta(hit.id());
523	>	eeRecHit.phi[eeRecHit.n] = getPhi(hit.id());
524	>	eeRecHit.perp[eeRecHit.n] = getPerp(hit.id());
525	>	eeRecHit.etVtx[eeRecHit.n] = getEt(pos,hit.energy());
526	>	eeRecHit.etaVtx[eeRecHit.n] = pos.eta();
527	>	eeRecHit.phiVtx[eeRecHit.n] = pos.phi();
528	>	eeRecHit.perpVtx[eeRecHit.n] = pos.rho();
529	>
530	>	}
531	>
532		eeRecHit.isjet[eeRecHit.n] = false;
533	<	if(!excludeJets_){
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(eeRecHit.eta[eeRecHit.n],eeRecHit.phi[eeRecHit.n],jet.eta(),jet.phi());
#	Line 274 | Line 540 | RecHitTreeProducer::analyze(const edm::E
540		}
541		eeRecHit.n++;
542		}
543	<
544	<	eeTree->Fill();
545	<	ebTree->Fill();
546	<
547	<	hbheTree->Fill();
548	<	hfTree->Fill();
543	>	}
544	>
545	>	if(doTowers_){
546	>
547	>	for(unsigned int i = 0; i < towers->size(); ++i){
548	>	const CaloTower & hit= (*towers)[i];
549	>	if (getEt(hit.id(),hit.energy())<towerPtMin_) continue;
550	>
551	>	myTowers.et[myTowers.n] = hit.p4(vtx).Et();
552	>	myTowers.eta[myTowers.n] = hit.p4(vtx).Eta();
553	>	myTowers.phi[myTowers.n] = hit.p4(vtx).Phi();
554	>	myTowers.emEt[myTowers.n] = hit.emEt(vtx);
555	>	myTowers.hadEt[myTowers.n] = hit.hadEt(vtx);
556	>
557	>	if (saveBothVtx_) {
558	>	myTowers.e[myTowers.n] = hit.energy();
559	>	myTowers.et[myTowers.n] = getEt(hit.id(),hit.energy());
560	>	myTowers.eta[myTowers.n] = getEta(hit.id());
561	>	myTowers.phi[myTowers.n] = getPhi(hit.id());
562	>	myTowers.isjet[myTowers.n] = false;
563	>	myTowers.etVtx[myTowers.n] = hit.p4(vtx).Et();
564	>	myTowers.etaVtx[myTowers.n] = hit.p4(vtx).Eta();
565	>	myTowers.emEtVtx[myTowers.n] = hit.emEt(vtx);
566	>	myTowers.hadEtVtx[myTowers.n] = hit.hadEt(vtx);
567	>	}
568	>
569	>	myTowers.isjet[myTowers.n] = false;
570	>
571	>	if(hit.ieta() > 29 && hit.energy() > hfTowerThreshold_) nHFtowerPlus++;
572	>	if(hit.ieta() < -29 && hit.energy() > hfTowerThreshold_) nHFtowerMinus++;
573	>
574	>	if(useJets_){
575	>	for(unsigned int j = 0 ; j < jets->size(); ++j){
576	>	const reco::Jet& jet = (*jets)[j];
577	>	double dr = reco::deltaR(myTowers.eta[myTowers.n],myTowers.phi[myTowers.n],jet.eta(),jet.phi());
578	>	if(dr < cone){ myTowers.isjet[myTowers.n] = true; }
579	>	}
580	>	}
581	>	myTowers.n++;
582	>	}
583	>
584	>	}
585	>
586	>	if(doBasicClusters_ && !doEbyEonly_){
587	>	for(unsigned int j = 0 ; j < jets->size(); ++j){
588	>	const reco::Jet& jet = (*jets)[j];
589	>	myBC.n = 0;
590	>	myBC.jtpt = jet.pt();
591	>	myBC.jteta = jet.eta();
592	>	myBC.jtphi = jet.phi();
593	>
594	>	for(unsigned int i = 0; i < bClusters->size(); ++i){
595	>	const reco::BasicCluster & bc= (*bClusters)[i];
596	>	double dr = reco::deltaR(bc.eta(),bc.phi(),jet.eta(),jet.phi());
597	>	if(dr < cone){
598	>	myBC.e[myBC.n] = bc.energy();
599	>	myBC.et[myBC.n] = bc.energy()*sin(bc.position().theta());
600	>	myBC.eta[myBC.n] = bc.eta();
601	>	myBC.phi[myBC.n] = bc.phi();
602	>	myBC.n++;
603	>	}
604	>	}
605	>	bcTree->Fill();
606	>	}
607	>	}
608	>
609	>	if(doCastor_){
610	>
611	>	edm::Handle<CastorRecHitCollection> casrechits;
612	>	try{ ev.getByLabel("castorreco",casrechits); }
613	>	catch(...) { edm::LogWarning(" CASTOR ") << " Cannot get Castor RecHits " << std::endl; }
614	>
615	>	int nhits = 0;
616	>	double energyCastor = 0;
617	>
618	>	if(casrechits.failedToGet()!=0 || !casrechits.isValid()) {
619	>	edm::LogWarning(" CASTOR ") << " Cannot read CastorRecHitCollection" << std::endl;
620	>	} else {
621	>	for(size_t i1 = 0; i1 < casrechits->size(); ++i1) {
622	>	const CastorRecHit & rh = (*casrechits)[i1];
623	>	HcalCastorDetId castorid = rh.id();
624	>	energyCastor += rh.energy();
625	>	if (nhits  < 224) {
626	>	castorRecHit.e[nhits] = rh.energy();
627	>	castorRecHit.iphi[nhits] = castorid.sector();
628	>	castorRecHit.depth[nhits] = castorid.module();
629	>	castorRecHit.phi[nhits] = getPhi(castorid);
630	>	}
631	>
632	>	nhits++;
633	>
634	>	} // end loop castor rechits
635	>	}
636	>
637	>	castorRecHit.n = nhits;
638	>	castorTree->Fill();
639	>	}
640	>
641	>	if(doZDCRecHit_){
642	>
643	>	edm::Handle<ZDCRecHitCollection> zdcrechits;
644	>
645	>	try{ ev.getByLabel("zdcreco",zdcrechits); }
646	>	catch(...) { edm::LogWarning(" ZDC ") << " Cannot get ZDC RecHits " << std::endl; }
647	>
648	>	int nhits = 0;
649	>
650	>	if (zdcrechits.failedToGet()!=0 || !zdcrechits.isValid()) {
651	>	edm::LogWarning(" ZDC ") << " Cannot read ZDCRecHitCollection" << std::endl;
652	>	} else {
653	>	for(size_t i1 = 0; i1 < zdcrechits->size(); ++i1) {
654	>	const ZDCRecHit & rh = (*zdcrechits)[i1];
655	>	HcalZDCDetId zdcid = rh.id();
656	>	if (nhits  < 18) {
657	>	zdcRecHit.e[nhits] = rh.energy();
658	>	zdcRecHit.zside[nhits] = zdcid.zside();
659	>	zdcRecHit.section[nhits] = zdcid.section();
660	>	zdcRecHit.channel[nhits] = zdcid.channel();
661	>	}
662	>
663	>	nhits++;
664	>
665	>	} // end loop zdc rechits
666	>	}
667	>
668	>	zdcRecHit.n = nhits;
669	>	zdcRecHitTree->Fill();
670	>
671	>	}
672	>
673	>	if(doZDCDigi_){
674	>
675	>	edm::Handle<ZDCDigiCollection> zdcdigis;
676	>
677	>	try{ ev.getByLabel("hcalDigis",zdcdigis); }
678	>	catch(...) { edm::LogWarning(" ZDC ") << " Cannot get ZDC Digis " << std::endl; }
679	>
680	>	int nhits = 0;
681	>
682	>	if (zdcdigis.failedToGet()!=0 || !zdcdigis.isValid()) {
683	>	edm::LogWarning(" ZDC ") << " Cannot read ZDCDigiCollection" << std::endl;
684	>	} else {
685	>	for(size_t i1 = 0; i1 < zdcdigis->size(); ++i1) {
686	>
687	>	const ZDCDataFrame & rh = (*zdcdigis)[i1];
688	>	HcalZDCDetId zdcid = rh.id();
689	>
690	>	if (nhits  < 18) {
691	>
692	>	int ts = 0;
693	>
694	>	for(int j1 = 0; j1 < rh.size(); j1++){
695	>
696	>	zdcDigi.chargefC[nhits][ts]= rh[ts].nominal_fC();
697	>	zdcDigi.adc[nhits][ts]= rh[ts].adc();
698	>	zdcDigi.zside[nhits][ts] = zdcid.zside();
699	>	zdcDigi.section[nhits][ts] = zdcid.section();
700	>	zdcDigi.channel[nhits][ts] = zdcid.channel();
701	>	zdcDigi.ts[nhits][ts] = ts;
702	>
703	>	ts++;
704	>	}
705	>
706	>	zdcDigi.nts=ts;
707	>
708	>	}
709	>
710	>	nhits++;
711	>
712	>	} // end loop zdc rechits
713	>	}
714	>
715	>	zdcDigi.n = nhits;
716	>	zdcDigiTree->Fill();
717	>
718	>	}
719	>
720	>	if(!doEbyEonly_){
721	>	towerTree->Fill();
722	>
723	>	eeTree->Fill();
724	>	ebTree->Fill();
725	>
726	>	hbheTree->Fill();
727	>	hfTree->Fill();
728	>
729	>	if (doFastJet_) {
730	>	bkgTree->Fill();
731	>	}
732	>	}
733	>
734	>	nt->Fill(nHFtowerPlus,nHFtowerMinus,nHFlongPlus,nHFlongMinus,nHFshortPlus,nHFshortMinus);
735
736		}
737
#	Line 289 | Line 741 | void
741		RecHitTreeProducer::beginJob()
742		{
743
744	<	hbheTree = fs->make<TTree>("hbhe","");
744	>	hbheTree = fs->make<TTree>("hbhe",versionTag);
745		hbheTree->Branch("n",&hbheRecHit.n,"n/I");
746		hbheTree->Branch("e",hbheRecHit.e,"e[n]/F");
747		hbheTree->Branch("et",hbheRecHit.et,"et[n]/F");
748		hbheTree->Branch("eta",hbheRecHit.eta,"eta[n]/F");
749		hbheTree->Branch("phi",hbheRecHit.phi,"phi[n]/F");
750	<	hbheTree->Branch("isjet",hbheRecHit.isjet,"isjet[n]/I");
751	<
752	<	hfTree = fs->make<TTree>("hf","");
750	>	hbheTree->Branch("perp",hbheRecHit.perp,"perp[n]/F");
751	>
752	>	hbheTree->Branch("isjet",hbheRecHit.isjet,"isjet[n]/O");
753	>	hbheTree->Branch("depth",hfRecHit.depth,"depth[n]/I");
754	>
755	>	hfTree = fs->make<TTree>("hf",versionTag);
756		hfTree->Branch("n",&hfRecHit.n,"n/I");
757		hfTree->Branch("e",hfRecHit.e,"e[n]/F");
758		hfTree->Branch("et",hfRecHit.et,"et[n]/F");
759		hfTree->Branch("eta",hfRecHit.eta,"eta[n]/F");
760		hfTree->Branch("phi",hfRecHit.phi,"phi[n]/F");
761	<	hfTree->Branch("isjet",hfRecHit.isjet,"isjet[n]/I");
761	>	hfTree->Branch("perp",hfRecHit.perp,"perp[n]/F");
762	>	hfTree->Branch("depth",hfRecHit.depth,"depth[n]/I");
763	>	hfTree->Branch("isjet",hfRecHit.isjet,"isjet[n]/O");
764
765	<	eeTree = fs->make<TTree>("ee","");
765	>	eeTree = fs->make<TTree>("ee",versionTag);
766		eeTree->Branch("n",&eeRecHit.n,"n/I");
767		eeTree->Branch("e",eeRecHit.e,"e[n]/F");
768		eeTree->Branch("et",eeRecHit.et,"et[n]/F");
769		eeTree->Branch("eta",eeRecHit.eta,"eta[n]/F");
770		eeTree->Branch("phi",eeRecHit.phi,"phi[n]/F");
771	<	eeTree->Branch("isjet",eeRecHit.isjet,"isjet[n]/I");
771	>	eeTree->Branch("perp",eeRecHit.perp,"perp[n]/F");
772	>
773	>	eeTree->Branch("isjet",eeRecHit.isjet,"isjet[n]/O");
774
775	<	ebTree = fs->make<TTree>("eb","");
775	>	ebTree = fs->make<TTree>("eb",versionTag);
776		ebTree->Branch("n",&ebRecHit.n,"n/I");
777		ebTree->Branch("e",ebRecHit.e,"e[n]/F");
778		ebTree->Branch("et",ebRecHit.et,"et[n]/F");
779		ebTree->Branch("eta",ebRecHit.eta,"eta[n]/F");
780		ebTree->Branch("phi",ebRecHit.phi,"phi[n]/F");
781	+	ebTree->Branch("perp",ebRecHit.perp,"perp[n]/F");
782	+
783		ebTree->Branch("isjet",ebRecHit.isjet,"isjet[n]/O");
784
785	+	towerTree = fs->make<TTree>("tower",versionTag);
786	+	towerTree->Branch("n",&myTowers.n,"n/I");
787	+	towerTree->Branch("e",myTowers.e,"e[n]/F");
788	+	towerTree->Branch("et",myTowers.et,"et[n]/F");
789	+	towerTree->Branch("eta",myTowers.eta,"eta[n]/F");
790	+	towerTree->Branch("phi",myTowers.phi,"phi[n]/F");
791	+	towerTree->Branch("isjet",myTowers.isjet,"isjet[n]/O");
792	+	towerTree->Branch("emEt",myTowers.emEt,"emEt[n]/F");
793	+	towerTree->Branch("hadEt",myTowers.hadEt,"hadEt[n]/F");
794	+
795	+
796	+	if(doCastor_){
797	+	castorTree = fs->make<TTree>("castor",versionTag);
798	+	castorTree->Branch("n",&castorRecHit.n,"n/I");
799	+	castorTree->Branch("e",castorRecHit.e,"e[n]/F");
800	+	castorTree->Branch("iphi",castorRecHit.iphi,"iphi[n]/I");
801	+	castorTree->Branch("phi",castorRecHit.phi,"phi[n]/F");
802	+	castorTree->Branch("depth",castorRecHit.depth,"depth[n]/I");
803	+	}
804	+
805	+	if(doZDCRecHit_){
806	+	zdcRecHitTree = fs->make<TTree>("zdcrechit",versionTag);
807	+	zdcRecHitTree->Branch("n",&zdcRecHit.n,"n/I");
808	+	zdcRecHitTree->Branch("e",zdcRecHit.e,"e[n]/F");
809	+	zdcRecHitTree->Branch("zside",zdcRecHit.zside,"zside[n]/I");
810	+	zdcRecHitTree->Branch("section",zdcRecHit.section,"section[n]/I");
811	+	zdcRecHitTree->Branch("channel",zdcRecHit.channel,"channel[n]/I");
812	+	}
813	+
814	+	if(doZDCDigi_){
815	+	zdcDigiTree = fs->make<TTree>("zdcrechit",versionTag);
816	+	zdcDigiTree->Branch("n",&zdcDigi.n,"n/I");
817	+	zdcDigiTree->Branch("nts",&zdcDigi.nts,"nts/I");
818	+	zdcDigiTree->Branch("zside",zdcDigi.zside,"zside[n][nts]/I");
819	+	zdcDigiTree->Branch("section",zdcDigi.section,"section[n][nts]/I");
820	+	zdcDigiTree->Branch("channel",zdcDigi.channel,"channel[n][nts]/I");
821	+	zdcDigiTree->Branch("adc",zdcDigi.adc,"adc[n][nts]/I");
822	+	zdcDigiTree->Branch("chargefC",zdcDigi.chargefC,"chargefC[n][nts]/F");
823	+	}
824	+
825	+	if (saveBothVtx_) {
826	+	towerTree->Branch("etVtx",myTowers.etVtx,"etvtx[n]/F");
827	+	towerTree->Branch("etaVtx",myTowers.etaVtx,"etavtx[n]/F");
828	+	towerTree->Branch("emEtVtx",myTowers.emEtVtx,"emEtVtx[n]/F");
829	+	towerTree->Branch("hadEtVtx",myTowers.hadEtVtx,"hadEtVtx[n]/F");
830	+	}
831	+
832	+	if(doBasicClusters_){
833	+	bcTree = fs->make<TTree>("bc",versionTag);
834	+	bcTree->Branch("n",&myBC.n,"n/I");
835	+	bcTree->Branch("e",myBC.e,"e[n]/F");
836	+	bcTree->Branch("et",myBC.et,"et[n]/F");
837	+	bcTree->Branch("eta",myBC.eta,"eta[n]/F");
838	+	bcTree->Branch("phi",myBC.phi,"phi[n]/F");
839	+	bcTree->Branch("jtpt",&myBC.jtpt,"jtpt/F");
840	+	bcTree->Branch("jteta",&myBC.jteta,"jteta/F");
841	+	bcTree->Branch("jtphi",&myBC.jtphi,"jtphi/F");
842	+	//     bcTree->Branch("isjet",bcRecHit.isjet,"isjet[n]/O");
843	+	}
844	+
845	+	if(doFastJet_){
846	+	bkgTree = fs->make<TTree>("bkg",versionTag);
847	+	bkgTree->Branch("n",&bkg.n,"n/I");
848	+	bkgTree->Branch("rho",bkg.rho,"rho[n]/F");
849	+	bkgTree->Branch("sigma",bkg.sigma,"sigma[n]/F");
850	+	}
851	+
852	+	nt = fs->make<TNtuple>("ntEvent","","nHFplus:nHFminus:nHFlongPlus:nHFlongMinus:nHFshortPlus:nHFshortMinus");
853	+
854		}
855
856		// ------------ method called once each job just after ending the event loop  ------------
#	Line 328 | Line 858 | void
858		RecHitTreeProducer::endJob() {
859		}
860
861	<	double RecHitTreeProducer::getEt(const DetId &id, double energy){
861	>	math::XYZPoint RecHitTreeProducer::getPosition(const DetId &id, reco::Vertex::Point vtx){
862	>	const GlobalPoint& pos=geo->getPosition(id);
863	>	math::XYZPoint posV(pos.x() - vtx.x(),pos.y() - vtx.y(),pos.z() - vtx.z());
864	>	return posV;
865	>	}
866	>
867	>	double RecHitTreeProducer::getEt(math::XYZPoint pos, double energy){
868	>	double et = energy*sin(pos.theta());
869	>	return et;
870	>	}
871	>
872	>	double RecHitTreeProducer::getEt(const DetId &id, double energy, reco::Vertex::Point vtx){
873		const GlobalPoint& pos=geo->getPosition(id);
874		double et = energy*sin(pos.theta());
875		return et;
876		}
877
878	<	double RecHitTreeProducer::getEta(const DetId &id){
878	>	double RecHitTreeProducer::getEta(const DetId &id, reco::Vertex::Point vtx){
879		const GlobalPoint& pos=geo->getPosition(id);
880		double et = pos.eta();
881		return et;
882		}
883
884	<	double RecHitTreeProducer::getPhi(const DetId &id){
884	>	double RecHitTreeProducer::getPhi(const DetId &id, reco::Vertex::Point vtx){
885		const GlobalPoint& pos=geo->getPosition(id);
886		double et = pos.phi();
887		return et;
888		}
889
890	+	double RecHitTreeProducer::getPerp(const DetId &id, reco::Vertex::Point vtx){
891	+	const GlobalPoint& pos=geo->getPosition(id);
892	+	double et = pos.perp();
893	+	return et;
894	+	}
895
896	+	reco::Vertex::Point RecHitTreeProducer::getVtx(const edm::Event& ev)
897	+	{
898	+	ev.getByLabel(VtxSrc_,vtxs);
899	+	int greatestvtx = 0;
900	+	int nVertex = vtxs->size();
901	+
902	+	for (unsigned int i = 0 ; i< vtxs->size(); ++i){
903	+	unsigned int daughter = (*vtxs)[i].tracksSize();
904	+	if( daughter > (*vtxs)[greatestvtx].tracksSize()) greatestvtx = i;
905	+	//cout <<"Vertex: "<< (*vtxs)[i].position().z()<<" "<<daughter<<endl;
906	+	}
907	+
908	+	if(nVertex<=0){
909	+	return reco::Vertex::Point(0,0,0);
910	+	}
911	+	return (*vtxs)[greatestvtx].position();
912	+	}
913
914		//define this as a plug-in
915		DEFINE_FWK_MODULE(RecHitTreeProducer);

Diff Legend

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