CMS 3D CMS Logo
EGEnergyCorrector.cc
Go to the documentation of this file.
1 #include <TFile.h>
2 #include "RecoEgamma/EgammaTools/interface/EGEnergyCorrector.h"
3 #include "CondFormats/EgammaObjects/interface/GBRForest.h"
4 #include "CondFormats/DataRecord/interface/GBRWrapperRcd.h"
5 #include "FWCore/Framework/interface/ESHandle.h"
6 #include "DataFormats/EgammaCandidates/interface/Photon.h"
7 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
8 #include "DataFormats/Math/interface/deltaPhi.h"
9 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
10 #include "DataFormats/VertexReco/interface/Vertex.h"
11 #include "RecoEgamma/EgammaTools/interface/EcalClusterLocal.h"
12 #include "Geometry/Records/interface/CaloGeometryRecord.h"
13 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
14 
15 using namespace reco;
16 
17 //--------------------------------------------------------------------------------------------------
18 EGEnergyCorrector::EGEnergyCorrector() :
19 fReadereb(nullptr),
20 fReaderebvariance(nullptr),
21 fReaderee(nullptr),
22 fReadereevariance(nullptr),
23 fIsInitialized(kFALSE),
24 fOwnsForests(kFALSE),
25 fVals(nullptr)
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(regweights + "_EBCorrection", readereb);
68  iSetup.get<GBRWrapperRcd>().get(regweights + "_EBUncertainty", readerebvar);
69  iSetup.get<GBRWrapperRcd>().get(regweights + "_EECorrection", readeree);
70  iSetup.get<GBRWrapperRcd>().get(regweights + "_EEUncertainty", 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  edm::ESHandle<CaloGeometry> caloGeometry;
226  es.get<CaloGeometryRecord>().get(caloGeometry);
227  egammaTools::localEcalClusterCoordsEB(*b,*caloGeometry,betacry,bphicry,bieta,biphi,bthetatilt,bphitilt);
228 
229  fVals[56] = bieta; //crystal ieta
230  fVals[57] = biphi; //crystal iphi
231  fVals[58] = bieta%5; //submodule boundary eta symmetry
232  fVals[59] = biphi%2; //submodule boundary phi symmetry
233  fVals[60] = (TMath::Abs(bieta)<=25)*(bieta%25) + (TMath::Abs(bieta)>25)*((bieta-25*TMath::Abs(bieta)/bieta)%20); //module boundary eta approximate symmetry
234  fVals[61] = biphi%20; //module boundary phi symmetry
235  fVals[62] = betacry; //local coordinates with respect to closest crystal center at nominal shower depth
236  fVals[63] = bphicry;
237 
238 
239  //2nd cluster (meaningful gap corrections for converted photons)
240  float bc2etacry, bc2phicry, bc2thetatilt, bc2phitilt;
241  int bc2ieta, bc2iphi;
242  if (hasbc2) egammaTools::localEcalClusterCoordsEB( *b2,*caloGeometry,bc2etacry,bc2phicry,
243  bc2ieta,bc2iphi,bc2thetatilt,bc2phitilt );
244 
245  fVals[64] = hasbc2 ? bc2ieta : 0.;
246  fVals[65] = hasbc2 ? bc2iphi : 0.;
247  fVals[66] = hasbc2 ? bc2ieta%5 : 0.;
248  fVals[67] = hasbc2 ? bc2iphi%2 : 0.;
249  fVals[68] = hasbc2 ? (TMath::Abs(bc2ieta)<=25)*(bc2ieta%25) + (TMath::Abs(bc2ieta)>25)*((bc2ieta-25*TMath::Abs(bc2ieta)/bc2ieta)%20) : 0.;
250  fVals[69] = hasbc2 ? bc2iphi%20 : 0.;
251  fVals[70] = hasbc2 ? bc2etacry : 0.;
252  fVals[71] = hasbc2 ? bc2phicry : 0.;
253 
254  fVals[72] = vtxcol.size();
255 
256  }
257  else {
258  fVals[0] = s->rawEnergy();
259  fVals[1] = p.r9();
260  fVals[2] = s->eta();
261  fVals[3] = s->phi();
262  fVals[4] = p.e5x5()/s->rawEnergy();
263  fVals[5] = s->etaWidth();
264  fVals[6] = s->phiWidth();
265  fVals[7] = vtxcol.size();
266  }
267 
268 
269  const Double_t varscale = 1.253;
270  Double_t den;
271  const GBRForest *reader;
272  const GBRForest *readervar;
273  if (isbarrel) {
274  den = s->rawEnergy();
275  reader = fReadereb;
276  readervar = fReaderebvariance;
277  }
278  else {
279  den = s->rawEnergy() + s->preshowerEnergy();
280  reader = fReaderee;
281  readervar = fReadereevariance;
282  }
283 
284  Double_t ecor = reader->GetResponse(fVals)*den;
285  Double_t ecorerr = readervar->GetResponse(fVals)*den*varscale;
286 
287  //printf("ecor = %5f, ecorerr = %5f\n",ecor,ecorerr);
288 
289  return std::pair<double,double>(ecor,ecorerr);
290 }
291 
292 
293 //--------------------------------------------------------------------------------------------------
294 std::pair<double,double> EGEnergyCorrector::CorrectedEnergyWithError(const GsfElectron &e, const reco::VertexCollection& vtxcol, EcalClusterLazyTools &clustertools, const edm::EventSetup &es) {
295 
296  //apply v2 regression to electrons
297  //mostly duplicated from photon function above //TODO, make common underlying function
298 
299  //protection, this doesn't work properly on non-egamma-seeded electrons
300  if (!e.ecalDrivenSeed()) return std::pair<double,double>(0.,0.);
301 
302 
303  const SuperClusterRef s = e.superCluster();
304  const CaloClusterPtr b = s->seed();
305 
306  CaloClusterPtr b2;
307  Double_t ebcmax = -99.;
308  for (reco::CaloCluster_iterator bit = s->clustersBegin(); bit!=s->clustersEnd(); ++bit) {
309  const CaloClusterPtr bc = *bit;
310  if (bc->energy() > ebcmax && bc !=b) {
311  b2 = bc;
312  ebcmax = bc->energy();
313  }
314  }
315 
316  CaloClusterPtr bclast;
317  Double_t ebcmin = 1e6;
318  for (reco::CaloCluster_iterator bit = s->clustersBegin(); bit!=s->clustersEnd(); ++bit) {
319  const CaloClusterPtr bc = *bit;
320  if (bc->energy() < ebcmin && bc !=b) {
321  bclast = bc;
322  ebcmin = bc->energy();
323  }
324  }
325 
326  CaloClusterPtr bclast2;
327  ebcmin = 1e6;
328  for (reco::CaloCluster_iterator bit = s->clustersBegin(); bit!=s->clustersEnd(); ++bit) {
329  const CaloClusterPtr bc = *bit;
330  if (bc->energy() < ebcmin && bc !=b && bc!=bclast) {
331  bclast2 = bc;
332  ebcmin = bc->energy();
333  }
334  }
335 
336  Bool_t isbarrel = b->hitsAndFractions().at(0).first.subdetId()==EcalBarrel;
337  Bool_t hasbc2 = b2.isNonnull() && b2->energy()>0.;
338  Bool_t hasbclast = bclast.isNonnull() && bclast->energy()>0.;
339  Bool_t hasbclast2 = bclast2.isNonnull() && bclast2->energy()>0.;
340 
341 
342  if (isbarrel) {
343  fVals[0] = s->rawEnergy();
344  fVals[1] = clustertools.e3x3(*b)/s->rawEnergy(); //r9
345  fVals[2] = s->eta();
346  fVals[3] = s->phi();
347  fVals[4] = clustertools.e5x5(*b)/s->rawEnergy();
348  fVals[5] = e.hcalOverEcal();
349  fVals[6] = s->etaWidth();
350  fVals[7] = s->phiWidth();
351 
352 
353  double bemax = clustertools.eMax(*b);
354  double be2nd = clustertools.e2nd(*b);
355  double betop = clustertools.eTop(*b);
356  double bebottom = clustertools.eBottom(*b);
357  double beleft = clustertools.eLeft(*b);
358  double beright = clustertools.eRight(*b);
359 
360 
361  fVals[8] = b->eta()-s->eta();
362  fVals[9] = reco::deltaPhi(b->phi(),s->phi());
363  fVals[10] = b->energy()/s->rawEnergy();
364  fVals[11] = clustertools.e3x3(*b)/b->energy();
365  fVals[12] = clustertools.e5x5(*b)/b->energy();
366  fVals[13] = sqrt(clustertools.localCovariances(*b)[0]);
367  fVals[14] = sqrt(clustertools.localCovariances(*b)[2]);
368  fVals[15] = clustertools.localCovariances(*b)[1];
369  fVals[16] = bemax/b->energy();
370  fVals[17] = log(be2nd/bemax);
371  fVals[18] = log(betop/bemax);
372  fVals[19] = log(bebottom/bemax);
373  fVals[20] = log(beleft/bemax);
374  fVals[21] = log(beright/bemax);
375  fVals[22] = (betop-bebottom)/(betop+bebottom);
376  fVals[23] = (beleft-beright)/(beleft+beright);
377 
378 
379  double bc2emax = hasbc2 ? clustertools.eMax(*b2) : 0.;
380  double bc2e2nd = hasbc2 ? clustertools.e2nd(*b2) : 0.;
381  double bc2etop = hasbc2 ? clustertools.eTop(*b2) : 0.;
382  double bc2ebottom = hasbc2 ? clustertools.eBottom(*b2) : 0.;
383  double bc2eleft = hasbc2 ? clustertools.eLeft(*b2) : 0.;
384  double bc2eright = hasbc2 ? clustertools.eRight(*b2) : 0.;
385 
386  fVals[24] = hasbc2 ? (b2->eta()-s->eta()) : 0.;
387  fVals[25] = hasbc2 ? reco::deltaPhi(b2->phi(),s->phi()) : 0.;
388  fVals[26] = hasbc2 ? b2->energy()/s->rawEnergy() : 0.;
389  fVals[27] = hasbc2 ? clustertools.e3x3(*b2)/b2->energy() : 0.;
390  fVals[28] = hasbc2 ? clustertools.e5x5(*b2)/b2->energy() : 0.;
391  fVals[29] = hasbc2 ? sqrt(clustertools.localCovariances(*b2)[0]) : 0.;
392  fVals[30] = hasbc2 ? sqrt(clustertools.localCovariances(*b2)[2]) : 0.;
393  fVals[31] = hasbc2 ? clustertools.localCovariances(*b)[1] : 0.;
394  fVals[32] = hasbc2 ? bc2emax/b2->energy() : 0.;
395  fVals[33] = hasbc2 ? log(bc2e2nd/bc2emax) : 0.;
396  fVals[34] = hasbc2 ? log(bc2etop/bc2emax) : 0.;
397  fVals[35] = hasbc2 ? log(bc2ebottom/bc2emax) : 0.;
398  fVals[36] = hasbc2 ? log(bc2eleft/bc2emax) : 0.;
399  fVals[37] = hasbc2 ? log(bc2eright/bc2emax) : 0.;
400  fVals[38] = hasbc2 ? (bc2etop-bc2ebottom)/(bc2etop+bc2ebottom) : 0.;
401  fVals[39] = hasbc2 ? (bc2eleft-bc2eright)/(bc2eleft+bc2eright) : 0.;
402 
403  fVals[40] = hasbclast ? (bclast->eta()-s->eta()) : 0.;
404  fVals[41] = hasbclast ? reco::deltaPhi(bclast->phi(),s->phi()) : 0.;
405  fVals[42] = hasbclast ? bclast->energy()/s->rawEnergy() : 0.;
406  fVals[43] = hasbclast ? clustertools.e3x3(*bclast)/bclast->energy() : 0.;
407  fVals[44] = hasbclast ? clustertools.e5x5(*bclast)/bclast->energy() : 0.;
408  fVals[45] = hasbclast ? sqrt(clustertools.localCovariances(*bclast)[0]) : 0.;
409  fVals[46] = hasbclast ? sqrt(clustertools.localCovariances(*bclast)[2]) : 0.;
410  fVals[47] = hasbclast ? clustertools.localCovariances(*bclast)[1] : 0.;
411 
412  fVals[48] = hasbclast2 ? (bclast2->eta()-s->eta()) : 0.;
413  fVals[49] = hasbclast2 ? reco::deltaPhi(bclast2->phi(),s->phi()) : 0.;
414  fVals[50] = hasbclast2 ? bclast2->energy()/s->rawEnergy() : 0.;
415  fVals[51] = hasbclast2 ? clustertools.e3x3(*bclast2)/bclast2->energy() : 0.;
416  fVals[52] = hasbclast2 ? clustertools.e5x5(*bclast2)/bclast2->energy() : 0.;
417  fVals[53] = hasbclast2 ? sqrt(clustertools.localCovariances(*bclast2)[0]) : 0.;
418  fVals[54] = hasbclast2 ? sqrt(clustertools.localCovariances(*bclast2)[2]) : 0.;
419  fVals[55] = hasbclast2 ? clustertools.localCovariances(*bclast2)[1] : 0.;
420 
421 
422  float betacry, bphicry, bthetatilt, bphitilt;
423  int bieta, biphi;
424  edm::ESHandle<CaloGeometry> caloGeometry;
425  es.get<CaloGeometryRecord>().get(caloGeometry);
426  egammaTools::localEcalClusterCoordsEB(*b,*caloGeometry,betacry,bphicry,bieta,biphi,bthetatilt,bphitilt);
427 
428  fVals[56] = bieta;
429  fVals[57] = biphi;
430  fVals[58] = bieta%5;
431  fVals[59] = biphi%2;
432  fVals[60] = (TMath::Abs(bieta)<=25)*(bieta%25) + (TMath::Abs(bieta)>25)*((bieta-25*TMath::Abs(bieta)/bieta)%20);
433  fVals[61] = biphi%20;
434  fVals[62] = betacry;
435  fVals[63] = bphicry;
436 
437  float bc2etacry, bc2phicry, bc2thetatilt, bc2phitilt;
438  int bc2ieta, bc2iphi;
439  if (hasbc2) egammaTools::localEcalClusterCoordsEB( *b2,*caloGeometry,bc2etacry,bc2phicry,
440  bc2ieta,bc2iphi,bc2thetatilt,bc2phitilt );
441 
442  fVals[64] = hasbc2 ? bc2ieta : 0.;
443  fVals[65] = hasbc2 ? bc2iphi : 0.;
444  fVals[66] = hasbc2 ? bc2ieta%5 : 0.;
445  fVals[67] = hasbc2 ? bc2iphi%2 : 0.;
446  fVals[68] = hasbc2 ? (TMath::Abs(bc2ieta)<=25)*(bc2ieta%25) + (TMath::Abs(bc2ieta)>25)*((bc2ieta-25*TMath::Abs(bc2ieta)/bc2ieta)%20) : 0.;
447  fVals[69] = hasbc2 ? bc2iphi%20 : 0.;
448  fVals[70] = hasbc2 ? bc2etacry : 0.;
449  fVals[71] = hasbc2 ? bc2phicry : 0.;
450 
451  fVals[72] = vtxcol.size();
452 
453  }
454  else {
455  fVals[0] = s->rawEnergy();
456  fVals[1] = clustertools.e3x3(*b)/s->rawEnergy(); //r9
457  fVals[2] = s->eta();
458  fVals[3] = s->phi();
459  fVals[4] = clustertools.e5x5(*b)/s->rawEnergy();
460  fVals[5] = s->etaWidth();
461  fVals[6] = s->phiWidth();
462  fVals[7] = vtxcol.size();
463  }
464 
465 
466  const Double_t varscale = 1.253;
467  Double_t den;
468  const GBRForest *reader;
469  const GBRForest *readervar;
470  if (isbarrel) {
471  den = s->rawEnergy();
472  reader = fReadereb;
473  readervar = fReaderebvariance;
474  }
475  else {
476  den = s->rawEnergy() + s->preshowerEnergy();
477  reader = fReaderee;
478  readervar = fReadereevariance;
479  }
480 
481  Double_t ecor = reader->GetResponse(fVals)*den;
482  Double_t ecorerr = readervar->GetResponse(fVals)*den*varscale;
483 
484  //printf("ecor = %5f, ecorerr = %5f\n",ecor,ecorerr);
485 
486  return std::pair<double,double>(ecor,ecorerr);
487 
488 }
489 
490 //--------------------------------------------------------------------------------------------------
491 std::pair<double,double> EGEnergyCorrector::CorrectedEnergyWithErrorV3(const Photon &p, const reco::VertexCollection& vtxcol, double rho, EcalClusterLazyTools &clustertools, const edm::EventSetup &es, bool applyRescale) {
492 
493  const SuperClusterRef s = p.superCluster();
494  const CaloClusterPtr b = s->seed(); //seed basic cluster
495 
496  Bool_t isbarrel = b->hitsAndFractions().at(0).first.subdetId()==EcalBarrel;
497 
498  //basic supercluster variables
499  fVals[0] = s->rawEnergy();
500  fVals[1] = s->eta();
501  fVals[2] = s->phi();
502  fVals[3] = p.r9();
503  fVals[4] = p.e5x5()/s->rawEnergy();
504  fVals[5] = s->etaWidth();
505  fVals[6] = s->phiWidth();
506  fVals[7] = s->clustersSize();
507  fVals[8] = p.hadTowOverEm();
508  fVals[9] = rho;
509  fVals[10] = vtxcol.size();
510 
511  //seed basic cluster variables
512  double bemax = clustertools.eMax(*b);
513  double be2nd = clustertools.e2nd(*b);
514  double betop = clustertools.eTop(*b);
515  double bebottom = clustertools.eBottom(*b);
516  double beleft = clustertools.eLeft(*b);
517  double beright = clustertools.eRight(*b);
518 
519  double be2x5max = clustertools.e2x5Max(*b);
520  double be2x5top = clustertools.e2x5Top(*b);
521  double be2x5bottom = clustertools.e2x5Bottom(*b);
522  double be2x5left = clustertools.e2x5Left(*b);
523  double be2x5right = clustertools.e2x5Right(*b);
524 
525  fVals[11] = b->eta()-s->eta();
526  fVals[12] = reco::deltaPhi(b->phi(),s->phi());
527  fVals[13] = b->energy()/s->rawEnergy();
528  fVals[14] = clustertools.e3x3(*b)/b->energy();
529  fVals[15] = clustertools.e5x5(*b)/b->energy();
530  fVals[16] = sqrt(clustertools.localCovariances(*b)[0]); //sigietaieta
531  fVals[17] = sqrt(clustertools.localCovariances(*b)[2]); //sigiphiiphi
532  fVals[18] = clustertools.localCovariances(*b)[1]; //sigietaiphi
533  fVals[19] = bemax/b->energy(); //crystal energy ratio gap variables
534  fVals[20] = be2nd/b->energy();
535  fVals[21] = betop/b->energy();
536  fVals[22] = bebottom/b->energy();
537  fVals[23] = beleft/b->energy();
538  fVals[24] = beright/b->energy();
539  fVals[25] = be2x5max/b->energy(); //crystal energy ratio gap variables
540  fVals[26] = be2x5top/b->energy();
541  fVals[27] = be2x5bottom/b->energy();
542  fVals[28] = be2x5left/b->energy();
543  fVals[29] = be2x5right/b->energy();
544 
545  if (isbarrel) {
546  //local coordinates and crystal indices (barrel only)
547 
548  //seed cluster
549  float betacry, bphicry, bthetatilt, bphitilt;
550  int bieta, biphi;
551  edm::ESHandle<CaloGeometry> caloGeometry;
552  es.get<CaloGeometryRecord>().get(caloGeometry);
553  egammaTools::localEcalClusterCoordsEB(*b,*caloGeometry,betacry,bphicry,bieta,biphi,bthetatilt,bphitilt);
554 
555  fVals[30] = bieta; //crystal ieta
556  fVals[31] = biphi; //crystal iphi
557  fVals[32] = bieta%5; //submodule boundary eta symmetry
558  fVals[33] = biphi%2; //submodule boundary phi symmetry
559  fVals[34] = (TMath::Abs(bieta)<=25)*(bieta%25) + (TMath::Abs(bieta)>25)*((bieta-25*TMath::Abs(bieta)/bieta)%20); //module boundary eta approximate symmetry
560  fVals[35] = biphi%20; //module boundary phi symmetry
561  fVals[36] = betacry; //local coordinates with respect to closest crystal center at nominal shower depth
562  fVals[37] = bphicry;
563 
564  }
565  else {
566  //preshower energy ratio (endcap only)
567  fVals[30] = s->preshowerEnergy()/s->rawEnergy();
568  }
569 
570 // if (isbarrel) {
571 // for (int i=0; i<38; ++i) printf("%i: %5f\n",i,fVals[i]);
572 // }
573 // else for (int i=0; i<31; ++i) printf("%i: %5f\n",i,fVals[i]);
574 
575  Double_t den;
576  const GBRForest *reader;
577  const GBRForest *readervar;
578  if (isbarrel) {
579  den = s->rawEnergy();
580  reader = fReadereb;
581  readervar = fReaderebvariance;
582  }
583  else {
584  den = s->rawEnergy() + s->preshowerEnergy();
585  reader = fReaderee;
586  readervar = fReadereevariance;
587  }
588 
589  Double_t ecor = reader->GetResponse(fVals)*den;
590 
591 
592  //apply shower shape rescaling - for Monte Carlo only, and only for calculation of energy uncertainty
593  if (applyRescale) {
594  if (isbarrel) {
595  fVals[3] = 1.0045*p.r9() +0.001; //r9
596  fVals[5] = 1.04302*s->etaWidth() - 0.000618; //etawidth
597  fVals[6] = 1.00002*s->phiWidth() - 0.000371; //phiwidth
598  fVals[14] = fVals[3]*s->rawEnergy()/b->energy(); //compute consistent e3x3/eseed after r9 rescaling
599  if (fVals[15]<=1.0) // rescale e5x5/eseed only if value is <=1.0, don't allow scaled values to exceed 1.0
600  fVals[15] = TMath::Min(1.0,1.0022*p.e5x5()/b->energy());
601 
602  fVals[4] = fVals[15]*b->energy()/s->rawEnergy(); // compute consistent e5x5()/rawEnergy() after e5x5/eseed resacling
603 
604  fVals[16] = 0.891832*sqrt(clustertools.localCovariances(*b)[0]) + 0.0009133; //sigietaieta
605  fVals[17] = 0.993*sqrt(clustertools.localCovariances(*b)[2]); //sigiphiiphi
606 
607  fVals[19] = 1.012*bemax/b->energy(); //crystal energy ratio gap variables
608  fVals[20] = 1.0*be2nd/b->energy();
609  fVals[21] = 0.94*betop/b->energy();
610  fVals[22] = 0.94*bebottom/b->energy();
611  fVals[23] = 0.94*beleft/b->energy();
612  fVals[24] = 0.94*beright/b->energy();
613  fVals[25] = 1.006*be2x5max/b->energy(); //crystal energy ratio gap variables
614  fVals[26] = 1.09*be2x5top/b->energy();
615  fVals[27] = 1.09*be2x5bottom/b->energy();
616  fVals[28] = 1.09*be2x5left/b->energy();
617  fVals[29] = 1.09*be2x5right/b->energy();
618 
619  }
620  else {
621  fVals[3] = 1.0086*p.r9() -0.0007; //r9
622  fVals[4] = TMath::Min(1.0,1.0022*p.e5x5()/s->rawEnergy()); //e5x5/rawenergy
623  fVals[5] = 0.903254*s->etaWidth() +0.001346; //etawidth
624  fVals[6] = 0.99992*s->phiWidth() + 4.8e-07; //phiwidth
625  fVals[13] = TMath::Min(1.0,1.0022*b->energy()/s->rawEnergy()); //eseed/rawenergy (practically equivalent to e5x5)
626 
627 
628  fVals[14] = fVals[3]*s->rawEnergy()/b->energy(); //compute consistent e3x3/eseed after r9 rescaling
629 
630  fVals[16] = 0.9947*sqrt(clustertools.localCovariances(*b)[0]) + 0.00003; //sigietaieta
631 
632  fVals[19] = 1.005*bemax/b->energy(); //crystal energy ratio gap variables
633  fVals[20] = 1.02*be2nd/b->energy();
634  fVals[21] = 0.96*betop/b->energy();
635  fVals[22] = 0.96*bebottom/b->energy();
636  fVals[23] = 0.96*beleft/b->energy();
637  fVals[24] = 0.96*beright/b->energy();
638  fVals[25] = 1.0075*be2x5max/b->energy(); //crystal energy ratio gap variables
639  fVals[26] = 1.13*be2x5top/b->energy();
640  fVals[27] = 1.13*be2x5bottom/b->energy();
641  fVals[28] = 1.13*be2x5left/b->energy();
642  fVals[29] = 1.13*be2x5right/b->energy();
643  }
644 
645  }
646 
647  Double_t ecorerr = readervar->GetResponse(fVals)*den;
648 
649  //printf("ecor = %5f, ecorerr = %5f\n",ecor,ecorerr);
650 
651  return std::pair<double,double>(ecor,ecorerr);
652 }
653 
654 //--------------------------------------------------------------------------------------------------
655 std::pair<double,double> EGEnergyCorrector::CorrectedEnergyWithErrorV3(const GsfElectron &e, const reco::VertexCollection& vtxcol, double rho, EcalClusterLazyTools &clustertools, const edm::EventSetup &es) {
656 
657  const SuperClusterRef s = e.superCluster();
658  const CaloClusterPtr b = s->seed(); //seed basic cluster
659 
660  Bool_t isbarrel = b->hitsAndFractions().at(0).first.subdetId()==EcalBarrel;
661 
662  //basic supercluster variables
663  fVals[0] = s->rawEnergy();
664  fVals[1] = s->eta();
665  fVals[2] = s->phi();
666  fVals[3] = clustertools.e3x3(*b)/s->rawEnergy(); //r9
667  fVals[4] = clustertools.e5x5(*b)/s->rawEnergy();
668  fVals[5] = s->etaWidth();
669  fVals[6] = s->phiWidth();
670  fVals[7] = s->clustersSize();
671  fVals[8] = e.hcalOverEcalBc();
672  fVals[9] = rho;
673  fVals[10] = vtxcol.size();
674 
675  //seed basic cluster variables
676  double bemax = clustertools.eMax(*b);
677  double be2nd = clustertools.e2nd(*b);
678  double betop = clustertools.eTop(*b);
679  double bebottom = clustertools.eBottom(*b);
680  double beleft = clustertools.eLeft(*b);
681  double beright = clustertools.eRight(*b);
682 
683  double be2x5max = clustertools.e2x5Max(*b);
684  double be2x5top = clustertools.e2x5Top(*b);
685  double be2x5bottom = clustertools.e2x5Bottom(*b);
686  double be2x5left = clustertools.e2x5Left(*b);
687  double be2x5right = clustertools.e2x5Right(*b);
688 
689  fVals[11] = b->eta()-s->eta();
690  fVals[12] = reco::deltaPhi(b->phi(),s->phi());
691  fVals[13] = b->energy()/s->rawEnergy();
692  fVals[14] = clustertools.e3x3(*b)/b->energy();
693  fVals[15] = clustertools.e5x5(*b)/b->energy();
694  fVals[16] = sqrt(clustertools.localCovariances(*b)[0]); //sigietaieta
695  fVals[17] = sqrt(clustertools.localCovariances(*b)[2]); //sigiphiiphi
696  fVals[18] = clustertools.localCovariances(*b)[1]; //sigietaiphi
697  fVals[19] = bemax/b->energy(); //crystal energy ratio gap variables
698  fVals[20] = be2nd/b->energy();
699  fVals[21] = betop/b->energy();
700  fVals[22] = bebottom/b->energy();
701  fVals[23] = beleft/b->energy();
702  fVals[24] = beright/b->energy();
703  fVals[25] = be2x5max/b->energy(); //crystal energy ratio gap variables
704  fVals[26] = be2x5top/b->energy();
705  fVals[27] = be2x5bottom/b->energy();
706  fVals[28] = be2x5left/b->energy();
707  fVals[29] = be2x5right/b->energy();
708 
709  if (isbarrel) {
710  //local coordinates and crystal indices (barrel only)
711 
712  //seed cluster
713  float betacry, bphicry, bthetatilt, bphitilt;
714  int bieta, biphi;
715  edm::ESHandle<CaloGeometry> caloGeometry;
716  es.get<CaloGeometryRecord>().get(caloGeometry);
717  egammaTools::localEcalClusterCoordsEB(*b,*caloGeometry,betacry,bphicry,bieta,biphi,bthetatilt,bphitilt);
718 
719  fVals[30] = bieta; //crystal ieta
720  fVals[31] = biphi; //crystal iphi
721  fVals[32] = bieta%5; //submodule boundary eta symmetry
722  fVals[33] = biphi%2; //submodule boundary phi symmetry
723  fVals[34] = (TMath::Abs(bieta)<=25)*(bieta%25) + (TMath::Abs(bieta)>25)*((bieta-25*TMath::Abs(bieta)/bieta)%20); //module boundary eta approximate symmetry
724  fVals[35] = biphi%20; //module boundary phi symmetry
725  fVals[36] = betacry; //local coordinates with respect to closest crystal center at nominal shower depth
726  fVals[37] = bphicry;
727 
728  }
729  else {
730  //preshower energy ratio (endcap only)
731  fVals[30] = s->preshowerEnergy()/s->rawEnergy();
732  }
733 
734  Double_t den;
735  const GBRForest *reader;
736  const GBRForest *readervar;
737  if (isbarrel) {
738  den = s->rawEnergy();
739  reader = fReadereb;
740  readervar = fReaderebvariance;
741  }
742  else {
743  den = s->rawEnergy() + s->preshowerEnergy();
744  reader = fReaderee;
745  readervar = fReadereevariance;
746  }
747 
748  Double_t ecor = reader->GetResponse(fVals)*den;
749  Double_t ecorerr = readervar->GetResponse(fVals)*den;
750 
751  //printf("ecor = %5f, ecorerr = %5f\n",ecor,ecorerr);
752 
753  return std::pair<double,double>(ecor,ecorerr);
754 }
755 
GBRForest::GetResponse
double GetResponse(const float *vector) const
Definition: GBRForest.h:53
reco::deltaPhi
constexpr double deltaPhi(double phi1, double phi2)
Definition: deltaPhi.h:22
egammaTools::localEcalClusterCoordsEB
void localEcalClusterCoordsEB(const reco::CaloCluster &bclus, const CaloGeometry &geom, float &etacry, float &phicry, int &ieta, int &iphi, float &thetatilt, float &phitilt)
Definition: EcalClusterLocal.cc:14
Photon
Definition: Photon.py:1
EGEnergyCorrector::Initialize
void Initialize(const edm::EventSetup &iSetup, std::string regweights, bool weightsFromDB=false)
Definition: EGEnergyCorrector.cc:47
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
reco::Photon::e5x5
float e5x5() const
Definition: Photon.h:234
EGEnergyCorrector::fReaderee
const GBRForest * fReaderee
Definition: EGEnergyCorrector.h:38
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
EGEnergyCorrector::fReadereb
const GBRForest * fReadereb
Definition: EGEnergyCorrector.h:36
EGEnergyCorrector::fReaderebvariance
const GBRForest * fReaderebvariance
Definition: EGEnergyCorrector.h:37
nullptr
#define nullptr
Definition: GCC11Compatibility.h:37
reco::VertexCollection
std::vector< Vertex > VertexCollection
collection of Vertex objects
Definition: VertexFwd.h:9
reco::CaloCluster::hitsAndFractions
const std::vector< std::pair< DetId, float > > & hitsAndFractions() const
Definition: CaloCluster.h:197
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
Min
T Min(T a, T b)
Definition: MathUtil.h:39
reco::Photon::superCluster
reco::SuperClusterRef superCluster() const override
Ref to SuperCluster.
reco::CaloCluster::eta
double eta() const
pseudorapidity of cluster centroid
Definition: CaloCluster.h:168
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
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
Abs
T Abs(T a)
Definition: MathUtil.h:49
reco::GsfElectron::hcalOverEcal
float hcalOverEcal() const
Definition: GsfElectron.h:448
reco::CaloCluster::energy
double energy() const
cluster energy
Definition: CaloCluster.h:126
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:221
reco::Photon::hadronicOverEm
float hadronicOverEm() const
the total hadronic over electromagnetic fraction
Definition: Photon.h:212
reco::GsfElectron::hcalOverEcalBc
float hcalOverEcalBc() const
Definition: GsfElectron.h:452
edm::Ptr::isNonnull
bool isNonnull() const
Checks for non-null.
Definition: Ptr.h:168
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:491
b
double b
Definition: hdecay.h:120
reco
fixed size matrix
Definition: AlignmentAlgorithmBase.h:44
EGEnergyCorrector::fReadereevariance
const GBRForest * fReadereevariance
Definition: EGEnergyCorrector.h:39
edm::EventSetup::get
T get() const
Definition: EventSetup.h:71
reco::Photon::r9
float r9() const
Definition: Photon.h:240
reco::GsfElectron::superCluster
SuperClusterRef superCluster() const override
reference to a SuperCluster
Definition: GsfElectron.h:185
reco::CaloCluster::phi
double phi() const
azimuthal angle of cluster centroid
Definition: CaloCluster.h:171
edm::ESHandle::product
T const * product() const
Definition: ESHandle.h:86
reco::GsfElectron::ecalDrivenSeed
bool ecalDrivenSeed() const
Definition: GsfElectron.h:188
Generated for CMSSW Reference Manual by doxygen 1.8.11