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.6 2010/10/18 16:13:37 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 "CommonTools/UtilAlgos/interface/TFileService.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "DataFormats/Math/interface/deltaR.h"

#include "TTree.h"

using namespace std;

static const int MAXJETS = 500;

struct etdr{
   double et;
   double dr;
};

struct JRA{

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


struct JRAV{

  float jtpt;
  float jtcorpt;
  float refpt;
  float refcorpt;
  float jteta;
  float refeta;
  float jtphi;
  float refphi;
 
};

//
// 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_;

  double matchR_;  
   double genPtMin_;
   double ptMin_;
   double emfMin_;
   double n90Min_;
   double n90hitMin_;

   edm::InputTag jetTag_;
  edm::InputTag matchTag_;

   JRA jra_;
   TTree* t;

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

   edm::Handle<reco::JetView> jets;
   edm::Handle<pat::JetCollection> patjets;
  edm::Handle<reco::JetView> matchedJets;

   edm::Service<TFileService> fs;

};

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)

{
   //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);

   jetTag_ = iConfig.getUntrackedParameter<edm::InputTag>("src",edm::InputTag("selectedPatJets"));
   matchTag_ = iConfig.getUntrackedParameter<edm::InputTag>("match",edm::InputTag("selectedPatJets"));
}


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);
   if(matchNew_)iEvent.getByLabel(matchTag_,matchedJets);

   std::vector<JRAV> jraV;

   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();

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

     if(matchNew_){
       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_){
           jv.refcorpt = -99;
           jv.refpt = match.pt();
           jv.refeta = match.eta();
           jv.refphi = match.phi();
         }
       }
     }
     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];
     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_.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("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("jtcorpt",jra_.jtcorpt,"jtcorpt[nref]/F");
   t->Branch("refpt",jra_.refpt,"refpt[nref]/F");
   t->Branch("refcorpt",jra_.refpt,"refcorpt[nref]/F");
   t->Branch("jteta",jra_.jteta,"jteta[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");


}

