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.3 2010/09/10 14:47:48 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 "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;
};

//
// class declaration
//

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

      // ----------member data ---------------------------

   bool usePat_;
   bool doMC_;

   double genJetPtMin_;


   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::Service<TFileService> fs;

};

//
// constants, enums and typedefs
//

//
// static data member definitions
//

//
// constructors and destructor
//
HiJetResponseAnalyzer::HiJetResponseAnalyzer(const edm::ParameterSet& iConfig)

{
   //now do what ever initialization is needed
   genJetPtMin_ = 0;

   usePat_ = iConfig.getUntrackedParameter<bool>("usePat",true);
   doMC_ = iConfig.getUntrackedParameter<bool>("doMC",true);

   jetTag_ = iConfig.getUntrackedParameter<edm::InputTag>("src",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);

   edm::Handle<pat::JetCollection> patjets;
   if(usePat_)iEvent.getByLabel(jetTag_,patjets);

   jra_.nref = 0;
   for(unsigned int j = 0 ; j < jets->size(); ++j){
      
      //      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();
      jra_.jtcorpt[jra_.nref] = jet.pt();

      if(usePat_){
         const pat::Jet& patjet = (*patjets)[j];

         jra_.jtpt[jra_.nref] = patjet.correctedJet("raw").pt();
         jra_.jtcorpt[jra_.nref] = patjet.pt();
         
         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++;
      
   }

   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("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.4
Committed:	Tue Sep 28 16:02:52 2010 UTC (14 years, 7 months ago) by yilmaz
Content type:	text/plain
Branch:	MAIN
Changes since 1.3:	+15 -6 lines
Log Message:	Full jra tree producable
#	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.4	// $Id: HiJetResponseAnalyzer.cc,v 1.3 2010/09/10 14:47:48 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
46
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			//
76			// class declaration
77			//
78
79			class HiJetResponseAnalyzer : public edm::EDAnalyzer {
80			public:
81			explicit HiJetResponseAnalyzer(const edm::ParameterSet&);
82			~HiJetResponseAnalyzer();
83
84
85			private:
86			virtual void beginJob() ;
87			virtual void analyze(const edm::Event&, const edm::EventSetup&);
88			virtual void endJob() ;
89
90			// ----------member data ---------------------------
91
92	yilmaz	1.2	bool usePat_;
93			bool doMC_;
94
95	yilmaz	1.1	double genJetPtMin_;
96
97	yilmaz	1.2
98			edm::InputTag jetTag_;
99
100	yilmaz	1.1	JRA jra_;
101			TTree* t;
102
103			edm::Handle<edm::GenHIEvent> mc;
104			edm::Handle<reco::Centrality> cent;
105
106			edm::Handle<reco::JetView> jets;
107	yilmaz	1.2	// edm::Handle<pat::JetCollection> jets;
108	yilmaz	1.1
109
110	yilmaz	1.2	edm::Service<TFileService> fs;
111	yilmaz	1.1
112			};
113
114			//
115			// constants, enums and typedefs
116			//
117
118			//
119			// static data member definitions
120			//
121
122			//
123			// constructors and destructor
124			//
125			HiJetResponseAnalyzer::HiJetResponseAnalyzer(const edm::ParameterSet& iConfig)
126
127			{
128			//now do what ever initialization is needed
129			genJetPtMin_ = 0;
130	yilmaz	1.2
131	yilmaz	1.4	usePat_ = iConfig.getUntrackedParameter<bool>("usePat",true);
132			doMC_ = iConfig.getUntrackedParameter<bool>("doMC",true);
133	yilmaz	1.2
134	yilmaz	1.4	jetTag_ = iConfig.getUntrackedParameter<edm::InputTag>("src",edm::InputTag("selectedPatJets"));
135	yilmaz	1.2
136	yilmaz	1.1	}
137
138
139			HiJetResponseAnalyzer::~HiJetResponseAnalyzer()
140			{
141
142			// do anything here that needs to be done at desctruction time
143			// (e.g. close files, deallocate resources etc.)
144
145			}
146
147
148			//
149			// member functions
150			//
151
152			// ------------ method called to for each event ------------
153			void
154			HiJetResponseAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
155			{
156			using namespace edm;
157
158	yilmaz	1.2	iEvent.getByLabel(jetTag_,jets);
159	yilmaz	1.1
160	yilmaz	1.4	edm::Handle<pat::JetCollection> patjets;
161			if(usePat_)iEvent.getByLabel(jetTag_,patjets);
162
163	yilmaz	1.2	jra_.nref = 0;
164	yilmaz	1.1	for(unsigned int j = 0 ; j < jets->size(); ++j){
165	yilmaz	1.2
166			// const pat::Jet& jet = (*jets)[j];
167			const reco::Jet jet = (*jets)[j];
168			jra_.jtpt[jra_.nref] = jet.pt();
169			jra_.jteta[jra_.nref] = jet.eta();
170			jra_.jtphi[jra_.nref] = jet.phi();
171	yilmaz	1.4	jra_.jtcorpt[jra_.nref] = jet.pt();
172
173	yilmaz	1.2	if(usePat_){
174	yilmaz	1.4	const pat::Jet& patjet = (*patjets)[j];
175
176			jra_.jtpt[jra_.nref] = patjet.correctedJet("raw").pt();
177			jra_.jtcorpt[jra_.nref] = patjet.pt();
178	yilmaz	1.2
179			if(doMC_ && patjet.genJet() != 0){
180			// if(jet.genJet()->pt() < genJetPtMin_) continue;
181			jra_.refpt[jra_.nref] = patjet.genJet()->pt();
182			jra_.refeta[jra_.nref] = patjet.genJet()->eta();
183			jra_.refphi[jra_.nref] = patjet.genJet()->phi();
184
185			}else{
186			jra_.refpt[jra_.nref] = -99;
187			jra_.refeta[jra_.nref] = -99;
188			jra_.refphi[jra_.nref] = -99;
189			}
190	yilmaz	1.1	}
191	yilmaz	1.2
192			jra_.nref++;
193
194	yilmaz	1.1	}
195
196			t->Fill();
197
198			}
199
200
201			// ------------ method called once each job just before starting event loop ------------
202			void
203			HiJetResponseAnalyzer::beginJob()
204			{
205	yilmaz	1.2
206			t= fs->make<TTree>("t","Jet Response Analyzer");
207			t->Branch("b",&jra_.b,"b/F");
208			t->Branch("hf",&jra_.hf,"hf/F");
209			t->Branch("nref",&jra_.nref,"nref/I");
210			t->Branch("jtpt",jra_.jtpt,"jtpt[nref]/F");
211	yilmaz	1.4	t->Branch("jtcorpt",jra_.jtcorpt,"jtcorpt[nref]/F");
212	yilmaz	1.2	t->Branch("refpt",jra_.refpt,"refpt[nref]/F");
213			t->Branch("jteta",jra_.jteta,"jteta[nref]/F");
214			t->Branch("refeta",jra_.refeta,"refeta[nref]/F");
215			t->Branch("jtphi",jra_.jtphi,"jtphi[nref]/F");
216			t->Branch("refphi",jra_.refphi,"refphi[nref]/F");
217			t->Branch("weight",&jra_.weight,"weight/F");
218			t->Branch("bin",&jra_.bin,"bin/I");
219
220
221
222	yilmaz	1.1	}
223
224			// ------------ method called once each job just after ending the event loop ------------
225			void
226			HiJetResponseAnalyzer::endJob() {
227			}
228
229			//define this as a plug-in
230			DEFINE_FWK_MODULE(HiJetResponseAnalyzer);