ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/CmsHi/JetAnalysis/src/RecHitComparison.cc
Revision: 1.10
Committed: Wed Mar 30 21:21:10 2011 UTC (14 years, 1 month ago) by yilmaz
Content type: text/plain
Branch: MAIN
CVS Tags: HiForest_V02_85, HiForest_V02_84, HiForest_V02_83, HiForest_V02_82, HiForest_V02_81, HiForest_V02_80, HiForest_V02_79, HiForest_V02_78, HiForest_V02_77, HiForest_V02_76, HiForest_V02_75, HiForest_V02_74, HiForest_V02_73, HiForest_V02_72, HiForest_V02_71, HiForest_V02_70, HiForest_V02_69, HiForest_V02_68, HiForest_V02_67, HiForest_V02_66, HiForest_V02_65, HiForest_V02_64, HiForest_V02_63, HiForest_V02_62, HiForest_V02_61, HiForest_V02_60, HiForest_V02_59, HiForest_V02_58, HiForest_V02_57, HiForest_V02_56, HiForest_V02_55, HiForest_V02_54, HiForest_V02_53, HiForest_V02_52, HiForest_V02_51, HiForest_V02_50, HiForest_V02_49, HiForest_V02_48, HiForest_V02_47, HiForest_V02_46, HiForest_V02_45, HiForest_V02_44, HiForest_V02_43, HiForest_V02_42, HiForest_V02_41, HiForest_V02_40, HiForest_V02_39, HiForest_V02_38, HiForest_V02_37, HiForest_V02_36, HiForest_V02_35, HiForest_V02_34, HiForest_V02_33, HiForest_V02_32, HiForest_V02_31, HiForest_V02_30, HiForest_V02_27, HiForest_V02_26, QM_2012, HiForest_V02_25, HiForest_V02_24, HiForest_V02_23, HiForest_V02_22, HiForest_V02_21, HiForest_V02_20, HiForest_V02_19, HiForest_V02_18, HiForest_V02_17, HiForest_V02_16, HiForest_V02_15, HiForest_V02_14, HiForest_V02_13, HiForest_V02_12, HiForest_V02_11, HiForest_V02_10, HiForest_V02_09, HiForest_V02_08, HiForest_V02_07, HiForest_V02_06, HiForest_V02_05, HiForest_V02_04, HiForest_V02_03, HiForest_V02_02, HiForest_V02_01, HiForest_V02_00, hi44X_02, hi413_03, hi441_1, hi441_0, hi413_11, hi413_10, hi413_09, hi413_08, hi413_07, hi413_06, hi413_05, hi413_04, hi413_02, hi39X_01, tag_d20110915, cmssw39x_base, cmssw39X_base, HEAD
Branch point for: branch_44x, cmssw39x_branch
Changes since 1.9: +26 -20 lines
Log Message:
jet optional

File Contents

