gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-2191-ga4742ac
a Global And Modular Bsm Inference Tool
Utils.hpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
18 
19 #pragma once
20 
21 #include <functional>
22 #include <memory>
23 #include <cfloat>
24 
25 #include "HEPUtils/MathUtils.h"
26 #include "HEPUtils/BinnedFn.h"
27 #include "HEPUtils/Event.h"
28 #include "HEPUtils/FastJet.h"
29 
30 namespace Gambit
31 {
32 
33  namespace ColliderBit
34  {
35 
36 
38  static const double GeV = 1, MeV = 1e-3, TeV = 1e3;
39 
40 
42  using HEPUtils::Event;
44  using HEPUtils::Particle;
46  using HEPUtils::Jet;
48  using HEPUtils::P4;
49 
51  using HEPUtils::add_quad;
52 
54  typedef std::vector<const HEPUtils::Particle*> ParticlePtrs;
55 
57  typedef std::vector<const HEPUtils::Jet*> JetPtrs;
58 
60 
61 
63  inline bool amIaJet(const HEPUtils::Jet *jet) { (void)jet; return true; }
64 
66  inline bool amIaBJet(const HEPUtils::Jet *jet) { return jet->btag(); }
67 
69  inline bool amIaJet(const HEPUtils::Particle *part) { (void)part; return false; }
70 
72  inline bool amIaBJet(const HEPUtils::Particle *part) { (void)part; return true; }
73 
75 
77 
78 
80  template <typename CONTAINER, typename RMFN>
81  inline void iremoveerase(CONTAINER& c, const RMFN& fn) {
82  auto newend = std::remove_if(c.begin(), c.end(), fn);
83  c.erase(newend, c.end());
84  }
85 
87  inline void ifilter_reject(ParticlePtrs& particles,
88  std::function<bool(const Particle*)> rejfn, bool do_delete=true) {
89  iremoveerase(particles, [&](const Particle* p) {
90  const bool rm = rejfn(p);
91  if (rm && do_delete) delete p;
92  return rm;
93  });
94  }
95 
97  inline void ifilter_select(ParticlePtrs& particles,
98  std::function<bool(const Particle*)> selfn, bool do_delete=true) {
99  ifilter_reject(particles, [&](const Particle* p) { return !selfn(p); }, do_delete);
100  }
101 
102 
105  inline ParticlePtrs filter_reject(const ParticlePtrs& particles,
106  std::function<bool(const Particle*)> rejfn, bool do_delete=true) {
107  ParticlePtrs rtn = particles;
108  ifilter_reject(rtn, rejfn, do_delete);
109  return rtn;
110  }
111 
113  inline ParticlePtrs filter_select(const ParticlePtrs& particles,
114  std::function<bool(const Particle*)> selfn, bool do_delete=true) {
115  return filter_reject(particles, [&](const Particle* p) { return !selfn(p); }, do_delete);
116  }
117 
119 
120 
121 
123 
124 
126  inline void ifilter_reject(JetPtrs& jets,
127  std::function<bool(const Jet*)> rejfn, bool do_delete=true) {
128  iremoveerase(jets, [&](const Jet* j) {
129  const bool rm = rejfn(j);
130  if (rm && do_delete) delete j;
131  return rm;
132  });
133  }
134 
136  inline void ifilter_select(JetPtrs& jets,
137  std::function<bool(const Jet*)> selfn, bool do_delete=true) {
138  ifilter_reject(jets, [&](const Jet* j) { return !selfn(j); }, do_delete);
139  }
140 
141 
144  inline JetPtrs filter_reject(const JetPtrs& jets,
145  std::function<bool(const Jet*)> rejfn, bool do_delete=true) {
146  JetPtrs rtn = jets;
147  ifilter_reject(rtn, rejfn, do_delete);
148  return rtn;
149  }
150 
152  inline JetPtrs filter_select(const JetPtrs& jets,
153  std::function<bool(const Jet*)> selfn, bool do_delete=true) {
154  return filter_reject(jets, [&](const Jet* j) { return !selfn(j); }, do_delete);
155  }
156 
158 
159 
160 
162 
163 
164 
166 
167 
169  bool random_bool(double eff);
170 
172  inline bool random_bool(const HEPUtils::BinnedFn1D<double>& effmap, double x) {
173  return random_bool( effmap.get_at(x) );
174  }
175 
177  inline bool random_bool(const HEPUtils::BinnedFn2D<double>& effmap, double x, double y) {
178  return random_bool( effmap.get_at(x, y) );
179  }
180 
182 
183 
185 
186 
188  void filtereff(std::vector<const HEPUtils::Particle*>& particles, double eff, bool do_delete=false);
189 
191  void filtereff(std::vector<const HEPUtils::Particle*>& particles, std::function<double(const HEPUtils::Particle*)> eff_fn, bool do_delete=false);
192 
194  void filtereff_pt(std::vector<const HEPUtils::Particle*>& particles, const HEPUtils::BinnedFn1D<double>& eff_pt, bool do_delete=false);
195 
197  void filtereff_etapt(std::vector<const HEPUtils::Particle*>& particles, const HEPUtils::BinnedFn2D<double>& eff_etapt, bool do_delete=false);
198 
200 
201 
203 
204 
207  inline bool has_tag(const HEPUtils::BinnedFn2D<double>& effmap, double eta, double pt) {
208  try {
209  return random_bool(effmap, fabs(eta), pt);
210  } catch (...) {
211  return false; // No tag if outside lookup range... be careful!
212  }
213  }
214 
216  inline std::map<const HEPUtils::Jet*,bool> generateBTagsMap(const std::vector<const HEPUtils::Jet*>& jets,
217  double bTagEff, double cMissTagEff, double otherMissTagEff,
218  double pTmin = 0., double absEtaMax = DBL_MAX)
219  {
220  std::map<const HEPUtils::Jet*,bool> bTagsMap;
221  for (const HEPUtils::Jet* j : jets)
222  {
223  bool genBTag = false;
224  if((j->pT() > pTmin) && (j->abseta() < absEtaMax))
225  {
226  if(j->btag())
227  {
228  if(random_bool(bTagEff)) { genBTag = true; }
229  }
230  else if(j->ctag())
231  {
232  if(random_bool(cMissTagEff)) { genBTag = true; }
233  }
234  else
235  {
236  if(random_bool(otherMissTagEff)) { genBTag = true; }
237  }
238  }
239  bTagsMap[j] = genBTag;
240  }
241  return bTagsMap;
242  }
243 
244 
245  template <typename NUM1, typename NUM2>
246  inline size_t binIndex(NUM1 val, const std::vector<NUM2>& binedges, bool allow_overflow=false) {
247  if (val < binedges.front()) return -1;
248  if (val >= binedges.back()) return allow_overflow ? int(binedges.size())-1 : -1;
249  return std::distance(binedges.begin(), --std::upper_bound(binedges.begin(), binedges.end(), val));
250  }
251 
252 
254  inline std::vector<double> mk_bin_values(const std::vector<double>& binEdgeValues) {
255  std::vector<double> results;
256  results.reserve(binEdgeValues.size()-1);
257  for (size_t i = 0; i < binEdgeValues.size()-1; ++i)
258  results.push_back( (binEdgeValues[i] + binEdgeValues[i+1])/2.0 );
259  return results;
260  }
262  inline std::vector<double> makeBinValues(const std::vector<double>& binEdgeValues) {
263  return mk_bin_values(binEdgeValues);
264  }
265 
266 
268  template <typename MOM>
269  inline std::vector<std::shared_ptr<HEPUtils::Jet>> get_jets(const std::vector<MOM*>& moms, double R,
270  double ptmin=0*GeV, FJNS::JetAlgorithm alg=FJNS::antikt_algorithm) {
271  // Make PseudoJets
272  std::vector<FJNS::PseudoJet> constituents;
273  for (const MOM* p : moms) constituents.push_back(HEPUtils::mk_pseudojet(*p));
274  // Run clustering
275  std::vector<FJNS::PseudoJet> jets = HEPUtils::get_jets(constituents, R, ptmin, alg);
276  // Make newly-allocated Jets
277  std::vector<std::shared_ptr<HEPUtils::Jet>> rtn;
278  for (const FJNS::PseudoJet& j : jets) rtn.push_back(std::make_shared<HEPUtils::Jet>(HEPUtils::mk_p4(j)));
279  return rtn;
280  }
281 
282 
284  inline bool object_in_cone(const HEPUtils::Event& e, const HEPUtils::P4& p4, double ptmin, double rmax, double rmin=0.05) {
285  for (const HEPUtils::Particle* p : e.visible_particles())
286  if (p->pT() > ptmin && HEPUtils::in_range(HEPUtils::deltaR_eta(p4, *p), rmin, rmax)) return true;
287  for (const HEPUtils::Jet* j : e.jets())
288  if (j->pT() > ptmin && HEPUtils::in_range(HEPUtils::deltaR_eta(p4, *j), rmin, rmax)) return true;
289  return false;
290  }
291 
292 
298  template<typename MOMPTRS1, typename MOMPTRS2>
299  void removeOverlap(MOMPTRS1& momstofilter, const MOMPTRS2& momstocmp, double deltaRMax, bool use_rapidity=false, double pTmax = DBL_MAX, double btageff = 0)
300  {
301  ifilter_reject(momstofilter, [&](const typename MOMPTRS1::value_type& mom1)
302  {
303  for (const typename MOMPTRS2::value_type& mom2 : momstocmp) {
304  const double dR = (use_rapidity) ? deltaR_rap(mom1->mom(), mom2->mom()) : deltaR_eta(mom1->mom(), mom2->mom());
305  if (dR < deltaRMax && mom1->pT() < pTmax && ( !amIaBJet(mom1) || !random_bool(btageff) ) ) return true;
306  }
307  return false;
308  }, false);
309  }
310 
317  template<typename MOMPTRS1, typename MOMPTRS2>
318  void removeOverlap(MOMPTRS1& momstofilter, const MOMPTRS2& momstocmp, double (*deltaRMax)(const double), bool use_rapidity=false, double pTmax = DBL_MAX, double btageff = 0)
319  {
320  ifilter_reject(momstofilter, [&](const typename MOMPTRS1::value_type& mom1)
321  {
322  for (const typename MOMPTRS2::value_type& mom2 : momstocmp) {
323  const double dR = (use_rapidity) ? deltaR_rap(mom1->mom(), mom2->mom()) : deltaR_eta(mom1->mom(), mom2->mom());
324  if (dR < deltaRMax(mom1->pT()) && mom1->pT() < pTmax && ( !amIaBJet(mom1) || !random_bool(btageff) ) ) return true;
325  }
326  return false;
327  }, false);
328  }
329 
334  template<typename MOMPTRS1>
335  void removeOverlapIfBjet(MOMPTRS1& momstofilter, std::vector<const HEPUtils::Jet*>& jets, double deltaRMax, bool use_rapidity=false, double pTmax = DBL_MAX)
336  {
337  ifilter_reject(momstofilter, [&](const typename MOMPTRS1::value_type& mom1)
338  {
339  for (const HEPUtils::Jet* jet : jets) {
340  const double dR = (use_rapidity) ? deltaR_rap(mom1->mom(), jet->mom()) : deltaR_eta(mom1->mom(), jet->mom());
341  if (dR < deltaRMax && mom1->pT() < pTmax && jet->btag() ) return true;
342  }
343  return false;
344  }, false);
345  }
346 
347 
349  template <typename CONTAINER, typename FN>
350  inline bool all_of(const CONTAINER& c, const FN& f) {
351  return std::all_of(std::begin(c), std::end(c), f);
352  }
353 
355  template <typename CONTAINER, typename FN>
356  inline bool any_of(const CONTAINER& c, const FN& f) {
357  return std::any_of(std::begin(c), std::end(c), f);
358  }
359 
361  template <typename CONTAINER, typename FN>
362  inline bool none_of(const CONTAINER& c, const FN& f) {
363  return std::none_of(std::begin(c), std::end(c), f);
364  }
365 
366 
368  std::vector<std::vector<const HEPUtils::Particle*>> getSFOSpairs(std::vector<const HEPUtils::Particle*> particles);
369 
371  std::vector<std::vector<const HEPUtils::Particle*>> getOSpairs(std::vector<const HEPUtils::Particle*> particles);
372 
374  std::vector<std::vector<const HEPUtils::Particle*>> getSSpairs(std::vector<const HEPUtils::Particle*> particles);
375 
376 
378 
379 
381  inline void sortBy(ParticlePtrs& particles, std::function<bool(const Particle*, const Particle*)> cmpfn) {
382  std::sort(particles.begin(), particles.end(), cmpfn);
383  }
384 
386  inline bool cmpParticlesByPt(const HEPUtils::Particle* lep1, const HEPUtils::Particle* lep2) { return lep1->pT() > lep2->pT(); }
387 
388  // Sort a particles list by decreasing pT
389  inline void sortByPt(ParticlePtrs& particles) { sortBy(particles, cmpParticlesByPt); }
390 
391 
393  inline void sortBy(JetPtrs& jets, std::function<bool(const Jet*, const Jet*)> cmpfn) {
394  std::sort(jets.begin(), jets.end(), cmpfn);
395  }
396 
398  inline bool cmpJetsByPt(const HEPUtils::Jet* jet1, const HEPUtils::Jet* jet2) { return jet1->pT() > jet2->pT(); }
399 
400  // Sort a jets list by decreasing pT
401  inline void sortByPt(JetPtrs& jets) { sortBy(jets, cmpJetsByPt); }
403 
404 
406 
407 
409  inline int countPt(const std::vector<const Particle*>& particles, double pTlim)
410  {
411  int sum = 0;
412  for (const Particle* p : particles)
413  {
414  if (p->pT() > pTlim) { sum++; }
415  }
416  return sum;
417  }
418 
420  inline int countPt(const std::vector<const Jet*>& jets, double pTlim)
421  {
422  int sum = 0;
423  for (const Jet* j : jets)
424  {
425  if (j->pT() > pTlim) { sum++; }
426  }
427  return sum;
428  }
429 
431 
432 
434 
435 
437  inline double scalarSumPt(const std::vector<const Particle*>& particles, double pTlim=0.)
438  {
439  double sum = 0.;
440  for (const Particle* p : particles)
441  {
442  if (p->pT() > pTlim) { sum += p->pT(); }
443  }
444  return sum;
445  }
446 
448  inline double scalarSumPt(const std::vector<const Jet*>& jets, double pTlim=0.)
449  {
450  double sum = 0.;
451  for (const Jet* j : jets)
452  {
453  if (j->pT() > pTlim) { sum += j->pT(); }
454  }
455  return sum;
456  }
457 
459 
460  }
461 
462 }
void filtereff_etapt(std::vector< const 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 removeOverlapIfBjet(MOMPTRS1 &momstofilter, std::vector< const HEPUtils::Jet *> &jets, double deltaRMax, bool use_rapidity=false, double pTmax=DBL_MAX)
Overlap removal for checking against b-jets – discard from first list if within deltaRMax of a b-jet...
Definition: Utils.hpp:335
bool has_tag(const HEPUtils::BinnedFn2D< double > &effmap, double eta, double pt)
Randomly get a tag result (can be anything) from a 2D |eta|-pT efficiency map.
Definition: Utils.hpp:207
bool object_in_cone(const HEPUtils::Event &e, const HEPUtils::P4 &p4, double ptmin, double rmax, double rmin=0.05)
Check if there&#39;s a physics object above ptmin in an annulus rmin..rmax around the given four-momentum...
Definition: Utils.hpp:284
void removeOverlap(MOMPTRS1 &momstofilter, const MOMPTRS2 &momstocmp, double deltaRMax, bool use_rapidity=false, double pTmax=DBL_MAX, double btageff=0)
Overlap removal – discard from first list if within deltaRMax of any from the second list Optional a...
Definition: Utils.hpp:299
double scalarSumPt(const std::vector< const Particle *> &particles, double pTlim=0.)
Scalar sum pT of particles with pT > pTlim (default pTlim = 0)
Definition: Utils.hpp:437
std::vector< const HEPUtils::Jet * > JetPtrs
Typedef for a vector of Jet pointers.
Definition: Utils.hpp:57
ParticlePtrs filter_reject(const ParticlePtrs &particles, std::function< bool(const Particle *)> rejfn, bool do_delete=true)
Filter a supplied particle vector by rejecting those which fail a supplied cut.
Definition: Utils.hpp:105
ParticlePtrs filter_select(const ParticlePtrs &particles, std::function< bool(const Particle *)> selfn, bool do_delete=true)
Filter a supplied particle vector by keeping those which pass a supplied cut.
Definition: Utils.hpp:113
FJNS::PseudoJet mk_pseudojet(const Vec4T &p)
Definition: Py8Utils.hpp:39
std::vector< double > makeBinValues(const std::vector< double > &binEdgeValues)
Alias.
Definition: Utils.hpp:262
HEPUtils::P4 mk_p4(const Vec4T &p)
Definition: Py8Utils.hpp:45
bool none_of(const CONTAINER &c, const FN &f)
Non-iterator version of std::none_of.
Definition: Utils.hpp:362
std::vector< std::vector< const HEPUtils::Particle * > > getOSpairs(std::vector< const HEPUtils::Particle *> particles)
Utility function for returning a collection of oppsosite-sign particle pairs.
Definition: Utils.cpp:86
std::vector< const HEPUtils::Particle * > ParticlePtrs
Typedef for a vector of Particle pointers.
Definition: Utils.hpp:54
std::map< const HEPUtils::Jet *, bool > generateBTagsMap(const std::vector< const HEPUtils::Jet *> &jets, double bTagEff, double cMissTagEff, double otherMissTagEff, double pTmin=0., double absEtaMax=DBL_MAX)
Return a map<Jet*,bool> containing a generated b-tag for every jet in the input vector.
Definition: Utils.hpp:216
std::vector< double > mk_bin_values(const std::vector< double > &binEdgeValues)
Make a vector of central bin values from a vector of bin edge values using linear interpolation...
Definition: Utils.hpp:254
bool cmpJetsByPt(const HEPUtils::Jet *jet1, const HEPUtils::Jet *jet2)
Comparison function to give a jets sorting order decreasing by pT.
Definition: Utils.hpp:398
int countPt(const std::vector< const Particle *> &particles, double pTlim)
Count number of particles that have pT > pTlim.
Definition: Utils.hpp:409
size_t binIndex(NUM1 val, const std::vector< NUM2 > &binedges, bool allow_overflow=false)
Definition: Utils.hpp:246
void ifilter_reject(ParticlePtrs &particles, std::function< bool(const Particle *)> rejfn, bool do_delete=true)
In-place filter a supplied particle vector by rejecting those which fail a supplied cut...
Definition: Utils.hpp:87
void sortByPt(ParticlePtrs &particles)
Definition: Utils.hpp:389
std::vector< std::vector< const HEPUtils::Particle * > > getSFOSpairs(std::vector< const HEPUtils::Particle *> particles)
Utility function for returning a collection of same-flavour, oppsosite-sign particle pairs...
Definition: Utils.cpp:69
void sortBy(ParticlePtrs &particles, std::function< bool(const Particle *, const Particle *)> cmpfn)
Particle-sorting function.
Definition: Utils.hpp:381
bool cmpParticlesByPt(const HEPUtils::Particle *lep1, const HEPUtils::Particle *lep2)
Comparison function to give a particles sorting order decreasing by pT.
Definition: Utils.hpp:386
hb_ModelParameters void
void filtereff_pt(std::vector< const HEPUtils::Particle *> &particles, const HEPUtils::BinnedFn1D< double > &eff_pt, bool do_delete=false)
Utility function for filtering a supplied particle vector by sampling wrt a binned 1D efficiency map ...
Definition: Utils.cpp:43
bool amIaJet(const HEPUtils::Jet *jet)
Identifier for jets true.
Definition: Utils.hpp:63
std::vector< std::shared_ptr< HEPUtils::Jet > > get_jets(const std::vector< MOM *> &moms, double R, double ptmin=0 *GeV, FJNS::JetAlgorithm alg=FJNS::antikt_algorithm)
Run jet clustering from any P4-compatible momentum type.
Definition: Utils.hpp:269
bool any_of(const CONTAINER &c, const FN &f)
Non-iterator version of std::any_of.
Definition: Utils.hpp:356
bool all_of(const CONTAINER &c, const FN &f)
Non-iterator version of std::all_of.
Definition: Utils.hpp:350
void iremoveerase(CONTAINER &c, const RMFN &fn)
Convenience combination of remove_if and erase.
Definition: Utils.hpp:81
std::vector< std::vector< const HEPUtils::Particle * > > getSSpairs(std::vector< const HEPUtils::Particle *> particles)
Utility function for returning a collection of same-sign particle pairs.
Definition: Utils.cpp:103
TODO: see if we can use this one:
Definition: Analysis.hpp:33
void ifilter_select(ParticlePtrs &particles, std::function< bool(const Particle *)> selfn, bool do_delete=true)
In-place filter a supplied particle vector by keeping those which pass a supplied cut...
Definition: Utils.hpp:97
void filtereff(std::vector< const 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
bool amIaBJet(const HEPUtils::Jet *jet)
Indentifier for b-jets true.
Definition: Utils.hpp:66
bool random_bool(double eff)
Return a random true/false at a success rate given by a number.
Definition: Utils.cpp:10