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.17 by yjlee, Thu Nov 3 16:12:42 2011 UTC

#	Line 12 | Line 12
12		*/
13		//
14		// Original Author:  Yetkin Yilmaz
15	+	// Modified: Frank Ma, Yen-Jie Lee
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 doHF_;
198	+	bool hasVtx_;
199
200		bool doFastJet_;
201
#	Line 206 | Line 225 | RecHitTreeProducer::RecHitTreeProducer(c
225		HcalRecHitHBHESrc_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHBHERecHitSrc",edm::InputTag("hbhereco"));
226		BCSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("BasicClusterSrc1",edm::InputTag("ecalRecHit","EcalRecHitsEB","RECO"));
227		TowerSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("towersSrc",edm::InputTag("towerMaker"));
228	<	JetSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("JetSrc",edm::InputTag("iterativeCone5CaloJets"));
228	>	VtxSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("vtxSrc",edm::InputTag("hiSelectedVertex"));
229	>	JetSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("JetSrc",edm::InputTag("iterativeConePu5CaloJets"));
230		useJets_ = iConfig.getUntrackedParameter<bool>("useJets",true);
231		doBasicClusters_ = iConfig.getUntrackedParameter<bool>("doBasicClusters",false);
232		doTowers_ = iConfig.getUntrackedParameter<bool>("doTowers",true);
233		doEcal_ = iConfig.getUntrackedParameter<bool>("doEcal",true);
234		doHcal_ = iConfig.getUntrackedParameter<bool>("doHcal",true);
235	+	doHF_ = iConfig.getUntrackedParameter<bool>("doHF",true);
236	+	hasVtx_ = iConfig.getUntrackedParameter<bool>("hasVtx",false);
237		doFastJet_ = iConfig.getUntrackedParameter<bool>("doFastJet",true);
238		FastJetTag_ = iConfig.getUntrackedParameter<edm::InputTag>("FastJetTag",edm::InputTag("kt4CaloJets"));
239		doEbyEonly_ = iConfig.getUntrackedParameter<bool>("doEbyEonly",false);
240		hfTowerThreshold_ = iConfig.getUntrackedParameter<double>("HFtowerMin",3.);
241	<	hfLongThreshold_ = iConfig.getUntrackedParameter<double>("HFlongMin",3.);
242	<	hfShortThreshold_ = iConfig.getUntrackedParameter<double>("HFshortMin",3.);
241	>	hfLongThreshold_ = iConfig.getUntrackedParameter<double>("HFlongMin",0.5);
242	>	hfShortThreshold_ = iConfig.getUntrackedParameter<double>("HFshortMin",0.85);
243	>	hbhePtMin_ = iConfig.getUntrackedParameter<double>("HBHETreePtMin",0);
244	>	hfPtMin_ = iConfig.getUntrackedParameter<double>("HFTreePtMin",0);
245	>	ebPtMin_ = iConfig.getUntrackedParameter<double>("EBTreePtMin",0);
246	>	eePtMin_ = iConfig.getUntrackedParameter<double>("EETreePtMin",0.);
247	>	towerPtMin_ = iConfig.getUntrackedParameter<double>("TowerTreePtMin",0.);
248		}
249
250
#	Line 247 | Line 274 | RecHitTreeProducer::analyze(const edm::E
274		myTowers.n = 0;
275		bkg.n = 0;
276
277	<	if(doEcal_){
278	<	ev.getByLabel(EBSrc_,ebHits);
279	<	ev.getByLabel(EESrc_,eeHits);
277	>	// get vertex
278	>	reco::Vertex::Point vtx;
279	>	if (hasVtx_) vtx = getVtx(ev);
280	>
281	>	if(doEcal_){
282	>	ev.getByLabel(EBSrc_,ebHits);
283	>	ev.getByLabel(EESrc_,eeHits);
284		}
285	+
286		if(doHcal_){
287	<	ev.getByLabel(HcalRecHitHFSrc_,hfHits);
288	<	ev.getByLabel(HcalRecHitHBHESrc_,hbheHits);
287	>	ev.getByLabel(HcalRecHitHBHESrc_,hbheHits);
288	>	}
289	>	if(doHF_){
290	>	ev.getByLabel(HcalRecHitHFSrc_,hfHits);
291		}
292	+
293		if(useJets_) {
294		ev.getByLabel(JetSrc_,jets);
295		}
#	Line 285 | Line 320 | RecHitTreeProducer::analyze(const edm::E
320		geo = pGeo.product();
321		}
322
323	<	int nHFlong = 0;
324	<	int nHFshort = 0;
325	<	int nHFtower = 0;
323	>	int nHFlongPlus = 0;
324	>	int nHFshortPlus = 0;
325	>	int nHFtowerPlus = 0;
326	>	int nHFlongMinus = 0;
327	>	int nHFshortMinus = 0;
328	>	int nHFtowerMinus = 0;
329	>
330	>
331
332	<	if(doHcal_){
332	>	if(doHF_){
333		for(unsigned int i = 0; i < hfHits->size(); ++i){
334		const HFRecHit & hit= (*hfHits)[i];
335		hfRecHit.e[hfRecHit.n] = hit.energy();
#	Line 298 | Line 338 | RecHitTreeProducer::analyze(const edm::E
338		hfRecHit.phi[hfRecHit.n] = getPhi(hit.id());
339		hfRecHit.isjet[hfRecHit.n] = false;
340		hfRecHit.depth[hfRecHit.n] = hit.id().depth();
341	<	if(hit.energy() > hfShortThreshold_ && hit.id().depth() != 1) nHFshort++;
342	<	if(hit.energy() > hfLongThreshold_ && hit.id().depth() == 1) nHFlong++;
341	>
342	>	if(hit.id().ieta() > 0){
343	>	if(hit.energy() > hfShortThreshold_ && hit.id().depth() != 1) nHFshortPlus++;
344	>	if(hit.energy() > hfLongThreshold_ && hit.id().depth() == 1) nHFlongPlus++;
345	>	}else{
346	>	if(hit.energy() > hfShortThreshold_ && hit.id().depth() != 1) nHFshortMinus++;
347	>	if(hit.energy() > hfLongThreshold_ && hit.id().depth() == 1) nHFlongMinus++;
348	>	}
349
350		if(useJets_){
351		for(unsigned int j = 0 ; j < jets->size(); ++j){
#	Line 308 | Line 354 | RecHitTreeProducer::analyze(const edm::E
354		if(dr < cone){ hfRecHit.isjet[hfRecHit.n] = true; }
355		}
356		}
357	<	hfRecHit.n++;
357	>	if (hfRecHit.et[hfRecHit.n]>=hfPtMin_) hfRecHit.n++;
358		}
359	<	if(!doEbyEonly_){
359	>
360	>	if(doHcal_ && !doEbyEonly_){
361		for(unsigned int i = 0; i < hbheHits->size(); ++i){
362		const HBHERecHit & hit= (*hbheHits)[i];
363	+	if (getEt(hit.id(),hit.energy())<hbhePtMin_) continue;
364		hbheRecHit.e[hbheRecHit.n] = hit.energy();
365		hbheRecHit.et[hbheRecHit.n] = getEt(hit.id(),hit.energy());
366		hbheRecHit.eta[hbheRecHit.n] = getEta(hit.id());
#	Line 332 | Line 380 | RecHitTreeProducer::analyze(const edm::E
380		if(doEcal_ && !doEbyEonly_){
381		for(unsigned int i = 0; i < ebHits->size(); ++i){
382		const EcalRecHit & hit= (*ebHits)[i];
383	+	if (getEt(hit.id(),hit.energy())<ebPtMin_) continue;
384		ebRecHit.e[ebRecHit.n] = hit.energy();
385		ebRecHit.et[ebRecHit.n] = getEt(hit.id(),hit.energy());
386		ebRecHit.eta[ebRecHit.n] = getEta(hit.id());
#	Line 349 | Line 398 | RecHitTreeProducer::analyze(const edm::E
398
399		for(unsigned int i = 0; i < eeHits->size(); ++i){
400		const EcalRecHit & hit= (*eeHits)[i];
401	+	if (getEt(hit.id(),hit.energy())<eePtMin_) continue;
402		eeRecHit.e[eeRecHit.n] = hit.energy();
403		eeRecHit.et[eeRecHit.n] = getEt(hit.id(),hit.energy());
404		eeRecHit.eta[eeRecHit.n] = getEta(hit.id());
#	Line 369 | Line 419 | RecHitTreeProducer::analyze(const edm::E
419
420		for(unsigned int i = 0; i < towers->size(); ++i){
421		const CaloTower & hit= (*towers)[i];
422	+	if (getEt(hit.id(),hit.energy())<towerPtMin_) continue;
423		myTowers.e[myTowers.n] = hit.energy();
424		myTowers.et[myTowers.n] = getEt(hit.id(),hit.energy());
425		myTowers.eta[myTowers.n] = getEta(hit.id());
426		myTowers.phi[myTowers.n] = getPhi(hit.id());
427		myTowers.isjet[myTowers.n] = false;
428	+	if (hasVtx_) {
429	+	myTowers.etvtx[myTowers.n] = hit.p4(vtx).Et();
430	+	myTowers.etavtx[myTowers.n] = hit.p4(vtx).Eta();
431	+	myTowers.emEtVtx[myTowers.n] = hit.emEt(vtx);
432	+	myTowers.hadEtVtx[myTowers.n] = hit.hadEt(vtx);
433	+	}
434
435	<	if(hit.ieta() > 29 && hit.energy() > hfTowerThreshold_) nHFtower++;
435	>	if(hit.ieta() > 29 && hit.energy() > hfTowerThreshold_) nHFtowerPlus++;
436	>	if(hit.ieta() < -29 && hit.energy() > hfTowerThreshold_) nHFtowerMinus++;
437
438		if(useJets_){
439		for(unsigned int j = 0 ; j < jets->size(); ++j){
#	Line 426 | Line 484 | RecHitTreeProducer::analyze(const edm::E
484		}
485		}
486
487	<	nt->Fill(nHFtower,nHFlong,nHFshort);
487	>	nt->Fill(nHFtowerPlus,nHFtowerMinus,nHFlongPlus,nHFlongMinus,nHFshortPlus,nHFshortMinus);
488
489		}
490
#	Line 436 | Line 494 | void
494		RecHitTreeProducer::beginJob()
495		{
496
497	<	hbheTree = fs->make<TTree>("hbhe","");
497	>	hbheTree = fs->make<TTree>("hbhe",versionTag);
498		hbheTree->Branch("n",&hbheRecHit.n,"n/I");
499		hbheTree->Branch("e",hbheRecHit.e,"e[n]/F");
500		hbheTree->Branch("et",hbheRecHit.et,"et[n]/F");
#	Line 444 | Line 502 | RecHitTreeProducer::beginJob()
502		hbheTree->Branch("phi",hbheRecHit.phi,"phi[n]/F");
503		hbheTree->Branch("isjet",hbheRecHit.isjet,"isjet[n]/O");
504
505	<	hfTree = fs->make<TTree>("hf","");
505	>	hfTree = fs->make<TTree>("hf",versionTag);
506		hfTree->Branch("n",&hfRecHit.n,"n/I");
507		hfTree->Branch("e",hfRecHit.e,"e[n]/F");
508		hfTree->Branch("et",hfRecHit.et,"et[n]/F");
#	Line 453 | Line 511 | RecHitTreeProducer::beginJob()
511		hfTree->Branch("depth",hfRecHit.depth,"depth[n]/I");
512		hfTree->Branch("isjet",hfRecHit.isjet,"isjet[n]/O");
513
514	<	eeTree = fs->make<TTree>("ee","");
514	>	eeTree = fs->make<TTree>("ee",versionTag);
515		eeTree->Branch("n",&eeRecHit.n,"n/I");
516		eeTree->Branch("e",eeRecHit.e,"e[n]/F");
517		eeTree->Branch("et",eeRecHit.et,"et[n]/F");
#	Line 461 | Line 519 | RecHitTreeProducer::beginJob()
519		eeTree->Branch("phi",eeRecHit.phi,"phi[n]/F");
520		eeTree->Branch("isjet",eeRecHit.isjet,"isjet[n]/O");
521
522	<	ebTree = fs->make<TTree>("eb","");
522	>	ebTree = fs->make<TTree>("eb",versionTag);
523		ebTree->Branch("n",&ebRecHit.n,"n/I");
524		ebTree->Branch("e",ebRecHit.e,"e[n]/F");
525		ebTree->Branch("et",ebRecHit.et,"et[n]/F");
#	Line 469 | Line 527 | RecHitTreeProducer::beginJob()
527		ebTree->Branch("phi",ebRecHit.phi,"phi[n]/F");
528		ebTree->Branch("isjet",ebRecHit.isjet,"isjet[n]/O");
529
530	<	towerTree = fs->make<TTree>("tower","");
530	>	towerTree = fs->make<TTree>("tower",versionTag);
531		towerTree->Branch("n",&myTowers.n,"n/I");
532		towerTree->Branch("e",myTowers.e,"e[n]/F");
533		towerTree->Branch("et",myTowers.et,"et[n]/F");
534		towerTree->Branch("eta",myTowers.eta,"eta[n]/F");
535		towerTree->Branch("phi",myTowers.phi,"phi[n]/F");
536		towerTree->Branch("isjet",myTowers.isjet,"isjet[n]/O");
537	+	if (hasVtx_) {
538	+	towerTree->Branch("etvtx",myTowers.etvtx,"etvtx[n]/F");
539	+	towerTree->Branch("etavtx",myTowers.etavtx,"etavtx[n]/F");
540	+	towerTree->Branch("emEtVtx",myTowers.emEtVtx,"emEtVtx[n]/F");
541	+	towerTree->Branch("hadEtVtx",myTowers.hadEtVtx,"hadEtVtx[n]/F");
542	+	}
543
544
545		if(doBasicClusters_){
546	<	bcTree = fs->make<TTree>("bc","");
546	>	bcTree = fs->make<TTree>("bc",versionTag);
547		bcTree->Branch("n",&myBC.n,"n/I");
548		bcTree->Branch("e",myBC.e,"e[n]/F");
549		bcTree->Branch("et",myBC.et,"et[n]/F");
#	Line 492 | Line 556 | RecHitTreeProducer::beginJob()
556		}
557
558		if(doFastJet_){
559	<	bkgTree = fs->make<TTree>("bkg","");
559	>	bkgTree = fs->make<TTree>("bkg",versionTag);
560		bkgTree->Branch("n",&bkg.n,"n/I");
561		bkgTree->Branch("rho",bkg.rho,"rho[n]/F");
562		bkgTree->Branch("sigma",bkg.sigma,"sigma[n]/F");
563		}
564
565	<	nt = fs->make<TNtuple>("ntEvent","","nHF:nHFlong:nHFshort");
565	>	nt = fs->make<TNtuple>("ntEvent","","nHFplus:nHFminus:nHFlongPlus:nHFlongMinus:nHFshortPlus:nHFshortMinus");
566
567		}
568
#	Line 525 | Line 589 | double RecHitTreeProducer::getPhi(const
589		return et;
590		}
591
592	<
592	>	reco::Vertex::Point RecHitTreeProducer::getVtx(const edm::Event& ev)
593	>	{
594	>	ev.getByLabel(VtxSrc_,vtxs);
595	>	int greatestvtx = 0;
596	>	int nVertex = vtxs->size();
597	>
598	>	for (unsigned int i = 0 ; i< vtxs->size(); ++i){
599	>	unsigned int daughter = (*vtxs)[i].tracksSize();
600	>	if( daughter > (*vtxs)[greatestvtx].tracksSize()) greatestvtx = i;
601	>	//cout <<"Vertex: "<< (*vtxs)[i].position().z()<<" "<<daughter<<endl;
602	>	}
603	>
604	>	if(nVertex<=0){
605	>	return reco::Vertex::Point(-999,-999,-999);
606	>	}
607	>	return (*vtxs)[greatestvtx].position();
608	>	}
609
610		//define this as a plug-in
611		DEFINE_FWK_MODULE(RecHitTreeProducer);

Diff Legend

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