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

# User Rev Content
1 yilmaz 1.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 yilmaz 1.2 // $Id: HiJetMatchAnalyzer.cc,v 1.1 2011/09/20 19:48:45 yilmaz Exp $
17 yilmaz 1.1 //
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 yilmaz 1.2 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 yilmaz 1.1 }
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);