rivet is hosted by Hepforge, IPPP Durham
Rivet  2.7.0
SmearingFunctions.hh
1 // -*- C++ -*-
2 #ifndef RIVET_SmearingFunctions_HH
3 #define RIVET_SmearingFunctions_HH
4 
5 #include "Rivet/Tools/MomentumSmearingFunctions.hh"
6 #include "Rivet/Tools/ParticleSmearingFunctions.hh"
7 #include "Rivet/Tools/JetSmearingFunctions.hh"
8 
9 namespace Rivet {
10 
11 
13 
14 
18  inline double ELECTRON_EFF_ATLAS_RUN1(const Particle& e) {
19  if (e.abseta() > 2.5) return 0;
20  if (e.pT() < 10*GeV) return 0;
21  return (e.abseta() < 1.5) ? 0.95 : 0.85;
22  }
23 
26  inline double ELECTRON_EFF_ATLAS_RUN2(const Particle& e) {
27  return ELECTRON_EFF_ATLAS_RUN1(e);
28  }
29 
30 
35  inline double ELECTRON_IDEFF_ATLAS_RUN2_LOOSE(const Particle& e) {
36 
37  // Manually symmetrised eta eff histogram
38  const static vector<double> edges_eta = { 0.0, 0.1, 0.8, 1.37, 1.52, 2.01, 2.37, 2.47 };
39  const static vector<double> effs_eta = { 0.950, 0.965, 0.955, 0.885, 0.950, 0.935, 0.90 };
40  // Et eff histogram (10-20 is a guess)
41  const static vector<double> edges_et = { 0, 10, 20, 25, 30, 35, 40, 45, 50, 60, 80 };
42  const static vector<double> effs_et = { 0.0, 0.90, 0.91, 0.92, 0.94, 0.95, 0.955, 0.965, 0.97, 0.98 };
43 
44  if (e.abseta() > 2.47) return 0.0; // no ID outside the tracker
45 
46  const int i_eta = binIndex(e.abseta(), edges_eta);
47  const int i_et = binIndex(e.Et()/GeV, edges_et, true);
48  const double eff = effs_et[i_et] * effs_eta[i_eta] / 0.95; //< norm factor as approximate double differential
49  return min(eff, 1.0);
50  }
51 
52 
54  inline double ELECTRON_IDEFF_ATLAS_RUN1_MEDIUM(const Particle& e) {
55 
56  const static vector<double> eta_edges_10 = {0.000, 0.049, 0.454, 1.107, 1.46, 1.790, 2.277, 2.500};
57  const static vector<double> eta_vals_10 = {0.730, 0.757, 0.780, 0.771, 0.77, 0.777, 0.778};
58 
59  const static vector<double> eta_edges_15 = {0.000, 0.053, 0.456, 1.102, 1.463, 1.783, 2.263, 2.500};
60  const static vector<double> eta_vals_15 = {0.780, 0.800, 0.819, 0.759, 0.749, 0.813, 0.829};
61 
62  const static vector<double> eta_edges_20 = {0.000, 0.065, 0.362, 0.719, 0.980, 1.289, 1.455, 1.681, 1.942, 2.239, 2.452, 2.500};
63  const static vector<double> eta_vals_20 = {0.794, 0.806, 0.816, 0.806, 0.797, 0.774, 0.764, 0.788, 0.793, 0.806, 0.825};
64 
65  const static vector<double> eta_edges_25 = {0.000, 0.077, 0.338, 0.742, 1.004, 1.265, 1.467, 1.692, 1.940, 2.227, 2.452, 2.500};
66  const static vector<double> eta_vals_25 = {0.833, 0.843, 0.853, 0.845, 0.839, 0.804, 0.790, 0.825, 0.830, 0.833, 0.839};
67 
68  const static vector<double> eta_edges_30 = {0.000, 0.077, 0.350, 0.707, 0.980, 1.289, 1.479, 1.681, 1.942, 2.239, 2.441, 2.500};
69  const static vector<double> eta_vals_30 = {0.863, 0.872, 0.881, 0.874, 0.870, 0.824, 0.808, 0.847, 0.845, 0.840, 0.842};
70 
71  const static vector<double> eta_edges_35 = {0.000, 0.058, 0.344, 0.700, 1.009, 1.270, 1.458, 1.685, 1.935, 2.231, 2.468, 2.500};
72  const static vector<double> eta_vals_35 = {0.878, 0.889, 0.901, 0.895, 0.893, 0.849, 0.835, 0.868, 0.863, 0.845, 0.832};
73 
74  const static vector<double> eta_edges_40 = {0.000, 0.047, 0.355, 0.699, 0.983, 1.280, 1.446, 1.694, 1.943, 2.227, 2.441, 2.500};
75  const static vector<double> eta_vals_40 = {0.894, 0.901, 0.909, 0.905, 0.904, 0.875, 0.868, 0.889, 0.876, 0.848, 0.827};
76 
77  const static vector<double> eta_edges_45 = {0.000, 0.058, 0.356, 0.712, 0.997, 1.282, 1.459, 1.686, 1.935, 2.220, 2.444, 2.500};
78  const static vector<double> eta_vals_45 = {0.900, 0.911, 0.923, 0.918, 0.917, 0.897, 0.891, 0.904, 0.894, 0.843, 0.796};
79 
80  const static vector<double> eta_edges_50 = {0.000, 0.059, 0.355, 0.711, 0.983, 1.280, 1.469, 1.682, 1.919, 2.227, 2.441, 2.500};
81  const static vector<double> eta_vals_50 = {0.903, 0.913, 0.923, 0.922, 0.923, 0.903, 0.898, 0.908, 0.895, 0.831, 0.774};
82 
83  const static vector<double> eta_edges_60 = {0.000, 0.053, 0.351, 0.720, 1.006, 1.291, 1.469, 1.696, 1.946, 2.243, 2.455, 2.500};
84  const static vector<double> eta_vals_60 = {0.903, 0.917, 0.928, 0.924, 0.927, 0.915, 0.911, 0.915, 0.899, 0.827, 0.760};
85 
86  const static vector<double> eta_edges_80 = {0.000, 0.053, 0.351, 0.720, 0.994, 1.292, 1.482, 1.708, 1.934, 2.220, 2.458, 2.500};
87  const static vector<double> eta_vals_80 = {0.936, 0.942, 0.952, 0.956, 0.956, 0.934, 0.931, 0.944, 0.933, 0.940, 0.948};
88 
89  const static vector<double> et_edges = { 10, 15, 20, 25, 30, 35, 40, 45, 50, 60, 80 };
90  const static vector< vector<double> > et_eta_edges = { eta_edges_10, eta_edges_15, eta_edges_20, eta_edges_25, eta_edges_30, eta_edges_35, eta_edges_40, eta_edges_45, eta_edges_50, eta_edges_60, eta_edges_80 };
91  const static vector< vector<double> > et_eta_vals = { eta_vals_10, eta_vals_15, eta_vals_20, eta_vals_25, eta_vals_30, eta_vals_35, eta_vals_40, eta_vals_45, eta_vals_50, eta_vals_60, eta_vals_80 };
92 
93  if (e.abseta() > 2.5 || e.Et() < 10*GeV) return 0.0;
94  const int i_et = binIndex(e.Et()/GeV, et_edges, true);
95  const int i_eta = binIndex(e.abseta(), et_eta_edges[i_et]);
96  return et_eta_vals[i_et][i_eta];
97  }
98 
101  inline double ELECTRON_IDEFF_ATLAS_RUN2_MEDIUM(const Particle& e) {
103  }
104 
105 
107  inline double ELECTRON_IDEFF_ATLAS_RUN1_TIGHT(const Particle& e) {
108 
109  const static vector<double> eta_edges_10 = {0.000, 0.049, 0.459, 1.100, 1.461, 1.789, 2.270, 2.500};
110  const static vector<double> eta_vals_10 = {0.581, 0.632, 0.668, 0.558, 0.548, 0.662, 0.690};
111 
112  const static vector<double> eta_edges_15 = {0.000, 0.053, 0.450, 1.096, 1.463, 1.783, 2.269, 2.500};
113  const static vector<double> eta_vals_15 = {0.630, 0.678, 0.714, 0.633, 0.616, 0.700, 0.733};
114 
115  const static vector<double> eta_edges_20 = {0.000, 0.065, 0.362, 0.719, 0.992, 1.277, 1.479, 1.692, 1.930, 2.227, 2.464, 2.500};
116  const static vector<double> eta_vals_20 = {0.653, 0.695, 0.735, 0.714, 0.688, 0.635, 0.625, 0.655, 0.680, 0.691, 0.674};
117 
118  const static vector<double> eta_edges_25 = {0.000, 0.077, 0.362, 0.719, 0.992, 1.300, 1.479, 1.692, 1.942, 2.227, 2.464, 2.500};
119  const static vector<double> eta_vals_25 = {0.692, 0.732, 0.768, 0.750, 0.726, 0.677, 0.667, 0.692, 0.710, 0.706, 0.679};
120 
121  const static vector<double> eta_edges_30 = {0.000, 0.053, 0.362, 0.719, 1.004, 1.277, 1.467, 1.681, 1.954, 2.239, 2.452, 2.500};
122  const static vector<double> eta_vals_30 = {0.724, 0.763, 0.804, 0.789, 0.762, 0.702, 0.690, 0.720, 0.731, 0.714, 0.681};
123 
124  const static vector<double> eta_edges_35 = {0.000, 0.044, 0.342, 0.711, 0.971, 1.280, 1.456, 1.683, 1.944, 2.218, 2.442, 2.500};
125  const static vector<double> eta_vals_35 = {0.736, 0.778, 0.824, 0.811, 0.784, 0.730, 0.718, 0.739, 0.743, 0.718, 0.678};
126 
127  const static vector<double> eta_edges_40 = {0.000, 0.047, 0.355, 0.699, 0.983, 1.268, 1.457, 1.671, 1.931, 2.204, 2.453, 2.500};
128  const static vector<double> eta_vals_40 = {0.741, 0.774, 0.823, 0.823, 0.802, 0.764, 0.756, 0.771, 0.771, 0.734, 0.684};
129 
130  const static vector<double> eta_edges_45 = {0.000, 0.056, 0.354, 0.711, 0.984, 1.280, 1.458, 1.684, 1.945, 2.207, 2.442, 2.500};
131  const static vector<double> eta_vals_45 = {0.758, 0.792, 0.841, 0.841, 0.823, 0.792, 0.786, 0.796, 0.794, 0.734, 0.663};
132 
133  const static vector<double> eta_edges_50 = {0.000, 0.059, 0.355, 0.699, 0.983, 1.268, 1.446, 1.682, 1.943, 2.216, 2.453, 2.500};
134  const static vector<double> eta_vals_50 = {0.771, 0.806, 0.855, 0.858, 0.843, 0.810, 0.800, 0.808, 0.802, 0.730, 0.653};
135 
136  const static vector<double> eta_edges_60 = {0.000, 0.050, 0.350, 0.707, 0.981, 1.278, 1.468, 1.694, 1.944, 2.242, 2.453, 2.500};
137  const static vector<double> eta_vals_60 = {0.773, 0.816, 0.866, 0.865, 0.853, 0.820, 0.812, 0.817, 0.804, 0.726, 0.645};
138 
139  const static vector<double> eta_edges_80 = {0.000, 0.051, 0.374, 0.720, 0.981, 1.279, 1.468, 1.707, 1.945, 2.207, 2.457, 2.500};
140  const static vector<double> eta_vals_80 = {0.819, 0.855, 0.899, 0.906, 0.900, 0.869, 0.865, 0.873, 0.869, 0.868, 0.859};
141 
142  const static vector<double> et_edges = { 10, 15, 20, 25, 30, 35, 40, 45, 50, 60, 80 };
143  const static vector< vector<double> > et_eta_edges = { eta_edges_10, eta_edges_15, eta_edges_20, eta_edges_25, eta_edges_30, eta_edges_35, eta_edges_40, eta_edges_45, eta_edges_50, eta_edges_60, eta_edges_80 };
144  const static vector< vector<double> > et_eta_vals = { eta_vals_10, eta_vals_15, eta_vals_20, eta_vals_25, eta_vals_30, eta_vals_35, eta_vals_40, eta_vals_45, eta_vals_50, eta_vals_60, eta_vals_80 };
145 
146  if (e.abseta() > 2.5 || e.Et() < 10*GeV) return 0.0;
147  const int i_et = binIndex(e.Et()/GeV, et_edges, true);
148  const int i_eta = binIndex(e.abseta(), et_eta_edges[i_et]);
149  return et_eta_vals[i_et][i_eta];
150  }
151 
154  inline double ELECTRON_IDEFF_ATLAS_RUN2_TIGHT(const Particle& e) {
156  }
157 
158 
159 
162  static const vector<double> edges_eta = {0., 2.5, 3.};
163  static const vector<double> edges_pt = {0., 0.1, 25.};
164  static const vector<double> e2s = {0.000, 0.015, 0.005,
165  0.005, 0.005, 0.005,
166  0.107, 0.107, 0.107};
167  static const vector<double> es = {0.00, 0.00, 0.05,
168  0.05, 0.05, 0.05,
169  2.08, 2.08, 2.08};
170  static const vector<double> cs = {0.00, 0.00, 0.25,
171  0.25, 0.25, 0.25,
172  0.00, 0.00, 0.00};
173 
174  const int i_eta = binIndex(e.abseta(), edges_eta, true);
175  const int i_pt = binIndex(e.pT()/GeV, edges_pt, true);
176  const int i = i_eta*edges_pt.size() + i_pt;
177 
178  // Calculate absolute resolution in GeV
179  const double c1 = sqr(e2s[i]), c2 = sqr(es[i]), c3 = sqr(cs[i]);
180  const double resolution = sqrt(c1*e.E2() + c2*e.E() + c3) * GeV;
181 
182  // normal_distribution<> d(e.E(), resolution);
183  // const double mass = e.mass2() > 0 ? e.mass() : 0; //< numerical carefulness...
184  // const double smeared_E = max(d(gen), mass); //< can't let the energy go below the mass!
185  // return Particle(e.pid(), FourMomentum::mkEtaPhiME(e.eta(), e.phi(), mass, smeared_E));
186  return Particle(e.pid(), P4_SMEAR_E_GAUSS(e, resolution));
187  }
188 
189 
193  return ELECTRON_SMEAR_ATLAS_RUN1(e);
194  }
195 
196 
198 
199 
200 
202  inline double ELECTRON_EFF_CMS_RUN1(const Particle& e) {
203  if (e.abseta() > 2.5) return 0;
204  if (e.pT() < 10*GeV) return 0;
205  return (e.abseta() < 1.5) ? 0.95 : 0.85;
206  }
207 
208 
211  inline double ELECTRON_EFF_CMS_RUN2(const Particle& e) {
212  return ELECTRON_EFF_CMS_RUN1(e);
213  }
214 
215 
223  // Calculate absolute resolution in GeV from functional form
224  double resolution = 0;
225  const double abseta = e.abseta();
226  if (e.pT() > 0.1*GeV && abseta < 2.5) { //< should be a given from efficiencies
227  if (abseta < 0.5) {
228  resolution = add_quad(0.06, 1.3e-3 * e.pT()/GeV) * GeV;
229  } else if (abseta < 1.5) {
230  resolution = add_quad(0.10, 1.7e-3 * e.pT()/GeV) * GeV;
231  } else { // still |eta| < 2.5
232  resolution = add_quad(0.25, 3.1e-3 * e.pT()/GeV) * GeV;
233  }
234  }
235 
236  // normal_distribution<> d(e.E(), resolution);
237  // const double mass = e.mass2() > 0 ? e.mass() : 0; //< numerical carefulness...
238  // const double smeared_E = max(d(gen), mass); //< can't let the energy go below the mass!
239  // return Particle(e.pid(), FourMomentum::mkEtaPhiME(e.eta(), e.phi(), mass, smeared_E));
240  return Particle(e.pid(), P4_SMEAR_E_GAUSS(e, resolution));
241  }
242 
246  return ELECTRON_SMEAR_CMS_RUN1(e);
247  }
248 
250 
251 
252 
254 
255 
259  inline double PHOTON_EFF_ATLAS_RUN1(const Particle& y) {
260  if (y.pT() < 10*GeV) return 0;
261  if (inRange(y.abseta(), 1.37, 1.52) || y.abseta() > 2.37) return 0;
262 
263  static const vector<double> edges_eta = {0., 0.6, 1.37, 1.52, 1.81, 2.37};
264  static const vector<double> edges_pt = {10., 15., 20., 25., 30., 35., 40., 45.,
265  50., 60., 80., 100., 125., 150., 175., 250.};
266  static const vector<double> effs = {0.53, 0.65, 0.73, 0.83, 0.86, 0.93, 0.94, 0.96,
267  0.97, 0.98, 0.98, 0.98, 0.98, 0.98, 0.98, 0.98,//
268  0.45, 0.57, 0.67, 0.74, 0.84, 0.87, 0.93, 0.94,
269  0.95, 0.96, 0.97, 0.98, 0.98, 0.99, 0.99, 0.99,//
270  0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00,
271  0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00,//
272  0.48, 0.56, 0.68, 0.76, 0.86, 0.90, 0.93, 0.95,
273  0.96, 0.97, 0.98, 0.99, 0.99, 1.00, 1.00, 1.00,//
274  0.50, 0.61, 0.74, 0.82, 0.88, 0.92, 0.94, 0.95,
275  0.96, 0.97, 0.98, 0.98, 0.98, 0.98, 0.99, 0.99};
276 
277  const int i_eta = binIndex(y.abseta(), edges_eta);
278  const int i_pt = binIndex(y.pT()/GeV, edges_pt, true);
279  const int i = i_eta*edges_pt.size() + i_pt;
280  const double eff = effs[i];
281  return eff;
282  }
283 
287  inline double PHOTON_EFF_ATLAS_RUN2(const Particle& y) {
288  if (y.pT() < 10*GeV) return 0;
289  if (inRange(y.abseta(), 1.37, 1.52) || y.abseta() > 2.37) return 0;
290 
291  static const vector<double> edges_eta = {0., 0.6, 1.37, 1.52, 1.81, 2.37};
292  static const vector<double> edges_pt = {10., 15., 20., 25., 30., 35., 40., 45.,
293  50., 60., 80., 100., 125., 150., 175., 250.};
294  static const vector<double> effs = {0.55, 0.70, 0.85, 0.89, 0.93, 0.95, 0.96, 0.96,
295  0.97, 0.97, 0.98, 0.97, 0.97, 0.97, 0.97, 0.97,//
296  0.47, 0.66, 0.79, 0.86, 0.89, 0.94, 0.96, 0.97,
297  0.97, 0.98, 0.97, 0.98, 0.98, 0.98, 0.98, 0.98,//
298  0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00,
299  0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00,//
300  0.54, 0.71, 0.84, 0.88, 0.92, 0.93, 0.94, 0.95,
301  0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96,//
302  0.61, 0.74, 0.83, 0.88, 0.91, 0.94, 0.95, 0.96,
303  0.97, 0.98, 0.98, 0.98, 0.98, 0.98, 0.98, 0.98};
304 
305  const int i_eta = binIndex(y.abseta(), edges_eta);
306  const int i_pt = binIndex(y.pT()/GeV, edges_pt, true);
307  const int i = i_eta*edges_pt.size() + i_pt;
308  const double eff = effs[i];
309  return eff;
310  }
311 
314  inline double PHOTON_EFF_CMS_RUN1(const Particle& y) {
315  if (y.pT() < 10*GeV || y.abseta() > 2.5) return 0;
316  return (y.abseta() < 1.5) ? 0.95 : 0.85;
317  }
318 
321  inline double PHOTON_EFF_CMS_RUN2(const Particle& y) {
322  return PHOTON_EFF_CMS_RUN1(y);
323  }
324 
325 
327  inline Particle PHOTON_SMEAR_ATLAS_RUN1(const Particle& y) { return y; }
328  inline Particle PHOTON_SMEAR_ATLAS_RUN2(const Particle& y) { return y; }
329  inline Particle PHOTON_SMEAR_CMS_RUN1(const Particle& y) { return y; }
330  inline Particle PHOTON_SMEAR_CMS_RUN2(const Particle& y) { return y; }
331 
333 
334 
335 
337 
338 
340  inline double MUON_EFF_ATLAS_RUN1(const Particle& m) {
341  if (m.abseta() > 2.7) return 0;
342  if (m.pT() < 10*GeV) return 0;
343  return (m.abseta() < 1.5) ? 0.95 : 0.85;
344  }
345 
350  inline double MUON_EFF_ATLAS_RUN2(const Particle& m) {
351  if (m.abseta() > 2.7) return 0;
352  static const vector<double> edges_pt = {0., 3.5, 4., 5., 6., 7., 8., 10.};
353  static const vector<double> effs = {0.00, 0.76, 0.94, 0.97, 0.98, 0.98, 0.98, 0.99};
354  const int i_pt = binIndex(m.pT()/GeV, edges_pt, true);
355  const double eff = effs[i_pt];
356  return eff;
357  }
358 
359 
362  static const vector<double> edges_eta = {0, 1.5, 2.5};
363  static const vector<double> edges_pt = {0, 0.1, 1.0, 10., 200.};
364  static const vector<double> res = {0., 0.03, 0.02, 0.03, 0.05,
365  0., 0.04, 0.03, 0.04, 0.05};
366 
367  const int i_eta = binIndex(m.abseta(), edges_eta);
368  const int i_pt = binIndex(m.pT()/GeV, edges_pt, true);
369  const int i = i_eta*edges_pt.size() + i_pt;
370 
371  const double resolution = res[i];
372 
373  // Smear by a Gaussian centered on the current pT, with width given by the resolution
374  // normal_distribution<> d(m.pT(), resolution*m.pT());
375  // const double smeared_pt = max(d(gen), 0.);
376  // const double mass = m.mass2() > 0 ? m.mass() : 0; //< numerical carefulness...
377  // return Particle(m.pid(), FourMomentum::mkEtaPhiMPt(m.eta(), m.phi(), mass, smeared_pt));
378  return Particle(m.pid(), P4_SMEAR_PT_GAUSS(m, resolution*m.pT()));
379  }
380 
384  return MUON_SMEAR_ATLAS_RUN1(m);
385  }
386 
387 
388 
389 
391  inline double MUON_EFF_CMS_RUN1(const Particle& m) {
392  if (m.abseta() > 2.4) return 0;
393  if (m.pT() < 10*GeV) return 0;
394  return 0.95 * (m.abseta() < 1.5 ? 1 : exp(0.5 - 5e-4*m.pT()/GeV));
395  }
396 
399  inline double MUON_EFF_CMS_RUN2(const Particle& m) {
400  return MUON_EFF_CMS_RUN1(m);
401  }
402 
403 
406  // Calculate fractional resolution
407  // for pT > 0.1 GeV, mom resolution = |eta| < 0.5 -> sqrt(0.01^2 + pt^2 * 2.0e-4^2)
408  // |eta| < 1.5 -> sqrt(0.02^2 + pt^2 * 3.0e-4^2)
409  // |eta| < 2.5 -> sqrt(0.05^2 + pt^2 * 2.6e-4^2)
410  double resolution = 0;
411  const double abseta = m.abseta();
412  if (m.pT() > 0.1*GeV && abseta < 2.5) {
413  if (abseta < 0.5) {
414  resolution = add_quad(0.01, 2.0e-4 * m.pT()/GeV);
415  } else if (abseta < 1.5) {
416  resolution = add_quad(0.02, 3.0e-4 * m.pT()/GeV);
417  } else { // still |eta| < 2.5... but isn't CMS' mu acceptance < 2.4?
418  resolution = add_quad(0.05, 2.6e-4 * m.pT()/GeV);
419  }
420  }
421 
422  // Smear by a Gaussian centered on the current pT, with width given by the resolution
423  // normal_distribution<> d(m.pT(), resolution*m.pT());
424  // const double smeared_pt = max(d(gen), 0.);
425  // const double mass = m.mass2() > 0 ? m.mass() : 0; //< numerical carefulness...
426  // return Particle(m.pid(), FourMomentum::mkEtaPhiMPt(m.eta(), m.phi(), mass, smeared_pt));
427  return Particle(m.pid(), P4_SMEAR_PT_GAUSS(m, resolution*m.pT()));
428  }
429 
433  return MUON_SMEAR_CMS_RUN1(m);
434  }
435 
437 
438 
439 
441 
442 
450  inline double TAU_EFF_ATLAS_RUN1(const Particle& t) {
451  if (t.abseta() > 2.5) return 0; //< hmm... mostly
452  double pThadvis = 0;
453  Particles chargedhadrons;
454  for (const Particle& p : t.children()) {
455  if (p.isHadron()) {
456  pThadvis += p.pT(); //< right definition? Paper is unclear
457  if (p.charge3() != 0 && p.abseta() < 2.5 && p.pT() > 1*GeV) chargedhadrons += p;
458  }
459  }
460  if (chargedhadrons.empty()) return 0; //< leptonic tau
461  if (pThadvis < 20*GeV) return 0; //< below threshold
462  if (pThadvis < 40*GeV) {
463  if (chargedhadrons.size() == 1) return (t.abspid() == PID::TAU) ? 0.56 : 1/20.;
464  if (chargedhadrons.size() == 3) return (t.abspid() == PID::TAU) ? 0.38 : 1/100.;
465  } else {
466  if (chargedhadrons.size() == 1) return (t.abspid() == PID::TAU) ? 0.56 : 1/25.;
467  if (chargedhadrons.size() == 3) return (t.abspid() == PID::TAU) ? 0.38 : 1/400.;
468  }
469  return 0;
470  }
471 
472 
478  inline double TAU_EFF_ATLAS_RUN2(const Particle& t) {
479  if (t.abseta() > 2.5) return 0; //< hmm... mostly
480  double pThadvis = 0;
481  Particles chargedhadrons;
482  for (const Particle& p : t.children()) {
483  if (p.isHadron()) {
484  pThadvis += p.pT(); //< right definition? Paper is unclear
485  if (p.charge3() != 0 && p.abseta() < 2.5 && p.pT() > 1*GeV) chargedhadrons += p;
486  }
487  }
488  if (chargedhadrons.empty()) return 0; //< leptonic tau
489  if (pThadvis < 20*GeV) return 0; //< below threshold
490  if (chargedhadrons.size() == 1) return (t.abspid() == PID::TAU) ? 0.55 : 1/50.;
491  if (chargedhadrons.size() == 3) return (t.abspid() == PID::TAU) ? 0.40 : 1/110.;
492  return 0;
493  }
494 
495 
499  // Const fractional resolution for now
500  static const double resolution = 0.03;
501 
502  // Smear by a Gaussian centered on 1 with width given by the (fractional) resolution
504  const double fsmear = max(randnorm(1., resolution), 0.);
505  const double mass = t.mass2() > 0 ? t.mass() : 0; //< numerical carefulness...
506  return Particle(t.pid(), FourMomentum::mkXYZM(t.px()*fsmear, t.py()*fsmear, t.pz()*fsmear, mass));
507  }
508 
509 
513  return TAU_SMEAR_ATLAS_RUN1(t);
514  }
515 
516 
520  inline double TAU_EFF_CMS_RUN2(const Particle& t) {
521  return (t.abspid() == PID::TAU) ? 0.6 : 0;
522  }
523 
527  inline double TAU_EFF_CMS_RUN1(const Particle& t) {
528  return TAU_EFF_CMS_RUN2(t);
529  }
530 
531 
535  return TAU_SMEAR_ATLAS_RUN1(t);
536  }
537 
538 
542  return TAU_SMEAR_CMS_RUN1(t);
543  }
544 
546 
547 
548 
550 
551 
553  inline double JET_BTAG_ATLAS_RUN1(const Jet& j) {
555  if (j.abseta() > 2.5) return 0;
556  const auto ftagsel = [&](const Particle& p){ return p.pT() > 5*GeV && deltaR(p,j) < 0.3; };
557  if (j.bTagged(ftagsel)) return 0.80*tanh(0.003*j.pT()/GeV)*(30/(1+0.0860*j.pT()/GeV));
558  if (j.cTagged(ftagsel)) return 0.20*tanh(0.020*j.pT()/GeV)*( 1/(1+0.0034*j.pT()/GeV));
559  return 0.002 + 7.3e-6*j.pT()/GeV;
560  }
562  inline double JET_BTAG_ATLAS_RUN2_MV2C20(const Jet& j) {
563  if (j.abseta() > 2.5) return 0;
564  if (j.bTagged(Cuts::pT > 5*GeV)) return 0.77;
565  if (j.cTagged(Cuts::pT > 5*GeV)) return 1/4.5;
566  return 1/140.;
567  }
569  inline double JET_BTAG_ATLAS_RUN2_MV2C10(const Jet& j) {
570  if (j.abseta() > 2.5) return 0;
571  if (j.bTagged(Cuts::pT > 5*GeV)) return 0.77;
572  if (j.cTagged(Cuts::pT > 5*GeV)) return 1/6.0;
573  return 1/134.;
574  }
575 
576 
578  inline Jet JET_SMEAR_ATLAS_RUN1(const Jet& j) {
579  // Jet energy resolution lookup
580  // Implemented by Matthias Danninger for GAMBIT, based roughly on
581  // https://atlas.web.cern.ch/Atlas/GROUPS/PHYSICS/CONFNOTES/ATLAS-CONF-2015-017/
582  // Parameterisation can be still improved, but eta dependence is minimal
584  static const vector<double> binedges_pt = {0., 50., 70., 100., 150., 200., 1000., 10000.};
585  static const vector<double> jer = {0.145, 0.115, 0.095, 0.075, 0.07, 0.05, 0.04, 0.04}; //< note overflow value
586  const int ipt = binIndex(j.pT()/GeV, binedges_pt, true);
587  if (ipt < 0) return j;
588  const double resolution = jer.at(ipt);
589 
590  // Smear by a Gaussian centered on 1 with width given by the (fractional) resolution
592  const double fsmear = max(randnorm(1., resolution), 0.);
593  const double mass = j.mass2() > 0 ? j.mass() : 0; //< numerical carefulness...
594  Jet rtn(FourMomentum::mkXYZM(j.px()*fsmear, j.py()*fsmear, j.pz()*fsmear, mass));
595  //if (deltaPhi(j, rtn) > 0.01) cout << "jdphi: " << deltaPhi(j, rtn) << endl;
596  return rtn;
597  }
598 
601  inline Jet JET_SMEAR_ATLAS_RUN2(const Jet& j) {
602  return JET_SMEAR_ATLAS_RUN1(j);
603  }
604 
607  inline Jet JET_SMEAR_CMS_RUN2(const Jet& j) {
608  return JET_SMEAR_ATLAS_RUN1(j);
609  }
610 
612 
613 
615 
616 
617  inline Vector3 MET_SMEAR_IDENTITY(const Vector3& met, double) { return met; }
618 
622  inline Vector3 MET_SMEAR_ATLAS_RUN1(const Vector3& met, double set) {
623  Vector3 smeared_met = met;
624 
625  // Linearity offset (Fig 14)
626  if (met.mod() < 25*GeV) smeared_met *= 1.05;
627  else if (met.mod() < 40*GeV) smeared_met *= (1.05 - (0.04/15)*(met.mod()/GeV - 25)); //< linear decrease
628  else smeared_met *= 1.01;
629 
630  // Smear by a Gaussian with width given by the resolution(sumEt) ~ 0.45 sqrt(sumEt) GeV
631  const double resolution = 0.45 * sqrt(set/GeV) * GeV;
632  const double metsmear = max(randnorm(smeared_met.mod(), resolution), 0.);
633  smeared_met = metsmear * smeared_met.unit();
634 
635  return smeared_met;
636  }
637 
640  inline Vector3 MET_SMEAR_ATLAS_RUN2(const Vector3& met, double set) {
641  return MET_SMEAR_ATLAS_RUN1(met, set);
642  }
643 
646  inline Vector3 MET_SMEAR_CMS_RUN1(const Vector3& met, double set) {
647  Vector3 smeared_met = met;
648 
649  // Calculate parallel and perpendicular resolutions and combine in quadrature (?)
650  const double resolution_x = (1.1 + 0.6*sqrt(set/GeV)) * GeV;
651  const double resolution_y = (1.4 + 0.6*sqrt(set/GeV)) * GeV;
652  const double resolution = sqrt(sqr(resolution_x) + sqr(resolution_y));
653 
654  // Smear by a Gaussian with width given by the resolution
655  const double metsmear = max(randnorm(smeared_met.mod(), resolution), 0.);
656  smeared_met = metsmear * smeared_met.unit();
657 
658  return smeared_met;
659  }
660 
663  inline Vector3 MET_SMEAR_CMS_RUN2(const Vector3& met, double set) {
664  Vector3 smeared_met = met;
665 
666  // Calculate parallel and perpendicular resolutions and combine in quadrature (?)
667  const double resolution_para = ( 2.0 + 0.64*sqrt(set/GeV)) * GeV;
668  const double resolution_perp = (-1.5 + 0.68*sqrt(set/GeV)) * GeV;
669  const double resolution = sqrt(sqr(resolution_para) + sqr(resolution_perp));
670 
671  // Smear by a Gaussian with width given by the resolution
672  const double metsmear = max(randnorm(smeared_met.mod(), resolution), 0.);
673  smeared_met = metsmear * smeared_met.unit();
674 
675  return smeared_met;
676  }
677 
679 
680 
682 
683 
685  inline double TRK_EFF_ATLAS_RUN1(const Particle& p) {
686  if (p.charge3() == 0) return 0;
687  if (p.abseta() > 2.5) return 0;
688  if (p.pT() < 0.1*GeV) return 0;
689 
690  if (p.abspid() == PID::ELECTRON) {
691  if (p.abseta() < 1.5) {
692  if (p.pT() < 1*GeV) return 0.73;
693  if (p.pT() < 100*GeV) return 0.95;
694  return 0.99;
695  } else {
696  if (p.pT() < 1*GeV) return 0.50;
697  if (p.pT() < 100*GeV) return 0.83;
698  else return 0.90;
699  }
700  } else if (p.abspid() == PID::MUON) {
701  if (p.abseta() < 1.5) {
702  return (p.pT() < 1*GeV) ? 0.75 : 0.99;
703  } else {
704  return (p.pT() < 1*GeV) ? 0.70 : 0.98;
705  }
706  } else { // charged hadrons
707  if (p.abseta() < 1.5) {
708  return (p.pT() < 1*GeV) ? 0.70 : 0.95;
709  } else {
710  return (p.pT() < 1*GeV) ? 0.60 : 0.85;
711  }
712  }
713  }
714 
717  inline double TRK_EFF_ATLAS_RUN2(const Particle& p) {
718  return TRK_EFF_ATLAS_RUN1(p);
719  }
720 
721 
723  inline double TRK_EFF_CMS_RUN1(const Particle& p) {
724  if (p.charge3() == 0) return 0;
725  if (p.abseta() > 2.5) return 0;
726  if (p.pT() < 0.1*GeV) return 0;
727 
728  if (p.abspid() == PID::ELECTRON) {
729  if (p.abseta() < 1.5) {
730  if (p.pT() < 1*GeV) return 0.73;
731  if (p.pT() < 100*GeV) return 0.95;
732  return 0.99;
733  } else {
734  if (p.pT() < 1*GeV) return 0.50;
735  if (p.pT() < 100*GeV) return 0.83;
736  else return 0.90;
737  }
738  } else if (p.abspid() == PID::MUON) {
739  if (p.abseta() < 1.5) {
740  return (p.pT() < 1*GeV) ? 0.75 : 0.99;
741  } else {
742  return (p.pT() < 1*GeV) ? 0.70 : 0.98;
743  }
744  } else { // charged hadrons
745  if (p.abseta() < 1.5) {
746  return (p.pT() < 1*GeV) ? 0.70 : 0.95;
747  } else {
748  return (p.pT() < 1*GeV) ? 0.60 : 0.85;
749  }
750  }
751  }
752 
755  inline double TRK_EFF_CMS_RUN2(const Particle& p) {
756  return TRK_EFF_CMS_RUN1(p);
757  }
758 
760 
761 
762 }
763 
764 #endif
Definition: ALICE_2010_I880049.cc:13
FourMomentum P4_SMEAR_E_GAUSS(const FourMomentum &p, double resolution)
Definition: MomentumSmearingFunctions.hh:46
double JET_BTAG_ATLAS_RUN2_MV2C10(const Jet &j)
Return the ATLAS Run 2 MC2c10 jet flavour tagging efficiency for the given Jet.
Definition: SmearingFunctions.hh:569
double MUON_EFF_CMS_RUN2(const Particle &m)
Definition: SmearingFunctions.hh:399
std::enable_if< std::is_arithmetic< NUM1 >::value &&std::is_arithmetic< NUM2 >::value, int >::type binIndex(NUM1 val, std::initializer_list< NUM2 > binedges, bool allow_overflow=false)
Return the bin index of the given value, val, given a vector of bin edges.
Definition: MathUtils.hh:356
Jet JET_SMEAR_CMS_RUN2(const Jet &j)
Definition: SmearingFunctions.hh:607
double pz() const
z component of momentum.
Definition: ParticleBase.hh:124
bool bTagged(const Cut &c=Cuts::open()) const
Does this jet have at least one b-tag (that passes an optional Cut)?
Definition: Jet.hh:110
double TAU_EFF_CMS_RUN2(const Particle &t)
Definition: SmearingFunctions.hh:520
double ELECTRON_IDEFF_ATLAS_RUN1_MEDIUM(const Particle &e)
ATLAS Run 1 &#39;medium&#39; electron identification/selection efficiency.
Definition: SmearingFunctions.hh:54
double TRK_EFF_CMS_RUN1(const Particle &p)
CMS Run 1 tracking efficiency.
Definition: SmearingFunctions.hh:723
Vector3 MET_SMEAR_CMS_RUN1(const Vector3 &met, double set)
Definition: SmearingFunctions.hh:646
double abseta() const
Get the directly (alias).
Definition: ParticleBase.hh:91
Vector3 MET_SMEAR_ATLAS_RUN2(const Vector3 &met, double set)
Definition: SmearingFunctions.hh:640
Particle TAU_SMEAR_CMS_RUN1(const Particle &t)
Definition: SmearingFunctions.hh:534
Particle PHOTON_SMEAR_ATLAS_RUN1(const Particle &y)
Definition: SmearingFunctions.hh:327
std::enable_if< std::is_arithmetic< NUM >::value, NUM >::type add_quad(NUM a, NUM b)
Named number-type addition in quadrature operation.
Definition: MathUtils.hh:201
double TRK_EFF_CMS_RUN2(const Particle &p)
Definition: SmearingFunctions.hh:755
PdgId pid() const
This Particle&#39;s PDG ID code.
Definition: Particle.hh:135
Vector3 MET_SMEAR_CMS_RUN2(const Vector3 &met, double set)
Definition: SmearingFunctions.hh:663
double E2() const
Get the energy-squared.
Definition: ParticleBase.hh:56
Particles children(const Cut &c=Cuts::OPEN) const
Get a list of the direct descendants from the current particle (with optional selection Cut) ...
Definition: Particle.cc:104
Particle representation, either from a HepMC::GenEvent or reconstructed.
Definition: Particle.hh:18
double mass2() const
Get the mass**2 directly.
Definition: ParticleBase.hh:82
double mass() const
Get the mass directly.
Definition: ParticleBase.hh:80
double PHOTON_EFF_CMS_RUN2(const Particle &y)
Definition: SmearingFunctions.hh:321
double E() const
Get the energy directly.
Definition: ParticleBase.hh:51
Particle ELECTRON_SMEAR_CMS_RUN2(const Particle &e)
Definition: SmearingFunctions.hh:245
Vector3 MET_SMEAR_ATLAS_RUN1(const Vector3 &met, double set)
ATLAS Run 1 ETmiss smearing.
Definition: SmearingFunctions.hh:622
double PHOTON_EFF_CMS_RUN1(const Particle &y)
Definition: SmearingFunctions.hh:314
double JET_BTAG_ATLAS_RUN1(const Jet &j)
Return the ATLAS Run 1 jet flavour tagging efficiency for the given Jet.
Definition: SmearingFunctions.hh:553
double Et() const
Get the directly.
Definition: ParticleBase.hh:75
double ELECTRON_IDEFF_ATLAS_RUN2_TIGHT(const Particle &e)
ATLAS Run 2 &#39;tight&#39; electron identification/selection efficiency.
Definition: SmearingFunctions.hh:154
double ELECTRON_EFF_CMS_RUN1(const Particle &e)
CMS Run 1 electron reconstruction efficiency.
Definition: SmearingFunctions.hh:202
Particle MUON_SMEAR_ATLAS_RUN2(const Particle &m)
Definition: SmearingFunctions.hh:383
Particle MUON_SMEAR_ATLAS_RUN1(const Particle &m)
ATLAS Run 1 muon reco smearing.
Definition: SmearingFunctions.hh:361
Particle ELECTRON_SMEAR_ATLAS_RUN2(const Particle &e)
Definition: SmearingFunctions.hh:192
double ELECTRON_EFF_ATLAS_RUN1(const Particle &e)
Definition: SmearingFunctions.hh:18
double TAU_EFF_ATLAS_RUN1(const Particle &t)
ATLAS Run 1 8 TeV tau efficiencies (medium working point)
Definition: SmearingFunctions.hh:450
std::enable_if< std::is_arithmetic< N1 >::value &&std::is_arithmetic< N2 >::value &&std::is_arithmetic< N3 >::value, bool >::type inRange(N1 value, N2 low, N3 high, RangeBoundary lowbound=CLOSED, RangeBoundary highbound=OPEN)
Determine if value is in the range low to high, for floating point numbers.
Definition: MathUtils.hh:103
Particle MUON_SMEAR_CMS_RUN1(const Particle &m)
CMS Run 1 muon reco smearing.
Definition: SmearingFunctions.hh:405
double TRK_EFF_ATLAS_RUN2(const Particle &p)
Definition: SmearingFunctions.hh:717
Particle ELECTRON_SMEAR_CMS_RUN1(const Particle &e)
CMS electron energy smearing, preserving direction.
Definition: SmearingFunctions.hh:222
double pT() const
Get the directly (alias).
Definition: ParticleBase.hh:63
double deltaR(double rap1, double phi1, double rap2, double phi2)
Definition: MathUtils.hh:597
int charge3() const
Three times the charge of this Particle (i.e. integer multiple of smallest quark charge).
Definition: Particle.hh:155
Jet JET_SMEAR_ATLAS_RUN2(const Jet &j)
Definition: SmearingFunctions.hh:601
double TAU_EFF_ATLAS_RUN2(const Particle &t)
ATLAS Run 2 13 TeV tau efficiencies (medium working point)
Definition: SmearingFunctions.hh:478
double max(const vector< double > &in, double errval=DBL_NAN)
Find the maximum value in the vector.
Definition: Utils.hh:465
Particle TAU_SMEAR_ATLAS_RUN2(const Particle &t)
Definition: SmearingFunctions.hh:512
double ELECTRON_EFF_ATLAS_RUN2(const Particle &e)
Definition: SmearingFunctions.hh:26
Particle TAU_SMEAR_CMS_RUN2(const Particle &t)
Definition: SmearingFunctions.hh:541
Vector3 unit() const
Synonym for unitVec.
Definition: Vector3.hh:105
double ELECTRON_IDEFF_ATLAS_RUN1_TIGHT(const Particle &e)
ATLAS Run 1 &#39;tight&#39; electron identification/selection efficiency.
Definition: SmearingFunctions.hh:107
double MUON_EFF_CMS_RUN1(const Particle &m)
CMS Run 1 muon reco efficiency.
Definition: SmearingFunctions.hh:391
double ELECTRON_EFF_CMS_RUN2(const Particle &e)
Definition: SmearingFunctions.hh:211
double px() const
x component of momentum.
Definition: ParticleBase.hh:120
double ELECTRON_IDEFF_ATLAS_RUN2_MEDIUM(const Particle &e)
ATLAS Run 2 &#39;medium&#39; electron identification/selection efficiency.
Definition: SmearingFunctions.hh:101
double TAU_EFF_CMS_RUN1(const Particle &t)
Definition: SmearingFunctions.hh:527
double PHOTON_EFF_ATLAS_RUN1(const Particle &y)
ATLAS Run 2 photon reco efficiency.
Definition: SmearingFunctions.hh:259
Representation of a clustered jet of particles.
Definition: Jet.hh:18
PdgId abspid() const
Absolute value of the PDG ID code.
Definition: Particle.hh:137
double JET_BTAG_ATLAS_RUN2_MV2C20(const Jet &j)
Return the ATLAS Run 2 MC2c20 jet flavour tagging efficiency for the given Jet.
Definition: SmearingFunctions.hh:562
Three-dimensional specialisation of Vector.
Definition: Vector3.hh:26
Particle ELECTRON_SMEAR_ATLAS_RUN1(const Particle &e)
ATLAS Run 1 electron reco smearing.
Definition: SmearingFunctions.hh:161
double ELECTRON_IDEFF_ATLAS_RUN2_LOOSE(const Particle &e)
ATLAS Run 2 &#39;loose&#39; electron identification/selection efficiency.
Definition: SmearingFunctions.hh:35
double MUON_EFF_ATLAS_RUN1(const Particle &m)
ATLAS Run 1 muon reco efficiency.
Definition: SmearingFunctions.hh:340
Particle TAU_SMEAR_ATLAS_RUN1(const Particle &t)
Definition: SmearingFunctions.hh:498
double TRK_EFF_ATLAS_RUN1(const Particle &p)
ATLAS Run 1 tracking efficiency.
Definition: SmearingFunctions.hh:685
bool cTagged(const Cut &c=Cuts::open()) const
Does this jet have at least one c-tag (that passes an optional Cut)?
Definition: Jet.hh:123
double PHOTON_EFF_ATLAS_RUN2(const Particle &y)
ATLAS Run 2 photon reco efficiency.
Definition: SmearingFunctions.hh:287
static FourMomentum mkXYZM(double px, double py, double pz, double mass)
Make a vector from (px,py,pz) coordinates and the mass.
Definition: Vector4.hh:776
double min(const vector< double > &in, double errval=DBL_NAN)
Find the minimum value in the vector.
Definition: Utils.hh:459
std::enable_if< std::is_arithmetic< NUM >::value, NUM >::type sqr(NUM a)
Named number-type squaring operation.
Definition: MathUtils.hh:189
double mod() const
Calculate the modulus of a vector. .
Definition: VectorN.hh:95
double py() const
y component of momentum.
Definition: ParticleBase.hh:122
FourMomentum P4_SMEAR_PT_GAUSS(const FourMomentum &p, double resolution)
Smear a FourMomentum&#39;s transverse momentum using a Gaussian of absolute width resolution.
Definition: MomentumSmearingFunctions.hh:53
Particle MUON_SMEAR_CMS_RUN2(const Particle &m)
Definition: SmearingFunctions.hh:432
double MUON_EFF_ATLAS_RUN2(const Particle &m)
ATLAS Run 2 muon reco efficiency.
Definition: SmearingFunctions.hh:350
double randnorm(double loc, double scale)
Return a Gaussian/normal sampled random number with the given mean and width.
Definition: Random.cc:47
Jet JET_SMEAR_ATLAS_RUN1(const Jet &j)
ATLAS Run 1 jet smearing.
Definition: SmearingFunctions.hh:578