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.9 by yilmaz, Tue Nov 9 22:54:45 2010 UTC vs.
Revision 1.16 by yjlee, Tue Oct 25 13:41:08 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		//
19		//
20
21	<
21	>	#define versionTag "v1"
22		// system include files
23		#include <memory>
24		#include <vector>
#	Line 60 | Line 61
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
#	Line 76 | Line 80 | struct MyRecHit{
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;
#	Line 102 | 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 123 | Line 133 | class RecHitTreeProducer : public edm::E
133
134		edm::Handle<reco::BasicClusterCollection> bClusters;
135		edm::Handle<CaloTowerCollection> towers;
136	<
136	>	edm::Handle<reco::VertexCollection> vtxs;
137
138		typedef vector<EcalRecHit>::const_iterator EcalIterator;
139
#	Line 156 | Line 166 | class RecHitTreeProducer : public edm::E
166		double hbheThreshold_;
167		double ebThreshold_;
168		double eeThreshold_;
169	<
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;
#	Line 167 | Line 183 | class RecHitTreeProducer : public edm::E
183		edm::InputTag EESrc_;
184		edm::InputTag BCSrc_;
185		edm::InputTag TowerSrc_;
186	+	edm::InputTag VtxSrc_;
187
188		edm::InputTag JetSrc_;
189
#	Line 177 | Line 194 | class RecHitTreeProducer : public edm::E
194		bool doTowers_;
195		bool doEcal_;
196		bool doHcal_;
197	+	bool hasVtx_;
198
199		bool doFastJet_;
200
#	Line 206 | Line 224 | RecHitTreeProducer::RecHitTreeProducer(c
224		HcalRecHitHBHESrc_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHBHERecHitSrc",edm::InputTag("hbhereco"));
225		BCSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("BasicClusterSrc1",edm::InputTag("ecalRecHit","EcalRecHitsEB","RECO"));
226		TowerSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("towersSrc",edm::InputTag("towerMaker"));
227	<	JetSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("JetSrc",edm::InputTag("iterativeCone5CaloJets"));
227	>	VtxSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("vtxSrc",edm::InputTag("hiSelectedVertex"));
228	>	JetSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("JetSrc",edm::InputTag("iterativeConePu5CaloJets"));
229		useJets_ = iConfig.getUntrackedParameter<bool>("useJets",true);
230		doBasicClusters_ = iConfig.getUntrackedParameter<bool>("doBasicClusters",false);
231		doTowers_ = iConfig.getUntrackedParameter<bool>("doTowers",true);
232		doEcal_ = iConfig.getUntrackedParameter<bool>("doEcal",true);
233		doHcal_ = iConfig.getUntrackedParameter<bool>("doHcal",true);
234	+	hasVtx_ = iConfig.getUntrackedParameter<bool>("hasVtx",false);
235		doFastJet_ = iConfig.getUntrackedParameter<bool>("doFastJet",true);
236		FastJetTag_ = iConfig.getUntrackedParameter<edm::InputTag>("FastJetTag",edm::InputTag("kt4CaloJets"));
237		doEbyEonly_ = iConfig.getUntrackedParameter<bool>("doEbyEonly",false);
238		hfTowerThreshold_ = iConfig.getUntrackedParameter<double>("HFtowerMin",3.);
239	<	hfLongThreshold_ = iConfig.getUntrackedParameter<double>("HFlongMin",3.);
240	<	hfShortThreshold_ = iConfig.getUntrackedParameter<double>("HFshortMin",3.);
239	>	hfLongThreshold_ = iConfig.getUntrackedParameter<double>("HFlongMin",0.5);
240	>	hfShortThreshold_ = iConfig.getUntrackedParameter<double>("HFshortMin",0.85);
241	>	hbhePtMin_ = iConfig.getUntrackedParameter<double>("HBHETreePtMin",0);
242	>	hfPtMin_ = iConfig.getUntrackedParameter<double>("HFTreePtMin",0);
243	>	ebPtMin_ = iConfig.getUntrackedParameter<double>("EBTreePtMin",0);
244	>	eePtMin_ = iConfig.getUntrackedParameter<double>("EETreePtMin",0.);
245	>	towerPtMin_ = iConfig.getUntrackedParameter<double>("TowerTreePtMin",0.);
246		}
247
248
#	Line 247 | Line 272 | RecHitTreeProducer::analyze(const edm::E
272		myTowers.n = 0;
273		bkg.n = 0;
274
275	+	// get vertex
276	+	reco::Vertex::Point vtx;
277	+	if (hasVtx_) vtx = getVtx(ev);
278	+
279		if(doEcal_){
280		ev.getByLabel(EBSrc_,ebHits);
281		ev.getByLabel(EESrc_,eeHits);
#	Line 285 | Line 314 | RecHitTreeProducer::analyze(const edm::E
314		geo = pGeo.product();
315		}
316
317	<	int nHFlong = 0;
318	<	int nHFshort = 0;
319	<	int nHFtower = 0;
317	>	int nHFlongPlus = 0;
318	>	int nHFshortPlus = 0;
319	>	int nHFtowerPlus = 0;
320	>	int nHFlongMinus = 0;
321	>	int nHFshortMinus = 0;
322	>	int nHFtowerMinus = 0;
323	>
324	>
325
326		if(doHcal_){
327		for(unsigned int i = 0; i < hfHits->size(); ++i){
#	Line 298 | Line 332 | RecHitTreeProducer::analyze(const edm::E
332		hfRecHit.phi[hfRecHit.n] = getPhi(hit.id());
333		hfRecHit.isjet[hfRecHit.n] = false;
334		hfRecHit.depth[hfRecHit.n] = hit.id().depth();
335	<	if(hit.energy() > hfShortThreshold_ && hit.id().depth() != 1) nHFshort++;
336	<	if(hit.energy() > hfLongThreshold_ && hit.id().depth() == 1) nHFlong++;
335	>
336	>	if(hit.id().ieta() > 0){
337	>	if(hit.energy() > hfShortThreshold_ && hit.id().depth() != 1) nHFshortPlus++;
338	>	if(hit.energy() > hfLongThreshold_ && hit.id().depth() == 1) nHFlongPlus++;
339	>	}else{
340	>	if(hit.energy() > hfShortThreshold_ && hit.id().depth() != 1) nHFshortMinus++;
341	>	if(hit.energy() > hfLongThreshold_ && hit.id().depth() == 1) nHFlongMinus++;
342	>	}
343
344		if(useJets_){
345		for(unsigned int j = 0 ; j < jets->size(); ++j){
#	Line 308 | Line 348 | RecHitTreeProducer::analyze(const edm::E
348		if(dr < cone){ hfRecHit.isjet[hfRecHit.n] = true; }
349		}
350		}
351	<	hfRecHit.n++;
351	>	if (hfRecHit.et[hfRecHit.n]>=hfPtMin_) hfRecHit.n++;
352		}
353		if(!doEbyEonly_){
354		for(unsigned int i = 0; i < hbheHits->size(); ++i){
355		const HBHERecHit & hit= (*hbheHits)[i];
356	+	if (getEt(hit.id(),hit.energy())<hbhePtMin_) continue;
357		hbheRecHit.e[hbheRecHit.n] = hit.energy();
358		hbheRecHit.et[hbheRecHit.n] = getEt(hit.id(),hit.energy());
359		hbheRecHit.eta[hbheRecHit.n] = getEta(hit.id());
#	Line 332 | Line 373 | RecHitTreeProducer::analyze(const edm::E
373		if(doEcal_ && !doEbyEonly_){
374		for(unsigned int i = 0; i < ebHits->size(); ++i){
375		const EcalRecHit & hit= (*ebHits)[i];
376	+	if (getEt(hit.id(),hit.energy())<ebPtMin_) continue;
377		ebRecHit.e[ebRecHit.n] = hit.energy();
378		ebRecHit.et[ebRecHit.n] = getEt(hit.id(),hit.energy());
379		ebRecHit.eta[ebRecHit.n] = getEta(hit.id());
#	Line 349 | Line 391 | RecHitTreeProducer::analyze(const edm::E
391
392		for(unsigned int i = 0; i < eeHits->size(); ++i){
393		const EcalRecHit & hit= (*eeHits)[i];
394	+	if (getEt(hit.id(),hit.energy())<eePtMin_) continue;
395		eeRecHit.e[eeRecHit.n] = hit.energy();
396		eeRecHit.et[eeRecHit.n] = getEt(hit.id(),hit.energy());
397		eeRecHit.eta[eeRecHit.n] = getEta(hit.id());
#	Line 369 | Line 412 | RecHitTreeProducer::analyze(const edm::E
412
413		for(unsigned int i = 0; i < towers->size(); ++i){
414		const CaloTower & hit= (*towers)[i];
415	+	if (getEt(hit.id(),hit.energy())<towerPtMin_) continue;
416		myTowers.e[myTowers.n] = hit.energy();
417		myTowers.et[myTowers.n] = getEt(hit.id(),hit.energy());
418		myTowers.eta[myTowers.n] = getEta(hit.id());
419		myTowers.phi[myTowers.n] = getPhi(hit.id());
420		myTowers.isjet[myTowers.n] = false;
421	+	if (hasVtx_) {
422	+	myTowers.etvtx[myTowers.n] = hit.p4(vtx).Et();
423	+	myTowers.etavtx[myTowers.n] = hit.p4(vtx).Eta();
424	+	myTowers.emEtVtx[myTowers.n] = hit.emEt(vtx);
425	+	myTowers.hadEtVtx[myTowers.n] = hit.hadEt(vtx);
426	+	}
427
428	<	if(hit.ieta() > 29 && hit.energy() > hfTowerThreshold_) nHFtower++;
428	>	if(hit.ieta() > 29 && hit.energy() > hfTowerThreshold_) nHFtowerPlus++;
429	>	if(hit.ieta() < -29 && hit.energy() > hfTowerThreshold_) nHFtowerMinus++;
430
431		if(useJets_){
432		for(unsigned int j = 0 ; j < jets->size(); ++j){
#	Line 426 | Line 477 | RecHitTreeProducer::analyze(const edm::E
477		}
478		}
479
480	<	nt->Fill(nHFtower,nHFlong,nHFshort);
480	>	nt->Fill(nHFtowerPlus,nHFtowerMinus,nHFlongPlus,nHFlongMinus,nHFshortPlus,nHFshortMinus);
481
482		}
483
#	Line 436 | Line 487 | void
487		RecHitTreeProducer::beginJob()
488		{
489
490	<	hbheTree = fs->make<TTree>("hbhe","");
490	>	hbheTree = fs->make<TTree>("hbhe",versionTag);
491		hbheTree->Branch("n",&hbheRecHit.n,"n/I");
492		hbheTree->Branch("e",hbheRecHit.e,"e[n]/F");
493		hbheTree->Branch("et",hbheRecHit.et,"et[n]/F");
#	Line 444 | Line 495 | RecHitTreeProducer::beginJob()
495		hbheTree->Branch("phi",hbheRecHit.phi,"phi[n]/F");
496		hbheTree->Branch("isjet",hbheRecHit.isjet,"isjet[n]/O");
497
498	<	hfTree = fs->make<TTree>("hf","");
498	>	hfTree = fs->make<TTree>("hf",versionTag);
499		hfTree->Branch("n",&hfRecHit.n,"n/I");
500		hfTree->Branch("e",hfRecHit.e,"e[n]/F");
501		hfTree->Branch("et",hfRecHit.et,"et[n]/F");
#	Line 453 | Line 504 | RecHitTreeProducer::beginJob()
504		hfTree->Branch("depth",hfRecHit.depth,"depth[n]/I");
505		hfTree->Branch("isjet",hfRecHit.isjet,"isjet[n]/O");
506
507	<	eeTree = fs->make<TTree>("ee","");
507	>	eeTree = fs->make<TTree>("ee",versionTag);
508		eeTree->Branch("n",&eeRecHit.n,"n/I");
509		eeTree->Branch("e",eeRecHit.e,"e[n]/F");
510		eeTree->Branch("et",eeRecHit.et,"et[n]/F");
#	Line 461 | Line 512 | RecHitTreeProducer::beginJob()
512		eeTree->Branch("phi",eeRecHit.phi,"phi[n]/F");
513		eeTree->Branch("isjet",eeRecHit.isjet,"isjet[n]/O");
514
515	<	ebTree = fs->make<TTree>("eb","");
515	>	ebTree = fs->make<TTree>("eb",versionTag);
516		ebTree->Branch("n",&ebRecHit.n,"n/I");
517		ebTree->Branch("e",ebRecHit.e,"e[n]/F");
518		ebTree->Branch("et",ebRecHit.et,"et[n]/F");
#	Line 469 | Line 520 | RecHitTreeProducer::beginJob()
520		ebTree->Branch("phi",ebRecHit.phi,"phi[n]/F");
521		ebTree->Branch("isjet",ebRecHit.isjet,"isjet[n]/O");
522
523	<	towerTree = fs->make<TTree>("tower","");
523	>	towerTree = fs->make<TTree>("tower",versionTag);
524		towerTree->Branch("n",&myTowers.n,"n/I");
525		towerTree->Branch("e",myTowers.e,"e[n]/F");
526		towerTree->Branch("et",myTowers.et,"et[n]/F");
527		towerTree->Branch("eta",myTowers.eta,"eta[n]/F");
528		towerTree->Branch("phi",myTowers.phi,"phi[n]/F");
529		towerTree->Branch("isjet",myTowers.isjet,"isjet[n]/O");
530	+	if (hasVtx_) {
531	+	towerTree->Branch("etvtx",myTowers.etvtx,"etvtx[n]/F");
532	+	towerTree->Branch("etavtx",myTowers.etavtx,"etavtx[n]/F");
533	+	towerTree->Branch("emEtVtx",myTowers.emEtVtx,"emEtVtx[n]/F");
534	+	towerTree->Branch("hadEtVtx",myTowers.hadEtVtx,"hadEtVtx[n]/F");
535	+	}
536
537
538		if(doBasicClusters_){
539	<	bcTree = fs->make<TTree>("bc","");
539	>	bcTree = fs->make<TTree>("bc",versionTag);
540		bcTree->Branch("n",&myBC.n,"n/I");
541		bcTree->Branch("e",myBC.e,"e[n]/F");
542		bcTree->Branch("et",myBC.et,"et[n]/F");
#	Line 492 | Line 549 | RecHitTreeProducer::beginJob()
549		}
550
551		if(doFastJet_){
552	<	bkgTree = fs->make<TTree>("bkg","");
552	>	bkgTree = fs->make<TTree>("bkg",versionTag);
553		bkgTree->Branch("n",&bkg.n,"n/I");
554		bkgTree->Branch("rho",bkg.rho,"rho[n]/F");
555		bkgTree->Branch("sigma",bkg.sigma,"sigma[n]/F");
556		}
557
558	<	nt = fs->make<TNtuple>("ntEvent","","nHF:nHFlong:nHFshort");
558	>	nt = fs->make<TNtuple>("ntEvent","","nHFplus:nHFminus:nHFlongPlus:nHFlongMinus:nHFshortPlus:nHFshortMinus");
559
560		}
561
#	Line 525 | Line 582 | double RecHitTreeProducer::getPhi(const
582		return et;
583		}
584
585	<
585	>	reco::Vertex::Point RecHitTreeProducer::getVtx(const edm::Event& ev)
586	>	{
587	>	ev.getByLabel(VtxSrc_,vtxs);
588	>	int greatestvtx = 0;
589	>	int nVertex = vtxs->size();
590	>
591	>	for (unsigned int i = 0 ; i< vtxs->size(); ++i){
592	>	unsigned int daughter = (*vtxs)[i].tracksSize();
593	>	if( daughter > (*vtxs)[greatestvtx].tracksSize()) greatestvtx = i;
594	>	//cout <<"Vertex: "<< (*vtxs)[i].position().z()<<" "<<daughter<<endl;
595	>	}
596	>
597	>	if(nVertex<=0){
598	>	return reco::Vertex::Point(-999,-999,-999);
599	>	}
600	>	return (*vtxs)[greatestvtx].position();
601	>	}
602
603		//define this as a plug-in
604		DEFINE_FWK_MODULE(RecHitTreeProducer);

Diff Legend

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