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