ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/CmsHi/JetAnalysis/src/HiJetMatchAnalyzer.cc
Revision: 1.2
Committed: Sun Apr 22 19:12:44 2012 UTC (13 years 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, HEAD
Branch point for: branch_44x
Changes since 1.1: +22 -16 lines
Log Message:
up

File Contents

# Content
1 // -*- C++ -*-
2 //
3 // Package: HiJetMatchAnalyzer
4 // Class: HiJetMatchAnalyzer
5 //
6 /**\class HiJetMatchAnalyzer HiJetMatchAnalyzer.cc CmsHi/HiJetMatchAnalyzer/src/HiJetMatchAnalyzer.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: HiJetMatchAnalyzer.cc,v 1.1 2011/09/20 19:48:45 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 HiJetMatchAnalyzer : public edm::EDAnalyzer {
113 public:
114 explicit HiJetMatchAnalyzer(const edm::ParameterSet&);
115 ~HiJetMatchAnalyzer();
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 HiJetMatchAnalyzer::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 HiJetMatchAnalyzer::HiJetMatchAnalyzer(const edm::ParameterSet& iConfig)
208
209 {
210
211 //now do what ever initialization is needed
212 matchR_ = iConfig.getUntrackedParameter<double>("matchR",0.25);
213
214 ptMin_ = iConfig.getUntrackedParameter<double>("ptMin",0);
215 genPtMin_ = iConfig.getUntrackedParameter<double>("genPtMin",20);
216 emfMin_ = iConfig.getUntrackedParameter<double>("emfMin",0.01);
217 n90Min_ = iConfig.getUntrackedParameter<double>("n90Min",1);
218 n90hitMin_ = iConfig.getUntrackedParameter<double>("n90hitMin",1);
219
220 filterJets_ = iConfig.getUntrackedParameter<bool>("filterJets",true);
221 diJetsOnly_ = iConfig.getUntrackedParameter<bool>("diJetsOnly",false);
222 matchDiJets_ = iConfig.getUntrackedParameter<bool>("matchDiJets",false);
223 matchPatGen_ = iConfig.getUntrackedParameter<bool>("matchPatGen",false);
224
225 matchNew_ = iConfig.getUntrackedParameter<bool>("matchNew",false);
226
227 usePat_ = iConfig.getUntrackedParameter<bool>("usePat",true);
228 doMC_ = iConfig.getUntrackedParameter<bool>("doMC",true);
229
230 sortJets_ = iConfig.getUntrackedParameter<bool>("sortJets",true);
231 correctJets_ = iConfig.getUntrackedParameter<bool>("correctJets",false);
232
233 correctMatchedJets_ = iConfig.getUntrackedParameter<std::vector<int> >("correctMatchedJets",std::vector<int>(0));
234 doMatchedFastJets_ = iConfig.getUntrackedParameter<std::vector<int> >("doMatchedFastJets",std::vector<int>(0));
235
236 jetTag_ = iConfig.getUntrackedParameter<edm::InputTag>("src",edm::InputTag("selectedPatJets"));
237 matchTag_ = iConfig.getUntrackedParameter<edm::InputTag>("match",edm::InputTag("selectedPatJets"));
238 matchTags_ = iConfig.getUntrackedParameter<std::vector<edm::InputTag> >("matches",std::vector<edm::InputTag>(0));
239
240 ktSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("ktSrc",edm::InputTag("kt4CaloJets"));
241 akSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("akSrc",edm::InputTag("ak5CaloJets"));
242 doFastJets_ = iConfig.getUntrackedParameter<bool>("doFastJets",false);
243
244 getFastJets_ = iConfig.getUntrackedParameter<bool>("getFastJets",false);
245
246 for(unsigned int i = 0; i < doMatchedFastJets_.size(); ++i){
247 getFastJets_ = getFastJets_ || (bool)doMatchedFastJets_[i];
248 }
249
250
251 if(correctJets_){
252
253 levels_ = iConfig.getUntrackedParameter<string>("corrLevels","L2Relative:L3Absolute");
254
255 algo_ = iConfig.getUntrackedParameter<string>("algo","IC5Calo");
256 tags_ = "";
257
258 string l[2] = {"L2Relative","L3Absolute"};
259
260 for(int i = 0; i <2; ++i){
261 edm::FileInPath fip("CondFormats/JetMETObjects/data/Spring10_"+l[i]+"_"+algo_+".txt");
262 tags_ += fip.fullPath();
263 if(i < 2 - 1)tags_ +=":";
264 }
265
266 jetCorrector_ = new FactorizedJetCorrector(levels_, tags_);
267 }
268
269
270
271
272 }
273
274
275 HiJetMatchAnalyzer::~HiJetMatchAnalyzer()
276 {
277
278 // do anything here that needs to be done at desctruction time
279 // (e.g. close files, deallocate resources etc.)
280
281 }
282
283
284 //
285 // member functions
286 //
287
288 // ------------ method called to for each event ------------
289 void
290 HiJetMatchAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
291 {
292 using namespace edm;
293
294 iEvent.getByLabel(jetTag_,jets);
295 if(usePat_)iEvent.getByLabel(jetTag_,patjets);
296 std::vector<JRAV> jraV;
297
298 if(getFastJets_){
299 iEvent.getByLabel(edm::InputTag(ktSrc_.label(),"rhos"),ktRhos);
300 iEvent.getByLabel(edm::InputTag(akSrc_.label(),"rhos"),akRhos);
301 }
302
303 for(unsigned int j = 0 ; j < jets->size(); ++j){
304 if(filterJets_ && !selectJet(j)) continue;
305 const reco::Jet& jet = (*jets)[j];
306 JRAV jv;
307 jv.jtpt = jet.pt();
308 jv.jteta = jet.eta();
309 jv.jtphi = jet.phi();
310 jv.jtrawpt = jet.pt();
311 jv.area = jet.jetArea();
312 jv.pu = jet.pileup();
313 jv.index = j;
314
315 double pt = jet.pt();
316 double ktRho = -1, akRho = -1;
317 if(getFastJets_){
318 ktRho = getRho(jv.jteta,*ktRhos);
319 akRho = getRho(jv.jteta,*akRhos);
320 }
321
322 jv.rho = akRho;
323
324 if(doFastJets_){
325 jv.pu = jet.jetArea()*akRho;
326 pt -= jv.pu;
327 }
328 jv.jtpt = pt;
329
330 if(correctJets_){
331 jetCorrector_->setJetEta(jet.eta());
332 jetCorrector_->setJetPt(pt);
333 // jetCorrector_->setJetE(jet.energy());
334 vector<float> corrs = jetCorrector_->getSubCorrections();
335 jv.l2 = corrs[0];
336 jv.l3 = corrs[1];
337 jv.jtpt = pt*jv.l2*jv.l3;
338 }
339
340
341 if(usePat_){
342 const pat::Jet& patjet = (*patjets)[j];
343
344 jv.jtrawpt = patjet.correctedJet("raw").pt();
345 jv.jtpt = patjet.pt();
346
347 if(doMC_ && matchPatGen_ && patjet.genJet() != 0){
348 if(patjet.genJet()->pt() < genPtMin_) continue;
349 jv.refpt = patjet.genJet()->pt();
350 jv.refeta = patjet.genJet()->eta();
351 jv.refphi = patjet.genJet()->phi();
352 }else{
353 jv.refpt = -99;
354 jv.refeta = -99;
355 jv.refphi = -99;
356 }
357 }
358 jraV.push_back(jv);
359 }
360
361 if(sortJets_){
362 std::sort(jraV.begin(),jraV.end(),comparePt);
363 }
364
365 for(unsigned int i = 0; i < jraV.size(); ++i){
366 JRAV& jv = jraV[i];
367 const reco::Jet& jet = (*jets)[jv.index];
368
369 if(matchNew_){
370 for(unsigned int im = 0; im < matchTags_.size(); ++im){
371 edm::Handle<reco::JetView> matchedJets;
372 iEvent.getByLabel(matchTags_[im],matchedJets);
373 jraMatch_[im].jtrawpt[i] = -99;
374 jraMatch_[im].jtpt[i] = -99;
375 jraMatch_[im].jteta[i] = -99;
376 jraMatch_[im].jtphi[i] = -99;
377 jraMatch_[im].area[i] = -99;
378 jraMatch_[im].l2[i] = -99;
379 jraMatch_[im].l3[i] = -99;
380 jraMatch_[im].pu[i] = -99;
381 for(unsigned int m = 0 ; m < matchedJets->size(); ++m){
382 const reco::Jet& match = (*matchedJets)[m];
383 double dr = reco::deltaR(jet.eta(),jet.phi(),match.eta(),match.phi());
384 if(dr < matchR_ && match.pt() > genPtMin_){
385 jraMatch_[im].jtpt[i] = match.pt();
386 jraMatch_[im].jtrawpt[i] = match.pt();
387 jraMatch_[im].jteta[i] = match.eta();
388 jraMatch_[im].jtphi[i] = match.phi();
389 jraMatch_[im].area[i] = match.jetArea();
390
391 double ktRhoM = -1, akRhoM = -1;
392 if(getFastJets_){
393 ktRhoM = getRho(jraMatch_[im].jteta[i],*ktRhos);
394 akRhoM = getRho(jraMatch_[im].jteta[i],*akRhos);
395 }
396
397 jraMatch_[im].rho[i] = akRhoM;
398 double pt = jraMatch_[im].jtpt[i];
399
400 if((bool)doMatchedFastJets_[im]){
401 jraMatch_[im].pu[i] = jraMatch_[im].area[i]*akRhoM;
402 pt -= jraMatch_[im].pu[i];
403 }
404
405 jraMatch_[im].jtpt[i] = pt;
406
407 if((bool)correctMatchedJets_[im]){
408
409 jetCorrector_->setJetEta(jraMatch_[im].jteta[i]);
410 jetCorrector_->setJetPt(pt);
411 // jetCorrector_->setJetE(match.energy());
412 vector<float> corrs = jetCorrector_->getSubCorrections();
413 jraMatch_[im].l2[i] = corrs[0];
414 // Should it be re-set for L3???
415 jraMatch_[im].l3[i] = corrs[1];
416 jraMatch_[im].jtpt[i] = pt*jraMatch_[im].l2[i]*jraMatch_[im].l3[i];
417 }
418
419 }
420 }
421 }
422 }
423
424 jra_.jtpt[i] = jv.jtpt;
425 jra_.jteta[i] = jv.jteta;
426 jra_.jtphi[i] = jv.jtphi;
427 jra_.jtrawpt[i] = jv.jtrawpt;
428 jra_.refpt[i] = jv.refpt;
429 jra_.refeta[i] = jv.refeta;
430 jra_.refphi[i] = jv.refphi;
431
432 jra_.area[i] = jv.area;
433 jra_.pu[i] = jv.pu;
434 jra_.rho[i] = jv.rho;
435 jra_.l2[i] = jv.l2;
436 jra_.l3[i] = jv.l3;
437
438 }
439 jra_.nref = jraV.size();
440
441 t->Fill();
442
443 }
444
445 // ------------ method called once each job just before starting event loop ------------
446
447
448 void
449 HiJetMatchAnalyzer::beginJob(){
450 t= fs->make<TTree>("t","Jet Response Analyzer");
451 t->Branch("nref",&jra_.nref,"nref/I");
452 t->Branch("jtpt",jra_.jtpt,"jtpt[nref]/F");
453 t->Branch("jteta",jra_.jteta,"jteta[nref]/F");
454 t->Branch("jtphi",jra_.jtphi,"jtphi[nref]/F");
455 t->Branch("jtrawpt",jra_.jtrawpt,"jtrawpt[nref]/F");
456
457 if(correctJets_){
458 t->Branch("l2",jra_.l2,"l2[nref]/F");
459 t->Branch("l3",jra_.l3,"l3[nref]/F");
460 }
461 t->Branch("area",jra_.area,"area[nref]/F");
462 t->Branch("pu",jra_.pu,"pu[nref]/F");
463 t->Branch("rho",jra_.rho,"rho[nref]/F");
464
465 t->Branch("refpt",jra_.refpt,"refpt[nref]/F");
466 t->Branch("refcorpt",jra_.refpt,"refcorpt[nref]/F");
467 t->Branch("refeta",jra_.refeta,"refeta[nref]/F");
468 t->Branch("refphi",jra_.refphi,"refphi[nref]/F");
469 t->Branch("weight",&jra_.weight,"weight/F");
470
471 jraMatch_.clear();
472 for(unsigned int im = 0; im < matchTags_.size(); ++im){
473 JRA jrm;
474 jraMatch_.push_back(jrm);
475 }
476
477 for(unsigned int im = 0; im < matchTags_.size(); ++im){
478 t->Branch(Form("jtpt%d",im),jraMatch_[im].jtpt,Form("jtpt%d[nref]/F",im));
479 t->Branch(Form("jteta%d",im),jraMatch_[im].jteta,Form("jteta%d[nref]/F",im));
480 t->Branch(Form("jtphi%d",im),jraMatch_[im].jtphi,Form("jtphi%d[nref]/F",im));
481 t->Branch(Form("jtrawpt%d",im),jraMatch_[im].jtrawpt,Form("jtrawpt%d[nref]/F",im));
482
483 if((bool)correctMatchedJets_[im]){
484 t->Branch(Form("l2_%d",im),jraMatch_[im].l2,Form("l2_%d[nref]/F",im));
485 t->Branch(Form("l3_%d",im),jraMatch_[im].l3,Form("l3_%d[nref]/F",im));
486 }
487
488 t->Branch(Form("area%d",im),jraMatch_[im].area,Form("area%d[nref]/F",im));
489 t->Branch(Form("pu%d",im),jraMatch_[im].pu,Form("pu%d[nref]/F",im));
490 t->Branch(Form("rho%d",im),jraMatch_[im].rho,Form("rho%d[nref]/F",im));
491
492 }
493
494 }
495
496 // ------------ method called once each job just after ending the event loop ------------
497 void
498 HiJetMatchAnalyzer::endJob() {
499 }
500
501 //define this as a plug-in
502 DEFINE_FWK_MODULE(HiJetMatchAnalyzer);