ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

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

Comparing UserCode/CmsHi/JetAnalysis/src/RecHitTreeProducer.cc (file contents):
Revision 1.1 by yilmaz, Sat Oct 16 13:46:13 2010 UTC vs.
Revision 1.14 by frankma, Thu Oct 13 10:40:53 2011 UTC

#	Line 12 | Line 12
12		*/
13		//
14		// Original Author:  Yetkin Yilmaz
15	+	// Modified: Frank Ma
16		//         Created:  Tue Sep  7 11:38:19 EDT 2010
17		// $Id$
18		//
#	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 100000
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	+	bool isjet[MAXHITS];
83	+	float etvtx[MAXHITS];
84	+	float etavtx[MAXHITS];
85	+	float emEtVtx[MAXHITS];
86	+	float hadEtVtx[MAXHITS];
87	+
88	+	float jtpt;
89	+	float jteta;
90	+	float jtphi;
91
92		};
93
94
95	+	struct MyBkg{
96	+	int n;
97	+	float rho[50];
98	+	float sigma[50];
99	+	};
100	+
101
102		//
103		// class declaration
#	Line 90 | Line 110 | class RecHitTreeProducer : public edm::E
110		double       getEt(const DetId &id, double energy);
111		double       getEta(const DetId &id);
112		double       getPhi(const DetId &id);
113	+	reco::Vertex::Point getVtx(const edm::Event& ev);
114	+
115
116
117		private:
#	Line 99 | Line 121 | class RecHitTreeProducer : public edm::E
121
122		// ----------member data ---------------------------
123
102	–
124		edm::Handle<reco::Centrality> cent;
125		edm::Handle<vector<double> > ktRhos;
126		edm::Handle<vector<double> > akRhos;
127
128	<	edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > jets1;
129	<	edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > jets2;
109	<
110	<	//   typedef edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > >::const_iterator EcalIterator;
111	<	typedef vector<EcalRecHit>::const_iterator EcalIterator;
128	>	edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > ebHits;
129	>	edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > eeHits;
130
131	<	edm::Handle<reco::CaloJetCollection> signalJets;
131	>	edm::Handle<HFRecHitCollection> hfHits;
132	>	edm::Handle<HBHERecHitCollection> hbheHits;
133
134	+	edm::Handle<reco::BasicClusterCollection> bClusters;
135	+	edm::Handle<CaloTowerCollection> towers;
136	+	edm::Handle<reco::VertexCollection> vtxs;
137	+
138	+	typedef vector<EcalRecHit>::const_iterator EcalIterator;
139	+
140	+	edm::Handle<reco::CaloJetCollection> jets;
141	+
142	+	edm::Handle<std::vector<double> > rhos;
143	+	edm::Handle<std::vector<double> > sigmas;
144	+
145		MyRecHit hbheRecHit;
146		MyRecHit hfRecHit;
147	+	MyRecHit ebRecHit;
148	+	MyRecHit eeRecHit;
149	+	MyRecHit myBC;
150	+	MyRecHit myTowers;
151	+	MyBkg bkg;
152
153	+	TNtuple* nt;
154		TTree* hbheTree;
155		TTree* hfTree;
156	<
157	<	double cone;
158	<
156	>	TTree* ebTree;
157	>	TTree* eeTree;
158	>	TTree* bcTree;
159	>	TTree* towerTree;
160	>	TTree* bkgTree;
161	>
162	>	double cone;
163	>	double hfTowerThreshold_;
164	>	double hfLongThreshold_;
165	>	double hfShortThreshold_;
166	>	double hbheThreshold_;
167	>	double ebThreshold_;
168	>	double eeThreshold_;
169	>
170	>	double hbhePtMin_;
171	>	double hfPtMin_;
172	>	double ebPtMin_;
173	>	double eePtMin_;
174	>	double towerPtMin_;
175	>
176		edm::Service<TFileService> fs;
177		const CentralityBins * cbins_;
178		const CaloGeometry *geo;
179
180		edm::InputTag HcalRecHitHFSrc_;
181		edm::InputTag HcalRecHitHBHESrc_;
182	+	edm::InputTag EBSrc_;
183	+	edm::InputTag EESrc_;
184	+	edm::InputTag BCSrc_;
185	+	edm::InputTag TowerSrc_;
186	+	edm::InputTag VtxSrc_;
187	+
188	+	edm::InputTag JetSrc_;
189	+
190	+	edm::InputTag FastJetTag_;
191
192	+	bool useJets_;
193	+	bool doBasicClusters_;
194	+	bool doTowers_;
195	+	bool doEcal_;
196	+	bool doHcal_;
197
198	+	bool doFastJet_;
199	+
200	+	bool doEbyEonly_;
201		};
202
203		//
#	Line 147 | Line 217 | RecHitTreeProducer::RecHitTreeProducer(c
217		cone(0.5)
218		{
219		//now do what ever initialization is needed
220	+	EBSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("EBRecHitSrc",edm::InputTag("ecalRecHit","EcalRecHitsEB"));
221	+	EESrc_ = iConfig.getUntrackedParameter<edm::InputTag>("EERecHitSrc",edm::InputTag("ecalRecHit","EcalRecHitsEE"));
222		HcalRecHitHFSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHFRecHitSrc",edm::InputTag("hfreco"));
223		HcalRecHitHBHESrc_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHBHERecHitSrc",edm::InputTag("hbhereco"));
224	+	BCSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("BasicClusterSrc1",edm::InputTag("ecalRecHit","EcalRecHitsEB","RECO"));
225	+	TowerSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("towersSrc",edm::InputTag("towerMaker"));
226	+	VtxSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("vtxSrc",edm::InputTag("hiSelectedVertex"));
227	+	JetSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("JetSrc",edm::InputTag("iterativeConePu5CaloJets"));
228	+	useJets_ = iConfig.getUntrackedParameter<bool>("useJets",true);
229	+	doBasicClusters_ = iConfig.getUntrackedParameter<bool>("doBasicClusters",false);
230	+	doTowers_ = iConfig.getUntrackedParameter<bool>("doTowers",true);
231	+	doEcal_ = iConfig.getUntrackedParameter<bool>("doEcal",true);
232	+	doHcal_ = iConfig.getUntrackedParameter<bool>("doHcal",true);
233	+	doFastJet_ = iConfig.getUntrackedParameter<bool>("doFastJet",true);
234	+	FastJetTag_ = iConfig.getUntrackedParameter<edm::InputTag>("FastJetTag",edm::InputTag("kt4CaloJets"));
235	+	doEbyEonly_ = iConfig.getUntrackedParameter<bool>("doEbyEonly",false);
236	+	hfTowerThreshold_ = iConfig.getUntrackedParameter<double>("HFtowerMin",3.);
237	+	hfLongThreshold_ = iConfig.getUntrackedParameter<double>("HFlongMin",0.5);
238	+	hfShortThreshold_ = iConfig.getUntrackedParameter<double>("HFshortMin",0.85);
239	+	hbhePtMin_ = iConfig.getUntrackedParameter<double>("HBHETreePtMin",0);
240	+	hfPtMin_ = iConfig.getUntrackedParameter<double>("HFTreePtMin",0);
241	+	ebPtMin_ = iConfig.getUntrackedParameter<double>("EBTreePtMin",0);
242	+	eePtMin_ = iConfig.getUntrackedParameter<double>("EETreePtMin",0.);
243	+	towerPtMin_ = iConfig.getUntrackedParameter<double>("TowerTreePtMin",0.);
244		}
245
246
#	Line 172 | Line 264 | RecHitTreeProducer::analyze(const edm::E
264
265		hfRecHit.n = 0;
266		hbheRecHit.n = 0;
267	<
268	<	edm::Handle<HFRecHitCollection> hfHits;
269	<	edm::Handle<HBHERecHitCollection> hbheHits;
270	<
267	>	ebRecHit.n = 0;
268	>	eeRecHit.n = 0;
269	>	myBC.n = 0;
270	>	myTowers.n = 0;
271	>	bkg.n = 0;
272	>
273	>	// get vertex
274	>	reco::Vertex::Point vtx = getVtx(ev);
275	>
276	>	if(doEcal_){
277	>	ev.getByLabel(EBSrc_,ebHits);
278	>	ev.getByLabel(EESrc_,eeHits);
279	>	}
280	>	if(doHcal_){
281		ev.getByLabel(HcalRecHitHFSrc_,hfHits);
282		ev.getByLabel(HcalRecHitHBHESrc_,hbheHits);
283	<
284	<
285	<	if(0 && !cbins_) cbins_ = getCentralityBinsFromDB(iSetup);
283	>	}
284	>	if(useJets_) {
285	>	ev.getByLabel(JetSrc_,jets);
286	>	}
287	>
288	>	if(doBasicClusters_){
289	>	ev.getByLabel(BCSrc_,bClusters);
290	>	}
291	>
292	>	if(doTowers_){
293	>	ev.getByLabel(TowerSrc_,towers);
294	>	}
295	>
296	>	if(doFastJet_){
297	>	ev.getByLabel(edm::InputTag(FastJetTag_.label(),"rhos",FastJetTag_.process()),rhos);
298	>	ev.getByLabel(edm::InputTag(FastJetTag_.label(),"sigmas",FastJetTag_.process()),sigmas);
299	>	bkg.n = rhos->size();
300	>	for(unsigned int i = 0; i < rhos->size(); ++i){
301	>	bkg.rho[i] = (*rhos)[i];
302	>	bkg.sigma[i] = (*sigmas)[i];
303	>	}
304	>	}
305	>
306	>	if(0 && !cbins_) cbins_ = getCentralityBinsFromDB(iSetup);
307
308		if(!geo){
309		edm::ESHandle<CaloGeometry> pGeo;
#	Line 188 | Line 311 | RecHitTreeProducer::analyze(const edm::E
311		geo = pGeo.product();
312		}
313
314	+	int nHFlongPlus = 0;
315	+	int nHFshortPlus = 0;
316	+	int nHFtowerPlus = 0;
317	+	int nHFlongMinus = 0;
318	+	int nHFshortMinus = 0;
319	+	int nHFtowerMinus = 0;
320	+
321	+
322	+
323	+	if(doHcal_){
324		for(unsigned int i = 0; i < hfHits->size(); ++i){
325		const HFRecHit & hit= (*hfHits)[i];
326		hfRecHit.e[hfRecHit.n] = hit.energy();
327		hfRecHit.et[hfRecHit.n] = getEt(hit.id(),hit.energy());
328		hfRecHit.eta[hfRecHit.n] = getEta(hit.id());
329		hfRecHit.phi[hfRecHit.n] = getPhi(hit.id());
330	<	//     hfRecHit.ieta[hfRecHit.n] = hit.id().ieta();
331	<	//     hfRecHit.iphi[hfRecHit.n] = hit.id().iphi();
199	<
200	<	hfRecHit.n++;
201	<	}
330	>	hfRecHit.isjet[hfRecHit.n] = false;
331	>	hfRecHit.depth[hfRecHit.n] = hit.id().depth();
332
333	+	if(hit.id().ieta() > 0){
334	+	if(hit.energy() > hfShortThreshold_ && hit.id().depth() != 1) nHFshortPlus++;
335	+	if(hit.energy() > hfLongThreshold_ && hit.id().depth() == 1) nHFlongPlus++;
336	+	}else{
337	+	if(hit.energy() > hfShortThreshold_ && hit.id().depth() != 1) nHFshortMinus++;
338	+	if(hit.energy() > hfLongThreshold_ && hit.id().depth() == 1) nHFlongMinus++;
339	+	}
340	+
341	+	if(useJets_){
342	+	for(unsigned int j = 0 ; j < jets->size(); ++j){
343	+	const reco::Jet& jet = (*jets)[j];
344	+	double dr = reco::deltaR(hfRecHit.eta[hfRecHit.n],hfRecHit.phi[hfRecHit.n],jet.eta(),jet.phi());
345	+	if(dr < cone){ hfRecHit.isjet[hfRecHit.n] = true; }
346	+	}
347	+	}
348	+	if (hfRecHit.et[hfRecHit.n]>=hfPtMin_) hfRecHit.n++;
349	+	}
350	+	if(!doEbyEonly_){
351		for(unsigned int i = 0; i < hbheHits->size(); ++i){
352		const HBHERecHit & hit= (*hbheHits)[i];
353	<	hbheRecHit.e[hfRecHit.n] = hit.energy();
354	<	hbheRecHit.et[hfRecHit.n] = getEt(hit.id(),hit.energy());
355	<	hbheRecHit.eta[hfRecHit.n] = getEta(hit.id());
356	<	hbheRecHit.phi[hfRecHit.n] = getPhi(hit.id());
357	<	//     hbheRecHit.ieta[hfRecHit.n] = hit.id().ieta();
358	<	//     hbheRecHit.iphi[hfRecHit.n] = hit.id().iphi();
359	<
353	>	if (getEt(hit.id(),hit.energy())<hbhePtMin_) continue;
354	>	hbheRecHit.e[hbheRecHit.n] = hit.energy();
355	>	hbheRecHit.et[hbheRecHit.n] = getEt(hit.id(),hit.energy());
356	>	hbheRecHit.eta[hbheRecHit.n] = getEta(hit.id());
357	>	hbheRecHit.phi[hbheRecHit.n] = getPhi(hit.id());
358	>	hbheRecHit.isjet[hbheRecHit.n] = false;
359	>	if(useJets_){
360	>	for(unsigned int j = 0 ; j < jets->size(); ++j){
361	>	const reco::Jet& jet = (*jets)[j];
362	>	double dr = reco::deltaR(hbheRecHit.eta[hbheRecHit.n],hbheRecHit.phi[hbheRecHit.n],jet.eta(),jet.phi());
363	>	if(dr < cone){ hbheRecHit.isjet[hbheRecHit.n] = true; }
364	>	}
365	>	}
366		hbheRecHit.n++;
367		}
368	+	}
369	+	}
370	+	if(doEcal_ && !doEbyEonly_){
371	+	for(unsigned int i = 0; i < ebHits->size(); ++i){
372	+	const EcalRecHit & hit= (*ebHits)[i];
373	+	if (getEt(hit.id(),hit.energy())<ebPtMin_) continue;
374	+	ebRecHit.e[ebRecHit.n] = hit.energy();
375	+	ebRecHit.et[ebRecHit.n] = getEt(hit.id(),hit.energy());
376	+	ebRecHit.eta[ebRecHit.n] = getEta(hit.id());
377	+	ebRecHit.phi[ebRecHit.n] = getPhi(hit.id());
378	+	ebRecHit.isjet[ebRecHit.n] = false;
379	+	if(useJets_){
380	+	for(unsigned int j = 0 ; j < jets->size(); ++j){
381	+	const reco::Jet& jet = (*jets)[j];
382	+	double dr = reco::deltaR(ebRecHit.eta[ebRecHit.n],ebRecHit.phi[ebRecHit.n],jet.eta(),jet.phi());
383	+	if(dr < cone){ ebRecHit.isjet[ebRecHit.n] = true; }
384	+	}
385	+	}
386	+	ebRecHit.n++;
387	+	}
388	+
389	+	for(unsigned int i = 0; i < eeHits->size(); ++i){
390	+	const EcalRecHit & hit= (*eeHits)[i];
391	+	if (getEt(hit.id(),hit.energy())<eePtMin_) continue;
392	+	eeRecHit.e[eeRecHit.n] = hit.energy();
393	+	eeRecHit.et[eeRecHit.n] = getEt(hit.id(),hit.energy());
394	+	eeRecHit.eta[eeRecHit.n] = getEta(hit.id());
395	+	eeRecHit.phi[eeRecHit.n] = getPhi(hit.id());
396	+	eeRecHit.isjet[eeRecHit.n] = false;
397	+	if(useJets_){
398	+	for(unsigned int j = 0 ; j < jets->size(); ++j){
399	+	const reco::Jet& jet = (*jets)[j];
400	+	double dr = reco::deltaR(eeRecHit.eta[eeRecHit.n],eeRecHit.phi[eeRecHit.n],jet.eta(),jet.phi());
401	+	if(dr < cone){ eeRecHit.isjet[eeRecHit.n] = true; }
402	+	}
403	+	}
404	+	eeRecHit.n++;
405	+	}
406	+	}
407
408	<	hbheTree->Fill();
216	<	hfTree->Fill();
408	>	if(doTowers_){
409
410	+	for(unsigned int i = 0; i < towers->size(); ++i){
411	+	const CaloTower & hit= (*towers)[i];
412	+	if (getEt(hit.id(),hit.energy())<towerPtMin_) continue;
413	+	myTowers.e[myTowers.n] = hit.energy();
414	+	myTowers.et[myTowers.n] = getEt(hit.id(),hit.energy());
415	+	myTowers.eta[myTowers.n] = getEta(hit.id());
416	+	myTowers.phi[myTowers.n] = getPhi(hit.id());
417	+	myTowers.isjet[myTowers.n] = false;
418	+	myTowers.etvtx[myTowers.n] = hit.p4(vtx).Et();
419	+	myTowers.etavtx[myTowers.n] = hit.p4(vtx).Eta();
420	+	myTowers.emEtVtx[myTowers.n] = hit.emEt(vtx);
421	+	myTowers.hadEtVtx[myTowers.n] = hit.hadEt(vtx);
422	+
423	+	if(hit.ieta() > 29 && hit.energy() > hfTowerThreshold_) nHFtowerPlus++;
424	+	if(hit.ieta() < -29 && hit.energy() > hfTowerThreshold_) nHFtowerMinus++;
425	+
426	+	if(useJets_){
427	+	for(unsigned int j = 0 ; j < jets->size(); ++j){
428	+	const reco::Jet& jet = (*jets)[j];
429	+	double dr = reco::deltaR(myTowers.eta[myTowers.n],myTowers.phi[myTowers.n],jet.eta(),jet.phi());
430	+	if(dr < cone){ myTowers.isjet[myTowers.n] = true; }
431	+	}
432	+	}
433	+	myTowers.n++;
434	+	}
435	+
436	+	}
437	+
438	+	if(doBasicClusters_ && !doEbyEonly_){
439	+	for(unsigned int j = 0 ; j < jets->size(); ++j){
440	+	const reco::Jet& jet = (*jets)[j];
441	+	myBC.n = 0;
442	+	myBC.jtpt = jet.pt();
443	+	myBC.jteta = jet.eta();
444	+	myBC.jtphi = jet.phi();
445	+
446	+	for(unsigned int i = 0; i < bClusters->size(); ++i){
447	+	const reco::BasicCluster & bc= (*bClusters)[i];
448	+	double dr = reco::deltaR(bc.eta(),bc.phi(),jet.eta(),jet.phi());
449	+	if(dr < cone){
450	+	myBC.e[myBC.n] = bc.energy();
451	+	myBC.et[myBC.n] = bc.energy()*sin(bc.position().theta());
452	+	myBC.eta[myBC.n] = bc.eta();
453	+	myBC.phi[myBC.n] = bc.phi();
454	+	myBC.n++;
455	+	}
456	+	}
457	+	bcTree->Fill();
458	+	}
459	+	}
460	+
461	+	if(!doEbyEonly_){
462	+	towerTree->Fill();
463	+
464	+	eeTree->Fill();
465	+	ebTree->Fill();
466	+
467	+	hbheTree->Fill();
468	+	hfTree->Fill();
469	+
470	+	if (doFastJet_) {
471	+	bkgTree->Fill();
472	+	}
473	+	}
474	+
475	+	nt->Fill(nHFtowerPlus,nHFtowerMinus,nHFlongPlus,nHFlongMinus,nHFshortPlus,nHFshortMinus);
476	+
477		}
478
479
#	Line 222 | Line 481 | RecHitTreeProducer::analyze(const edm::E
481		void
482		RecHitTreeProducer::beginJob()
483		{
484	<
484	>
485		hbheTree = fs->make<TTree>("hbhe","");
486		hbheTree->Branch("n",&hbheRecHit.n,"n/I");
487		hbheTree->Branch("e",hbheRecHit.e,"e[n]/F");
488		hbheTree->Branch("et",hbheRecHit.et,"et[n]/F");
489		hbheTree->Branch("eta",hbheRecHit.eta,"eta[n]/F");
490		hbheTree->Branch("phi",hbheRecHit.phi,"phi[n]/F");
491	<
491	>	hbheTree->Branch("isjet",hbheRecHit.isjet,"isjet[n]/O");
492	>
493		hfTree = fs->make<TTree>("hf","");
494		hfTree->Branch("n",&hfRecHit.n,"n/I");
495		hfTree->Branch("e",hfRecHit.e,"e[n]/F");
496		hfTree->Branch("et",hfRecHit.et,"et[n]/F");
497		hfTree->Branch("eta",hfRecHit.eta,"eta[n]/F");
498		hfTree->Branch("phi",hfRecHit.phi,"phi[n]/F");
499	+	hfTree->Branch("depth",hfRecHit.depth,"depth[n]/I");
500	+	hfTree->Branch("isjet",hfRecHit.isjet,"isjet[n]/O");
501	+
502	+	eeTree = fs->make<TTree>("ee","");
503	+	eeTree->Branch("n",&eeRecHit.n,"n/I");
504	+	eeTree->Branch("e",eeRecHit.e,"e[n]/F");
505	+	eeTree->Branch("et",eeRecHit.et,"et[n]/F");
506	+	eeTree->Branch("eta",eeRecHit.eta,"eta[n]/F");
507	+	eeTree->Branch("phi",eeRecHit.phi,"phi[n]/F");
508	+	eeTree->Branch("isjet",eeRecHit.isjet,"isjet[n]/O");
509	+
510	+	ebTree = fs->make<TTree>("eb","");
511	+	ebTree->Branch("n",&ebRecHit.n,"n/I");
512	+	ebTree->Branch("e",ebRecHit.e,"e[n]/F");
513	+	ebTree->Branch("et",ebRecHit.et,"et[n]/F");
514	+	ebTree->Branch("eta",ebRecHit.eta,"eta[n]/F");
515	+	ebTree->Branch("phi",ebRecHit.phi,"phi[n]/F");
516	+	ebTree->Branch("isjet",ebRecHit.isjet,"isjet[n]/O");
517	+
518	+	towerTree = fs->make<TTree>("tower","");
519	+	towerTree->Branch("n",&myTowers.n,"n/I");
520	+	towerTree->Branch("e",myTowers.e,"e[n]/F");
521	+	towerTree->Branch("et",myTowers.et,"et[n]/F");
522	+	towerTree->Branch("eta",myTowers.eta,"eta[n]/F");
523	+	towerTree->Branch("phi",myTowers.phi,"phi[n]/F");
524	+	towerTree->Branch("isjet",myTowers.isjet,"isjet[n]/O");
525	+	towerTree->Branch("etvtx",myTowers.etvtx,"etvtx[n]/F");
526	+	towerTree->Branch("etavtx",myTowers.etavtx,"etavtx[n]/F");
527	+	towerTree->Branch("emEtVtx",myTowers.emEtVtx,"emEtVtx[n]/F");
528	+	towerTree->Branch("hadEtVtx",myTowers.hadEtVtx,"hadEtVtx[n]/F");
529	+
530	+
531	+	if(doBasicClusters_){
532	+	bcTree = fs->make<TTree>("bc","");
533	+	bcTree->Branch("n",&myBC.n,"n/I");
534	+	bcTree->Branch("e",myBC.e,"e[n]/F");
535	+	bcTree->Branch("et",myBC.et,"et[n]/F");
536	+	bcTree->Branch("eta",myBC.eta,"eta[n]/F");
537	+	bcTree->Branch("phi",myBC.phi,"phi[n]/F");
538	+	bcTree->Branch("jtpt",&myBC.jtpt,"jtpt/F");
539	+	bcTree->Branch("jteta",&myBC.jteta,"jteta/F");
540	+	bcTree->Branch("jtphi",&myBC.jtphi,"jtphi/F");
541	+	//     bcTree->Branch("isjet",bcRecHit.isjet,"isjet[n]/O");
542	+	}
543	+
544	+	if(doFastJet_){
545	+	bkgTree = fs->make<TTree>("bkg","");
546	+	bkgTree->Branch("n",&bkg.n,"n/I");
547	+	bkgTree->Branch("rho",bkg.rho,"rho[n]/F");
548	+	bkgTree->Branch("sigma",bkg.sigma,"sigma[n]/F");
549	+	}
550	+
551	+	nt = fs->make<TNtuple>("ntEvent","","nHFplus:nHFminus:nHFlongPlus:nHFlongMinus:nHFshortPlus:nHFshortMinus");
552
553		}
554
#	Line 249 | Line 562 | double RecHitTreeProducer::getEt(const D
562		double et = energy*sin(pos.theta());
563		return et;
564		}
565	<
565	>
566		double RecHitTreeProducer::getEta(const DetId &id){
567		const GlobalPoint& pos=geo->getPosition(id);
568		double et = pos.eta();
#	Line 262 | Line 575 | double RecHitTreeProducer::getPhi(const
575		return et;
576		}
577
578	<
578	>	reco::Vertex::Point RecHitTreeProducer::getVtx(const edm::Event& ev)
579	>	{
580	>	ev.getByLabel(VtxSrc_,vtxs);
581	>	int greatestvtx = 0;
582	>	int nVertex = vtxs->size();
583	>
584	>	for (unsigned int i = 0 ; i< vtxs->size(); ++i){
585	>	unsigned int daughter = (*vtxs)[i].tracksSize();
586	>	if( daughter > (*vtxs)[greatestvtx].tracksSize()) greatestvtx = i;
587	>	//cout <<"Vertex: "<< (*vtxs)[i].position().z()<<" "<<daughter<<endl;
588	>	}
589	>
590	>	if(nVertex<=0){
591	>	return reco::Vertex::Point(-999,-999,-999);
592	>	}
593	>	return (*vtxs)[greatestvtx].position();
594	>	}
595
596		//define this as a plug-in
597		DEFINE_FWK_MODULE(RecHitTreeProducer);

Diff Legend

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