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.12 by yilmaz, Thu Nov 11 18:37:16 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	<
71	>	int depth[MAXHITS];
72		int n;
73
74		float e[MAXHITS];
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;
125	<
115	<	edm::Handle<reco::CaloJetCollection> signalJets;
124	>	edm::Handle<reco::BasicClusterCollection> bClusters;
125	>	edm::Handle<CaloTowerCollection> towers;
126
127	+	typedef vector<EcalRecHit>::const_iterator EcalIterator;
128	+
129	+	edm::Handle<reco::CaloJetCollection> jets;
130	+
131	+	edm::Handle<std::vector<double> > rhos;
132	+	edm::Handle<std::vector<double> > sigmas;
133	+
134		MyRecHit hbheRecHit;
135		MyRecHit hfRecHit;
136		MyRecHit ebRecHit;
137		MyRecHit eeRecHit;
138	+	MyRecHit myBC;
139	+	MyRecHit myTowers;
140	+	MyBkg bkg;
141
142	+	TNtuple* nt;
143		TTree* hbheTree;
144		TTree* hfTree;
145		TTree* ebTree;
146		TTree* eeTree;
147	+	TTree* bcTree;
148	+	TTree* towerTree;
149	+	TTree* bkgTree;
150	+
151		double cone;
152	+	double hfTowerThreshold_;
153	+	double hfLongThreshold_;
154	+	double hfShortThreshold_;
155	+	double hbheThreshold_;
156	+	double ebThreshold_;
157	+	double eeThreshold_;
158
159		edm::Service<TFileService> fs;
160		const CentralityBins * cbins_;
#	Line 133 | Line 164 | class RecHitTreeProducer : public edm::E
164		edm::InputTag HcalRecHitHBHESrc_;
165		edm::InputTag EBSrc_;
166		edm::InputTag EESrc_;
167	+	edm::InputTag BCSrc_;
168	+	edm::InputTag TowerSrc_;
169	+
170	+	edm::InputTag JetSrc_;
171	+
172	+	edm::InputTag FastJetTag_;
173
174	+	bool useJets_;
175	+	bool doBasicClusters_;
176	+	bool doTowers_;
177	+	bool doEcal_;
178	+	bool doHcal_;
179
180	+	bool doFastJet_;
181	+
182	+	bool doEbyEonly_;
183		};
184
185		//
#	Line 154 | Line 199 | RecHitTreeProducer::RecHitTreeProducer(c
199		cone(0.5)
200		{
201		//now do what ever initialization is needed
157	–
202		EBSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("EBRecHitSrc",edm::InputTag("ecalRecHit","EcalRecHitsEB"));
203		EESrc_ = iConfig.getUntrackedParameter<edm::InputTag>("EERecHitSrc",edm::InputTag("ecalRecHit","EcalRecHitsEE"));
204		HcalRecHitHFSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHFRecHitSrc",edm::InputTag("hfreco"));
205		HcalRecHitHBHESrc_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHBHERecHitSrc",edm::InputTag("hbhereco"));
206	+	BCSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("BasicClusterSrc1",edm::InputTag("ecalRecHit","EcalRecHitsEB","RECO"));
207	+	TowerSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("towersSrc",edm::InputTag("towerMaker"));
208	+	JetSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("JetSrc",edm::InputTag("iterativeConePu5CaloJets"));
209	+	useJets_ = iConfig.getUntrackedParameter<bool>("useJets",true);
210	+	doBasicClusters_ = iConfig.getUntrackedParameter<bool>("doBasicClusters",false);
211	+	doTowers_ = iConfig.getUntrackedParameter<bool>("doTowers",true);
212	+	doEcal_ = iConfig.getUntrackedParameter<bool>("doEcal",true);
213	+	doHcal_ = iConfig.getUntrackedParameter<bool>("doHcal",true);
214	+	doFastJet_ = iConfig.getUntrackedParameter<bool>("doFastJet",true);
215	+	FastJetTag_ = iConfig.getUntrackedParameter<edm::InputTag>("FastJetTag",edm::InputTag("kt4CaloJets"));
216	+	doEbyEonly_ = iConfig.getUntrackedParameter<bool>("doEbyEonly",false);
217	+	hfTowerThreshold_ = iConfig.getUntrackedParameter<double>("HFtowerMin",3.);
218	+	hfLongThreshold_ = iConfig.getUntrackedParameter<double>("HFlongMin",0.5);
219	+	hfShortThreshold_ = iConfig.getUntrackedParameter<double>("HFshortMin",0.85);
220		}
221
222
#	Line 184 | Line 242 | RecHitTreeProducer::analyze(const edm::E
242		hbheRecHit.n = 0;
243		ebRecHit.n = 0;
244		eeRecHit.n = 0;
245	+	myBC.n = 0;
246	+	myTowers.n = 0;
247	+	bkg.n = 0;
248
249	+	if(doEcal_){
250		ev.getByLabel(EBSrc_,ebHits);
251		ev.getByLabel(EESrc_,eeHits);
252	<
252	>	}
253	>	if(doHcal_){
254		ev.getByLabel(HcalRecHitHFSrc_,hfHits);
255		ev.getByLabel(HcalRecHitHBHESrc_,hbheHits);
256	<
257	<	if(0 && !cbins_) cbins_ = getCentralityBinsFromDB(iSetup);
256	>	}
257	>	if(useJets_) {
258	>	ev.getByLabel(JetSrc_,jets);
259	>	}
260	>
261	>	if(doBasicClusters_){
262	>	ev.getByLabel(BCSrc_,bClusters);
263	>	}
264	>
265	>	if(doTowers_){
266	>	ev.getByLabel(TowerSrc_,towers);
267	>	}
268	>
269	>	if(doFastJet_){
270	>	ev.getByLabel(edm::InputTag(FastJetTag_.label(),"rhos",FastJetTag_.process()),rhos);
271	>	ev.getByLabel(edm::InputTag(FastJetTag_.label(),"sigmas",FastJetTag_.process()),sigmas);
272	>	bkg.n = rhos->size();
273	>	for(unsigned int i = 0; i < rhos->size(); ++i){
274	>	bkg.rho[i] = (*rhos)[i];
275	>	bkg.sigma[i] = (*sigmas)[i];
276	>	}
277	>	}
278	>
279	>	if(0 && !cbins_) cbins_ = getCentralityBinsFromDB(iSetup);
280
281		if(!geo){
282		edm::ESHandle<CaloGeometry> pGeo;
#	Line 199 | Line 284 | RecHitTreeProducer::analyze(const edm::E
284		geo = pGeo.product();
285		}
286
287	+	int nHFlongPlus = 0;
288	+	int nHFshortPlus = 0;
289	+	int nHFtowerPlus = 0;
290	+	int nHFlongMinus = 0;
291	+	int nHFshortMinus = 0;
292	+	int nHFtowerMinus = 0;
293	+
294	+
295	+
296	+	if(doHcal_){
297		for(unsigned int i = 0; i < hfHits->size(); ++i){
298		const HFRecHit & hit= (*hfHits)[i];
299		hfRecHit.e[hfRecHit.n] = hit.energy();
300		hfRecHit.et[hfRecHit.n] = getEt(hit.id(),hit.energy());
301		hfRecHit.eta[hfRecHit.n] = getEta(hit.id());
302		hfRecHit.phi[hfRecHit.n] = getPhi(hit.id());
303	+	hfRecHit.isjet[hfRecHit.n] = false;
304	+	hfRecHit.depth[hfRecHit.n] = hit.id().depth();
305	+
306	+	if(hit.id().ieta() > 0){
307	+	if(hit.energy() > hfShortThreshold_ && hit.id().depth() != 1) nHFshortPlus++;
308	+	if(hit.energy() > hfLongThreshold_ && hit.id().depth() == 1) nHFlongPlus++;
309	+	}else{
310	+	if(hit.energy() > hfShortThreshold_ && hit.id().depth() != 1) nHFshortMinus++;
311	+	if(hit.energy() > hfLongThreshold_ && hit.id().depth() == 1) nHFlongMinus++;
312	+	}
313	+
314	+	if(useJets_){
315	+	for(unsigned int j = 0 ; j < jets->size(); ++j){
316	+	const reco::Jet& jet = (*jets)[j];
317	+	double dr = reco::deltaR(hfRecHit.eta[hfRecHit.n],hfRecHit.phi[hfRecHit.n],jet.eta(),jet.phi());
318	+	if(dr < cone){ hfRecHit.isjet[hfRecHit.n] = true; }
319	+	}
320	+	}
321		hfRecHit.n++;
322		}
323	<
323	>	if(!doEbyEonly_){
324		for(unsigned int i = 0; i < hbheHits->size(); ++i){
325		const HBHERecHit & hit= (*hbheHits)[i];
326	<	hbheRecHit.e[hfRecHit.n] = hit.energy();
327	<	hbheRecHit.et[hfRecHit.n] = getEt(hit.id(),hit.energy());
328	<	hbheRecHit.eta[hfRecHit.n] = getEta(hit.id());
329	<	hbheRecHit.phi[hfRecHit.n] = getPhi(hit.id());
326	>	hbheRecHit.e[hbheRecHit.n] = hit.energy();
327	>	hbheRecHit.et[hbheRecHit.n] = getEt(hit.id(),hit.energy());
328	>	hbheRecHit.eta[hbheRecHit.n] = getEta(hit.id());
329	>	hbheRecHit.phi[hbheRecHit.n] = getPhi(hit.id());
330	>	hbheRecHit.isjet[hbheRecHit.n] = false;
331	>	if(useJets_){
332	>	for(unsigned int j = 0 ; j < jets->size(); ++j){
333	>	const reco::Jet& jet = (*jets)[j];
334	>	double dr = reco::deltaR(hbheRecHit.eta[hbheRecHit.n],hbheRecHit.phi[hbheRecHit.n],jet.eta(),jet.phi());
335	>	if(dr < cone){ hbheRecHit.isjet[hbheRecHit.n] = true; }
336	>	}
337	>	}
338		hbheRecHit.n++;
339		}
340	<
340	>	}
341	>	}
342	>	if(doEcal_ && !doEbyEonly_){
343		for(unsigned int i = 0; i < ebHits->size(); ++i){
344		const EcalRecHit & hit= (*ebHits)[i];
345		ebRecHit.e[ebRecHit.n] = hit.energy();
346		ebRecHit.et[ebRecHit.n] = getEt(hit.id(),hit.energy());
347		ebRecHit.eta[ebRecHit.n] = getEta(hit.id());
348		ebRecHit.phi[ebRecHit.n] = getPhi(hit.id());
349	+	ebRecHit.isjet[ebRecHit.n] = false;
350	+	if(useJets_){
351	+	for(unsigned int j = 0 ; j < jets->size(); ++j){
352	+	const reco::Jet& jet = (*jets)[j];
353	+	double dr = reco::deltaR(ebRecHit.eta[ebRecHit.n],ebRecHit.phi[ebRecHit.n],jet.eta(),jet.phi());
354	+	if(dr < cone){ ebRecHit.isjet[ebRecHit.n] = true; }
355	+	}
356	+	}
357		ebRecHit.n++;
358		}
359	<
359	>
360		for(unsigned int i = 0; i < eeHits->size(); ++i){
361		const EcalRecHit & hit= (*eeHits)[i];
362		eeRecHit.e[eeRecHit.n] = hit.energy();
363		eeRecHit.et[eeRecHit.n] = getEt(hit.id(),hit.energy());
364		eeRecHit.eta[eeRecHit.n] = getEta(hit.id());
365		eeRecHit.phi[eeRecHit.n] = getPhi(hit.id());
366	+	eeRecHit.isjet[eeRecHit.n] = false;
367	+	if(useJets_){
368	+	for(unsigned int j = 0 ; j < jets->size(); ++j){
369	+	const reco::Jet& jet = (*jets)[j];
370	+	double dr = reco::deltaR(eeRecHit.eta[eeRecHit.n],eeRecHit.phi[eeRecHit.n],jet.eta(),jet.phi());
371	+	if(dr < cone){ eeRecHit.isjet[eeRecHit.n] = true; }
372	+	}
373	+	}
374		eeRecHit.n++;
375		}
376	+	}
377	+
378	+	if(doTowers_){
379	+
380	+	for(unsigned int i = 0; i < towers->size(); ++i){
381	+	const CaloTower & hit= (*towers)[i];
382	+	myTowers.e[myTowers.n] = hit.energy();
383	+	myTowers.et[myTowers.n] = getEt(hit.id(),hit.energy());
384	+	myTowers.eta[myTowers.n] = getEta(hit.id());
385	+	myTowers.phi[myTowers.n] = getPhi(hit.id());
386	+	myTowers.isjet[myTowers.n] = false;
387	+
388	+	if(hit.ieta() > 29 && hit.energy() > hfTowerThreshold_) nHFtowerPlus++;
389	+	if(hit.ieta() < -29 && hit.energy() > hfTowerThreshold_) nHFtowerMinus++;
390	+
391	+	if(useJets_){
392	+	for(unsigned int j = 0 ; j < jets->size(); ++j){
393	+	const reco::Jet& jet = (*jets)[j];
394	+	double dr = reco::deltaR(myTowers.eta[myTowers.n],myTowers.phi[myTowers.n],jet.eta(),jet.phi());
395	+	if(dr < cone){ myTowers.isjet[myTowers.n] = true; }
396	+	}
397	+	}
398	+	myTowers.n++;
399	+	}
400	+
401	+	}
402
403	<	eeTree->Fill();
404	<	ebTree->Fill();
403	>	if(doBasicClusters_ && !doEbyEonly_){
404	>	for(unsigned int j = 0 ; j < jets->size(); ++j){
405	>	const reco::Jet& jet = (*jets)[j];
406	>	myBC.n = 0;
407	>	myBC.jtpt = jet.pt();
408	>	myBC.jteta = jet.eta();
409	>	myBC.jtphi = jet.phi();
410	>
411	>	for(unsigned int i = 0; i < bClusters->size(); ++i){
412	>	const reco::BasicCluster & bc= (*bClusters)[i];
413	>	double dr = reco::deltaR(bc.eta(),bc.phi(),jet.eta(),jet.phi());
414	>	if(dr < cone){
415	>	myBC.e[myBC.n] = bc.energy();
416	>	myBC.et[myBC.n] = bc.energy()*sin(bc.position().theta());
417	>	myBC.eta[myBC.n] = bc.eta();
418	>	myBC.phi[myBC.n] = bc.phi();
419	>	myBC.n++;
420	>	}
421	>	}
422	>	bcTree->Fill();
423	>	}
424	>	}
425
426	<	hbheTree->Fill();
427	<	hfTree->Fill();
426	>	if(!doEbyEonly_){
427	>	towerTree->Fill();
428	>
429	>	eeTree->Fill();
430	>	ebTree->Fill();
431	>
432	>	hbheTree->Fill();
433	>	hfTree->Fill();
434	>
435	>	if (doFastJet_) {
436	>	bkgTree->Fill();
437	>	}
438	>	}
439
440	+	nt->Fill(nHFtowerPlus,nHFtowerMinus,nHFlongPlus,nHFlongMinus,nHFshortPlus,nHFshortMinus);
441	+
442		}
443
444
#	Line 248 | Line 446 | RecHitTreeProducer::analyze(const edm::E
446		void
447		RecHitTreeProducer::beginJob()
448		{
449	<
449	>
450		hbheTree = fs->make<TTree>("hbhe","");
451		hbheTree->Branch("n",&hbheRecHit.n,"n/I");
452		hbheTree->Branch("e",hbheRecHit.e,"e[n]/F");
453		hbheTree->Branch("et",hbheRecHit.et,"et[n]/F");
454		hbheTree->Branch("eta",hbheRecHit.eta,"eta[n]/F");
455		hbheTree->Branch("phi",hbheRecHit.phi,"phi[n]/F");
456	<
456	>	hbheTree->Branch("isjet",hbheRecHit.isjet,"isjet[n]/O");
457	>
458		hfTree = fs->make<TTree>("hf","");
459		hfTree->Branch("n",&hfRecHit.n,"n/I");
460		hfTree->Branch("e",hfRecHit.e,"e[n]/F");
461		hfTree->Branch("et",hfRecHit.et,"et[n]/F");
462		hfTree->Branch("eta",hfRecHit.eta,"eta[n]/F");
463		hfTree->Branch("phi",hfRecHit.phi,"phi[n]/F");
464	+	hfTree->Branch("depth",hfRecHit.depth,"depth[n]/I");
465	+	hfTree->Branch("isjet",hfRecHit.isjet,"isjet[n]/O");
466
467		eeTree = fs->make<TTree>("ee","");
468		eeTree->Branch("n",&eeRecHit.n,"n/I");
#	Line 269 | Line 470 | RecHitTreeProducer::beginJob()
470		eeTree->Branch("et",eeRecHit.et,"et[n]/F");
471		eeTree->Branch("eta",eeRecHit.eta,"eta[n]/F");
472		eeTree->Branch("phi",eeRecHit.phi,"phi[n]/F");
473	<
473	>	eeTree->Branch("isjet",eeRecHit.isjet,"isjet[n]/O");
474	>
475		ebTree = fs->make<TTree>("eb","");
476		ebTree->Branch("n",&ebRecHit.n,"n/I");
477		ebTree->Branch("e",ebRecHit.e,"e[n]/F");
478		ebTree->Branch("et",ebRecHit.et,"et[n]/F");
479		ebTree->Branch("eta",ebRecHit.eta,"eta[n]/F");
480		ebTree->Branch("phi",ebRecHit.phi,"phi[n]/F");
481	+	ebTree->Branch("isjet",ebRecHit.isjet,"isjet[n]/O");
482	+
483	+	towerTree = fs->make<TTree>("tower","");
484	+	towerTree->Branch("n",&myTowers.n,"n/I");
485	+	towerTree->Branch("e",myTowers.e,"e[n]/F");
486	+	towerTree->Branch("et",myTowers.et,"et[n]/F");
487	+	towerTree->Branch("eta",myTowers.eta,"eta[n]/F");
488	+	towerTree->Branch("phi",myTowers.phi,"phi[n]/F");
489	+	towerTree->Branch("isjet",myTowers.isjet,"isjet[n]/O");
490	+
491	+
492	+	if(doBasicClusters_){
493	+	bcTree = fs->make<TTree>("bc","");
494	+	bcTree->Branch("n",&myBC.n,"n/I");
495	+	bcTree->Branch("e",myBC.e,"e[n]/F");
496	+	bcTree->Branch("et",myBC.et,"et[n]/F");
497	+	bcTree->Branch("eta",myBC.eta,"eta[n]/F");
498	+	bcTree->Branch("phi",myBC.phi,"phi[n]/F");
499	+	bcTree->Branch("jtpt",&myBC.jtpt,"jtpt/F");
500	+	bcTree->Branch("jteta",&myBC.jteta,"jteta/F");
501	+	bcTree->Branch("jtphi",&myBC.jtphi,"jtphi/F");
502	+	//     bcTree->Branch("isjet",bcRecHit.isjet,"isjet[n]/O");
503	+	}
504	+
505	+	if(doFastJet_){
506	+	bkgTree = fs->make<TTree>("bkg","");
507	+	bkgTree->Branch("n",&bkg.n,"n/I");
508	+	bkgTree->Branch("rho",bkg.rho,"rho[n]/F");
509	+	bkgTree->Branch("sigma",bkg.sigma,"sigma[n]/F");
510	+	}
511	+
512	+	nt = fs->make<TNtuple>("ntEvent","","nHFplus:nHFminus:nHFlongPlus:nHFlongMinus:nHFshortPlus:nHFshortMinus");
513
514		}
515
#	Line 289 | Line 523 | double RecHitTreeProducer::getEt(const D
523		double et = energy*sin(pos.theta());
524		return et;
525		}
526	<
526	>
527		double RecHitTreeProducer::getEta(const DetId &id){
528		const GlobalPoint& pos=geo->getPosition(id);
529		double et = pos.eta();

Diff Legend

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