ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/CmsHi/JetAnalysis/src/RecHitComparison.cc
Revision: 1.3
Committed: Wed Oct 20 14:11:46 2010 UTC (14 years, 6 months ago) by yilmaz
Content type: text/plain
Branch: MAIN
Changes since 1.2: +5 -3 lines
Log Message:
upadte

File Contents

# User Rev Content
1 yilmaz 1.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 yilmaz 1.3 // $Id: RecHitComparison.cc,v 1.2 2010/10/18 16:13:37 yilmaz Exp $
17 yilmaz 1.1 //
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/Centrality.h"
52     #include "SimDataFormats/HiGenData/interface/GenHIEvent.h"
53     #include "DataFormats/PatCandidates/interface/Jet.h"
54     #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
55 yilmaz 1.2 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
56     #include "DataFormats/HcalDetId/interface/HcalDetId.h"
57 yilmaz 1.1
58     #include "TNtuple.h"
59    
60     using namespace std;
61    
62     //
63     // class declaration
64     //
65    
66     class RecHitComparison : public edm::EDAnalyzer {
67     public:
68     explicit RecHitComparison(const edm::ParameterSet&);
69     ~RecHitComparison();
70    
71    
72     private:
73     virtual void beginJob() ;
74     virtual void analyze(const edm::Event&, const edm::EventSetup&);
75     virtual void endJob() ;
76    
77     // ----------member data ---------------------------
78    
79    
80     edm::Handle<reco::Centrality> cent;
81     edm::Handle<vector<double> > ktRhos;
82     edm::Handle<vector<double> > akRhos;
83    
84 yilmaz 1.2 edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > ebHits1;
85     edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > ebHits2;
86     edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > eeHits1;
87     edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > eeHits2;
88    
89     edm::Handle<HFRecHitCollection> hfHits1;
90     edm::Handle<HFRecHitCollection> hfHits2;
91     edm::Handle<HBHERecHitCollection> hbheHits1;
92     edm::Handle<HBHERecHitCollection> hbheHits2;
93 yilmaz 1.1
94     typedef vector<EcalRecHit>::const_iterator EcalIterator;
95 yilmaz 1.2 typedef vector<HFRecHit>::const_iterator HFIterator;
96     typedef vector<HBHERecHit>::const_iterator HBHEIterator;
97 yilmaz 1.1
98     edm::Handle<reco::CaloJetCollection> signalJets;
99    
100 yilmaz 1.2 edm::InputTag HcalRecHitHFSrc1_;
101     edm::InputTag HcalRecHitHFSrc2_;
102     edm::InputTag HcalRecHitHBHESrc1_;
103     edm::InputTag HcalRecHitHBHESrc2_;
104     edm::InputTag EBSrc1_;
105     edm::InputTag EBSrc2_;
106     edm::InputTag EESrc1_;
107     edm::InputTag EESrc2_;
108 yilmaz 1.3 edm::InputTag signalTag_;
109 yilmaz 1.2
110     TNtuple* ntEB;
111     TNtuple* ntEE;
112     TNtuple* ntHBHE;
113     TNtuple* ntHF;
114 yilmaz 1.1 TNtuple* ntjet;
115    
116     double cone;
117    
118     edm::Service<TFileService> fs;
119     const CentralityBins * cbins_;
120     const CaloGeometry *geo;
121     };
122    
123     //
124     // constants, enums and typedefs
125     //
126    
127     //
128     // static data member definitions
129     //
130    
131     //
132     // constructors and destructor
133     //
134     RecHitComparison::RecHitComparison(const edm::ParameterSet& iConfig) :
135     cbins_(0),
136     geo(0),
137     cone(0.5)
138     {
139     //now do what ever initialization is needed
140 yilmaz 1.3 signalTag_ = iConfig.getUntrackedParameter<edm::InputTag>("signalJets",edm::InputTag("iterativeCone5CaloJets","","SIGNAL"));
141    
142 yilmaz 1.2 HcalRecHitHFSrc1_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHFRecHitSrc1",edm::InputTag("hfreco"));
143     HcalRecHitHFSrc2_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHFRecHitSrc2",edm::InputTag("hfreco"));
144     HcalRecHitHBHESrc1_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHBHERecHitSrc1",edm::InputTag("hbhereco"));
145     HcalRecHitHBHESrc2_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHBHERecHitSrc2",edm::InputTag("hbhereco"));
146     EBSrc1_ = iConfig.getUntrackedParameter<edm::InputTag>("EBRecHitSrc1",edm::InputTag("ecalRecHit","EcalRecHitsEB","RECO"));
147     EBSrc2_ = iConfig.getUntrackedParameter<edm::InputTag>("EBRecHitSrc2",edm::InputTag("ecalRecHit","EcalRecHitsEB","SIGNALONLY"));
148     EESrc1_ = iConfig.getUntrackedParameter<edm::InputTag>("EERecHitSrc1",edm::InputTag("ecalRecHit","EcalRecHitsEE","RECO"));
149     EESrc2_ = iConfig.getUntrackedParameter<edm::InputTag>("EERecHitSrc2",edm::InputTag("ecalRecHit","EcalRecHitsEE","SIGNALONLY"));
150 yilmaz 1.1
151     }
152    
153    
154     RecHitComparison::~RecHitComparison()
155     {
156    
157     // do anything here that needs to be done at desctruction time
158     // (e.g. close files, deallocate resources etc.)
159    
160     }
161    
162    
163     //
164     // member functions
165     //
166    
167     // ------------ method called to for each event ------------
168     void
169     RecHitComparison::analyze(const edm::Event& ev, const edm::EventSetup& iSetup)
170     {
171    
172     if(!cbins_) cbins_ = getCentralityBinsFromDB(iSetup);
173     if(!geo){
174     edm::ESHandle<CaloGeometry> pGeo;
175     iSetup.get<CaloGeometryRecord>().get(pGeo);
176     geo = pGeo.product();
177     }
178    
179     edm::InputTag centTag("hiCentrality","","RECO");
180    
181     using namespace edm;
182    
183     ev.getByLabel(centTag,cent);
184 yilmaz 1.2 ev.getByLabel(EBSrc1_,ebHits1);
185     ev.getByLabel(EBSrc1_,ebHits2);
186 yilmaz 1.3 ev.getByLabel(signalTag_,signalJets);
187 yilmaz 1.1
188 yilmaz 1.2 ev.getByLabel(HcalRecHitHFSrc1_,hfHits1);
189     ev.getByLabel(HcalRecHitHFSrc2_,hfHits2);
190     ev.getByLabel(HcalRecHitHBHESrc1_,hbheHits1);
191     ev.getByLabel(HcalRecHitHBHESrc2_,hbheHits2);
192    
193     ev.getByLabel(EESrc1_,eeHits1);
194     ev.getByLabel(EESrc2_,eeHits2);
195    
196 yilmaz 1.1 double hf = cent->EtHFhitSum();
197     double sumet = cent->EtMidRapiditySum();
198     int run = ev.id().run();
199     run = 1;
200     int bin = cbins_->getBin(hf);
201     int margin = 0;
202    
203     vector<double> fFull;
204     vector<double> f05;
205     vector<double> f1;
206     vector<double> f15;
207     vector<double> f2;
208     vector<double> f25;
209     vector<double> f3;
210    
211     int njets = signalJets->size();
212     fFull.reserve(njets);
213     f05.reserve(njets);
214     f1.reserve(njets);
215     f15.reserve(njets);
216     f2.reserve(njets);
217     f25.reserve(njets);
218     f3.reserve(njets);
219    
220     for(unsigned int j1 = 0 ; j1 < signalJets->size(); ++j1){
221     fFull.push_back(0);
222     f05.push_back(0);
223     f1.push_back(0);
224     f15.push_back(0);
225     f2.push_back(0);
226     f25.push_back(0);
227     f3.push_back(0);
228     }
229    
230 yilmaz 1.2 for(unsigned int j1 = 0 ; j1 < ebHits1->size(); ++j1){
231 yilmaz 1.1
232 yilmaz 1.2 const EcalRecHit& jet1 = (*ebHits1)[j1];
233 yilmaz 1.1 double e1 = jet1.energy();
234    
235     const GlobalPoint& pos1=geo->getPosition(jet1.id());
236     double eta1 = pos1.eta();
237     double phi1 = pos1.phi();
238     double et1 = e1*sin(pos1.theta());
239    
240     double drjet = -1;
241     double jetpt = -1;
242     bool isjet = false;
243     int matchedJet = -1;
244    
245     for(unsigned int j = 0 ; j < signalJets->size(); ++j){
246     const reco::CaloJet & jet = (*signalJets)[j];
247     double dr = reco::deltaR(eta1,phi1,jet.eta(),jet.phi());
248     if(dr < cone){
249     jetpt = jet.pt();
250     drjet = dr;
251     isjet = true;
252     matchedJet = j;
253     fFull[j] += et1;
254    
255     if(et1 > 0.5){
256     f05[j] += et1;
257     }
258     if(et1 > 1.){
259     f1[j] += et1;
260     }
261     if(et1 > 1.5){
262     f15[j] += et1;
263     }
264     if(et1 > 2){
265     f2[j] += et1;
266     }
267     if(et1 > 2.5){
268     f25[j] += et1;
269     }
270     if(et1 > 3.){
271     f3[j] += et1;
272     }
273     }
274     }
275    
276     GlobalPoint pos2;
277     double e2 = -1;
278 yilmaz 1.2 EcalIterator hitit = ebHits2->find(jet1.id());
279     if(hitit != ebHits2->end()){
280 yilmaz 1.1 e2 = hitit->energy();
281     pos2=geo->getPosition(hitit->id());
282     }else{
283     e2 = 0;
284     pos2 = pos1;
285     }
286    
287     double eta2 = pos2.eta();
288     double phi2 = pos2.eta();
289     double et2 = e2*sin(pos2.theta());
290 yilmaz 1.2 if(isjet) ntEB->Fill(e1,et1,e2,et2,eta2,phi2,sumet,hf,bin,jetpt,drjet);
291     }
292    
293     for(unsigned int i = 0; i < eeHits1->size(); ++i){
294     const EcalRecHit & jet1= (*eeHits1)[i];
295     double e1 = jet1.energy();
296     const GlobalPoint& pos1=geo->getPosition(jet1.id());
297     double eta1 = pos1.eta();
298     double phi1 = pos1.phi();
299     double et1 = e1*sin(pos1.theta());
300     double drjet = -1;
301     double jetpt = -1;
302     bool isjet = false;
303     int matchedJet = -1;
304     for(unsigned int j = 0 ; j < signalJets->size(); ++j){
305     const reco::CaloJet & jet = (*signalJets)[j];
306     double dr = reco::deltaR(eta1,phi1,jet.eta(),jet.phi());
307     if(dr < cone){
308     jetpt = jet.pt();
309     drjet = dr;
310     isjet = true;
311     matchedJet = j;
312     }
313     }
314     GlobalPoint pos2;
315     double e2 = -1;
316     EcalIterator hitit = eeHits2->find(jet1.id());
317     if(hitit != eeHits2->end()){
318     e2 = hitit->energy();
319     pos2=geo->getPosition(hitit->id());
320     }else{
321     e2 = 0;
322     pos2 = pos1;
323     }
324     double eta2 = pos2.eta();
325     double phi2 = pos2.eta();
326     double et2 = e2*sin(pos2.theta());
327     if(isjet) ntEE->Fill(e1,et1,e2,et2,eta2,phi2,sumet,hf,bin,jetpt,drjet);
328     }
329    
330     for(unsigned int i = 0; i < hbheHits1->size(); ++i){
331     const HBHERecHit & jet1= (*hbheHits1)[i];
332     double e1 = jet1.energy();
333     const GlobalPoint& pos1=geo->getPosition(jet1.id());
334     double eta1 = pos1.eta();
335     double phi1 = pos1.phi();
336     double et1 = e1*sin(pos1.theta());
337     double drjet = -1;
338     double jetpt = -1;
339     bool isjet = false;
340     int matchedJet = -1;
341     for(unsigned int j = 0 ; j < signalJets->size(); ++j){
342     const reco::CaloJet & jet = (*signalJets)[j];
343     double dr = reco::deltaR(eta1,phi1,jet.eta(),jet.phi());
344     if(dr < cone){
345     jetpt = jet.pt();
346     drjet = dr;
347     isjet = true;
348     matchedJet = j;
349     }
350     }
351     GlobalPoint pos2;
352     double e2 = -1;
353     HBHEIterator hitit = hbheHits2->find(jet1.id());
354     if(hitit != hbheHits2->end()){
355     e2 = hitit->energy();
356     pos2=geo->getPosition(hitit->id());
357     }else{
358     e2 = 0;
359     pos2 = pos1;
360     }
361     double eta2 = pos2.eta();
362     double phi2 = pos2.eta();
363     double et2 = e2*sin(pos2.theta());
364     if(isjet) ntHBHE->Fill(e1,et1,e2,et2,eta2,phi2,sumet,hf,bin,jetpt,drjet);
365     }
366    
367     for(unsigned int i = 0; i < hfHits1->size(); ++i){
368     const HFRecHit & jet1= (*hfHits1)[i];
369     double e1 = jet1.energy();
370     const GlobalPoint& pos1=geo->getPosition(jet1.id());
371     double eta1 = pos1.eta();
372     double phi1 = pos1.phi();
373     double et1 = e1*sin(pos1.theta());
374     double drjet = -1;
375     double jetpt = -1;
376     bool isjet = false;
377     int matchedJet = -1;
378     for(unsigned int j = 0 ; j < signalJets->size(); ++j){
379     const reco::CaloJet & jet = (*signalJets)[j];
380     double dr = reco::deltaR(eta1,phi1,jet.eta(),jet.phi());
381     if(dr < cone){
382     jetpt = jet.pt();
383     drjet = dr;
384     isjet = true;
385     matchedJet = j;
386     }
387     }
388     GlobalPoint pos2;
389     double e2 = -1;
390     HFIterator hitit = hfHits2->find(jet1.id());
391     if(hitit != hfHits2->end()){
392     e2 = hitit->energy();
393     pos2=geo->getPosition(hitit->id());
394     }else{
395     e2 = 0;
396     pos2 = pos1;
397     }
398     double eta2 = pos2.eta();
399     double phi2 = pos2.eta();
400     double et2 = e2*sin(pos2.theta());
401     if(isjet) ntHF->Fill(e1,et1,e2,et2,eta2,phi2,sumet,hf,bin,jetpt,drjet);
402 yilmaz 1.1 }
403    
404     for(unsigned int j1 = 0 ; j1 < signalJets->size(); ++j1){
405     const reco::CaloJet & jet = (*signalJets)[j1];
406     double em = (jet.emEnergyInEB() + jet.emEnergyInEE()) * sin(jet.theta());
407     double emf = jet.emEnergyFraction();
408     double pt = jet.pt();
409     double eta = jet.eta();
410     ntjet->Fill(bin,pt,eta,fFull[j1],f05[j1],f1[j1],f15[j1],f2[j1],f25[j1],f3[j1],em,emf);
411     }
412    
413     }
414    
415    
416     // ------------ method called once each job just before starting event loop ------------
417     void
418     RecHitComparison::beginJob()
419     {
420 yilmaz 1.2 ntEB = fs->make<TNtuple>("ntEB","","e1:et1:e2:et2:eta:phi:sumet:hf:bin:ptjet:drjet");
421     ntEE = fs->make<TNtuple>("ntEE","","e1:et1:e2:et2:eta:phi:sumet:hf:bin:ptjet:drjet");
422     ntHBHE = fs->make<TNtuple>("ntHBHE","","e1:et1:e2:et2:eta:phi:sumet:hf:bin:ptjet:drjet");
423     ntHF = fs->make<TNtuple>("ntHF","","e1:et1:e2:et2:eta:phi:sumet:hf:bin:ptjet:drjet");
424 yilmaz 1.1 ntjet = fs->make<TNtuple>("ntjet","","bin:pt:eta:ethit:f05:f1:f15:f2:f25:f3:em:emf");
425    
426     }
427    
428     // ------------ method called once each job just after ending the event loop ------------
429     void
430     RecHitComparison::endJob() {
431     }
432    
433     //define this as a plug-in
434     DEFINE_FWK_MODULE(RecHitComparison);