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