[ViewVC] Diff of: cvsroot/UserCode/CmsHi/JetAnalysis/src/HiJetResponseAnalyzer.cc

Comparing UserCode/CmsHi/JetAnalysis/src/HiJetResponseAnalyzer.cc (file contents):
Revision 1.2 by yilmaz, Fri Sep 10 12:54:47 2010 UTC vs.
Revision 1.10 by yilmaz, Thu Jan 20 19:41:10 2011 UTC

+#include "DataFormats/JetReco/interface/CaloJetCollection.h"
+#include "DataFormats/JetReco/interface/GenJetCollection.h"
-+
+#include "CondFormats/JetMETObjects/interface/FactorizedJetCorrector.h"
-+
+#include "JetMETCorrections/Objects/interface/JetCorrectionsRecord.h"
-+
+#include "CondFormats/JetMETObjects/interface/JetCorrectionUncertainty.h"
-+
+#include "CommonTools/UtilAlgos/interface/TFileService.h"
+#include "FWCore/ServiceRegistry/interface/Service.h"
-+
+#include "DataFormats/Math/interface/deltaR.h"
-+
+#include "FWCore/ParameterSet/interface/FileInPath.h"
-+
+#include "FWCore/Framework/interface/ESWatcher.h"
-<
->
+#include "CmsHi/JetAnalysis/interface/RhoGetter.h"
+#include "TTree.h"
+using namespace std;
-+
+using namespace edm;
+static const int MAXJETS = 500;
+   double dr;
+};
-<
+struct JRA{
->
+class JRA{
-+
+public:
+   int nref;
+   int bin;
+   float b;
+   float hf;
+   float jtpt[MAXJETS];
-+
+   float jtcorpt[MAXJETS];
+   float refpt[MAXJETS];
+   float jteta[MAXJETS];
+   float refeta[MAXJETS];
+   float jtphi[MAXJETS];
+   float refphi[MAXJETS];
-+
+   float l2[MAXJETS];
-+
+   float l3[MAXJETS];
-+
+   float area[MAXJETS];
-+
+   float pu[MAXJETS];
-+
+   float rho[MAXJETS];
+   float weight;
+};
-+
+struct JRAV{
-+
+  int index;
-+
+  float jtpt;
-+
+  float jtcorpt;
-+
+  float refpt;
-+
+  float refcorpt;
-+
+  float jteta;
-+
+  float refeta;
-+
+  float jtphi;
-+
+  float refphi;
-+
+   float l2;
-+
+   float l3;
-+
+   float area;
-+
+   float pu;
-+
+   float rho;
-+
-+
+};
-+
+//
+// class declaration
+//
-+
+bool compareCorPt(JRAV a, JRAV b) {return a.jtcorpt > b.jtcorpt;}
-+
+bool comparePt(JRAV a, JRAV b) {return a.jtpt > b.jtpt;}
+class HiJetResponseAnalyzer : public edm::EDAnalyzer {
+   public:
+      virtual void beginJob() ;
+      virtual void analyze(const edm::Event&, const edm::EventSetup&);
+      virtual void endJob() ;
-<
->
+   bool selectJet(int i);
+      // ----------member data ---------------------------
-<
+   bool usePat_;
-<
+   bool doMC_;
-<
-<
+   double genJetPtMin_;
->
+  bool usePat_;
->
+  bool doMC_;
->
+  bool filterJets_;
->
+  bool diJetsOnly_;
->
+  bool matchDiJets_;
->
+  bool matchPatGen_;
->
+  bool matchNew_;
->
+  bool sortJets_;
->
+  bool correctJets_;
->
+   bool getFastJets_;
->
->
+  double matchR_;
->
+   double genPtMin_;
->
+   double ptMin_;
->
+   double emfMin_;
->
+   double n90Min_;
->
+   double n90hitMin_;
->
->
+  edm::InputTag jetTag_;
->
+  edm::InputTag matchTag_;
->
+  std::vector<edm::InputTag> matchTags_;
->
->
+  JRA jra_;
->
+  std::vector<JRA> jraMatch_;
-–
-–
+   edm::InputTag jetTag_;
-–
-–
+   JRA jra_;
+   TTree* t;
+   edm::Handle<edm::GenHIEvent> mc;
+   edm::Handle<reco::Centrality> cent;
+   edm::Handle<reco::JetView> jets;
-<
+   //   edm::Handle<pat::JetCollection> jets;
->
+   edm::Handle<pat::JetCollection> patjets;
->
->
+   FactorizedJetCorrector* jetCorrector_;
->
+   edm::ESWatcher<JetCorrectionsRecord> watchJetCorrectionsRecord_;
-+
+   std::string tags_;
-+
+   std::string levels_;
-+
+   std::string algo_;
+   edm::Service<TFileService> fs;
-+
+   edm::Handle<vector<double> > ktRhos;
-+
+   edm::Handle<vector<double> > akRhos;
-+
+   bool doFastJets_;
-+
-+
+   vector<int> doMatchedFastJets_;
-+
+   vector<int> correctMatchedJets_;
-+
-+
+   InputTag ktSrc_;
-+
+   InputTag akSrc_;
-+
-+
+};
-+
+bool HiJetResponseAnalyzer::selectJet(int i){
-+
+   const reco::Jet& jet = (*jets)[i];
-+
+   if(usePat_){
-+
+      const pat::Jet& patjet = (*patjets)[i];
-+
+      if(patjet.emEnergyFraction() <= emfMin_) return false;
-+
+      if(patjet.jetID().n90Hits <= n90hitMin_) return false;
-+
+      if(doMC_){
-+
-+
+      }
-+
-+
+   }
-+
-+
+   return true;
-+
+}
-+
-+
+//
+// constants, enums and typedefs
+//
+HiJetResponseAnalyzer::HiJetResponseAnalyzer(const edm::ParameterSet& iConfig)
-+
-+
+   levels_ = iConfig.getUntrackedParameter<string>("corrLevels","L2Relative:L3Absolute");
-+
-+
+   algo_ = iConfig.getUntrackedParameter<string>("algo","IC5Calo");
-+
+   tags_ = "";
-+
-+
+   string l[2] = {"L2Relative","L3Absolute"};
-+
-+
+   for(int i = 0; i <2; ++i){
-+
+      edm::FileInPath fip("CondFormats/JetMETObjects/data/Spring10_"+l[i]+"_"+algo_+".txt");
-+
+      tags_ += fip.fullPath();
-+
+      if(i < 2 - 1)tags_ +=":";
-+
+   }
-+
-+
+   jetCorrector_ = new FactorizedJetCorrector(levels_, tags_);
-+
+   //now do what ever initialization is needed
-<
+   genJetPtMin_ = 0;
->
+  matchR_ = iConfig.getUntrackedParameter<double>("matchR",0.25);
->
->
+   ptMin_ = iConfig.getUntrackedParameter<double>("ptMin",0);
->
+   genPtMin_ = iConfig.getUntrackedParameter<double>("genPtMin",20);
->
+   emfMin_ = iConfig.getUntrackedParameter<double>("emfMin",0.01);
->
+   n90Min_ = iConfig.getUntrackedParameter<double>("n90Min",1);
->
+   n90hitMin_ = iConfig.getUntrackedParameter<double>("n90hitMin",1);
->
->
+   filterJets_ = iConfig.getUntrackedParameter<bool>("filterJets",true);
->
+   diJetsOnly_ = iConfig.getUntrackedParameter<bool>("diJetsOnly",false);
->
+   matchDiJets_ = iConfig.getUntrackedParameter<bool>("matchDiJets",false);
->
+   matchPatGen_ = iConfig.getUntrackedParameter<bool>("matchPatGen",false);
->
->
+   matchNew_ = iConfig.getUntrackedParameter<bool>("matchNew",false);
->
->
+   usePat_ = iConfig.getUntrackedParameter<bool>("usePat",true);
->
+   doMC_ = iConfig.getUntrackedParameter<bool>("doMC",true);
->
->
+   sortJets_ = iConfig.getUntrackedParameter<bool>("sortJets",true);
->
+   correctJets_ = iConfig.getUntrackedParameter<bool>("correctJets",false);
-<
+   usePat_ = false;
-<
+   doMC_ = false;
->
+   correctMatchedJets_ = iConfig.getUntrackedParameter<std::vector<int> >("correctMatchedJets",std::vector<int>(0));
->
+   doMatchedFastJets_ = iConfig.getUntrackedParameter<std::vector<int> >("doMatchedFastJets",std::vector<int>(0));
->
->
+   jetTag_ = iConfig.getUntrackedParameter<edm::InputTag>("src",edm::InputTag("selectedPatJets"));
->
+   matchTag_ = iConfig.getUntrackedParameter<edm::InputTag>("match",edm::InputTag("selectedPatJets"));
->
+   matchTags_ = iConfig.getUntrackedParameter<std::vector<edm::InputTag> >("matches",std::vector<edm::InputTag>(0));
->
->
+   ktSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("ktSrc",edm::InputTag("kt4CaloJets"));
->
+   akSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("akSrc",edm::InputTag("ak5CaloJets"));
->
+   doFastJets_ = iConfig.getUntrackedParameter<bool>("doFastJets",false);
->
->
+   getFastJets_ = true;
->
+   for(unsigned int i = 0; i < doMatchedFastJets_.size(); ++i){
->
+      getFastJets_ = getFastJets_ || (bool)doMatchedFastJets_[i];
->
+   }
-–
+   jetTag_ = edm::InputTag("selectedPatJets");
+   using namespace edm;
+   iEvent.getByLabel(jetTag_,jets);
-<
+   cout<<"a"<<endl;
->
+   if(usePat_)iEvent.getByLabel(jetTag_,patjets);
->
+   std::vector<JRAV> jraV;
->
->
+   if(getFastJets_){
->
+      iEvent.getByLabel(edm::InputTag(ktSrc_.label(),"rhos"),ktRhos);
->
+      iEvent.getByLabel(edm::InputTag(akSrc_.label(),"rhos"),akRhos);
->
+   }
-–
+   jra_.nref = 0;
+   for(unsigned int j = 0 ; j < jets->size(); ++j){
-<
+      cout<<"b"<<endl;
-<
-<
+      //      const pat::Jet& jet = (*jets)[j];
-<
+      const reco::Jet jet = (*jets)[j];
-<
+      jra_.jtpt[jra_.nref] = jet.pt();
-<
+      jra_.jteta[jra_.nref] = jet.eta();
-<
+      jra_.jtphi[jra_.nref] = jet.phi();
-<
+      cout<<"c"<<endl;
-<
-<
+      if(usePat_){
-<
+         const pat::Jet& patjet = (*jets)[j];
-<
-<
+         if(doMC_ && patjet.genJet() != 0){
-<
+            //      if(jet.genJet()->pt() < genJetPtMin_) continue;
-<
+            jra_.refpt[jra_.nref] = patjet.genJet()->pt();
-<
+            jra_.refeta[jra_.nref] = patjet.genJet()->eta();
-<
+            jra_.refphi[jra_.nref] = patjet.genJet()->phi();
-<
-<
+         }else{
-<
+            jra_.refpt[jra_.nref] = -99;
-<
+            jra_.refeta[jra_.nref] = -99;
-<
+            jra_.refphi[jra_.nref] = -99;
-<
+         }
-<
+      }
-<
-<
+      jra_.nref++;
-<
->
+     if(filterJets_ && !selectJet(j)) continue;
->
+     const reco::Jet& jet = (*jets)[j];
->
+     JRAV jv;
->
+     jv.jtpt = jet.pt();
->
+     jv.jteta = jet.eta();
->
+     jv.jtphi = jet.phi();
->
+     jv.jtcorpt = jet.pt();
->
+     jv.area = jet.jetArea();
->
+     jv.pu  = jet.pileup();
->
+     jv.index = j;
->
->
+     double pt = jet.pt();
->
+     double ktRho = getRho(jv.jteta,*ktRhos);
->
+     double akRho = getRho(jv.jteta,*akRhos);
->
->
+     jv.rho = ktRho;
->
+     if(correctJets_){
->
+        if(doFastJets_){
->
+           jv.pu = jet.jetArea()*ktRho;
->
+           pt -= jv.pu;
->
+        }
->
->
+        jetCorrector_->setJetEta(jet.eta());
->
+        jetCorrector_->setJetPt(pt);
->
+        //      jetCorrector_->setJetE(jet.energy());
->
+        vector<float> corrs = jetCorrector_->getSubCorrections();
->
+        jv.l2 = corrs[0];
->
+        jv.l3 = corrs[1];
->
+        jv.jtcorpt = pt*jv.l2*jv.l3;
->
->
+     }
->
+     if(usePat_){
->
+       const pat::Jet& patjet = (*patjets)[j];
->
->
+       jv.jtpt = patjet.correctedJet("raw").pt();
->
+       jv.jtcorpt = patjet.pt();
->
->
+       if(doMC_ && matchPatGen_ && patjet.genJet() != 0){
->
+         if(patjet.genJet()->pt() < genPtMin_) continue;
->
+         jv.refpt = patjet.genJet()->pt();
->
+         jv.refeta = patjet.genJet()->eta();
->
+         jv.refphi = patjet.genJet()->phi();
->
+       }else{
->
+         jv.refpt = -99;
->
+         jv.refeta = -99;
->
+         jv.refphi = -99;
->
+       }
->
+     }
->
+     jraV.push_back(jv);
->
+   }
->
->
+   if(sortJets_){
->
+     if(usePat_ || correctJets_) std::sort(jraV.begin(),jraV.end(),compareCorPt);
->
+     else std::sort(jraV.begin(),jraV.end(),comparePt);
-+
+   for(unsigned int i = 0; i < jraV.size(); ++i){
-+
+     JRAV& jv = jraV[i];
-+
+     const reco::Jet& jet = (*jets)[jv.index];
-+
-+
+     if(matchNew_){
-+
+       for(unsigned int im = 0; im < matchTags_.size(); ++im){
-+
+          edm::Handle<reco::JetView> matchedJets;
-+
+          iEvent.getByLabel(matchTags_[im],matchedJets);
-+
+         jraMatch_[im].jtcorpt[i] = -99;
-+
+         jraMatch_[im].jtpt[i] = -99;
-+
+         jraMatch_[im].jteta[i] = -99;
-+
+         jraMatch_[im].jtphi[i] = -99;
-+
+         jraMatch_[im].area[i] = -99;
-+
+         jraMatch_[im].l2[i] = -99;
-+
+         jraMatch_[im].l3[i] = -99;
-+
+         jraMatch_[im].pu[i] = -99;
-+
+         for(unsigned int m = 0 ; m < matchedJets->size(); ++m){
-+
+           const reco::Jet& match = (*matchedJets)[m];
-+
+           double dr = reco::deltaR(jet.eta(),jet.phi(),match.eta(),match.phi());
-+
+           if(dr < matchR_ && match.pt() > genPtMin_){
-+
+             jraMatch_[im].jtcorpt[i] = -99;
-+
+             jraMatch_[im].jtpt[i] = match.pt();
-+
+             jraMatch_[im].jteta[i] = match.eta();
-+
+             jraMatch_[im].jtphi[i] = match.phi();
-+
+             jraMatch_[im].area[i] = match.jetArea();
-+
-+
+             double  ktRho = getRho(jraMatch_[im].jteta[i],*ktRhos);
-+
+             double akRho  = getRho(jraMatch_[im].jteta[i],*akRhos);
-+
-+
+             jraMatch_[im].rho[i] = ktRho;
-+
-+
+             if((bool)correctMatchedJets_[im]){
-+
-+
+                double ktRho;
-+
+                double pt = jraMatch_[im].jtpt[i];
-+
-+
+                if((bool)doMatchedFastJets_[im]){
-+
+                   jraMatch_[im].pu[i] = jraMatch_[im].area[i]*ktRho;
-+
+                   pt -= jraMatch_[im].pu[i];
-+
+                }
-+
+                jetCorrector_->setJetEta(jraMatch_[im].jteta[i]);
-+
+                jetCorrector_->setJetPt(pt);
-+
+                //              jetCorrector_->setJetE(match.energy());
-+
+                vector<float> corrs = jetCorrector_->getSubCorrections();
-+
+                jraMatch_[im].l2[i] = corrs[0];
-+
+                // Should it be re-set for L3???
-+
+                jraMatch_[im].l3[i] = corrs[1];
-+
+                jraMatch_[im].jtcorpt[i] = pt*jraMatch_[im].l2[i]*jraMatch_[im].l3[i];
-+
-+
+             }
-+
-+
-+
+           }
-+
+         }
-+
+       }
-+
+     }
-+
-+
+     jra_.jtpt[i] = jv.jtpt;
-+
+     jra_.jteta[i] = jv.jteta;
-+
+     jra_.jtphi[i] = jv.jtphi;
-+
+     jra_.jtcorpt[i] = jv.jtcorpt;
-+
+     jra_.refpt[i] = jv.refpt;
-+
+     jra_.refeta[i] = jv.refeta;
-+
+     jra_.refphi[i] = jv.refphi;
-+
-+
+     jra_.area[i] = jv.area;
-+
+     jra_.pu[i] = jv.pu;
-+
+     jra_.rho[i] = jv.rho;
-+
+     jra_.l2[i] = jv.l2;
-+
+     jra_.l3[i] = jv.l3;
-+
-+
+   }
-+
+   jra_.nref = jraV.size();
-+
+   t->Fill();
-–
+// ------------ method called once each job just before starting event loop  ------------
-+
-+
+void
-<
+HiJetResponseAnalyzer::beginJob()
-<
+{
->
+HiJetResponseAnalyzer::beginJob(){
->
+  t= fs->make<TTree>("t","Jet Response Analyzer");
->
+  t->Branch("nref",&jra_.nref,"nref/I");
->
+  t->Branch("jtpt",jra_.jtpt,"jtpt[nref]/F");
->
+  t->Branch("jteta",jra_.jteta,"jteta[nref]/F");
->
+  t->Branch("jtphi",jra_.jtphi,"jtphi[nref]/F");
->
->
+  if(correctJets_){
->
+     t->Branch("jtcorpt",jra_.jtcorpt,"jtcorpt[nref]/F");
->
+     t->Branch("l2",jra_.l2,"l2[nref]/F");
->
+     t->Branch("l3",jra_.l3,"l3[nref]/F");
->
+  }
->
+  t->Branch("area",jra_.area,"area[nref]/F");
->
+  t->Branch("pu",jra_.pu,"pu[nref]/F");
->
+  t->Branch("rho",jra_.rho,"rho[nref]/F");
-<
+   t= fs->make<TTree>("t","Jet Response Analyzer");
-<
+   t->Branch("b",&jra_.b,"b/F");
-<
+   t->Branch("hf",&jra_.hf,"hf/F");
-<
+   t->Branch("nref",&jra_.nref,"nref/I");
-<
+   t->Branch("jtpt",jra_.jtpt,"jtpt[nref]/F");
-<
+   t->Branch("refpt",jra_.refpt,"refpt[nref]/F");
-<
+   t->Branch("jteta",jra_.jteta,"jteta[nref]/F");
->
+  t->Branch("refpt",jra_.refpt,"refpt[nref]/F");
->
+  t->Branch("refcorpt",jra_.refpt,"refcorpt[nref]/F");
+   t->Branch("refeta",jra_.refeta,"refeta[nref]/F");
-–
+   t->Branch("jtphi",jra_.jtphi,"jtphi[nref]/F");
+   t->Branch("refphi",jra_.refphi,"refphi[nref]/F");
+   t->Branch("weight",&jra_.weight,"weight/F");
-–
+   t->Branch("bin",&jra_.bin,"bin/I");
-+
+   jraMatch_.clear();
-+
+   for(unsigned int im = 0; im < matchTags_.size(); ++im){
-+
+     JRA jrm;
-+
+     jraMatch_.push_back(jrm);
-+
+   }
-+
+   for(unsigned int im = 0; im < matchTags_.size(); ++im){
-+
+     t->Branch(Form("jtpt%d",im),jraMatch_[im].jtpt,Form("jtpt%d[nref]/F",im));
-+
+     t->Branch(Form("jteta%d",im),jraMatch_[im].jteta,Form("jteta%d[nref]/F",im));
-+
+     t->Branch(Form("jtphi%d",im),jraMatch_[im].jtphi,Form("jtphi%d[nref]/F",im));
-+
-+
+     if((bool)correctMatchedJets_[im]){
-+
+        t->Branch(Form("jtcorpt%d",im),jraMatch_[im].jtcorpt,Form("jtcorpt%d[nref]/F",im));
-+
+        t->Branch(Form("l2_%d",im),jraMatch_[im].l2,Form("l2_%d[nref]/F",im));
-+
+        t->Branch(Form("l3_%d",im),jraMatch_[im].l3,Form("l3_%d[nref]/F",im));
-+
+     }
-+
-+
+     t->Branch(Form("area%d",im),jraMatch_[im].area,Form("area%d[nref]/F",im));
-+
+     t->Branch(Form("pu%d",im),jraMatch_[im].pu,Form("pu%d[nref]/F",im));
-+
+     t->Branch(Form("rho%d",im),jraMatch_[im].rho,Form("rho%d[nref]/F",im));
-+
+   }
-+
+// ------------ method called once each job just after ending the event loop  ------------

Diff Legend

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

Comparing UserCode/CmsHi/JetAnalysis/src/HiJetResponseAnalyzer.cc (file contents): Revision 1.2 by yilmaz, Fri Sep 10 12:54:47 2010 UTC vs. Revision 1.10 by yilmaz, Thu Jan 20 19:41:10 2011 UTC

Diff Legend

Comparing UserCode/CmsHi/JetAnalysis/src/HiJetResponseAnalyzer.cc (file contents):
Revision 1.2 by yilmaz, Fri Sep 10 12:54:47 2010 UTC vs.
Revision 1.10 by yilmaz, Thu Jan 20 19:41:10 2011 UTC