CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
EGEnergyCorrector.cc
Go to the documentation of this file.
1 // $Id: EGEnergyCorrector.cc,v 1.11 2012/06/17 20:19:52 bendavid Exp $
2 
3 #include <TFile.h>
4 #include "../interface/EGEnergyCorrector.h"
5 #include "CondFormats/EgammaObjects/interface/GBRForest.h"
6 #include "CondFormats/DataRecord/interface/GBRWrapperRcd.h"
7 #include "FWCore/Framework/interface/ESHandle.h"
8 #include "DataFormats/EgammaCandidates/interface/Photon.h"
9 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
10 #include "DataFormats/Math/interface/deltaPhi.h"
11 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterLazyTools.h"
12 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
13 #include "DataFormats/VertexReco/interface/Vertex.h"
14 
15 using namespace reco;
16 
17 //--------------------------------------------------------------------------------------------------
18 EGEnergyCorrector::EGEnergyCorrector() :
19 fReadereb(0),
20 fReaderebvariance(0),
21 fReaderee(0),
22 fReadereevariance(0),
23 fIsInitialized(kFALSE),
24 fOwnsForests(kFALSE),
25 fVals(0)
26 {
27  // Constructor.
28 }
29 
30 
31 //--------------------------------------------------------------------------------------------------
32 EGEnergyCorrector::~EGEnergyCorrector()
33 {
34 
35  if (fVals) delete [] fVals;
36 
37  if (fOwnsForests) {
38  if (fReadereb) delete fReadereb;
39  if (fReaderebvariance) delete fReaderebvariance;
40  if (fReaderee) delete fReaderee;
41  if (fReadereevariance) delete fReadereevariance;
42  }
43 
44 }
45 
46 //--------------------------------------------------------------------------------------------------
47 void EGEnergyCorrector::Initialize(const edm::EventSetup &iSetup, std::string regweights, bool weightsFromDB) {
48  fIsInitialized = kTRUE;
49 
50  if (fVals) delete [] fVals;
51  if (fOwnsForests) {
52  if (fReadereb) delete fReadereb;
53  if (fReaderebvariance) delete fReaderebvariance;
54  if (fReaderee) delete fReaderee;
55  if (fReadereevariance) delete fReadereevariance;
56  }
57 
58  fVals = new Float_t[73];
59 
60  if (weightsFromDB) { //weights from event setup
61 
62  edm::ESHandle<GBRForest> readereb;
63  edm::ESHandle<GBRForest> readerebvar;
64  edm::ESHandle<GBRForest> readeree;
65  edm::ESHandle<GBRForest> readereevar;
66 
67  iSetup.get<GBRWrapperRcd>().get(std::string(TString::Format("%s_EBCorrection",regweights.c_str())),readereb);
68  iSetup.get<GBRWrapperRcd>().get(std::string(TString::Format("%s_EBUncertainty",regweights.c_str())),readerebvar);
69  iSetup.get<GBRWrapperRcd>().get(std::string(TString::Format("%s_EECorrection",regweights.c_str())),readeree);
70  iSetup.get<GBRWrapperRcd>().get(std::string(TString::Format("%s_EEUncertainty",regweights.c_str())),readereevar);
71 
72  fReadereb = readereb.product();
73  fReaderebvariance = readerebvar.product();
74  fReaderee = readeree.product();
75  fReadereevariance = readereevar.product();
76 
77  }
78  else { //weights from root file
79  fOwnsForests = kTRUE;
80 
81  TFile *fgbr = TFile::Open(regweights.c_str(),"READ");
82  fReadereb = (GBRForest*)fgbr->Get("EBCorrection");
83  fReaderebvariance = (GBRForest*)fgbr->Get("EBUncertainty");
84  fReaderee = (GBRForest*)fgbr->Get("EECorrection");
85  fReadereevariance = (GBRForest*)fgbr->Get("EEUncertainty");
86  fgbr->Close();
87  }
88 
89 }
90 
91 //--------------------------------------------------------------------------------------------------
92 std::pair<double,double> EGEnergyCorrector::CorrectedEnergyWithError(const Photon &p, const reco::VertexCollection& vtxcol, EcalClusterLazyTools &clustertools, const edm::EventSetup &es) {
93 
94  const SuperClusterRef s = p.superCluster();
95  const CaloClusterPtr b = s->seed(); //seed basic cluster
96 
97  //highest energy basic cluster excluding seed basic cluster
98  CaloClusterPtr b2;
99  Double_t ebcmax = -99.;
100  for (reco::CaloCluster_iterator bit = s->clustersBegin(); bit!=s->clustersEnd(); ++bit) {
101  const CaloClusterPtr bc = *bit;
102  if (bc->energy() > ebcmax && bc !=b) {
103  b2 = bc;
104  ebcmax = bc->energy();
105  }
106  }
107 
108  //lowest energy basic cluster excluding seed (for pileup mitigation)
109  CaloClusterPtr bclast;
110  Double_t ebcmin = 1e6;
111  for (reco::CaloCluster_iterator bit = s->clustersBegin(); bit!=s->clustersEnd(); ++bit) {
112  const CaloClusterPtr bc = *bit;
113  if (bc->energy() < ebcmin && bc !=b) {
114  bclast = bc;
115  ebcmin = bc->energy();
116  }
117  }
118 
119  //2nd lowest energy basic cluster excluding seed (for pileup mitigation)
120  CaloClusterPtr bclast2;
121  ebcmin = 1e6;
122  for (reco::CaloCluster_iterator bit = s->clustersBegin(); bit!=s->clustersEnd(); ++bit) {
123  const CaloClusterPtr bc = *bit;
124  if (bc->energy() < ebcmin && bc !=b && bc!=bclast) {
125  bclast2 = bc;
126  ebcmin = bc->energy();
127  }
128  }
129 
130  Bool_t isbarrel = b->hitsAndFractions().at(0).first.subdetId()==EcalBarrel;
131  Bool_t hasbc2 = b2.isNonnull() && b2->energy()>0.;
132  Bool_t hasbclast = bclast.isNonnull() && bclast->energy()>0.;
133  Bool_t hasbclast2 = bclast2.isNonnull() && bclast2->energy()>0.;
134 
135 
136  if (isbarrel) {
137 
138  //basic supercluster variables
139  fVals[0] = s->rawEnergy();
140  fVals[1] = p.r9();
141  fVals[2] = s->eta();
142  fVals[3] = s->phi();
143  fVals[4] = p.e5x5()/s->rawEnergy();
144  fVals[5] = p.hadronicOverEm();
145  fVals[6] = s->etaWidth();
146  fVals[7] = s->phiWidth();
147 
148 
149  //seed basic cluster variables
150  double bemax = clustertools.eMax(*b);
151  double be2nd = clustertools.e2nd(*b);
152  double betop = clustertools.eTop(*b);
153  double bebottom = clustertools.eBottom(*b);
154  double beleft = clustertools.eLeft(*b);
155  double beright = clustertools.eRight(*b);
156 
157 
158  fVals[8] = b->eta()-s->eta();
159  fVals[9] = reco::deltaPhi(b->phi(),s->phi());
160  fVals[10] = b->energy()/s->rawEnergy();
161  fVals[11] = clustertools.e3x3(*b)/b->energy();
162  fVals[12] = clustertools.e5x5(*b)/b->energy();
163  fVals[13] = sqrt(clustertools.localCovariances(*b)[0]); //sigietaieta
164  fVals[14] = sqrt(clustertools.localCovariances(*b)[2]); //sigiphiiphi
165  fVals[15] = clustertools.localCovariances(*b)[1]; //sigietaiphi
166  fVals[16] = bemax/b->energy(); //crystal energy ratio gap variables
167  fVals[17] = log(be2nd/bemax);
168  fVals[18] = log(betop/bemax);
169  fVals[19] = log(bebottom/bemax);
170  fVals[20] = log(beleft/bemax);
171  fVals[21] = log(beright/bemax);
172  fVals[22] = (betop-bebottom)/(betop+bebottom);
173  fVals[23] = (beleft-beright)/(beleft+beright);
174 
175 
176  double bc2emax = hasbc2 ? clustertools.eMax(*b2) : 0.;
177  double bc2e2nd = hasbc2 ? clustertools.e2nd(*b2) : 0.;
178  double bc2etop = hasbc2 ? clustertools.eTop(*b2) : 0.;
179  double bc2ebottom = hasbc2 ? clustertools.eBottom(*b2) : 0.;
180  double bc2eleft = hasbc2 ? clustertools.eLeft(*b2) : 0.;
181  double bc2eright = hasbc2 ? clustertools.eRight(*b2) : 0.;
182 
183  fVals[24] = hasbc2 ? (b2->eta()-s->eta()) : 0.;
184  fVals[25] = hasbc2 ? reco::deltaPhi(b2->phi(),s->phi()) : 0.;
185  fVals[26] = hasbc2 ? b2->energy()/s->rawEnergy() : 0.;
186  fVals[27] = hasbc2 ? clustertools.e3x3(*b2)/b2->energy() : 0.;
187  fVals[28] = hasbc2 ? clustertools.e5x5(*b2)/b2->energy() : 0.;
188  fVals[29] = hasbc2 ? sqrt(clustertools.localCovariances(*b2)[0]) : 0.;
189  fVals[30] = hasbc2 ? sqrt(clustertools.localCovariances(*b2)[2]) : 0.;
190  fVals[31] = hasbc2 ? clustertools.localCovariances(*b)[1] : 0.;
191  fVals[32] = hasbc2 ? bc2emax/b2->energy() : 0.;
192  fVals[33] = hasbc2 ? log(bc2e2nd/bc2emax) : 0.;
193  fVals[34] = hasbc2 ? log(bc2etop/bc2emax) : 0.;
194  fVals[35] = hasbc2 ? log(bc2ebottom/bc2emax) : 0.;
195  fVals[36] = hasbc2 ? log(bc2eleft/bc2emax) : 0.;
196  fVals[37] = hasbc2 ? log(bc2eright/bc2emax) : 0.;
197  fVals[38] = hasbc2 ? (bc2etop-bc2ebottom)/(bc2etop+bc2ebottom) : 0.;
198  fVals[39] = hasbc2 ? (bc2eleft-bc2eright)/(bc2eleft+bc2eright) : 0.;
199 
200  fVals[40] = hasbclast ? (bclast->eta()-s->eta()) : 0.;
201  fVals[41] = hasbclast ? reco::deltaPhi(bclast->phi(),s->phi()) : 0.;
202  fVals[42] = hasbclast ? bclast->energy()/s->rawEnergy() : 0.;
203  fVals[43] = hasbclast ? clustertools.e3x3(*bclast)/bclast->energy() : 0.;
204  fVals[44] = hasbclast ? clustertools.e5x5(*bclast)/bclast->energy() : 0.;
205  fVals[45] = hasbclast ? sqrt(clustertools.localCovariances(*bclast)[0]) : 0.;
206  fVals[46] = hasbclast ? sqrt(clustertools.localCovariances(*bclast)[2]) : 0.;
207  fVals[47] = hasbclast ? clustertools.localCovariances(*bclast)[1] : 0.;
208 
209  fVals[48] = hasbclast2 ? (bclast2->eta()-s->eta()) : 0.;
210  fVals[49] = hasbclast2 ? reco::deltaPhi(bclast2->phi(),s->phi()) : 0.;
211  fVals[50] = hasbclast2 ? bclast2->energy()/s->rawEnergy() : 0.;
212  fVals[51] = hasbclast2 ? clustertools.e3x3(*bclast2)/bclast2->energy() : 0.;
213  fVals[52] = hasbclast2 ? clustertools.e5x5(*bclast2)/bclast2->energy() : 0.;
214  fVals[53] = hasbclast2 ? sqrt(clustertools.localCovariances(*bclast2)[0]) : 0.;
215  fVals[54] = hasbclast2 ? sqrt(clustertools.localCovariances(*bclast2)[2]) : 0.;
216  fVals[55] = hasbclast2 ? clustertools.localCovariances(*bclast2)[1] : 0.;
217 
218 
219  //local coordinates and crystal indices
220 
221 
222  //seed cluster
223  float betacry, bphicry, bthetatilt, bphitilt;
224  int bieta, biphi;
225  _ecalLocal.localCoordsEB(*b,es,betacry,bphicry,bieta,biphi,bthetatilt,bphitilt);
226 
227  fVals[56] = bieta; //crystal ieta
228  fVals[57] = biphi; //crystal iphi
229  fVals[58] = bieta%5; //submodule boundary eta symmetry
230  fVals[59] = biphi%2; //submodule boundary phi symmetry
231  fVals[60] = (TMath::Abs(bieta)<=25)*(bieta%25) + (TMath::Abs(bieta)>25)*((bieta-25*TMath::Abs(bieta)/bieta)%20); //module boundary eta approximate symmetry
232  fVals[61] = biphi%20; //module boundary phi symmetry
233  fVals[62] = betacry; //local coordinates with respect to closest crystal center at nominal shower depth
234  fVals[63] = bphicry;
235 
236 
237  //2nd cluster (meaningful gap corrections for converted photons)
238  float bc2etacry, bc2phicry, bc2thetatilt, bc2phitilt;
239  int bc2ieta, bc2iphi;
240  if (hasbc2) _ecalLocal.localCoordsEB(*b2,es,bc2etacry,bc2phicry,bc2ieta,bc2iphi,bc2thetatilt,bc2phitilt);
241 
242  fVals[64] = hasbc2 ? bc2ieta : 0.;
243  fVals[65] = hasbc2 ? bc2iphi : 0.;
244  fVals[66] = hasbc2 ? bc2ieta%5 : 0.;
245  fVals[67] = hasbc2 ? bc2iphi%2 : 0.;
246  fVals[68] = hasbc2 ? (TMath::Abs(bc2ieta)<=25)*(bc2ieta%25) + (TMath::Abs(bc2ieta)>25)*((bc2ieta-25*TMath::Abs(bc2ieta)/bc2ieta)%20) : 0.;
247  fVals[69] = hasbc2 ? bc2iphi%20 : 0.;
248  fVals[70] = hasbc2 ? bc2etacry : 0.;
249  fVals[71] = hasbc2 ? bc2phicry : 0.;
250 
251  fVals[72] = vtxcol.size();
252 
253  }
254  else {
255  fVals[0] = s->rawEnergy();
256  fVals[1] = p.r9();
257  fVals[2] = s->eta();
258  fVals[3] = s->phi();
259  fVals[4] = p.e5x5()/s->rawEnergy();
260  fVals[5] = s->etaWidth();
261  fVals[6] = s->phiWidth();
262  fVals[7] = vtxcol.size();
263  }
264 
265 
266  const Double_t varscale = 1.253;
267  Double_t den;
268  const GBRForest *reader;
269  const GBRForest *readervar;
270  if (isbarrel) {
271  den = s->rawEnergy();
272  reader = fReadereb;
273  readervar = fReaderebvariance;
274  }
275  else {
276  den = s->rawEnergy() + s->preshowerEnergy();
277  reader = fReaderee;
278  readervar = fReadereevariance;
279  }
280 
281  Double_t ecor = reader->GetResponse(fVals)*den;
282  Double_t ecorerr = readervar->GetResponse(fVals)*den*varscale;
283 
284  //printf("ecor = %5f, ecorerr = %5f\n",ecor,ecorerr);
285 
286  return std::pair<double,double>(ecor,ecorerr);
287 }
288 
289 
290 //--------------------------------------------------------------------------------------------------
291 std::pair<double,double> EGEnergyCorrector::CorrectedEnergyWithError(const GsfElectron &e, const reco::VertexCollection& vtxcol, EcalClusterLazyTools &clustertools, const edm::EventSetup &es) {
292 
293  //apply v2 regression to electrons
294  //mostly duplicated from photon function above //TODO, make common underlying function
295 
296  //protection, this doesn't work properly on non-egamma-seeded electrons
297  if (!e.ecalDrivenSeed()) return std::pair<double,double>(0.,0.);
298 
299 
300  const SuperClusterRef s = e.superCluster();
301  const CaloClusterPtr b = s->seed();
302 
303  CaloClusterPtr b2;
304  Double_t ebcmax = -99.;
305  for (reco::CaloCluster_iterator bit = s->clustersBegin(); bit!=s->clustersEnd(); ++bit) {
306  const CaloClusterPtr bc = *bit;
307  if (bc->energy() > ebcmax && bc !=b) {
308  b2 = bc;
309  ebcmax = bc->energy();
310  }
311  }
312 
313  CaloClusterPtr bclast;
314  Double_t ebcmin = 1e6;
315  for (reco::CaloCluster_iterator bit = s->clustersBegin(); bit!=s->clustersEnd(); ++bit) {
316  const CaloClusterPtr bc = *bit;
317  if (bc->energy() < ebcmin && bc !=b) {
318  bclast = bc;
319  ebcmin = bc->energy();
320  }
321  }
322 
323  CaloClusterPtr bclast2;
324  ebcmin = 1e6;
325  for (reco::CaloCluster_iterator bit = s->clustersBegin(); bit!=s->clustersEnd(); ++bit) {
326  const CaloClusterPtr bc = *bit;
327  if (bc->energy() < ebcmin && bc !=b && bc!=bclast) {
328  bclast2 = bc;
329  ebcmin = bc->energy();
330  }
331  }
332 
333  Bool_t isbarrel = b->hitsAndFractions().at(0).first.subdetId()==EcalBarrel;
334  Bool_t hasbc2 = b2.isNonnull() && b2->energy()>0.;
335  Bool_t hasbclast = bclast.isNonnull() && bclast->energy()>0.;
336  Bool_t hasbclast2 = bclast2.isNonnull() && bclast2->energy()>0.;
337 
338 
339  if (isbarrel) {
340  fVals[0] = s->rawEnergy();
341  fVals[1] = clustertools.e3x3(*b)/s->rawEnergy(); //r9
342  fVals[2] = s->eta();
343  fVals[3] = s->phi();
344  fVals[4] = clustertools.e5x5(*b)/s->rawEnergy();
345  fVals[5] = e.hcalOverEcal();
346  fVals[6] = s->etaWidth();
347  fVals[7] = s->phiWidth();
348 
349 
350  double bemax = clustertools.eMax(*b);
351  double be2nd = clustertools.e2nd(*b);
352  double betop = clustertools.eTop(*b);
353  double bebottom = clustertools.eBottom(*b);
354  double beleft = clustertools.eLeft(*b);
355  double beright = clustertools.eRight(*b);
356 
357 
358  fVals[8] = b->eta()-s->eta();
359  fVals[9] = reco::deltaPhi(b->phi(),s->phi());
360  fVals[10] = b->energy()/s->rawEnergy();
361  fVals[11] = clustertools.e3x3(*b)/b->energy();
362  fVals[12] = clustertools.e5x5(*b)/b->energy();
363  fVals[13] = sqrt(clustertools.localCovariances(*b)[0]);
364  fVals[14] = sqrt(clustertools.localCovariances(*b)[2]);
365  fVals[15] = clustertools.localCovariances(*b)[1];
366  fVals[16] = bemax/b->energy();
367  fVals[17] = log(be2nd/bemax);
368  fVals[18] = log(betop/bemax);
369  fVals[19] = log(bebottom/bemax);
370  fVals[20] = log(beleft/bemax);
371  fVals[21] = log(beright/bemax);
372  fVals[22] = (betop-bebottom)/(betop+bebottom);
373  fVals[23] = (beleft-beright)/(beleft+beright);
374 
375 
376  double bc2emax = hasbc2 ? clustertools.eMax(*b2) : 0.;
377  double bc2e2nd = hasbc2 ? clustertools.e2nd(*b2) : 0.;
378  double bc2etop = hasbc2 ? clustertools.eTop(*b2) : 0.;
379  double bc2ebottom = hasbc2 ? clustertools.eBottom(*b2) : 0.;
380  double bc2eleft = hasbc2 ? clustertools.eLeft(*b2) : 0.;
381  double bc2eright = hasbc2 ? clustertools.eRight(*b2) : 0.;
382 
383  fVals[24] = hasbc2 ? (b2->eta()-s->eta()) : 0.;
384  fVals[25] = hasbc2 ? reco::deltaPhi(b2->phi(),s->phi()) : 0.;
385  fVals[26] = hasbc2 ? b2->energy()/s->rawEnergy() : 0.;
386  fVals[27] = hasbc2 ? clustertools.e3x3(*b2)/b2->energy() : 0.;
387  fVals[28] = hasbc2 ? clustertools.e5x5(*b2)/b2->energy() : 0.;
388  fVals[29] = hasbc2 ? sqrt(clustertools.localCovariances(*b2)[0]) : 0.;
389  fVals[30] = hasbc2 ? sqrt(clustertools.localCovariances(*b2)[2]) : 0.;
390  fVals[31] = hasbc2 ? clustertools.localCovariances(*b)[1] : 0.;
391  fVals[32] = hasbc2 ? bc2emax/b2->energy() : 0.;
392  fVals[33] = hasbc2 ? log(bc2e2nd/bc2emax) : 0.;
393  fVals[34] = hasbc2 ? log(bc2etop/bc2emax) : 0.;
394  fVals[35] = hasbc2 ? log(bc2ebottom/bc2emax) : 0.;
395  fVals[36] = hasbc2 ? log(bc2eleft/bc2emax) : 0.;
396  fVals[37] = hasbc2 ? log(bc2eright/bc2emax) : 0.;
397  fVals[38] = hasbc2 ? (bc2etop-bc2ebottom)/(bc2etop+bc2ebottom) : 0.;
398  fVals[39] = hasbc2 ? (bc2eleft-bc2eright)/(bc2eleft+bc2eright) : 0.;
399 
400  fVals[40] = hasbclast ? (bclast->eta()-s->eta()) : 0.;
401  fVals[41] = hasbclast ? reco::deltaPhi(bclast->phi(),s->phi()) : 0.;
402  fVals[42] = hasbclast ? bclast->energy()/s->rawEnergy() : 0.;
403  fVals[43] = hasbclast ? clustertools.e3x3(*bclast)/bclast->energy() : 0.;
404  fVals[44] = hasbclast ? clustertools.e5x5(*bclast)/bclast->energy() : 0.;
405  fVals[45] = hasbclast ? sqrt(clustertools.localCovariances(*bclast)[0]) : 0.;
406  fVals[46] = hasbclast ? sqrt(clustertools.localCovariances(*bclast)[2]) : 0.;
407  fVals[47] = hasbclast ? clustertools.localCovariances(*bclast)[1] : 0.;
408 
409  fVals[48] = hasbclast2 ? (bclast2->eta()-s->eta()) : 0.;
410  fVals[49] = hasbclast2 ? reco::deltaPhi(bclast2->phi(),s->phi()) : 0.;
411  fVals[50] = hasbclast2 ? bclast2->energy()/s->rawEnergy() : 0.;
412  fVals[51] = hasbclast2 ? clustertools.e3x3(*bclast2)/bclast2->energy() : 0.;
413  fVals[52] = hasbclast2 ? clustertools.e5x5(*bclast2)/bclast2->energy() : 0.;
414  fVals[53] = hasbclast2 ? sqrt(clustertools.localCovariances(*bclast2)[0]) : 0.;
415  fVals[54] = hasbclast2 ? sqrt(clustertools.localCovariances(*bclast2)[2]) : 0.;
416  fVals[55] = hasbclast2 ? clustertools.localCovariances(*bclast2)[1] : 0.;
417 
418 
419  float betacry, bphicry, bthetatilt, bphitilt;
420  int bieta, biphi;
421  _ecalLocal.localCoordsEB(*b,es,betacry,bphicry,bieta,biphi,bthetatilt,bphitilt);
422 
423  fVals[56] = bieta;
424  fVals[57] = biphi;
425  fVals[58] = bieta%5;
426  fVals[59] = biphi%2;
427  fVals[60] = (TMath::Abs(bieta)<=25)*(bieta%25) + (TMath::Abs(bieta)>25)*((bieta-25*TMath::Abs(bieta)/bieta)%20);
428  fVals[61] = biphi%20;
429  fVals[62] = betacry;
430  fVals[63] = bphicry;
431 
432  float bc2etacry, bc2phicry, bc2thetatilt, bc2phitilt;
433  int bc2ieta, bc2iphi;
434  if (hasbc2) _ecalLocal.localCoordsEB(*b2,es,bc2etacry,bc2phicry,bc2ieta,bc2iphi,bc2thetatilt,bc2phitilt);
435 
436  fVals[64] = hasbc2 ? bc2ieta : 0.;
437  fVals[65] = hasbc2 ? bc2iphi : 0.;
438  fVals[66] = hasbc2 ? bc2ieta%5 : 0.;
439  fVals[67] = hasbc2 ? bc2iphi%2 : 0.;
440  fVals[68] = hasbc2 ? (TMath::Abs(bc2ieta)<=25)*(bc2ieta%25) + (TMath::Abs(bc2ieta)>25)*((bc2ieta-25*TMath::Abs(bc2ieta)/bc2ieta)%20) : 0.;
441  fVals[69] = hasbc2 ? bc2iphi%20 : 0.;
442  fVals[70] = hasbc2 ? bc2etacry : 0.;
443  fVals[71] = hasbc2 ? bc2phicry : 0.;
444 
445  fVals[72] = vtxcol.size();
446 
447  }
448  else {
449  fVals[0] = s->rawEnergy();
450  fVals[1] = clustertools.e3x3(*b)/s->rawEnergy(); //r9
451  fVals[2] = s->eta();
452  fVals[3] = s->phi();
453  fVals[4] = clustertools.e5x5(*b)/s->rawEnergy();
454  fVals[5] = s->etaWidth();
455  fVals[6] = s->phiWidth();
456  fVals[7] = vtxcol.size();
457  }
458 
459 
460  const Double_t varscale = 1.253;
461  Double_t den;
462  const GBRForest *reader;
463  const GBRForest *readervar;
464  if (isbarrel) {
465  den = s->rawEnergy();
466  reader = fReadereb;
467  readervar = fReaderebvariance;
468  }
469  else {
470  den = s->rawEnergy() + s->preshowerEnergy();
471  reader = fReaderee;
472  readervar = fReadereevariance;
473  }
474 
475  Double_t ecor = reader->GetResponse(fVals)*den;
476  Double_t ecorerr = readervar->GetResponse(fVals)*den*varscale;
477 
478  //printf("ecor = %5f, ecorerr = %5f\n",ecor,ecorerr);
479 
480  return std::pair<double,double>(ecor,ecorerr);
481 
482 }
483 
484 //--------------------------------------------------------------------------------------------------
485 std::pair<double,double> EGEnergyCorrector::CorrectedEnergyWithErrorV3(const Photon &p, const reco::VertexCollection& vtxcol, double rho, EcalClusterLazyTools &clustertools, const edm::EventSetup &es, bool applyRescale) {
486 
487  const SuperClusterRef s = p.superCluster();
488  const CaloClusterPtr b = s->seed(); //seed basic cluster
489 
490  Bool_t isbarrel = b->hitsAndFractions().at(0).first.subdetId()==EcalBarrel;
491 
492  //basic supercluster variables
493  fVals[0] = s->rawEnergy();
494  fVals[1] = s->eta();
495  fVals[2] = s->phi();
496  fVals[3] = p.r9();
497  fVals[4] = p.e5x5()/s->rawEnergy();
498  fVals[5] = s->etaWidth();
499  fVals[6] = s->phiWidth();
500  fVals[7] = s->clustersSize();
501  fVals[8] = p.hadTowOverEm();
502  fVals[9] = rho;
503  fVals[10] = vtxcol.size();
504 
505  //seed basic cluster variables
506  double bemax = clustertools.eMax(*b);
507  double be2nd = clustertools.e2nd(*b);
508  double betop = clustertools.eTop(*b);
509  double bebottom = clustertools.eBottom(*b);
510  double beleft = clustertools.eLeft(*b);
511  double beright = clustertools.eRight(*b);
512 
513  double be2x5max = clustertools.e2x5Max(*b);
514  double be2x5top = clustertools.e2x5Top(*b);
515  double be2x5bottom = clustertools.e2x5Bottom(*b);
516  double be2x5left = clustertools.e2x5Left(*b);
517  double be2x5right = clustertools.e2x5Right(*b);
518 
519  fVals[11] = b->eta()-s->eta();
520  fVals[12] = reco::deltaPhi(b->phi(),s->phi());
521  fVals[13] = b->energy()/s->rawEnergy();
522  fVals[14] = clustertools.e3x3(*b)/b->energy();
523  fVals[15] = clustertools.e5x5(*b)/b->energy();
524  fVals[16] = sqrt(clustertools.localCovariances(*b)[0]); //sigietaieta
525  fVals[17] = sqrt(clustertools.localCovariances(*b)[2]); //sigiphiiphi
526  fVals[18] = clustertools.localCovariances(*b)[1]; //sigietaiphi
527  fVals[19] = bemax/b->energy(); //crystal energy ratio gap variables
528  fVals[20] = be2nd/b->energy();
529  fVals[21] = betop/b->energy();
530  fVals[22] = bebottom/b->energy();
531  fVals[23] = beleft/b->energy();
532  fVals[24] = beright/b->energy();
533  fVals[25] = be2x5max/b->energy(); //crystal energy ratio gap variables
534  fVals[26] = be2x5top/b->energy();
535  fVals[27] = be2x5bottom/b->energy();
536  fVals[28] = be2x5left/b->energy();
537  fVals[29] = be2x5right/b->energy();
538 
539  if (isbarrel) {
540  //local coordinates and crystal indices (barrel only)
541 
542  //seed cluster
543  float betacry, bphicry, bthetatilt, bphitilt;
544  int bieta, biphi;
545  _ecalLocal.localCoordsEB(*b,es,betacry,bphicry,bieta,biphi,bthetatilt,bphitilt);
546 
547  fVals[30] = bieta; //crystal ieta
548  fVals[31] = biphi; //crystal iphi
549  fVals[32] = bieta%5; //submodule boundary eta symmetry
550  fVals[33] = biphi%2; //submodule boundary phi symmetry
551  fVals[34] = (TMath::Abs(bieta)<=25)*(bieta%25) + (TMath::Abs(bieta)>25)*((bieta-25*TMath::Abs(bieta)/bieta)%20); //module boundary eta approximate symmetry
552  fVals[35] = biphi%20; //module boundary phi symmetry
553  fVals[36] = betacry; //local coordinates with respect to closest crystal center at nominal shower depth
554  fVals[37] = bphicry;
555 
556  }
557  else {
558  //preshower energy ratio (endcap only)
559  fVals[30] = s->preshowerEnergy()/s->rawEnergy();
560  }
561 
562 // if (isbarrel) {
563 // for (int i=0; i<38; ++i) printf("%i: %5f\n",i,fVals[i]);
564 // }
565 // else for (int i=0; i<31; ++i) printf("%i: %5f\n",i,fVals[i]);
566 
567  Double_t den;
568  const GBRForest *reader;
569  const GBRForest *readervar;
570  if (isbarrel) {
571  den = s->rawEnergy();
572  reader = fReadereb;
573  readervar = fReaderebvariance;
574  }
575  else {
576  den = s->rawEnergy() + s->preshowerEnergy();
577  reader = fReaderee;
578  readervar = fReadereevariance;
579  }
580 
581  Double_t ecor = reader->GetResponse(fVals)*den;
582 
583 
584  //apply shower shape rescaling - for Monte Carlo only, and only for calculation of energy uncertainty
585  if (applyRescale) {
586  if (isbarrel) {
587  fVals[3] = 1.0045*p.r9() +0.001; //r9
588  fVals[5] = 1.04302*s->etaWidth() - 0.000618; //etawidth
589  fVals[6] = 1.00002*s->phiWidth() - 0.000371; //phiwidth
590  fVals[14] = fVals[3]*s->rawEnergy()/b->energy(); //compute consistent e3x3/eseed after r9 rescaling
591  if (fVals[15]<=1.0) // rescale e5x5/eseed only if value is <=1.0, don't allow scaled values to exceed 1.0
592  fVals[15] = TMath::Min(1.0,1.0022*p.e5x5()/b->energy());
593 
594  fVals[4] = fVals[15]*b->energy()/s->rawEnergy(); // compute consistent e5x5()/rawEnergy() after e5x5/eseed resacling
595 
596  fVals[16] = 0.891832*sqrt(clustertools.localCovariances(*b)[0]) + 0.0009133; //sigietaieta
597  fVals[17] = 0.993*sqrt(clustertools.localCovariances(*b)[2]); //sigiphiiphi
598 
599  fVals[19] = 1.012*bemax/b->energy(); //crystal energy ratio gap variables
600  fVals[20] = 1.0*be2nd/b->energy();
601  fVals[21] = 0.94*betop/b->energy();
602  fVals[22] = 0.94*bebottom/b->energy();
603  fVals[23] = 0.94*beleft/b->energy();
604  fVals[24] = 0.94*beright/b->energy();
605  fVals[25] = 1.006*be2x5max/b->energy(); //crystal energy ratio gap variables
606  fVals[26] = 1.09*be2x5top/b->energy();
607  fVals[27] = 1.09*be2x5bottom/b->energy();
608  fVals[28] = 1.09*be2x5left/b->energy();
609  fVals[29] = 1.09*be2x5right/b->energy();
610 
611  }
612  else {
613  fVals[3] = 1.0086*p.r9() -0.0007; //r9
614  fVals[4] = TMath::Min(1.0,1.0022*p.e5x5()/s->rawEnergy()); //e5x5/rawenergy
615  fVals[5] = 0.903254*s->etaWidth() +0.001346; //etawidth
616  fVals[6] = 0.99992*s->phiWidth() + 4.8e-07; //phiwidth
617  fVals[13] = TMath::Min(1.0,1.0022*b->energy()/s->rawEnergy()); //eseed/rawenergy (practically equivalent to e5x5)
618 
619 
620  fVals[14] = fVals[3]*s->rawEnergy()/b->energy(); //compute consistent e3x3/eseed after r9 rescaling
621 
622  fVals[16] = 0.9947*sqrt(clustertools.localCovariances(*b)[0]) + 0.00003; //sigietaieta
623 
624  fVals[19] = 1.005*bemax/b->energy(); //crystal energy ratio gap variables
625  fVals[20] = 1.02*be2nd/b->energy();
626  fVals[21] = 0.96*betop/b->energy();
627  fVals[22] = 0.96*bebottom/b->energy();
628  fVals[23] = 0.96*beleft/b->energy();
629  fVals[24] = 0.96*beright/b->energy();
630  fVals[25] = 1.0075*be2x5max/b->energy(); //crystal energy ratio gap variables
631  fVals[26] = 1.13*be2x5top/b->energy();
632  fVals[27] = 1.13*be2x5bottom/b->energy();
633  fVals[28] = 1.13*be2x5left/b->energy();
634  fVals[29] = 1.13*be2x5right/b->energy();
635  }
636 
637  }
638 
639  Double_t ecorerr = readervar->GetResponse(fVals)*den;
640 
641  //printf("ecor = %5f, ecorerr = %5f\n",ecor,ecorerr);
642 
643  return std::pair<double,double>(ecor,ecorerr);
644 }
645 
646 //--------------------------------------------------------------------------------------------------
647 std::pair<double,double> EGEnergyCorrector::CorrectedEnergyWithErrorV3(const GsfElectron &e, const reco::VertexCollection& vtxcol, double rho, EcalClusterLazyTools &clustertools, const edm::EventSetup &es) {
648 
649  const SuperClusterRef s = e.superCluster();
650  const CaloClusterPtr b = s->seed(); //seed basic cluster
651 
652  Bool_t isbarrel = b->hitsAndFractions().at(0).first.subdetId()==EcalBarrel;
653 
654  //basic supercluster variables
655  fVals[0] = s->rawEnergy();
656  fVals[1] = s->eta();
657  fVals[2] = s->phi();
658  fVals[3] = clustertools.e3x3(*b)/s->rawEnergy(); //r9
659  fVals[4] = clustertools.e5x5(*b)/s->rawEnergy();
660  fVals[5] = s->etaWidth();
661  fVals[6] = s->phiWidth();
662  fVals[7] = s->clustersSize();
663  fVals[8] = e.hcalOverEcalBc();
664  fVals[9] = rho;
665  fVals[10] = vtxcol.size();
666 
667  //seed basic cluster variables
668  double bemax = clustertools.eMax(*b);
669  double be2nd = clustertools.e2nd(*b);
670  double betop = clustertools.eTop(*b);
671  double bebottom = clustertools.eBottom(*b);
672  double beleft = clustertools.eLeft(*b);
673  double beright = clustertools.eRight(*b);
674 
675  double be2x5max = clustertools.e2x5Max(*b);
676  double be2x5top = clustertools.e2x5Top(*b);
677  double be2x5bottom = clustertools.e2x5Bottom(*b);
678  double be2x5left = clustertools.e2x5Left(*b);
679  double be2x5right = clustertools.e2x5Right(*b);
680 
681  fVals[11] = b->eta()-s->eta();
682  fVals[12] = reco::deltaPhi(b->phi(),s->phi());
683  fVals[13] = b->energy()/s->rawEnergy();
684  fVals[14] = clustertools.e3x3(*b)/b->energy();
685  fVals[15] = clustertools.e5x5(*b)/b->energy();
686  fVals[16] = sqrt(clustertools.localCovariances(*b)[0]); //sigietaieta
687  fVals[17] = sqrt(clustertools.localCovariances(*b)[2]); //sigiphiiphi
688  fVals[18] = clustertools.localCovariances(*b)[1]; //sigietaiphi
689  fVals[19] = bemax/b->energy(); //crystal energy ratio gap variables
690  fVals[20] = be2nd/b->energy();
691  fVals[21] = betop/b->energy();
692  fVals[22] = bebottom/b->energy();
693  fVals[23] = beleft/b->energy();
694  fVals[24] = beright/b->energy();
695  fVals[25] = be2x5max/b->energy(); //crystal energy ratio gap variables
696  fVals[26] = be2x5top/b->energy();
697  fVals[27] = be2x5bottom/b->energy();
698  fVals[28] = be2x5left/b->energy();
699  fVals[29] = be2x5right/b->energy();
700 
701  if (isbarrel) {
702  //local coordinates and crystal indices (barrel only)
703 
704  //seed cluster
705  float betacry, bphicry, bthetatilt, bphitilt;
706  int bieta, biphi;
707  _ecalLocal.localCoordsEB(*b,es,betacry,bphicry,bieta,biphi,bthetatilt,bphitilt);
708 
709  fVals[30] = bieta; //crystal ieta
710  fVals[31] = biphi; //crystal iphi
711  fVals[32] = bieta%5; //submodule boundary eta symmetry
712  fVals[33] = biphi%2; //submodule boundary phi symmetry
713  fVals[34] = (TMath::Abs(bieta)<=25)*(bieta%25) + (TMath::Abs(bieta)>25)*((bieta-25*TMath::Abs(bieta)/bieta)%20); //module boundary eta approximate symmetry
714  fVals[35] = biphi%20; //module boundary phi symmetry
715  fVals[36] = betacry; //local coordinates with respect to closest crystal center at nominal shower depth
716  fVals[37] = bphicry;
717 
718  }
719  else {
720  //preshower energy ratio (endcap only)
721  fVals[30] = s->preshowerEnergy()/s->rawEnergy();
722  }
723 
724  Double_t den;
725  const GBRForest *reader;
726  const GBRForest *readervar;
727  if (isbarrel) {
728  den = s->rawEnergy();
729  reader = fReadereb;
730  readervar = fReaderebvariance;
731  }
732  else {
733  den = s->rawEnergy() + s->preshowerEnergy();
734  reader = fReaderee;
735  readervar = fReadereevariance;
736  }
737 
738  Double_t ecor = reader->GetResponse(fVals)*den;
739  Double_t ecorerr = readervar->GetResponse(fVals)*den;
740 
741  //printf("ecor = %5f, ecorerr = %5f\n",ecor,ecorerr);
742 
743  return std::pair<double,double>(ecor,ecorerr);
744 }
745 
EcalClusterLazyTools::localCovariances
std::vector< float > localCovariances(const reco::BasicCluster &cluster, float w0=4.7)
Definition: EcalClusterLazyTools.cc:471
GBRForest::GetResponse
double GetResponse(const float *vector) const
Definition: GBRForest.h:51
EcalClusterLazyTools::eBottom
float eBottom(const reco::BasicCluster &cluster)
Definition: EcalClusterLazyTools.cc:385
EGEnergyCorrector::Initialize
void Initialize(const edm::EventSetup &iSetup, std::string regweights, bool weightsFromDB=false)
Definition: EGEnergyCorrector.cc:47
EcalClusterLazyTools::e2nd
float e2nd(const reco::BasicCluster &cluster)
Definition: EcalClusterLazyTools.cc:407
reco::Photon::e5x5
float e5x5() const
Definition: Photon.h:185
siStripFEDMonitor_P5_cff.Min
tuple Min
Definition: siStripFEDMonitor_P5_cff.py:55
EGEnergyCorrector::fReaderee
const GBRForest * fReaderee
Definition: EGEnergyCorrector.h:40
EGEnergyCorrector::fReadereb
const GBRForest * fReadereb
Definition: EGEnergyCorrector.h:38
EGEnergyCorrector::fReaderebvariance
const GBRForest * fReaderebvariance
Definition: EGEnergyCorrector.h:39
EcalClusterLazyTools::eRight
float eRight(const reco::BasicCluster &cluster)
Definition: EcalClusterLazyTools.cc:363
rho
Definition: DDAxes.h:10
create_public_lumi_plots.log
log
Definition: create_public_lumi_plots.py:1108
EcalClusterLazyTools::e3x3
float e3x3(const reco::BasicCluster &cluster)
Definition: EcalClusterLazyTools.cc:264
EcalClusterLazyTools::eMax
float eMax(const reco::BasicCluster &cluster)
Definition: EcalClusterLazyTools.cc:396
reco::Photon::superCluster
reco::SuperClusterRef superCluster() const
Ref to SuperCluster.
Definition: Photon.cc:59
reco::VertexCollection
std::vector< Vertex > VertexCollection
collection of Vertex objects
Definition: VertexFwd.h:9
dt_dqm_sourceclient_common_cff.reco
tuple reco
Definition: dt_dqm_sourceclient_common_cff.py:105
EGEnergyCorrector::CorrectedEnergyWithError
std::pair< double, double > CorrectedEnergyWithError(const reco::Photon &p, const reco::VertexCollection &vtxcol, EcalClusterLazyTools &clustertools, const edm::EventSetup &es)
Definition: EGEnergyCorrector.cc:92
edm::Ptr::isNonnull
bool isNonnull() const
Checks for non-null.
Definition: Ptr.h:152
EcalClusterLazyTools::e2x5Max
float e2x5Max(const reco::BasicCluster &cluster)
Definition: EcalClusterLazyTools.cc:341
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:46
reco::GsfElectron::superCluster
virtual SuperClusterRef superCluster() const
reference to a SuperCluster
Definition: GsfElectron.h:168
EcalClusterLazyTools::e2x5Right
float e2x5Right(const reco::BasicCluster &cluster)
Definition: EcalClusterLazyTools.cc:297
reco::GsfElectron::hcalOverEcal
float hcalOverEcal() const
Definition: GsfElectron.h:390
reco::Photon::hadTowOverEm
float hadTowOverEm() const
the ration of hadronic energy in towers behind the BCs in the SC and the SC energy ...
Definition: Photon.h:174
reco::Photon::hadronicOverEm
float hadronicOverEm() const
the total hadronic over electromagnetic fraction
Definition: Photon.h:167
reco::GsfElectron::hcalOverEcalBc
float hcalOverEcalBc() const
Definition: GsfElectron.h:394
EcalClusterLazyTools::eTop
float eTop(const reco::BasicCluster &cluster)
Definition: EcalClusterLazyTools.cc:374
EcalClusterLazyTools::e2x5Bottom
float e2x5Bottom(const reco::BasicCluster &cluster)
Definition: EcalClusterLazyTools.cc:330
reco::deltaPhi
double deltaPhi(double phi1, double phi2)
Definition: deltaPhi.h:12
EGEnergyCorrector::CorrectedEnergyWithErrorV3
std::pair< double, double > CorrectedEnergyWithErrorV3(const reco::Photon &p, const reco::VertexCollection &vtxcol, double rho, EcalClusterLazyTools &clustertools, const edm::EventSetup &es, bool applyRescale=false)
Definition: EGEnergyCorrector.cc:485
EcalClusterLazyTools::e2x5Left
float e2x5Left(const reco::BasicCluster &cluster)
Definition: EcalClusterLazyTools.cc:308
EcalClusterLazyTools::e5x5
float e5x5(const reco::BasicCluster &cluster)
Definition: EcalClusterLazyTools.cc:286
EcalClusterLazyTools::eLeft
float eLeft(const reco::BasicCluster &cluster)
Definition: EcalClusterLazyTools.cc:352
edm::EventSetup::get
const T & get() const
Definition: EventSetup.h:55
edm::ESHandle::product
T const * product() const
Definition: ESHandle.h:62
b
double b
Definition: hdecay.h:120
EcalClusterLocal::localCoordsEB
void localCoordsEB(const reco::CaloCluster &bclus, const edm::EventSetup &es, float &etacry, float &phicry, int &ieta, int &iphi, float &thetatilt, float &phitilt) const
Definition: EcalClusterLocal.cc:22
AlCaHLTBitMon_ParallelJobs.p
tuple p
Definition: AlCaHLTBitMon_ParallelJobs.py:152
EGEnergyCorrector::fReadereevariance
const GBRForest * fReadereevariance
Definition: EGEnergyCorrector.h:41
EGEnergyCorrector::_ecalLocal
EcalClusterLocal _ecalLocal
Definition: EGEnergyCorrector.h:47
reco::Photon::r9
float r9() const
Definition: Photon.h:191
EcalClusterLazyTools::e2x5Top
float e2x5Top(const reco::BasicCluster &cluster)
Definition: EcalClusterLazyTools.cc:319
reco::GsfElectron::ecalDrivenSeed
bool ecalDrivenSeed() const
Definition: GsfElectron.h:172
Generated for CMSSW Reference Manual by doxygen 1.8.5