ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/csander/JetResolutionFromMC/MCResolutions/src/MCResolutions.cc
Revision: 1.16
Committed: Thu Jul 26 09:48:00 2012 UTC (12 years, 9 months ago) by kaussen
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 1.15: +44 -174 lines
Log Message: 
Some code reduction

File Contents

# User Rev Content
1 csander 1.1 // -*- C++ -*-
2     //
3     // Package: MCResolutions
4     // Class: MCResolutions
5 csander 1.2 //
6 csander 1.1 /**\class MCResolutions MCResolutions.cc JetResolutionFromMC/MCResolutions/src/MCResolutions.cc
7 csander 1.2
8 csander 1.1 Description: [one line class summary]
9 csander 1.2
10 csander 1.1 Implementation:
11 csander 1.10 [Notes on implementation]
12     */
13 csander 1.1 //
14     // Original Author: Christian Sander,,,
15     // Created: Wed Oct 6 18:22:21 CEST 2010
16 kaussen 1.16 // $Id: MCResolutions.cc,v 1.15 2012/07/25 09:44:33 kaussen Exp $
17 csander 1.1 //
18     //
19    
20 csander 1.2 #include "JetResolutionFromMC/MCResolutions/interface/MCResolutions.h"
21 csander 1.1
22     //
23     // constructors and destructor
24     //
25 csander 1.3 MCResolutions::MCResolutions(const edm::ParameterSet& iConfig) {
26 csander 1.10 //now do what ever initialization is needed
27     _jetTag = iConfig.getParameter<edm::InputTag> ("jetTag");
28     _muonTag = iConfig.getParameter<edm::InputTag> ("muonTag");
29     _genJetTag = iConfig.getParameter<edm::InputTag> ("genJetTag");
30     _weightName = iConfig.getParameter<edm::InputTag> ("weightName");
31     _EBRecHits = iConfig.getParameter<edm::InputTag> ("EBRecHits");
32     _EERecHits = iConfig.getParameter<edm::InputTag> ("EERecHits");
33     _jetMultPtCut = iConfig.getParameter<double> ("jetMultPtCut");
34     _jetMultEtaCut = iConfig.getParameter<double> ("jetMultEtaCut");
35     _deltaPhiDiJet = iConfig.getParameter<double> ("deltaPhiDiJet");
36     _absCut3rdJet = iConfig.getParameter<double> ("absCut3rdJet");
37     _relCut3rdJet = iConfig.getParameter<double> ("relCut3rdJet");
38     _deltaRMatch = iConfig.getParameter<double> ("deltaRMatch");
39     _deltaRMatchVeto = iConfig.getParameter<double> ("deltaRMatchVeto");
40     _absPtVeto = iConfig.getParameter<double> ("absPtVeto");
41     _relPtVeto = iConfig.getParameter<double> ("relPtVeto");
42     _deltaRDeadECal = iConfig.getParameter<double> ("deltaRDeadECal");
43     _GenJetPtCut = iConfig.getParameter<double> ("GenJetPtCut");
44     _maskedEcalChannelStatusThreshold = iConfig.getParameter<int> ("maskedEcalChannelStatusThreshold");
45     _fileName = iConfig.getParameter<std::string> ("fileName");
46    
47     hfile = new TFile(_fileName.c_str(), "RECREATE", "Jet response in pT and eta bins");
48     //hfile->mkdir(_dirName.c_str(), _dirName.c_str());
49    
50     PtBinEdges.push_back(0);
51     PtBinEdges.push_back(20);
52     PtBinEdges.push_back(30);
53     PtBinEdges.push_back(50);
54     PtBinEdges.push_back(80);
55     PtBinEdges.push_back(120);
56     PtBinEdges.push_back(170);
57     PtBinEdges.push_back(230);
58     PtBinEdges.push_back(300);
59     PtBinEdges.push_back(380);
60     PtBinEdges.push_back(470);
61     PtBinEdges.push_back(570);
62     PtBinEdges.push_back(680);
63     PtBinEdges.push_back(800);
64     PtBinEdges.push_back(1000);
65     PtBinEdges.push_back(1300);
66     PtBinEdges.push_back(1700);
67     PtBinEdges.push_back(2200);
68 csander 1.14 PtBinEdges.push_back(2800);
69     PtBinEdges.push_back(3500);
70 csander 1.10
71     EtaBinEdges.push_back(0.0);
72 csander 1.12 EtaBinEdges.push_back(0.3);
73 csander 1.10 EtaBinEdges.push_back(0.5);
74 csander 1.12 EtaBinEdges.push_back(0.8);
75 csander 1.10 EtaBinEdges.push_back(1.1);
76 csander 1.12 EtaBinEdges.push_back(1.4);
77 csander 1.10 EtaBinEdges.push_back(1.7);
78 csander 1.12 EtaBinEdges.push_back(2.0);
79 csander 1.10 EtaBinEdges.push_back(2.3);
80 csander 1.12 EtaBinEdges.push_back(2.8);
81 csander 1.10 EtaBinEdges.push_back(3.2);
82 csander 1.12 EtaBinEdges.push_back(4.1);
83 csander 1.10 EtaBinEdges.push_back(5.0);
84    
85 csander 1.12 // EtaBinEdges.push_back(0.0);
86     // EtaBinEdges.push_back(0.5);
87     // EtaBinEdges.push_back(1.1);
88     // EtaBinEdges.push_back(1.7);
89     // EtaBinEdges.push_back(2.3);
90     // EtaBinEdges.push_back(3.2);
91     // EtaBinEdges.push_back(5.0);
92    
93 csander 1.10 //// Array of histograms for jet resolutions (all jet multiplicities)
94 kaussen 1.16 ResizeHistoVector(h_tot_DiJet_JetResPt_Pt);
95     ResizeHistoVector(h_b_DiJet_JetResPt_Pt);
96     ResizeHistoVector(h_nob_DiJet_JetResPt_Pt);
97     ResizeHistoVector(h_dead_DiJet_JetResPt_Pt);
98     ResizeHistoVector(h_deadb_DiJet_JetResPt_Pt);
99 csander 1.10
100     //// Array of histograms for jet resolutions (all jet multiplicities)
101 kaussen 1.16 ResizeHistoVector(h_tot_JetAll_JetResPt_Pt);
102     ResizeHistoVector(h_b_JetAll_JetResPt_Pt);
103     ResizeHistoVector(h_nob_JetAll_JetResPt_Pt);
104     ResizeHistoVector(h_dead_JetAll_JetResPt_Pt);
105     ResizeHistoVector(h_deadb_JetAll_JetResPt_Pt);
106 csander 1.13
107     //// Array of histograms for jet resolutions (first jet)
108 kaussen 1.16 ResizeHistoVector(h_tot_Jet1_JetResPt_Pt);
109     ResizeHistoVector(h_b_Jet1_JetResPt_Pt);
110     ResizeHistoVector(h_nob_Jet1_JetResPt_Pt);
111     ResizeHistoVector(h_dead_Jet1_JetResPt_Pt);
112     ResizeHistoVector(h_deadb_Jet1_JetResPt_Pt);
113 csander 1.13
114     //// Array of histograms for jet resolutions (second jet)
115 kaussen 1.16 ResizeHistoVector(h_tot_Jet2_JetResPt_Pt);
116     ResizeHistoVector(h_b_Jet2_JetResPt_Pt);
117     ResizeHistoVector(h_nob_Jet2_JetResPt_Pt);
118     ResizeHistoVector(h_dead_Jet2_JetResPt_Pt);
119     ResizeHistoVector(h_deadb_Jet2_JetResPt_Pt);
120 csander 1.13
121 kaussen 1.15 //// Array of histograms for jet resolutions (third jet)
122 kaussen 1.16 ResizeHistoVector(h_tot_Jet3_JetResPt_Pt);
123     ResizeHistoVector(h_b_Jet3_JetResPt_Pt);
124     ResizeHistoVector(h_nob_Jet3_JetResPt_Pt);
125     ResizeHistoVector(h_dead_Jet3_JetResPt_Pt);
126     ResizeHistoVector(h_deadb_Jet3_JetResPt_Pt);
127 kaussen 1.15
128     //// Array of histograms for jet resolutions (fourth jet)
129 kaussen 1.16 ResizeHistoVector(h_tot_Jet4_JetResPt_Pt);
130     ResizeHistoVector(h_b_Jet4_JetResPt_Pt);
131     ResizeHistoVector(h_nob_Jet4_JetResPt_Pt);
132     ResizeHistoVector(h_dead_Jet4_JetResPt_Pt);
133     ResizeHistoVector(h_deadb_Jet4_JetResPt_Pt);
134 kaussen 1.15
135     //// Array of histograms for jet resolutions (fifth+ jet)
136 kaussen 1.16 ResizeHistoVector(h_tot_Jet5p_JetResPt_Pt);
137     ResizeHistoVector(h_b_Jet5p_JetResPt_Pt);
138     ResizeHistoVector(h_nob_Jet5p_JetResPt_Pt);
139     ResizeHistoVector(h_dead_Jet5p_JetResPt_Pt);
140     ResizeHistoVector(h_deadb_Jet5p_JetResPt_Pt);
141 csander 1.1
142     }
143    
144 csander 1.3 MCResolutions::~MCResolutions() {
145 csander 1.2
146 csander 1.10 // do anything here that needs to be done at desctruction time
147     // (e.g. close files, deallocate resources etc.)
148 csander 1.1
149     }
150    
151     //
152     // member functions
153     //
154    
155     // ------------ method called to for each event ------------
156 csander 1.10 void MCResolutions::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
157     using namespace std;
158 csander 1.5
159 csander 1.10 // Event setup
160     envSet(iSetup);
161 csander 1.9
162 csander 1.10 getChannelStatusMaps();
163 csander 1.9
164 csander 1.10 //Weight
165     edm::Handle<double> event_weight;
166     bool findWeight = iEvent.getByLabel(_weightName, event_weight);
167     weight = (event_weight.isValid() ? (*event_weight) : 1.0);
168     if (!findWeight) {
169     cout << "Weight not found!" << endl;
170     }
171    
172     //GenJets
173     edm::Handle<edm::View<reco::GenJet> > Jets_gen;
174     iEvent.getByLabel(_genJetTag, Jets_gen);
175    
176     //RecoJets
177     edm::Handle<edm::View<pat::Jet> > Jets_rec;
178     iEvent.getByLabel(_jetTag, Jets_rec);
179    
180     edm::Handle<reco::GenParticleCollection> genParticles;
181     iEvent.getByLabel("genParticles", genParticles);
182    
183 kaussen 1.15 //// Do DiJet selection
184 csander 1.10 bool isDiJet = false;
185     const reco::GenJet* gj1 = 0;
186     const reco::GenJet* gj2 = 0;
187     const reco::GenJet* gj3 = 0;
188     int k = 0;
189     for (edm::View<reco::GenJet>::const_iterator it = Jets_gen->begin(); it != Jets_gen->end(); ++it) {
190     ++k;
191     if (k == 1)
192     gj1 = &(*it);
193     if (k == 2)
194     gj2 = &(*it);
195     if (k == 3) {
196     gj3 = &(*it);
197     break;
198     }
199     }
200     if (k > 1) {
201     if (deltaPhi(gj1->phi(), gj2->phi()) > _deltaPhiDiJet) {
202     if (k == 2)
203     isDiJet = true;
204     if (k == 3) {
205     //if (gj3->pt()<(_relCut3rdJet*(gj1->pt()+gj2->pt())/2)){
206     if (gj3->pt() < (_relCut3rdJet * gj1->pt())) {
207     isDiJet = true;
208 csander 1.4 }
209 csander 1.10 }
210     }
211     }
212    
213 csander 1.13 //// Calculate jet multiplicity
214 csander 1.10 double NGenJet = 0;
215     for (edm::View<reco::GenJet>::const_iterator it = Jets_gen-> begin(); it != Jets_gen-> end(); ++it) {
216     if (it->pt() > _jetMultPtCut && abs(it->eta()) < _jetMultEtaCut)
217     ++NGenJet;
218     }
219    
220     int JetNo = 0;
221     for (edm::View<reco::GenJet>::const_iterator it = Jets_gen->begin(); it != Jets_gen->end(); ++it) {
222     ++JetNo;
223    
224     if (it->pt() < _GenJetPtCut)
225     continue;
226 csander 1.13
227 csander 1.10 //// First look if there is no significant GenJet near the tested GenJet
228     double dRgenjet = 999.;
229     double PtdRmin = 0;
230     for (edm::View<reco::GenJet>::const_iterator kt = Jets_gen-> begin(); kt != Jets_gen->end(); ++kt) {
231     if (&(*it) != &(*kt))
232 csander 1.3 continue;
233 csander 1.10 double dR = deltaR(*it, *kt);
234     if (dR < dRgenjet) {
235     dRgenjet = dR;
236     PtdRmin = kt->pt();
237     }
238     }
239     if (dRgenjet < _deltaRMatchVeto && PtdRmin / it->pt() < _relPtVeto)
240     continue;
241     const pat::Jet* matchedJet = 0;
242     const pat::Jet* nearestJet = 0;
243     double dRmatched = 999.;
244     double dRnearest = 999.;
245     for (edm::View<pat::Jet>::const_iterator jt = Jets_rec-> begin(); jt != Jets_rec->end(); ++jt) {
246     double dR = deltaR(*it, *jt);
247     if (dR < dRmatched) {
248     nearestJet = matchedJet;
249     dRnearest = dRmatched;
250     matchedJet = &(*jt);
251     dRmatched = dR;
252     } else if (dR < dRnearest) {
253     nearestJet = &(*jt);
254     dRnearest = dR;
255     }
256     }
257 csander 1.13
258 csander 1.10 //// look if there is no further significant CaloJet near the genJet
259     if (dRmatched < _deltaRMatch && (nearestJet == 0 || dRnearest > _deltaRMatchVeto || (nearestJet->pt()
260     < _absPtVeto && nearestJet->pt() / matchedJet->pt() < _relPtVeto))) {
261     double res = matchedJet->pt() / it->pt();
262 csander 1.13 h_tot_JetAll_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
263 csander 1.10 if (isDiJet && JetNo < 3) {
264     h_tot_DiJet_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
265     }
266 csander 1.13 if (JetNo == 1) {
267     h_tot_Jet1_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
268     }
269     if (JetNo == 2) {
270     h_tot_Jet2_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
271     }
272 kaussen 1.15 if (JetNo == 3) {
273     h_tot_Jet3_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
274     }
275     if (JetNo == 4) {
276     h_tot_Jet4_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
277     }
278     if (JetNo > 4) {
279     h_tot_Jet5p_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
280 csander 1.13 }
281    
282     //// Use algorithmic matching for heavy flavour ID
283 csander 1.10 bool bTag = false;
284 csander 1.13 if (fabs(matchedJet->partonFlavour()) == 4 || fabs(matchedJet->partonFlavour()) == 5) {
285     bTag = true;
286 csander 1.10 }
287 csander 1.6
288 csander 1.10 if (bTag) {
289 csander 1.13 h_b_JetAll_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
290 csander 1.10 if (isDiJet && JetNo < 3) {
291     h_b_DiJet_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
292     }
293 csander 1.13 if (JetNo == 1) {
294     h_b_Jet1_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
295     }
296     if (JetNo == 2) {
297     h_b_Jet2_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
298     }
299 kaussen 1.15 if (JetNo == 3) {
300     h_b_Jet3_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
301     }
302     if (JetNo == 4) {
303     h_b_Jet4_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
304     }
305     if (JetNo > 4) {
306     h_b_Jet5p_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
307 csander 1.13 }
308 csander 1.10 } else {
309 csander 1.13 h_nob_JetAll_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
310 csander 1.10 if (isDiJet && JetNo < 3) {
311     h_nob_DiJet_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
312 csander 1.6 }
313 csander 1.13 if (JetNo == 1) {
314     h_nob_Jet1_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
315     }
316     if (JetNo == 2) {
317     h_nob_Jet2_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
318     }
319 kaussen 1.15 if (JetNo == 3) {
320     h_nob_Jet3_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
321     }
322     if (JetNo == 4) {
323     h_nob_Jet4_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
324     }
325     if (JetNo > 4) {
326     h_nob_Jet5p_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
327 csander 1.13 }
328 csander 1.10 }
329    
330     //// look if matched jet points into direction of masked ECAL cluster
331     for (std::map<DetId, std::vector<double> >::iterator kt = EcalAllDeadChannelsValMap.begin(); kt
332     != EcalAllDeadChannelsValMap.end(); ++kt) {
333     math::PtEtaPhiMLorentzVectorD Evec(100., kt->second.at(0), kt->second.at(1), 0.);
334     double dR_dead = deltaR(*matchedJet, Evec);
335     if (dR_dead < _deltaRDeadECal) {
336 csander 1.13 h_dead_JetAll_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
337 csander 1.10 if (isDiJet && JetNo < 3) {
338     h_dead_DiJet_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
339 csander 1.13 }
340     if (JetNo == 1) {
341     h_dead_Jet1_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
342     }
343     if (JetNo == 2) {
344     h_dead_Jet2_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
345     }
346 kaussen 1.15 if (JetNo == 3) {
347     h_dead_Jet3_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
348     }
349     if (JetNo == 4) {
350     h_dead_Jet4_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
351     }
352     if (JetNo > 4) {
353     h_dead_Jet5p_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
354 csander 1.13 }
355     if (bTag) {
356     h_deadb_JetAll_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
357     if (isDiJet && JetNo < 3) {
358 csander 1.10 h_deadb_DiJet_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
359 csander 1.13 }
360     if (JetNo == 1) {
361     h_deadb_Jet1_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
362     }
363     if (JetNo == 2) {
364     h_deadb_Jet2_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
365     }
366 kaussen 1.15 if (JetNo == 3) {
367     h_deadb_Jet3_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
368     }
369     if (JetNo == 4) {
370     h_deadb_Jet4_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
371     }
372     if (JetNo > 4) {
373     h_deadb_Jet5p_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
374 csander 1.10 }
375     }
376     break;
377 csander 1.3 }
378 csander 1.10 }
379     }
380     }
381 csander 1.1
382     }
383    
384     // ------------ method called once each job just before starting event loop ------------
385 csander 1.10 void MCResolutions::beginJob() {
386 csander 1.3
387 csander 1.10 for (unsigned int i_pt = 0; i_pt < PtBinEdges.size() - 1; ++i_pt) {
388     for (unsigned int i_eta = 0; i_eta < EtaBinEdges.size() - 1; ++i_eta) {
389     char hname[100];
390     //// Book histograms (all jet multiplicities)
391 csander 1.13 sprintf(hname, "h_tot_JetAll_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
392     h_tot_JetAll_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
393     h_tot_JetAll_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
394     sprintf(hname, "h_b_JetAll_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
395     h_b_JetAll_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
396     h_b_JetAll_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
397     sprintf(hname, "h_nob_JetAll_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
398     h_nob_JetAll_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
399     h_nob_JetAll_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
400     sprintf(hname, "h_dead_JetAll_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
401     h_dead_JetAll_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
402     h_dead_JetAll_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
403     sprintf(hname, "h_deadb_JetAll_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
404     h_deadb_JetAll_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
405     h_deadb_JetAll_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
406 csander 1.10 //// Book histograms (DiJets)
407     sprintf(hname, "h_tot_DiJet_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
408     h_tot_DiJet_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
409     h_tot_DiJet_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
410     sprintf(hname, "h_b_DiJet_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
411     h_b_DiJet_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
412     h_b_DiJet_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
413     sprintf(hname, "h_nob_DiJet_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
414     h_nob_DiJet_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
415     h_nob_DiJet_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
416     sprintf(hname, "h_dead_DiJet_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
417     h_dead_DiJet_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
418     h_dead_DiJet_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
419     sprintf(hname, "h_deadb_DiJet_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
420     h_deadb_DiJet_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
421     h_deadb_DiJet_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
422 csander 1.13 //// Book histograms (leading jet)
423     sprintf(hname, "h_tot_Jet1_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
424     h_tot_Jet1_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
425     h_tot_Jet1_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
426     sprintf(hname, "h_b_Jet1_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
427     h_b_Jet1_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
428     h_b_Jet1_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
429     sprintf(hname, "h_nob_Jet1_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
430     h_nob_Jet1_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
431     h_nob_Jet1_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
432     sprintf(hname, "h_dead_Jet1_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
433     h_dead_Jet1_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
434     h_dead_Jet1_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
435     sprintf(hname, "h_deadb_Jet1_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
436     h_deadb_Jet1_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
437     h_deadb_Jet1_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
438     //// Book histograms (second jet)
439     sprintf(hname, "h_tot_Jet2_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
440     h_tot_Jet2_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
441     h_tot_Jet2_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
442     sprintf(hname, "h_b_Jet2_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
443     h_b_Jet2_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
444     h_b_Jet2_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
445     sprintf(hname, "h_nob_Jet2_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
446     h_nob_Jet2_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
447     h_nob_Jet2_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
448     sprintf(hname, "h_dead_Jet2_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
449     h_dead_Jet2_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
450     h_dead_Jet2_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
451     sprintf(hname, "h_deadb_Jet2_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
452     h_deadb_Jet2_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
453     h_deadb_Jet2_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
454 kaussen 1.15 //// Book histograms (third jet)
455     sprintf(hname, "h_tot_Jet3_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
456     h_tot_Jet3_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
457     h_tot_Jet3_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
458     sprintf(hname, "h_b_Jet3_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
459     h_b_Jet3_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
460     h_b_Jet3_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
461     sprintf(hname, "h_nob_Jet3_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
462     h_nob_Jet3_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
463     h_nob_Jet3_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
464     sprintf(hname, "h_dead_Jet3_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
465     h_dead_Jet3_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
466     h_dead_Jet3_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
467     sprintf(hname, "h_deadb_Jet3_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
468     h_deadb_Jet3_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
469     h_deadb_Jet3_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
470     //// Book histograms (fourth jet)
471     sprintf(hname, "h_tot_Jet4_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
472     h_tot_Jet4_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
473     h_tot_Jet4_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
474     sprintf(hname, "h_b_Jet4_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
475     h_b_Jet4_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
476     h_b_Jet4_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
477     sprintf(hname, "h_nob_Jet4_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
478     h_nob_Jet4_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
479     h_nob_Jet4_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
480     sprintf(hname, "h_dead_Jet4_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
481     h_dead_Jet4_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
482     h_dead_Jet4_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
483     sprintf(hname, "h_deadb_Jet4_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
484     h_deadb_Jet4_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
485     h_deadb_Jet4_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
486     //// Book histograms (>= fifth jet)
487     sprintf(hname, "h_tot_Jet5p_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
488     h_tot_Jet5p_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
489     h_tot_Jet5p_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
490     sprintf(hname, "h_b_Jet5p_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
491     h_b_Jet5p_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
492     h_b_Jet5p_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
493     sprintf(hname, "h_nob_Jet5p_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
494     h_nob_Jet5p_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
495     h_nob_Jet5p_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
496     sprintf(hname, "h_dead_Jet5p_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
497     h_dead_Jet5p_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
498     h_dead_Jet5p_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
499     sprintf(hname, "h_deadb_Jet5p_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
500     h_deadb_Jet5p_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
501     h_deadb_Jet5p_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
502 csander 1.10 }
503     }
504    
505     mapsReady = false;
506 csander 1.3 }
507 csander 1.1
508     // ------------ method called once each job just after ending the event loop ------------
509 csander 1.10 void MCResolutions::endJob() {
510 csander 1.5
511 csander 1.10 hfile->cd();
512     //hfile->cd(_dirName.c_str());
513     // Save all objects in this file
514     for (unsigned int i_pt = 0; i_pt < PtBinEdges.size() - 1; ++i_pt) {
515     for (unsigned int i_eta = 0; i_eta < EtaBinEdges.size() - 1; ++i_eta) {
516     // total
517 csander 1.13 hfile->WriteTObject(h_tot_JetAll_JetResPt_Pt.at(i_eta).at(i_pt));
518 csander 1.10 hfile->WriteTObject(h_tot_DiJet_JetResPt_Pt.at(i_eta).at(i_pt));
519 csander 1.13 hfile->WriteTObject(h_tot_Jet1_JetResPt_Pt.at(i_eta).at(i_pt));
520     hfile->WriteTObject(h_tot_Jet2_JetResPt_Pt.at(i_eta).at(i_pt));
521 kaussen 1.15 hfile->WriteTObject(h_tot_Jet3_JetResPt_Pt.at(i_eta).at(i_pt));
522     hfile->WriteTObject(h_tot_Jet4_JetResPt_Pt.at(i_eta).at(i_pt));
523     hfile->WriteTObject(h_tot_Jet5p_JetResPt_Pt.at(i_eta).at(i_pt));
524 csander 1.10 // with btag
525 csander 1.13 hfile->WriteTObject(h_b_JetAll_JetResPt_Pt.at(i_eta).at(i_pt));
526 csander 1.10 hfile->WriteTObject(h_b_DiJet_JetResPt_Pt.at(i_eta).at(i_pt));
527 csander 1.13 hfile->WriteTObject(h_b_Jet1_JetResPt_Pt.at(i_eta).at(i_pt));
528     hfile->WriteTObject(h_b_Jet2_JetResPt_Pt.at(i_eta).at(i_pt));
529 kaussen 1.15 hfile->WriteTObject(h_b_Jet3_JetResPt_Pt.at(i_eta).at(i_pt));
530     hfile->WriteTObject(h_b_Jet4_JetResPt_Pt.at(i_eta).at(i_pt));
531     hfile->WriteTObject(h_b_Jet5p_JetResPt_Pt.at(i_eta).at(i_pt));
532 csander 1.10 // without btag
533 csander 1.13 hfile->WriteTObject(h_nob_JetAll_JetResPt_Pt.at(i_eta).at(i_pt));
534 csander 1.10 hfile->WriteTObject(h_nob_DiJet_JetResPt_Pt.at(i_eta).at(i_pt));
535 csander 1.13 hfile->WriteTObject(h_nob_Jet1_JetResPt_Pt.at(i_eta).at(i_pt));
536     hfile->WriteTObject(h_nob_Jet2_JetResPt_Pt.at(i_eta).at(i_pt));
537 kaussen 1.15 hfile->WriteTObject(h_nob_Jet3_JetResPt_Pt.at(i_eta).at(i_pt));
538     hfile->WriteTObject(h_nob_Jet4_JetResPt_Pt.at(i_eta).at(i_pt));
539     hfile->WriteTObject(h_nob_Jet5p_JetResPt_Pt.at(i_eta).at(i_pt));
540 csander 1.10 // in direction of dead ECAL cells
541 csander 1.13 hfile->WriteTObject(h_dead_JetAll_JetResPt_Pt.at(i_eta).at(i_pt));
542 csander 1.10 hfile->WriteTObject(h_dead_DiJet_JetResPt_Pt.at(i_eta).at(i_pt));
543 csander 1.13 hfile->WriteTObject(h_dead_Jet1_JetResPt_Pt.at(i_eta).at(i_pt));
544     hfile->WriteTObject(h_dead_Jet2_JetResPt_Pt.at(i_eta).at(i_pt));
545 kaussen 1.15 hfile->WriteTObject(h_dead_Jet3_JetResPt_Pt.at(i_eta).at(i_pt));
546     hfile->WriteTObject(h_dead_Jet4_JetResPt_Pt.at(i_eta).at(i_pt));
547     hfile->WriteTObject(h_dead_Jet5p_JetResPt_Pt.at(i_eta).at(i_pt));
548 csander 1.10 // in direction of dead ECAL cells and with b tag
549 csander 1.13 hfile->WriteTObject(h_deadb_JetAll_JetResPt_Pt.at(i_eta).at(i_pt));
550 csander 1.10 hfile->WriteTObject(h_deadb_DiJet_JetResPt_Pt.at(i_eta).at(i_pt));
551 csander 1.13 hfile->WriteTObject(h_deadb_Jet1_JetResPt_Pt.at(i_eta).at(i_pt));
552     hfile->WriteTObject(h_deadb_Jet2_JetResPt_Pt.at(i_eta).at(i_pt));
553 kaussen 1.15 hfile->WriteTObject(h_deadb_Jet3_JetResPt_Pt.at(i_eta).at(i_pt));
554     hfile->WriteTObject(h_deadb_Jet4_JetResPt_Pt.at(i_eta).at(i_pt));
555     hfile->WriteTObject(h_deadb_Jet5p_JetResPt_Pt.at(i_eta).at(i_pt));
556 csander 1.10 }
557     }
558    
559     hfile->cd();
560     hfile->WriteObject(&PtBinEdges, "PtBinEdges");
561     hfile->WriteObject(&EtaBinEdges, "EtaBinEdges");
562     //hfile->ls();
563 csander 1.5
564 csander 1.10 // Close the file.
565     hfile->Close();
566 csander 1.5
567     }
568 csander 1.3
569     int MCResolutions::PtBin(const double& pt) {
570 csander 1.10 int i_pt = -1;
571     for (std::vector<double>::const_iterator it = PtBinEdges.begin(); it != PtBinEdges.end(); ++it) {
572     if ((*it) > pt)
573     break;
574     ++i_pt;
575     }
576     if (i_pt < 0)
577     i_pt = 0;
578     if (i_pt > (int) PtBinEdges.size() - 2)
579     i_pt = (int) PtBinEdges.size() - 2;
580 csander 1.3
581 csander 1.10 return i_pt;
582 csander 1.3 }
583    
584     int MCResolutions::EtaBin(const double& eta) {
585 csander 1.10 int i_eta = -1;
586     for (std::vector<double>::const_iterator it = EtaBinEdges.begin(); it != EtaBinEdges.end(); ++it) {
587     if ((*it) > fabs(eta))
588     break;
589     ++i_eta;
590     }
591     if (i_eta < 0)
592     i_eta = 0;
593     if (i_eta > (int) EtaBinEdges.size() - 2)
594     i_eta = (int) EtaBinEdges.size() - 2;
595     return i_eta;
596 csander 1.3 }
597    
598     void MCResolutions::envSet(const edm::EventSetup& iSetup) {
599    
600 csander 1.10 ecalScale.setEventSetup(iSetup);
601 csander 1.3
602 csander 1.10 iSetup.get<EcalChannelStatusRcd> ().get(ecalStatus);
603     iSetup.get<CaloGeometryRecord> ().get(geometry);
604 csander 1.3
605 csander 1.10 if (!ecalStatus.isValid())
606     throw "Failed to get ECAL channel status!";
607     if (!geometry.isValid())
608     throw "Failed to get the geometry!";
609 csander 1.3
610     }
611    
612     int MCResolutions::getChannelStatusMaps() {
613    
614 csander 1.10 if (mapsReady)
615     return -1;
616 csander 1.3
617 csander 1.10 EcalAllDeadChannelsValMap.clear();
618 csander 1.3
619 csander 1.10 // Loop over EB ...
620     for (int ieta = -85; ieta <= 85; ieta++) {
621     for (int iphi = 0; iphi <= 360; iphi++) {
622     if (!EBDetId::validDetId(ieta, iphi))
623     continue;
624    
625     const EBDetId detid = EBDetId(ieta, iphi, EBDetId::ETAPHIMODE);
626     EcalChannelStatus::const_iterator chit = ecalStatus->find(detid);
627     // refer https://twiki.cern.ch/twiki/bin/viewauth/CMS/EcalChannelStatus
628 csander 1.11 int status = (chit != ecalStatus->end()) ? chit->getStatusCode() & 0x1F : -1;
629 csander 1.10 //std::cout << ieta << " " << iphi << " " << status << std:: endl;
630     const CaloSubdetectorGeometry* subGeom = geometry->getSubdetectorGeometry(detid);
631     const CaloCellGeometry* cellGeom = subGeom->getGeometry(detid);
632     double eta = cellGeom->getPosition().eta();
633     double phi = cellGeom->getPosition().phi();
634     double theta = cellGeom->getPosition().theta();
635    
636     if (status >= _maskedEcalChannelStatusThreshold) {
637     std::vector<double> valVec;
638     valVec.push_back(eta);
639     valVec.push_back(phi);
640     valVec.push_back(theta);
641     //std::cout << eta << " " << phi << std::endl;
642     EcalAllDeadChannelsValMap.insert(std::make_pair(detid, valVec));
643     }
644     } // end loop iphi
645     } // end loop ieta
646    
647     // Loop over EE detid
648     for (int ix = 0; ix <= 100; ix++) {
649     for (int iy = 0; iy <= 100; iy++) {
650     for (int iz = -1; iz <= 1; iz++) {
651     if (iz == 0)
652     continue;
653     if (!EEDetId::validDetId(ix, iy, iz))
654     continue;
655    
656     const EEDetId detid = EEDetId(ix, iy, iz, EEDetId::XYMODE);
657     EcalChannelStatus::const_iterator chit = ecalStatus->find(detid);
658 csander 1.11 int status = (chit != ecalStatus->end()) ? chit->getStatusCode() & 0x1F : -1;
659 csander 1.10
660     const CaloSubdetectorGeometry* subGeom = geometry->getSubdetectorGeometry(detid);
661     const CaloCellGeometry* cellGeom = subGeom->getGeometry(detid);
662     double eta = cellGeom->getPosition().eta();
663     double phi = cellGeom->getPosition().phi();
664 csander 1.3 double theta = cellGeom->getPosition().theta();
665    
666 csander 1.10 if (status >= _maskedEcalChannelStatusThreshold) {
667     std::vector<double> valVec;
668     valVec.push_back(eta);
669     valVec.push_back(phi);
670     valVec.push_back(theta);
671     //std::cout << eta << " " << phi << std::endl;
672     EcalAllDeadChannelsValMap.insert(std::make_pair(detid, valVec));
673 csander 1.3 }
674 csander 1.10 } // end loop iz
675     } // end loop iy
676     } // end loop ix
677 csander 1.3
678 csander 1.10 mapsReady = true;
679     return 1;
680 csander 1.3 }
681 csander 1.1
682 kaussen 1.16 void MCResolutions::ResizeHistoVector(std::vector<std::vector<TH1F*> > &histoVector) {
683    
684     histoVector.resize(EtaBinEdges.size() - 1);
685     for (std::vector<std::vector<TH1F*> >::iterator it = histoVector.begin(); it != histoVector.end(); ++it) {
686     it->resize(PtBinEdges.size() - 1);
687     }
688     }
689    
690 csander 1.1 //define this as a plug-in
691 csander 1.10 DEFINE_FWK_MODULE( MCResolutions);