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.6 by yilmaz, Mon Nov 1 21:48:31 2010 UTC vs.
Revision 1.13 by frankma, Thu Sep 22 08:51:47 2011 UTC

#	Line 68 | Line 68 | using namespace std;
68		#define MAXHITS 1000000
69
70		struct MyRecHit{
71	<
71	>	int depth[MAXHITS];
72		int n;
73
74		float e[MAXHITS];
#	Line 84 | Line 84 | struct MyRecHit{
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 118 | Line 124 | class RecHitTreeProducer : public edm::E
124		edm::Handle<reco::BasicClusterCollection> bClusters;
125		edm::Handle<CaloTowerCollection> towers;
126
121	–
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;
#	Line 129 | Line 137 | class RecHitTreeProducer : public edm::E
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	<
152	>	double hfTowerThreshold_;
153	>	double hfLongThreshold_;
154	>	double hfShortThreshold_;
155	>	double hbheThreshold_;
156	>	double ebThreshold_;
157	>	double eeThreshold_;
158	>
159	>	double hbhePtMin_;
160	>	double hfPtMin_;
161	>	double ebPtMin_;
162	>	double eePtMin_;
163	>	double towerPtMin_;
164	>
165		edm::Service<TFileService> fs;
166		const CentralityBins * cbins_;
167		const CaloGeometry *geo;
#	Line 148 | Line 171 | class RecHitTreeProducer : public edm::E
171		edm::InputTag EBSrc_;
172		edm::InputTag EESrc_;
173		edm::InputTag BCSrc_;
151	–
174		edm::InputTag TowerSrc_;
175
176		edm::InputTag JetSrc_;
177
178	+	edm::InputTag FastJetTag_;
179	+
180		bool useJets_;
181		bool doBasicClusters_;
182	+	bool doTowers_;
183	+	bool doEcal_;
184	+	bool doHcal_;
185	+
186	+	bool doFastJet_;
187
188	+	bool doEbyEonly_;
189		};
190
191		//
#	Line 181 | Line 211 | RecHitTreeProducer::RecHitTreeProducer(c
211		HcalRecHitHBHESrc_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHBHERecHitSrc",edm::InputTag("hbhereco"));
212		BCSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("BasicClusterSrc1",edm::InputTag("ecalRecHit","EcalRecHitsEB","RECO"));
213		TowerSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("towersSrc",edm::InputTag("towerMaker"));
214	<	JetSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("JetSrc",edm::InputTag("iterativeCone5CaloJets"));
215	<	useJets_ = iConfig.getUntrackedParameter<bool>("useJets",false);
214	>	JetSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("JetSrc",edm::InputTag("iterativeConePu5CaloJets"));
215	>	useJets_ = iConfig.getUntrackedParameter<bool>("useJets",true);
216		doBasicClusters_ = iConfig.getUntrackedParameter<bool>("doBasicClusters",false);
217	<
217	>	doTowers_ = iConfig.getUntrackedParameter<bool>("doTowers",true);
218	>	doEcal_ = iConfig.getUntrackedParameter<bool>("doEcal",true);
219	>	doHcal_ = iConfig.getUntrackedParameter<bool>("doHcal",true);
220	>	doFastJet_ = iConfig.getUntrackedParameter<bool>("doFastJet",true);
221	>	FastJetTag_ = iConfig.getUntrackedParameter<edm::InputTag>("FastJetTag",edm::InputTag("kt4CaloJets"));
222	>	doEbyEonly_ = iConfig.getUntrackedParameter<bool>("doEbyEonly",false);
223	>	hfTowerThreshold_ = iConfig.getUntrackedParameter<double>("HFtowerMin",3.);
224	>	hfLongThreshold_ = iConfig.getUntrackedParameter<double>("HFlongMin",0.5);
225	>	hfShortThreshold_ = iConfig.getUntrackedParameter<double>("HFshortMin",0.85);
226	>	hbhePtMin_ = iConfig.getUntrackedParameter<double>("HBHETreePtMin",0);
227	>	hfPtMin_ = iConfig.getUntrackedParameter<double>("HFTreePtMin",0);
228	>	ebPtMin_ = iConfig.getUntrackedParameter<double>("EBTreePtMin",0);
229	>	eePtMin_ = iConfig.getUntrackedParameter<double>("EETreePtMin",0.);
230	>	towerPtMin_ = iConfig.getUntrackedParameter<double>("TowerTreePtMin",0.);
231		}
232
233
#	Line 210 | Line 253 | RecHitTreeProducer::analyze(const edm::E
253		hbheRecHit.n = 0;
254		ebRecHit.n = 0;
255		eeRecHit.n = 0;
256	+	myBC.n = 0;
257	+	myTowers.n = 0;
258	+	bkg.n = 0;
259
260	+	if(doEcal_){
261		ev.getByLabel(EBSrc_,ebHits);
262		ev.getByLabel(EESrc_,eeHits);
263	<
263	>	}
264	>	if(doHcal_){
265		ev.getByLabel(HcalRecHitHFSrc_,hfHits);
266		ev.getByLabel(HcalRecHitHBHESrc_,hbheHits);
267	<
267	>	}
268		if(useJets_) {
269		ev.getByLabel(JetSrc_,jets);
270		}
#	Line 225 | Line 273 | RecHitTreeProducer::analyze(const edm::E
273		ev.getByLabel(BCSrc_,bClusters);
274		}
275
276	+	if(doTowers_){
277	+	ev.getByLabel(TowerSrc_,towers);
278	+	}
279
280	+	if(doFastJet_){
281	+	ev.getByLabel(edm::InputTag(FastJetTag_.label(),"rhos",FastJetTag_.process()),rhos);
282	+	ev.getByLabel(edm::InputTag(FastJetTag_.label(),"sigmas",FastJetTag_.process()),sigmas);
283	+	bkg.n = rhos->size();
284	+	for(unsigned int i = 0; i < rhos->size(); ++i){
285	+	bkg.rho[i] = (*rhos)[i];
286	+	bkg.sigma[i] = (*sigmas)[i];
287	+	}
288	+	}
289
290		if(0 && !cbins_) cbins_ = getCentralityBinsFromDB(iSetup);
291
#	Line 234 | Line 294 | RecHitTreeProducer::analyze(const edm::E
294		iSetup.get<CaloGeometryRecord>().get(pGeo);
295		geo = pGeo.product();
296		}
297	+
298	+	int nHFlongPlus = 0;
299	+	int nHFshortPlus = 0;
300	+	int nHFtowerPlus = 0;
301	+	int nHFlongMinus = 0;
302	+	int nHFshortMinus = 0;
303	+	int nHFtowerMinus = 0;
304	+
305	+
306
307	+	if(doHcal_){
308		for(unsigned int i = 0; i < hfHits->size(); ++i){
309		const HFRecHit & hit= (*hfHits)[i];
310		hfRecHit.e[hfRecHit.n] = hit.energy();
#	Line 242 | Line 312 | RecHitTreeProducer::analyze(const edm::E
312		hfRecHit.eta[hfRecHit.n] = getEta(hit.id());
313		hfRecHit.phi[hfRecHit.n] = getPhi(hit.id());
314		hfRecHit.isjet[hfRecHit.n] = false;
315	<	if(!useJets_){
315	>	hfRecHit.depth[hfRecHit.n] = hit.id().depth();
316	>
317	>	if(hit.id().ieta() > 0){
318	>	if(hit.energy() > hfShortThreshold_ && hit.id().depth() != 1) nHFshortPlus++;
319	>	if(hit.energy() > hfLongThreshold_ && hit.id().depth() == 1) nHFlongPlus++;
320	>	}else{
321	>	if(hit.energy() > hfShortThreshold_ && hit.id().depth() != 1) nHFshortMinus++;
322	>	if(hit.energy() > hfLongThreshold_ && hit.id().depth() == 1) nHFlongMinus++;
323	>	}
324	>
325	>	if(useJets_){
326		for(unsigned int j = 0 ; j < jets->size(); ++j){
327		const reco::Jet& jet = (*jets)[j];
328		double dr = reco::deltaR(hfRecHit.eta[hfRecHit.n],hfRecHit.phi[hfRecHit.n],jet.eta(),jet.phi());
329		if(dr < cone){ hfRecHit.isjet[hfRecHit.n] = true; }
330		}
331		}
332	<	hfRecHit.n++;
332	>	if (hfRecHit.et[hfRecHit.n]>=hfPtMin_) hfRecHit.n++;
333		}
334	<
334	>	if(!doEbyEonly_){
335		for(unsigned int i = 0; i < hbheHits->size(); ++i){
336		const HBHERecHit & hit= (*hbheHits)[i];
337	+	if (getEt(hit.id(),hit.energy())<hbhePtMin_) continue;
338		hbheRecHit.e[hbheRecHit.n] = hit.energy();
339		hbheRecHit.et[hbheRecHit.n] = getEt(hit.id(),hit.energy());
340		hbheRecHit.eta[hbheRecHit.n] = getEta(hit.id());
341		hbheRecHit.phi[hbheRecHit.n] = getPhi(hit.id());
342		hbheRecHit.isjet[hbheRecHit.n] = false;
343	<	if(!useJets_){
343	>	if(useJets_){
344		for(unsigned int j = 0 ; j < jets->size(); ++j){
345		const reco::Jet& jet = (*jets)[j];
346		double dr = reco::deltaR(hbheRecHit.eta[hbheRecHit.n],hbheRecHit.phi[hbheRecHit.n],jet.eta(),jet.phi());
#	Line 268 | Line 349 | RecHitTreeProducer::analyze(const edm::E
349		}
350		hbheRecHit.n++;
351		}
352	<
352	>	}
353	>	}
354	>	if(doEcal_ && !doEbyEonly_){
355		for(unsigned int i = 0; i < ebHits->size(); ++i){
356		const EcalRecHit & hit= (*ebHits)[i];
357	+	if (getEt(hit.id(),hit.energy())<ebPtMin_) continue;
358		ebRecHit.e[ebRecHit.n] = hit.energy();
359		ebRecHit.et[ebRecHit.n] = getEt(hit.id(),hit.energy());
360		ebRecHit.eta[ebRecHit.n] = getEta(hit.id());
361		ebRecHit.phi[ebRecHit.n] = getPhi(hit.id());
362		ebRecHit.isjet[ebRecHit.n] = false;
363	<	if(!useJets_){
363	>	if(useJets_){
364		for(unsigned int j = 0 ; j < jets->size(); ++j){
365		const reco::Jet& jet = (*jets)[j];
366		double dr = reco::deltaR(ebRecHit.eta[ebRecHit.n],ebRecHit.phi[ebRecHit.n],jet.eta(),jet.phi());
#	Line 288 | Line 372 | RecHitTreeProducer::analyze(const edm::E
372
373		for(unsigned int i = 0; i < eeHits->size(); ++i){
374		const EcalRecHit & hit= (*eeHits)[i];
375	+	if (getEt(hit.id(),hit.energy())<eePtMin_) continue;
376		eeRecHit.e[eeRecHit.n] = hit.energy();
377		eeRecHit.et[eeRecHit.n] = getEt(hit.id(),hit.energy());
378		eeRecHit.eta[eeRecHit.n] = getEta(hit.id());
379		eeRecHit.phi[eeRecHit.n] = getPhi(hit.id());
380		eeRecHit.isjet[eeRecHit.n] = false;
381	<	if(!useJets_){
381	>	if(useJets_){
382		for(unsigned int j = 0 ; j < jets->size(); ++j){
383		const reco::Jet& jet = (*jets)[j];
384		double dr = reco::deltaR(eeRecHit.eta[eeRecHit.n],eeRecHit.phi[eeRecHit.n],jet.eta(),jet.phi());
#	Line 302 | Line 387 | RecHitTreeProducer::analyze(const edm::E
387		}
388		eeRecHit.n++;
389		}
390	+	}
391
392	+	if(doTowers_){
393	+
394	+	for(unsigned int i = 0; i < towers->size(); ++i){
395	+	const CaloTower & hit= (*towers)[i];
396	+	if (getEt(hit.id(),hit.energy())<towerPtMin_) continue;
397	+	myTowers.e[myTowers.n] = hit.energy();
398	+	myTowers.et[myTowers.n] = getEt(hit.id(),hit.energy());
399	+	myTowers.eta[myTowers.n] = getEta(hit.id());
400	+	myTowers.phi[myTowers.n] = getPhi(hit.id());
401	+	myTowers.isjet[myTowers.n] = false;
402	+
403	+	if(hit.ieta() > 29 && hit.energy() > hfTowerThreshold_) nHFtowerPlus++;
404	+	if(hit.ieta() < -29 && hit.energy() > hfTowerThreshold_) nHFtowerMinus++;
405	+
406	+	if(useJets_){
407	+	for(unsigned int j = 0 ; j < jets->size(); ++j){
408	+	const reco::Jet& jet = (*jets)[j];
409	+	double dr = reco::deltaR(myTowers.eta[myTowers.n],myTowers.phi[myTowers.n],jet.eta(),jet.phi());
410	+	if(dr < cone){ myTowers.isjet[myTowers.n] = true; }
411	+	}
412	+	}
413	+	myTowers.n++;
414	+	}
415	+
416	+	}
417
418	<	if(doBasicClusters_){
418	>	if(doBasicClusters_ && !doEbyEonly_){
419		for(unsigned int j = 0 ; j < jets->size(); ++j){
420		const reco::Jet& jet = (*jets)[j];
421		myBC.n = 0;
#	Line 326 | Line 437 | RecHitTreeProducer::analyze(const edm::E
437		bcTree->Fill();
438		}
439		}
440	<
441	<	eeTree->Fill();
442	<	ebTree->Fill();
443	<
444	<	hbheTree->Fill();
445	<	hfTree->Fill();
440	>
441	>	if(!doEbyEonly_){
442	>	towerTree->Fill();
443	>
444	>	eeTree->Fill();
445	>	ebTree->Fill();
446	>
447	>	hbheTree->Fill();
448	>	hfTree->Fill();
449	>
450	>	if (doFastJet_) {
451	>	bkgTree->Fill();
452	>	}
453	>	}
454	>
455	>	nt->Fill(nHFtowerPlus,nHFtowerMinus,nHFlongPlus,nHFlongMinus,nHFshortPlus,nHFshortMinus);
456
457		}
458
#	Line 355 | Line 476 | RecHitTreeProducer::beginJob()
476		hfTree->Branch("et",hfRecHit.et,"et[n]/F");
477		hfTree->Branch("eta",hfRecHit.eta,"eta[n]/F");
478		hfTree->Branch("phi",hfRecHit.phi,"phi[n]/F");
479	+	hfTree->Branch("depth",hfRecHit.depth,"depth[n]/I");
480		hfTree->Branch("isjet",hfRecHit.isjet,"isjet[n]/O");
481
482		eeTree = fs->make<TTree>("ee","");
#	Line 373 | Line 495 | RecHitTreeProducer::beginJob()
495		ebTree->Branch("phi",ebRecHit.phi,"phi[n]/F");
496		ebTree->Branch("isjet",ebRecHit.isjet,"isjet[n]/O");
497
498	+	towerTree = fs->make<TTree>("tower","");
499	+	towerTree->Branch("n",&myTowers.n,"n/I");
500	+	towerTree->Branch("e",myTowers.e,"e[n]/F");
501	+	towerTree->Branch("et",myTowers.et,"et[n]/F");
502	+	towerTree->Branch("eta",myTowers.eta,"eta[n]/F");
503	+	towerTree->Branch("phi",myTowers.phi,"phi[n]/F");
504	+	towerTree->Branch("isjet",myTowers.isjet,"isjet[n]/O");
505	+
506	+
507		if(doBasicClusters_){
508		bcTree = fs->make<TTree>("bc","");
509		bcTree->Branch("n",&myBC.n,"n/I");
#	Line 386 | Line 517 | RecHitTreeProducer::beginJob()
517		//     bcTree->Branch("isjet",bcRecHit.isjet,"isjet[n]/O");
518		}
519
520	<
520	>	if(doFastJet_){
521	>	bkgTree = fs->make<TTree>("bkg","");
522	>	bkgTree->Branch("n",&bkg.n,"n/I");
523	>	bkgTree->Branch("rho",bkg.rho,"rho[n]/F");
524	>	bkgTree->Branch("sigma",bkg.sigma,"sigma[n]/F");
525	>	}
526	>
527	>	nt = fs->make<TNtuple>("ntEvent","","nHFplus:nHFminus:nHFlongPlus:nHFlongMinus:nHFshortPlus:nHFshortMinus");
528
529		}
530

Diff Legend

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