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.1 by yilmaz, Sat Oct 16 13:46:13 2010 UTC vs.
Revision 1.6 by yilmaz, Mon Nov 1 21:48:31 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
#	Line 74 | Line 75 | struct MyRecHit{
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
#	Line 99 | Line 105 | class RecHitTreeProducer : public edm::E
105
106		// ----------member data ---------------------------
107
102	–
108		edm::Handle<reco::Centrality> cent;
109		edm::Handle<vector<double> > ktRhos;
110		edm::Handle<vector<double> > akRhos;
111
112	<	edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > jets1;
113	<	edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > jets2;
112	>	edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > ebHits;
113	>	edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > eeHits;
114
115	<	//   typedef edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > >::const_iterator EcalIterator;
116	<	typedef vector<EcalRecHit>::const_iterator EcalIterator;
115	>	edm::Handle<HFRecHitCollection> hfHits;
116	>	edm::Handle<HBHERecHitCollection> hbheHits;
117
118	<	edm::Handle<reco::CaloJetCollection> signalJets;
118	>	edm::Handle<reco::BasicClusterCollection> bClusters;
119	>	edm::Handle<CaloTowerCollection> towers;
120
121	+
122	+	typedef vector<EcalRecHit>::const_iterator EcalIterator;
123	+
124	+	edm::Handle<reco::CaloJetCollection> jets;
125	+
126		MyRecHit hbheRecHit;
127		MyRecHit hfRecHit;
128	+	MyRecHit ebRecHit;
129	+	MyRecHit eeRecHit;
130	+	MyRecHit myBC;
131	+	MyRecHit myTowers;
132
133		TTree* hbheTree;
134		TTree* hfTree;
135	+	TTree* ebTree;
136	+	TTree* eeTree;
137	+	TTree* bcTree;
138	+	TTree* towerTree;
139
140	<	double cone;
140	>	double cone;
141
142		edm::Service<TFileService> fs;
143		const CentralityBins * cbins_;
#	Line 126 | Line 145 | class RecHitTreeProducer : public edm::E
145
146		edm::InputTag HcalRecHitHFSrc_;
147		edm::InputTag HcalRecHitHBHESrc_;
148	+	edm::InputTag EBSrc_;
149	+	edm::InputTag EESrc_;
150	+	edm::InputTag BCSrc_;
151	+
152	+	edm::InputTag TowerSrc_;
153	+
154	+	edm::InputTag JetSrc_;
155
156	+	bool useJets_;
157	+	bool doBasicClusters_;
158
159		};
160
#	Line 147 | Line 175 | RecHitTreeProducer::RecHitTreeProducer(c
175		cone(0.5)
176		{
177		//now do what ever initialization is needed
178	+	EBSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("EBRecHitSrc",edm::InputTag("ecalRecHit","EcalRecHitsEB"));
179	+	EESrc_ = iConfig.getUntrackedParameter<edm::InputTag>("EERecHitSrc",edm::InputTag("ecalRecHit","EcalRecHitsEE"));
180		HcalRecHitHFSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHFRecHitSrc",edm::InputTag("hfreco"));
181		HcalRecHitHBHESrc_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHBHERecHitSrc",edm::InputTag("hbhereco"));
182	+	BCSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("BasicClusterSrc1",edm::InputTag("ecalRecHit","EcalRecHitsEB","RECO"));
183	+	TowerSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("towersSrc",edm::InputTag("towerMaker"));
184	+	JetSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("JetSrc",edm::InputTag("iterativeCone5CaloJets"));
185	+	useJets_ = iConfig.getUntrackedParameter<bool>("useJets",false);
186	+	doBasicClusters_ = iConfig.getUntrackedParameter<bool>("doBasicClusters",false);
187	+
188		}
189
190
#	Line 172 | Line 208 | RecHitTreeProducer::analyze(const edm::E
208
209		hfRecHit.n = 0;
210		hbheRecHit.n = 0;
211	+	ebRecHit.n = 0;
212	+	eeRecHit.n = 0;
213
214	<	edm::Handle<HFRecHitCollection> hfHits;
215	<	edm::Handle<HBHERecHitCollection> hbheHits;
214	>	ev.getByLabel(EBSrc_,ebHits);
215	>	ev.getByLabel(EESrc_,eeHits);
216
217		ev.getByLabel(HcalRecHitHFSrc_,hfHits);
218		ev.getByLabel(HcalRecHitHBHESrc_,hbheHits);
219
220	+	if(useJets_) {
221	+	ev.getByLabel(JetSrc_,jets);
222	+	}
223	+
224	+	if(doBasicClusters_){
225	+	ev.getByLabel(BCSrc_,bClusters);
226	+	}
227
228	<	if(0 && !cbins_) cbins_ = getCentralityBinsFromDB(iSetup);
228	>
229	>
230	>	if(0 && !cbins_) cbins_ = getCentralityBinsFromDB(iSetup);
231
232		if(!geo){
233		edm::ESHandle<CaloGeometry> pGeo;
234		iSetup.get<CaloGeometryRecord>().get(pGeo);
235		geo = pGeo.product();
236		}
237	<
237	>
238		for(unsigned int i = 0; i < hfHits->size(); ++i){
239		const HFRecHit & hit= (*hfHits)[i];
240		hfRecHit.e[hfRecHit.n] = hit.energy();
241		hfRecHit.et[hfRecHit.n] = getEt(hit.id(),hit.energy());
242		hfRecHit.eta[hfRecHit.n] = getEta(hit.id());
243		hfRecHit.phi[hfRecHit.n] = getPhi(hit.id());
244	<	//     hfRecHit.ieta[hfRecHit.n] = hit.id().ieta();
245	<	//     hfRecHit.iphi[hfRecHit.n] = hit.id().iphi();
246	<
244	>	hfRecHit.isjet[hfRecHit.n] = false;
245	>	if(!useJets_){
246	>	for(unsigned int j = 0 ; j < jets->size(); ++j){
247	>	const reco::Jet& jet = (*jets)[j];
248	>	double dr = reco::deltaR(hfRecHit.eta[hfRecHit.n],hfRecHit.phi[hfRecHit.n],jet.eta(),jet.phi());
249	>	if(dr < cone){ hfRecHit.isjet[hfRecHit.n] = true; }
250	>	}
251	>	}
252		hfRecHit.n++;
253		}
254	<
254	>
255		for(unsigned int i = 0; i < hbheHits->size(); ++i){
256		const HBHERecHit & hit= (*hbheHits)[i];
257	<	hbheRecHit.e[hfRecHit.n] = hit.energy();
258	<	hbheRecHit.et[hfRecHit.n] = getEt(hit.id(),hit.energy());
259	<	hbheRecHit.eta[hfRecHit.n] = getEta(hit.id());
260	<	hbheRecHit.phi[hfRecHit.n] = getPhi(hit.id());
261	<	//     hbheRecHit.ieta[hfRecHit.n] = hit.id().ieta();
262	<	//     hbheRecHit.iphi[hfRecHit.n] = hit.id().iphi();
263	<
257	>	hbheRecHit.e[hbheRecHit.n] = hit.energy();
258	>	hbheRecHit.et[hbheRecHit.n] = getEt(hit.id(),hit.energy());
259	>	hbheRecHit.eta[hbheRecHit.n] = getEta(hit.id());
260	>	hbheRecHit.phi[hbheRecHit.n] = getPhi(hit.id());
261	>	hbheRecHit.isjet[hbheRecHit.n] = false;
262	>	if(!useJets_){
263	>	for(unsigned int j = 0 ; j < jets->size(); ++j){
264	>	const reco::Jet& jet = (*jets)[j];
265	>	double dr = reco::deltaR(hbheRecHit.eta[hbheRecHit.n],hbheRecHit.phi[hbheRecHit.n],jet.eta(),jet.phi());
266	>	if(dr < cone){ hbheRecHit.isjet[hbheRecHit.n] = true; }
267	>	}
268	>	}
269		hbheRecHit.n++;
270		}
271	+
272	+	for(unsigned int i = 0; i < ebHits->size(); ++i){
273	+	const EcalRecHit & hit= (*ebHits)[i];
274	+	ebRecHit.e[ebRecHit.n] = hit.energy();
275	+	ebRecHit.et[ebRecHit.n] = getEt(hit.id(),hit.energy());
276	+	ebRecHit.eta[ebRecHit.n] = getEta(hit.id());
277	+	ebRecHit.phi[ebRecHit.n] = getPhi(hit.id());
278	+	ebRecHit.isjet[ebRecHit.n] = false;
279	+	if(!useJets_){
280	+	for(unsigned int j = 0 ; j < jets->size(); ++j){
281	+	const reco::Jet& jet = (*jets)[j];
282	+	double dr = reco::deltaR(ebRecHit.eta[ebRecHit.n],ebRecHit.phi[ebRecHit.n],jet.eta(),jet.phi());
283	+	if(dr < cone){ ebRecHit.isjet[ebRecHit.n] = true; }
284	+	}
285	+	}
286	+	ebRecHit.n++;
287	+	}
288	+
289	+	for(unsigned int i = 0; i < eeHits->size(); ++i){
290	+	const EcalRecHit & hit= (*eeHits)[i];
291	+	eeRecHit.e[eeRecHit.n] = hit.energy();
292	+	eeRecHit.et[eeRecHit.n] = getEt(hit.id(),hit.energy());
293	+	eeRecHit.eta[eeRecHit.n] = getEta(hit.id());
294	+	eeRecHit.phi[eeRecHit.n] = getPhi(hit.id());
295	+	eeRecHit.isjet[eeRecHit.n] = false;
296	+	if(!useJets_){
297	+	for(unsigned int j = 0 ; j < jets->size(); ++j){
298	+	const reco::Jet& jet = (*jets)[j];
299	+	double dr = reco::deltaR(eeRecHit.eta[eeRecHit.n],eeRecHit.phi[eeRecHit.n],jet.eta(),jet.phi());
300	+	if(dr < cone){ eeRecHit.isjet[eeRecHit.n] = true; }
301	+	}
302	+	}
303	+	eeRecHit.n++;
304	+	}
305	+
306
307	+	if(doBasicClusters_){
308	+	for(unsigned int j = 0 ; j < jets->size(); ++j){
309	+	const reco::Jet& jet = (*jets)[j];
310	+	myBC.n = 0;
311	+	myBC.jtpt = jet.pt();
312	+	myBC.jteta = jet.eta();
313	+	myBC.jtphi = jet.phi();
314	+
315	+	for(unsigned int i = 0; i < bClusters->size(); ++i){
316	+	const reco::BasicCluster & bc= (*bClusters)[i];
317	+	double dr = reco::deltaR(bc.eta(),bc.phi(),jet.eta(),jet.phi());
318	+	if(dr < cone){
319	+	myBC.e[myBC.n] = bc.energy();
320	+	myBC.et[myBC.n] = bc.energy()*sin(bc.position().theta());
321	+	myBC.eta[myBC.n] = bc.eta();
322	+	myBC.phi[myBC.n] = bc.phi();
323	+	myBC.n++;
324	+	}
325	+	}
326	+	bcTree->Fill();
327	+	}
328	+	}
329	+
330	+	eeTree->Fill();
331	+	ebTree->Fill();
332	+
333		hbheTree->Fill();
334		hfTree->Fill();
335	<
335	>
336		}
337
338
#	Line 222 | Line 340 | RecHitTreeProducer::analyze(const edm::E
340		void
341		RecHitTreeProducer::beginJob()
342		{
343	<
343	>
344		hbheTree = fs->make<TTree>("hbhe","");
345		hbheTree->Branch("n",&hbheRecHit.n,"n/I");
346		hbheTree->Branch("e",hbheRecHit.e,"e[n]/F");
347		hbheTree->Branch("et",hbheRecHit.et,"et[n]/F");
348		hbheTree->Branch("eta",hbheRecHit.eta,"eta[n]/F");
349		hbheTree->Branch("phi",hbheRecHit.phi,"phi[n]/F");
350	<
350	>	hbheTree->Branch("isjet",hbheRecHit.isjet,"isjet[n]/O");
351	>
352		hfTree = fs->make<TTree>("hf","");
353		hfTree->Branch("n",&hfRecHit.n,"n/I");
354		hfTree->Branch("e",hfRecHit.e,"e[n]/F");
355		hfTree->Branch("et",hfRecHit.et,"et[n]/F");
356		hfTree->Branch("eta",hfRecHit.eta,"eta[n]/F");
357		hfTree->Branch("phi",hfRecHit.phi,"phi[n]/F");
358	+	hfTree->Branch("isjet",hfRecHit.isjet,"isjet[n]/O");
359	+
360	+	eeTree = fs->make<TTree>("ee","");
361	+	eeTree->Branch("n",&eeRecHit.n,"n/I");
362	+	eeTree->Branch("e",eeRecHit.e,"e[n]/F");
363	+	eeTree->Branch("et",eeRecHit.et,"et[n]/F");
364	+	eeTree->Branch("eta",eeRecHit.eta,"eta[n]/F");
365	+	eeTree->Branch("phi",eeRecHit.phi,"phi[n]/F");
366	+	eeTree->Branch("isjet",eeRecHit.isjet,"isjet[n]/O");
367	+
368	+	ebTree = fs->make<TTree>("eb","");
369	+	ebTree->Branch("n",&ebRecHit.n,"n/I");
370	+	ebTree->Branch("e",ebRecHit.e,"e[n]/F");
371	+	ebTree->Branch("et",ebRecHit.et,"et[n]/F");
372	+	ebTree->Branch("eta",ebRecHit.eta,"eta[n]/F");
373	+	ebTree->Branch("phi",ebRecHit.phi,"phi[n]/F");
374	+	ebTree->Branch("isjet",ebRecHit.isjet,"isjet[n]/O");
375	+
376	+	if(doBasicClusters_){
377	+	bcTree = fs->make<TTree>("bc","");
378	+	bcTree->Branch("n",&myBC.n,"n/I");
379	+	bcTree->Branch("e",myBC.e,"e[n]/F");
380	+	bcTree->Branch("et",myBC.et,"et[n]/F");
381	+	bcTree->Branch("eta",myBC.eta,"eta[n]/F");
382	+	bcTree->Branch("phi",myBC.phi,"phi[n]/F");
383	+	bcTree->Branch("jtpt",&myBC.jtpt,"jtpt/F");
384	+	bcTree->Branch("jteta",&myBC.jteta,"jteta/F");
385	+	bcTree->Branch("jtphi",&myBC.jtphi,"jtphi/F");
386	+	//     bcTree->Branch("isjet",bcRecHit.isjet,"isjet[n]/O");
387	+	}
388	+
389	+
390
391		}
392
#	Line 249 | Line 400 | double RecHitTreeProducer::getEt(const D
400		double et = energy*sin(pos.theta());
401		return et;
402		}
403	<
403	>
404		double RecHitTreeProducer::getEta(const DetId &id){
405		const GlobalPoint& pos=geo->getPosition(id);
406		double et = pos.eta();

Diff Legend

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