ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/CmsHi/JetAnalysis/src/HiJetResponseAnalyzer.cc
Revision: 1.12
Committed: Tue Sep 20 19:48:45 2011 UTC (13 years, 7 months ago) by yilmaz
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 1.11: +1 -1 lines
State: FILE REMOVED
Error occurred while calculating annotation data.
Log Message:
rename the matching analyzer

File Contents

# Content
1 // -*- C++ -*-
2 //
3 // Package: HiJetResponseAnalyzer
4 // Class: HiJetResponseAnalyzer
5 //
6 /**\class HiJetResponseAnalyzer HiJetResponseAnalyzer.cc CmsHi/HiJetResponseAnalyzer/src/HiJetResponseAnalyzer.cc
7
8 Description: [one line class summary]
9
10 Implementation:
11 [Notes on implementation]
12 */
13 //
14 // Original Author: Yetkin Yilmaz
15 // Created: Thu Sep 9 10:38:59 EDT 2010
16 // $Id: HiJetResponseAnalyzer.cc,v 1.11 2011/04/01 12:09:17 yilmaz Exp $
17 //
18 //
19
20
21 // system include files
22 #include <memory>
23 #include <iostream>
24
25 // user include files
26 #include "FWCore/Framework/interface/Frameworkfwd.h"
27 #include "FWCore/Framework/interface/EDAnalyzer.h"
28
29 #include "FWCore/Framework/interface/Event.h"
30 #include "FWCore/Framework/interface/MakerMacros.h"
31
32 #include "FWCore/ParameterSet/interface/ParameterSet.h"
33
34 #include "FWCore/Utilities/interface/InputTag.h"
35 #include "DataFormats/HeavyIonEvent/interface/CentralityBins.h"
36 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
37 #include "DataFormats/HeavyIonEvent/interface/Centrality.h"
38 #include "DataFormats/PatCandidates/interface/Jet.h"
39 #include "SimDataFormats/HiGenData/interface/GenHIEvent.h"
40 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
41 #include "DataFormats/JetReco/interface/GenJetCollection.h"
42
43 #include "CondFormats/JetMETObjects/interface/FactorizedJetCorrector.h"
44 #include "JetMETCorrections/Objects/interface/JetCorrectionsRecord.h"
45 #include "CondFormats/JetMETObjects/interface/JetCorrectionUncertainty.h"
46
47 #include "CommonTools/UtilAlgos/interface/TFileService.h"
48 #include "FWCore/ServiceRegistry/interface/Service.h"
49 #include "DataFormats/Math/interface/deltaR.h"
50 #include "FWCore/ParameterSet/interface/FileInPath.h"
51 #include "FWCore/Framework/interface/ESWatcher.h"
52
53 #include "CmsHi/JetAnalysis/interface/RhoGetter.h"
54 #include "TTree.h"
55
56 using namespace std;
57 using namespace edm;
58
59 static const int MAXJETS = 500;
60
61 struct etdr{
62 double et;
63 double dr;
64 };
65
66 class JRA{
67
68 public:
69 int nref;
70 int bin;
71 float b;
72 float hf;
73 float jtpt[MAXJETS];
74 float jtrawpt[MAXJETS];
75 float refpt[MAXJETS];
76 float jteta[MAXJETS];
77 float refeta[MAXJETS];
78 float jtphi[MAXJETS];
79 float refphi[MAXJETS];
80 float l2[MAXJETS];
81 float l3[MAXJETS];
82 float area[MAXJETS];
83 float pu[MAXJETS];
84 float rho[MAXJETS];
85
86 float weight;
87 };
88
89 struct JRAV{
90 int index;
91 float jtpt;
92 float jtrawpt;
93 float refpt;
94 float refcorpt;
95 float jteta;
96 float refeta;
97 float jtphi;
98 float refphi;
99 float l2;
100 float l3;
101 float area;
102 float pu;
103 float rho;
104
105 };
106
107 //
108 // class declaration
109 //
110 bool comparePt(JRAV a, JRAV b) {return a.jtpt > b.jtpt;}
111
112 class HiJetResponseAnalyzer : public edm::EDAnalyzer {
113 public:
114 explicit HiJetResponseAnalyzer(const edm::ParameterSet&);
115 ~HiJetResponseAnalyzer();
116
117
118 private:
119 virtual void beginJob() ;
120 virtual void analyze(const edm::Event&, const edm::EventSetup&);
121 virtual void endJob() ;
122 bool selectJet(int i);
123 // ----------member data ---------------------------
124
125 bool usePat_;
126 bool doMC_;
127 bool filterJets_;
128 bool diJetsOnly_;
129 bool matchDiJets_;
130 bool matchPatGen_;
131 bool matchNew_;
132 bool sortJets_;
133 bool correctJets_;
134 bool getFastJets_;
135
136 double matchR_;
137 double genPtMin_;
138 double ptMin_;
139 double emfMin_;
140 double n90Min_;
141 double n90hitMin_;
142
143 edm::InputTag jetTag_;
144 edm::InputTag matchTag_;
145 std::vector<edm::InputTag> matchTags_;
146
147 JRA jra_;
148 std::vector<JRA> jraMatch_;
149
150 TTree* t;
151
152 edm::Handle<edm::GenHIEvent> mc;
153 edm::Handle<reco::Centrality> cent;
154
155 edm::Handle<reco::JetView> jets;
156 edm::Handle<pat::JetCollection> patjets;
157
158 FactorizedJetCorrector* jetCorrector_;
159 edm::ESWatcher<JetCorrectionsRecord> watchJetCorrectionsRecord_;
160
161 std::string tags_;
162 std::string levels_;
163 std::string algo_;
164
165 edm::Service<TFileService> fs;
166
167 edm::Handle<vector<double> > ktRhos;
168 edm::Handle<vector<double> > akRhos;
169 bool doFastJets_;
170
171 vector<int> doMatchedFastJets_;
172 vector<int> correctMatchedJets_;
173
174 InputTag ktSrc_;
175 InputTag akSrc_;
176
177
178 };
179
180 bool HiJetResponseAnalyzer::selectJet(int i){
181 const reco::Jet& jet = (*jets)[i];
182 if(usePat_){
183 const pat::Jet& patjet = (*patjets)[i];
184 if(patjet.emEnergyFraction() <= emfMin_) return false;
185 if(patjet.jetID().n90Hits <= n90hitMin_) return false;
186 if(doMC_){
187
188 }
189
190 }
191
192 return true;
193 }
194
195
196 //
197 // constants, enums and typedefs
198 //
199
200 //
201 // static data member definitions
202 //
203
204 //
205 // constructors and destructor
206 //
207 HiJetResponseAnalyzer::HiJetResponseAnalyzer(const edm::ParameterSet& iConfig)
208
209 {
210
211 levels_ = iConfig.getUntrackedParameter<string>("corrLevels","L2Relative:L3Absolute");
212
213 algo_ = iConfig.getUntrackedParameter<string>("algo","IC5Calo");
214 tags_ = "";
215
216 string l[2] = {"L2Relative","L3Absolute"};
217
218 for(int i = 0; i <2; ++i){
219 edm::FileInPath fip("CondFormats/JetMETObjects/data/Spring10_"+l[i]+"_"+algo_+".txt");
220 tags_ += fip.fullPath();
221 if(i < 2 - 1)tags_ +=":";
222 }
223
224 jetCorrector_ = new FactorizedJetCorrector(levels_, tags_);
225
226 //now do what ever initialization is needed
227 matchR_ = iConfig.getUntrackedParameter<double>("matchR",0.25);
228
229 ptMin_ = iConfig.getUntrackedParameter<double>("ptMin",0);
230 genPtMin_ = iConfig.getUntrackedParameter<double>("genPtMin",20);
231 emfMin_ = iConfig.getUntrackedParameter<double>("emfMin",0.01);
232 n90Min_ = iConfig.getUntrackedParameter<double>("n90Min",1);
233 n90hitMin_ = iConfig.getUntrackedParameter<double>("n90hitMin",1);
234
235 filterJets_ = iConfig.getUntrackedParameter<bool>("filterJets",true);
236 diJetsOnly_ = iConfig.getUntrackedParameter<bool>("diJetsOnly",false);
237 matchDiJets_ = iConfig.getUntrackedParameter<bool>("matchDiJets",false);
238 matchPatGen_ = iConfig.getUntrackedParameter<bool>("matchPatGen",false);
239
240 matchNew_ = iConfig.getUntrackedParameter<bool>("matchNew",false);
241
242 usePat_ = iConfig.getUntrackedParameter<bool>("usePat",true);
243 doMC_ = iConfig.getUntrackedParameter<bool>("doMC",true);
244
245 sortJets_ = iConfig.getUntrackedParameter<bool>("sortJets",true);
246 correctJets_ = iConfig.getUntrackedParameter<bool>("correctJets",false);
247
248 correctMatchedJets_ = iConfig.getUntrackedParameter<std::vector<int> >("correctMatchedJets",std::vector<int>(0));
249 doMatchedFastJets_ = iConfig.getUntrackedParameter<std::vector<int> >("doMatchedFastJets",std::vector<int>(0));
250
251 jetTag_ = iConfig.getUntrackedParameter<edm::InputTag>("src",edm::InputTag("selectedPatJets"));
252 matchTag_ = iConfig.getUntrackedParameter<edm::InputTag>("match",edm::InputTag("selectedPatJets"));
253 matchTags_ = iConfig.getUntrackedParameter<std::vector<edm::InputTag> >("matches",std::vector<edm::InputTag>(0));
254
255 ktSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("ktSrc",edm::InputTag("kt4CaloJets"));
256 akSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("akSrc",edm::InputTag("ak5CaloJets"));
257 doFastJets_ = iConfig.getUntrackedParameter<bool>("doFastJets",false);
258
259 getFastJets_ = iConfig.getUntrackedParameter<bool>("getFastJets",false);
260
261 for(unsigned int i = 0; i < doMatchedFastJets_.size(); ++i){
262 getFastJets_ = getFastJets_ || (bool)doMatchedFastJets_[i];
263 }
264
265
266 }
267
268
269 HiJetResponseAnalyzer::~HiJetResponseAnalyzer()
270 {
271
272 // do anything here that needs to be done at desctruction time
273 // (e.g. close files, deallocate resources etc.)
274
275 }
276
277
278 //
279 // member functions
280 //
281
282 // ------------ method called to for each event ------------
283 void
284 HiJetResponseAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
285 {
286 using namespace edm;
287
288 iEvent.getByLabel(jetTag_,jets);
289 if(usePat_)iEvent.getByLabel(jetTag_,patjets);
290 std::vector<JRAV> jraV;
291
292 if(getFastJets_){
293 iEvent.getByLabel(edm::InputTag(ktSrc_.label(),"rhos"),ktRhos);
294 iEvent.getByLabel(edm::InputTag(akSrc_.label(),"rhos"),akRhos);
295 }
296
297 for(unsigned int j = 0 ; j < jets->size(); ++j){
298 if(filterJets_ && !selectJet(j)) continue;
299 const reco::Jet& jet = (*jets)[j];
300 JRAV jv;
301 jv.jtpt = jet.pt();
302 jv.jteta = jet.eta();
303 jv.jtphi = jet.phi();
304 jv.jtrawpt = jet.pt();
305 jv.area = jet.jetArea();
306 jv.pu = jet.pileup();
307 jv.index = j;
308
309 double pt = jet.pt();
310 double ktRho = -1, akRho = -1;
311 if(getFastJets_){
312 ktRho = getRho(jv.jteta,*ktRhos);
313 akRho = getRho(jv.jteta,*akRhos);
314 }
315
316 jv.rho = akRho;
317
318 if(doFastJets_){
319 jv.pu = jet.jetArea()*akRho;
320 pt -= jv.pu;
321 }
322 jv.jtpt = pt;
323
324 if(correctJets_){
325 jetCorrector_->setJetEta(jet.eta());
326 jetCorrector_->setJetPt(pt);
327 // jetCorrector_->setJetE(jet.energy());
328 vector<float> corrs = jetCorrector_->getSubCorrections();
329 jv.l2 = corrs[0];
330 jv.l3 = corrs[1];
331 jv.jtpt = pt*jv.l2*jv.l3;
332 }
333
334
335 if(usePat_){
336 const pat::Jet& patjet = (*patjets)[j];
337
338 jv.jtrawpt = patjet.correctedJet("raw").pt();
339 jv.jtpt = patjet.pt();
340
341 if(doMC_ && matchPatGen_ && patjet.genJet() != 0){
342 if(patjet.genJet()->pt() < genPtMin_) continue;
343 jv.refpt = patjet.genJet()->pt();
344 jv.refeta = patjet.genJet()->eta();
345 jv.refphi = patjet.genJet()->phi();
346 }else{
347 jv.refpt = -99;
348 jv.refeta = -99;
349 jv.refphi = -99;
350 }
351 }
352 jraV.push_back(jv);
353 }
354
355 if(sortJets_){
356 std::sort(jraV.begin(),jraV.end(),comparePt);
357 }
358
359 for(unsigned int i = 0; i < jraV.size(); ++i){
360 JRAV& jv = jraV[i];
361 const reco::Jet& jet = (*jets)[jv.index];
362
363 if(matchNew_){
364 for(unsigned int im = 0; im < matchTags_.size(); ++im){
365 edm::Handle<reco::JetView> matchedJets;
366 iEvent.getByLabel(matchTags_[im],matchedJets);
367 jraMatch_[im].jtrawpt[i] = -99;
368 jraMatch_[im].jtpt[i] = -99;
369 jraMatch_[im].jteta[i] = -99;
370 jraMatch_[im].jtphi[i] = -99;
371 jraMatch_[im].area[i] = -99;
372 jraMatch_[im].l2[i] = -99;
373 jraMatch_[im].l3[i] = -99;
374 jraMatch_[im].pu[i] = -99;
375 for(unsigned int m = 0 ; m < matchedJets->size(); ++m){
376 const reco::Jet& match = (*matchedJets)[m];
377 double dr = reco::deltaR(jet.eta(),jet.phi(),match.eta(),match.phi());
378 if(dr < matchR_ && match.pt() > genPtMin_){
379 jraMatch_[im].jtpt[i] = match.pt();
380 jraMatch_[im].jtrawpt[i] = match.pt();
381 jraMatch_[im].jteta[i] = match.eta();
382 jraMatch_[im].jtphi[i] = match.phi();
383 jraMatch_[im].area[i] = match.jetArea();
384
385 double ktRhoM = -1, akRhoM = -1;
386 if(getFastJets_){
387 ktRhoM = getRho(jraMatch_[im].jteta[i],*ktRhos);
388 akRhoM = getRho(jraMatch_[im].jteta[i],*akRhos);
389 }
390
391 jraMatch_[im].rho[i] = akRhoM;
392 double pt = jraMatch_[im].jtpt[i];
393
394 if((bool)doMatchedFastJets_[im]){
395 jraMatch_[im].pu[i] = jraMatch_[im].area[i]*akRhoM;
396 pt -= jraMatch_[im].pu[i];
397 }
398
399 jraMatch_[im].jtpt[i] = pt;
400
401 if((bool)correctMatchedJets_[im]){
402
403 jetCorrector_->setJetEta(jraMatch_[im].jteta[i]);
404 jetCorrector_->setJetPt(pt);
405 // jetCorrector_->setJetE(match.energy());
406 vector<float> corrs = jetCorrector_->getSubCorrections();
407 jraMatch_[im].l2[i] = corrs[0];
408 // Should it be re-set for L3???
409 jraMatch_[im].l3[i] = corrs[1];
410 jraMatch_[im].jtpt[i] = pt*jraMatch_[im].l2[i]*jraMatch_[im].l3[i];
411 }
412
413 }
414 }
415 }
416 }
417
418 jra_.jtpt[i] = jv.jtpt;
419 jra_.jteta[i] = jv.jteta;
420 jra_.jtphi[i] = jv.jtphi;
421 jra_.jtrawpt[i] = jv.jtrawpt;
422 jra_.refpt[i] = jv.refpt;
423 jra_.refeta[i] = jv.refeta;
424 jra_.refphi[i] = jv.refphi;
425
426 jra_.area[i] = jv.area;
427 jra_.pu[i] = jv.pu;
428 jra_.rho[i] = jv.rho;
429 jra_.l2[i] = jv.l2;
430 jra_.l3[i] = jv.l3;
431
432 }
433 jra_.nref = jraV.size();
434
435 t->Fill();
436
437 }
438
439 // ------------ method called once each job just before starting event loop ------------
440
441
442 void
443 HiJetResponseAnalyzer::beginJob(){
444 t= fs->make<TTree>("t","Jet Response Analyzer");
445 t->Branch("nref",&jra_.nref,"nref/I");
446 t->Branch("jtpt",jra_.jtpt,"jtpt[nref]/F");
447 t->Branch("jteta",jra_.jteta,"jteta[nref]/F");
448 t->Branch("jtphi",jra_.jtphi,"jtphi[nref]/F");
449 t->Branch("jtrawpt",jra_.jtrawpt,"jtrawpt[nref]/F");
450
451 if(correctJets_){
452 t->Branch("l2",jra_.l2,"l2[nref]/F");
453 t->Branch("l3",jra_.l3,"l3[nref]/F");
454 }
455 t->Branch("area",jra_.area,"area[nref]/F");
456 t->Branch("pu",jra_.pu,"pu[nref]/F");
457 t->Branch("rho",jra_.rho,"rho[nref]/F");
458
459 t->Branch("refpt",jra_.refpt,"refpt[nref]/F");
460 t->Branch("refcorpt",jra_.refpt,"refcorpt[nref]/F");
461 t->Branch("refeta",jra_.refeta,"refeta[nref]/F");
462 t->Branch("refphi",jra_.refphi,"refphi[nref]/F");
463 t->Branch("weight",&jra_.weight,"weight/F");
464
465 jraMatch_.clear();
466 for(unsigned int im = 0; im < matchTags_.size(); ++im){
467 JRA jrm;
468 jraMatch_.push_back(jrm);
469 }
470
471 for(unsigned int im = 0; im < matchTags_.size(); ++im){
472 t->Branch(Form("jtpt%d",im),jraMatch_[im].jtpt,Form("jtpt%d[nref]/F",im));
473 t->Branch(Form("jteta%d",im),jraMatch_[im].jteta,Form("jteta%d[nref]/F",im));
474 t->Branch(Form("jtphi%d",im),jraMatch_[im].jtphi,Form("jtphi%d[nref]/F",im));
475 t->Branch(Form("jtrawpt%d",im),jraMatch_[im].jtrawpt,Form("jtrawpt%d[nref]/F",im));
476
477 if((bool)correctMatchedJets_[im]){
478 t->Branch(Form("l2_%d",im),jraMatch_[im].l2,Form("l2_%d[nref]/F",im));
479 t->Branch(Form("l3_%d",im),jraMatch_[im].l3,Form("l3_%d[nref]/F",im));
480 }
481
482 t->Branch(Form("area%d",im),jraMatch_[im].area,Form("area%d[nref]/F",im));
483 t->Branch(Form("pu%d",im),jraMatch_[im].pu,Form("pu%d[nref]/F",im));
484 t->Branch(Form("rho%d",im),jraMatch_[im].rho,Form("rho%d[nref]/F",im));
485
486 }
487
488 }
489
490 // ------------ method called once each job just after ending the event loop ------------
491 void
492 HiJetResponseAnalyzer::endJob() {
493 }
494
495 //define this as a plug-in
496 DEFINE_FWK_MODULE(HiJetResponseAnalyzer);