gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-252-gf9a3f78
a Global And Modular Bsm Inference Tool
CMSEfficiencies.hpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
19 
20 
21 #pragma once
22 
23 #include <cfloat>
24 
27 
28 #include "HEPUtils/MathUtils.h"
29 #include "HEPUtils/BinnedFn.h"
30 #include "HEPUtils/Event.h"
31 #include <iomanip>
32 #include <algorithm>
33 
34 
35 namespace Gambit
36 {
37 
38  namespace ColliderBit
39  {
40 
41 
43  namespace CMS
44  {
45 
48 
50  inline void applyElectronTrackingEff(std::vector<HEPUtils::Particle*>& electrons) {
51  static HEPUtils::BinnedFn2D<double> _elTrackEff2d({{0, 1.5, 2.5, DBL_MAX}}, //< |eta|
52  {{0, 0.1, 1.0, DBL_MAX}}, //< pT
53  {{0., 0.70, 0.95, // |eta| 0.1-1.5
54  0., 0.60, 0.85, // |eta| 1.5-2.5
55  0., 0., 0.}}); // |eta| > 2.5
56  filtereff_etapt(electrons, _elTrackEff2d);
57  }
58 
59 
63  inline void applyElectronEff(std::vector<HEPUtils::Particle*>& electrons) {
64  static HEPUtils::BinnedFn2D<double> _elEff2d({{0, 1.5, 2.5, DBL_MAX}}, //< |eta|
65  {{0, 10., DBL_MAX}}, //< pT
66  {{0., 0.95, // |eta| 0.1-1.5
67  0., 0.85, // |eta| 1.5-2.5
68  0., 0.}}); // |eta| > 2.5
69  filtereff_etapt(electrons, _elEff2d);
70  }
71 
72 
75  inline void applyMuonTrackEff(std::vector<HEPUtils::Particle*>& muons) {
76  static HEPUtils::BinnedFn2D<double> _muTrackEff2d({{0, 1.5, 2.5, DBL_MAX}}, //< |eta|
77  {{0, 0.1, 1.0, DBL_MAX}}, //< pT
78  {{0, 0.75, 0.99, // |eta| 0.1-1.5
79  0, 0.70, 0.98, // |eta| 1.5-2.5
80  0, 0, 0}}); // |eta| > 2.5
81  filtereff_etapt(muons, _muTrackEff2d);
82  }
83 
84 
86  inline void applyMuonEff(std::vector<HEPUtils::Particle*>& muons) {
87  if(muons.empty()) return;
88  auto keptMuonsEnd = std::remove_if(muons.begin(), muons.end(),
89  [](const HEPUtils::Particle* p) {
90  bool rm(p->abseta() > 2.4 || p->pT() < 10);
91  if (!rm)
92  {
93  const double eff = 0.95 * (p->pT() <= 1.0e3 ? 1 : exp(0.5 - 5e-4*p->pT()));
94  rm = !random_bool(eff);
95  }
96  return rm;
97  } );
98  muons.erase(keptMuonsEnd, muons.end());
99  }
100 
101 
104  inline void applyTauEfficiency(std::vector<HEPUtils::Particle*>& taus) {
105  filtereff(taus, 0.6);
106  }
107 
108 
113  inline void smearElectronEnergy(std::vector<HEPUtils::Particle*>& electrons) {
114 
115  // Now loop over the electrons and smear the 4-vectors
116  for (HEPUtils::Particle* e : electrons) {
117 
118  // Calculate resolution
119  // for pT > 0.1 GeV, E resolution = |eta| < 0.5 -> sqrt(0.06^2 + pt^2 * 1.3e-3^2)
120  // |eta| < 1.5 -> sqrt(0.10^2 + pt^2 * 1.7e-3^2)
121  // |eta| < 2.5 -> sqrt(0.25^2 + pt^2 * 3.1e-3^2)
122  double resolution = 0;
123  const double abseta = e->abseta();
124  if (e->pT() > 0.1 && abseta < 2.5) {
125  if (abseta < 0.5) {
126  resolution = HEPUtils::add_quad(0.06, 1.3e-3 * e->pT());
127  } else if (abseta < 1.5) {
128  resolution = HEPUtils::add_quad(0.10, 1.7e-3 * e->pT());
129  } else { // still |eta| < 2.5
130  resolution = HEPUtils::add_quad(0.25, 3.1e-3 * e->pT());
131  }
132  }
133 
134  // Smear by a Gaussian centered on the current energy, with width given by the resolution
135  if (resolution > 0) {
136  std::normal_distribution<> d(e->E(), resolution);
137  double smeared_E = d(Random::rng());
138  if (smeared_E < 0) smeared_E = 0;
139  // double smeared_pt = smeared_E/cosh(e->eta()); ///< @todo Should be cosh(|eta|)?
140  // std::cout << "BEFORE eta " << electron->eta() << std::std::endl;
141  e->set_mom(HEPUtils::P4::mkEtaPhiME(e->eta(), e->phi(), e->mass(), smeared_E));
142  // std::cout << "AFTER eta " << electron->eta() << std::std::endl;
143  }
144  }
145  }
146 
147 
152  inline void smearMuonMomentum(std::vector<HEPUtils::Particle*>& muons) {
153 
154  // Now loop over the muons and smear the 4-vectors
155  for (HEPUtils::Particle* p : muons) {
156 
157  // Calculate resolution
158  // for pT > 0.1 GeV, mom resolution = |eta| < 0.5 -> sqrt(0.01^2 + pt^2 * 2.0e-4^2)
159  // |eta| < 1.5 -> sqrt(0.02^2 + pt^2 * 3.0e-4^2)
160  // |eta| < 2.5 -> sqrt(0.05^2 + pt^2 * 2.6e-4^2)
161  double resolution = 0;
162  const double abseta = p->abseta();
163  if (p->pT() > 0.1 && abseta < 2.5) {
164  if (abseta < 0.5) {
165  resolution = HEPUtils::add_quad(0.01, 2.0e-4 * p->pT());
166  } else if (abseta < 1.5) {
167  resolution = HEPUtils::add_quad(0.02, 3.0e-4 * p->pT());
168  } else { // still |eta| < 2.5... but isn't CMS' mu acceptance < 2.4?
169  resolution = HEPUtils::add_quad(0.05, 2.6e-4 * p->pT());
170  }
171  }
172 
173  // Smear by a Gaussian centered on the current pT, with width given by the resolution
174  std::normal_distribution<> d(p->pT(), resolution*p->pT());
175  double smeared_pt = d(Random::rng());
176  if (smeared_pt < 0) smeared_pt = 0;
177  // const double smeared_E = smeared_pt*cosh(mu->eta()); ///< @todo Should be cosh(|eta|)?
178  // std::cout << "Muon pt " << mu_pt << " smeared " << smeared_pt << std::endl;
179  p->set_mom(HEPUtils::P4::mkEtaPhiMPt(p->eta(), p->phi(), p->mass(), smeared_pt));
180  }
181  }
182 
183 
190  inline void smearJets(std::vector<HEPUtils::Jet*>& jets) {
191 
192  // Const resolution for now
197  //const double resolution = 0.03;
198 
199  // Matthias jet smearing implemented roughly from https://arxiv.org/pdf/1607.03663.pdf
200  // Parameterisation can be still improved as functional form is given
201  // Pileup of <mu>=25 is taken, as JER depends strongly on mu
202  // CMS does not include information about JER at eta>1.3
203  const std::vector<double> binedges_eta = {0,10.};
204  const std::vector<double> binedges_pt = {0,20,30,40,50.,70.,100.,150.,200.,1000.,10000.};
205  const std::vector<double> JetsJER = {0.3,0.2,0.16,0.145,0.12,0.1,0.09,0.08,0.06,0.05};
206  static HEPUtils::BinnedFn2D<double> _resJets2D(binedges_eta,binedges_pt,JetsJER);
207 
208  // Now loop over the jets and smear the 4-vectors
209  for (HEPUtils::Jet* jet : jets) {
210  const double resolution = _resJets2D.get_at(jet->abseta(), jet->pT());
211  std::normal_distribution<> d(1., resolution);
212  // Smear by a Gaussian centered on 1 with width given by the (fractional) resolution
213  double smear_factor = d(Random::rng());
214  jet->set_mom(HEPUtils::P4::mkXYZM(jet->mom().px()*smear_factor, jet->mom().py()*smear_factor, jet->mom().pz()*smear_factor, jet->mass()));
215  }
216  }
217 
218 
225  inline void smearTaus(std::vector<HEPUtils::Particle*>& taus) {
226 
227  // Const resolution for now
228  const double resolution = 0.03;
229 
230  // Now loop over the jets and smear the 4-vectors
231  std::normal_distribution<> d(1., resolution);
232  for (HEPUtils::Particle* p : taus) {
233  // Smear by a Gaussian centered on 1 with width given by the (fractional) resolution
234  double smear_factor = d(Random::rng());
236  p->set_mom(HEPUtils::P4::mkXYZM(p->mom().px()*smear_factor, p->mom().py()*smear_factor, p->mom().pz()*smear_factor, p->mass()));
237  }
238  }
239 
242  inline void applyCSVv2MediumBtagEff(std::vector<const HEPUtils::Jet*>& bjets) {
243  if (bjets.empty()) return;
244 
245  const static std::vector<double> binedges_et = {25., 40., 60., 80., 100., 150., 200., 250., 300., 400., 500.,DBL_MAX };
246  const static std::vector<double> bineffs_et = {0.58, 0.61, 0.63, 0.64, 0.65, 0.62,0.6, 0.58, 0.56, 0.52, 0.48};
247  const static HEPUtils::BinnedFn1D<double> _eff_et(binedges_et, bineffs_et);
248 
249  auto keptBjetsEnd = std::remove_if(bjets.begin(), bjets.end(),
250  [](const HEPUtils::Jet* bjet) {
251  const double bjet_pt = bjet->pT();
252  const double bjet_aeta = bjet->abseta();
253  if (bjet_aeta > 2.4 || bjet_pt < 25) return true;
254  const double eff = _eff_et.get_at(bjet_pt);
255  return random_bool(1-eff);
256  } );
257  bjets.erase(keptBjetsEnd, bjets.end());
258  }
259 
260  inline void applyCSVv2MediumBtagEff(std::vector<HEPUtils::Jet*>& bjets) {
261  applyCSVv2MediumBtagEff(reinterpret_cast<std::vector<const HEPUtils::Jet*>&>(bjets));
262  }
263 
266  inline void applyCSVv2LooseBtagEff(std::vector<const HEPUtils::Jet*>& bjets) {
267  if (bjets.empty()) return;
268 
269  const static std::vector<double> binedges_et = {25., 40., 60., 80., 100., 150., 200., 250., 300., 400., 500.,DBL_MAX };
270  const static std::vector<double> bineffs_et = {0.78, 0.80, 0.82, 0.83, 0.84, 0.825, 0.82, 0.81, 0.8, 0.795, 0.79};
271  const static HEPUtils::BinnedFn1D<double> _eff_et(binedges_et, bineffs_et);
272 
273  auto keptBjetsEnd = std::remove_if(bjets.begin(), bjets.end(),
274  [](const HEPUtils::Jet* bjet) {
275  const double bjet_pt = bjet->pT();
276  const double bjet_aeta = bjet->abseta();
277  if (bjet_aeta > 2.4 || bjet_pt < 25) return true;
278  const double eff = _eff_et.get_at(bjet_pt);
279  return random_bool(1-eff);
280  } );
281  bjets.erase(keptBjetsEnd, bjets.end());
282  }
283 
284  inline void applyCSVv2LooseBtagEff(std::vector<HEPUtils::Jet*>& bjets) {
285  applyCSVv2LooseBtagEff(reinterpret_cast<std::vector<const HEPUtils::Jet*>&>(bjets));
286  }
287 
288 
290  inline void applyBtagMisId(double mis_id_prob, std::vector<const HEPUtils::Jet*>& jets, std::vector<const HEPUtils::Jet*>& bjets) {
291  if (jets.empty()) return;
292  for (const HEPUtils::Jet* jet : jets) {
293  // Only apply misidentification rate for non-b-jets
294  if (!jet->btag() && random_bool(mis_id_prob)) bjets.push_back(jet);
295  }
296  }
297 
298  inline void applyBtagMisId(double mis_id_prob, std::vector<HEPUtils::Jet*>& jets, std::vector<HEPUtils::Jet*>& bjets) {
299  applyBtagMisId(mis_id_prob, reinterpret_cast<std::vector<const HEPUtils::Jet*>&>(jets), reinterpret_cast<std::vector<const HEPUtils::Jet*>&>(bjets));
300  }
301 
302 
305  inline void applyCSVv2LooseBtagMisId(std::vector<const HEPUtils::Jet*>& jets, std::vector<const HEPUtils::Jet*>& bjets) {
306  if (jets.empty()) return;
307  // For now we apply the (pT-averaged) light-flavour misidentification rate to all jets.
308  // Realistically, the rate should be higher for c-jets.
309  const static double mis_id_prob = 0.089;
310  applyBtagMisId(mis_id_prob, jets, bjets);
311  }
312 
313  inline void applyCSVv2LooseBtagMisId(std::vector<HEPUtils::Jet*>& jets, std::vector<HEPUtils::Jet*>& bjets) {
314  applyCSVv2LooseBtagMisId(reinterpret_cast<std::vector<const HEPUtils::Jet*>&>(jets), reinterpret_cast<std::vector<const HEPUtils::Jet*>&>(bjets));
315  }
316 
317 
319  inline void applyCSVv2LooseBtagEffAndMisId(std::vector<const HEPUtils::Jet*>& jets, std::vector<const HEPUtils::Jet*>& bjets) {
320  if (jets.empty() && bjets.empty()) return;
321  // Apply b-tag efficiency
322  applyCSVv2LooseBtagEff(bjets);
323  // Apply misidentification rate to the non-b-jets in the jets vector
324  applyCSVv2LooseBtagMisId(jets, bjets);
325  }
326 
327  inline void applyCSVv2LooseBtagEffAndMisId(std::vector<HEPUtils::Jet*>& jets, std::vector<HEPUtils::Jet*>& bjets) {
328  applyCSVv2LooseBtagEffAndMisId(reinterpret_cast<std::vector<const HEPUtils::Jet*>&>(jets), reinterpret_cast<std::vector<const HEPUtils::Jet*>&>(bjets));
329  }
330 
331 
334  inline void applyCSVv2MediumBtagMisId(std::vector<const HEPUtils::Jet*>& jets, std::vector<const HEPUtils::Jet*>& bjets) {
335  if (jets.empty()) return;
336  // For now we apply the (pT-averaged) light-flavour misidentification rate to all jets.
337  // Realistically, the rate should be higher for c-jets.
338  const static double mis_id_prob = 0.009;
339  applyBtagMisId(mis_id_prob, jets, bjets);
340  }
341 
342  inline void applyCSVv2MediumBtagMisId(std::vector<HEPUtils::Jet*>& jets, std::vector<HEPUtils::Jet*>& bjets) {
343  applyCSVv2MediumBtagMisId(reinterpret_cast<std::vector<const HEPUtils::Jet*>&>(jets), reinterpret_cast<std::vector<const HEPUtils::Jet*>&>(bjets));
344  }
345 
346 
348  inline void applyCSVv2MediumBtagEffAndMisId(std::vector<const HEPUtils::Jet*>& jets, std::vector<const HEPUtils::Jet*>& bjets) {
349  if (jets.empty() && bjets.empty()) return;
350  // Apply b-tag efficiency
352  // Apply misidentification rate to the non-b-jets in the jets vector
353  applyCSVv2MediumBtagMisId(jets, bjets);
354  }
355 
356  inline void applyCSVv2MediumBtagEffAndMisId(std::vector<HEPUtils::Jet*>& jets, std::vector<HEPUtils::Jet*>& bjets) {
357  applyCSVv2MediumBtagEffAndMisId(reinterpret_cast<std::vector<const HEPUtils::Jet*>&>(jets), reinterpret_cast<std::vector<const HEPUtils::Jet*>&>(bjets));
358  }
359 
361 
362  }
363  }
364 }
void applyElectronEff(std::vector< HEPUtils::Particle *> &electrons)
Randomly filter the supplied particle list by parameterised electron efficiency.
void filtereff(std::vector< HEPUtils::Particle *> &particles, double eff, bool do_delete=false)
Utility function for filtering a supplied particle vector by sampling wrt an efficiency scalar...
Definition: Utils.cpp:16
void applyCSVv2LooseBtagMisId(std::vector< const HEPUtils::Jet *> &jets, std::vector< const HEPUtils::Jet *> &bjets)
Apply b-tag misidentification rate for CSVv2 loose WP.
void applyCSVv2LooseBtagEffAndMisId(std::vector< const HEPUtils::Jet *> &jets, std::vector< const HEPUtils::Jet *> &bjets)
Apply both b-tag efficiency and misidentification rate for CSVv2 loose WP.
void smearTaus(std::vector< HEPUtils::Particle *> &taus)
Randomly smear the supplied hadronic taus&#39; momenta by parameterised resolutions.
void smearJets(std::vector< HEPUtils::Jet *> &jets)
Randomly smear the supplied jets&#39; momenta by parameterised resolutions.
void smearMuonMomentum(std::vector< HEPUtils::Particle *> &muons)
Randomly smear the supplied muons&#39; momenta by parameterised resolutions.
void smearElectronEnergy(std::vector< HEPUtils::Particle *> &electrons)
Randomly smear the supplied electrons&#39; momenta by parameterised resolutions.
void applyElectronTrackingEff(std::vector< HEPUtils::Particle *> &electrons)
Randomly filter the supplied particle list by parameterised electron tracking efficiency.
void applyMuonTrackEff(std::vector< HEPUtils::Particle *> &muons)
Randomly filter the supplied particle list by parameterised muon tracking efficiency.
void applyCSVv2MediumBtagEff(std::vector< const HEPUtils::Jet *> &bjets)
Apply efficiency function to CSVv2 medium WP b-tagged jets.
void applyCSVv2MediumBtagEffAndMisId(std::vector< const HEPUtils::Jet *> &jets, std::vector< const HEPUtils::Jet *> &bjets)
Apply both b-tag efficiency and misidentification rate for CSVv2 medium WP.
A threadsafe interface to the STL random number generators.
void filtereff_etapt(std::vector< HEPUtils::Particle *> &particles, const HEPUtils::BinnedFn2D< double > &eff_etapt, bool do_delete=false)
Utility function for filtering a supplied particle vector by sampling wrt a binned 2D efficiency map ...
Definition: Utils.cpp:56
void applyCSVv2MediumBtagMisId(std::vector< const HEPUtils::Jet *> &jets, std::vector< const HEPUtils::Jet *> &bjets)
Apply b-tag misidentification rate for CSVv2 medium WP.
void applyCSVv2LooseBtagEff(std::vector< const HEPUtils::Jet *> &bjets)
Apply efficiency function to CSVv2 loose WP b-tagged jets.
void applyBtagMisId(double mis_id_prob, std::vector< const HEPUtils::Jet *> &jets, std::vector< const HEPUtils::Jet *> &bjets)
Apply user-specified b-tag misidentification rate (flat)
static Utils::threadsafe_rng & rng()
Return a threadsafe wrapper for the chosen RNG engine (to be passed to e.g.
TODO: see if we can use this one:
Definition: Analysis.hpp:33
void applyMuonEff(std::vector< HEPUtils::Particle *> &muons)
Randomly filter the supplied particle list by parameterised muon efficiency.
void applyTauEfficiency(std::vector< HEPUtils::Particle *> &taus)
Randomly filter the supplied particle list by parameterised tau efficiency.
bool random_bool(double eff)
Return a random true/false at a success rate given by a number.
Definition: Utils.cpp:10