ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/CmsHi/JetAnalysis/src/HiJetMatchAnalyzer.cc
Revision: 1.1
Committed: Tue Sep 20 19:48:45 2011 UTC (13 years, 7 months ago) by yilmaz
Content type: text/plain
Branch: MAIN
CVS Tags: hi44X_02, hi413_03, hi413_02, hi441_1, hi441_0, hi413_11, hi413_10, hi413_09
Log Message:
rename the matching analyzer

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     // $Id: HiJetMatchAnalyzer.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 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     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     HiJetMatchAnalyzer::~HiJetMatchAnalyzer()
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     HiJetMatchAnalyzer::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     HiJetMatchAnalyzer::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     HiJetMatchAnalyzer::endJob() {
493     }
494    
495     //define this as a plug-in
496     DEFINE_FWK_MODULE(HiJetMatchAnalyzer);