# Content
1 // -*- C++ -*-
2 //
3 // Package: RecHitComparison
4 // Class: RecHitComparison
5 //
6 /**\class RecHitComparison RecHitComparison.cc CmsHi/RecHitComparison/src/RecHitComparison.cc
7
8 Description: [one line class summary]
9
10 Implementation:
11 [Notes on implementation]
12 */
13 //
14 // Original Author: Yetkin Yilmaz
15 // Created: Tue Sep 7 11:38:19 EDT 2010
16 // $Id: RecHitComparison.cc,v 1.9 2011/03/30 12:13:22 yilmaz Exp $
17 //
18 //
19
20
21 // system include files
22 #include <memory>
23 #include <vector>
24 #include <iostream>
25
26 // user include files
27 #include "FWCore/Framework/interface/Frameworkfwd.h"
28 #include "FWCore/Framework/interface/EDAnalyzer.h"
29
30 #include "FWCore/Framework/interface/Event.h"
31 #include "FWCore/Framework/interface/MakerMacros.h"
32 #include "FWCore/Framework/interface/ESHandle.h"
33 #include "FWCore/ParameterSet/interface/ParameterSet.h"
34
35 #include "DataFormats/DetId/interface/DetId.h"
36 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
37 #include "Geometry/Records/interface/IdealGeometryRecord.h"
38 #include "Geometry/Records/interface/CaloGeometryRecord.h"
39 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
40
41 #include "CommonTools/UtilAlgos/interface/TFileService.h"
42 #include "FWCore/ServiceRegistry/interface/Service.h"
43
44 #include "DataFormats/Math/interface/deltaR.h"
45
46 #include "DataFormats/Common/interface/Handle.h"
47 #include "DataFormats/FWLite/interface/ChainEvent.h"
48 #include "FWCore/Utilities/interface/InputTag.h"
49 #include "DataFormats/HeavyIonEvent/interface/CentralityBins.h"
50 #include "DataFormats/CaloTowers/interface/CaloTower.h"
51 #include "DataFormats/HeavyIonEvent/interface/CentralityProvider.h"
52 #include "SimDataFormats/HiGenData/interface/GenHIEvent.h"
53 #include "DataFormats/PatCandidates/interface/Jet.h"
54 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
55 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
56 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
57
58 #include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
59
60 #include "TNtuple.h"
61
62 using namespace std;
63
64 //
65 // class declaration
66 //
67
68 class RecHitComparison : public edm::EDAnalyzer {
69 public:
70 explicit RecHitComparison(const edm::ParameterSet&);
71 ~RecHitComparison();
72
73
74 private:
75 virtual void beginJob() ;
76 virtual void analyze(const edm::Event&, const edm::EventSetup&);
77 virtual void endJob() ;
78
79 // ----------member data ---------------------------
80 edm::Handle<vector<double> > ktRhos;
81 edm::Handle<vector<double> > akRhos;
82
83 edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > ebHits1;
84 edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > ebHits2;
85 edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > eeHits1;
86 edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > eeHits2;
87
88 edm::Handle<HFRecHitCollection> hfHits1;
89 edm::Handle<HFRecHitCollection> hfHits2;
90 edm::Handle<HBHERecHitCollection> hbheHits1;
91 edm::Handle<HBHERecHitCollection> hbheHits2;
92
93 edm::Handle<reco::BasicClusterCollection> bClusters1;
94 edm::Handle<reco::BasicClusterCollection> bClusters2;
95
96
97 typedef vector<EcalRecHit>::const_iterator EcalIterator;
98 typedef vector<HFRecHit>::const_iterator HFIterator;
99 typedef vector<HBHERecHit>::const_iterator HBHEIterator;
100
101 edm::Handle<reco::CaloJetCollection> signalJets;
102
103 edm::InputTag HcalRecHitHFSrc1_;
104 edm::InputTag HcalRecHitHFSrc2_;
105 edm::InputTag HcalRecHitHBHESrc1_;
106 edm::InputTag HcalRecHitHBHESrc2_;
107 edm::InputTag EBSrc1_;
108 edm::InputTag EBSrc2_;
109 edm::InputTag EESrc1_;
110 edm::InputTag EESrc2_;
111
112 edm::InputTag BCSrc1_;
113 edm::InputTag BCSrc2_;
114
115 edm::InputTag signalTag_;
116
117 TNtuple* ntBC;
118 TNtuple* ntEB;
119 TNtuple* ntEE;
120 TNtuple* ntHBHE;
121 TNtuple* ntHF;
122 TNtuple* ntjet;
123
124 double cone;
125 bool jetsOnly_;
126 bool doBasicClusters_;
127 bool doJetCone_;
128
129 edm::Service<TFileService> fs;
130 CentralityProvider* centrality_;
131 const CaloGeometry *geo;
132 };
133
134 //
135 // constants, enums and typedefs
136 //
137
138 //
139 // static data member definitions
140 //
141
142 //
143 // constructors and destructor
144 //
145 RecHitComparison::RecHitComparison(const edm::ParameterSet& iConfig) :
146 centrality_(0),
147 geo(0),
148 cone(0.5)
149 {
150 //now do what ever initialization is needed
151 jetsOnly_ = iConfig.getUntrackedParameter<bool>("jetsOnly",false);
152 doBasicClusters_ = iConfig.getUntrackedParameter<bool>("doBasicClusters",false);
153
154 doJetCone_ = iConfig.getUntrackedParameter<bool>("doJetCone",false);
155 signalTag_ = iConfig.getUntrackedParameter<edm::InputTag>("signalJets",edm::InputTag("iterativeCone5CaloJets","","SIGNAL"));
156
157 if(!doJetCone_) jetsOnly_ = 0;
158
159 HcalRecHitHFSrc1_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHFRecHitSrc1",edm::InputTag("hfreco"));
160 HcalRecHitHFSrc2_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHFRecHitSrc2",edm::InputTag("hfreco"));
161 HcalRecHitHBHESrc1_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHBHERecHitSrc1",edm::InputTag("hbhereco"));
162 HcalRecHitHBHESrc2_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHBHERecHitSrc2",edm::InputTag("hbhereco"));
163 EBSrc1_ = iConfig.getUntrackedParameter<edm::InputTag>("EBRecHitSrc1",edm::InputTag("ecalRecHit","EcalRecHitsEB","RECOBKG"));
164 EBSrc2_ = iConfig.getUntrackedParameter<edm::InputTag>("EBRecHitSrc2",edm::InputTag("ecalRecHit","EcalRecHitsEB","S"));
165 EESrc1_ = iConfig.getUntrackedParameter<edm::InputTag>("EERecHitSrc1",edm::InputTag("ecalRecHit","EcalRecHitsEE","RECO"));
166 EESrc2_ = iConfig.getUntrackedParameter<edm::InputTag>("EERecHitSrc2",edm::InputTag("ecalRecHit","EcalRecHitsEE","SIGNALONLY"));
167 BCSrc1_ = iConfig.getUntrackedParameter<edm::InputTag>("BasicClusterSrc1",edm::InputTag("ecalRecHit","EcalRecHitsEB","RECO"));
168 BCSrc2_ = iConfig.getUntrackedParameter<edm::InputTag>("BasicClusterSrc2",edm::InputTag("ecalRecHit","EcalRecHitsEB","SIGNALONLY"));
169
170 }
171
172
173 RecHitComparison::~RecHitComparison()
174 {
175
176 // do anything here that needs to be done at desctruction time
177 // (e.g. close files, deallocate resources etc.)
178
179 }
180
181
182 //
183 // member functions
184 //
185
186 // ------------ method called to for each event ------------
187 void
188 RecHitComparison::analyze(const edm::Event& ev, const edm::EventSetup& iSetup)
189 {
190 if(!centrality_) centrality_ = new CentralityProvider(iSetup);
191 if(!geo){
192 edm::ESHandle<CaloGeometry> pGeo;
193 iSetup.get<CaloGeometryRecord>().get(pGeo);
194 geo = pGeo.product();
195 }
196
197 using namespace edm;
198 ev.getByLabel(EBSrc1_,ebHits1);
199 ev.getByLabel(EBSrc2_,ebHits2);
200
201 if(doJetCone_) ev.getByLabel(signalTag_,signalJets);
202
203 ev.getByLabel(HcalRecHitHFSrc1_,hfHits1);
204 ev.getByLabel(HcalRecHitHFSrc2_,hfHits2);
205 ev.getByLabel(HcalRecHitHBHESrc1_,hbheHits1);
206 ev.getByLabel(HcalRecHitHBHESrc2_,hbheHits2);
207 ev.getByLabel(EESrc1_,eeHits1);
208 ev.getByLabel(EESrc2_,eeHits2);
209
210 if(doBasicClusters_){
211 ev.getByLabel(BCSrc1_,bClusters1);
212 ev.getByLabel(BCSrc2_,bClusters2);
213 }
214
215 centrality_->newEvent(ev,iSetup);
216
217 double hf = centrality_->centralityValue();
218 int bin = centrality_->getBin();
219
220 vector<double> fFull;
221 vector<double> f05;
222 vector<double> f1;
223 vector<double> f15;
224 vector<double> f2;
225 vector<double> f25;
226 vector<double> f3;
227
228 int njets = 0;
229
230 if(doJetCone_) njets = signalJets->size();
231 fFull.reserve(njets);
232 f05.reserve(njets);
233 f1.reserve(njets);
234 f15.reserve(njets);
235 f2.reserve(njets);
236 f25.reserve(njets);
237 f3.reserve(njets);
238
239 if(doJetCone_){
240 for(unsigned int j1 = 0 ; j1 < signalJets->size(); ++j1){
241 fFull.push_back(0);
242 f05.push_back(0);
243 f1.push_back(0);
244 f15.push_back(0);
245 f2.push_back(0);
246 f25.push_back(0);
247 f3.push_back(0);
248 }
249 }
250
251 for(unsigned int j1 = 0 ; j1 < ebHits1->size(); ++j1){
252
253 const EcalRecHit& jet1 = (*ebHits1)[j1];
254 double e1 = jet1.energy();
255
256 const GlobalPoint& pos1=geo->getPosition(jet1.id());
257 double eta1 = pos1.eta();
258 double phi1 = pos1.phi();
259 double et1 = e1*sin(pos1.theta());
260
261 double drjet = -1;
262 double jetpt = -1;
263 bool isjet = false;
264 int matchedJet = -1;
265
266 if(doJetCone_){
267 for(unsigned int j = 0 ; j < signalJets->size(); ++j){
268 const reco::CaloJet & jet = (*signalJets)[j];
269 double dr = reco::deltaR(eta1,phi1,jet.eta(),jet.phi());
270 if(dr < cone){
271 jetpt = jet.pt();
272 drjet = dr;
273 isjet = true;
274 matchedJet = j;
275 fFull[j] += et1;
276
277 if(et1 > 0.5){
278 f05[j] += et1;
279 }
280 if(et1 > 1.){
281 f1[j] += et1;
282 }
283 if(et1 > 1.5){
284 f15[j] += et1;
285 }
286 if(et1 > 2){
287 f2[j] += et1;
288 }
289 if(et1 > 2.5){
290 f25[j] += et1;
291 }
292 if(et1 > 3.){
293 f3[j] += et1;
294 }
295 }
296 }
297 }
298
299 GlobalPoint pos2;
300 double e2 = -1;
301 EcalIterator hitit = ebHits2->find(jet1.id());
302 if(hitit != ebHits2->end()){
303 e2 = hitit->energy();
304 pos2=geo->getPosition(hitit->id());
305 }else{
306 e2 = 0;
307 pos2 = pos1;
308 }
309
310 double eta2 = pos2.eta();
311 double phi2 = pos2.eta();
312 double et2 = e2*sin(pos2.theta());
313 if(!jetsOnly_ || isjet) ntEB->Fill(e1,et1,e2,et2,eta2,phi2,hf,bin,jetpt,drjet);
314 }
315
316 for(unsigned int i = 0; i < eeHits1->size(); ++i){
317 const EcalRecHit & jet1= (*eeHits1)[i];
318 double e1 = jet1.energy();
319 const GlobalPoint& pos1=geo->getPosition(jet1.id());
320 double eta1 = pos1.eta();
321 double phi1 = pos1.phi();
322 double et1 = e1*sin(pos1.theta());
323 double drjet = -1;
324 double jetpt = -1;
325 bool isjet = false;
326 int matchedJet = -1;
327 if(doJetCone_){
328 for(unsigned int j = 0 ; j < signalJets->size(); ++j){
329 const reco::CaloJet & jet = (*signalJets)[j];
330 double dr = reco::deltaR(eta1,phi1,jet.eta(),jet.phi());
331 if(dr < cone){
332 jetpt = jet.pt();
333 drjet = dr;
334 isjet = true;
335 matchedJet = j;
336 }
337 }
338 }
339
340 GlobalPoint pos2;
341 double e2 = -1;
342 EcalIterator hitit = eeHits2->find(jet1.id());
343 if(hitit != eeHits2->end()){
344 e2 = hitit->energy();
345 pos2=geo->getPosition(hitit->id());
346 }else{
347 e2 = 0;
348 pos2 = pos1;
349 }
350 double eta2 = pos2.eta();
351 double phi2 = pos2.eta();
352 double et2 = e2*sin(pos2.theta());
353 if(!jetsOnly_ || isjet) ntEE->Fill(e1,et1,e2,et2,eta2,phi2,hf,bin,jetpt,drjet);
354 }
355
356 for(unsigned int i = 0; i < hbheHits1->size(); ++i){
357 const HBHERecHit & jet1= (*hbheHits1)[i];
358 double e1 = jet1.energy();
359 const GlobalPoint& pos1=geo->getPosition(jet1.id());
360 double eta1 = pos1.eta();
361 double phi1 = pos1.phi();
362 double et1 = e1*sin(pos1.theta());
363 double drjet = -1;
364 double jetpt = -1;
365 bool isjet = false;
366 int matchedJet = -1;
367 if(doJetCone_){
368 for(unsigned int j = 0 ; j < signalJets->size(); ++j){
369 const reco::CaloJet & jet = (*signalJets)[j];
370 double dr = reco::deltaR(eta1,phi1,jet.eta(),jet.phi());
371 if(dr < cone){
372 jetpt = jet.pt();
373 drjet = dr;
374 isjet = true;
375 matchedJet = j;
376 }
377 }
378 }
379
380 GlobalPoint pos2;
381 double e2 = -1;
382 HBHEIterator hitit = hbheHits2->find(jet1.id());
383 if(hitit != hbheHits2->end()){
384 e2 = hitit->energy();
385 pos2=geo->getPosition(hitit->id());
386 }else{
387 e2 = 0;
388 pos2 = pos1;
389 }
390 double eta2 = pos2.eta();
391 double phi2 = pos2.eta();
392 double et2 = e2*sin(pos2.theta());
393 if(!jetsOnly_ || isjet) ntHBHE->Fill(e1,et1,e2,et2,eta2,phi2,hf,bin,jetpt,drjet);
394 }
395
396 for(unsigned int i = 0; i < hfHits1->size(); ++i){
397 const HFRecHit & jet1= (*hfHits1)[i];
398 double e1 = jet1.energy();
399 const GlobalPoint& pos1=geo->getPosition(jet1.id());
400 double eta1 = pos1.eta();
401 double phi1 = pos1.phi();
402 double et1 = e1*sin(pos1.theta());
403 double drjet = -1;
404 double jetpt = -1;
405 bool isjet = false;
406 int matchedJet = -1;
407 if(doJetCone_){
408 for(unsigned int j = 0 ; j < signalJets->size(); ++j){
409 const reco::CaloJet & jet = (*signalJets)[j];
410 double dr = reco::deltaR(eta1,phi1,jet.eta(),jet.phi());
411 if(dr < cone){
412 jetpt = jet.pt();
413 drjet = dr;
414 isjet = true;
415 matchedJet = j;
416 }
417 }
418 }
419 GlobalPoint pos2;
420 double e2 = -1;
421 HFIterator hitit = hfHits2->find(jet1.id());
422 if(hitit != hfHits2->end()){
423 e2 = hitit->energy();
424 pos2=geo->getPosition(hitit->id());
425 }else{
426 e2 = 0;
427 pos2 = pos1;
428 }
429 double eta2 = pos2.eta();
430 double phi2 = pos2.eta();
431 double et2 = e2*sin(pos2.theta());
432 if(!jetsOnly_ || isjet) ntHF->Fill(e1,et1,e2,et2,eta2,phi2,hf,bin,jetpt,drjet);
433 }
434
435 if(doJetCone_){
436 for(unsigned int j1 = 0 ; j1 < signalJets->size(); ++j1){
437 const reco::CaloJet & jet = (*signalJets)[j1];
438 double em = (jet.emEnergyInEB() + jet.emEnergyInEE()) * sin(jet.theta());
439 double emf = jet.emEnergyFraction();
440 double pt = jet.pt();
441 double eta = jet.eta();
442 ntjet->Fill(bin,pt,eta,fFull[j1],f05[j1],f1[j1],f15[j1],f2[j1],f25[j1],f3[j1],em,emf);
443 }
444 }
445
446 }
447
448
449 // ------------ method called once each job just before starting event loop ------------
450 void
451 RecHitComparison::beginJob()
452 {
453 ntEB = fs->make<TNtuple>("ntEB","","e1:et1:e2:et2:eta:phi:hf:bin:ptjet:drjet");
454 ntEE = fs->make<TNtuple>("ntEE","","e1:et1:e2:et2:eta:phi:hf:bin:ptjet:drjet");
455 ntHBHE = fs->make<TNtuple>("ntHBHE","","e1:et1:e2:et2:eta:phi:hf:bin:ptjet:drjet");
456 ntHF = fs->make<TNtuple>("ntHF","","e1:et1:e2:et2:eta:phi:hf:bin:ptjet:drjet");
457
458 ntBC = fs->make<TNtuple>("ntBC","","e1:et1:e2:et2:eta:phi:hf:bin:ptjet:drjet");
459
460 ntjet = fs->make<TNtuple>("ntjet","","bin:pt:eta:ethit:f05:f1:f15:f2:f25:f3:em:emf");
461
462 }
463
464 // ------------ method called once each job just after ending the event loop ------------
465 void
466 RecHitComparison::endJob() {
467 }
468
469 //define this as a plug-in
470 DEFINE_FWK_MODULE(RecHitComparison);