ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/CmsHi/JetAnalysis/src/RecHitComparison.cc
Revision: 1.2
Committed: Mon Oct 18 16:13:37 2010 UTC (14 years, 6 months ago) by yilmaz
Content type: text/plain
Branch: MAIN
Changes since 1.1: +165 -20 lines
Log Message:
compilable analyzers

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