ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/csander/JetResolutionFromMC/MCResolutions/src/MCResolutions.cc
Revision: 1.5
Committed: Fri Oct 22 12:03:23 2010 UTC (14 years, 6 months ago) by csander
Content type: text/plain
Branch: MAIN
Changes since 1.4: +524 -54 lines
Log Message:
Added E, pt, and Eta vector to output

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     [Notes on implementation]
12     */
13     //
14     // Original Author: Christian Sander,,,
15     // Created: Wed Oct 6 18:22:21 CEST 2010
16 csander 1.5 // $Id: MCResolutions.cc,v 1.4 2010/10/11 18:01:08 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.2 //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 csander 1.3 _bTag = iConfig.getParameter<edm::InputTag> ("bTag");
32 csander 1.2 _EBRecHits = iConfig.getParameter<edm::InputTag> ("EBRecHits");
33     _EERecHits = iConfig.getParameter<edm::InputTag> ("EERecHits");
34 csander 1.5 _jetMultPtCut = iConfig.getParameter<double> ("jetMultPtCut");
35     _jetMultEtaCut = iConfig.getParameter<double> ("jetMultEtaCut");
36 csander 1.2 _deltaPhiDiJet = iConfig.getParameter<double> ("deltaPhiDiJet");
37     _absCut3rdJet = iConfig.getParameter<double> ("absCut3rdJet");
38     _relCut3rdJet = iConfig.getParameter<double> ("relCut3rdJet");
39     _deltaRMatch = iConfig.getParameter<double> ("deltaRMatch");
40     _deltaRMatchVeto = iConfig.getParameter<double> ("deltaRMatchVeto");
41     _absPtVeto = iConfig.getParameter<double> ("absPtVeto");
42     _relPtVeto = iConfig.getParameter<double> ("relPtVeto");
43     _deltaRDeadECal = iConfig.getParameter<double> ("deltaRDeadECal");
44 csander 1.3 _GenJetPtCut = iConfig.getParameter<double> ("GenJetPtCut");
45     _bTagCut = iConfig.getParameter<double> ("bTagCut");
46     _bTagDeltaR = iConfig.getParameter<double> ("bTagDelatR");
47     _maskedEcalChannelStatusThreshold = iConfig.getParameter<int> ("maskedEcalChannelStatusThreshold");
48 csander 1.5 _deadMethod = iConfig.getParameter<std::string> ("deadMethod");
49     _fileName = iConfig.getParameter<std::string> ("fileName");
50     _dirName = iConfig.getParameter<std::string> ("dirName");
51    
52     hfile = new TFile(_fileName.c_str(), "RECREATE", "Jet response in pT and eta bins");
53     hfile->mkdir(_dirName.c_str(), _dirName.c_str());
54 csander 1.3
55     EBinEdges.push_back(0);
56     EBinEdges.push_back(20);
57     EBinEdges.push_back(30);
58     EBinEdges.push_back(50);
59     EBinEdges.push_back(80);
60     EBinEdges.push_back(120);
61     EBinEdges.push_back(170);
62     EBinEdges.push_back(230);
63     EBinEdges.push_back(300);
64     EBinEdges.push_back(380);
65     EBinEdges.push_back(470);
66     EBinEdges.push_back(570);
67     EBinEdges.push_back(680);
68     EBinEdges.push_back(800);
69     EBinEdges.push_back(1000);
70     EBinEdges.push_back(1300);
71     EBinEdges.push_back(1700);
72     EBinEdges.push_back(2200);
73    
74     PtBinEdges.push_back(0);
75     PtBinEdges.push_back(20);
76     PtBinEdges.push_back(30);
77     PtBinEdges.push_back(50);
78     PtBinEdges.push_back(80);
79     PtBinEdges.push_back(120);
80     PtBinEdges.push_back(170);
81     PtBinEdges.push_back(230);
82     PtBinEdges.push_back(300);
83     PtBinEdges.push_back(380);
84     PtBinEdges.push_back(470);
85     PtBinEdges.push_back(570);
86     PtBinEdges.push_back(680);
87     PtBinEdges.push_back(800);
88     PtBinEdges.push_back(1000);
89     PtBinEdges.push_back(1300);
90     PtBinEdges.push_back(1700);
91     PtBinEdges.push_back(2200);
92    
93     EtaBinEdges.push_back(0.0);
94     EtaBinEdges.push_back(1.4);
95 csander 1.5 //EtaBinEdges.push_back(2.6);
96 csander 1.3 EtaBinEdges.push_back(3.2);
97 csander 1.5 //EtaBinEdges.push_back(4.5);
98 csander 1.3 EtaBinEdges.push_back(5.0);
99    
100 csander 1.5 //// Array of histograms for jet resolutions (all jet multiplicities)
101     h_tot_NJetAll_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
102     h_b_NJetAll_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
103     h_dead_NJetAll_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
104     h_deadb_NJetAll_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
105     for (std::vector<std::vector<TH1F*> >::iterator it = h_tot_NJetAll_JetResPt_Pt.begin(); it != h_tot_NJetAll_JetResPt_Pt.end(); ++it) {
106     it->resize(PtBinEdges.size() - 1);
107     }
108     for (std::vector<std::vector<TH1F*> >::iterator it = h_b_NJetAll_JetResPt_Pt.begin(); it != h_b_NJetAll_JetResPt_Pt.end(); ++it) {
109 csander 1.3 it->resize(PtBinEdges.size() - 1);
110     }
111 csander 1.5 for (std::vector<std::vector<TH1F*> >::iterator it = h_dead_NJetAll_JetResPt_Pt.begin(); it != h_dead_NJetAll_JetResPt_Pt.end(); ++it) {
112 csander 1.3 it->resize(PtBinEdges.size() - 1);
113     }
114 csander 1.5 for (std::vector<std::vector<TH1F*> >::iterator it = h_deadb_NJetAll_JetResPt_Pt.begin(); it != h_deadb_NJetAll_JetResPt_Pt.end(); ++it) {
115 csander 1.3 it->resize(PtBinEdges.size() - 1);
116     }
117 csander 1.5 h_tot_NJetAll_JetResPt_E.resize(EtaBinEdges.size() - 1);
118     h_b_NJetAll_JetResPt_E.resize(EtaBinEdges.size() - 1);
119     h_dead_NJetAll_JetResPt_E.resize(EtaBinEdges.size() - 1);
120     h_deadb_NJetAll_JetResPt_E.resize(EtaBinEdges.size() - 1);
121     for (std::vector<std::vector<TH1F*> >::iterator it = h_tot_NJetAll_JetResPt_E.begin(); it != h_tot_NJetAll_JetResPt_E.end(); ++it) {
122     it->resize(EBinEdges.size() - 1);
123     }
124     for (std::vector<std::vector<TH1F*> >::iterator it = h_b_NJetAll_JetResPt_E.begin(); it != h_b_NJetAll_JetResPt_E.end(); ++it) {
125     it->resize(EBinEdges.size() - 1);
126     }
127     for (std::vector<std::vector<TH1F*> >::iterator it = h_dead_NJetAll_JetResPt_E.begin(); it != h_dead_NJetAll_JetResPt_E.end(); ++it) {
128     it->resize(EBinEdges.size() - 1);
129     }
130     for (std::vector<std::vector<TH1F*> >::iterator it = h_deadb_NJetAll_JetResPt_E.begin(); it != h_deadb_NJetAll_JetResPt_E.end(); ++it) {
131     it->resize(EBinEdges.size() - 1);
132     }
133 csander 1.3
134 csander 1.5 //// Array of histograms for jet resolutions (jet multiplicities = 2)
135     h_tot_NJet2_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
136     h_b_NJet2_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
137     h_dead_NJet2_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
138     h_deadb_NJet2_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
139     for (std::vector<std::vector<TH1F*> >::iterator it = h_tot_NJet2_JetResPt_Pt.begin(); it != h_tot_NJet2_JetResPt_Pt.end(); ++it) {
140     it->resize(PtBinEdges.size() - 1);
141     }
142     for (std::vector<std::vector<TH1F*> >::iterator it = h_b_NJet2_JetResPt_Pt.begin(); it != h_b_NJet2_JetResPt_Pt.end(); ++it) {
143     it->resize(PtBinEdges.size() - 1);
144     }
145     for (std::vector<std::vector<TH1F*> >::iterator it = h_dead_NJet2_JetResPt_Pt.begin(); it != h_dead_NJet2_JetResPt_Pt.end(); ++it) {
146     it->resize(PtBinEdges.size() - 1);
147     }
148     for (std::vector<std::vector<TH1F*> >::iterator it = h_deadb_NJet2_JetResPt_Pt.begin(); it != h_deadb_NJet2_JetResPt_Pt.end(); ++it) {
149     it->resize(PtBinEdges.size() - 1);
150     }
151     h_tot_NJet2_JetResPt_E.resize(EtaBinEdges.size() - 1);
152     h_b_NJet2_JetResPt_E.resize(EtaBinEdges.size() - 1);
153     h_dead_NJet2_JetResPt_E.resize(EtaBinEdges.size() - 1);
154     h_deadb_NJet2_JetResPt_E.resize(EtaBinEdges.size() - 1);
155     for (std::vector<std::vector<TH1F*> >::iterator it = h_tot_NJet2_JetResPt_E.begin(); it != h_tot_NJet2_JetResPt_E.end(); ++it) {
156     it->resize(EBinEdges.size() - 1);
157     }
158     for (std::vector<std::vector<TH1F*> >::iterator it = h_b_NJet2_JetResPt_E.begin(); it != h_b_NJet2_JetResPt_E.end(); ++it) {
159     it->resize(EBinEdges.size() - 1);
160     }
161     for (std::vector<std::vector<TH1F*> >::iterator it = h_dead_NJet2_JetResPt_E.begin(); it != h_dead_NJet2_JetResPt_E.end(); ++it) {
162     it->resize(EBinEdges.size() - 1);
163     }
164     for (std::vector<std::vector<TH1F*> >::iterator it = h_deadb_NJet2_JetResPt_E.begin(); it != h_deadb_NJet2_JetResPt_E.end(); ++it) {
165     it->resize(EBinEdges.size() - 1);
166     }
167 csander 1.3
168 csander 1.5 //// Array of histograms for jet resolutions (jet multiplicities = 3)
169     h_tot_NJet3_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
170     h_b_NJet3_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
171     h_dead_NJet3_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
172     h_deadb_NJet3_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
173     for (std::vector<std::vector<TH1F*> >::iterator it = h_tot_NJet3_JetResPt_Pt.begin(); it != h_tot_NJet3_JetResPt_Pt.end(); ++it) {
174     it->resize(PtBinEdges.size() - 1);
175     }
176     for (std::vector<std::vector<TH1F*> >::iterator it = h_b_NJet3_JetResPt_Pt.begin(); it != h_b_NJet3_JetResPt_Pt.end(); ++it) {
177     it->resize(PtBinEdges.size() - 1);
178     }
179     for (std::vector<std::vector<TH1F*> >::iterator it = h_dead_NJet3_JetResPt_Pt.begin(); it != h_dead_NJet3_JetResPt_Pt.end(); ++it) {
180     it->resize(PtBinEdges.size() - 1);
181     }
182     for (std::vector<std::vector<TH1F*> >::iterator it = h_deadb_NJet3_JetResPt_Pt.begin(); it != h_deadb_NJet3_JetResPt_Pt.end(); ++it) {
183     it->resize(PtBinEdges.size() - 1);
184     }
185     h_tot_NJet3_JetResPt_E.resize(EtaBinEdges.size() - 1);
186     h_b_NJet3_JetResPt_E.resize(EtaBinEdges.size() - 1);
187     h_dead_NJet3_JetResPt_E.resize(EtaBinEdges.size() - 1);
188     h_deadb_NJet3_JetResPt_E.resize(EtaBinEdges.size() - 1);
189     for (std::vector<std::vector<TH1F*> >::iterator it = h_tot_NJet3_JetResPt_E.begin(); it != h_tot_NJet3_JetResPt_E.end(); ++it) {
190     it->resize(EBinEdges.size() - 1);
191     }
192     for (std::vector<std::vector<TH1F*> >::iterator it = h_b_NJet3_JetResPt_E.begin(); it != h_b_NJet3_JetResPt_E.end(); ++it) {
193 csander 1.3 it->resize(EBinEdges.size() - 1);
194     }
195 csander 1.5 for (std::vector<std::vector<TH1F*> >::iterator it = h_dead_NJet3_JetResPt_E.begin(); it != h_dead_NJet3_JetResPt_E.end(); ++it) {
196 csander 1.3 it->resize(EBinEdges.size() - 1);
197     }
198 csander 1.5 for (std::vector<std::vector<TH1F*> >::iterator it = h_deadb_NJet3_JetResPt_E.begin(); it != h_deadb_NJet3_JetResPt_E.end(); ++it) {
199     it->resize(EBinEdges.size() - 1);
200     }
201    
202     //// Array of histograms for jet resolutions (jet multiplicities = 4)
203     h_tot_NJet4_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
204     h_b_NJet4_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
205     h_dead_NJet4_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
206     h_deadb_NJet4_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
207     for (std::vector<std::vector<TH1F*> >::iterator it = h_tot_NJet4_JetResPt_Pt.begin(); it != h_tot_NJet4_JetResPt_Pt.end(); ++it) {
208     it->resize(PtBinEdges.size() - 1);
209     }
210     for (std::vector<std::vector<TH1F*> >::iterator it = h_b_NJet4_JetResPt_Pt.begin(); it != h_b_NJet4_JetResPt_Pt.end(); ++it) {
211     it->resize(PtBinEdges.size() - 1);
212     }
213     for (std::vector<std::vector<TH1F*> >::iterator it = h_dead_NJet4_JetResPt_Pt.begin(); it != h_dead_NJet4_JetResPt_Pt.end(); ++it) {
214     it->resize(PtBinEdges.size() - 1);
215     }
216     for (std::vector<std::vector<TH1F*> >::iterator it = h_deadb_NJet4_JetResPt_Pt.begin(); it != h_deadb_NJet4_JetResPt_Pt.end(); ++it) {
217     it->resize(PtBinEdges.size() - 1);
218     }
219     h_tot_NJet4_JetResPt_E.resize(EtaBinEdges.size() - 1);
220     h_b_NJet4_JetResPt_E.resize(EtaBinEdges.size() - 1);
221     h_dead_NJet4_JetResPt_E.resize(EtaBinEdges.size() - 1);
222     h_deadb_NJet4_JetResPt_E.resize(EtaBinEdges.size() - 1);
223     for (std::vector<std::vector<TH1F*> >::iterator it = h_tot_NJet4_JetResPt_E.begin(); it != h_tot_NJet4_JetResPt_E.end(); ++it) {
224     it->resize(EBinEdges.size() - 1);
225     }
226     for (std::vector<std::vector<TH1F*> >::iterator it = h_b_NJet4_JetResPt_E.begin(); it != h_b_NJet4_JetResPt_E.end(); ++it) {
227     it->resize(EBinEdges.size() - 1);
228     }
229     for (std::vector<std::vector<TH1F*> >::iterator it = h_dead_NJet4_JetResPt_E.begin(); it != h_dead_NJet4_JetResPt_E.end(); ++it) {
230     it->resize(EBinEdges.size() - 1);
231     }
232     for (std::vector<std::vector<TH1F*> >::iterator it = h_deadb_NJet4_JetResPt_E.begin(); it != h_deadb_NJet4_JetResPt_E.end(); ++it) {
233     it->resize(EBinEdges.size() - 1);
234     }
235    
236     //// Array of histograms for jet resolutions (jet multiplicities >= 5)
237     h_tot_NJet5p_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
238     h_b_NJet5p_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
239     h_dead_NJet5p_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
240     h_deadb_NJet5p_JetResPt_Pt.resize(EtaBinEdges.size() - 1);
241     for (std::vector<std::vector<TH1F*> >::iterator it = h_tot_NJet5p_JetResPt_Pt.begin(); it != h_tot_NJet5p_JetResPt_Pt.end(); ++it) {
242     it->resize(PtBinEdges.size() - 1);
243     }
244     for (std::vector<std::vector<TH1F*> >::iterator it = h_b_NJet5p_JetResPt_Pt.begin(); it != h_b_NJet5p_JetResPt_Pt.end(); ++it) {
245     it->resize(PtBinEdges.size() - 1);
246     }
247     for (std::vector<std::vector<TH1F*> >::iterator it = h_dead_NJet5p_JetResPt_Pt.begin(); it != h_dead_NJet5p_JetResPt_Pt.end(); ++it) {
248     it->resize(PtBinEdges.size() - 1);
249     }
250     for (std::vector<std::vector<TH1F*> >::iterator it = h_deadb_NJet5p_JetResPt_Pt.begin(); it != h_deadb_NJet5p_JetResPt_Pt.end(); ++it) {
251     it->resize(PtBinEdges.size() - 1);
252     }
253     h_tot_NJet5p_JetResPt_E.resize(EtaBinEdges.size() - 1);
254     h_b_NJet5p_JetResPt_E.resize(EtaBinEdges.size() - 1);
255     h_dead_NJet5p_JetResPt_E.resize(EtaBinEdges.size() - 1);
256     h_deadb_NJet5p_JetResPt_E.resize(EtaBinEdges.size() - 1);
257     for (std::vector<std::vector<TH1F*> >::iterator it = h_tot_NJet5p_JetResPt_E.begin(); it != h_tot_NJet5p_JetResPt_E.end(); ++it) {
258     it->resize(EBinEdges.size() - 1);
259     }
260     for (std::vector<std::vector<TH1F*> >::iterator it = h_b_NJet5p_JetResPt_E.begin(); it != h_b_NJet5p_JetResPt_E.end(); ++it) {
261     it->resize(EBinEdges.size() - 1);
262     }
263     for (std::vector<std::vector<TH1F*> >::iterator it = h_dead_NJet5p_JetResPt_E.begin(); it != h_dead_NJet5p_JetResPt_E.end(); ++it) {
264     it->resize(EBinEdges.size() - 1);
265     }
266     for (std::vector<std::vector<TH1F*> >::iterator it = h_deadb_NJet5p_JetResPt_E.begin(); it != h_deadb_NJet5p_JetResPt_E.end(); ++it) {
267 csander 1.3 it->resize(EBinEdges.size() - 1);
268     }
269 csander 1.1
270     }
271    
272 csander 1.3 MCResolutions::~MCResolutions() {
273 csander 1.2
274     // do anything here that needs to be done at desctruction time
275     // (e.g. close files, deallocate resources etc.)
276 csander 1.1
277     }
278    
279     //
280     // member functions
281     //
282    
283     // ------------ method called to for each event ------------
284     void
285 csander 1.3 MCResolutions::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
286     using namespace std;
287    
288     // Event setup
289     envSet(iSetup);
290    
291     getChannelStatusMaps();
292 csander 1.1
293 csander 1.2 //Weight
294     edm::Handle<double> event_weight;
295 csander 1.5 bool findWeight = iEvent.getByLabel(_weightName, event_weight);
296 csander 1.2 weight = (event_weight.isValid() ? (*event_weight) : 1.0);
297 csander 1.5 if (!findWeight) {
298     cout << "Weight not found!" << endl;
299     }
300 csander 1.2
301     //GenJets
302     edm::Handle<edm::View<reco::GenJet> > Jets_gen;
303 csander 1.3 iEvent.getByLabel(_genJetTag, Jets_gen);
304 csander 1.2
305     //PATJets
306     edm::Handle<edm::View<reco::CaloJet> > Jets_rec;
307 csander 1.3 iEvent.getByLabel(_jetTag, Jets_rec);
308    
309     //softMuonBJetTags
310     edm::Handle<reco::JetTagCollection> bTagHandle;
311     iEvent.getByLabel("trackCountingHighPurBJetTags", bTagHandle);
312     const reco::JetTagCollection & bTags = *(bTagHandle.product());
313    
314 csander 1.5 //Ullas dead Ecal tagger
315     edm::InputTag ecalAnomalousFilterTag("EcalAnomalousEventFilter","anomalousECALVariables","MCResolutions");
316     Handle<AnomalousECALVariables> anomalousECALvarsHandle;
317     AnomalousECALVariables anomalousECALvars;
318     if (_deadMethod == "Ulla") {
319     iEvent.getByLabel(ecalAnomalousFilterTag, anomalousECALvarsHandle);
320     if (anomalousECALvarsHandle.isValid()) {
321     anomalousECALvars = *anomalousECALvarsHandle;
322     } else {
323     cout << "anomalous ECAL Vars not valid/found" << endl;
324     }
325     }
326    
327     //// Get jet multiplicity
328     double NGenJet = 0;
329 csander 1.3 for (edm::View<reco::GenJet>::const_iterator it = Jets_gen->
330     begin();
331 csander 1.5 it != Jets_gen->
332     end();
333     ++it) {
334     if (it->pt()>_jetMultPtCut && abs(it->eta()) < _jetMultEtaCut)
335     ++NGenJet;
336     }
337    
338     for (edm::View<reco::GenJet>
339     ::const_iterator it = Jets_gen->
340     begin();
341 csander 1.3 it != Jets_gen->end();
342     ++it) {
343 csander 1.4 if (it->pt()
344     <_GenJetPtCut)
345     continue;
346     //// First look if there is no significant GenJet near the tested GenJet
347     double dRgenjet = 999.;
348     double PtdRmin = 0;
349     for (edm::View<reco::GenJet>::const_iterator kt = Jets_gen->
350     begin();
351     kt != Jets_gen->end();
352     ++kt) {
353     if (&(*it) != &(*kt))
354     continue;
355     double dR = deltaR (*it, *kt);
356     if (dR < dRgenjet ) {
357     dRgenjet = dR;
358     PtdRmin = kt->pt();
359     }
360     }
361     if (dRgenjet < _deltaRMatchVeto && PtdRmin/it->pt()
362     < _relPtVeto)
363 csander 1.3 continue;
364     const reco::CaloJet* matchedJet = 0;
365     const reco::CaloJet* nearestJet = 0;
366     double dRmatched = 999.;
367     double dRnearest = 999.;
368 csander 1.4 for (edm::View<reco::CaloJet>
369     ::const_iterator jt = Jets_rec->
370     begin();
371 csander 1.3 jt != Jets_rec->end();
372     ++jt) {
373     double dR = deltaR(*it, *jt);
374     if (dR < dRmatched) {
375     nearestJet = matchedJet;
376     dRnearest = dRmatched;
377     matchedJet = &(*jt);
378     dRmatched = dR;
379     } else if (dR < dRnearest) {
380     nearestJet = &(*jt);
381     dRnearest = dR;
382     }
383     }
384     //// look if there is no further significant CaloJet near the genJet
385     if (dRmatched < _deltaRMatch && ( nearestJet == 0 || dRnearest > _deltaRMatchVeto || (nearestJet->pt()
386     < _absPtVeto && nearestJet->pt()/matchedJet->pt() < _relPtVeto)
387     )
388     ) {
389     //cout << dRmatched << endl;
390     double res = matchedJet->pt() / it->pt();
391 csander 1.5 h_tot_NJetAll_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
392     h_tot_NJetAll_JetResPt_E.at(EtaBin(it->eta())).at(EBin(it->energy()))->Fill(res, weight);
393     if (NGenJet == 2) {
394     h_tot_NJet2_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
395     h_tot_NJet2_JetResPt_E.at(EtaBin(it->eta())).at(EBin(it->energy()))->Fill(res, weight);
396     }
397     if (NGenJet == 3) {
398     h_tot_NJet3_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
399     h_tot_NJet3_JetResPt_E.at(EtaBin(it->eta())).at(EBin(it->energy()))->Fill(res, weight);
400     }
401     if (NGenJet == 4) {
402     h_tot_NJet4_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
403     h_tot_NJet4_JetResPt_E.at(EtaBin(it->eta())).at(EBin(it->energy()))->Fill(res, weight);
404     }
405     if (NGenJet >= 5) {
406     h_tot_NJet5p_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
407     h_tot_NJet5p_JetResPt_E.at(EtaBin(it->eta())).at(EBin(it->energy()))->Fill(res, weight);
408     }
409 csander 1.3 //// look if matched jet can be matched to a b-tagged jet
410 csander 1.5 bool bTag = false;
411 csander 1.3 for (int i = 0; i != (int) bTags.size(); ++i) {
412     double dR_btag = deltaR(*matchedJet, *(bTags[i].first));
413     if (dR_btag < _bTagDeltaR && bTags[i].second > _bTagCut) {
414     //cout << "matched jet has b tag discriminator = "<<bTags[i].second << endl;
415 csander 1.5 h_b_NJetAll_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
416     h_b_NJetAll_JetResPt_E.at(EtaBin(it->eta())).at(EBin(it->energy()))->Fill(res, weight);
417     if (NGenJet == 2) {
418     h_b_NJet2_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
419     h_b_NJet2_JetResPt_E.at(EtaBin(it->eta())).at(EBin(it->energy()))->Fill(res, weight);
420     }
421     if (NGenJet == 3) {
422     h_b_NJet3_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
423     h_b_NJet3_JetResPt_E.at(EtaBin(it->eta())).at(EBin(it->energy()))->Fill(res, weight);
424     }
425     if (NGenJet == 4) {
426     h_b_NJet4_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
427     h_b_NJet4_JetResPt_E.at(EtaBin(it->eta())).at(EBin(it->energy()))->Fill(res, weight);
428     }
429     if (NGenJet >= 5) {
430     h_b_NJet5p_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
431     h_b_NJet5p_JetResPt_E.at(EtaBin(it->eta())).at(EBin(it->energy()))->Fill(res, weight);
432     }
433     bTag = true;
434 csander 1.3 break;
435     }
436     }
437 csander 1.5 if (_deadMethod == "dR") {
438     //// look if matched jet points into direction of masked ECAL cluster
439     for (std::map<DetId, std::vector<double> >::iterator kt = EcalAllDeadChannelsValMap.begin(); kt != EcalAllDeadChannelsValMap.end(); ++kt) {
440     math::PtEtaPhiMLorentzVectorD Evec(100., kt->second.at(0), kt->second.at(1), 0.);
441     double dR_dead = deltaR(*matchedJet, Evec);
442     if (dR_dead < _deltaRDeadECal) {
443     //cout << "jet is pointing to dead Ecal cluster (eta, phi)= "<< it->eta() << ", " << it->phi() << endl;
444     h_dead_NJetAll_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
445     h_dead_NJetAll_JetResPt_E.at(EtaBin(it->eta())).at(EBin(it->energy()))->Fill(res, weight);
446     if (bTag) {
447     h_deadb_NJetAll_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
448     h_deadb_NJetAll_JetResPt_E.at(EtaBin(it->eta())).at(EBin(it->energy()))->Fill(res, weight);
449     }
450     if (NGenJet == 2) {
451     h_dead_NJet2_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
452     h_dead_NJet2_JetResPt_E.at(EtaBin(it->eta())).at(EBin(it->energy()))->Fill(res, weight);
453     if (bTag) {
454     h_deadb_NJet2_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
455     h_deadb_NJet2_JetResPt_E.at(EtaBin(it->eta())).at(EBin(it->energy()))->Fill(res, weight);
456     }
457     }
458     if (NGenJet == 3) {
459     h_dead_NJet3_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
460     h_dead_NJet3_JetResPt_E.at(EtaBin(it->eta())).at(EBin(it->energy()))->Fill(res, weight);
461     if (bTag) {
462     h_deadb_NJet3_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
463     h_deadb_NJet3_JetResPt_E.at(EtaBin(it->eta())).at(EBin(it->energy()))->Fill(res, weight);
464     }
465     }
466     if (NGenJet == 4) {
467     h_dead_NJet4_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
468     h_dead_NJet4_JetResPt_E.at(EtaBin(it->eta())).at(EBin(it->energy()))->Fill(res, weight);
469     if (bTag) {
470     h_deadb_NJet4_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
471     h_deadb_NJet4_JetResPt_E.at(EtaBin(it->eta())).at(EBin(it->energy()))->Fill(res, weight);
472     }
473     }
474     if (NGenJet >= 5) {
475     h_dead_NJet5p_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
476     h_dead_NJet5p_JetResPt_E.at(EtaBin(it->eta())).at(EBin(it->energy()))->Fill(res, weight);
477     if (bTag) {
478     h_deadb_NJet5p_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
479     h_deadb_NJet5p_JetResPt_E.at(EtaBin(it->eta())).at(EBin(it->energy()))->Fill(res, weight);
480     }
481     }
482     break;
483     }
484     }
485     }
486     if (_deadMethod == "Boundary") {
487     if (anomalousECALvars.isEcalNoise()) {
488     h_dead_NJetAll_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
489     h_dead_NJetAll_JetResPt_E.at(EtaBin(it->eta())).at(EBin(it->energy()))->Fill(res, weight);
490     if (bTag) {
491     h_deadb_NJetAll_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
492     h_deadb_NJetAll_JetResPt_E.at(EtaBin(it->eta())).at(EBin(it->energy()))->Fill(res, weight);
493     }
494     if (NGenJet == 2) {
495     h_dead_NJet2_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
496     h_dead_NJet2_JetResPt_E.at(EtaBin(it->eta())).at(EBin(it->energy()))->Fill(res, weight);
497     if (bTag) {
498     h_deadb_NJet2_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
499     h_deadb_NJet2_JetResPt_E.at(EtaBin(it->eta())).at(EBin(it->energy()))->Fill(res, weight);
500     }
501     }
502     if (NGenJet == 3) {
503     h_dead_NJet3_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
504     h_dead_NJet3_JetResPt_E.at(EtaBin(it->eta())).at(EBin(it->energy()))->Fill(res, weight);
505     if (bTag) {
506     h_deadb_NJet3_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
507     h_deadb_NJet3_JetResPt_E.at(EtaBin(it->eta())).at(EBin(it->energy()))->Fill(res, weight);
508     }
509     }
510     if (NGenJet == 4) {
511     h_dead_NJet4_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
512     h_dead_NJet4_JetResPt_E.at(EtaBin(it->eta())).at(EBin(it->energy()))->Fill(res, weight);
513     if (bTag) {
514     h_deadb_NJet4_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
515     h_deadb_NJet4_JetResPt_E.at(EtaBin(it->eta())).at(EBin(it->energy()))->Fill(res, weight);
516     }
517     }
518     if (NGenJet >= 5) {
519     h_dead_NJet5p_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
520     h_dead_NJet5p_JetResPt_E.at(EtaBin(it->eta())).at(EBin(it->energy()))->Fill(res, weight);
521     if (bTag) {
522     h_deadb_NJet5p_JetResPt_Pt.at(EtaBin(it->eta())).at(PtBin(it->pt()))->Fill(res, weight);
523     h_deadb_NJet5p_JetResPt_E.at(EtaBin(it->eta())).at(EBin(it->energy()))->Fill(res, weight);
524     }
525     }
526 csander 1.3 }
527     }
528     }
529     }
530 csander 1.1
531     }
532    
533    
534     // ------------ method called once each job just before starting event loop ------------
535 csander 1.2 void
536 csander 1.3 MCResolutions::beginJob() {
537    
538     for (unsigned int i_pt = 0; i_pt < PtBinEdges.size() - 1; ++i_pt) {
539     for (unsigned int i_eta = 0; i_eta < EtaBinEdges.size() - 1; ++i_eta) {
540     char hname[100];
541 csander 1.5 //// Book histograms (all jet multiplicities)
542     sprintf(hname, "h_tot_NJetAll_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
543     h_tot_NJetAll_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
544     h_tot_NJetAll_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
545     //hfile->Append(h_tot_NJetAll_JetResPt_Pt.at(i_eta).at(i_pt));
546     sprintf(hname, "h_b_NJetAll_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
547     h_b_NJetAll_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
548     h_b_NJetAll_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
549     sprintf(hname, "h_dead_NJetAll_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
550     h_dead_NJetAll_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
551     h_dead_NJetAll_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
552     sprintf(hname, "h_deadb_NJetAll_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
553     h_deadb_NJetAll_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
554     h_deadb_NJetAll_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
555     //// Book histograms (jet multiplicities = 2)
556     sprintf(hname, "h_tot_NJet2_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
557     h_tot_NJet2_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
558     h_tot_NJet2_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
559     sprintf(hname, "h_b_NJet2_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
560     h_b_NJet2_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
561     h_b_NJet2_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
562     sprintf(hname, "h_dead_NJet2_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
563     h_dead_NJet2_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
564     h_dead_NJet2_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
565     sprintf(hname, "h_deadb_NJet2_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
566     h_deadb_NJet2_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
567     h_deadb_NJet2_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
568     //// Book histograms (jet multiplicities = 3)
569     sprintf(hname, "h_tot_NJet3_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
570     h_tot_NJet3_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
571     h_tot_NJet3_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
572     sprintf(hname, "h_b_NJet3_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
573     h_b_NJet3_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
574     h_b_NJet3_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
575     sprintf(hname, "h_dead_NJet3_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
576     h_dead_NJet3_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
577     h_dead_NJet3_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
578     sprintf(hname, "h_deadb_NJet3_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
579     h_deadb_NJet3_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
580     h_deadb_NJet3_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
581     //// Book histograms (jet multiplicities = 4)
582     sprintf(hname, "h_tot_NJet4_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
583     h_tot_NJet4_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
584     h_tot_NJet4_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
585     sprintf(hname, "h_b_NJet4_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
586     h_b_NJet4_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
587     h_b_NJet4_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
588     sprintf(hname, "h_dead_NJet4_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
589     h_dead_NJet4_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
590     h_dead_NJet4_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
591     sprintf(hname, "h_deadb_NJet4_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
592     h_deadb_NJet4_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
593     h_deadb_NJet4_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
594     //// Book histograms (jet multiplicities >= 5)
595     sprintf(hname, "h_tot_NJet5p_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
596     h_tot_NJet5p_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
597     h_tot_NJet5p_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
598     sprintf(hname, "h_b_NJet5p_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
599     h_b_NJet5p_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
600     h_b_NJet5p_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
601     sprintf(hname, "h_dead_NJet5p_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
602     h_dead_NJet5p_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
603     h_dead_NJet5p_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
604     sprintf(hname, "h_deadb_NJet5p_ResponsePt_Pt%i_Eta%i", i_pt, i_eta);
605     h_deadb_NJet5p_JetResPt_Pt.at(i_eta).at(i_pt) = new TH1F(hname, hname, 300, 0., 3.);
606     h_deadb_NJet5p_JetResPt_Pt.at(i_eta).at(i_pt)->Sumw2();
607 csander 1.3 }
608     }
609    
610     for (unsigned int i_e = 0; i_e < EBinEdges.size() - 1; ++i_e) {
611     for (unsigned int i_eta = 0; i_eta < EtaBinEdges.size() - 1; ++i_eta) {
612     char hname[100];
613 csander 1.5 //// Book histograms (all jet multiplicities)
614     sprintf(hname, "h_tot_NJetAll_ResponsePt_E%i_Eta%i", i_e, i_eta);
615     h_tot_NJetAll_JetResPt_E.at(i_eta).at(i_e) = new TH1F(hname, hname, 300, 0., 3.);
616     h_tot_NJetAll_JetResPt_E.at(i_eta).at(i_e)->Sumw2();
617     sprintf(hname, "h_b_NJetAll_ResponsePt_E%i_Eta%i", i_e, i_eta);
618     h_b_NJetAll_JetResPt_E.at(i_eta).at(i_e) = new TH1F(hname, hname, 300, 0., 3.);
619     h_b_NJetAll_JetResPt_E.at(i_eta).at(i_e)->Sumw2();
620     sprintf(hname, "h_dead_NJetAll_ResponsePt_E%i_Eta%i", i_e, i_eta);
621     h_dead_NJetAll_JetResPt_E.at(i_eta).at(i_e) = new TH1F(hname, hname, 300, 0., 3.);
622     h_dead_NJetAll_JetResPt_E.at(i_eta).at(i_e)->Sumw2();
623     sprintf(hname, "h_deadb_NJetAll_ResponsePt_E%i_Eta%i", i_e, i_eta);
624     h_deadb_NJetAll_JetResPt_E.at(i_eta).at(i_e) = new TH1F(hname, hname, 300, 0., 3.);
625     h_deadb_NJetAll_JetResPt_E.at(i_eta).at(i_e)->Sumw2();
626     //// Book histograms (jet multiplicities = 2)
627     sprintf(hname, "h_tot_NJet2_ResponsePt_E%i_Eta%i", i_e, i_eta);
628     h_tot_NJet2_JetResPt_E.at(i_eta).at(i_e) = new TH1F(hname, hname, 300, 0., 3.);
629     h_tot_NJet2_JetResPt_E.at(i_eta).at(i_e)->Sumw2();
630     sprintf(hname, "h_b_NJet2_ResponsePt_E%i_Eta%i", i_e, i_eta);
631     h_b_NJet2_JetResPt_E.at(i_eta).at(i_e) = new TH1F(hname, hname, 300, 0., 3.);
632     h_b_NJet2_JetResPt_E.at(i_eta).at(i_e)->Sumw2();
633     sprintf(hname, "h_dead_NJet2_ResponsePt_E%i_Eta%i", i_e, i_eta);
634     h_dead_NJet2_JetResPt_E.at(i_eta).at(i_e) = new TH1F(hname, hname, 300, 0., 3.);
635     h_dead_NJet2_JetResPt_E.at(i_eta).at(i_e)->Sumw2();
636     sprintf(hname, "h_deadb_NJet2_ResponsePt_E%i_Eta%i", i_e, i_eta);
637     h_deadb_NJet2_JetResPt_E.at(i_eta).at(i_e) = new TH1F(hname, hname, 300, 0., 3.);
638     h_deadb_NJet2_JetResPt_E.at(i_eta).at(i_e)->Sumw2();
639     //// Book histograms (jet multiplicities = 3)
640     sprintf(hname, "h_tot_NJet3_ResponsePt_E%i_Eta%i", i_e, i_eta);
641     h_tot_NJet3_JetResPt_E.at(i_eta).at(i_e) = new TH1F(hname, hname, 300, 0., 3.);
642     h_tot_NJet3_JetResPt_E.at(i_eta).at(i_e)->Sumw2();
643     sprintf(hname, "h_b_NJet3_ResponsePt_E%i_Eta%i", i_e, i_eta);
644     h_b_NJet3_JetResPt_E.at(i_eta).at(i_e) = new TH1F(hname, hname, 300, 0., 3.);
645     h_b_NJet3_JetResPt_E.at(i_eta).at(i_e)->Sumw2();
646     sprintf(hname, "h_dead_NJet3_ResponsePt_E%i_Eta%i", i_e, i_eta);
647     h_dead_NJet3_JetResPt_E.at(i_eta).at(i_e) = new TH1F(hname, hname, 300, 0., 3.);
648     h_dead_NJet3_JetResPt_E.at(i_eta).at(i_e)->Sumw2();
649     sprintf(hname, "h_deadb_NJet3_ResponsePt_E%i_Eta%i", i_e, i_eta);
650     h_deadb_NJet3_JetResPt_E.at(i_eta).at(i_e) = new TH1F(hname, hname, 300, 0., 3.);
651     h_deadb_NJet3_JetResPt_E.at(i_eta).at(i_e)->Sumw2();
652     //// Book histograms (jet multiplicities = 4)
653     sprintf(hname, "h_tot_NJet4_ResponsePt_E%i_Eta%i", i_e, i_eta);
654     h_tot_NJet4_JetResPt_E.at(i_eta).at(i_e) = new TH1F(hname, hname, 300, 0., 3.);
655     h_tot_NJet4_JetResPt_E.at(i_eta).at(i_e)->Sumw2();
656     sprintf(hname, "h_b_NJet4_ResponsePt_E%i_Eta%i", i_e, i_eta);
657     h_b_NJet4_JetResPt_E.at(i_eta).at(i_e) = new TH1F(hname, hname, 300, 0., 3.);
658     h_b_NJet4_JetResPt_E.at(i_eta).at(i_e)->Sumw2();
659     sprintf(hname, "h_dead_NJet4_ResponsePt_E%i_Eta%i", i_e, i_eta);
660     h_dead_NJet4_JetResPt_E.at(i_eta).at(i_e) = new TH1F(hname, hname, 300, 0., 3.);
661     h_dead_NJet4_JetResPt_E.at(i_eta).at(i_e)->Sumw2();
662     sprintf(hname, "h_deadb_NJet4_ResponsePt_E%i_Eta%i", i_e, i_eta);
663     h_deadb_NJet4_JetResPt_E.at(i_eta).at(i_e) = new TH1F(hname, hname, 300, 0., 3.);
664     h_deadb_NJet4_JetResPt_E.at(i_eta).at(i_e)->Sumw2();
665     //// Book histograms (jet multiplicities = 5p)
666     sprintf(hname, "h_tot_NJet5p_ResponsePt_E%i_Eta%i", i_e, i_eta);
667     h_tot_NJet5p_JetResPt_E.at(i_eta).at(i_e) = new TH1F(hname, hname, 300, 0., 3.);
668     h_tot_NJet5p_JetResPt_E.at(i_eta).at(i_e)->Sumw2();
669     sprintf(hname, "h_b_NJet5p_ResponsePt_E%i_Eta%i", i_e, i_eta);
670     h_b_NJet5p_JetResPt_E.at(i_eta).at(i_e) = new TH1F(hname, hname, 300, 0., 3.);
671     h_b_NJet5p_JetResPt_E.at(i_eta).at(i_e)->Sumw2();
672     sprintf(hname, "h_dead_NJet5p_ResponsePt_E%i_Eta%i", i_e, i_eta);
673     h_dead_NJet5p_JetResPt_E.at(i_eta).at(i_e) = new TH1F(hname, hname, 300, 0., 3.);
674     h_dead_NJet5p_JetResPt_E.at(i_eta).at(i_e)->Sumw2();
675     sprintf(hname, "h_deadb_NJet5p_ResponsePt_E%i_Eta%i", i_e, i_eta);
676     h_deadb_NJet5p_JetResPt_E.at(i_eta).at(i_e) = new TH1F(hname, hname, 300, 0., 3.);
677     h_deadb_NJet5p_JetResPt_E.at(i_eta).at(i_e)->Sumw2();
678 csander 1.3 }
679     }
680    
681     mapsReady = false;
682     }
683 csander 1.1
684     // ------------ method called once each job just after ending the event loop ------------
685 csander 1.2 void
686 csander 1.5 MCResolutions::endJob() {
687    
688     hfile->cd();
689     hfile->cd(_dirName.c_str());
690     // Save all objects in this file
691     for (unsigned int i_pt = 0; i_pt < PtBinEdges.size() - 1; ++i_pt) {
692     for (unsigned int i_eta = 0; i_eta < EtaBinEdges.size() - 1; ++i_eta) {
693     // total
694     h_tot_NJetAll_JetResPt_Pt.at(i_eta).at(i_pt)->Write();
695     h_tot_NJet2_JetResPt_Pt.at(i_eta).at(i_pt)->Write();
696     h_tot_NJet3_JetResPt_Pt.at(i_eta).at(i_pt)->Write();
697     h_tot_NJet4_JetResPt_Pt.at(i_eta).at(i_pt)->Write();
698     h_tot_NJet5p_JetResPt_Pt.at(i_eta).at(i_pt)->Write();
699     // with btag
700     h_b_NJetAll_JetResPt_Pt.at(i_eta).at(i_pt)->Write();
701     h_b_NJet2_JetResPt_Pt.at(i_eta).at(i_pt)->Write();
702     h_b_NJet3_JetResPt_Pt.at(i_eta).at(i_pt)->Write();
703     h_b_NJet4_JetResPt_Pt.at(i_eta).at(i_pt)->Write();
704     h_b_NJet5p_JetResPt_Pt.at(i_eta).at(i_pt)->Write();
705     // in direction of dead ECAL cells
706     h_dead_NJetAll_JetResPt_Pt.at(i_eta).at(i_pt)->Write();
707     h_dead_NJet2_JetResPt_Pt.at(i_eta).at(i_pt)->Write();
708     h_dead_NJet3_JetResPt_Pt.at(i_eta).at(i_pt)->Write();
709     h_dead_NJet4_JetResPt_Pt.at(i_eta).at(i_pt)->Write();
710     h_dead_NJet5p_JetResPt_Pt.at(i_eta).at(i_pt)->Write();
711     // in direction of dead ECAL cells and with b tag
712     h_deadb_NJetAll_JetResPt_Pt.at(i_eta).at(i_pt)->Write();
713     h_deadb_NJet2_JetResPt_Pt.at(i_eta).at(i_pt)->Write();
714     h_deadb_NJet3_JetResPt_Pt.at(i_eta).at(i_pt)->Write();
715     h_deadb_NJet4_JetResPt_Pt.at(i_eta).at(i_pt)->Write();
716     h_deadb_NJet5p_JetResPt_Pt.at(i_eta).at(i_pt)->Write();
717     }
718     }
719     for (unsigned int i_e = 0; i_e < EBinEdges.size() - 1; ++i_e) {
720     for (unsigned int i_eta = 0; i_eta < EtaBinEdges.size() - 1; ++i_eta) {
721     // total
722     h_tot_NJetAll_JetResPt_E.at(i_eta).at(i_e)->Write();
723     h_tot_NJet2_JetResPt_E.at(i_eta).at(i_e)->Write();
724     h_tot_NJet3_JetResPt_E.at(i_eta).at(i_e)->Write();
725     h_tot_NJet4_JetResPt_E.at(i_eta).at(i_e)->Write();
726     h_tot_NJet5p_JetResPt_E.at(i_eta).at(i_e)->Write();
727     // with btag
728     h_b_NJetAll_JetResPt_E.at(i_eta).at(i_e)->Write();
729     h_b_NJet2_JetResPt_E.at(i_eta).at(i_e)->Write();
730     h_b_NJet3_JetResPt_E.at(i_eta).at(i_e)->Write();
731     h_b_NJet4_JetResPt_E.at(i_eta).at(i_e)->Write();
732     h_b_NJet5p_JetResPt_E.at(i_eta).at(i_e)->Write();
733     // in direction of dead ECAL cells
734     h_dead_NJetAll_JetResPt_E.at(i_eta).at(i_e)->Write();
735     h_dead_NJet2_JetResPt_E.at(i_eta).at(i_e)->Write();
736     h_dead_NJet3_JetResPt_E.at(i_eta).at(i_e)->Write();
737     h_dead_NJet4_JetResPt_E.at(i_eta).at(i_e)->Write();
738     h_dead_NJet5p_JetResPt_E.at(i_eta).at(i_e)->Write();
739     // in direction of dead ECAL cells and with b tag
740     h_deadb_NJetAll_JetResPt_E.at(i_eta).at(i_e)->Write();
741     h_deadb_NJet2_JetResPt_E.at(i_eta).at(i_e)->Write();
742     h_deadb_NJet3_JetResPt_E.at(i_eta).at(i_e)->Write();
743     h_deadb_NJet4_JetResPt_E.at(i_eta).at(i_e)->Write();
744     h_deadb_NJet5p_JetResPt_E.at(i_eta).at(i_e)->Write();
745     }
746     }
747    
748     hfile->cd();
749     hfile->WriteObject(&EBinEdges, "EBinEdges");
750     hfile->WriteObject(&PtBinEdges, "PtBinEdges");
751     hfile->WriteObject(&EtaBinEdges, "EtaBinEdges");
752    
753     // Close the file.
754     hfile->Close();
755    
756     }
757 csander 1.3
758     int MCResolutions::EBin(const double& e) {
759     int i_e = -1;
760     for (std::vector<double>::const_iterator it = EBinEdges.begin(); it != EBinEdges.end(); ++it) {
761     if ((*it) > e)
762     break;
763     ++i_e;
764     }
765     if (i_e < 0)
766     i_e = 0;
767     if (i_e > (int)EBinEdges.size() - 2)
768     i_e = (int)EBinEdges.size() - 2;
769    
770     return i_e;
771     }
772    
773     int MCResolutions::PtBin(const double& pt) {
774     int i_pt = -1;
775     for (std::vector<double>::const_iterator it = PtBinEdges.begin(); it != PtBinEdges.end(); ++it) {
776     if ((*it) > pt)
777     break;
778     ++i_pt;
779     }
780     if (i_pt < 0)
781     i_pt = 0;
782     if (i_pt > (int)PtBinEdges.size() - 2)
783     i_pt = (int)PtBinEdges.size() - 2;
784    
785     return i_pt;
786     }
787    
788     int MCResolutions::EtaBin(const double& eta) {
789     int i_eta = -1;
790     for (std::vector<double>::const_iterator it = EtaBinEdges.begin(); it != EtaBinEdges.end(); ++it) {
791     if ((*it) > fabs(eta))
792     break;
793     ++i_eta;
794     }
795     if (i_eta < 0)
796     i_eta = 0;
797     if (i_eta > (int)EtaBinEdges.size() - 2)
798     i_eta = (int)EtaBinEdges.size() - 2;
799     return i_eta;
800     }
801    
802     void MCResolutions::envSet(const edm::EventSetup& iSetup) {
803    
804     ecalScale.setEventSetup( iSetup );
805    
806     iSetup.get<EcalChannelStatusRcd> ().get(ecalStatus);
807     iSetup.get<CaloGeometryRecord> ().get(geometry);
808    
809     if ( !ecalStatus.isValid() )
810     throw "Failed to get ECAL channel status!";
811     if ( !geometry.isValid() )
812     throw "Failed to get the geometry!";
813    
814     }
815    
816     int MCResolutions::getChannelStatusMaps() {
817    
818     if( mapsReady )
819     return -1;
820    
821     EcalAllDeadChannelsValMap.clear();
822    
823     // Loop over EB ...
824     for( int ieta=-85; ieta<=85; ieta++ ) {
825     for( int iphi=0; iphi<=360; iphi++ ) {
826     if(! EBDetId::validDetId( ieta, iphi ) )
827     continue;
828    
829     const EBDetId detid = EBDetId( ieta, iphi, EBDetId::ETAPHIMODE );
830     EcalChannelStatus::const_iterator chit = ecalStatus->find( detid );
831     // refer https://twiki.cern.ch/twiki/bin/viewauth/CMS/EcalChannelStatus
832     int status = ( chit != ecalStatus->end() ) ? chit->getStatusCode() : -1;
833     //std::cout << ieta << " " << iphi << " " << status << std:: endl;
834     const CaloSubdetectorGeometry* subGeom = geometry->getSubdetectorGeometry (detid);
835     const CaloCellGeometry* cellGeom = subGeom->getGeometry (detid);
836     double eta = cellGeom->getPosition ().eta ();
837     double phi = cellGeom->getPosition ().phi ();
838     double theta = cellGeom->getPosition().theta();
839    
840     if(status >= _maskedEcalChannelStatusThreshold) {
841     std::vector<double> valVec;
842     valVec.push_back(eta);
843     valVec.push_back(phi);
844     valVec.push_back(theta);
845     //std::cout << eta << " " << phi << std::endl;
846     EcalAllDeadChannelsValMap.insert( std::make_pair(detid, valVec) );
847     }
848     } // end loop iphi
849     } // end loop ieta
850    
851     // Loop over EE detid
852     for( int ix=0; ix<=100; ix++ ) {
853     for( int iy=0; iy<=100; iy++ ) {
854     for( int iz=-1; iz<=1; iz++ ) {
855     if(iz==0)
856     continue;
857     if(! EEDetId::validDetId( ix, iy, iz ) )
858     continue;
859    
860     const EEDetId detid = EEDetId( ix, iy, iz, EEDetId::XYMODE );
861     EcalChannelStatus::const_iterator chit = ecalStatus->find( detid );
862     int status = ( chit != ecalStatus->end() ) ? chit->getStatusCode() : -1;
863    
864     const CaloSubdetectorGeometry* subGeom = geometry->getSubdetectorGeometry (detid);
865     const CaloCellGeometry* cellGeom = subGeom->getGeometry (detid);
866     double eta = cellGeom->getPosition ().eta () ;
867     double phi = cellGeom->getPosition ().phi () ;
868     double theta = cellGeom->getPosition().theta();
869    
870     if(status >= _maskedEcalChannelStatusThreshold) {
871     std::vector<double> valVec;
872     valVec.push_back(eta);
873     valVec.push_back(phi);
874     valVec.push_back(theta);
875     //std::cout << eta << " " << phi << std::endl;
876     EcalAllDeadChannelsValMap.insert( std::make_pair(detid, valVec) );
877     }
878     } // end loop iz
879     } // end loop iy
880     } // end loop ix
881    
882     mapsReady = true;
883     return 1;
884     }
885 csander 1.1
886     //define this as a plug-in
887     DEFINE_FWK_MODULE(MCResolutions);