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

Diff Legend

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