// ------------ 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.7
Committed:	Thu Oct 21 11:26:40 2010 UTC (14 years, 6 months ago) by yilmaz
Content type:	text/plain
Branch:	MAIN
Changes since 1.6:	+8 -9 lines
Log Message:	Bug Fixx
#	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.7	// $Id: HiJetResponseAnalyzer.cc,v 1.6 2010/10/18 16:13:37 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.2	#include "CommonTools/UtilAlgos/interface/TFileService.h"
44			#include "FWCore/ServiceRegistry/interface/Service.h"
45	yilmaz	1.6	#include "DataFormats/Math/interface/deltaR.h"
46	yilmaz	1.2
47	yilmaz	1.1	#include "TTree.h"
48
49			using namespace std;
50
51			static const int MAXJETS = 500;
52
53			struct etdr{
54			double et;
55			double dr;
56			};
57
58			struct JRA{
59
60			int nref;
61			int bin;
62			float b;
63			float hf;
64			float jtpt[MAXJETS];
65	yilmaz	1.4	float jtcorpt[MAXJETS];
66	yilmaz	1.1	float refpt[MAXJETS];
67			float jteta[MAXJETS];
68			float refeta[MAXJETS];
69			float jtphi[MAXJETS];
70			float refphi[MAXJETS];
71
72			float weight;
73			};
74
75	yilmaz	1.6
76			struct JRAV{
77
78			float jtpt;
79			float jtcorpt;
80			float refpt;
81			float refcorpt;
82			float jteta;
83			float refeta;
84			float jtphi;
85			float refphi;
86
87			};
88
89	yilmaz	1.1	//
90			// class declaration
91			//
92	yilmaz	1.6	bool compareCorPt(JRAV a, JRAV b) {return a.jtcorpt > b.jtcorpt;}
93			bool comparePt(JRAV a, JRAV b) {return a.jtpt > b.jtpt;}
94	yilmaz	1.1
95			class HiJetResponseAnalyzer : public edm::EDAnalyzer {
96			public:
97			explicit HiJetResponseAnalyzer(const edm::ParameterSet&);
98			~HiJetResponseAnalyzer();
99
100
101			private:
102			virtual void beginJob() ;
103			virtual void analyze(const edm::Event&, const edm::EventSetup&);
104			virtual void endJob() ;
105	yilmaz	1.5	bool selectJet(int i);
106	yilmaz	1.1	// ----------member data ---------------------------
107
108	yilmaz	1.6	bool usePat_;
109			bool doMC_;
110			bool filterJets_;
111			bool diJetsOnly_;
112			bool matchDiJets_;
113			bool matchPatGen_;
114			bool matchNew_;
115			bool sortJets_;
116			bool correctJets_;
117	yilmaz	1.5
118	yilmaz	1.6	double matchR_;
119	yilmaz	1.5	double genPtMin_;
120			double ptMin_;
121			double emfMin_;
122			double n90Min_;
123			double n90hitMin_;
124	yilmaz	1.2
125			edm::InputTag jetTag_;
126	yilmaz	1.6	edm::InputTag matchTag_;
127	yilmaz	1.2
128	yilmaz	1.1	JRA jra_;
129			TTree* t;
130
131			edm::Handle<edm::GenHIEvent> mc;
132			edm::Handle<reco::Centrality> cent;
133
134			edm::Handle<reco::JetView> jets;
135	yilmaz	1.5	edm::Handle<pat::JetCollection> patjets;
136	yilmaz	1.6	edm::Handle<reco::JetView> matchedJets;
137	yilmaz	1.1
138	yilmaz	1.2	edm::Service<TFileService> fs;
139	yilmaz	1.1
140			};
141
142	yilmaz	1.5	bool HiJetResponseAnalyzer::selectJet(int i){
143			const reco::Jet& jet = (*jets)[i];
144			if(usePat_){
145			const pat::Jet& patjet = (*patjets)[i];
146			if(patjet.emEnergyFraction() <= emfMin_) return false;
147			if(patjet.jetID().n90Hits <= n90hitMin_) return false;
148			if(doMC_){
149
150			}
151
152			}
153
154			return true;
155			}
156
157
158	yilmaz	1.1	//
159			// constants, enums and typedefs
160			//
161
162			//
163			// static data member definitions
164			//
165
166			//
167			// constructors and destructor
168			//
169			HiJetResponseAnalyzer::HiJetResponseAnalyzer(const edm::ParameterSet& iConfig)
170
171			{
172			//now do what ever initialization is needed
173	yilmaz	1.6	matchR_ = iConfig.getUntrackedParameter<double>("matchR",0.25);
174	yilmaz	1.2
175	yilmaz	1.5	ptMin_ = iConfig.getUntrackedParameter<double>("ptMin",0);
176	yilmaz	1.6	genPtMin_ = iConfig.getUntrackedParameter<double>("genPtMin",20);
177			emfMin_ = iConfig.getUntrackedParameter<double>("emfMin",0.01);
178			n90Min_ = iConfig.getUntrackedParameter<double>("n90Min",1);
179			n90hitMin_ = iConfig.getUntrackedParameter<double>("n90hitMin",1);
180	yilmaz	1.5
181			filterJets_ = iConfig.getUntrackedParameter<bool>("filterJets",true);
182	yilmaz	1.6	diJetsOnly_ = iConfig.getUntrackedParameter<bool>("diJetsOnly",false);
183			matchDiJets_ = iConfig.getUntrackedParameter<bool>("matchDiJets",false);
184			matchPatGen_ = iConfig.getUntrackedParameter<bool>("matchPatGen",false);
185
186			matchNew_ = iConfig.getUntrackedParameter<bool>("matchNew",false);
187
188	yilmaz	1.4	usePat_ = iConfig.getUntrackedParameter<bool>("usePat",true);
189			doMC_ = iConfig.getUntrackedParameter<bool>("doMC",true);
190	yilmaz	1.6
191			sortJets_ = iConfig.getUntrackedParameter<bool>("sortJets",true);
192			correctJets_ = iConfig.getUntrackedParameter<bool>("correctJets",false);
193
194	yilmaz	1.4	jetTag_ = iConfig.getUntrackedParameter<edm::InputTag>("src",edm::InputTag("selectedPatJets"));
195	yilmaz	1.6	matchTag_ = iConfig.getUntrackedParameter<edm::InputTag>("match",edm::InputTag("selectedPatJets"));
196	yilmaz	1.1	}
197
198
199			HiJetResponseAnalyzer::~HiJetResponseAnalyzer()
200			{
201
202			// do anything here that needs to be done at desctruction time
203			// (e.g. close files, deallocate resources etc.)
204
205			}
206
207
208			//
209			// member functions
210			//
211
212			// ------------ method called to for each event ------------
213			void
214			HiJetResponseAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
215			{
216			using namespace edm;
217
218	yilmaz	1.2	iEvent.getByLabel(jetTag_,jets);
219	yilmaz	1.4	if(usePat_)iEvent.getByLabel(jetTag_,patjets);
220	yilmaz	1.6	if(matchNew_)iEvent.getByLabel(matchTag_,matchedJets);
221
222			std::vector<JRAV> jraV;
223	yilmaz	1.4
224	yilmaz	1.1	for(unsigned int j = 0 ; j < jets->size(); ++j){
225	yilmaz	1.6	if(filterJets_ && !selectJet(j)) continue;
226			const reco::Jet& jet = (*jets)[j];
227			JRAV jv;
228			jv.jtpt = jet.pt();
229			jv.jteta = jet.eta();
230			jv.jtphi = jet.phi();
231			jv.jtcorpt = jet.pt();
232
233			if(usePat_){
234			const pat::Jet& patjet = (*patjets)[j];
235
236			jv.jtpt = patjet.correctedJet("raw").pt();
237			jv.jtcorpt = patjet.pt();
238
239			if(doMC_ && matchPatGen_ && patjet.genJet() != 0){
240			if(patjet.genJet()->pt() < genPtMin_) continue;
241			jv.refpt = patjet.genJet()->pt();
242			jv.refeta = patjet.genJet()->eta();
243			jv.refphi = patjet.genJet()->phi();
244
245			}else{
246			jv.refpt = -99;
247			jv.refeta = -99;
248			jv.refphi = -99;
249			}
250			}
251
252			if(matchNew_){
253			for(unsigned int m = 0 ; m < matchedJets->size(); ++m){
254			const reco::Jet& match = (*matchedJets)[m];
255			double dr = reco::deltaR(jet.eta(),jet.phi(),match.eta(),match.phi());
256			if(dr < matchR_){
257			jv.refcorpt = -99;
258			jv.refpt = match.pt();
259			jv.refeta = match.eta();
260			jv.refphi = match.phi();
261			}
262			}
263			}
264			jraV.push_back(jv);
265			}
266
267			if(sortJets_){
268			if(usePat_ \|\| correctJets_) std::sort(jraV.begin(),jraV.end(),compareCorPt);
269			else std::sort(jraV.begin(),jraV.end(),comparePt);
270			}
271	yilmaz	1.4
272	yilmaz	1.6	for(unsigned int i = 0; i < jraV.size(); ++i){
273			JRAV& jv = jraV[i];
274	yilmaz	1.7	jra_.jtpt[i] = jv.jtpt;
275			jra_.jteta[i] = jv.jteta;
276			jra_.jtphi[i] = jv.jtphi;
277			jra_.jtcorpt[i] = jv.jtcorpt;
278			jra_.refpt[i] = jv.refpt;
279			jra_.refeta[i] = jv.refeta;
280			jra_.refphi[i] = jv.refphi;
281	yilmaz	1.1	}
282	yilmaz	1.6	jra_.nref = jraV.size();
283	yilmaz	1.1
284			t->Fill();
285
286			}
287
288
289			// ------------ method called once each job just before starting event loop ------------
290			void
291			HiJetResponseAnalyzer::beginJob()
292			{
293	yilmaz	1.2
294			t= fs->make<TTree>("t","Jet Response Analyzer");
295			t->Branch("b",&jra_.b,"b/F");
296			t->Branch("hf",&jra_.hf,"hf/F");
297			t->Branch("nref",&jra_.nref,"nref/I");
298			t->Branch("jtpt",jra_.jtpt,"jtpt[nref]/F");
299	yilmaz	1.4	t->Branch("jtcorpt",jra_.jtcorpt,"jtcorpt[nref]/F");
300	yilmaz	1.2	t->Branch("refpt",jra_.refpt,"refpt[nref]/F");
301	yilmaz	1.6	t->Branch("refcorpt",jra_.refpt,"refcorpt[nref]/F");
302	yilmaz	1.2	t->Branch("jteta",jra_.jteta,"jteta[nref]/F");
303			t->Branch("refeta",jra_.refeta,"refeta[nref]/F");
304			t->Branch("jtphi",jra_.jtphi,"jtphi[nref]/F");
305			t->Branch("refphi",jra_.refphi,"refphi[nref]/F");
306			t->Branch("weight",&jra_.weight,"weight/F");
307			t->Branch("bin",&jra_.bin,"bin/I");
308
309
310
311	yilmaz	1.1	}
312
313			// ------------ method called once each job just after ending the event loop ------------
314			void
315			HiJetResponseAnalyzer::endJob() {
316			}
317
318			//define this as a plug-in
319			DEFINE_FWK_MODULE(HiJetResponseAnalyzer);