CMS 3D CMS Logo

EnergyScaleCorrection.cc
Go to the documentation of this file.
2 
5 
6 #include <cstdlib>
7 #include <cmath>
8 #include <iomanip>
9 #include <algorithm>
10 #include <sstream>
11 #include <iterator>
12 
13 EnergyScaleCorrection::EnergyScaleCorrection(const std::string& correctionFileName, unsigned int genSeed)
14  : smearingType_(ECALELF) {
15  if (!correctionFileName.empty()) {
16  std::string filename = correctionFileName + "_scales.dat";
17  readScalesFromFile(filename);
18  if (scales_.empty()) {
19  throw cms::Exception("EnergyScaleCorrection") << "scale correction map empty";
20  }
21  }
22 
23  if (!correctionFileName.empty()) {
24  std::string filename = correctionFileName + "_smearings.dat";
25  readSmearingsFromFile(filename);
26  if (smearings_.empty()) {
27  throw cms::Exception("EnergyScaleCorrection") << "smearing correction map empty";
28  }
29  }
30 }
31 
33  double et,
34  double eta,
35  double r9,
36  unsigned int gainSeed,
37  std::bitset<kErrNrBits> uncBitMask) const {
38  const ScaleCorrection* scaleCorr = getScaleCorr(runNumber, et, eta, r9, gainSeed);
39  if (scaleCorr != nullptr)
40  return scaleCorr->scale();
41  else
42  return kDefaultScaleVal_;
43 }
44 
46  double et,
47  double eta,
48  double r9,
49  unsigned int gainSeed,
50  std::bitset<kErrNrBits> uncBitMask) const {
51  const ScaleCorrection* scaleCorr = getScaleCorr(runNumber, et, eta, r9, gainSeed);
52  if (scaleCorr != nullptr)
53  return scaleCorr->scaleErr(uncBitMask);
54  else
55  return 0.;
56 }
57 
59  int runnr, double et, double eta, double r9, unsigned int gainSeed, ParamSmear par, float nSigma) const {
60  if (par == kRho)
61  return smearingSigma(runnr, et, eta, r9, gainSeed, nSigma, 0.);
62  if (par == kPhi)
63  return smearingSigma(runnr, et, eta, r9, gainSeed, 0., nSigma);
64  return smearingSigma(runnr, et, eta, r9, gainSeed, 0., 0.);
65 }
66 
68  int runnr, double et, double eta, double r9, unsigned int gainSeed, float nrSigmaRho, float nrSigmaPhi) const {
69  const SmearCorrection* smearCorr = getSmearCorr(runnr, et, eta, r9, gainSeed);
70 
71  if (smearCorr != nullptr)
72  return smearCorr->sigma(nrSigmaRho, nrSigmaPhi);
73  else
74  return kDefaultSmearVal_;
75 }
76 
78  unsigned int runnr, double et, double eta, double r9, unsigned int gainSeed) const {
79  // buld the category based on the values of the object
80  CorrectionCategory category(runnr, et, eta, r9, gainSeed);
81  auto result = scales_.find(category);
82 
83  if (result == scales_.end()) {
84  edm::LogInfo("EnergyScaleCorrection")
85  << "Scale category not found: " << category << " Returning uncorrected value.";
86  return nullptr;
87  }
88 
89  //validate the result, just to be sure
90  if (!result->first.inCategory(runnr, et, eta, r9, gainSeed)) {
91  throw cms::Exception("LogicError") << " error found scale category " << result->first
92  << " that does not contain run " << runnr << " et " << et << " eta " << eta
93  << " r9 " << r9 << " gain seed " << gainSeed;
94  }
95  return &result->second;
96 }
97 
99  unsigned int runnr, double et, double eta, double r9, unsigned int gainSeed) const {
100  // buld the category based on the values of the object
101  CorrectionCategory category(runnr, et, eta, r9, gainSeed);
102  auto result = smearings_.find(category);
103 
104  if (result == smearings_.end()) {
105  edm::LogInfo("EnergyScaleCorrection")
106  << "Smear category not found: " << category << " Returning uncorrected value.";
107  return nullptr;
108  }
109 
110  //validate the result, just to be sure
111  if (!result->first.inCategory(runnr, et, eta, r9, gainSeed)) {
112  throw cms::Exception("LogicError") << " error found smear category " << result->first
113  << " that does not contain run " << runnr << " et " << et << " eta " << eta
114  << " r9 " << r9 << " gain seed " << gainSeed;
115  }
116  return &result->second;
117 }
118 
120  int runMin,
121  int runMax,
122  double energyScale,
123  double energyScaleErrStat,
124  double energyScaleErrSyst,
125  double energyScaleErrGain) {
126  CorrectionCategory cat(category, runMin, runMax); // build the category from the string
127  auto result = scales_.equal_range(cat);
128  if (result.first != result.second) {
129  throw cms::Exception("ConfigError") << "Category already defined! " << cat;
130  }
131 
132  ScaleCorrection corr(energyScale, energyScaleErrStat, energyScaleErrSyst, energyScaleErrGain);
133  scales_.insert(result.first, {cat, corr}); //use a hint where to insert
134 }
135 
137  int runMax,
138  double etaMin,
139  double etaMax,
140  double r9Min,
141  double r9Max,
142  double etMin,
143  double etMax,
144  unsigned int gain,
145  double energyScale,
146  double energyScaleErrStat,
147  double energyScaleErrSyst,
148  double energyScaleErrGain) {
149  CorrectionCategory cat(runMin, runMax, etaMin, etaMax, r9Min, r9Max, etMin, etMax, gain);
150 
151  auto result = scales_.equal_range(cat);
152  if (result.first != result.second) {
153  throw cms::Exception("ConfigError") << "Category already defined! " << cat;
154  }
155 
156  ScaleCorrection corr(energyScale, energyScaleErrStat, energyScaleErrSyst, energyScaleErrGain);
157  scales_.insert(result.first, {cat, corr});
158 }
159 
161  int runMin,
162  int runMax,
163  double rho,
164  double errRho,
165  double phi,
166  double errPhi,
167  double eMean,
168  double errEMean) {
169  CorrectionCategory cat(category);
170  auto res = smearings_.equal_range(cat);
171 
172  if (res.first != res.second) {
173  throw cms::Exception("EnergyScaleCorrection") << "Smearing category already defined " << cat;
174  }
175 
176  SmearCorrection corr(rho, errRho, phi, errPhi, eMean, errEMean);
177  smearings_.insert(res.first, {cat, corr}); //use a hint from res
178 }
179 
181  if (value <= 1) {
183  } else {
185  }
186 }
187 
189  std::ifstream file(edm::FileInPath(filename).fullPath().c_str());
190 
191  if (!file.good()) {
192  throw cms::Exception("EnergyScaleCorrection") << "file " << filename << " not readable.";
193  }
194 
195  int runMin, runMax;
196  std::string category, region2;
197  double energyScale, energyScaleErr, energyScaleErrStat, energyScaleErrSyst, energyScaleErrGain;
198 
199  double etaMin;
200  double etaMax;
201  double r9Min;
202  double r9Max;
203  double etMin;
204  double etMax;
205  unsigned int gain;
206 
207  // TO count the #columns in that txt file and decide based on that the version to read
209  std::stringstream stream;
210  getline(file, line);
211  stream.clear();
212  stream << line;
213 
214  int ncols = std::distance(std::istream_iterator<std::string>(stream), std::istream_iterator<std::string>());
215 
216  file.seekg(0, std::ios::beg);
217 
218  if (ncols == 9) {
219  for (file >> category; file.good(); file >> category) {
220  file >> region2 >> runMin >> runMax >> energyScale >> energyScaleErr >> energyScaleErrStat >>
221  energyScaleErrSyst >> energyScaleErrGain;
222  addScale(category, runMin, runMax, energyScale, energyScaleErrStat, energyScaleErrSyst, energyScaleErrGain);
223  }
224  } else {
225  if (file.peek() == 'r')
226  file.ignore(1000, 10);
227 
228  for (file >> runMin; file.good(); file >> runMin) {
229  file >> runMax >> etaMin >> etaMax >> r9Min >> r9Max >> etMin >> etMax >> gain >> energyScale >> energyScaleErr;
230  file.ignore(1000, 10);
231  energyScaleErrStat = energyScaleErr;
232  energyScaleErrSyst = 0;
233  energyScaleErrGain = 0;
234 
235  addScale(runMin,
236  runMax,
237  etaMin,
238  etaMax,
239  r9Min,
240  r9Max,
241  etMin,
242  etMax,
243  gain,
244  energyScale,
245  energyScaleErrStat,
246  energyScaleErrSyst,
247  energyScaleErrGain);
248  }
249  }
250 
251  file.close();
252  return;
253 }
254 
255 //also more or less untouched function from the orginal package
257  std::ifstream file(edm::FileInPath(filename).fullPath().c_str());
258  if (!file.good()) {
259  throw cms::Exception("EnergyScaleCorrection") << "file " << filename << " not readable";
260  }
261 
262  int runMin = 0;
263  int runMax = 900000;
264  int unused = 0;
265  std::string category, region2;
266  double rho, phi, eMean, errRho, errPhi, errEMean;
267  double etaMin, etaMax, r9Min, r9Max;
268  std::string phiString, errPhiString;
269 
270  while (file.peek() != EOF && file.good()) {
271  if (file.peek() == 10) { // 10 = \n
272  file.get();
273  continue;
274  }
275 
276  if (file.peek() == 35) { // 35 = #
277  file.ignore(1000, 10); // ignore the rest of the line until \n
278  continue;
279  }
280 
281  if (smearingType_ == UNKNOWN) { // trying to guess: not recommended
282  throw cms::Exception("ConfigError") << "unknown smearing type";
283 
284  } else if (smearingType_ == GLOBE) {
285  file >> category >> unused >> etaMin >> etaMax >> r9Min >> r9Max >> runMin >> runMax >> eMean >> errEMean >>
286  rho >> errRho >> phi >> errPhi;
287 
288  addSmearing(category, runMin, runMax, rho, errRho, phi, errPhi, eMean, errEMean);
289 
290  } else if (smearingType_ == ECALELF) {
291  file >> category >> eMean >> errEMean >> rho >> errRho >> phiString >> errPhiString;
292 
293  if (phiString == "M_PI_2")
294  phi = M_PI_2;
295  else
296  phi = std::stod(phiString);
297 
298  if (errPhiString == "M_PI_2")
299  errPhi = M_PI_2;
300  else
301  errPhi = std::stod(errPhiString);
302 
303  addSmearing(category, runMin, runMax, rho, errRho, phi, errPhi, eMean, errEMean);
304 
305  } else {
306  file >> category >> rho >> phi;
307  errRho = errPhi = eMean = errEMean = 0;
308  addSmearing(category, runMin, runMax, rho, errRho, phi, errPhi, eMean, errEMean);
309  }
310  }
311 
312  file.close();
313  return;
314 }
315 
316 std::ostream& EnergyScaleCorrection::ScaleCorrection::print(std::ostream& os) const {
317  os << "( " << scale_ << " +/- " << scaleErrStat_ << " +/- " << scaleErrSyst_ << " +/- " << scaleErrGain_ << ")";
318  return os;
319 }
320 
321 float EnergyScaleCorrection::ScaleCorrection::scaleErr(const std::bitset<kErrNrBits>& uncBitMask) const {
322  double totErr(0);
323  auto pow2 = [](const double& x) { return x * x; };
324 
325  if (uncBitMask.test(kErrStatBitNr))
326  totErr += pow2(scaleErrStat_);
327  if (uncBitMask.test(kErrSystBitNr))
328  totErr += pow2(scaleErrSyst_);
329  if (uncBitMask.test(kErrGainBitNr))
330  totErr += pow2(scaleErrGain_);
331 
332  return std::sqrt(totErr);
333 }
334 
335 std::ostream& EnergyScaleCorrection::SmearCorrection::print(std::ostream& os) const {
336  os << rho_ << " +/- " << rhoErr_ << "\t" << phi_ << " +/- " << phiErr_ << "\t" << eMean_ << " +/- " << eMeanErr_;
337  return os;
338 }
339 
340 //here be dragons
341 //this function is nasty and needs to be replaced
343  : runMin_(runnrMin),
344  runMax_(runnrMax),
345  etaMin_(0),
346  etaMax_(3),
347  r9Min_(-1),
348  r9Max_(999),
349  etMin_(0),
350  etMax_(9999999),
351  gain_(0) {
352  size_t p1, p2; // boundary
353 
354  // eta region
355  p1 = category.find("absEta_");
356  if (category.find("absEta_0_1") != std::string::npos) {
357  etaMin_ = 0;
358  etaMax_ = 1;
359  } else if (category.find("absEta_1_1.4442") != std::string::npos) {
360  etaMin_ = 1;
361  etaMax_ = 1.479;
362  } else if (category.find("absEta_1.566_2") != std::string::npos) {
363  etaMin_ = 1.479;
364  etaMax_ = 2;
365  } else if (category.find("absEta_2_2.5") != std::string::npos) {
366  etaMin_ = 2;
367  etaMax_ = 3;
368  } else {
369  if (p1 != std::string::npos) {
370  p1 = category.find("_", p1);
371  p2 = category.find("_", p1 + 1);
372  etaMin_ = std::stof(category.substr(p1 + 1, p2 - p1 - 1));
373  p1 = p2;
374  p2 = category.find("-", p1);
375  etaMax_ = std::stof(category.substr(p1 + 1, p2 - p1 - 1));
376  }
377  }
378 
379  if (category.find("EBlowEta") != std::string::npos) {
380  etaMin_ = 0;
381  etaMax_ = 1;
382  };
383  if (category.find("EBhighEta") != std::string::npos) {
384  etaMin_ = 1;
385  etaMax_ = 1.479;
386  };
387  if (category.find("EElowEta") != std::string::npos) {
388  etaMin_ = 1.479;
389  etaMax_ = 2;
390  };
391  if (category.find("EEhighEta") != std::string::npos) {
392  etaMin_ = 2;
393  etaMax_ = 7;
394  };
395 
396  // Et region
397  p1 = category.find("-Et_");
398 
399  if (p1 != std::string::npos) {
400  p1 = category.find("_", p1);
401  p2 = category.find("_", p1 + 1);
402  etMin_ = std::stof(category.substr(p1 + 1, p2 - p1 - 1));
403  p1 = p2;
404  p2 = category.find("-", p1);
405  etMax_ = std::stof(category.substr(p1 + 1, p2 - p1 - 1));
406  }
407 
408  if (category.find("gold") != std::string::npos || category.find("Gold") != std::string::npos ||
409  category.find("highR9") != std::string::npos) {
410  r9Min_ = 0.94;
412  } else if (category.find("bad") != std::string::npos || category.find("Bad") != std::string::npos ||
413  category.find("lowR9") != std::string::npos) {
414  r9Min_ = -1;
415  r9Max_ = 0.94;
416  };
417  // R9 region
418  p1 = category.find("-R9");
419  if (p1 != std::string::npos) {
420  p1 = category.find("_", p1);
421  p2 = category.find("_", p1 + 1);
422  r9Min_ = std::stof(category.substr(p1 + 1, p2 - p1 - 1));
423  // If there is one value, just set lower bound
424  if (p2 != std::string::npos) {
425  p1 = p2;
426  p2 = category.find("-", p1);
427  r9Max_ = std::stof(category.substr(p1 + 1, p2 - p1 - 1));
428  if (r9Max_ >= 1.0)
430  }
431  }
432  //------------------------------
433  p1 = category.find("gainEle_"); // Position of first character
434  if (p1 != std::string::npos) {
435  p1 += 8; // Position of character after _
436  p2 = category.find("-", p1); // Position of - or end of string
437  gain_ = std::stoul(category.substr(p1, p2 - p1), nullptr);
438  }
439  //so turns out the code does an inclusive X<=Y<=Z search for bins
440  //which is what we want for run numbers
441  //however then the problem is when we get a value exactly at the bin boundary
442  //for the et/eta/r9 which then gives multiple bins
443  //so we just decrement the maxValues ever so slightly to ensure that they are different
444  //from the next bins min value
445  etMax_ = std::nextafterf(etMax_, std::numeric_limits<float>::min());
446  etaMax_ = std::nextafterf(etaMax_, std::numeric_limits<float>::min());
447  r9Max_ = std::nextafterf(r9Max_, std::numeric_limits<float>::min());
448 }
449 
452  unsigned int runMax,
453  float etaMin,
454  float etaMax,
455  float r9Min,
456  float r9Max,
457  float etMin,
458  float etMax,
459  unsigned int gainSeed)
460  : runMin_(runMin),
461  runMax_(runMax),
462  etaMin_(etaMin),
463  etaMax_(etaMax),
464  r9Min_(r9Min),
465  r9Max_(r9Max),
466  etMin_(etMin),
467  etMax_(etMax),
468  gain_(gainSeed) {
473  etMax_ = std::nextafterf(etMax_, std::numeric_limits<float>::min());
474  etaMax_ = std::nextafterf(etaMax_, std::numeric_limits<float>::min());
475  r9Max_ = std::nextafterf(r9Max_, std::numeric_limits<float>::min());
476 };
477 
479  const unsigned int runnr, const float et, const float eta, const float r9, const unsigned int gainSeed) const {
480  return runnr >= runMin_ && runnr <= runMax_ && et >= etMin_ && et <= etMax_ && eta >= etaMin_ && eta <= etaMax_ &&
481  r9 >= r9Min_ && r9 <= r9Max_ && (gain_ == 0 || gainSeed == gain_);
482 }
483 
485  if (runMin_ < b.runMin_ && runMax_ < b.runMax_)
486  return true;
487  if (runMax_ > b.runMax_ && runMin_ > b.runMin_)
488  return false;
489 
490  if (etaMin_ < b.etaMin_ && etaMax_ < b.etaMax_)
491  return true;
492  if (etaMax_ > b.etaMax_ && etaMin_ > b.etaMin_)
493  return false;
494 
495  if (r9Min_ < b.r9Min_ && r9Max_ < b.r9Max_)
496  return true;
497  if (r9Max_ > b.r9Max_ && r9Min_ > b.r9Min_)
498  return false;
499 
500  if (etMin_ < b.etMin_ && etMax_ < b.etMax_)
501  return true;
502  if (etMax_ > b.etMax_ && etMin_ > b.etMin_)
503  return false;
504 
505  if (gain_ == 0 || b.gain_ == 0)
506  return false; // if corrections are not categorized in gain then default gain value should always return false in order to have a match with the category
507  if (gain_ < b.gain_)
508  return true;
509  else
510  return false;
511  return false;
512 }
513 
514 std::ostream& EnergyScaleCorrection::CorrectionCategory::print(std::ostream& os) const {
515  os << runMin_ << " " << runMax_ << "\t" << etaMin_ << " " << etaMax_ << "\t" << r9Min_ << " " << r9Max_ << "\t"
516  << etMin_ << " " << etMax_ << "\t" << gain_;
517  return os;
518 }
void readSmearingsFromFile(const std::string &filename)
static constexpr float kDefaultSmearVal_
void addScale(const std::string &category, int runMin, int runMax, double deltaP, double errDeltaP, double errSystDeltaP, double errDeltaPGain)
CorrectionCategory(const std::string &category, int runnrMin=0, int runnrMax=999999)
#define M_PI_2
std::ostream & print(std::ostream &os) const
Definition: Electron.h:6
bool inCategory(const unsigned int runnr, const float et, const float eta, const float r9, const unsigned int gainSeed) const
void readScalesFromFile(const std::string &filename)
def cat(path)
Definition: eostools.py:401
const ScaleCorrection * getScaleCorr(unsigned int runnr, double et, double eta, double r9, unsigned int gainSeed) const
T sqrt(T t)
Definition: SSEVec.h:18
void setSmearingType(FileFormat value)
float smearingSigma(int runnr, double et, double eta, double r9, unsigned int gainSeed, ParamSmear par, float nSigma=0.) const
static constexpr float kDefaultScaleVal_
const SmearCorrection * getSmearCorr(unsigned int runnr, double et, double eta, double r9, unsigned int gainSeed) const
unsigned int gain_
12, 6, 1, 61 (double gain switch)
std::map< CorrectionCategory, SmearCorrection > smearings_
Definition: value.py:1
T min(T a, T b)
Definition: MathUtil.h:58
double p2[4]
Definition: TauolaWrapper.h:90
JetCorrectorParameters corr
Definition: classes.h:5
float scaleErr(const std::bitset< kErrNrBits > &uncBitMask) const
float sigma(const float et, const float nrSigmaRho=0., const float nrSigmaPhi=0.) const
void addSmearing(const std::string &category, int runMin, int runMax, double rho, double errRho, double phi, double errPhi, double eMean, double errEMean)
double b
Definition: hdecay.h:120
et
define resolution functions of each parameter
double p1[4]
Definition: TauolaWrapper.h:89
std::ostream & print(std::ostream &os) const
std::map< CorrectionCategory, ScaleCorrection > scales_
static const char gain_[]
bool operator<(const CorrectionCategory &b) const
std::ostream & print(std::ostream &os) const
float scaleCorrUncert(unsigned int runnr, double et, double eta, double r9, unsigned int gainSeed, std::bitset< kErrNrBits > uncBitMask=kErrNone) const
float scaleCorr(unsigned int runnr, double et, double eta, double r9, unsigned int gainSeed=12, std::bitset< kErrNrBits > uncBitMask=kErrNone) const