ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/csander/JetResolutionFromMC/MCResolutions/src/MCResolutions.cc
Revision: 1.15
Committed: Wed Jul 25 09:44:33 2012 UTC (12 years, 9 months ago) by kaussen
Content type: text/plain
Branch: MAIN
Changes since 1.14: +175 -53 lines
Log Message: 
Jet multiplicity bins extended from 1,2,3plus to 1,2,3,4,5plus

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.15 // $Id: MCResolutions.cc,v 1.14 2011/08/08 16:35:55 csander 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     h_tot_DiJet_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
95     h_b_DiJet_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
96     h_nob_DiJet_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
97     h_dead_DiJet_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
98     h_deadb_DiJet_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
99     for (std::vector<std::vector<TH1F*> >::iterator it = h_tot_DiJet_JetResPt_Pt.begin(); it
100     != h_tot_DiJet_JetResPt_Pt.end(); ++it) {
101     it->resize(PtBinEdges.size() - 1);
102     }
103     for (std::vector<std::vector<TH1F*> >::iterator it = h_b_DiJet_JetResPt_Pt.begin(); it
104     != h_b_DiJet_JetResPt_Pt.end(); ++it) {
105     it->resize(PtBinEdges.size() - 1);
106     }
107     for (std::vector<std::vector<TH1F*> >::iterator it = h_nob_DiJet_JetResPt_Pt.begin(); it
108     != h_nob_DiJet_JetResPt_Pt.end(); ++it) {
109     it->resize(PtBinEdges.size() - 1);
110     }
111     for (std::vector<std::vector<TH1F*> >::iterator it = h_dead_DiJet_JetResPt_Pt.begin(); it
112     != h_dead_DiJet_JetResPt_Pt.end(); ++it) {
113     it->resize(PtBinEdges.size() - 1);
114     }
115     for (std::vector<std::vector<TH1F*> >::iterator it = h_deadb_DiJet_JetResPt_Pt.begin(); it
116     != h_deadb_DiJet_JetResPt_Pt.end(); ++it) {
117     it->resize(PtBinEdges.size() - 1);
118     }
119    
120     //// Array of histograms for jet resolutions (all jet multiplicities)
121 csander 1.13 h_tot_JetAll_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
122     h_b_JetAll_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
123     h_nob_JetAll_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
124     h_dead_JetAll_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
125     h_deadb_JetAll_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
126     for (std::vector<std::vector<TH1F*> >::iterator it = h_tot_JetAll_JetResPt_Pt.begin(); it
127     != h_tot_JetAll_JetResPt_Pt.end(); ++it) {
128     it->resize(PtBinEdges.size() - 1);
129     }
130     for (std::vector<std::vector<TH1F*> >::iterator it = h_b_JetAll_JetResPt_Pt.begin(); it
131     != h_b_JetAll_JetResPt_Pt.end(); ++it) {
132     it->resize(PtBinEdges.size() - 1);
133     }
134     for (std::vector<std::vector<TH1F*> >::iterator it = h_nob_JetAll_JetResPt_Pt.begin(); it
135     != h_nob_JetAll_JetResPt_Pt.end(); ++it) {
136     it->resize(PtBinEdges.size() - 1);
137     }
138     for (std::vector<std::vector<TH1F*> >::iterator it = h_dead_JetAll_JetResPt_Pt.begin(); it
139     != h_dead_JetAll_JetResPt_Pt.end(); ++it) {
140     it->resize(PtBinEdges.size() - 1);
141     }
142     for (std::vector<std::vector<TH1F*> >::iterator it = h_deadb_JetAll_JetResPt_Pt.begin(); it
143     != h_deadb_JetAll_JetResPt_Pt.end(); ++it) {
144     it->resize(PtBinEdges.size() - 1);
145     }
146    
147     //// Array of histograms for jet resolutions (first jet)
148     h_tot_Jet1_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
149     h_b_Jet1_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
150     h_nob_Jet1_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
151     h_dead_Jet1_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
152     h_deadb_Jet1_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
153     for (std::vector<std::vector<TH1F*> >::iterator it = h_tot_Jet1_JetResPt_Pt.begin(); it
154     != h_tot_Jet1_JetResPt_Pt.end(); ++it) {
155     it->resize(PtBinEdges.size() - 1);
156     }
157     for (std::vector<std::vector<TH1F*> >::iterator it = h_b_Jet1_JetResPt_Pt.begin(); it != h_b_Jet1_JetResPt_Pt.end(); ++it) {
158     it->resize(PtBinEdges.size() - 1);
159     }
160     for (std::vector<std::vector<TH1F*> >::iterator it = h_nob_Jet1_JetResPt_Pt.begin(); it
161     != h_nob_Jet1_JetResPt_Pt.end(); ++it) {
162     it->resize(PtBinEdges.size() - 1);
163     }
164     for (std::vector<std::vector<TH1F*> >::iterator it = h_dead_Jet1_JetResPt_Pt.begin(); it
165     != h_dead_Jet1_JetResPt_Pt.end(); ++it) {
166     it->resize(PtBinEdges.size() - 1);
167     }
168     for (std::vector<std::vector<TH1F*> >::iterator it = h_deadb_Jet1_JetResPt_Pt.begin(); it
169     != h_deadb_Jet1_JetResPt_Pt.end(); ++it) {
170     it->resize(PtBinEdges.size() - 1);
171     }
172    
173     //// Array of histograms for jet resolutions (second jet)
174     h_tot_Jet2_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
175     h_b_Jet2_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
176     h_nob_Jet2_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
177     h_dead_Jet2_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
178     h_deadb_Jet2_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
179     for (std::vector<std::vector<TH1F*> >::iterator it = h_tot_Jet2_JetResPt_Pt.begin(); it
180     != h_tot_Jet2_JetResPt_Pt.end(); ++it) {
181     it->resize(PtBinEdges.size() - 1);
182     }
183     for (std::vector<std::vector<TH1F*> >::iterator it = h_b_Jet2_JetResPt_Pt.begin(); it != h_b_Jet2_JetResPt_Pt.end(); ++it) {
184     it->resize(PtBinEdges.size() - 1);
185     }
186     for (std::vector<std::vector<TH1F*> >::iterator it = h_nob_Jet2_JetResPt_Pt.begin(); it
187     != h_nob_Jet2_JetResPt_Pt.end(); ++it) {
188     it->resize(PtBinEdges.size() - 1);
189     }
190     for (std::vector<std::vector<TH1F*> >::iterator it = h_dead_Jet2_JetResPt_Pt.begin(); it
191     != h_dead_Jet2_JetResPt_Pt.end(); ++it) {
192     it->resize(PtBinEdges.size() - 1);
193     }
194     for (std::vector<std::vector<TH1F*> >::iterator it = h_deadb_Jet2_JetResPt_Pt.begin(); it
195     != h_deadb_Jet2_JetResPt_Pt.end(); ++it) {
196     it->resize(PtBinEdges.size() - 1);
197     }
198    
199 kaussen 1.15 //// Array of histograms for jet resolutions (third jet)
200     h_tot_Jet3_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
201     h_b_Jet3_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
202     h_nob_Jet3_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
203     h_dead_Jet3_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
204     h_deadb_Jet3_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
205     for (std::vector<std::vector<TH1F*> >::iterator it = h_tot_Jet3_JetResPt_Pt.begin(); it
206     != h_tot_Jet3_JetResPt_Pt.end(); ++it) {
207 csander 1.13 it->resize(PtBinEdges.size() - 1);
208     }
209 kaussen 1.15 for (std::vector<std::vector<TH1F*> >::iterator it = h_b_Jet3_JetResPt_Pt.begin(); it
210     != h_b_Jet3_JetResPt_Pt.end(); ++it) {
211 csander 1.13 it->resize(PtBinEdges.size() - 1);
212     }
213 kaussen 1.15 for (std::vector<std::vector<TH1F*> >::iterator it = h_nob_Jet3_JetResPt_Pt.begin(); it
214     != h_nob_Jet3_JetResPt_Pt.end(); ++it) {
215 csander 1.13 it->resize(PtBinEdges.size() - 1);
216     }
217 kaussen 1.15 for (std::vector<std::vector<TH1F*> >::iterator it = h_dead_Jet3_JetResPt_Pt.begin(); it
218     != h_dead_Jet3_JetResPt_Pt.end(); ++it) {
219 csander 1.13 it->resize(PtBinEdges.size() - 1);
220     }
221 kaussen 1.15 for (std::vector<std::vector<TH1F*> >::iterator it = h_deadb_Jet3_JetResPt_Pt.begin(); it
222     != h_deadb_Jet3_JetResPt_Pt.end(); ++it) {
223     it->resize(PtBinEdges.size() - 1);
224     }
225    
226     //// Array of histograms for jet resolutions (fourth jet)
227     h_tot_Jet4_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
228     h_b_Jet4_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
229     h_nob_Jet4_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
230     h_dead_Jet4_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
231     h_deadb_Jet4_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
232     for (std::vector<std::vector<TH1F*> >::iterator it = h_tot_Jet4_JetResPt_Pt.begin(); it
233     != h_tot_Jet4_JetResPt_Pt.end(); ++it) {
234     it->resize(PtBinEdges.size() - 1);
235     }
236     for (std::vector<std::vector<TH1F*> >::iterator it = h_b_Jet4_JetResPt_Pt.begin(); it
237     != h_b_Jet4_JetResPt_Pt.end(); ++it) {
238     it->resize(PtBinEdges.size() - 1);
239     }
240     for (std::vector<std::vector<TH1F*> >::iterator it = h_nob_Jet4_JetResPt_Pt.begin(); it
241     != h_nob_Jet4_JetResPt_Pt.end(); ++it) {
242     it->resize(PtBinEdges.size() - 1);
243     }
244     for (std::vector<std::vector<TH1F*> >::iterator it = h_dead_Jet4_JetResPt_Pt.begin(); it
245     != h_dead_Jet4_JetResPt_Pt.end(); ++it) {
246     it->resize(PtBinEdges.size() - 1);
247     }
248     for (std::vector<std::vector<TH1F*> >::iterator it = h_deadb_Jet4_JetResPt_Pt.begin(); it
249     != h_deadb_Jet4_JetResPt_Pt.end(); ++it) {
250     it->resize(PtBinEdges.size() - 1);
251     }
252    
253     //// Array of histograms for jet resolutions (fifth+ jet)
254     h_tot_Jet5p_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
255     h_b_Jet5p_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
256     h_nob_Jet5p_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
257     h_dead_Jet5p_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
258     h_deadb_Jet5p_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
259     for (std::vector<std::vector<TH1F*> >::iterator it = h_tot_Jet5p_JetResPt_Pt.begin(); it
260     != h_tot_Jet5p_JetResPt_Pt.end(); ++it) {
261     it->resize(PtBinEdges.size() - 1);
262     }
263     for (std::vector<std::vector<TH1F*> >::iterator it = h_b_Jet5p_JetResPt_Pt.begin(); it
264     != h_b_Jet5p_JetResPt_Pt.end(); ++it) {
265     it->resize(PtBinEdges.size() - 1);
266     }
267     for (std::vector<std::vector<TH1F*> >::iterator it = h_nob_Jet5p_JetResPt_Pt.begin(); it
268     != h_nob_Jet5p_JetResPt_Pt.end(); ++it) {
269     it->resize(PtBinEdges.size() - 1);
270     }
271     for (std::vector<std::vector<TH1F*> >::iterator it = h_dead_Jet5p_JetResPt_Pt.begin(); it
272     != h_dead_Jet5p_JetResPt_Pt.end(); ++it) {
273     it->resize(PtBinEdges.size() - 1);
274     }
275     for (std::vector<std::vector<TH1F*> >::iterator it = h_deadb_Jet5p_JetResPt_Pt.begin(); it
276     != h_deadb_Jet5p_JetResPt_Pt.end(); ++it) {
277 csander 1.13 it->resize(PtBinEdges.size() - 1);
278     }
279 csander 1.1
280     }
281    
282 csander 1.3 MCResolutions::~MCResolutions() {
283 csander 1.2
284 csander 1.10 // do anything here that needs to be done at desctruction time
285     // (e.g. close files, deallocate resources etc.)
286 csander 1.1
287     }
288    
289     //
290     // member functions
291     //
292    
293     // ------------ method called to for each event ------------
294 csander 1.10 void MCResolutions::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
295     using namespace std;
296 csander 1.5
297 csander 1.10 // Event setup
298     envSet(iSetup);
299 csander 1.9
300 csander 1.10 getChannelStatusMaps();
301 csander 1.9
302 csander 1.10 //Weight
303     edm::Handle<double> event_weight;
304     bool findWeight = iEvent.getByLabel(_weightName, event_weight);
305     weight = (event_weight.isValid() ? (*event_weight) : 1.0);
306     if (!findWeight) {
307     cout << "Weight not found!" << endl;
308     }
309    
310     //GenJets
311     edm::Handle<edm::View<reco::GenJet> > Jets_gen;
312     iEvent.getByLabel(_genJetTag, Jets_gen);
313    
314     //RecoJets
315     edm::Handle<edm::View<pat::Jet> > Jets_rec;
316     iEvent.getByLabel(_jetTag, Jets_rec);
317    
318     edm::Handle<reco::GenParticleCollection> genParticles;
319     iEvent.getByLabel("genParticles", genParticles);
320    
321 kaussen 1.15 //// Do DiJet selection
322 csander 1.10 bool isDiJet = false;
323     const reco::GenJet* gj1 = 0;
324     const reco::GenJet* gj2 = 0;
325     const reco::GenJet* gj3 = 0;
326     int k = 0;
327     for (edm::View<reco::GenJet>::const_iterator it = Jets_gen->begin(); it != Jets_gen->end(); ++it) {
328     ++k;
329     if (k == 1)
330     gj1 = &(*it);
331     if (k == 2)
332     gj2 = &(*it);
333     if (k == 3) {
334     gj3 = &(*it);
335     break;
336     }
337     }
338     if (k > 1) {
339     if (deltaPhi(gj1->phi(), gj2->phi()) > _deltaPhiDiJet) {
340     if (k == 2)
341     isDiJet = true;
342     if (k == 3) {
343     //if (gj3->pt()<(_relCut3rdJet*(gj1->pt()+gj2->pt())/2)){
344     if (gj3->pt() < (_relCut3rdJet * gj1->pt())) {
345     isDiJet = true;
346 csander 1.4 }
347 csander 1.10 }
348     }
349     }
350    
351 csander 1.13 //// Calculate jet multiplicity
352 csander 1.10 double NGenJet = 0;
353     for (edm::View<reco::GenJet>::const_iterator it = Jets_gen-> begin(); it != Jets_gen-> end(); ++it) {
354     if (it->pt() > _jetMultPtCut && abs(it->eta()) < _jetMultEtaCut)
355     ++NGenJet;
356     }
357    
358     int JetNo = 0;
359     for (edm::View<reco::GenJet>::const_iterator it = Jets_gen->begin(); it != Jets_gen->end(); ++it) {
360     ++JetNo;
361    
362     if (it->pt() < _GenJetPtCut)
363     continue;
364 csander 1.13
365 csander 1.10 //// First look if there is no significant GenJet near the tested GenJet
366     double dRgenjet = 999.;
367     double PtdRmin = 0;
368     for (edm::View<reco::GenJet>::const_iterator kt = Jets_gen-> begin(); kt != Jets_gen->end(); ++kt) {
369     if (&(*it) != &(*kt))
370 csander 1.3 continue;
371 csander 1.10 double dR = deltaR(*it, *kt);
372     if (dR < dRgenjet) {
373     dRgenjet = dR;
374     PtdRmin = kt->pt();
375     }
376     }
377     if (dRgenjet < _deltaRMatchVeto && PtdRmin / it->pt() < _relPtVeto)
378     continue;
379     const pat::Jet* matchedJet = 0;
380     const pat::Jet* nearestJet = 0;
381     double dRmatched = 999.;
382     double dRnearest = 999.;
383     for (edm::View<pat::Jet>::const_iterator jt = Jets_rec-> begin(); jt != Jets_rec->end(); ++jt) {
384     double dR = deltaR(*it, *jt);
385     if (dR < dRmatched) {
386     nearestJet = matchedJet;
387     dRnearest = dRmatched;
388     matchedJet = &(*jt);
389     dRmatched = dR;
390     } else if (dR < dRnearest) {
391     nearestJet = &(*jt);
392     dRnearest = dR;
393     }
394     }
395 csander 1.13
396 csander 1.10 //// look if there is no further significant CaloJet near the genJet
397     if (dRmatched < _deltaRMatch && (nearestJet == 0 || dRnearest > _deltaRMatchVeto || (nearestJet->pt()
398     < _absPtVeto && nearestJet->pt() / matchedJet->pt() < _relPtVeto))) {
399     double res = matchedJet->pt() / it->pt();
400 csander 1.13 h_tot_JetAll_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
401 csander 1.10 if (isDiJet && JetNo < 3) {
402     h_tot_DiJet_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
403     }
404 csander 1.13 if (JetNo == 1) {
405     h_tot_Jet1_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
406     }
407     if (JetNo == 2) {
408     h_tot_Jet2_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
409     }
410 kaussen 1.15 if (JetNo == 3) {
411     h_tot_Jet3_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
412     }
413     if (JetNo == 4) {
414     h_tot_Jet4_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
415     }
416     if (JetNo > 4) {
417     h_tot_Jet5p_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
418 csander 1.13 }
419    
420     //// Use algorithmic matching for heavy flavour ID
421 csander 1.10 bool bTag = false;
422 csander 1.13 if (fabs(matchedJet->partonFlavour()) == 4 || fabs(matchedJet->partonFlavour()) == 5) {
423     bTag = true;
424 csander 1.10 }
425 csander 1.6
426 csander 1.10 if (bTag) {
427 csander 1.13 h_b_JetAll_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
428 csander 1.10 if (isDiJet && JetNo < 3) {
429     h_b_DiJet_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
430     }
431 csander 1.13 if (JetNo == 1) {
432     h_b_Jet1_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
433     }
434     if (JetNo == 2) {
435     h_b_Jet2_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
436     }
437 kaussen 1.15 if (JetNo == 3) {
438     h_b_Jet3_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
439     }
440     if (JetNo == 4) {
441     h_b_Jet4_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
442     }
443     if (JetNo > 4) {
444     h_b_Jet5p_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
445 csander 1.13 }
446 csander 1.10 } else {
447 csander 1.13 h_nob_JetAll_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
448 csander 1.10 if (isDiJet && JetNo < 3) {
449     h_nob_DiJet_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
450 csander 1.6 }
451 csander 1.13 if (JetNo == 1) {
452     h_nob_Jet1_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
453     }
454     if (JetNo == 2) {
455     h_nob_Jet2_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
456     }
457 kaussen 1.15 if (JetNo == 3) {
458     h_nob_Jet3_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
459     }
460     if (JetNo == 4) {
461     h_nob_Jet4_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
462     }
463     if (JetNo > 4) {
464     h_nob_Jet5p_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
465 csander 1.13 }
466 csander 1.10 }
467    
468     //// look if matched jet points into direction of masked ECAL cluster
469     for (std::map<DetId, std::vector<double> >::iterator kt = EcalAllDeadChannelsValMap.begin(); kt
470     != EcalAllDeadChannelsValMap.end(); ++kt) {
471     math::PtEtaPhiMLorentzVectorD Evec(100., kt->second.at(0), kt->second.at(1), 0.);
472     double dR_dead = deltaR(*matchedJet, Evec);
473     if (dR_dead < _deltaRDeadECal) {
474 csander 1.13 h_dead_JetAll_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
475 csander 1.10 if (isDiJet && JetNo < 3) {
476     h_dead_DiJet_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
477 csander 1.13 }
478     if (JetNo == 1) {
479     h_dead_Jet1_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
480     }
481     if (JetNo == 2) {
482     h_dead_Jet2_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
483     }
484 kaussen 1.15 if (JetNo == 3) {
485     h_dead_Jet3_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
486     }
487     if (JetNo == 4) {
488     h_dead_Jet4_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
489     }
490     if (JetNo > 4) {
491     h_dead_Jet5p_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
492 csander 1.13 }
493     if (bTag) {
494     h_deadb_JetAll_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
495     if (isDiJet && JetNo < 3) {
496 csander 1.10 h_deadb_DiJet_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
497 csander 1.13 }
498     if (JetNo == 1) {
499     h_deadb_Jet1_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
500     }
501     if (JetNo == 2) {
502     h_deadb_Jet2_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
503     }
504 kaussen 1.15 if (JetNo == 3) {
505     h_deadb_Jet3_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
506     }
507     if (JetNo == 4) {
508     h_deadb_Jet4_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
509     }
510     if (JetNo > 4) {
511     h_deadb_Jet5p_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
512 csander 1.10 }
513     }
514     break;
515 csander 1.3 }
516 csander 1.10 }
517     }
518     }
519 csander 1.1
520     }
521    
522     // ------------ method called once each job just before starting event loop ------------
523 csander 1.10 void MCResolutions::beginJob() {
524 csander 1.3
525 csander 1.10 for (unsigned int i_pt = 0; i_pt < PtBinEdges.size() - 1; ++i_pt) {
526     for (unsigned int i_eta = 0; i_eta < EtaBinEdges.size() - 1; ++i_eta) {
527     char hname[100];
528     //// Book histograms (all jet multiplicities)
529 csander 1.13 sprintf(hname, "h_tot_JetAll_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
530     h_tot_JetAll_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
531     h_tot_JetAll_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
532     sprintf(hname, "h_b_JetAll_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
533     h_b_JetAll_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
534     h_b_JetAll_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
535     sprintf(hname, "h_nob_JetAll_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
536     h_nob_JetAll_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
537     h_nob_JetAll_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
538     sprintf(hname, "h_dead_JetAll_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
539     h_dead_JetAll_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
540     h_dead_JetAll_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
541     sprintf(hname, "h_deadb_JetAll_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
542     h_deadb_JetAll_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
543     h_deadb_JetAll_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
544 csander 1.10 //// Book histograms (DiJets)
545     sprintf(hname, "h_tot_DiJet_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
546     h_tot_DiJet_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
547     h_tot_DiJet_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
548     sprintf(hname, "h_b_DiJet_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
549     h_b_DiJet_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
550     h_b_DiJet_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
551     sprintf(hname, "h_nob_DiJet_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
552     h_nob_DiJet_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
553     h_nob_DiJet_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
554     sprintf(hname, "h_dead_DiJet_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
555     h_dead_DiJet_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
556     h_dead_DiJet_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
557     sprintf(hname, "h_deadb_DiJet_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
558     h_deadb_DiJet_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
559     h_deadb_DiJet_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
560 csander 1.13 //// Book histograms (leading jet)
561     sprintf(hname, "h_tot_Jet1_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
562     h_tot_Jet1_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
563     h_tot_Jet1_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
564     sprintf(hname, "h_b_Jet1_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
565     h_b_Jet1_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
566     h_b_Jet1_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
567     sprintf(hname, "h_nob_Jet1_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
568     h_nob_Jet1_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
569     h_nob_Jet1_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
570     sprintf(hname, "h_dead_Jet1_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
571     h_dead_Jet1_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
572     h_dead_Jet1_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
573     sprintf(hname, "h_deadb_Jet1_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
574     h_deadb_Jet1_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
575     h_deadb_Jet1_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
576     //// Book histograms (second jet)
577     sprintf(hname, "h_tot_Jet2_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
578     h_tot_Jet2_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
579     h_tot_Jet2_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
580     sprintf(hname, "h_b_Jet2_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
581     h_b_Jet2_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
582     h_b_Jet2_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
583     sprintf(hname, "h_nob_Jet2_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
584     h_nob_Jet2_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
585     h_nob_Jet2_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
586     sprintf(hname, "h_dead_Jet2_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
587     h_dead_Jet2_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
588     h_dead_Jet2_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
589     sprintf(hname, "h_deadb_Jet2_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
590     h_deadb_Jet2_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
591     h_deadb_Jet2_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
592 kaussen 1.15 //// Book histograms (third jet)
593     sprintf(hname, "h_tot_Jet3_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
594     h_tot_Jet3_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
595     h_tot_Jet3_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
596     sprintf(hname, "h_b_Jet3_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
597     h_b_Jet3_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
598     h_b_Jet3_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
599     sprintf(hname, "h_nob_Jet3_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
600     h_nob_Jet3_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
601     h_nob_Jet3_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
602     sprintf(hname, "h_dead_Jet3_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
603     h_dead_Jet3_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
604     h_dead_Jet3_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
605     sprintf(hname, "h_deadb_Jet3_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
606     h_deadb_Jet3_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
607     h_deadb_Jet3_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
608     //// Book histograms (fourth jet)
609     sprintf(hname, "h_tot_Jet4_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
610     h_tot_Jet4_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
611     h_tot_Jet4_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
612     sprintf(hname, "h_b_Jet4_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
613     h_b_Jet4_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
614     h_b_Jet4_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
615     sprintf(hname, "h_nob_Jet4_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
616     h_nob_Jet4_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
617     h_nob_Jet4_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
618     sprintf(hname, "h_dead_Jet4_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
619     h_dead_Jet4_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
620     h_dead_Jet4_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
621     sprintf(hname, "h_deadb_Jet4_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
622     h_deadb_Jet4_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
623     h_deadb_Jet4_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
624     //// Book histograms (>= fifth jet)
625     sprintf(hname, "h_tot_Jet5p_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
626     h_tot_Jet5p_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
627     h_tot_Jet5p_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
628     sprintf(hname, "h_b_Jet5p_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
629     h_b_Jet5p_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
630     h_b_Jet5p_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
631     sprintf(hname, "h_nob_Jet5p_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
632     h_nob_Jet5p_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
633     h_nob_Jet5p_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
634     sprintf(hname, "h_dead_Jet5p_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
635     h_dead_Jet5p_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
636     h_dead_Jet5p_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
637     sprintf(hname, "h_deadb_Jet5p_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
638     h_deadb_Jet5p_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
639     h_deadb_Jet5p_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
640 csander 1.10 }
641     }
642    
643     mapsReady = false;
644 csander 1.3 }
645 csander 1.1
646     // ------------ method called once each job just after ending the event loop ------------
647 csander 1.10 void MCResolutions::endJob() {
648 csander 1.5
649 csander 1.10 hfile->cd();
650     //hfile->cd(_dirName.c_str());
651     // Save all objects in this file
652     for (unsigned int i_pt = 0; i_pt < PtBinEdges.size() - 1; ++i_pt) {
653     for (unsigned int i_eta = 0; i_eta < EtaBinEdges.size() - 1; ++i_eta) {
654     // total
655 csander 1.13 hfile->WriteTObject(h_tot_JetAll_JetResPt_Pt.at(i_eta).at(i_pt));
656 csander 1.10 hfile->WriteTObject(h_tot_DiJet_JetResPt_Pt.at(i_eta).at(i_pt));
657 csander 1.13 hfile->WriteTObject(h_tot_Jet1_JetResPt_Pt.at(i_eta).at(i_pt));
658     hfile->WriteTObject(h_tot_Jet2_JetResPt_Pt.at(i_eta).at(i_pt));
659 kaussen 1.15 hfile->WriteTObject(h_tot_Jet3_JetResPt_Pt.at(i_eta).at(i_pt));
660     hfile->WriteTObject(h_tot_Jet4_JetResPt_Pt.at(i_eta).at(i_pt));
661     hfile->WriteTObject(h_tot_Jet5p_JetResPt_Pt.at(i_eta).at(i_pt));
662 csander 1.10 // with btag
663 csander 1.13 hfile->WriteTObject(h_b_JetAll_JetResPt_Pt.at(i_eta).at(i_pt));
664 csander 1.10 hfile->WriteTObject(h_b_DiJet_JetResPt_Pt.at(i_eta).at(i_pt));
665 csander 1.13 hfile->WriteTObject(h_b_Jet1_JetResPt_Pt.at(i_eta).at(i_pt));
666     hfile->WriteTObject(h_b_Jet2_JetResPt_Pt.at(i_eta).at(i_pt));
667 kaussen 1.15 hfile->WriteTObject(h_b_Jet3_JetResPt_Pt.at(i_eta).at(i_pt));
668     hfile->WriteTObject(h_b_Jet4_JetResPt_Pt.at(i_eta).at(i_pt));
669     hfile->WriteTObject(h_b_Jet5p_JetResPt_Pt.at(i_eta).at(i_pt));
670 csander 1.10 // without btag
671 csander 1.13 hfile->WriteTObject(h_nob_JetAll_JetResPt_Pt.at(i_eta).at(i_pt));
672 csander 1.10 hfile->WriteTObject(h_nob_DiJet_JetResPt_Pt.at(i_eta).at(i_pt));
673 csander 1.13 hfile->WriteTObject(h_nob_Jet1_JetResPt_Pt.at(i_eta).at(i_pt));
674     hfile->WriteTObject(h_nob_Jet2_JetResPt_Pt.at(i_eta).at(i_pt));
675 kaussen 1.15 hfile->WriteTObject(h_nob_Jet3_JetResPt_Pt.at(i_eta).at(i_pt));
676     hfile->WriteTObject(h_nob_Jet4_JetResPt_Pt.at(i_eta).at(i_pt));
677     hfile->WriteTObject(h_nob_Jet5p_JetResPt_Pt.at(i_eta).at(i_pt));
678 csander 1.10 // in direction of dead ECAL cells
679 csander 1.13 hfile->WriteTObject(h_dead_JetAll_JetResPt_Pt.at(i_eta).at(i_pt));
680 csander 1.10 hfile->WriteTObject(h_dead_DiJet_JetResPt_Pt.at(i_eta).at(i_pt));
681 csander 1.13 hfile->WriteTObject(h_dead_Jet1_JetResPt_Pt.at(i_eta).at(i_pt));
682     hfile->WriteTObject(h_dead_Jet2_JetResPt_Pt.at(i_eta).at(i_pt));
683 kaussen 1.15 hfile->WriteTObject(h_dead_Jet3_JetResPt_Pt.at(i_eta).at(i_pt));
684     hfile->WriteTObject(h_dead_Jet4_JetResPt_Pt.at(i_eta).at(i_pt));
685     hfile->WriteTObject(h_dead_Jet5p_JetResPt_Pt.at(i_eta).at(i_pt));
686 csander 1.10 // in direction of dead ECAL cells and with b tag
687 csander 1.13 hfile->WriteTObject(h_deadb_JetAll_JetResPt_Pt.at(i_eta).at(i_pt));
688 csander 1.10 hfile->WriteTObject(h_deadb_DiJet_JetResPt_Pt.at(i_eta).at(i_pt));
689 csander 1.13 hfile->WriteTObject(h_deadb_Jet1_JetResPt_Pt.at(i_eta).at(i_pt));
690     hfile->WriteTObject(h_deadb_Jet2_JetResPt_Pt.at(i_eta).at(i_pt));
691 kaussen 1.15 hfile->WriteTObject(h_deadb_Jet3_JetResPt_Pt.at(i_eta).at(i_pt));
692     hfile->WriteTObject(h_deadb_Jet4_JetResPt_Pt.at(i_eta).at(i_pt));
693     hfile->WriteTObject(h_deadb_Jet5p_JetResPt_Pt.at(i_eta).at(i_pt));
694 csander 1.10 }
695     }
696    
697     hfile->cd();
698     hfile->WriteObject(&PtBinEdges, "PtBinEdges");
699     hfile->WriteObject(&EtaBinEdges, "EtaBinEdges");
700     //hfile->ls();
701 csander 1.5
702 csander 1.10 // Close the file.
703     hfile->Close();
704 csander 1.5
705     }
706 csander 1.3
707     int MCResolutions::PtBin(const double& pt) {
708 csander 1.10 int i_pt = -1;
709     for (std::vector<double>::const_iterator it = PtBinEdges.begin(); it != PtBinEdges.end(); ++it) {
710     if ((*it) > pt)
711     break;
712     ++i_pt;
713     }
714     if (i_pt < 0)
715     i_pt = 0;
716     if (i_pt > (int) PtBinEdges.size() - 2)
717     i_pt = (int) PtBinEdges.size() - 2;
718 csander 1.3
719 csander 1.10 return i_pt;
720 csander 1.3 }
721    
722     int MCResolutions::EtaBin(const double& eta) {
723 csander 1.10 int i_eta = -1;
724     for (std::vector<double>::const_iterator it = EtaBinEdges.begin(); it != EtaBinEdges.end(); ++it) {
725     if ((*it) > fabs(eta))
726     break;
727     ++i_eta;
728     }
729     if (i_eta < 0)
730     i_eta = 0;
731     if (i_eta > (int) EtaBinEdges.size() - 2)
732     i_eta = (int) EtaBinEdges.size() - 2;
733     return i_eta;
734 csander 1.3 }
735    
736     void MCResolutions::envSet(const edm::EventSetup& iSetup) {
737    
738 csander 1.10 ecalScale.setEventSetup(iSetup);
739 csander 1.3
740 csander 1.10 iSetup.get<EcalChannelStatusRcd> ().get(ecalStatus);
741     iSetup.get<CaloGeometryRecord> ().get(geometry);
742 csander 1.3
743 csander 1.10 if (!ecalStatus.isValid())
744     throw "Failed to get ECAL channel status!";
745     if (!geometry.isValid())
746     throw "Failed to get the geometry!";
747 csander 1.3
748     }
749    
750     int MCResolutions::getChannelStatusMaps() {
751    
752 csander 1.10 if (mapsReady)
753     return -1;
754 csander 1.3
755 csander 1.10 EcalAllDeadChannelsValMap.clear();
756 csander 1.3
757 csander 1.10 // Loop over EB ...
758     for (int ieta = -85; ieta <= 85; ieta++) {
759     for (int iphi = 0; iphi <= 360; iphi++) {
760     if (!EBDetId::validDetId(ieta, iphi))
761     continue;
762    
763     const EBDetId detid = EBDetId(ieta, iphi, EBDetId::ETAPHIMODE);
764     EcalChannelStatus::const_iterator chit = ecalStatus->find(detid);
765     // refer https://twiki.cern.ch/twiki/bin/viewauth/CMS/EcalChannelStatus
766 csander 1.11 int status = (chit != ecalStatus->end()) ? chit->getStatusCode() & 0x1F : -1;
767 csander 1.10 //std::cout << ieta << " " << iphi << " " << status << std:: endl;
768     const CaloSubdetectorGeometry* subGeom = geometry->getSubdetectorGeometry(detid);
769     const CaloCellGeometry* cellGeom = subGeom->getGeometry(detid);
770     double eta = cellGeom->getPosition().eta();
771     double phi = cellGeom->getPosition().phi();
772     double theta = cellGeom->getPosition().theta();
773    
774     if (status >= _maskedEcalChannelStatusThreshold) {
775     std::vector<double> valVec;
776     valVec.push_back(eta);
777     valVec.push_back(phi);
778     valVec.push_back(theta);
779     //std::cout << eta << " " << phi << std::endl;
780     EcalAllDeadChannelsValMap.insert(std::make_pair(detid, valVec));
781     }
782     } // end loop iphi
783     } // end loop ieta
784    
785     // Loop over EE detid
786     for (int ix = 0; ix <= 100; ix++) {
787     for (int iy = 0; iy <= 100; iy++) {
788     for (int iz = -1; iz <= 1; iz++) {
789     if (iz == 0)
790     continue;
791     if (!EEDetId::validDetId(ix, iy, iz))
792     continue;
793    
794     const EEDetId detid = EEDetId(ix, iy, iz, EEDetId::XYMODE);
795     EcalChannelStatus::const_iterator chit = ecalStatus->find(detid);
796 csander 1.11 int status = (chit != ecalStatus->end()) ? chit->getStatusCode() & 0x1F : -1;
797 csander 1.10
798     const CaloSubdetectorGeometry* subGeom = geometry->getSubdetectorGeometry(detid);
799     const CaloCellGeometry* cellGeom = subGeom->getGeometry(detid);
800     double eta = cellGeom->getPosition().eta();
801     double phi = cellGeom->getPosition().phi();
802 csander 1.3 double theta = cellGeom->getPosition().theta();
803    
804 csander 1.10 if (status >= _maskedEcalChannelStatusThreshold) {
805     std::vector<double> valVec;
806     valVec.push_back(eta);
807     valVec.push_back(phi);
808     valVec.push_back(theta);
809     //std::cout << eta << " " << phi << std::endl;
810     EcalAllDeadChannelsValMap.insert(std::make_pair(detid, valVec));
811 csander 1.3 }
812 csander 1.10 } // end loop iz
813     } // end loop iy
814     } // end loop ix
815 csander 1.3
816 csander 1.10 mapsReady = true;
817     return 1;
818 csander 1.3 }
819 csander 1.1
820     //define this as a plug-in
821 csander 1.10 DEFINE_FWK_MODULE( MCResolutions);