CMS 3D CMS Logo

ElectronEnergyCorrector.cc
Go to the documentation of this file.
1 
7 #include "TMath.h"
8 
9 
10 /****************************************************************************
11  *
12  * Classification based eta corrections for the ecal cluster energy
13  *
14  * \author Federico Ferri - INFN Milano, Bicocca university
15  * \author Ivica Puljak - FESB, Split
16  * \author Stephanie Baffioni - Laboratoire Leprince-Ringuet - École polytechnique, CNRS/IN2P3
17  *
18  ****************************************************************************/
19 
21  {
23  double energyError = 999. ;
24  double ecalEnergy = electron.correctedEcalEnergy() ;
25  double eleEta = electron.superCluster()->eta() ;
26  double brem = electron.superCluster()->etaWidth()/electron.superCluster()->phiWidth() ;
27  energyError = uncertainty.computeElectronEnergyUncertainty(electron.classification(),eleEta,brem,ecalEnergy) ;
28  electron.setCorrectedEcalEnergyError(energyError) ;
29  }
30 
31 float energyError( float E, float * par )
32  { return sqrt( pow(par[0]/sqrt(E),2) + pow(par[1]/E,2) + pow(par[2],2) ) ; }
33 
35  {
36  double error = 999. ;
37  double ecalEnergy = electron.correctedEcalEnergy() ;
38 
39  if (electron.isEB())
40  {
41  float parEB[3] = { 5.24e-02, 2.01e-01, 1.00e-02} ;
42  error = ecalEnergy*energyError(ecalEnergy,parEB) ;
43  }
44  else if (electron.isEE())
45  {
46  float parEE[3] = { 1.46e-01, 9.21e-01, 1.94e-03} ;
47  error = ecalEnergy*energyError(ecalEnergy,parEE) ;
48  }
49  else
50  { edm::LogWarning("ElectronEnergyCorrector::simpleParameterizationUncertainty")<<"nor barrel neither endcap electron !" ; }
51 
52  electron.setCorrectedEcalEnergyError(error) ;
53  }
54 
57  {
58  if (electron.isEcalEnergyCorrected())
59  {
60  edm::LogWarning("ElectronEnergyCorrector::classBasedParameterizationEnergy")<<"already done" ;
61  return ;
62  }
63 
65  if ( (elClass <= reco::GsfElectron::UNKNOWN) ||
66  (elClass>reco::GsfElectron::GAP) )
67  {
68  edm::LogWarning("ElectronEnergyCorrector::classBasedParameterizationEnergy")<<"unexpected classification" ;
69  return ;
70  }
71 
72  // new corrections from N. Chanon et al., taken from EcalClusterCorrectionObjectSpecific.cc
73  float corr = 1.;
74  float corr2 = 1.;
75  float energy = electron.superCluster()->energy() ;
76 
77  //int subdet = electron.superCluster()->seed()->hitsAndFractions()[0].first.subdetId();
78 
79  if (electron.isEB())
80  {
81  float cetacorr = fEta(electron.superCluster()->rawEnergy(), electron.superCluster()->eta(), 0)/electron.superCluster()->rawEnergy();
82  energy = electron.superCluster()->rawEnergy()*cetacorr; //previously in CMSSW
83  //energy = superCluster.rawEnergy()*fEta(e5x5, superCluster.seed()->eta(), 0)/e5x5;
84  }
85  else if (electron.isEE())
86  {
87  energy = electron.superCluster()->rawEnergy()+electron.superCluster()->preshowerEnergy();
88  }
89  else
90  { edm::LogWarning("ElectronEnergyCorrector::classBasedParameterizationEnergy")<<"nor barrel neither endcap electron !" ; }
91 
92  corr = fBremEta(electron.superCluster()->phiWidth()/electron.superCluster()->etaWidth(), electron.superCluster()->eta(), 0,elClass);
93 
94  float et = energy*TMath::Sin(2*TMath::ATan(TMath::Exp(-electron.superCluster()->eta())))/corr;
95 
96  if (electron.isEB()) { corr2 = corr * fEt(et, 0,elClass) ; }
97  else if (electron.isEE()) { corr2 = corr * fEnergy(energy/corr, 1,elClass) ; }
98  else { edm::LogWarning("ElectronEnergyCorrector::classBasedParameterizationEnergy")<<"nor barrel neither endcap electron !" ; }
99 
100  float newEnergy = energy/corr2;
101 
102  // cracks
103  double crackcor = 1. ;
105  cIt = electron.superCluster()->clustersBegin() ;
106  cIt != electron.superCluster()->clustersEnd() ;
107  ++cIt )
108  {
109  const reco::CaloClusterPtr cc = *cIt ;
110  crackcor *=
111  ( electron.superCluster()->rawEnergy()
112  + cc->energy()*(crackCorrectionFunction_->getValue(*cc)-1.) )
113  / electron.superCluster()->rawEnergy() ;
114  }
115  newEnergy *= crackcor ;
116 
117  // register final value
118  electron.setCorrectedEcalEnergy(newEnergy) ;
119 
120  }
121 
122 // main correction function
123 // new corrections: taken from EcalClusterCorrectionObjectSpecific.cc (N. Chanon et al.)
124 // this is to prepare for class based corrections, for the time being the parameters are the same as for the SC corrections
125 // code fully duplicated here, to be improved; electron means algorithm==0 and mode==0
126 
127 float ElectronEnergyCorrector::fEta(float energy, float eta, int algorithm) const
128 {
129 
130  // corrections for electrons
131  if (algorithm!=0) {
132  edm::LogWarning("ElectronEnergyCorrector::fEta")<<"algorithm should be 0 for electrons !" ;
133  return energy;
134  }
135  //std::cout << "fEta function" << std::endl;
136 
137  // this correction is setup only for EB
138  if ( algorithm != 0 ) return energy;
139 
140  float ieta = fabs(eta)*(5/0.087);
141  // bypass the DB reading for the time being
142  //float p0 = (params_->params())[0]; // should be 40.2198
143  //float p1 = (params_->params())[1]; // should be -3.03103e-6
144  float p0 = 40.2198 ;
145  float p1 = -3.03103e-6 ;
146 
147  //std::cout << "ieta=" << ieta << std::endl;
148 
149  float correctedEnergy = energy;
150  if ( ieta < p0 ) correctedEnergy = energy;
151  else correctedEnergy = energy/(1.0 + p1*(ieta-p0)*(ieta-p0));
152  //std::cout << "ECEC fEta = " << correctedEnergy << std::endl;
153  return correctedEnergy;
154 
155  }
156 
158  ( float sigmaPhiSigmaEta, float eta, int algorithm, reco::GsfElectron::Classification cl ) const
159  {
160  // corrections for electrons
161  if (algorithm!=0)
162  {
163  edm::LogWarning("ElectronEnergyCorrector::fBremEta")<<"algorithm should be 0 for electrons !" ;
164  return 1. ;
165  }
166 
167  const float etaCrackMin = 1.44 ;
168  const float etaCrackMax = 1.56 ;
169 
170  //STD
171  const int nBinsEta = 14 ;
172  float leftEta[nBinsEta] = { 0.02, 0.25, 0.46, 0.81, 0.91, 1.01, 1.16, etaCrackMax, 1.653, 1.8, 2.0, 2.2, 2.3, 2.4 } ;
173  float rightEta[nBinsEta] = { 0.25, 0.42, 0.77, 0.91, 1.01, 1.13, etaCrackMin, 1.653, 1.8, 2.0, 2.2, 2.3, 2.4, 2.5 } ;
174 
175  // eta = 0
176  if ( TMath::Abs(eta) < leftEta[0] ) { eta = 0.02 ; }
177 
178  // outside acceptance
179  if ( TMath::Abs(eta) >= rightEta[nBinsEta-1] ) { eta = 2.49 ; } //if (DBG) std::cout << " WARNING [applyScCorrections]: TMath::Abs(eta) >= rightEta[nBinsEta-1] " << std::endl;}
180 
181  int tmpEta = -1 ;
182  for (int iEta = 0; iEta < nBinsEta; ++iEta)
183  {
184  if ( leftEta[iEta] <= TMath::Abs(eta) && TMath::Abs(eta) <rightEta[iEta] )
185  { tmpEta = iEta ; }
186  }
187 
188  float xcorr[nBinsEta][reco::GsfElectron::GAP+1] =
189  {
190  { 1.00227, 1.00227, 1.00227, 1.00227, 1.00227 },
191  { 1.00252, 1.00252, 1.00252, 1.00252, 1.00252 },
192  { 1.00225, 1.00225, 1.00225, 1.00225, 1.00225 },
193  { 1.00159, 1.00159, 1.00159, 1.00159, 1.00159 },
194  { 0.999475, 0.999475, 0.999475, 0.999475, 0.999475 },
195  { 0.997203, 0.997203, 0.997203, 0.997203, 0.997203 },
196  { 0.993886, 0.993886, 0.993886, 0.993886, 0.993886 },
197  { 0.971262, 0.971262, 0.971262, 0.971262, 0.971262 },
198  { 0.975922, 0.975922, 0.975922, 0.975922, 0.975922 },
199  { 0.979087, 0.979087, 0.979087, 0.979087, 0.979087 },
200  { 0.98495, 0.98495, 0.98495, 0.98495, 0.98495 },
201  { 0.98781, 0.98781, 0.98781, 0.98781, 0.98781 },
202  { 0.989546, 0.989546, 0.989546, 0.989546, 0.989546 },
203  { 0.989638, 0.989638, 0.989638, 0.989638, 0.989638 }
204  } ;
205 
206  float par0[nBinsEta][reco::GsfElectron::GAP+1] =
207  {
208  { 0.994949, 1.00718, 1.00718, 1.00556, 1.00718 },
209  { 1.009, 1.00713, 1.00713, 1.00248, 1.00713 },
210  { 0.999395, 1.00641, 1.00641, 1.00293, 1.00641 },
211  { 0.988662, 1.00761, 1.001, 0.99972, 1.00761 },
212  { 0.998443, 1.00682, 1.001, 1.00282, 1.00682 },
213  { 1.00285, 1.0073, 1.001, 1.00396, 1.0073 },
214  { 0.993053, 1.00462, 1.01341, 1.00184, 1.00462 },
215  { 1.10561, 0.972798, 1.02835, 0.995218, 0.972798 },
216  { 0.893741, 0.981672, 0.98982, 1.01712, 0.981672 },
217  { 0.911123, 0.98251, 1.03466, 1.00824, 0.98251 },
218  { 0.981931, 0.986123, 0.954295, 1.0202, 0.986123 },
219  { 0.905634, 0.990124, 0.928934, 0.998492, 0.990124 },
220  { 0.919343, 0.990187, 0.967526, 0.963923, 0.990187 },
221  { 0.844783, 0.99372, 0.923808, 0.953001, 0.99372 }
222  } ;
223 
224  float par1[nBinsEta][reco::GsfElectron::GAP+1] =
225  {
226  { 0.0111034, -0.00187886, -0.00187886, -0.00289304, -0.00187886 },
227  { -0.00969012, -0.00227574, -0.00227574, -0.00182187, -0.00227574 },
228  { 0.00389454, -0.00259935, -0.00259935, -0.00211059, -0.00259935 },
229  { 0.017095, -0.00433692, -0.00302335, -0.00241385, -0.00433692 },
230  { -0.00049009, -0.00551324, -0.00302335, -0.00532352, -0.00551324 },
231  { -0.00252723, -0.00799669, -0.00302335, -0.00823109, -0.00799669 },
232  { 0.00332567, -0.00870057, -0.0170581, -0.011482, -0.00870057 },
233  { -0.213285, -0.000771577, -0.036007, -0.0187526, -0.000771577 },
234  { 0.1741, -0.00202028, -0.0233995, -0.0302066, -0.00202028 },
235  { 0.152794, 0.00441308, -0.0468563, -0.0158817, 0.00441308 },
236  { 0.0351465, 0.00832913, 0.0358028, -0.0233262, 0.00832913 },
237  { 0.185781, 0.00742879, 0.08858, 0.00568078, 0.00742879 },
238  { 0.153088, 0.0094608, 0.0489979, 0.0491897, 0.0094608 },
239  { 0.296681, 0.00560406, 0.106492, 0.0652007, 0.00560406 }
240  } ;
241 
242  float par2[nBinsEta][reco::GsfElectron::GAP+1] =
243  {
244  { -0.00330844, 0, 0, 5.62441e-05, 0 },
245  { 0.00329373, 0, 0, -0.000113883, 0 },
246  { -0.00104661, 0, 0, -0.000152794, 0 },
247  { -0.0060409, 0, -0.000257724, -0.000202099, 0 },
248  { -0.000742866, 0, -0.000257724, -2.06003e-05, 0 },
249  { -0.00205425, 0, -0.000257724, 3.84179e-05, 0 },
250  { -0.00350757, 0, 0.000590483, 0.000323723, 0 },
251  { 0.0794596, -0.00276696, 0.00205854, 0.000356716, -0.00276696 },
252  { -0.092436, -0.00471028, 0.00062096, 0.00088347, -0.00471028 },
253  { -0.0855029, -0.00809139, 0.00284102, -0.00366903, -0.00809139 },
254  { -0.0306209, -0.00944584, -0.0145892, -0.00176969, -0.00944584 },
255  { -0.0996414, -0.00960462, -0.0328264, -0.00983844, -0.00960462 },
256  { -0.0784107, -0.010172, -0.0256722, -0.0215133, -0.010172 },
257  { -0.145815, -0.00943169, -0.0414525, -0.027087, -0.00943169 }
258  } ;
259 
260  float sigmaPhiSigmaEtaMin[reco::GsfElectron::GAP+1] = { 0.8, 0.8, 0.8, 0.8, 0.8 } ;
261  float sigmaPhiSigmaEtaMax[reco::GsfElectron::GAP+1] = { 5., 5., 5., 5., 5. } ;
262  float sigmaPhiSigmaEtaFit[reco::GsfElectron::GAP+1] = { 1.2, 1.2, 1.2, 1.2, 1.2 } ;
263 
264  // extra protections
265  // fix sigmaPhiSigmaEta boundaries
266  if ( sigmaPhiSigmaEta < sigmaPhiSigmaEtaMin[cl] )
267  { sigmaPhiSigmaEta = sigmaPhiSigmaEtaMin[cl] ; }
268  if ( sigmaPhiSigmaEta > sigmaPhiSigmaEtaMax[cl] )
269  { sigmaPhiSigmaEta = sigmaPhiSigmaEtaMax[cl] ; }
270 
271  // In eta cracks/gaps
272  if ( tmpEta == -1 ) // need to interpolate
273  {
274  float tmpInter = 1 ;
275  for ( int iEta = 0 ; iEta < nBinsEta-1 ; ++iEta )
276  {
277  if ( rightEta[iEta] <= TMath::Abs(eta) && TMath::Abs(eta) <leftEta[iEta+1] )
278  {
279  if ( sigmaPhiSigmaEta >= sigmaPhiSigmaEtaFit[cl] )
280  {
281  tmpInter =
282  ( par0[iEta][cl] +
283  sigmaPhiSigmaEta*par1[iEta][cl] +
284  sigmaPhiSigmaEta*sigmaPhiSigmaEta*par2[iEta][cl] +
285  par0[iEta+1][cl] +
286  sigmaPhiSigmaEta*par1[iEta+1][cl] +
287  sigmaPhiSigmaEta*sigmaPhiSigmaEta*par2[iEta+1][cl] ) / 2. ;
288  }
289  else tmpInter = (xcorr[iEta][cl] + xcorr[iEta+1][cl])/2. ;
290  }
291  }
292  return tmpInter ;
293  }
294 
295  if (sigmaPhiSigmaEta >= sigmaPhiSigmaEtaFit[cl])
296  { return par0[tmpEta][cl] + sigmaPhiSigmaEta*par1[tmpEta][cl] + sigmaPhiSigmaEta*sigmaPhiSigmaEta*par2[tmpEta][cl] ; }
297  else
298  { return xcorr[tmpEta][cl] ; }
299 
300  return 1. ;
301  }
302 
304  {
305  if (algorithm==0) //Electrons EB
306  {
307  const float parClassIndep[5] = { 0.97213, 0.999528, 5.61192e-06, 0.0143269, -17.1776 } ;
308  const float par[reco::GsfElectron::GAP+1][5] =
309  {
310  { 0.974327, 0.996127, 5.99401e-05, 0.159813, -3.80392 },
311  { 0.97213, 0.999528, 5.61192e-06, 0.0143269, -17.1776 },
312  { 0.940666, 0.988894, 0.00017474, 0.25603, -4.58153 },
313  { 0.969526, 0.98572, 0.000193842, 4.21548, -1.37159 },
314  { parClassIndep[0], parClassIndep[1], parClassIndep[2], parClassIndep[3], parClassIndep[4] }
315  } ;
316  if ( ET > 200 ) { ET = 200 ; }
317  if ( ET > 100 ) { return (parClassIndep[1]+ET*parClassIndep[2]) * (1-parClassIndep[3]*exp(ET/parClassIndep[4])) ; }
318  if ( ET >= 10 ) { return (par[cl][1]+ET*par[cl][2]) * (1-par[cl][3]*exp(ET/par[cl][4])) ; }
319  if ( ET >=5 ) { return par[cl][0] ; }
320  return 1. ;
321  }
322  else if (algorithm==1) //Electrons EE
323  {
324  float par[reco::GsfElectron::GAP+1][5] =
325  {
326  { 0.930081, 0.996683, 3.54079e-05, 0.0460187, -23.2461 },
327  { 0.930081, 0.996683, 3.54079e-05, 0.0460187, -23.2461 },
328  { 0.930081, 0.996683, 3.54079e-05, 0.0460187, -23.2461 },
329  { 0.930081, 0.996683, 3.54079e-05, 0.0460187, -23.2461 },
330  { 0.930081, 0.996683, 3.54079e-05, 0.0460187, -23.2461 }
331  };
332  if ( ET > 200 ) { ET = 200 ; }
333  if ( ET < 5 ) { return 1. ; }
334  if ( 5 <= ET && ET < 10 ) { return par[cl][0] ; }
335  if ( 10 <= ET && ET <= 200 ) { return ( par[cl][1] + ET*par[cl][2])*(1-par[cl][3]*exp(ET/par[cl][4])) ; }
336  }
337  else
338  { edm::LogWarning("ElectronEnergyCorrector::fEt")<<"algorithm should be 0 or 1 for electrons !" ; }
339  return 1. ;
340  }
341 
343  {
344  if (algorithm==0) // Electrons EB
345  { return 1. ; }
346  else if (algorithm==1) // Electrons EE
347  {
348  float par0[reco::GsfElectron::GAP+1] = { 400, 400, 400, 400, 400 } ;
349  float par1[reco::GsfElectron::GAP+1] = { 0.999545, 0.982475, 0.986217, 0.996763, 0.982475 } ;
350  float par2[reco::GsfElectron::GAP+1] = { 1.26568e-05, 4.95413e-05, 5.02161e-05, 2.8485e-06, 4.95413e-05 } ;
351  float par3[reco::GsfElectron::GAP+1] = { 0.0696757, 0.16886, 0.115317, 0.12098, 0.16886 } ;
352  float par4[reco::GsfElectron::GAP+1] = { -54.3468, -30.1517, -26.3805, -62.0538, -30.1517 } ;
353 
354  if ( E > par0[cl] ) { E = par0[cl] ; }
355  if ( E < 0 ) { return 1. ; }
356  if ( 0 <= E && E <= par0[cl] ) { return (par1[cl] + E*par2[cl] )*(1- par3[cl]*exp(E/par4[cl] )) ; }
357  }
358  else
359  { edm::LogWarning("ElectronEnergyCorrector::fEnergy")<<"algorithm should be 0 or 1 for electrons !" ; }
360  return 1.;
361  }
362 
363 
364 //==========================================================================
365 //
366 //==========================================================================
367 
368 double ElectronEnergyCorrector::fEtaBarrelGood( double scEta ) const
369  {
370  // f(eta) for the first 3 classes (0, 10 and 20) (estimated from 1Mevt single e sample)
371  // Ivica's new corrections 01/06
372  float p0 = 1.00149e+00 ;
373  float p1 = -2.06622e-03 ;
374  float p2 = -1.08793e-02 ;
375  float p3 = 1.54392e-02 ;
376  float p4 = -1.02056e-02 ;
377  double x = (double) std::abs(scEta) ;
378  return p0 + p1*x + p2*x*x + p3*x*x*x + p4*x*x*x*x ;
379  }
380 
381 double ElectronEnergyCorrector::fEtaBarrelBad(double scEta) const
382  {
383  // f(eta) for the class = 30 (estimated from 1Mevt single e sample)
384  // Ivica's new corrections 01/06
385  float p0 = 9.99063e-01;
386  float p1 = -2.63341e-02;
387  float p2 = 5.16054e-02;
388  float p3 = -4.95976e-02;
389  float p4 = 3.62304e-03;
390  double x = (double) std::abs(scEta) ;
391  return p0 + p1*x + p2*x*x + p3*x*x*x + p4*x*x*x*x ;
392  }
393 
394 double ElectronEnergyCorrector::fEtaEndcapGood( double scEta ) const
395  {
396  // f(eta) for the first 3 classes (100, 110 and 120)
397  // Ivica's new corrections 01/06
398  float p0 = -8.51093e-01 ;
399  float p1 = 3.54266e+00 ;
400  float p2 = -2.59288e+00 ;
401  float p3 = 8.58945e-01 ;
402  float p4 = -1.07844e-01 ;
403  double x = (double) std::abs(scEta) ;
404  return p0 + p1*x + p2*x*x + p3*x*x*x + p4*x*x*x*x ;
405  }
406 
407 double ElectronEnergyCorrector::fEtaEndcapBad( double scEta ) const
408  {
409  // f(eta) for the class = 130-134
410  // Ivica's new corrections 01/06
411  float p0 = -4.25221e+00 ;
412  float p1 = 1.01936e+01 ;
413  float p2 = -7.48247e+00 ;
414  float p3 = 2.45520e+00 ;
415  float p4 = -3.02872e-01 ;
416  double x = (double) std::abs(scEta);
417  return p0 + p1*x + p2*x*x + p3*x*x*x + p4*x*x*x*x ;
418  }
419 
double fEtaBarrelGood(double scEta) const
bool isEcalEnergyCorrected() const
Definition: GsfElectron.h:845
bool isEE() const
Definition: GsfElectron.h:357
bool isEB() const
Definition: GsfElectron.h:356
double fEtaEndcapBad(double scEta) const
float fBremEta(float sigmaPhiSigmaEta, float eta, int algorithm, reco::GsfElectron::Classification cl) const
return((rh^lh)&mask)
virtual float getValue(const reco::BasicCluster &, const EcalRecHitCollection &) const =0
void setCorrectedEcalEnergyError(float newEnergyError)
Definition: GsfElectron.cc:179
void simpleParameterizationUncertainty(reco::GsfElectron &)
T sqrt(T t)
Definition: SSEVec.h:18
double p4[4]
Definition: TauolaWrapper.h:92
T Abs(T a)
Definition: MathUtil.h:49
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
double p2[4]
Definition: TauolaWrapper.h:90
JetCorrectorParameters corr
Definition: classes.h:5
float correctedEcalEnergy() const
Definition: GsfElectron.h:846
float fEt(float et, int algorithm, reco::GsfElectron::Classification cl) const
void classBasedParameterizationEnergy(reco::GsfElectron &, const reco::BeamSpot &bs)
double fEtaEndcapGood(double scEta) const
Classification classification() const
Definition: GsfElectron.h:768
double computeElectronEnergyUncertainty(reco::GsfElectron::Classification c, double eta, double brem, double energy)
et
define resolution functions of each parameter
float fEnergy(float e, int algorithm, reco::GsfElectron::Classification cl) const
void setCorrectedEcalEnergy(float newEnergy)
Definition: GsfElectron.cc:182
double p1[4]
Definition: TauolaWrapper.h:89
double fEtaBarrelBad(double scEta) const
SuperClusterRef superCluster() const override
reference to a SuperCluster
Definition: GsfElectron.h:185
#define ET
float fEta(float energy, float eta, int algorithm) const
void classBasedParameterizationUncertainty(reco::GsfElectron &)
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
EcalClusterFunctionBaseClass * crackCorrectionFunction_
float energyError(float E, float *par)
double p3[4]
Definition: TauolaWrapper.h:91