CMS 3D CMS Logo

ElectronEnergyCalibrator.cc
Go to the documentation of this file.
2 
3 #include <CLHEP/Random/RandGaussQ.h>
4 #include <CLHEP/Random/RandFlat.h>
5 #include <CLHEP/Random/Random.h>
9 
10 /****************************************************************************
11  *
12  * Propagate SC calibration from Zee fit to the electrons
13  *
14  ****************************************************************************/
15 
16 #include <string>
17 #include <vector>
18 #include <fstream>
19 #include <sstream>
20 #include <iostream>
21 
22 
24 {
25  if ( !isMC_ ) // DATA
26  {
27  if ( verbose_ )
28  {
29  std::cout << "[ElectronEnergyCalibrator] Initialization in DATA mode" << std::endl;
30  }
31 
32  std::ifstream fin(pathData_.c_str());
33 
34  if (!fin){
35  throw cms::Exception("Configuration")
36  << "[ElectronEnergyCalibrator] Cannot open the file "
37  << pathData_ << "\n It is not found, missed or corrupted" ;
38  } else
39  {
40  if ( verbose_ )
41  {
42  std::cout << "[ElectronEnergyCalibrator] File "
43  << pathData_ << " succesfully opened" << std::endl;
44  }
45 
46  std::string s;
47  std::vector<std::string> selements;
48  std::string delimiter = ",";
49  nCorrValRaw = 0;
50 
51  while ( !fin.eof() )
52  {
53  getline(fin, s);
54  if ( !s.empty() )
55  {
56  splitString(s, selements, delimiter);
67 
68  nCorrValRaw++;
69 
70  selements.clear();
71  }
72  }
73 
74  fin.close();
75 
76  if ( verbose_ )
77  {
78  std::cout << "[ElectronEnergyCalibrator] File closed" << std::endl;
79  }
80 
81  }
82  // linearity corrections data
84  {
85  std::ifstream finlin(pathLinData_.c_str());
86 
87  if (!finlin)
88  {
89  throw cms::Exception("Configuration")
90  << "[ElectronEnergyCalibrator] Cannot open the file "<< pathLinData_ << "\n It is not found, missed or corrupted" ;
91  } else
92  {
93  if (verbose_)
94  {
95  std::cout<<"[ElectronEnergyCalibrator] File with Linearity Corrections "<<pathLinData_<<" succesfully opened"<<std::endl;
96  }
97 
98  std::string s;
99  std::vector<std::string> selements;
100  std::string delimiter = ",";
101  nLinCorrValRaw = 0;
102 
103  while ( !finlin.eof() )
104  {
105  getline(finlin, s);
106  if ( !s.empty() )
107  {
108  splitString(s, selements, delimiter);
109 
118 
119  nLinCorrValRaw++;
120 
121  selements.clear();
122  }
123  }
124 
125  finlin.close();
126 
127  if (verbose_)
128  {
129  std::cout<<"[ElectronEnergyCalibrator] File closed"<<std::endl;
130  }
131  }
132  }
133  } else // MC
134  {
135  if ( verbose_ )
136  {
137  std::cout << "[ElectronEnergyCalibrator] Initialization in MC mode" << std::endl;
138  }
139  }
140 }
141 
142 void ElectronEnergyCalibrator::splitString(const std::string &fullstr, std::vector<std::string> &elements, const std::string &delimiter)
143 {
144  std::string::size_type lastpos = fullstr.find_first_not_of(delimiter, 0);
145  std::string::size_type pos = fullstr.find_first_of(delimiter, lastpos);
146 
147  while ( ( std::string::npos != pos ) || ( std::string::npos != lastpos ) )
148  {
149  elements.push_back(fullstr.substr(lastpos, pos-lastpos));
150  lastpos = fullstr.find_first_not_of(delimiter, pos);
151  pos = fullstr.find_first_of(delimiter, lastpos);
152  }
153 }
154 
156 {
157  std::istringstream stm;
158  double val = 0;
159  stm.str(str);
160  stm >> val;
161  return val;
162 }
163 
165 {
166  double scale = 1.0;
167  double dsigMC=0.;
168  double corrMC=0.;
169  double run_ = electron.getRunNumber();
170  bool isEB = electron.isEB();
171  double eta = electron.getEta();
172  double r9 = electron.getR9();
173 
174  switch ( correctionsType_ )
175  {
176  case 1:
177  if ( verbose_ )
178  {
179  std::cout << "[ElectronEnergyCalibrator] Using regression energy for calibration" << std::endl;
180  }
181  newEnergy_ = electron.getRegEnergy();
182  newEnergyError_ = electron.getRegEnergyError();
183  break;
184  case 2:
185  if( verbose_ )
186  {
187  std::cout << "[ElectronEnergyCalibrator] Using scale corrections for new regression" << std::endl;
188  }
189  newEnergy_ = electron.getRegEnergy();
190  newEnergyError_ = electron.getRegEnergyError();
191  break;
192  case 3:
193  if ( verbose_ )
194  {
195  std::cout << "[ElectronEnergyCalibrator] Using standard ecal energy for calibration" << std::endl;
196  }
197  newEnergy_ = electron.getSCEnergy();
198  newEnergyError_ = electron.getSCEnergyError();
199  break;
200  }
201 
203  if ( !rng.isAvailable() )
204  {
205  throw cms::Exception("Configuration")
206  << "XXXXXXX requires the RandomNumberGeneratorService\n"
207  "which is not present in the configuration file. You must add the service\n"
208  "in the configuration file or remove the modules that require it.";
209  }
210 
211  if (!isMC_ )
212  {
213  for ( int i=0; i < nCorrValRaw; i++ )
214  {
215  if ( (run_ >= corrValArray[i].nRunMin ) && ( run_ <= corrValArray[i].nRunMax ) )
216  {
217  if ( isEB )
218  {
219  if ( fabs(eta) < 1 )
220  {
221  if ( r9<0.94 )
222  {
223  scale = corrValArray[i].corrCat0;
224  } else
225  {
226  scale = corrValArray[i].corrCat1;
227  }
228  } else
229  {
230  if ( r9<0.94 )
231  {
232  scale = corrValArray[i].corrCat2;
233  } else
234  {
235  scale = corrValArray[i].corrCat3;
236  }
237  }
238  } else
239  {
240  if ( fabs(eta) < 2 )
241  {
242  if ( r9<0.94 )
243  {
244  scale = corrValArray[i].corrCat4;
245  } else
246  {
247  scale = corrValArray[i].corrCat5;
248  }
249  } else
250  {
251  if ( r9<0.94 )
252  {
253  scale = corrValArray[i].corrCat6;
254  } else
255  {
256  scale = corrValArray[i].corrCat7;
257  }
258  }
259  }
260  }
261  }
263  }
264 
265  switch ( correctionsType_ )
266  {
267  case 1:
268  // Implementation of the MC smearing for regression energy type 1
269  if ( dataset_ == "Summer12_DR53X_HCP2012" || dataset_ == "Moriond2013" )
270  {
271  if ( !isMC_ )
272  {
273  if ( run_ <= 203002 )
274  {
275  if ( isEB && fabs(eta) < 1 && r9 < 0.94 ) dsigMC = 0.0103;
276  if ( isEB && fabs(eta) < 1 && r9 >= 0.94 ) dsigMC = 0.0090;
277  if ( isEB && fabs(eta) >= 1 && r9 < 0.94 ) dsigMC = 0.0190;
278  if ( isEB && fabs(eta) >= 1 && r9 >= 0.94 ) dsigMC = 0.0156;
279  if ( !isEB && fabs(eta) < 2 && r9 < 0.94 ) dsigMC = 0.0269;
280  if ( !isEB && fabs(eta) < 2 && r9 >= 0.94 ) dsigMC = 0.0287;
281  if ( !isEB && fabs(eta) >= 2 && r9 < 0.94 ) dsigMC = 0.0364;
282  if ( !isEB && fabs(eta) >= 2 && r9 >= 0.94 ) dsigMC = 0.0321;
283  } else
284  {
285  if ( isEB && fabs(eta) < 1 && r9 < 0.94 ) dsigMC = 0.0109;
286  if ( isEB && fabs(eta) < 1 && r9 >= 0.94 ) dsigMC = 0.0099;
287  if ( isEB && fabs(eta) >= 1 && r9 < 0.94 ) dsigMC = 0.0182;
288  if ( isEB && fabs(eta) >= 1 && r9 >= 0.94 ) dsigMC = 0.0200;
289  if ( !isEB && fabs(eta) < 2 && r9 < 0.94 ) dsigMC = 0.0282;
290  if ( !isEB && fabs(eta) < 2 && r9 >= 0.94 ) dsigMC = 0.0309;
291  if ( !isEB && fabs(eta) >= 2 && r9 < 0.94 ) dsigMC = 0.0386;
292  if ( !isEB && fabs(eta) >= 2 && r9 >= 0.94 ) dsigMC = 0.0359;
293  }
294  } else
295  {
296  CLHEP::RandFlat flatRandom(rng->getEngine(streamID));
297  double rn = flatRandom.fire();
298  if ( rn > lumiRatio_ )
299  {
300  if ( isEB && fabs(eta) < 1 && r9 < 0.94 ) dsigMC = 0.0109;
301  if ( isEB && fabs(eta) < 1 && r9 >= 0.94 ) dsigMC = 0.0099;
302  if ( isEB && fabs(eta) >= 1 && r9 < 0.94 ) dsigMC = 0.0182;
303  if ( isEB && fabs(eta) >= 1 && r9 >= 0.94 ) dsigMC = 0.0200;
304  if ( !isEB && fabs(eta) < 2 && r9 < 0.94 ) dsigMC = 0.0282;
305  if ( !isEB && fabs(eta) < 2 && r9 >= 0.94 ) dsigMC = 0.0309;
306  if ( !isEB && fabs(eta) >= 2 && r9 < 0.94 ) dsigMC = 0.0386;
307  if ( !isEB && fabs(eta) >= 2 && r9 >= 0.94 ) dsigMC = 0.0359;
308  } else
309  {
310  if ( isEB && fabs(eta) < 1 && r9 < 0.94 ) dsigMC = 0.0103;
311  if ( isEB && fabs(eta) < 1 && r9 >= 0.94 ) dsigMC = 0.0090;
312  if ( isEB && fabs(eta) >= 1 && r9 < 0.94 ) dsigMC = 0.0190;
313  if ( isEB && fabs(eta) >= 1 && r9 >= 0.94 ) dsigMC = 0.0156;
314  if ( !isEB && fabs(eta) < 2 && r9 < 0.94 ) dsigMC = 0.0269;
315  if ( !isEB && fabs(eta) < 2 && r9 >= 0.94 ) dsigMC = 0.0287;
316  if ( !isEB && fabs(eta) >= 2 && r9 < 0.94 ) dsigMC = 0.0364;
317  if ( !isEB && fabs(eta) >= 2 && r9 >= 0.94 ) dsigMC = 0.0321;
318  }
319  if ( lumiRatio_ == 0.0 )
320  {
321  if ( isEB && fabs(eta) < 1 && r9 < 0.94 ) dsigMC = 0.0103;
322  if ( isEB && fabs(eta) < 1 && r9 >= 0.94 ) dsigMC = 0.0090;
323  if ( isEB && fabs(eta) >= 1 && r9 < 0.94 ) dsigMC = 0.0190;
324  if ( isEB && fabs(eta) >= 1 && r9 >= 0.94 ) dsigMC = 0.0156;
325  if ( !isEB && fabs(eta) < 2 && r9 < 0.94 ) dsigMC = 0.0269;
326  if ( !isEB && fabs(eta) < 2 && r9 >= 0.94 ) dsigMC = 0.0287;
327  if ( !isEB && fabs(eta) >= 2 && r9 < 0.94 ) dsigMC = 0.0364;
328  if ( !isEB && fabs(eta) >= 2 && r9 >= 0.94 ) dsigMC = 0.0321;
329  }
330  if ( lumiRatio_ == 1.0 )
331  {
332  if ( isEB && fabs(eta) < 1 && r9 < 0.94 ) dsigMC = 0.0109;
333  if ( isEB && fabs(eta) < 1 && r9 >= 0.94 ) dsigMC = 0.0099;
334  if ( isEB && fabs(eta) >= 1 && r9 < 0.94 ) dsigMC = 0.0182;
335  if ( isEB && fabs(eta) >= 1 && r9 >= 0.94 ) dsigMC = 0.0200;
336  if ( !isEB && fabs(eta) < 2 && r9 < 0.94 ) dsigMC = 0.0282;
337  if ( !isEB && fabs(eta) < 2 && r9 >= 0.94 ) dsigMC = 0.0309;
338  if ( !isEB && fabs(eta) >= 2 && r9 < 0.94 ) dsigMC = 0.0386;
339  if ( !isEB && fabs(eta) >= 2 && r9 >= 0.94 ) dsigMC = 0.0359;
340  }
341  }
342  }
343  break;
344 
345  case 2:
346  // Implementation of the MC smearing for regression energy type 2
347  if ( dataset_ == "Summer12_LegacyPaper" || dataset_ == "22Jan2013ReReco" )
348  {
349  if ( isEB && fabs(eta) < 1 && r9 < 0.94 ) dsigMC = 0.0094;
350  if ( isEB && fabs(eta) < 1 && r9 >= 0.94 ) dsigMC = 0.0092;
351  if ( isEB && fabs(eta) >= 1 && r9 < 0.94 ) dsigMC = 0.0182;
352  if ( isEB && fabs(eta) >= 1 && r9 >= 0.94 ) dsigMC = 0.0139;
353  if ( !isEB && fabs(eta) < 2 && r9 < 0.94 ) dsigMC = 0.0220;
354  if ( !isEB && fabs(eta) < 2 && r9 >= 0.94 ) dsigMC = 0.0229;
355  if ( !isEB && fabs(eta) >= 2 && r9 < 0.94 ) dsigMC = 0.0290;
356  if ( !isEB && fabs(eta) >= 2 && r9 >= 0.94 ) dsigMC = 0.0234;
357  }
358  break;
359 
360  case 3:
361  // standard SC energy scale corrections implementation
362  if ( dataset_ == "Summer11" || dataset_ == "ReReco" )
363  { // values from https://indico.cern.ch/conferenceDisplay.py?confId=146386
364  if ( isEB && fabs(eta) < 1 && r9 < 0.94 ) dsigMC = 0.01;
365  if ( isEB && fabs(eta) < 1 && r9 >= 0.94 ) dsigMC = 0.0099;
366  if ( isEB && fabs(eta) >= 1 && r9 < 0.94 ) dsigMC = 0.0217;
367  if ( isEB && fabs(eta) >= 1 && r9 >= 0.94 ) dsigMC = 0.0157;
368  if ( !isEB && fabs(eta) < 2 && r9 < 0.94 ) dsigMC = 0.0326;
369  if ( !isEB && fabs(eta) < 2 && r9 >= 0.94 ) dsigMC = 0.0330;
370  if ( !isEB && fabs(eta) >= 2 && r9 < 0.94 ) dsigMC = 0.0331;
371  if ( !isEB && fabs(eta) >= 2 && r9 >= 0.94 ) dsigMC = 0.0378;
372  } else if ( dataset_ == "Fall11" || dataset_ == "Jan16ReReco" )
373  { // values from https://hypernews.cern.ch/HyperNews/CMS/get/higgs2g/634.html, consistant with Jan16ReReco corrections
374  if ( isEB && fabs(eta) < 1 && r9 < 0.94 ) dsigMC = 0.0096;
375  if ( isEB && fabs(eta) < 1 && r9 >= 0.94 ) dsigMC = 0.0074;
376  if ( isEB && fabs(eta) >= 1 && r9 < 0.94 ) dsigMC = 0.0196;
377  if ( isEB && fabs(eta) >= 1 && r9 >= 0.94 ) dsigMC = 0.0141;
378  if ( !isEB && fabs(eta) < 2 && r9 < 0.94 ) dsigMC = 0.0279;
379  if ( !isEB && fabs(eta) < 2 && r9 >= 0.94 ) dsigMC = 0.0268;
380  if ( !isEB && fabs(eta) >= 2 && r9 < 0.94 ) dsigMC = 0.0301;
381  if ( !isEB && fabs(eta) >= 2 && r9 >= 0.94 ) dsigMC = 0.0293;
382  } else if ( dataset_ == "Summer12" || dataset_ == "ICHEP2012" )
383  { // new values from https://twiki.cern.ch/twiki/pub/CMS/EcalEnergyResolutionWithZee/oriented-ICHEP-scales_resolution.pdf
384  if ( isEB && fabs(eta) < 1 && r9 < 0.94 ) dsigMC = 0.0119;
385  if ( isEB && fabs(eta) < 1 && r9 >= 0.94 ) dsigMC = 0.0107;
386  if ( isEB && fabs(eta) >= 1 && r9 < 0.94 ) dsigMC = 0.0240;
387  if ( isEB && fabs(eta) >= 1 && r9 >= 0.94 ) dsigMC = 0.0149;
388  if ( !isEB && fabs(eta) < 2 && r9 < 0.94 ) dsigMC = 0.0330;
389  if ( !isEB && fabs(eta) < 2 && r9 >= 0.94 ) dsigMC = 0.0375;
390  if ( !isEB && fabs(eta) >= 2 && r9 < 0.94 ) dsigMC = 0.0602;
391  if ( !isEB && fabs(eta) >= 2 && r9 >= 0.94 ) dsigMC = 0.0607;
392  } else if ( dataset_ == "Summer12_DR53X_HCP2012" || dataset_ == "Moriond2013" )
393  {
394  if ( isEB && fabs(eta) < 1 && r9 < 0.94 ) dsigMC = 0.0099;
395  if ( isEB && fabs(eta) < 1 && r9 >= 0.94 ) dsigMC = 0.0103;
396  if ( isEB && fabs(eta) >= 1 && r9 < 0.94 ) dsigMC = 0.0219;
397  if ( isEB && fabs(eta) >= 1 && r9 >= 0.94 ) dsigMC = 0.0158;
398  if ( !isEB && fabs(eta) < 2 && r9 < 0.94 ) dsigMC = 0.0222;
399  if ( !isEB && fabs(eta) < 2 && r9 >= 0.94 ) dsigMC = 0.0298;
400  if ( !isEB && fabs(eta) >= 2 && r9 < 0.94 ) dsigMC = 0.0318;
401  if ( !isEB && fabs(eta) >= 2 && r9 >= 0.94 ) dsigMC = 0.0302;
402  }
403  break;
404  }
405 
406  if ( isMC_ )
407  {
408  CLHEP::RandGaussQ gaussDistribution(rng->getEngine(streamID), 1.,dsigMC);
409  corrMC = gaussDistribution.fire();
410  if ( verbose_ )
411  {
412  std::cout << "[ElectronEnergyCalibrator] unsmeared energy " << newEnergy_ << std::endl;
413  }
414  if ( synchronization_ )
415  {
416  std::cout << "[ElectronEnergyCalibrator] "
417  << "======================= SYNCRONIZATION MODE! ======================="
418  << std::endl;
419  newEnergy_ = newEnergy_*(1+dsigMC);
420  } else
421  {
422  newEnergy_ = newEnergy_*corrMC;
423  }
424  if ( verbose_ )
425  {
426  std::cout << "[ElectronEnergyCalibrator] smeared energy " << newEnergy_ << std::endl;
427  }
428  }
429 
430  // correct energy error for MC and for data as error is obtained from (ideal) MC parametrisation
431  if ( updateEnergyErrors_ )
432  {
434  }
435  if ( verbose_ )
436  {
437  std::cout << "[ElectronEnergyCalibrator] initial energy "
438  << electron.getRegEnergy() << " recalibrated energy " << newEnergy_ << std::endl;
439  }
440  if ( verbose_ )
441  {
442  std::cout << "[ElectronEnergyCalibrator] initial energy error "
443  << electron.getRegEnergyError() << " recalibrated energy error "
444  << newEnergyError_ << std::endl;
445  }
446 
447  electron.setNewEnergy(newEnergy_);
449 }
450 
452 {
454  {
455  bool isEB = electron.isEB();
456  double eta = electron.getEta();
457  double theta = 2*atan(exp(-eta));
458  double p = electron.getCombinedMomentum();
459  double pt = p * fabs(sin(theta));
460  int classification = electron.getElClass();
461  double linscale = 0.;
462 
463  for (int i=0; i < nLinCorrValRaw; i++)
464  {
465  if ((pt >= linCorrValArray[i].ptMin) && (pt <= linCorrValArray[i].ptMax))
466  {
467  if (isEB)
468  {
469  if (fabs(eta) < 1)
470  {
471  if (classification<2)
472  {
473  linscale = linCorrValArray[i].corrCat0;
474  } else
475  {
476  linscale = linCorrValArray[i].corrCat3;
477  }
478  } else
479  {
480  if (classification<2)
481  {
482  linscale = linCorrValArray[i].corrCat1;
483  } else
484  {
485  linscale = linCorrValArray[i].corrCat4;
486  }
487  }
488  } else // !isEB
489  {
490  if (classification<2)
491  {
492  linscale = linCorrValArray[i].corrCat2;
493  } else
494  {
495  linscale = linCorrValArray[i].corrCat5;
496  }
497  }
498  }
499  }
500  double newP = p/(1.+linscale);
501  if (verbose_)
502  {
503  std::cout << "[ElectronEnergyCalibrator] Applying a linearity correction of " << 1./(1.+linscale) << " to " << pt << " GeV in pt" << std::endl;
504  }
505  electron.setCombinedMomentum(newP);
506  if (verbose_)
507  {
508  std::cout << "[ElectronEnergyCalibrator] calibrated transverse momentum " << pt << " GeV recalibrated for linearity to momentum " << electron.getCombinedMomentum()*fabs(sin(theta)) << " GeV" << std::endl;
509  }
510  }
511 }
correctionValues corrValArray[100]
int getElClass() const
double getSCEnergy() const
void setNewEnergyError(double newEnergyError)
double getRegEnergy() const
void setCombinedMomentum(double combinedMomentum)
void setNewEnergy(double newEnergy)
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
Geom::Theta< T > theta() const
void calibrate(SimpleElectron &electron, edm::StreamID const &)
virtual CLHEP::HepRandomEngine & getEngine(StreamID const &)=0
Use this engine in event methods.
uint16_t size_type
double stringToDouble(const std::string &str)
double getCombinedMomentum() const
unsigned int getRunNumber() const
T sqrt(T t)
Definition: SSEVec.h:18
bool isAvailable() const
Definition: Service.h:46
double getEta() const
float getR9() const
void correctLinearity(SimpleElectron &electron)
linearityCorrectionValues linCorrValArray[100]
double getRegEnergyError() const
double getSCEnergyError() const
bool isEB() const
void splitString(const std::string &fullstr, std::vector< std::string > &elements, const std::string &delimiter)