ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/CmsHi/JetAnalysis/src/HiJetResponseAnalyzer.cc
(Generate patch)

Comparing UserCode/CmsHi/JetAnalysis/src/HiJetResponseAnalyzer.cc (file contents):
Revision 1.9 by yilmaz, Sun Oct 24 14:39:32 2010 UTC vs.
Revision 1.10 by yilmaz, Thu Jan 20 19:41:10 2011 UTC

# Line 40 | Line 40
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  
# Line 55 | Line 63 | struct etdr{
63     double dr;
64   };
65  
66 < struct JRA{
66 > class JRA{
67  
68 + public:
69     int nref;
70     int bin;
71     float b;
# Line 68 | Line 77 | struct JRA{
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   };
# Line 82 | Line 96 | struct JRAV{
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  
# Line 113 | Line 132 | class HiJetResponseAnalyzer : public edm
132    bool matchNew_;
133    bool sortJets_;
134    bool correctJets_;
135 +   bool getFastJets_;
136  
137    double matchR_;  
138     double genPtMin_;
# Line 135 | Line 155 | class HiJetResponseAnalyzer : public edm
155  
156     edm::Handle<reco::JetView> jets;
157     edm::Handle<pat::JetCollection> patjets;
158 <  edm::Handle<reco::JetView> matchedJets;
158 >
159 >   FactorizedJetCorrector* jetCorrector_;
160 >   edm::ESWatcher<JetCorrectionsRecord> watchJetCorrectionsRecord_;
161 >
162 >   std::string tags_;
163 >   std::string levels_;
164 >   std::string algo_;
165  
166     edm::Service<TFileService> fs;
167  
168 +   edm::Handle<vector<double> > ktRhos;
169 +   edm::Handle<vector<double> > akRhos;
170 +   bool doFastJets_;
171 +
172 +   vector<int> doMatchedFastJets_;
173 +   vector<int> correctMatchedJets_;
174 +  
175 +   InputTag ktSrc_;
176 +   InputTag akSrc_;
177 +
178 +
179   };
180  
181   bool HiJetResponseAnalyzer::selectJet(int i){
# Line 171 | Line 208 | bool HiJetResponseAnalyzer::selectJet(in
208   HiJetResponseAnalyzer::HiJetResponseAnalyzer(const edm::ParameterSet& iConfig)
209  
210   {
211 +
212 +   levels_ = iConfig.getUntrackedParameter<string>("corrLevels","L2Relative:L3Absolute");
213 +
214 +   algo_ = iConfig.getUntrackedParameter<string>("algo","IC5Calo");
215 +   tags_ = "";
216 +
217 +   string l[2] = {"L2Relative","L3Absolute"};
218 +
219 +   for(int i = 0; i <2; ++i){
220 +      edm::FileInPath fip("CondFormats/JetMETObjects/data/Spring10_"+l[i]+"_"+algo_+".txt");
221 +      tags_ += fip.fullPath();
222 +      if(i < 2 - 1)tags_ +=":";
223 +   }
224 +
225 +   jetCorrector_ = new FactorizedJetCorrector(levels_, tags_);
226 +
227     //now do what ever initialization is needed
228    matchR_ = iConfig.getUntrackedParameter<double>("matchR",0.25);
229  
# Line 193 | Line 246 | HiJetResponseAnalyzer::HiJetResponseAnal
246     sortJets_ = iConfig.getUntrackedParameter<bool>("sortJets",true);
247     correctJets_ = iConfig.getUntrackedParameter<bool>("correctJets",false);
248  
249 +   correctMatchedJets_ = iConfig.getUntrackedParameter<std::vector<int> >("correctMatchedJets",std::vector<int>(0));
250 +   doMatchedFastJets_ = iConfig.getUntrackedParameter<std::vector<int> >("doMatchedFastJets",std::vector<int>(0));
251 +
252     jetTag_ = iConfig.getUntrackedParameter<edm::InputTag>("src",edm::InputTag("selectedPatJets"));
253     matchTag_ = iConfig.getUntrackedParameter<edm::InputTag>("match",edm::InputTag("selectedPatJets"));
254     matchTags_ = iConfig.getUntrackedParameter<std::vector<edm::InputTag> >("matches",std::vector<edm::InputTag>(0));
255 +
256 +   ktSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("ktSrc",edm::InputTag("kt4CaloJets"));
257 +   akSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("akSrc",edm::InputTag("ak5CaloJets"));
258 +   doFastJets_ = iConfig.getUntrackedParameter<bool>("doFastJets",false);
259 +
260 +   getFastJets_ = true;
261 +   for(unsigned int i = 0; i < doMatchedFastJets_.size(); ++i){
262 +      getFastJets_ = getFastJets_ || (bool)doMatchedFastJets_[i];
263 +   }
264 +
265 +
266   }
267  
268  
# Line 222 | Line 289 | HiJetResponseAnalyzer::analyze(const edm
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];
# Line 230 | Line 302 | HiJetResponseAnalyzer::analyze(const edm
302       jv.jteta = jet.eta();
303       jv.jtphi = jet.phi();
304       jv.jtcorpt = 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 = getRho(jv.jteta,*ktRhos);
311 +     double akRho = getRho(jv.jteta,*akRhos);
312 +
313 +     jv.rho = ktRho;
314 +     if(correctJets_){
315 +        if(doFastJets_){
316 +           jv.pu = jet.jetArea()*ktRho;
317 +           pt -= jv.pu;
318 +        }
319 +
320 +        jetCorrector_->setJetEta(jet.eta());
321 +        jetCorrector_->setJetPt(pt);
322 +        //      jetCorrector_->setJetE(jet.energy());
323 +        vector<float> corrs = jetCorrector_->getSubCorrections();
324 +        jv.l2 = corrs[0];
325 +        jv.l3 = corrs[1];
326 +        jv.jtcorpt = pt*jv.l2*jv.l3;
327 +
328 +     }
329       if(usePat_){
330         const pat::Jet& patjet = (*patjets)[j];
331  
# Line 248 | Line 343 | HiJetResponseAnalyzer::analyze(const edm
343           jv.refphi = -99;
344         }
345       }
251 –
346       jraV.push_back(jv);
347     }
348  
# Line 263 | Line 357 | HiJetResponseAnalyzer::analyze(const edm
357      
358       if(matchNew_){
359         for(unsigned int im = 0; im < matchTags_.size(); ++im){
360 <         iEvent.getByLabel(matchTags_[im],matchedJets);
360 >          edm::Handle<reco::JetView> matchedJets;
361 >          iEvent.getByLabel(matchTags_[im],matchedJets);
362 >         jraMatch_[im].jtcorpt[i] = -99;
363 >         jraMatch_[im].jtpt[i] = -99;
364 >         jraMatch_[im].jteta[i] = -99;
365 >         jraMatch_[im].jtphi[i] = -99;
366 >         jraMatch_[im].area[i] = -99;
367 >         jraMatch_[im].l2[i] = -99;
368 >         jraMatch_[im].l3[i] = -99;
369 >         jraMatch_[im].pu[i] = -99;
370           for(unsigned int m = 0 ; m < matchedJets->size(); ++m){
371             const reco::Jet& match = (*matchedJets)[m];
372             double dr = reco::deltaR(jet.eta(),jet.phi(),match.eta(),match.phi());
373 <           if(dr < matchR_){
373 >           if(dr < matchR_ && match.pt() > genPtMin_){
374               jraMatch_[im].jtcorpt[i] = -99;
375               jraMatch_[im].jtpt[i] = match.pt();
376               jraMatch_[im].jteta[i] = match.eta();
377               jraMatch_[im].jtphi[i] = match.phi();
378 +             jraMatch_[im].area[i] = match.jetArea();
379 +
380 +             double  ktRho = getRho(jraMatch_[im].jteta[i],*ktRhos);
381 +             double akRho  = getRho(jraMatch_[im].jteta[i],*akRhos);
382 +
383 +             jraMatch_[im].rho[i] = ktRho;
384 +
385 +             if((bool)correctMatchedJets_[im]){
386 +
387 +                double ktRho;
388 +                double pt = jraMatch_[im].jtpt[i];
389 +
390 +                if((bool)doMatchedFastJets_[im]){
391 +                   jraMatch_[im].pu[i] = jraMatch_[im].area[i]*ktRho;
392 +                   pt -= jraMatch_[im].pu[i];
393 +                }
394 +                jetCorrector_->setJetEta(jraMatch_[im].jteta[i]);
395 +                jetCorrector_->setJetPt(pt);
396 +                //              jetCorrector_->setJetE(match.energy());
397 +                vector<float> corrs = jetCorrector_->getSubCorrections();
398 +                jraMatch_[im].l2[i] = corrs[0];
399 +                // Should it be re-set for L3???
400 +                jraMatch_[im].l3[i] = corrs[1];
401 +                jraMatch_[im].jtcorpt[i] = pt*jraMatch_[im].l2[i]*jraMatch_[im].l3[i];
402 +
403 +             }
404 +
405 +
406             }
407           }
408         }
# Line 284 | Line 415 | HiJetResponseAnalyzer::analyze(const edm
415       jra_.refpt[i] = jv.refpt;
416       jra_.refeta[i] = jv.refeta;
417       jra_.refphi[i] = jv.refphi;
418 +
419 +     jra_.area[i] = jv.area;
420 +     jra_.pu[i] = jv.pu;
421 +     jra_.rho[i] = jv.rho;
422 +     jra_.l2[i] = jv.l2;
423 +     jra_.l3[i] = jv.l3;
424 +
425     }
426     jra_.nref = jraV.size();
427    
# Line 299 | Line 437 | HiJetResponseAnalyzer::beginJob(){
437    t= fs->make<TTree>("t","Jet Response Analyzer");
438    t->Branch("nref",&jra_.nref,"nref/I");
439    t->Branch("jtpt",jra_.jtpt,"jtpt[nref]/F");
302 –  t->Branch("jtcorpt",jra_.jtcorpt,"jtcorpt[nref]/F");
440    t->Branch("jteta",jra_.jteta,"jteta[nref]/F");
441    t->Branch("jtphi",jra_.jtphi,"jtphi[nref]/F");
442 +
443 +  if(correctJets_){
444 +     t->Branch("jtcorpt",jra_.jtcorpt,"jtcorpt[nref]/F");
445 +     t->Branch("l2",jra_.l2,"l2[nref]/F");
446 +     t->Branch("l3",jra_.l3,"l3[nref]/F");
447 +  }
448 +  t->Branch("area",jra_.area,"area[nref]/F");
449 +  t->Branch("pu",jra_.pu,"pu[nref]/F");
450 +  t->Branch("rho",jra_.rho,"rho[nref]/F");
451 +
452    t->Branch("refpt",jra_.refpt,"refpt[nref]/F");
453    t->Branch("refcorpt",jra_.refpt,"refcorpt[nref]/F");
454     t->Branch("refeta",jra_.refeta,"refeta[nref]/F");
455     t->Branch("refphi",jra_.refphi,"refphi[nref]/F");
456     t->Branch("weight",&jra_.weight,"weight/F");
457 +
458 +   jraMatch_.clear();
459     for(unsigned int im = 0; im < matchTags_.size(); ++im){
460       JRA jrm;
461       jraMatch_.push_back(jrm);
462 +   }
463 +
464 +   for(unsigned int im = 0; im < matchTags_.size(); ++im){
465       t->Branch(Form("jtpt%d",im),jraMatch_[im].jtpt,Form("jtpt%d[nref]/F",im));
314 –     t->Branch(Form("jtcorpt%d",im),jraMatch_[im].jtcorpt,Form("jtcorpt%d[nref]/F",im));
466       t->Branch(Form("jteta%d",im),jraMatch_[im].jteta,Form("jteta%d[nref]/F",im));
467       t->Branch(Form("jtphi%d",im),jraMatch_[im].jtphi,Form("jtphi%d[nref]/F",im));
468 +
469 +     if((bool)correctMatchedJets_[im]){
470 +        t->Branch(Form("jtcorpt%d",im),jraMatch_[im].jtcorpt,Form("jtcorpt%d[nref]/F",im));
471 +        t->Branch(Form("l2_%d",im),jraMatch_[im].l2,Form("l2_%d[nref]/F",im));
472 +        t->Branch(Form("l3_%d",im),jraMatch_[im].l3,Form("l3_%d[nref]/F",im));
473 +     }
474 +    
475 +     t->Branch(Form("area%d",im),jraMatch_[im].area,Form("area%d[nref]/F",im));
476 +     t->Branch(Form("pu%d",im),jraMatch_[im].pu,Form("pu%d[nref]/F",im));
477 +     t->Branch(Form("rho%d",im),jraMatch_[im].rho,Form("rho%d[nref]/F",im));
478 +
479     }
480    
481   }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines