JetAnalysis/src/HiJetResponseAnalyzer.cc

// -*- C++ -*-
//
// Package:    HiJetResponseAnalyzer
// Class:      HiJetResponseAnalyzer
// 
/**\class HiJetResponseAnalyzer HiJetResponseAnalyzer.cc CmsHi/HiJetResponseAnalyzer/src/HiJetResponseAnalyzer.cc

 Description: [one line class summary]

 Implementation:
     [Notes on implementation]
*/
//
// Original Author:  Yetkin Yilmaz
//         Created:  Thu Sep  9 10:38:59 EDT 2010
// $Id: HiJetResponseAnalyzer.cc,v 1.9 2010/10/24 14:39:32 yilmaz Exp $
//
//


// system include files
#include <memory>
#include <iostream>

// user include files
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/EDAnalyzer.h"

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"

#include "FWCore/ParameterSet/interface/ParameterSet.h"

#include "FWCore/Utilities/interface/InputTag.h"
#include "DataFormats/HeavyIonEvent/interface/CentralityBins.h"
#include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
#include "DataFormats/HeavyIonEvent/interface/Centrality.h"
#include "DataFormats/PatCandidates/interface/Jet.h"
#include "SimDataFormats/HiGenData/interface/GenHIEvent.h"
#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;

struct etdr{
   double et;
   double dr;
};

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:
      explicit HiJetResponseAnalyzer(const edm::ParameterSet&);
      ~HiJetResponseAnalyzer();


   private:
      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_;
  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_;

   TTree* t;

   edm::Handle<edm::GenHIEvent> mc;
   edm::Handle<reco::Centrality> cent;

   edm::Handle<reco::JetView> 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
//

//
// static data member definitions
//

//
// constructors and destructor
//
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
  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);

   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];
   }


}


HiJetResponseAnalyzer::~HiJetResponseAnalyzer()
{
 
   // do anything here that needs to be done at desctruction time
   // (e.g. close files, deallocate resources etc.)

}


//
// member functions
//

// ------------ method called to for each event  ------------
void
HiJetResponseAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
{
   using namespace edm;

   iEvent.getByLabel(jetTag_,jets);
   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);
   }

   for(unsigned int j = 0 ; j < jets->size(); ++j){
     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(){
  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->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("refphi",jra_.refphi,"refphi[nref]/F");
   t->Branch("weight",&jra_.weight,"weight/F");

   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  ------------
void 
HiJetResponseAnalyzer::endJob() {
}

//define this as a plug-in
DEFINE_FWK_MODULE(HiJetResponseAnalyzer);
Revision:	1.10
Committed:	Thu Jan 20 19:41:10 2011 UTC (14 years, 3 months ago) by yilmaz
Content type:	text/plain
Branch:	MAIN
Changes since 1.9:	+170 -8 lines
Log Message:	more efficient
#	User	Rev	Content
1	yilmaz	1.1	// -- C++ --
2			//
3			// Package: HiJetResponseAnalyzer
4			// Class: HiJetResponseAnalyzer
5			//
6			/**\class HiJetResponseAnalyzer HiJetResponseAnalyzer.cc CmsHi/HiJetResponseAnalyzer/src/HiJetResponseAnalyzer.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.10	// $Id: HiJetResponseAnalyzer.cc,v 1.9 2010/10/24 14:39:32 yilmaz Exp $
17	yilmaz	1.1	//
18			//
19
20
21			// system include files
22			#include <memory>
23	yilmaz	1.2	#include <iostream>
24	yilmaz	1.1
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	yilmaz	1.10	#include "CondFormats/JetMETObjects/interface/FactorizedJetCorrector.h"
44			#include "JetMETCorrections/Objects/interface/JetCorrectionsRecord.h"
45			#include "CondFormats/JetMETObjects/interface/JetCorrectionUncertainty.h"
46
47	yilmaz	1.2	#include "CommonTools/UtilAlgos/interface/TFileService.h"
48			#include "FWCore/ServiceRegistry/interface/Service.h"
49	yilmaz	1.6	#include "DataFormats/Math/interface/deltaR.h"
50	yilmaz	1.10	#include "FWCore/ParameterSet/interface/FileInPath.h"
51			#include "FWCore/Framework/interface/ESWatcher.h"
52	yilmaz	1.2
53	yilmaz	1.10	#include "CmsHi/JetAnalysis/interface/RhoGetter.h"
54	yilmaz	1.1	#include "TTree.h"
55
56			using namespace std;
57	yilmaz	1.10	using namespace edm;
58	yilmaz	1.1
59			static const int MAXJETS = 500;
60
61			struct etdr{
62			double et;
63			double dr;
64			};
65
66	yilmaz	1.10	class JRA{
67	yilmaz	1.1
68	yilmaz	1.10	public:
69	yilmaz	1.1	int nref;
70			int bin;
71			float b;
72			float hf;
73			float jtpt[MAXJETS];
74	yilmaz	1.4	float jtcorpt[MAXJETS];
75	yilmaz	1.1	float refpt[MAXJETS];
76			float jteta[MAXJETS];
77			float refeta[MAXJETS];
78			float jtphi[MAXJETS];
79			float refphi[MAXJETS];
80	yilmaz	1.10	float l2[MAXJETS];
81			float l3[MAXJETS];
82			float area[MAXJETS];
83			float pu[MAXJETS];
84			float rho[MAXJETS];
85	yilmaz	1.1
86			float weight;
87			};
88
89	yilmaz	1.6	struct JRAV{
90	yilmaz	1.8	int index;
91	yilmaz	1.6	float jtpt;
92			float jtcorpt;
93			float refpt;
94			float refcorpt;
95			float jteta;
96			float refeta;
97			float jtphi;
98			float refphi;
99	yilmaz	1.10	float l2;
100			float l3;
101			float area;
102			float pu;
103			float rho;
104	yilmaz	1.6
105			};
106
107	yilmaz	1.1	//
108			// class declaration
109			//
110	yilmaz	1.6	bool compareCorPt(JRAV a, JRAV b) {return a.jtcorpt > b.jtcorpt;}
111			bool comparePt(JRAV a, JRAV b) {return a.jtpt > b.jtpt;}
112	yilmaz	1.1
113			class HiJetResponseAnalyzer : public edm::EDAnalyzer {
114			public:
115			explicit HiJetResponseAnalyzer(const edm::ParameterSet&);
116			~HiJetResponseAnalyzer();
117
118
119			private:
120			virtual void beginJob() ;
121			virtual void analyze(const edm::Event&, const edm::EventSetup&);
122			virtual void endJob() ;
123	yilmaz	1.5	bool selectJet(int i);
124	yilmaz	1.1	// ----------member data ---------------------------
125
126	yilmaz	1.6	bool usePat_;
127			bool doMC_;
128			bool filterJets_;
129			bool diJetsOnly_;
130			bool matchDiJets_;
131			bool matchPatGen_;
132			bool matchNew_;
133			bool sortJets_;
134			bool correctJets_;
135	yilmaz	1.10	bool getFastJets_;
136	yilmaz	1.5
137	yilmaz	1.6	double matchR_;
138	yilmaz	1.5	double genPtMin_;
139			double ptMin_;
140			double emfMin_;
141			double n90Min_;
142			double n90hitMin_;
143	yilmaz	1.2
144	yilmaz	1.8	edm::InputTag jetTag_;
145	yilmaz	1.6	edm::InputTag matchTag_;
146	yilmaz	1.8	std::vector<edm::InputTag> matchTags_;
147
148			JRA jra_;
149			std::vector<JRA> jraMatch_;
150	yilmaz	1.2
151	yilmaz	1.1	TTree* t;
152
153			edm::Handle<edm::GenHIEvent> mc;
154			edm::Handle<reco::Centrality> cent;
155
156			edm::Handle<reco::JetView> jets;
157	yilmaz	1.5	edm::Handle<pat::JetCollection> patjets;
158	yilmaz	1.10
159			FactorizedJetCorrector* jetCorrector_;
160			edm::ESWatcher<JetCorrectionsRecord> watchJetCorrectionsRecord_;
161
162			std::string tags_;
163			std::string levels_;
164			std::string algo_;
165	yilmaz	1.1
166	yilmaz	1.2	edm::Service<TFileService> fs;
167	yilmaz	1.1
168	yilmaz	1.10	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	yilmaz	1.1	};
180
181	yilmaz	1.5	bool HiJetResponseAnalyzer::selectJet(int i){
182			const reco::Jet& jet = (*jets)[i];
183			if(usePat_){
184			const pat::Jet& patjet = (*patjets)[i];
185			if(patjet.emEnergyFraction() <= emfMin_) return false;
186			if(patjet.jetID().n90Hits <= n90hitMin_) return false;
187			if(doMC_){
188
189			}
190
191			}
192
193			return true;
194			}
195
196
197	yilmaz	1.1	//
198			// constants, enums and typedefs
199			//
200
201			//
202			// static data member definitions
203			//
204
205			//
206			// constructors and destructor
207			//
208			HiJetResponseAnalyzer::HiJetResponseAnalyzer(const edm::ParameterSet& iConfig)
209
210			{
211	yilmaz	1.10
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	yilmaz	1.1	//now do what ever initialization is needed
228	yilmaz	1.6	matchR_ = iConfig.getUntrackedParameter<double>("matchR",0.25);
229	yilmaz	1.2
230	yilmaz	1.5	ptMin_ = iConfig.getUntrackedParameter<double>("ptMin",0);
231	yilmaz	1.6	genPtMin_ = iConfig.getUntrackedParameter<double>("genPtMin",20);
232			emfMin_ = iConfig.getUntrackedParameter<double>("emfMin",0.01);
233			n90Min_ = iConfig.getUntrackedParameter<double>("n90Min",1);
234			n90hitMin_ = iConfig.getUntrackedParameter<double>("n90hitMin",1);
235	yilmaz	1.5
236			filterJets_ = iConfig.getUntrackedParameter<bool>("filterJets",true);
237	yilmaz	1.6	diJetsOnly_ = iConfig.getUntrackedParameter<bool>("diJetsOnly",false);
238			matchDiJets_ = iConfig.getUntrackedParameter<bool>("matchDiJets",false);
239			matchPatGen_ = iConfig.getUntrackedParameter<bool>("matchPatGen",false);
240
241			matchNew_ = iConfig.getUntrackedParameter<bool>("matchNew",false);
242
243	yilmaz	1.4	usePat_ = iConfig.getUntrackedParameter<bool>("usePat",true);
244			doMC_ = iConfig.getUntrackedParameter<bool>("doMC",true);
245	yilmaz	1.6
246			sortJets_ = iConfig.getUntrackedParameter<bool>("sortJets",true);
247			correctJets_ = iConfig.getUntrackedParameter<bool>("correctJets",false);
248
249	yilmaz	1.10	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	yilmaz	1.4	jetTag_ = iConfig.getUntrackedParameter<edm::InputTag>("src",edm::InputTag("selectedPatJets"));
253	yilmaz	1.6	matchTag_ = iConfig.getUntrackedParameter<edm::InputTag>("match",edm::InputTag("selectedPatJets"));
254	yilmaz	1.8	matchTags_ = iConfig.getUntrackedParameter<std::vector<edm::InputTag> >("matches",std::vector<edm::InputTag>(0));
255	yilmaz	1.10
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	yilmaz	1.1	}
267
268
269			HiJetResponseAnalyzer::~HiJetResponseAnalyzer()
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			HiJetResponseAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
285			{
286			using namespace edm;
287
288	yilmaz	1.2	iEvent.getByLabel(jetTag_,jets);
289	yilmaz	1.4	if(usePat_)iEvent.getByLabel(jetTag_,patjets);
290	yilmaz	1.6	std::vector<JRAV> jraV;
291	yilmaz	1.4
292	yilmaz	1.10	if(getFastJets_){
293			iEvent.getByLabel(edm::InputTag(ktSrc_.label(),"rhos"),ktRhos);
294			iEvent.getByLabel(edm::InputTag(akSrc_.label(),"rhos"),akRhos);
295			}
296
297	yilmaz	1.1	for(unsigned int j = 0 ; j < jets->size(); ++j){
298	yilmaz	1.6	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.jtcorpt = jet.pt();
305	yilmaz	1.10	jv.area = jet.jetArea();
306			jv.pu = jet.pileup();
307	yilmaz	1.8	jv.index = j;
308	yilmaz	1.10
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 = ptjv.l2jv.l3;
327
328			}
329	yilmaz	1.6	if(usePat_){
330			const pat::Jet& patjet = (*patjets)[j];
331
332			jv.jtpt = patjet.correctedJet("raw").pt();
333			jv.jtcorpt = patjet.pt();
334
335			if(doMC_ && matchPatGen_ && patjet.genJet() != 0){
336			if(patjet.genJet()->pt() < genPtMin_) continue;
337			jv.refpt = patjet.genJet()->pt();
338			jv.refeta = patjet.genJet()->eta();
339			jv.refphi = patjet.genJet()->phi();
340			}else{
341			jv.refpt = -99;
342			jv.refeta = -99;
343			jv.refphi = -99;
344			}
345			}
346			jraV.push_back(jv);
347			}
348
349			if(sortJets_){
350			if(usePat_ \|\| correctJets_) std::sort(jraV.begin(),jraV.end(),compareCorPt);
351			else std::sort(jraV.begin(),jraV.end(),comparePt);
352			}
353	yilmaz	1.4
354	yilmaz	1.6	for(unsigned int i = 0; i < jraV.size(); ++i){
355			JRAV& jv = jraV[i];
356	yilmaz	1.8	const reco::Jet& jet = (*jets)[jv.index];
357
358			if(matchNew_){
359			for(unsigned int im = 0; im < matchTags_.size(); ++im){
360	yilmaz	1.10	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	yilmaz	1.8	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	yilmaz	1.10	if(dr < matchR_ && match.pt() > genPtMin_){
374	yilmaz	1.8	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	yilmaz	1.10	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] = ptjraMatch_[im].l2[i]jraMatch_[im].l3[i];
402
403			}
404
405
406	yilmaz	1.8	}
407			}
408			}
409			}
410
411	yilmaz	1.7	jra_.jtpt[i] = jv.jtpt;
412			jra_.jteta[i] = jv.jteta;
413			jra_.jtphi[i] = jv.jtphi;
414			jra_.jtcorpt[i] = jv.jtcorpt;
415			jra_.refpt[i] = jv.refpt;
416			jra_.refeta[i] = jv.refeta;
417			jra_.refphi[i] = jv.refphi;
418	yilmaz	1.10
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	yilmaz	1.1	}
426	yilmaz	1.6	jra_.nref = jraV.size();
427	yilmaz	1.8
428	yilmaz	1.1	t->Fill();
429
430			}
431
432	yilmaz	1.8	// ------------ method called once each job just before starting event loop ------------
433
434	yilmaz	1.1
435			void
436	yilmaz	1.8	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");
440			t->Branch("jteta",jra_.jteta,"jteta[nref]/F");
441			t->Branch("jtphi",jra_.jtphi,"jtphi[nref]/F");
442	yilmaz	1.10
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	yilmaz	1.8	t->Branch("refpt",jra_.refpt,"refpt[nref]/F");
453			t->Branch("refcorpt",jra_.refpt,"refcorpt[nref]/F");
454	yilmaz	1.2	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	yilmaz	1.10
458			jraMatch_.clear();
459	yilmaz	1.8	for(unsigned int im = 0; im < matchTags_.size(); ++im){
460			JRA jrm;
461			jraMatch_.push_back(jrm);
462	yilmaz	1.10	}
463
464			for(unsigned int im = 0; im < matchTags_.size(); ++im){
465	yilmaz	1.8	t->Branch(Form("jtpt%d",im),jraMatch_[im].jtpt,Form("jtpt%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	yilmaz	1.10
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	yilmaz	1.8	}
480
481	yilmaz	1.1	}
482
483			// ------------ method called once each job just after ending the event loop ------------
484			void
485			HiJetResponseAnalyzer::endJob() {
486			}
487
488			//define this as a plug-in
489			DEFINE_FWK_MODULE(HiJetResponseAnalyzer);