CMS 3D CMS Logo

BTagCalibrationReader.cc
Go to the documentation of this file.
3 
4 
6 {
7  friend class BTagCalibrationReader;
8 
9 public:
10  struct TmpEntry {
11  float etaMin;
12  float etaMax;
13  float ptMin;
14  float ptMax;
15  float discrMin;
16  float discrMax;
17  TF1 func;
18  };
19 
20 private:
22  const std::string & sysType,
23  const std::vector<std::string> & otherSysTypes={});
24 
25  void load(const BTagCalibration & c,
27  std::string measurementType);
28 
29  double eval(BTagEntry::JetFlavor jf,
30  float eta,
31  float pt,
32  float discr) const;
33 
34  double eval_auto_bounds(const std::string & sys,
36  float eta,
37  float pt,
38  float discr) const;
39 
40  std::pair<float, float> min_max_pt(BTagEntry::JetFlavor jf,
41  float eta,
42  float discr) const;
43 
44  std::pair<float, float> min_max_eta(BTagEntry::JetFlavor jf,
45  float discr) const;
46 
49  std::vector<std::vector<TmpEntry> > tmpData_; // first index: jetFlavor
50  std::vector<bool> useAbsEta_; // first index: jetFlavor
51  std::map<std::string, std::shared_ptr<BTagCalibrationReaderImpl>> otherSysTypeReaders_;
52 };
53 
54 
57  const std::string & sysType,
58  const std::vector<std::string> & otherSysTypes):
59  op_(op),
60  sysType_(sysType),
61  tmpData_(3),
62  useAbsEta_(3, true)
63 {
64  for (const std::string & ost : otherSysTypes) {
65  if (otherSysTypeReaders_.count(ost)) {
66  throw cms::Exception("BTagCalibrationReader")
67  << "Every otherSysType should only be given once. Duplicate: "
68  << ost;
69  }
70  otherSysTypeReaders_[ost] = std::unique_ptr<BTagCalibrationReaderImpl>(
71  new BTagCalibrationReaderImpl(op, ost)
72  );
73  }
74 }
75 
77  const BTagCalibration & c,
79  std::string measurementType)
80 {
81  if (!tmpData_[jf].empty()) {
82  throw cms::Exception("BTagCalibrationReader")
83  << "Data for this jet-flavor is already loaded: "
84  << jf;
85  }
86 
87  BTagEntry::Parameters params(op_, measurementType, sysType_);
88  const std::vector<BTagEntry> &entries = c.getEntries(params);
89 
90  for (const auto &be : entries) {
91  if (be.params.jetFlavor != jf) {
92  continue;
93  }
94 
95  TmpEntry te;
96  te.etaMin = be.params.etaMin;
97  te.etaMax = be.params.etaMax;
98  te.ptMin = be.params.ptMin;
99  te.ptMax = be.params.ptMax;
100  te.discrMin = be.params.discrMin;
101  te.discrMax = be.params.discrMax;
102 
103  if (op_ == BTagEntry::OP_RESHAPING) {
104  te.func = TF1("", be.formula.c_str(),
105  be.params.discrMin, be.params.discrMax);
106  } else {
107  te.func = TF1("", be.formula.c_str(),
108  be.params.ptMin, be.params.ptMax);
109  }
110 
111  tmpData_[be.params.jetFlavor].push_back(te);
112  if (te.etaMin < 0) {
113  useAbsEta_[be.params.jetFlavor] = false;
114  }
115  }
116 
117  for (auto & p : otherSysTypeReaders_) {
118  p.second->load(c, jf, measurementType);
119  }
120 }
121 
124  float eta,
125  float pt,
126  float discr) const
127 {
128  bool use_discr = (op_ == BTagEntry::OP_RESHAPING);
129  if (useAbsEta_[jf] && eta < 0) {
130  eta = -eta;
131  }
132 
133  // search linearly through eta, pt and discr ranges and eval
134  // future: find some clever data structure based on intervals
135  const auto &entries = tmpData_.at(jf);
136  for (unsigned i=0; i<entries.size(); ++i) {
137  const auto &e = entries.at(i);
138  if (
139  e.etaMin <= eta && eta <= e.etaMax // find eta
140  && e.ptMin < pt && pt <= e.ptMax // check pt
141  ){
142  if (use_discr) { // discr. reshaping?
143  if (e.discrMin <= discr && discr < e.discrMax) { // check discr
144  return e.func.Eval(discr);
145  }
146  } else {
147  return e.func.Eval(pt);
148  }
149  }
150  }
151 
152  return 0.; // default value
153 }
154 
156  const std::string & sys,
158  float eta,
159  float pt,
160  float discr) const
161 {
162  auto sf_bounds_eta = min_max_eta(jf, discr);
163  bool eta_is_out_of_bounds = false;
164 
165  if (sf_bounds_eta.first < 0) sf_bounds_eta.first = -sf_bounds_eta.second;
166 
167  if (useAbsEta_[jf] && eta < 0) {
168  eta = -eta;
169  }
170 
171  if (eta <= sf_bounds_eta.first || eta > sf_bounds_eta.second ) {
172  eta_is_out_of_bounds = true;
173  }
174 
175  if (eta_is_out_of_bounds) {
176  return 1.;
177  }
178 
179 
180  auto sf_bounds = min_max_pt(jf, eta, discr);
181  float pt_for_eval = pt;
182  bool is_out_of_bounds = false;
183 
184  if (pt <= sf_bounds.first) {
185  pt_for_eval = sf_bounds.first + .0001;
186  is_out_of_bounds = true;
187  } else if (pt > sf_bounds.second) {
188  pt_for_eval = sf_bounds.second - .0001;
189  is_out_of_bounds = true;
190  }
191 
192  // get central SF (and maybe return)
193  double sf = eval(jf, eta, pt_for_eval, discr);
194  if (sys == sysType_) {
195  return sf;
196  }
197 
198  // get sys SF (and maybe return)
199  if (!otherSysTypeReaders_.count(sys)) {
200  throw cms::Exception("BTagCalibrationReader")
201  << "sysType not available (maybe not loaded?): "
202  << sys;
203  }
204  double sf_err = otherSysTypeReaders_.at(sys)->eval(jf, eta, pt_for_eval, discr);
205  if (!is_out_of_bounds) {
206  return sf_err;
207  }
208 
209  // double uncertainty on out-of-bounds and return
210  sf_err = sf + 2*(sf_err - sf);
211  return sf_err;
212 }
213 
216  float eta,
217  float discr) const
218 {
219  bool use_discr = (op_ == BTagEntry::OP_RESHAPING);
220  if (useAbsEta_[jf] && eta < 0) {
221  eta = -eta;
222  }
223 
224  const auto &entries = tmpData_.at(jf);
225  float min_pt = -1., max_pt = -1.;
226  for (const auto & e: entries) {
227  if (
228  e.etaMin <= eta && eta <=e.etaMax // find eta
229  ){
230  if (min_pt < 0.) { // init
231  min_pt = e.ptMin;
232  max_pt = e.ptMax;
233  continue;
234  }
235 
236  if (use_discr) { // discr. reshaping?
237  if (e.discrMin <= discr && discr < e.discrMax) { // check discr
238  min_pt = min_pt < e.ptMin ? min_pt : e.ptMin;
239  max_pt = max_pt > e.ptMax ? max_pt : e.ptMax;
240  }
241  } else {
242  min_pt = min_pt < e.ptMin ? min_pt : e.ptMin;
243  max_pt = max_pt > e.ptMax ? max_pt : e.ptMax;
244  }
245  }
246  }
247 
248  return std::make_pair(min_pt, max_pt);
249 }
250 
253  float discr) const
254 {
255  bool use_discr = (op_ == BTagEntry::OP_RESHAPING);
256 
257  const auto &entries = tmpData_.at(jf);
258  float min_eta = 0., max_eta = 0.;
259  for (const auto & e: entries) {
260 
261  if (use_discr) { // discr. reshaping?
262  if (e.discrMin <= discr && discr < e.discrMax) { // check discr
263  min_eta = min_eta < e.etaMin ? min_eta : e.etaMin;
264  max_eta = max_eta > e.etaMax ? max_eta : e.etaMax;
265  }
266  } else {
267  min_eta = min_eta < e.etaMin ? min_eta : e.etaMin;
268  max_eta = max_eta > e.etaMax ? max_eta : e.etaMax;
269  }
270  }
271 
272 
273  return std::make_pair(min_eta, max_eta);
274 }
275 
276 
278  const std::string & sysType,
279  const std::vector<std::string> & otherSysTypes):
280  pimpl(new BTagCalibrationReaderImpl(op, sysType, otherSysTypes)) {}
281 
284  const std::string & measurementType)
285 {
286  pimpl->load(c, jf, measurementType);
287 }
288 
290  float eta,
291  float pt,
292  float discr) const
293 {
294  return pimpl->eval(jf, eta, pt, discr);
295 }
296 
299  float eta,
300  float pt,
301  float discr) const
302 {
303  return pimpl->eval_auto_bounds(sys, jf, eta, pt, discr);
304 }
305 
307  float eta,
308  float discr) const
309 {
310  return pimpl->min_max_pt(jf, eta, discr);
311 }
312 
313 
std::vector< std::vector< TmpEntry > > tmpData_
const std::vector< BTagEntry > & getEntries(const BTagEntry::Parameters &par) const
double eval_auto_bounds(const std::string &sys, BTagEntry::JetFlavor jf, float eta, float pt, float discr) const
std::shared_ptr< BTagCalibrationReaderImpl > pimpl
OperatingPoint
Definition: BTagEntry.h:27
double eval(BTagEntry::JetFlavor jf, float eta, float pt, float discr=0.) const
BTagCalibrationReaderImpl(BTagEntry::OperatingPoint op, const std::string &sysType, const std::vector< std::string > &otherSysTypes={})
void load(const BTagCalibration &c, BTagEntry::JetFlavor jf, const std::string &measurementType="comb")
std::pair< float, float > min_max_pt(BTagEntry::JetFlavor jf, float eta, float discr=0.) const
double eval_auto_bounds(const std::string &sys, BTagEntry::JetFlavor jf, float eta, float pt, float discr=0.) const
std::pair< float, float > min_max_eta(BTagEntry::JetFlavor jf, float discr) const
std::map< std::string, std::shared_ptr< BTagCalibrationReaderImpl > > otherSysTypeReaders_
std::pair< float, float > min_max_pt(BTagEntry::JetFlavor jf, float eta, float discr) const
double eval(BTagEntry::JetFlavor jf, float eta, float pt, float discr) const
void load(const BTagCalibration &c, BTagEntry::JetFlavor jf, std::string measurementType)