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.2 by yilmaz, Fri Sep 10 12:54:47 2010 UTC vs.
Revision 1.8 by yilmaz, Thu Oct 21 16:24:04 2010 UTC

# Line 42 | Line 42
42  
43   #include "CommonTools/UtilAlgos/interface/TFileService.h"
44   #include "FWCore/ServiceRegistry/interface/Service.h"
45 <
45 > #include "DataFormats/Math/interface/deltaR.h"
46  
47   #include "TTree.h"
48  
# Line 62 | 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 71 | Line 72 | struct JRA{
72     float weight;
73   };
74  
75 + struct JRAV{
76 +  int index;
77 +  float jtpt;
78 +  float jtcorpt;
79 +  float refpt;
80 +  float refcorpt;
81 +  float jteta;
82 +  float refeta;
83 +  float jtphi;
84 +  float refphi;
85 +
86 + };
87 +
88   //
89   // class declaration
90   //
91 + bool compareCorPt(JRAV a, JRAV b) {return a.jtcorpt > b.jtcorpt;}
92 + bool comparePt(JRAV a, JRAV b) {return a.jtpt > b.jtpt;}
93  
94   class HiJetResponseAnalyzer : public edm::EDAnalyzer {
95     public:
# Line 85 | Line 101 | class HiJetResponseAnalyzer : public edm
101        virtual void beginJob() ;
102        virtual void analyze(const edm::Event&, const edm::EventSetup&);
103        virtual void endJob() ;
104 <
104 >   bool selectJet(int i);
105        // ----------member data ---------------------------
106  
107 <   bool usePat_;
108 <   bool doMC_;
107 >  bool usePat_;
108 >  bool doMC_;
109 >  bool filterJets_;
110 >  bool diJetsOnly_;
111 >  bool matchDiJets_;
112 >  bool matchPatGen_;
113 >  bool matchNew_;
114 >  bool sortJets_;
115 >  bool correctJets_;
116 >
117 >  double matchR_;  
118 >   double genPtMin_;
119 >   double ptMin_;
120 >   double emfMin_;
121 >   double n90Min_;
122 >   double n90hitMin_;
123 >
124 >  edm::InputTag jetTag_;
125 >  edm::InputTag matchTag_;
126 >  std::vector<edm::InputTag> matchTags_;
127 >  
128 >  JRA jra_;
129 >  std::vector<JRA> jraMatch_;
130  
94   double genJetPtMin_;
95
96
97   edm::InputTag jetTag_;
98
99   JRA jra_;
131     TTree* t;
132  
133     edm::Handle<edm::GenHIEvent> mc;
134     edm::Handle<reco::Centrality> cent;
135  
136     edm::Handle<reco::JetView> jets;
137 <   //   edm::Handle<pat::JetCollection> jets;
138 <
137 >   edm::Handle<pat::JetCollection> patjets;
138 >  edm::Handle<reco::JetView> matchedJets;
139  
140     edm::Service<TFileService> fs;
141  
142   };
143  
144 + bool HiJetResponseAnalyzer::selectJet(int i){
145 +   const reco::Jet& jet = (*jets)[i];
146 +   if(usePat_){
147 +      const pat::Jet& patjet = (*patjets)[i];
148 +      if(patjet.emEnergyFraction() <= emfMin_) return false;
149 +      if(patjet.jetID().n90Hits <= n90hitMin_) return false;
150 +      if(doMC_){
151 +        
152 +      }
153 +
154 +   }
155 +
156 +   return true;
157 + }
158 +
159 +
160   //
161   // constants, enums and typedefs
162   //
# Line 125 | Line 172 | HiJetResponseAnalyzer::HiJetResponseAnal
172  
173   {
174     //now do what ever initialization is needed
175 <   genJetPtMin_ = 0;
129 <
130 <   usePat_ = false;
131 <   doMC_ = false;
132 <
133 <   jetTag_ = edm::InputTag("selectedPatJets");
175 >  matchR_ = iConfig.getUntrackedParameter<double>("matchR",0.25);
176  
177 +   ptMin_ = iConfig.getUntrackedParameter<double>("ptMin",0);
178 +   genPtMin_ = iConfig.getUntrackedParameter<double>("genPtMin",20);
179 +   emfMin_ = iConfig.getUntrackedParameter<double>("emfMin",0.01);
180 +   n90Min_ = iConfig.getUntrackedParameter<double>("n90Min",1);
181 +   n90hitMin_ = iConfig.getUntrackedParameter<double>("n90hitMin",1);
182 +
183 +   filterJets_ = iConfig.getUntrackedParameter<bool>("filterJets",true);
184 +   diJetsOnly_ = iConfig.getUntrackedParameter<bool>("diJetsOnly",false);
185 +   matchDiJets_ = iConfig.getUntrackedParameter<bool>("matchDiJets",false);
186 +   matchPatGen_ = iConfig.getUntrackedParameter<bool>("matchPatGen",false);
187 +
188 +   matchNew_ = iConfig.getUntrackedParameter<bool>("matchNew",false);
189 +
190 +   usePat_ = iConfig.getUntrackedParameter<bool>("usePat",true);
191 +   doMC_ = iConfig.getUntrackedParameter<bool>("doMC",true);
192 +
193 +   sortJets_ = iConfig.getUntrackedParameter<bool>("sortJets",true);
194 +   correctJets_ = iConfig.getUntrackedParameter<bool>("correctJets",false);
195 +
196 +   jetTag_ = iConfig.getUntrackedParameter<edm::InputTag>("src",edm::InputTag("selectedPatJets"));
197 +   matchTag_ = iConfig.getUntrackedParameter<edm::InputTag>("match",edm::InputTag("selectedPatJets"));
198 +   matchTags_ = iConfig.getUntrackedParameter<std::vector<edm::InputTag> >("matches",std::vector<edm::InputTag>(0));
199   }
200  
201  
# Line 155 | Line 219 | HiJetResponseAnalyzer::analyze(const edm
219     using namespace edm;
220  
221     iEvent.getByLabel(jetTag_,jets);
222 <   cout<<"a"<<endl;
222 >   if(usePat_)iEvent.getByLabel(jetTag_,patjets);
223 >   std::vector<JRAV> jraV;
224  
160   jra_.nref = 0;
225     for(unsigned int j = 0 ; j < jets->size(); ++j){
226 <      cout<<"b"<<endl;
227 <      
228 <      //      const pat::Jet& jet = (*jets)[j];
229 <      const reco::Jet jet = (*jets)[j];      
230 <      jra_.jtpt[jra_.nref] = jet.pt();
231 <      jra_.jteta[jra_.nref] = jet.eta();
232 <      jra_.jtphi[jra_.nref] = jet.phi();
233 <      cout<<"c"<<endl;
234 <      
235 <      if(usePat_){
236 <         const pat::Jet& patjet = (*jets)[j];
237 <        
238 <         if(doMC_ && patjet.genJet() != 0){        
239 <            //      if(jet.genJet()->pt() < genJetPtMin_) continue;
240 <            jra_.refpt[jra_.nref] = patjet.genJet()->pt();
241 <            jra_.refeta[jra_.nref] = patjet.genJet()->eta();
242 <            jra_.refphi[jra_.nref] = patjet.genJet()->phi();
243 <                            
244 <         }else{
245 <            jra_.refpt[jra_.nref] = -99;
246 <            jra_.refeta[jra_.nref] = -99;
247 <            jra_.refphi[jra_.nref] = -99;
248 <         }
249 <      }
250 <      
251 <      jra_.nref++;
252 <      
226 >     if(filterJets_ && !selectJet(j)) continue;
227 >     const reco::Jet& jet = (*jets)[j];
228 >     JRAV jv;
229 >     jv.jtpt = jet.pt();
230 >     jv.jteta = jet.eta();
231 >     jv.jtphi = jet.phi();
232 >     jv.jtcorpt = jet.pt();
233 >     jv.index = j;
234 >     if(usePat_){
235 >       const pat::Jet& patjet = (*patjets)[j];
236 >
237 >       jv.jtpt = patjet.correctedJet("raw").pt();
238 >       jv.jtcorpt = patjet.pt();
239 >
240 >       if(doMC_ && matchPatGen_ && patjet.genJet() != 0){
241 >         if(patjet.genJet()->pt() < genPtMin_) continue;
242 >         jv.refpt = patjet.genJet()->pt();
243 >         jv.refeta = patjet.genJet()->eta();
244 >         jv.refphi = patjet.genJet()->phi();
245 >       }else{
246 >         jv.refpt = -99;
247 >         jv.refeta = -99;
248 >         jv.refphi = -99;
249 >       }
250 >     }
251 >
252 >     jraV.push_back(jv);
253     }
254  
255 +   if(sortJets_){
256 +     if(usePat_ || correctJets_) std::sort(jraV.begin(),jraV.end(),compareCorPt);
257 +     else std::sort(jraV.begin(),jraV.end(),comparePt);
258 +   }
259 +
260 +   for(unsigned int i = 0; i < jraV.size(); ++i){
261 +     JRAV& jv = jraV[i];
262 +     const reco::Jet& jet = (*jets)[jv.index];
263 +    
264 +     if(matchNew_){
265 +       for(unsigned int im = 0; im < matchTags_.size(); ++im){
266 +         iEvent.getByLabel(matchTags_[im],matchedJets);
267 +         for(unsigned int m = 0 ; m < matchedJets->size(); ++m){
268 +           const reco::Jet& match = (*matchedJets)[m];
269 +           double dr = reco::deltaR(jet.eta(),jet.phi(),match.eta(),match.phi());
270 +           if(dr < matchR_){
271 +             jraMatch_[im].jtcorpt[i] = -99;
272 +             jraMatch_[im].jtpt[i] = match.pt();
273 +             jraMatch_[im].jteta[i] = match.eta();
274 +             jraMatch_[im].jtphi[i] = match.phi();
275 +           }
276 +         }
277 +       }
278 +     }
279 +    
280 +     jra_.jtpt[i] = jv.jtpt;
281 +     jra_.jteta[i] = jv.jteta;
282 +     jra_.jtphi[i] = jv.jtphi;
283 +     jra_.jtcorpt[i] = jv.jtcorpt;
284 +     jra_.refpt[i] = jv.refpt;
285 +     jra_.refeta[i] = jv.refeta;
286 +     jra_.refphi[i] = jv.refphi;
287 +   }
288 +   jra_.nref = jraV.size();
289 +  
290     t->Fill();
291  
292   }
293  
195
294   // ------------ method called once each job just before starting event loop  ------------
197 void
198 HiJetResponseAnalyzer::beginJob()
199 {
295  
296 <   t= fs->make<TTree>("t","Jet Response Analyzer");
297 <   t->Branch("b",&jra_.b,"b/F");
298 <   t->Branch("hf",&jra_.hf,"hf/F");
299 <   t->Branch("nref",&jra_.nref,"nref/I");
300 <   t->Branch("jtpt",jra_.jtpt,"jtpt[nref]/F");
301 <   t->Branch("refpt",jra_.refpt,"refpt[nref]/F");
302 <   t->Branch("jteta",jra_.jteta,"jteta[nref]/F");
296 >
297 > void
298 > HiJetResponseAnalyzer::beginJob(){
299 >  t= fs->make<TTree>("t","Jet Response Analyzer");
300 >  t->Branch("b",&jra_.b,"b/F");
301 >  t->Branch("hf",&jra_.hf,"hf/F");
302 >  t->Branch("nref",&jra_.nref,"nref/I");
303 >  t->Branch("jtpt",jra_.jtpt,"jtpt[nref]/F");
304 >  t->Branch("jtcorpt",jra_.jtcorpt,"jtcorpt[nref]/F");
305 >  t->Branch("jteta",jra_.jteta,"jteta[nref]/F");
306 >  t->Branch("jtphi",jra_.jtphi,"jtphi[nref]/F");
307 >  t->Branch("refpt",jra_.refpt,"refpt[nref]/F");
308 >  t->Branch("refcorpt",jra_.refpt,"refcorpt[nref]/F");
309     t->Branch("refeta",jra_.refeta,"refeta[nref]/F");
209   t->Branch("jtphi",jra_.jtphi,"jtphi[nref]/F");
310     t->Branch("refphi",jra_.refphi,"refphi[nref]/F");
311     t->Branch("weight",&jra_.weight,"weight/F");
312     t->Branch("bin",&jra_.bin,"bin/I");
313 <
314 <
315 <
313 >   for(unsigned int im = 0; im < matchTags_.size(); ++im){
314 >     JRA jrm;
315 >     jraMatch_.push_back(jrm);
316 >     t->Branch(Form("jtpt%d",im),jraMatch_[im].jtpt,Form("jtpt%d[nref]/F",im));
317 >     t->Branch(Form("jtcorpt%d",im),jraMatch_[im].jtcorpt,Form("jtcorpt%d[nref]/F",im));
318 >     t->Branch(Form("jteta%d",im),jraMatch_[im].jteta,Form("jteta%d[nref]/F",im));
319 >     t->Branch(Form("jtphi%d",im),jraMatch_[im].jtphi,Form("jtphi%d[nref]/F",im));
320 >   }
321 >  
322   }
323  
324   // ------------ method called once each job just after ending the event loop  ------------

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines