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.2 by yilmaz, Mon Oct 18 16:13:37 2010 UTC vs.
Revision 1.7 by yilmaz, Tue Nov 2 12:58:33 2010 UTC

#	Line 57 | Line 57
57		#include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
58		#include "DataFormats/HcalDetId/interface/HcalDetId.h"
59		#include "DataFormats/EcalDetId/interface/EcalDetIdCollections.h"
60	+	#include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
61		#include "DataFormats/DetId/interface/DetId.h"
62
63
#	Line 64 | Line 65
65
66		using namespace std;
67
68	<	#define MAXHITS 100000
68	>	#define MAXHITS 1000000
69
70		struct MyRecHit{
71
#	Line 74 | Line 75 | struct MyRecHit{
75		float et[MAXHITS];
76		float eta[MAXHITS];
77		float phi[MAXHITS];
78	+	bool isjet[MAXHITS];
79	+
80	+	float jtpt;
81	+	float jteta;
82	+	float jtphi;
83
84		};
85
86
87	+	struct MyBkg{
88	+	int n;
89	+	float rho[50];
90	+	float sigma[50];
91	+	};
92	+
93
94		//
95		// class declaration
#	Line 99 | Line 111 | class RecHitTreeProducer : public edm::E
111
112		// ----------member data ---------------------------
113
102	–
114		edm::Handle<reco::Centrality> cent;
115		edm::Handle<vector<double> > ktRhos;
116		edm::Handle<vector<double> > akRhos;
#	Line 110 | Line 121 | class RecHitTreeProducer : public edm::E
121		edm::Handle<HFRecHitCollection> hfHits;
122		edm::Handle<HBHERecHitCollection> hbheHits;
123
124	<	typedef vector<EcalRecHit>::const_iterator EcalIterator;
124	>	edm::Handle<reco::BasicClusterCollection> bClusters;
125	>	edm::Handle<CaloTowerCollection> towers;
126	>
127
128	<	edm::Handle<reco::CaloJetCollection> signalJets;
128	>	typedef vector<EcalRecHit>::const_iterator EcalIterator;
129	>
130	>	edm::Handle<reco::CaloJetCollection> jets;
131
132	+	edm::Handle<std::vector<double> > rhos;
133	+	edm::Handle<std::vector<double> > sigmas;
134	+
135		MyRecHit hbheRecHit;
136		MyRecHit hfRecHit;
137		MyRecHit ebRecHit;
138		MyRecHit eeRecHit;
139	+	MyRecHit myBC;
140	+	MyRecHit myTowers;
141	+	MyBkg bkg;
142	+
143
144		TTree* hbheTree;
145		TTree* hfTree;
146		TTree* ebTree;
147		TTree* eeTree;
148	+	TTree* bcTree;
149	+	TTree* towerTree;
150	+	TTree* bkgTree;
151	+
152		double cone;
153
154		edm::Service<TFileService> fs;
#	Line 133 | Line 159 | class RecHitTreeProducer : public edm::E
159		edm::InputTag HcalRecHitHBHESrc_;
160		edm::InputTag EBSrc_;
161		edm::InputTag EESrc_;
162	+	edm::InputTag BCSrc_;
163	+	edm::InputTag TowerSrc_;
164	+
165	+	edm::InputTag JetSrc_;
166	+
167	+	edm::InputTag FastJetTag_;
168
169	+	bool useJets_;
170	+	bool doBasicClusters_;
171	+	bool doTowers_;
172	+	bool doEcal_;
173	+	bool doHcal_;
174
175	+	bool doFastJet_;
176		};
177
178		//
#	Line 154 | Line 192 | RecHitTreeProducer::RecHitTreeProducer(c
192		cone(0.5)
193		{
194		//now do what ever initialization is needed
157	–
195		EBSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("EBRecHitSrc",edm::InputTag("ecalRecHit","EcalRecHitsEB"));
196		EESrc_ = iConfig.getUntrackedParameter<edm::InputTag>("EERecHitSrc",edm::InputTag("ecalRecHit","EcalRecHitsEE"));
197		HcalRecHitHFSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHFRecHitSrc",edm::InputTag("hfreco"));
198		HcalRecHitHBHESrc_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHBHERecHitSrc",edm::InputTag("hbhereco"));
199	+	BCSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("BasicClusterSrc1",edm::InputTag("ecalRecHit","EcalRecHitsEB","RECO"));
200	+	TowerSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("towersSrc",edm::InputTag("towerMaker"));
201	+	JetSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("JetSrc",edm::InputTag("iterativeCone5CaloJets"));
202	+	useJets_ = iConfig.getUntrackedParameter<bool>("useJets",true);
203	+	doBasicClusters_ = iConfig.getUntrackedParameter<bool>("doBasicClusters",false);
204	+	doTowers_ = iConfig.getUntrackedParameter<bool>("doTowers",true);
205	+	doEcal_ = iConfig.getUntrackedParameter<bool>("doEcal",true);
206	+	doHcal_ = iConfig.getUntrackedParameter<bool>("doHcal",true);
207	+	doFastJet_ = iConfig.getUntrackedParameter<bool>("doFastJet",true);
208	+	FastJetTag_ = iConfig.getUntrackedParameter<edm::InputTag>("FastJetTag",edm::InputTag("kt4CaloJets"));
209	+
210		}
211
212
#	Line 184 | Line 232 | RecHitTreeProducer::analyze(const edm::E
232		hbheRecHit.n = 0;
233		ebRecHit.n = 0;
234		eeRecHit.n = 0;
235	+	myBC.n = 0;
236	+	myTowers.n = 0;
237	+	bkg.n = 0;
238
239	+	if(doEcal_){
240		ev.getByLabel(EBSrc_,ebHits);
241		ev.getByLabel(EESrc_,eeHits);
242	<
242	>	}
243	>	if(doHcal_){
244		ev.getByLabel(HcalRecHitHFSrc_,hfHits);
245		ev.getByLabel(HcalRecHitHBHESrc_,hbheHits);
246	<
247	<	if(0 && !cbins_) cbins_ = getCentralityBinsFromDB(iSetup);
246	>	}
247	>	if(useJets_) {
248	>	ev.getByLabel(JetSrc_,jets);
249	>	}
250	>
251	>	if(doBasicClusters_){
252	>	ev.getByLabel(BCSrc_,bClusters);
253	>	}
254	>
255	>	if(doTowers_){
256	>	ev.getByLabel(TowerSrc_,towers);
257	>	}
258	>
259	>	if(doFastJet_){
260	>	ev.getByLabel(edm::InputTag(FastJetTag_.label(),"rhos",FastJetTag_.process()),rhos);
261	>	ev.getByLabel(edm::InputTag(FastJetTag_.label(),"sigmas",FastJetTag_.process()),sigmas);
262	>	bkg.n = rhos->size();
263	>	for(unsigned int i = 0; i < rhos->size(); ++i){
264	>	bkg.rho[i] = (*rhos)[i];
265	>	bkg.sigma[i] = (*sigmas)[i];
266	>	}
267	>	}
268	>
269	>	if(0 && !cbins_) cbins_ = getCentralityBinsFromDB(iSetup);
270
271		if(!geo){
272		edm::ESHandle<CaloGeometry> pGeo;
273		iSetup.get<CaloGeometryRecord>().get(pGeo);
274		geo = pGeo.product();
275		}
276	<
276	>
277	>	if(doHcal_){
278		for(unsigned int i = 0; i < hfHits->size(); ++i){
279		const HFRecHit & hit= (*hfHits)[i];
280		hfRecHit.e[hfRecHit.n] = hit.energy();
281		hfRecHit.et[hfRecHit.n] = getEt(hit.id(),hit.energy());
282		hfRecHit.eta[hfRecHit.n] = getEta(hit.id());
283		hfRecHit.phi[hfRecHit.n] = getPhi(hit.id());
284	+	hfRecHit.isjet[hfRecHit.n] = false;
285	+	if(useJets_){
286	+	for(unsigned int j = 0 ; j < jets->size(); ++j){
287	+	const reco::Jet& jet = (*jets)[j];
288	+	double dr = reco::deltaR(hfRecHit.eta[hfRecHit.n],hfRecHit.phi[hfRecHit.n],jet.eta(),jet.phi());
289	+	if(dr < cone){ hfRecHit.isjet[hfRecHit.n] = true; }
290	+	}
291	+	}
292		hfRecHit.n++;
293		}
294	<
294	>
295		for(unsigned int i = 0; i < hbheHits->size(); ++i){
296		const HBHERecHit & hit= (*hbheHits)[i];
297	<	hbheRecHit.e[hfRecHit.n] = hit.energy();
298	<	hbheRecHit.et[hfRecHit.n] = getEt(hit.id(),hit.energy());
299	<	hbheRecHit.eta[hfRecHit.n] = getEta(hit.id());
300	<	hbheRecHit.phi[hfRecHit.n] = getPhi(hit.id());
297	>	hbheRecHit.e[hbheRecHit.n] = hit.energy();
298	>	hbheRecHit.et[hbheRecHit.n] = getEt(hit.id(),hit.energy());
299	>	hbheRecHit.eta[hbheRecHit.n] = getEta(hit.id());
300	>	hbheRecHit.phi[hbheRecHit.n] = getPhi(hit.id());
301	>	hbheRecHit.isjet[hbheRecHit.n] = false;
302	>	if(useJets_){
303	>	for(unsigned int j = 0 ; j < jets->size(); ++j){
304	>	const reco::Jet& jet = (*jets)[j];
305	>	double dr = reco::deltaR(hbheRecHit.eta[hbheRecHit.n],hbheRecHit.phi[hbheRecHit.n],jet.eta(),jet.phi());
306	>	if(dr < cone){ hbheRecHit.isjet[hbheRecHit.n] = true; }
307	>	}
308	>	}
309		hbheRecHit.n++;
310		}
311	<
311	>	}
312	>	if(doEcal_){
313		for(unsigned int i = 0; i < ebHits->size(); ++i){
314		const EcalRecHit & hit= (*ebHits)[i];
315		ebRecHit.e[ebRecHit.n] = hit.energy();
316		ebRecHit.et[ebRecHit.n] = getEt(hit.id(),hit.energy());
317		ebRecHit.eta[ebRecHit.n] = getEta(hit.id());
318		ebRecHit.phi[ebRecHit.n] = getPhi(hit.id());
319	+	ebRecHit.isjet[ebRecHit.n] = false;
320	+	if(useJets_){
321	+	for(unsigned int j = 0 ; j < jets->size(); ++j){
322	+	const reco::Jet& jet = (*jets)[j];
323	+	double dr = reco::deltaR(ebRecHit.eta[ebRecHit.n],ebRecHit.phi[ebRecHit.n],jet.eta(),jet.phi());
324	+	if(dr < cone){ ebRecHit.isjet[ebRecHit.n] = true; }
325	+	}
326	+	}
327		ebRecHit.n++;
328		}
329	<
329	>
330		for(unsigned int i = 0; i < eeHits->size(); ++i){
331		const EcalRecHit & hit= (*eeHits)[i];
332		eeRecHit.e[eeRecHit.n] = hit.energy();
333		eeRecHit.et[eeRecHit.n] = getEt(hit.id(),hit.energy());
334		eeRecHit.eta[eeRecHit.n] = getEta(hit.id());
335		eeRecHit.phi[eeRecHit.n] = getPhi(hit.id());
336	+	eeRecHit.isjet[eeRecHit.n] = false;
337	+	if(useJets_){
338	+	for(unsigned int j = 0 ; j < jets->size(); ++j){
339	+	const reco::Jet& jet = (*jets)[j];
340	+	double dr = reco::deltaR(eeRecHit.eta[eeRecHit.n],eeRecHit.phi[eeRecHit.n],jet.eta(),jet.phi());
341	+	if(dr < cone){ eeRecHit.isjet[eeRecHit.n] = true; }
342	+	}
343	+	}
344		eeRecHit.n++;
345		}
346	+	}
347	+
348	+	if(doTowers_){
349
350	+	for(unsigned int i = 0; i < towers->size(); ++i){
351	+	const CaloTower & hit= (*towers)[i];
352	+	myTowers.e[myTowers.n] = hit.energy();
353	+	myTowers.et[myTowers.n] = getEt(hit.id(),hit.energy());
354	+	myTowers.eta[myTowers.n] = getEta(hit.id());
355	+	myTowers.phi[myTowers.n] = getPhi(hit.id());
356	+	myTowers.isjet[myTowers.n] = false;
357	+	if(useJets_){
358	+	for(unsigned int j = 0 ; j < jets->size(); ++j){
359	+	const reco::Jet& jet = (*jets)[j];
360	+	double dr = reco::deltaR(myTowers.eta[myTowers.n],myTowers.phi[myTowers.n],jet.eta(),jet.phi());
361	+	if(dr < cone){ myTowers.isjet[myTowers.n] = true; }
362	+	}
363	+	}
364	+	myTowers.n++;
365	+	}
366	+
367	+	}
368	+
369	+	if(doBasicClusters_){
370	+	for(unsigned int j = 0 ; j < jets->size(); ++j){
371	+	const reco::Jet& jet = (*jets)[j];
372	+	myBC.n = 0;
373	+	myBC.jtpt = jet.pt();
374	+	myBC.jteta = jet.eta();
375	+	myBC.jtphi = jet.phi();
376	+
377	+	for(unsigned int i = 0; i < bClusters->size(); ++i){
378	+	const reco::BasicCluster & bc= (*bClusters)[i];
379	+	double dr = reco::deltaR(bc.eta(),bc.phi(),jet.eta(),jet.phi());
380	+	if(dr < cone){
381	+	myBC.e[myBC.n] = bc.energy();
382	+	myBC.et[myBC.n] = bc.energy()*sin(bc.position().theta());
383	+	myBC.eta[myBC.n] = bc.eta();
384	+	myBC.phi[myBC.n] = bc.phi();
385	+	myBC.n++;
386	+	}
387	+	}
388	+	bcTree->Fill();
389	+	}
390	+	}
391	+
392	+	towerTree->Fill();
393	+
394		eeTree->Fill();
395		ebTree->Fill();
396	<
396	>
397		hbheTree->Fill();
398		hfTree->Fill();
399	<
399	>	bkgTree->Fill();
400		}
401
402
#	Line 248 | Line 404 | RecHitTreeProducer::analyze(const edm::E
404		void
405		RecHitTreeProducer::beginJob()
406		{
407	<
407	>
408		hbheTree = fs->make<TTree>("hbhe","");
409		hbheTree->Branch("n",&hbheRecHit.n,"n/I");
410		hbheTree->Branch("e",hbheRecHit.e,"e[n]/F");
411		hbheTree->Branch("et",hbheRecHit.et,"et[n]/F");
412		hbheTree->Branch("eta",hbheRecHit.eta,"eta[n]/F");
413		hbheTree->Branch("phi",hbheRecHit.phi,"phi[n]/F");
414	<
414	>	hbheTree->Branch("isjet",hbheRecHit.isjet,"isjet[n]/O");
415	>
416		hfTree = fs->make<TTree>("hf","");
417		hfTree->Branch("n",&hfRecHit.n,"n/I");
418		hfTree->Branch("e",hfRecHit.e,"e[n]/F");
419		hfTree->Branch("et",hfRecHit.et,"et[n]/F");
420		hfTree->Branch("eta",hfRecHit.eta,"eta[n]/F");
421		hfTree->Branch("phi",hfRecHit.phi,"phi[n]/F");
422	+	hfTree->Branch("isjet",hfRecHit.isjet,"isjet[n]/O");
423
424		eeTree = fs->make<TTree>("ee","");
425		eeTree->Branch("n",&eeRecHit.n,"n/I");
#	Line 269 | Line 427 | RecHitTreeProducer::beginJob()
427		eeTree->Branch("et",eeRecHit.et,"et[n]/F");
428		eeTree->Branch("eta",eeRecHit.eta,"eta[n]/F");
429		eeTree->Branch("phi",eeRecHit.phi,"phi[n]/F");
430	<
430	>	eeTree->Branch("isjet",eeRecHit.isjet,"isjet[n]/O");
431	>
432		ebTree = fs->make<TTree>("eb","");
433		ebTree->Branch("n",&ebRecHit.n,"n/I");
434		ebTree->Branch("e",ebRecHit.e,"e[n]/F");
435		ebTree->Branch("et",ebRecHit.et,"et[n]/F");
436		ebTree->Branch("eta",ebRecHit.eta,"eta[n]/F");
437		ebTree->Branch("phi",ebRecHit.phi,"phi[n]/F");
438	+	ebTree->Branch("isjet",ebRecHit.isjet,"isjet[n]/O");
439	+
440	+	towerTree = fs->make<TTree>("tower","");
441	+	towerTree->Branch("n",&myTowers.n,"n/I");
442	+	towerTree->Branch("e",myTowers.e,"e[n]/F");
443	+	towerTree->Branch("et",myTowers.et,"et[n]/F");
444	+	towerTree->Branch("eta",myTowers.eta,"eta[n]/F");
445	+	towerTree->Branch("phi",myTowers.phi,"phi[n]/F");
446	+	towerTree->Branch("isjet",myTowers.isjet,"isjet[n]/O");
447	+
448	+
449	+	if(doBasicClusters_){
450	+	bcTree = fs->make<TTree>("bc","");
451	+	bcTree->Branch("n",&myBC.n,"n/I");
452	+	bcTree->Branch("e",myBC.e,"e[n]/F");
453	+	bcTree->Branch("et",myBC.et,"et[n]/F");
454	+	bcTree->Branch("eta",myBC.eta,"eta[n]/F");
455	+	bcTree->Branch("phi",myBC.phi,"phi[n]/F");
456	+	bcTree->Branch("jtpt",&myBC.jtpt,"jtpt/F");
457	+	bcTree->Branch("jteta",&myBC.jteta,"jteta/F");
458	+	bcTree->Branch("jtphi",&myBC.jtphi,"jtphi/F");
459	+	//     bcTree->Branch("isjet",bcRecHit.isjet,"isjet[n]/O");
460	+	}
461	+
462	+	if(doFastJet_){
463	+	bkgTree = fs->make<TTree>("bkg","");
464	+	bkgTree->Branch("n",&bkg.n,"n/I");
465	+	bkgTree->Branch("rho",bkg.rho,"rho[n]/F");
466	+	bkgTree->Branch("sigma",bkg.sigma,"sigma[n]/F");
467	+	}
468
469		}
470
#	Line 289 | Line 478 | double RecHitTreeProducer::getEt(const D
478		double et = energy*sin(pos.theta());
479		return et;
480		}
481	<
481	>
482		double RecHitTreeProducer::getEta(const DetId &id){
483		const GlobalPoint& pos=geo->getPosition(id);
484		double et = pos.eta();

Diff Legend

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