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.1 by yilmaz, Thu Sep 9 15:05:01 2010 UTC vs.
Revision 1.7 by yilmaz, Thu Oct 21 11:26:40 2010 UTC

# Line 20 | Line 20
20  
21   // system include files
22   #include <memory>
23 + #include <iostream>
24  
25   // user include files
26   #include "FWCore/Framework/interface/Frameworkfwd.h"
# Line 39 | Line 40
40   #include "DataFormats/JetReco/interface/CaloJetCollection.h"
41   #include "DataFormats/JetReco/interface/GenJetCollection.h"
42  
43 + #include "CommonTools/UtilAlgos/interface/TFileService.h"
44 + #include "FWCore/ServiceRegistry/interface/Service.h"
45 + #include "DataFormats/Math/interface/deltaR.h"
46 +
47   #include "TTree.h"
48  
49   using namespace std;
# Line 57 | Line 62 | struct JRA{
62     float b;
63     float hf;
64     float jtpt[MAXJETS];
65 +   float jtcorpt[MAXJETS];
66     float refpt[MAXJETS];
67     float jteta[MAXJETS];
68     float refeta[MAXJETS];
# Line 66 | Line 72 | struct JRA{
72     float weight;
73   };
74  
75 +
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   //
90   // class declaration
91   //
92 + bool compareCorPt(JRAV a, JRAV b) {return a.jtcorpt > b.jtcorpt;}
93 + bool comparePt(JRAV a, JRAV b) {return a.jtpt > b.jtpt;}
94  
95   class HiJetResponseAnalyzer : public edm::EDAnalyzer {
96     public:
# Line 80 | Line 102 | class HiJetResponseAnalyzer : public edm
102        virtual void beginJob() ;
103        virtual void analyze(const edm::Event&, const edm::EventSetup&);
104        virtual void endJob() ;
105 <
105 >   bool selectJet(int i);
106        // ----------member data ---------------------------
107  
108 <   double genJetPtMin_;
108 >  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 >
118 >  double matchR_;  
119 >   double genPtMin_;
120 >   double ptMin_;
121 >   double emfMin_;
122 >   double n90Min_;
123 >   double n90hitMin_;
124 >
125 >   edm::InputTag jetTag_;
126 >  edm::InputTag matchTag_;
127  
128     JRA jra_;
129     TTree* t;
# Line 92 | Line 132 | class HiJetResponseAnalyzer : public edm
132     edm::Handle<reco::Centrality> cent;
133  
134     edm::Handle<reco::JetView> jets;
135 +   edm::Handle<pat::JetCollection> patjets;
136 +  edm::Handle<reco::JetView> matchedJets;
137  
138 +   edm::Service<TFileService> fs;
139  
140 + };
141  
142 + 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  
99 };
157  
158   //
159   // constants, enums and typedefs
# Line 113 | Line 170 | HiJetResponseAnalyzer::HiJetResponseAnal
170  
171   {
172     //now do what ever initialization is needed
173 <   genJetPtMin_ = 0;
173 >  matchR_ = iConfig.getUntrackedParameter<double>("matchR",0.25);
174 >
175 >   ptMin_ = iConfig.getUntrackedParameter<double>("ptMin",0);
176 >   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 >
181 >   filterJets_ = iConfig.getUntrackedParameter<bool>("filterJets",true);
182 >   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 >   usePat_ = iConfig.getUntrackedParameter<bool>("usePat",true);
189 >   doMC_ = iConfig.getUntrackedParameter<bool>("doMC",true);
190 >
191 >   sortJets_ = iConfig.getUntrackedParameter<bool>("sortJets",true);
192 >   correctJets_ = iConfig.getUntrackedParameter<bool>("correctJets",false);
193 >
194 >   jetTag_ = iConfig.getUntrackedParameter<edm::InputTag>("src",edm::InputTag("selectedPatJets"));
195 >   matchTag_ = iConfig.getUntrackedParameter<edm::InputTag>("match",edm::InputTag("selectedPatJets"));
196   }
197  
198  
# Line 136 | Line 215 | HiJetResponseAnalyzer::analyze(const edm
215   {
216     using namespace edm;
217  
218 +   iEvent.getByLabel(jetTag_,jets);
219 +   if(usePat_)iEvent.getByLabel(jetTag_,patjets);
220 +   if(matchNew_)iEvent.getByLabel(matchTag_,matchedJets);
221  
222 <   for(unsigned int j = 0 ; j < jets->size(); ++j){
222 >   std::vector<JRAV> jraV;
223  
224 <      const pat::Jet jet = (pat::Jet)((*jets)[j]);
225 <      
226 <      if(jet.genJet() != 0){
227 <
228 <         if(jet.genJet()->pt() < genJetPtMin_) continue;
229 <         jra_.refpt[jra_.nref] = jet.genJet()->pt();
230 <         jra_.refeta[jra_.nref] = jet.genJet()->eta();
231 <         jra_.refphi[jra_.nref] = jet.genJet()->phi();
232 <          
233 <         jra_.jtpt[jra_.nref] = jet.pt();
234 <         jra_.jteta[jra_.nref] = jet.eta();
235 <         jra_.jtphi[jra_.nref] = jet.phi();
236 <          
237 <         jra_.nref++;
238 <          
239 <      }
224 >   for(unsigned int j = 0 ; j < jets->size(); ++j){
225 >     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 <   t->Fill();
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  
272 +   for(unsigned int i = 0; i < jraV.size(); ++i){
273 +     JRAV& jv = jraV[i];
274 +     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 +   }
282 +   jra_.nref = jraV.size();
283  
284 +   t->Fill();
285  
164 #ifdef THIS_IS_AN_EVENT_EXAMPLE
165   Handle<ExampleData> pIn;
166   iEvent.getByLabel("example",pIn);
167 #endif
168  
169 #ifdef THIS_IS_AN_EVENTSETUP_EXAMPLE
170   ESHandle<SetupData> pSetup;
171   iSetup.get<SetupRecord>().get(pSetup);
172 #endif
286   }
287  
288  
# Line 177 | Line 290 | HiJetResponseAnalyzer::analyze(const edm
290   void
291   HiJetResponseAnalyzer::beginJob()
292   {
293 +
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 +   t->Branch("jtcorpt",jra_.jtcorpt,"jtcorpt[nref]/F");
300 +   t->Branch("refpt",jra_.refpt,"refpt[nref]/F");
301 +   t->Branch("refcorpt",jra_.refpt,"refcorpt[nref]/F");
302 +   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   }
312  
313   // ------------ method called once each job just after ending the event loop  ------------

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines