CMS 3D CMS Logo

EnergyScaleCorrection_class.cc
Go to the documentation of this file.
3 #include <RooDataSet.h>
4 #include <RooArgSet.h>
5 
6 #define NGENMAX 10000
7 
8 // for exit(0)
9 #include <cstdlib>
10 #include <cfloat>
11 // for setw
12 #include <iomanip>
13 //for istreamstring
14 #include <sstream>
15 
16 //#define DEBUG
17 //#define PEDANTIC_OUTPUT
18 
20  doScale(false), doSmearings(false),
21  smearingType_(ECALELF)
22 {
23 
24  if(!correctionFileName.empty()) {
25  std::string filename = correctionFileName+"_scales.dat";
26  ReadFromFile(filename);
27  if(scales.empty()) {
28  std::cerr << "[ERROR] scale correction map empty" << std::endl;
29  exit(1);
30  }
31  }
32 
33  if(!correctionFileName.empty()) {
34  std::string filename = correctionFileName+"_smearings.dat";
35  ReadSmearingFromFile(filename);
36  if(smearings.empty()) {
37  std::cerr << "[ERROR] smearing correction map empty" << std::endl;
38  exit(1);
39  }
40  }
41 
42  return;
43 }
44 
46 {
47  return;
48 }
49 
50 
51 
53  double R9Ele, double etaSCEle, double EtEle) const
54 {
55  float correction = 1;
56  if(doScale) correction *= getScaleOffset(runNumber, isEBEle, R9Ele, etaSCEle, EtEle);
57 
58  return correction;
59 }
60 
61 
62 
64  double R9Ele, double etaSCEle, double EtEle) const
65 {
66  float statUncert = getScaleStatUncertainty(runNumber, isEBEle, R9Ele, etaSCEle, EtEle);
67  float systUncert = getScaleSystUncertainty(runNumber, isEBEle, R9Ele, etaSCEle, EtEle);
68 
69  return sqrt(statUncert * statUncert + systUncert * systUncert);
70 }
71 
72 
74  double R9Ele, double etaSCEle, double EtEle) const
75 {
76  return getScaleCorrection(runNumber, isEBEle, R9Ele, etaSCEle, EtEle).scale_err;
77 }
78 
80  double R9Ele, double etaSCEle, double EtEle) const
81 {
82  return getScaleCorrection(runNumber, isEBEle, R9Ele, etaSCEle, EtEle).scale_err_syst;
83 }
84 
85 
86 correctionValue_class EnergyScaleCorrection_class::getScaleCorrection(unsigned int runNumber, bool isEBEle, double R9Ele, double etaSCEle, double EtEle) const
87 {
88  // buld the category based on the values of the object
89  correctionCategory_class category(runNumber, etaSCEle, R9Ele, EtEle);
90  correction_map_t::const_iterator corr_itr = scales.find(category); // find the correction value in the map that associates the category to the correction
91 
92  if(corr_itr == scales.end()) { // if not in the standard classes, add it in the list of not defined classes
93 
94  // this part is commented because it makes the method not constant
95  // if(scales_not_defined.count(category) == 0) {
96  // correctionValue_class corr;
97  // scales_not_defined[category] = corr;
98  // }
99  // corr_itr = scales_not_defined.find(category);
101  std::cout << "[ERROR] Scale category not found: " << std::endl;
102  std::cout << category << std::endl;
103  std::cout << "Returning uncorrected value." << std::endl;
104  // exit(1);
105  correctionValue_class nocorr;
106  std::cout << nocorr << std::endl;
107  return nocorr;
108  }
109 
110 #ifdef DEBUG
111  std::cout << "[DEBUG] Checking scale correction for category: " << category << std::endl;
112  std::cout << "[DEBUG] Correction is: " << corr_itr->second
113  << std::endl
114  << " given for category " << corr_itr->first << std::endl;;
115 #endif
116  return corr_itr->second;
117 }
118 
119 float EnergyScaleCorrection_class::getScaleOffset(unsigned int runNumber, bool isEBEle, double R9Ele, double etaSCEle, double EtEle) const
120 {
121  // buld the category based on the values of the object
122  correctionCategory_class category(runNumber, etaSCEle, R9Ele, EtEle);
123  correction_map_t::const_iterator corr_itr = scales.find(category); // find the correction value in the map that associates the category to the correction
124 
125  if(corr_itr == scales.end()) { // if not in the standard classes, add it in the list of not defined classes
126 
127  // this part is commented because it makes the method not constant
128  // if(scales_not_defined.count(category) == 0) {
129  // correctionValue_class corr;
130  // scales_not_defined[category] = corr;
131  // }
132  // corr_itr = scales_not_defined.find(category);
134  std::cout << "[ERROR] Scale offset category not found: " << std::endl;
135  std::cout << category << std::endl;
136  std::cout << "Returning uncorrected value." << std::endl;
137  // exit(1);
138  correctionValue_class nocorr;
139  return nocorr.scale;
140  }
141 
142 #ifdef DEBUG
143  std::cout << "[DEBUG] Checking scale offset correction for category: " << category << std::endl;
144  std::cout << "[DEBUG] Correction is: " << corr_itr->second
145  << std::endl
146  << " given for category " << corr_itr->first << std::endl;;
147 #endif
148 
149  return corr_itr->second.scale;
150 }
151 
152 
159 {
160 #ifdef PEDANTIC_OUTPUT
161  std::cout << "[STATUS] Reading energy scale correction values from file: " << filename << std::endl;
162 #endif
163  std::cout << "[STATUS] Reading energy scale correction values from file: " << filename << std::endl;
164 
165  //std::ifstream Ccufile(edm::FileInPath(Ccufilename).fullPath().c_str(),std::ios::in);
166  std::ifstream f_in(edm::FileInPath(filename).fullPath().c_str());
167 
168  if(!f_in.good()) {
169  std::cerr << "[ERROR] file " << filename << " not readable" << std::endl;
170  exit(1);
171  return;
172  }
173 
174  int runMin, runMax;
175  TString category, region2;
176  double deltaP, err_deltaP, err_deltaP_stat, err_deltaP_syst;
177 
178 
179  for(f_in >> category; f_in.good(); f_in >> category) {
180  f_in >> region2
181  >> runMin >> runMax
182  >> deltaP >> err_deltaP >> err_deltaP_stat >> err_deltaP_syst;
183 
184  AddScale(category, runMin, runMax, deltaP, err_deltaP_stat, err_deltaP_syst);
185  }
186 
187  f_in.close();
188 
189  return;
190 }
191 
192 // this method adds the correction values read from the txt file to the map
193 void EnergyScaleCorrection_class::AddScale(TString category_, int runMin_, int runMax_, double deltaP_, double err_deltaP_, double err_syst_deltaP)
194 {
195 
196  correctionCategory_class cat(category_); // build the category from the string
197  cat.runmin = runMin_;
198  cat.runmax = runMax_;
199 
200  // the following check is not needed, can be removed
201  if(scales.count(cat) != 0) {
202  std::cerr << "[ERROR] Category already defined!" << std::endl;
203  std::cerr << " Adding category: " << cat << std::endl;
204  std::cerr << " Defined category: " << scales[cat] << std::endl;
205  exit(1);
206  }
207 
208  correctionValue_class corr; // define the correction values
209  corr.scale = deltaP_;
210  corr.scale_err = err_deltaP_;
211  corr.scale_err_syst = err_syst_deltaP;
212  scales[cat] = corr;
213 
214 #ifdef PEDANTIC_OUTPUT
215  std::cout << "[INFO:scale correction] " << cat << corr << std::endl;
216 #endif
217  return;
218 }
219 
220 //============================== SMEARING
221 void EnergyScaleCorrection_class::AddSmearing(TString category_, int runMin_, int runMax_,
222  double rho, double err_rho, double phi, double err_phi,
223  double Emean, double err_Emean)
224 {
225 
226  correctionCategory_class cat(category_);
227  cat.runmin = (runMin_ < 0) ? 0 : runMin_;
228  cat.runmax = runMax_;
229 
230  if(smearings.count(cat) != 0) {
231  std::cerr << "[ERROR] Smearing category already defined!" << std::endl;
232  std::cerr << " Adding category: " << cat << std::endl;
233  std::cerr << " Defined category: " << smearings[cat] << std::endl;
234  exit(1);
235  }
236 
238  corr.rho = rho;
239  corr.rho_err = err_rho;
240  corr.phi = phi;
241  corr.phi_err = err_phi;
242  corr.Emean = Emean;
243  corr.Emean_err = err_Emean;
244  smearings[cat] = corr;
245 
246 #ifdef PEDANTIC_OUTPUT
247 #ifndef CMSSW
248  std::cout << "[INFO:smearings] " << cat << corr << std::endl;
249 #else
250  edm::LogInfo("[INFO:smearings] ") << cat << corr;
251 #endif
252 #endif
253 
254  return;
255 }
256 
273 {
274 #ifdef PEDANTIC_OUTPUT
275  std::cout << "[STATUS] Reading smearing values from file: " << filename << std::endl;
276 #endif
277  //edm::FileInPath(Ccufilename).fullPath().c_str(),std::ios::in); .fullPath().c_str()
278  std::ifstream f_in(edm::FileInPath(filename).fullPath().c_str());
279  if(!f_in.good()) {
280  std::cerr << "[ERROR] file " << filename << " not readable" << std::endl;
281  exit(1);
282  }
283 
284  int runMin = 0, runMax = 900000;
285  int unused = 0;
286  TString category, region2;
287  //double smearing, err_smearing;
288  double rho, phi, Emean, err_rho, err_phi, err_Emean;
289  double etaMin, etaMax, r9Min, r9Max;
290  std::string phi_string, err_phi_string;
291 
292 
293  while(f_in.peek() != EOF && f_in.good()) {
294  if(f_in.peek() == 10) { // 10 = \n
295  f_in.get();
296  continue;
297  }
298 
299  if(f_in.peek() == 35) { // 35 = #
300  f_in.ignore(1000, 10); // ignore the rest of the line until \n
301  continue;
302  }
303 
304  if(smearingType_ == UNKNOWN) { // trying to guess: not recommended
305  std::cerr << "[ERROR] Not implemented" << std::endl;
306  assert(false);
307 
308  } else if(smearingType_ == GLOBE) {
309  f_in >> category >> unused >> etaMin >> etaMax >> r9Min >> r9Max >> runMin >> runMax >>
310  Emean >> err_Emean >>
311  rho >> err_rho >> phi >> err_phi;
312 
313  AddSmearing(category, runMin, runMax, rho, err_rho, phi, err_phi, Emean, err_Emean);
314 
315  } else if(smearingType_ == ECALELF) {
316  f_in >> category >>
317  Emean >> err_Emean >>
318  rho >> err_rho >> phi_string >> err_phi_string;
319 #ifdef DEBUG
320  std::cout << category
321  //<< "\t" << etaMin << "\t" << etaMax << "\t" << r9Min << "\t" << r9Max << "\t" << runMin << "\t" << runMax
322  << "\tEmean=" << Emean << "\t"
323  << rho << "\t" << err_rho << "\tphi_string="
324  << phi_string << "#\terr_phi_string=" << err_phi_string << std::endl;
325 #endif
326 
327  if(phi_string=="M_PI_2") phi=M_PI_2;
328  else phi = std::stod(phi_string);
329 
330  if(err_phi_string=="M_PI_2") err_phi=M_PI_2;
331  else err_phi = std::stod(err_phi_string);
332 
333 
334  AddSmearing(category, runMin, runMax, rho, err_rho, phi, err_phi, Emean, err_Emean);
335 
336  } else {
337  f_in >> category >> rho >> phi;
338  err_rho = err_phi = Emean = err_Emean = 0;
339  AddSmearing(category, runMin, runMax, rho, err_rho, phi, err_phi, Emean, err_Emean);
340  }
341 #ifdef DEBUG
342  std::cout << category << "\t" << etaMin << "\t" << etaMax << "\t" << r9Min << "\t" << r9Max << "\t" << runMin << "\t" << runMax << "\tEmean=" << Emean << "\t" << rho << "\t" << phi << std::endl;
343 #endif
344  }
345 
346  f_in.close();
347  // runCorrection_itr=runMin_map.begin();
348 
349 
350  return;
351 }
352 
353 
354 
355 float EnergyScaleCorrection_class::getSmearingSigma(int runNumber, bool isEBEle, float R9Ele, float etaSCEle, float EtEle, paramSmear_t par, float nSigma) const
356 {
357  if (par == kRho) return getSmearingSigma(runNumber, isEBEle, R9Ele, etaSCEle, EtEle, nSigma, 0.);
358  if (par == kPhi) return getSmearingSigma(runNumber, isEBEle, R9Ele, etaSCEle, EtEle, 0., nSigma);
359  return getSmearingSigma(runNumber, isEBEle, R9Ele, etaSCEle, EtEle, 0., 0.);
360 }
361 
362 float EnergyScaleCorrection_class::getSmearingSigma(int runNumber, bool isEBEle, float R9Ele, float etaSCEle, float EtEle, float nSigma_rho, float nSigma_phi) const
363 {
364 
365  correctionCategory_class category(runNumber, etaSCEle, R9Ele, EtEle);
366  correction_map_t::const_iterator corr_itr = smearings.find(category);
367  if(corr_itr == smearings.end()) { // if not in the standard classes, add it in the list of not defined classes
368  // the following commented part makes the method non const
369  // if(smearings_not_defined.count(category) == 0) {
370  // correctionValue_class corr;
371  // smearings_not_defined[category] = corr;
372  // }
373  corr_itr = smearings_not_defined.find(category);
374  std::cerr << "[WARNING] Smearing category not found: " << std::endl;
375  std::cerr << category << std::endl;
376  // exit(1);
377  }
378 
379 #ifdef DEBUG
380  std::cout << "[DEBUG] Checking smearing correction for category: " << category << std::endl;
381  std::cout << "[DEBUG] Correction is: " << corr_itr->second
382  << std::endl
383  << " given for category " << corr_itr->first;
384 #endif
385 
386  double rho = corr_itr->second.rho + corr_itr->second.rho_err * nSigma_rho;
387  double phi = corr_itr->second.phi + corr_itr->second.phi_err * nSigma_phi;
388 
389  double constTerm = rho * sin(phi);
390  double alpha = rho * corr_itr->second.Emean * cos( phi);
391 
392  return sqrt(constTerm * constTerm + alpha * alpha / EtEle);
393 
394 }
395 
396 float EnergyScaleCorrection_class::getSmearingRho(int runNumber, bool isEBEle, float R9Ele, float etaSCEle, float EtEle) const
397 {
398 
399  correctionCategory_class category(runNumber, etaSCEle, R9Ele, EtEle);
400  correction_map_t::const_iterator corr_itr = smearings.find(category);
401  if(corr_itr == smearings.end()) { // if not in the standard classes, add it in the list of not defined classes
402  // if(smearings_not_defined.count(category) == 0) {
403  // correctionValue_class corr;
404  // smearings_not_defined[category] = corr;
405  // }
406  corr_itr = smearings_not_defined.find(category);
407  }
408 
409  return corr_itr->second.rho;
410 }
411 
413 {
414  if(runmin < b.runmin && runmax < b.runmax) return true;
415  if(runmax > b.runmax && runmin > b.runmin) return false;
416 
417  if(etamin < b.etamin && etamax < b.etamax) return true;
418  if(etamax > b.etamax && etamin > b.etamin) return false;
419 
420  if(r9min < b.r9min && r9max < b.r9max) return true;
421  if(r9max > b.r9max && r9min > b.r9min) return false;
422 
423  if(etmin < b.etmin && etmax < b.etmax) return true;
424  if(etmax > b.etmax && etmin > b.etmin) return false;
425  return false;
426 
427 }
428 
430 {
431  std::string category(category_.Data());
432 #ifdef DEBUG
433  std::cout << "[DEBUG] correctionClass defined for category: " << category << std::endl;
434 #endif
435  // default values (corresponding to an existing category -- the worst one)
436  runmin = 0;
437  runmax = 999999;
438  etamin = 2;
439  etamax = 7;
440  r9min = -1;
441  r9max = 0.94;
442  etmin = -1;
443  etmax = 99999.;
444 
445  size_t p1, p2;// boundary
446  // eta region
447  p1 = category.find("absEta_");
448  if(category.find("absEta_0_1") != std::string::npos) {
449  etamin = 0;
450  etamax = 1;
451  } else if(category.find("absEta_1_1.4442") != std::string::npos) {
452  etamin = 1;
453  etamax = 1.479;
454  } //etamax=1.4442; }
455  else if(category.find("absEta_1.566_2") != std::string::npos) {
456  etamin = 1.479;
457  etamax = 2;
458  } //etamin=1.566; etamax=2;}
459  else if(category.find("absEta_2_2.5") != std::string::npos) {
460  etamin = 2;
461  etamax = 3;
462  } else {
463  if(p1 != std::string::npos) {
464  p1 = category.find("_", p1);
465  p2 = category.find("_", p1 + 1);
466  etamin = TString(category.substr(p1 + 1, p2 - p1 - 1)).Atof();
467  p1 = p2;
468  p2 = category.find("-", p1);
469  etamax = TString(category.substr(p1 + 1, p2 - p1 - 1)).Atof();
470  }
471  }
472 
473  if(category.find("EBlowEta") != std::string::npos) {
474  etamin = 0;
475  etamax = 1;
476  };
477  if(category.find("EBhighEta") != std::string::npos) {
478  etamin = 1;
479  etamax = 1.479;
480  };
481  if(category.find("EElowEta") != std::string::npos) {
482  etamin = 1.479;
483  etamax = 2;
484  };
485  if(category.find("EEhighEta") != std::string::npos) {
486  etamin = 2;
487  etamax = 7;
488  };
489 
490  // Et region
491  p1 = category.find("-Et_");
492  // p2 = p1 + 1;//this is needed only for the printouts
493  // std::cout << "p1 = " << p1 << "\t" << std::string::npos << "\t" << category.size() << std::endl;
494  //std::cout << etmin << "\t" << etmax << "\t" << category.substr(p1+1, p2-p1-1) << "\t" << p1 << "\t" << p2 << std::endl;
495  if(p1 != std::string::npos) {
496  p1 = category.find("_", p1);
497  p2 = category.find("_", p1 + 1);
498  etmin = TString(category.substr(p1 + 1, p2 - p1 - 1)).Atof();
499  p1 = p2;
500  p2 = category.find("-", p1);
501  etmax = TString(category.substr(p1 + 1, p2 - p1 - 1)).Atof();
502  // std::cout << etmin << "\t" << etmax << "\t" << category.substr(p1 + 1, p2 - p1 - 1) << std::endl;
503  }
504 
505  if(category.find("gold") != std::string::npos ||
506  category.find("Gold") != std::string::npos ||
507  category.find("highR9") != std::string::npos) {
508  r9min = 0.94;
509  r9max = FLT_MAX;
510  } else if(category.find("bad") != std::string::npos ||
511  category.find("Bad") != std::string::npos ||
512  category.find("lowR9") != std::string::npos
513  ) {
514  r9min = -1;
515  r9max = 0.94;
516  };
517 }
float etmin
min Et value for the bin
float alpha
Definition: AMPTWrapper.h:95
float getSmearingSigma(int runNumber, bool isEBEle, float R9Ele, float etaSCEle, float EtEle, paramSmear_t par, float nSigma=0.) const
float r9max
max R9 value for the bin
void AddScale(TString category_, int runMin_, int runMax_, double deltaP_, double err_deltaP_, double err_syst_deltaP)
float getScaleSystUncertainty(unsigned int runNumber, bool isEBEle, double R9Ele, double etaSCEle, double EtEle) const
float ScaleCorrectionUncertainty(unsigned int runNumber, bool isEBEle, double R9Ele, double etaSCEle, double EtEle) const
method to get scale correction uncertainties: it&#39;s stat+syst in eta x R9 categories ...
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
void ReadFromFile(TString filename)
category "runNumber" runMin runMax deltaP err_deltaP_per_bin err_deltaP_stat err_deltaP_syst ...
#define M_PI_2
float ScaleCorrection(unsigned int runNumber, bool isEBEle, double R9Ele, double etaSCEle, double EtEle) const
method to get energy scale corrections
correctionValue_class getScaleCorrection(unsigned int runNumber, bool isEBEle, double R9Ele, double etaSCEle, double EtEle) const
returns the correction value class
correctionCategory_class(TString category_)
constructor with name of the category according to ElectronCategory_class
float etmax
max Et value for the bin
def cat(path)
Definition: eostools.py:401
bool operator<(const correctionCategory_class &b) const
for ordering of the categories
void AddSmearing(TString category_, int runMin_, int runMax_, double rho, double err_rho, double phi, double err_phi, double Emean, double err_Emean)
T sqrt(T t)
Definition: SSEVec.h:18
~EnergyScaleCorrection_class(void)
dummy constructor needed in ElectronEnergyCalibratorRun2
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
float getScaleOffset(unsigned int runNumber, bool isEBEle, double R9Ele, double etaSCEle, double EtEle) const
double p2[4]
Definition: TauolaWrapper.h:90
float getSmearingRho(int runNumber, bool isEBEle, float R9Ele, float etaSCEle, float EtEle) const
public for sigmaE estimate
JetCorrectorParameters corr
Definition: classes.h:5
float etamin
min eta value for the bin
double b
Definition: hdecay.h:120
float r9min
min R9 vaule for the bin
float getScaleStatUncertainty(unsigned int runNumber, bool isEBEle, double R9Ele, double etaSCEle, double EtEle) const
void ReadSmearingFromFile(TString filename)
File structure: category constTerm alpha;.
double p1[4]
Definition: TauolaWrapper.h:89
float etamax
max eta value for the bin