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.17 by yjlee, Thu Nov 3 16:12:42 2011 UTC vs.
Revision 1.22 by yilmaz, Wed Jan 16 17:22:56 2013 UTC

#	Line 69 | Line 69
69
70		using namespace std;
71
72	<	#define MAXHITS 1000000
72	>	#define MAXHITS 100000
73
74		struct MyRecHit{
75		int depth[MAXHITS];
76		int n;
77
78	+	int ieta[MAXHITS];
79	+	int iphi[MAXHITS];
80	+
81		float e[MAXHITS];
82		float et[MAXHITS];
83		float eta[MAXHITS];
84		float phi[MAXHITS];
85	+	float perp[MAXHITS];
86	+	float emEt[MAXHITS];
87	+	float hadEt[MAXHITS];
88	+
89		bool isjet[MAXHITS];
90	<	float etvtx[MAXHITS];
91	<	float etavtx[MAXHITS];
90	>	float etVtx[MAXHITS];
91	>	float etaVtx[MAXHITS];
92	>	float phiVtx[MAXHITS];
93	>	float perpVtx[MAXHITS];
94		float emEtVtx[MAXHITS];
95		float hadEtVtx[MAXHITS];
96
#	Line 107 | Line 116 | class RecHitTreeProducer : public edm::E
116		public:
117		explicit RecHitTreeProducer(const edm::ParameterSet&);
118		~RecHitTreeProducer();
119	<	double       getEt(const DetId &id, double energy);
120	<	double       getEta(const DetId &id);
121	<	double       getPhi(const DetId &id);
119	>	double       getEt(const DetId &id, double energy, reco::Vertex::Point vtx = reco::Vertex::Point(0,0,0));
120	>	double       getEta(const DetId &id, reco::Vertex::Point vtx = reco::Vertex::Point(0,0,0));
121	>	double       getPhi(const DetId &id, reco::Vertex::Point vtx = reco::Vertex::Point(0,0,0));
122	>	double       getPerp(const DetId &id, reco::Vertex::Point vtx = reco::Vertex::Point(0,0,0));
123	>
124	>	math::XYZPoint getPosition(const DetId &id, reco::Vertex::Point vtx = reco::Vertex::Point(0,0,0));
125	>	double getEt(math::XYZPoint pos, double energy);
126	>
127		reco::Vertex::Point getVtx(const edm::Event& ev);
128
129
#	Line 148 | Line 162 | class RecHitTreeProducer : public edm::E
162		MyRecHit eeRecHit;
163		MyRecHit myBC;
164		MyRecHit myTowers;
165	+	MyRecHit castorRecHit;
166	+
167		MyBkg bkg;
168
169		TNtuple* nt;
#	Line 158 | Line 174 | class RecHitTreeProducer : public edm::E
174		TTree* bcTree;
175		TTree* towerTree;
176		TTree* bkgTree;
177	+	TTree* castorTree;
178
179		double cone;
180		double hfTowerThreshold_;
#	Line 195 | Line 212 | class RecHitTreeProducer : public edm::E
212		bool doEcal_;
213		bool doHcal_;
214		bool doHF_;
215	+	bool doCastor_;
216	+
217		bool hasVtx_;
218	+	bool saveBothVtx_;
219
220		bool doFastJet_;
221
222		bool doEbyEonly_;
223	+
224		};
225
226		//
#	Line 233 | Line 254 | RecHitTreeProducer::RecHitTreeProducer(c
254		doEcal_ = iConfig.getUntrackedParameter<bool>("doEcal",true);
255		doHcal_ = iConfig.getUntrackedParameter<bool>("doHcal",true);
256		doHF_ = iConfig.getUntrackedParameter<bool>("doHF",true);
257	<	hasVtx_ = iConfig.getUntrackedParameter<bool>("hasVtx",false);
257	>	doCastor_ = iConfig.getUntrackedParameter<bool>("doCASTOR",true);
258	>
259	>	hasVtx_ = iConfig.getUntrackedParameter<bool>("hasVtx",true);
260	>	saveBothVtx_ = iConfig.getUntrackedParameter<bool>("saveBothVtx",false);
261	>
262		doFastJet_ = iConfig.getUntrackedParameter<bool>("doFastJet",true);
263		FastJetTag_ = iConfig.getUntrackedParameter<edm::InputTag>("FastJetTag",edm::InputTag("kt4CaloJets"));
264		doEbyEonly_ = iConfig.getUntrackedParameter<bool>("doEbyEonly",false);
#	Line 245 | Line 270 | RecHitTreeProducer::RecHitTreeProducer(c
270		ebPtMin_ = iConfig.getUntrackedParameter<double>("EBTreePtMin",0);
271		eePtMin_ = iConfig.getUntrackedParameter<double>("EETreePtMin",0.);
272		towerPtMin_ = iConfig.getUntrackedParameter<double>("TowerTreePtMin",0.);
273	+
274	+
275		}
276
277
#	Line 275 | Line 302 | RecHitTreeProducer::analyze(const edm::E
302		bkg.n = 0;
303
304		// get vertex
305	<	reco::Vertex::Point vtx;
305	>	reco::Vertex::Point vtx(0,0,0);
306		if (hasVtx_) vtx = getVtx(ev);
307
308		if(doEcal_){
#	Line 326 | Line 353 | RecHitTreeProducer::analyze(const edm::E
353		int nHFlongMinus = 0;
354		int nHFshortMinus = 0;
355		int nHFtowerMinus = 0;
356	<
330	<
331	<
356	>
357		if(doHF_){
358		for(unsigned int i = 0; i < hfHits->size(); ++i){
359		const HFRecHit & hit= (*hfHits)[i];
360		hfRecHit.e[hfRecHit.n] = hit.energy();
361	<	hfRecHit.et[hfRecHit.n] = getEt(hit.id(),hit.energy());
362	<	hfRecHit.eta[hfRecHit.n] = getEta(hit.id());
363	<	hfRecHit.phi[hfRecHit.n] = getPhi(hit.id());
361	>	math::XYZPoint pos = getPosition(hit.id(),vtx);
362	>
363	>	if(!saveBothVtx_){
364	>	hfRecHit.et[hfRecHit.n] = getEt(pos,hit.energy());
365	>	hfRecHit.eta[hfRecHit.n] = pos.eta();
366	>	hfRecHit.phi[hfRecHit.n] = pos.phi();
367	>	hfRecHit.perp[hfRecHit.n] = pos.rho();
368	>	}else{
369	>	hfRecHit.et[hfRecHit.n] = getEt(hit.id(),hit.energy());
370	>	hfRecHit.eta[hfRecHit.n] = getEta(hit.id());
371	>	hfRecHit.phi[hfRecHit.n] = getPhi(hit.id());
372	>	hfRecHit.perp[hfRecHit.n] = getPerp(hit.id());
373	>
374	>	hfRecHit.etVtx[hfRecHit.n] = getEt(pos,hit.energy());
375	>	hfRecHit.etaVtx[hfRecHit.n] = pos.eta();
376	>	hfRecHit.phiVtx[hfRecHit.n] = pos.phi();
377	>	hfRecHit.perpVtx[hfRecHit.n] = pos.rho();
378	>
379	>	}
380	>
381		hfRecHit.isjet[hfRecHit.n] = false;
382		hfRecHit.depth[hfRecHit.n] = hit.id().depth();
383
#	Line 361 | Line 403 | RecHitTreeProducer::analyze(const edm::E
403		for(unsigned int i = 0; i < hbheHits->size(); ++i){
404		const HBHERecHit & hit= (*hbheHits)[i];
405		if (getEt(hit.id(),hit.energy())<hbhePtMin_) continue;
406	+
407		hbheRecHit.e[hbheRecHit.n] = hit.energy();
408	<	hbheRecHit.et[hbheRecHit.n] = getEt(hit.id(),hit.energy());
409	<	hbheRecHit.eta[hbheRecHit.n] = getEta(hit.id());
410	<	hbheRecHit.phi[hbheRecHit.n] = getPhi(hit.id());
411	<	hbheRecHit.isjet[hbheRecHit.n] = false;
408	>	math::XYZPoint pos = getPosition(hit.id(),vtx);
409	>
410	>	if(!saveBothVtx_){
411	>	hbheRecHit.et[hbheRecHit.n] = getEt(pos,hit.energy());
412	>	hbheRecHit.eta[hbheRecHit.n] = pos.eta();
413	>	hbheRecHit.phi[hbheRecHit.n] = pos.phi();
414	>	hbheRecHit.perp[hbheRecHit.n] = pos.rho();
415	>	}else{
416	>	hbheRecHit.et[hbheRecHit.n] = getEt(hit.id(),hit.energy());
417	>	hbheRecHit.eta[hbheRecHit.n] = getEta(hit.id());
418	>	hbheRecHit.phi[hbheRecHit.n] = getPhi(hit.id());
419	>	hbheRecHit.perp[hbheRecHit.n] = getPerp(hit.id());
420	>
421	>	hbheRecHit.etVtx[hbheRecHit.n] = getEt(pos,hit.energy());
422	>	hbheRecHit.etaVtx[hbheRecHit.n] = pos.eta();
423	>	hbheRecHit.phiVtx[hbheRecHit.n] = pos.phi();
424	>	hbheRecHit.perpVtx[hbheRecHit.n] = pos.rho();
425	>
426	>	}
427	>
428	>	hbheRecHit.isjet[hbheRecHit.n] = false;
429	>	hbheRecHit.depth[hbheRecHit.n] = hit.id().depth();
430	>
431		if(useJets_){
432		for(unsigned int j = 0 ; j < jets->size(); ++j){
433		const reco::Jet& jet = (*jets)[j];
#	Line 381 | Line 443 | RecHitTreeProducer::analyze(const edm::E
443		for(unsigned int i = 0; i < ebHits->size(); ++i){
444		const EcalRecHit & hit= (*ebHits)[i];
445		if (getEt(hit.id(),hit.energy())<ebPtMin_) continue;
446	+
447		ebRecHit.e[ebRecHit.n] = hit.energy();
448	<	ebRecHit.et[ebRecHit.n] = getEt(hit.id(),hit.energy());
449	<	ebRecHit.eta[ebRecHit.n] = getEta(hit.id());
450	<	ebRecHit.phi[ebRecHit.n] = getPhi(hit.id());
448	>	math::XYZPoint pos = getPosition(hit.id(),vtx);
449	>
450	>	if(!saveBothVtx_){
451	>	ebRecHit.et[ebRecHit.n] = getEt(pos,hit.energy());
452	>	ebRecHit.eta[ebRecHit.n] = pos.eta();
453	>	ebRecHit.phi[ebRecHit.n] = pos.phi();
454	>	ebRecHit.perp[ebRecHit.n] = pos.rho();
455	>	}else{
456	>	ebRecHit.et[ebRecHit.n] = getEt(hit.id(),hit.energy());
457	>	ebRecHit.eta[ebRecHit.n] = getEta(hit.id());
458	>	ebRecHit.phi[ebRecHit.n] = getPhi(hit.id());
459	>	ebRecHit.perp[ebRecHit.n] = getPerp(hit.id());
460	>	ebRecHit.etVtx[ebRecHit.n] = getEt(pos,hit.energy());
461	>	ebRecHit.etaVtx[ebRecHit.n] = pos.eta();
462	>	ebRecHit.phiVtx[ebRecHit.n] = pos.phi();
463	>	ebRecHit.perpVtx[ebRecHit.n] = pos.rho();
464	>
465	>	}
466	>
467		ebRecHit.isjet[ebRecHit.n] = false;
468		if(useJets_){
469		for(unsigned int j = 0 ; j < jets->size(); ++j){
#	Line 399 | Line 478 | RecHitTreeProducer::analyze(const edm::E
478		for(unsigned int i = 0; i < eeHits->size(); ++i){
479		const EcalRecHit & hit= (*eeHits)[i];
480		if (getEt(hit.id(),hit.energy())<eePtMin_) continue;
481	+
482		eeRecHit.e[eeRecHit.n] = hit.energy();
483	<	eeRecHit.et[eeRecHit.n] = getEt(hit.id(),hit.energy());
484	<	eeRecHit.eta[eeRecHit.n] = getEta(hit.id());
485	<	eeRecHit.phi[eeRecHit.n] = getPhi(hit.id());
483	>	math::XYZPoint pos = getPosition(hit.id(),vtx);
484	>
485	>	if(!saveBothVtx_){
486	>	eeRecHit.et[eeRecHit.n] = getEt(pos,hit.energy());
487	>	eeRecHit.eta[eeRecHit.n] = pos.eta();
488	>	eeRecHit.phi[eeRecHit.n] = pos.phi();
489	>	eeRecHit.perp[eeRecHit.n] = pos.rho();
490	>	}else{
491	>	eeRecHit.et[eeRecHit.n] = getEt(hit.id(),hit.energy());
492	>	eeRecHit.eta[eeRecHit.n] = getEta(hit.id());
493	>	eeRecHit.phi[eeRecHit.n] = getPhi(hit.id());
494	>	eeRecHit.perp[eeRecHit.n] = getPerp(hit.id());
495	>	eeRecHit.etVtx[eeRecHit.n] = getEt(pos,hit.energy());
496	>	eeRecHit.etaVtx[eeRecHit.n] = pos.eta();
497	>	eeRecHit.phiVtx[eeRecHit.n] = pos.phi();
498	>	eeRecHit.perpVtx[eeRecHit.n] = pos.rho();
499	>
500	>	}
501	>
502		eeRecHit.isjet[eeRecHit.n] = false;
503	+
504		if(useJets_){
505		for(unsigned int j = 0 ; j < jets->size(); ++j){
506		const reco::Jet& jet = (*jets)[j];
#	Line 420 | Line 517 | RecHitTreeProducer::analyze(const edm::E
517		for(unsigned int i = 0; i < towers->size(); ++i){
518		const CaloTower & hit= (*towers)[i];
519		if (getEt(hit.id(),hit.energy())<towerPtMin_) continue;
520	<	myTowers.e[myTowers.n] = hit.energy();
521	<	myTowers.et[myTowers.n] = getEt(hit.id(),hit.energy());
522	<	myTowers.eta[myTowers.n] = getEta(hit.id());
523	<	myTowers.phi[myTowers.n] = getPhi(hit.id());
524	<	myTowers.isjet[myTowers.n] = false;
525	<	if (hasVtx_) {
526	<	myTowers.etvtx[myTowers.n] = hit.p4(vtx).Et();
527	<	myTowers.etavtx[myTowers.n] = hit.p4(vtx).Eta();
520	>
521	>	myTowers.et[myTowers.n] = hit.p4(vtx).Et();
522	>	myTowers.eta[myTowers.n] = hit.p4(vtx).Eta();
523	>	myTowers.phi[myTowers.n] = hit.p4(vtx).Phi();
524	>	myTowers.emEt[myTowers.n] = hit.emEt(vtx);
525	>	myTowers.hadEt[myTowers.n] = hit.hadEt(vtx);
526	>
527	>	if (saveBothVtx_) {
528	>	myTowers.e[myTowers.n] = hit.energy();
529	>	myTowers.et[myTowers.n] = getEt(hit.id(),hit.energy());
530	>	myTowers.eta[myTowers.n] = getEta(hit.id());
531	>	myTowers.phi[myTowers.n] = getPhi(hit.id());
532	>	myTowers.isjet[myTowers.n] = false;
533	>	myTowers.etVtx[myTowers.n] = hit.p4(vtx).Et();
534	>	myTowers.etaVtx[myTowers.n] = hit.p4(vtx).Eta();
535		myTowers.emEtVtx[myTowers.n] = hit.emEt(vtx);
536		myTowers.hadEtVtx[myTowers.n] = hit.hadEt(vtx);
537		}
538
539	+	myTowers.isjet[myTowers.n] = false;
540	+
541		if(hit.ieta() > 29 && hit.energy() > hfTowerThreshold_) nHFtowerPlus++;
542		if(hit.ieta() < -29 && hit.energy() > hfTowerThreshold_) nHFtowerMinus++;
543
#	Line 470 | Line 576 | RecHitTreeProducer::analyze(const edm::E
576		}
577		}
578
579	+	if(doCastor_){
580	+
581	+	edm::Handle<CastorRecHitCollection> casrechits;
582	+	try{ ev.getByLabel("castorreco",casrechits); }
583	+	catch(...) { edm::LogWarning(" CASTOR ") << " Cannot get Castor RecHits " << std::endl; }
584	+
585	+	int nhits = 0;
586	+	double energyCastor = 0;
587	+
588	+	if(casrechits.failedToGet()!=0 || !casrechits.isValid()) {
589	+	edm::LogWarning(" CASTOR ") << " Cannot read CastorRecHitCollection" << std::endl;
590	+	} else {
591	+	for(size_t i1 = 0; i1 < casrechits->size(); ++i1) {
592	+	const CastorRecHit & rh = (*casrechits)[i1];
593	+	HcalCastorDetId castorid = rh.id();
594	+	energyCastor += rh.energy();
595	+	if (nhits  < 224) {
596	+	castorRecHit.e[nhits] = rh.energy();
597	+	castorRecHit.iphi[nhits] = castorid.sector();
598	+	castorRecHit.depth[nhits] = castorid.module();
599	+	castorRecHit.phi[nhits] = getPhi(castorid);
600	+	}
601	+
602	+	nhits++;
603	+
604	+	} // end loop castor rechits
605	+	}
606	+
607	+	castorRecHit.n = nhits;
608	+	castorTree->Fill();
609	+	}
610	+
611	+
612		if(!doEbyEonly_){
613		towerTree->Fill();
614
#	Line 500 | Line 639 | RecHitTreeProducer::beginJob()
639		hbheTree->Branch("et",hbheRecHit.et,"et[n]/F");
640		hbheTree->Branch("eta",hbheRecHit.eta,"eta[n]/F");
641		hbheTree->Branch("phi",hbheRecHit.phi,"phi[n]/F");
642	+	hbheTree->Branch("perp",hbheRecHit.perp,"perp[n]/F");
643	+
644		hbheTree->Branch("isjet",hbheRecHit.isjet,"isjet[n]/O");
645	<
645	>	hbheTree->Branch("depth",hfRecHit.depth,"depth[n]/I");
646	>
647		hfTree = fs->make<TTree>("hf",versionTag);
648		hfTree->Branch("n",&hfRecHit.n,"n/I");
649		hfTree->Branch("e",hfRecHit.e,"e[n]/F");
650		hfTree->Branch("et",hfRecHit.et,"et[n]/F");
651		hfTree->Branch("eta",hfRecHit.eta,"eta[n]/F");
652		hfTree->Branch("phi",hfRecHit.phi,"phi[n]/F");
653	+	hfTree->Branch("perp",hfRecHit.perp,"perp[n]/F");
654		hfTree->Branch("depth",hfRecHit.depth,"depth[n]/I");
655		hfTree->Branch("isjet",hfRecHit.isjet,"isjet[n]/O");
656
#	Line 517 | Line 660 | RecHitTreeProducer::beginJob()
660		eeTree->Branch("et",eeRecHit.et,"et[n]/F");
661		eeTree->Branch("eta",eeRecHit.eta,"eta[n]/F");
662		eeTree->Branch("phi",eeRecHit.phi,"phi[n]/F");
663	+	eeTree->Branch("perp",eeRecHit.perp,"perp[n]/F");
664	+
665		eeTree->Branch("isjet",eeRecHit.isjet,"isjet[n]/O");
666
667		ebTree = fs->make<TTree>("eb",versionTag);
#	Line 525 | Line 670 | RecHitTreeProducer::beginJob()
670		ebTree->Branch("et",ebRecHit.et,"et[n]/F");
671		ebTree->Branch("eta",ebRecHit.eta,"eta[n]/F");
672		ebTree->Branch("phi",ebRecHit.phi,"phi[n]/F");
673	+	ebTree->Branch("perp",ebRecHit.perp,"perp[n]/F");
674	+
675		ebTree->Branch("isjet",ebRecHit.isjet,"isjet[n]/O");
676
677		towerTree = fs->make<TTree>("tower",versionTag);
#	Line 534 | Line 681 | RecHitTreeProducer::beginJob()
681		towerTree->Branch("eta",myTowers.eta,"eta[n]/F");
682		towerTree->Branch("phi",myTowers.phi,"phi[n]/F");
683		towerTree->Branch("isjet",myTowers.isjet,"isjet[n]/O");
684	<	if (hasVtx_) {
685	<	towerTree->Branch("etvtx",myTowers.etvtx,"etvtx[n]/F");
686	<	towerTree->Branch("etavtx",myTowers.etavtx,"etavtx[n]/F");
684	>	towerTree->Branch("emEt",myTowers.emEt,"emEt[n]/F");
685	>	towerTree->Branch("hadEt",myTowers.hadEt,"hadEt[n]/F");
686	>
687	>
688	>	if(doCastor_){
689	>	castorTree = fs->make<TTree>("castor",versionTag);
690	>	castorTree->Branch("n",&castorRecHit.n,"n/I");
691	>	castorTree->Branch("e",castorRecHit.e,"e[n]/F");
692	>	castorTree->Branch("iphi",castorRecHit.iphi,"iphi[n]/I");
693	>	castorTree->Branch("phi",castorRecHit.phi,"phi[n]/F");
694	>	castorTree->Branch("depth",castorRecHit.depth,"depth[n]/I");
695	>	}
696	>
697	>	if (saveBothVtx_) {
698	>	towerTree->Branch("etVtx",myTowers.etVtx,"etvtx[n]/F");
699	>	towerTree->Branch("etaVtx",myTowers.etaVtx,"etavtx[n]/F");
700		towerTree->Branch("emEtVtx",myTowers.emEtVtx,"emEtVtx[n]/F");
701		towerTree->Branch("hadEtVtx",myTowers.hadEtVtx,"hadEtVtx[n]/F");
702		}
#	Line 571 | Line 731 | void
731		RecHitTreeProducer::endJob() {
732		}
733
734	<	double RecHitTreeProducer::getEt(const DetId &id, double energy){
734	>	math::XYZPoint RecHitTreeProducer::getPosition(const DetId &id, reco::Vertex::Point vtx){
735		const GlobalPoint& pos=geo->getPosition(id);
736	+	math::XYZPoint posV(pos.x() - vtx.x(),pos.y() - vtx.y(),pos.z() - vtx.z());
737	+	return posV;
738	+	}
739	+
740	+	double RecHitTreeProducer::getEt(math::XYZPoint pos, double energy){
741		double et = energy*sin(pos.theta());
742		return et;
743		}
744
745	<	double RecHitTreeProducer::getEta(const DetId &id){
745	>	double RecHitTreeProducer::getEt(const DetId &id, double energy, reco::Vertex::Point vtx){
746	>	const GlobalPoint& pos=geo->getPosition(id);
747	>	double et = energy*sin(pos.theta());
748	>	return et;
749	>	}
750	>
751	>	double RecHitTreeProducer::getEta(const DetId &id, reco::Vertex::Point vtx){
752		const GlobalPoint& pos=geo->getPosition(id);
753		double et = pos.eta();
754		return et;
755		}
756
757	<	double RecHitTreeProducer::getPhi(const DetId &id){
757	>	double RecHitTreeProducer::getPhi(const DetId &id, reco::Vertex::Point vtx){
758		const GlobalPoint& pos=geo->getPosition(id);
759		double et = pos.phi();
760		return et;
761		}
762
763	+	double RecHitTreeProducer::getPerp(const DetId &id, reco::Vertex::Point vtx){
764	+	const GlobalPoint& pos=geo->getPosition(id);
765	+	double et = pos.perp();
766	+	return et;
767	+	}
768	+
769		reco::Vertex::Point RecHitTreeProducer::getVtx(const edm::Event& ev)
770		{
771		ev.getByLabel(VtxSrc_,vtxs);
#	Line 602 | Line 779 | reco::Vertex::Point RecHitTreeProducer::
779		}
780
781		if(nVertex<=0){
782	<	return reco::Vertex::Point(-999,-999,-999);
782	>	return reco::Vertex::Point(0,0,0);
783		}
784		return (*vtxs)[greatestvtx].position();
785		}

Diff Legend

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