Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-09-18 09:31:17

0001 // -*- C++ -*-
0002 #ifndef RIVET_ExptSmearingFunctions_HH
0003 #define RIVET_ExptSmearingFunctions_HH
0004 
0005 #include "Rivet/Tools/MomentumSmearingFunctions.hh"
0006 #include "Rivet/Tools/ParticleSmearingFunctions.hh"
0007 #include "Rivet/Tools/JetSmearingFunctions.hh"
0008 
0009 namespace Rivet {
0010 
0011 
0012   /// @defgroup smearing Detector smearing & efficiency functions
0013   ///
0014   /// A family of smearing functions must include the following names and return types:
0015   ///
0016   ///   ParticleEffFn <PARTICLE>_EFF_<FAMILY>_<IDCLASS>
0017   ///   ParticleSmearFn <PARTICLE>_SMEAR_<FAMILY>
0018   ///   ParticleEffFn TRK_EFF_<FAMILY>
0019   ///   ParticleSmearFn TRK_SMEAR_<FAMILY>
0020   ///   JetEffFn JET_BTAG_<FAMILY>_<TAGGER>
0021   ///   JetSmearFn JET_SMEAR_<FAMILY>
0022   ///   METSmearParamsFn MET_SMEARPARAMS_<FAMILY>
0023   ///   METSmearFn MET_SMEAR_<FAMILY>
0024   ///
0025   /// where @c <PARTICLE> is all of @c ELECTRON, @c MUON, @c PHOTON, and @c TAU,
0026   /// @c <IDCLASS> is all of @c LOOSE, @c MEDIUM, and @c TIGHT. The family name
0027   /// and tagger name are free (within the constraints of function naming. Providing
0028   /// a full set may require use of alias/wrapper functions and/or the predefined
0029   /// eff/smearing @c IDENTITY functions.
0030   ///
0031   /// @{
0032 
0033 
0034   /// @defgroup smearing_elec Experiment-specific electron efficiency and smearing functions
0035   /// @{
0036 
0037   /// ATLAS Run 1 electron reconstruction efficiency
0038   /// @todo Include reco eff (but no e/y discrimination) in forward region
0039   inline double ELECTRON_RECOEFF_ATLAS_RUN1(const Particle& e) {
0040     if (e.abspid() != PID::ELECTRON) return 0;
0041     if (e.abseta() > 2.5) return 0;
0042     if (e.pT() < 10*GeV) return 0;
0043     return (e.abseta() < 1.5) ? 0.95 : 0.85;
0044   }
0045 
0046   /// ATLAS Run 2 electron reco efficiency
0047   ///
0048   /// Based on https://arxiv.org/pdf/1902.04655.pdf Fig 2
0049   inline double ELECTRON_RECOEFF_ATLAS_RUN2(const Particle& e) {
0050     if (e.abspid() != PID::ELECTRON) return 0;
0051     const double et = e.Et();
0052     if (e.abseta() > 2.5 || e.Et() < 2*GeV) return 0;
0053     if (et > 25*GeV) return 0.97;
0054     if (et > 10*GeV) return 0.92 + (et/GeV-10)/15.*0.05;
0055     if (et >  6*GeV) return 0.85 + (et/GeV-6)/4.*0.07;
0056     if (et >  5*GeV) return 0.70 + (et/GeV-5)/1.*0.15;
0057     if (et >  2*GeV) return 0.00 + (et/GeV-2)/3.*0.70;
0058     return 0;
0059   }
0060 
0061 
0062   /// @brief ATLAS Run 2 'loose' electron reco+identification efficiency
0063   ///
0064   /// Values read from Fig 3 of ATL-PHYS-PUB-2015-041
0065   ///
0066   /// @todo What about faking by jets or non-electrons?
0067   inline double ELECTRON_EFF_ATLAS_RUN2_LOOSE(const Particle& e) {
0068     if (e.abspid() != PID::ELECTRON) return 0;
0069 
0070     // Manually symmetrised eta eff histogram
0071     const static vector<double> edges_eta = { 0.0,   0.1,   0.8,   1.37,  1.52,  2.01,  2.37,  2.47 };
0072     const static vector<double> effs_eta  = { 0.950, 0.965, 0.955, 0.885, 0.950, 0.935, 0.90 };
0073     // Et eff histogram (10-20 is a guess)
0074     const static vector<double> edges_et = { 0,   10,   20,   25,   30,   35,   40,    45,    50,   60,  80 };
0075     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, 0.98 };//Extra value extrapolated for overflow
0076 
0077     if (e.abseta() > 2.47) return 0.0; // no ID outside the tracker
0078 
0079     const int i_eta = binIndex(e.abseta(), edges_eta);
0080     const int i_et = binIndex(e.Et()/GeV, edges_et, true);
0081     const double eff = effs_et[i_et] * effs_eta[i_eta] / 0.95; //< norm factor as approximate double differential
0082     return min(eff, 1.0) * ELECTRON_RECOEFF_ATLAS_RUN2(e);
0083   }
0084 
0085   /// @brief Pretend that ATLAS Run 1 loose was the same as in Run 2
0086   inline double ELECTRON_EFF_ATLAS_RUN1_LOOSE(const Particle& e) {
0087     return ELECTRON_EFF_ATLAS_RUN2_LOOSE(e);
0088   }
0089 
0090 
0091   /// @brief ATLAS Run 1 'medium' electron reco+identification efficiency
0092   inline double ELECTRON_EFF_ATLAS_RUN1_MEDIUM(const Particle& e) {
0093     if (e.abspid() != PID::ELECTRON) return 0;
0094 
0095     const static vector<double> eta_edges_10 = {0.000, 0.049, 0.454, 1.107, 1.46, 1.790, 2.277, 2.500};
0096     const static vector<double> eta_vals_10  = {0.730, 0.757, 0.780, 0.771, 0.77, 0.777, 0.778};
0097 
0098     const static vector<double> eta_edges_15 = {0.000, 0.053, 0.456, 1.102, 1.463, 1.783, 2.263, 2.500};
0099     const static vector<double> eta_vals_15  = {0.780, 0.800, 0.819, 0.759, 0.749, 0.813, 0.829};
0100 
0101     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};
0102     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};
0103 
0104     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};
0105     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};
0106 
0107     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};
0108     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};
0109 
0110     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};
0111     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};
0112 
0113     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};
0114     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};
0115 
0116     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};
0117     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};
0118 
0119     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};
0120     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};
0121 
0122     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};
0123     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};
0124 
0125     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};
0126     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};
0127 
0128     const static vector<double> et_edges = { 10, 15, 20, 25, 30, 35, 40, 45, 50, 60, 80 };
0129     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 };
0130     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 };
0131 
0132     if (e.abseta() > 2.5 || e.Et() < 10*GeV) return 0.0;
0133     const int i_et = binIndex(e.Et()/GeV, et_edges, true);
0134     const int i_eta = binIndex(e.abseta(), et_eta_edges[i_et]);
0135     return et_eta_vals[i_et][i_eta] * ELECTRON_RECOEFF_ATLAS_RUN1(e);
0136   }
0137 
0138   /// @brief ATLAS Run 2 'medium' electron reco+identification efficiency
0139   ///
0140   /// ~1% increase over Run 1 informed by Fig 1 in https://cds.cern.ch/record/2157687/files/ATLAS-CONF-2016-024.pdf
0141   inline double ELECTRON_EFF_ATLAS_RUN2_MEDIUM(const Particle& e) {
0142     if (e.abspid() != PID::ELECTRON) return 0;
0143     return 1.01 * ELECTRON_EFF_ATLAS_RUN1_MEDIUM(e);
0144   }
0145 
0146 
0147   /// @brief ATLAS Run 1 'tight' electron reco+identification efficiency
0148   inline double ELECTRON_EFF_ATLAS_RUN1_TIGHT(const Particle& e) {
0149     if (e.abspid() != PID::ELECTRON) return 0;
0150 
0151     const static vector<double> eta_edges_10 = {0.000, 0.049, 0.459, 1.100, 1.461, 1.789, 2.270, 2.500};
0152     const static vector<double> eta_vals_10  = {0.581, 0.632, 0.668, 0.558, 0.548, 0.662, 0.690};
0153 
0154     const static vector<double> eta_edges_15 = {0.000, 0.053, 0.450, 1.096, 1.463, 1.783, 2.269, 2.500};
0155     const static vector<double> eta_vals_15 =  {0.630, 0.678, 0.714, 0.633, 0.616, 0.700, 0.733};
0156 
0157     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};
0158     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};
0159 
0160     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};
0161     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};
0162 
0163     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};
0164     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};
0165 
0166     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};
0167     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};
0168 
0169     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};
0170     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};
0171 
0172     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};
0173     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};
0174 
0175     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};
0176     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};
0177 
0178     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};
0179     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};
0180 
0181     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};
0182     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};
0183 
0184     const static vector<double> et_edges = { 10, 15, 20, 25, 30, 35, 40, 45, 50, 60, 80 };
0185     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 };
0186     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 };
0187 
0188     if (e.abseta() > 2.5 || e.Et() < 10*GeV) return 0.0;
0189     const int i_et = binIndex(e.Et()/GeV, et_edges, true);
0190     const int i_eta = binIndex(e.abseta(), et_eta_edges[i_et]);
0191     return et_eta_vals[i_et][i_eta] * ELECTRON_RECOEFF_ATLAS_RUN1(e);
0192   }
0193 
0194   /// @brief ATLAS Run 2 'tight' electron reco+identification efficiency
0195   ///
0196   /// ~1% increase over Run 1 informed by Fig 1 in https://cds.cern.ch/record/2157687/files/ATLAS-CONF-2016-024.pdf
0197   inline double ELECTRON_EFF_ATLAS_RUN2_TIGHT(const Particle& e) {
0198     if (e.abspid() != PID::ELECTRON) return 0;
0199     const static vector<double> et_edges = { /* 10, 15, */ 20, 25, 30, 35, 40, 45, 50, 60, 80 };
0200     const static vector<double> et_effs = { 0.785, 0.805, 0.820, 0.830, 0.840, 0.850, 0.875, 0.910, 0.910 }; //Extra value extrapolated for overflow
0201     const static vector<double> eta_edges = {0.000, 0.051, 0.374, 0.720, 0.981, 1.279, 1.468, 1.707, 1.945, 2.207, 2.457, 2.500}; // from ET > 80 bin
0202     const static vector<double> eta_refs =    {0.819, 0.855, 0.899, 0.906, 0.900, 0.869, 0.865, 0.873, 0.869, 0.868, 0.859};
0203     if (e.abseta() > 2.5 || e.Et() < 20*GeV) return 0.0;
0204     const int i_et = binIndex(e.Et()/GeV, et_edges, true);
0205     const int i_eta = binIndex(e.abseta(), eta_edges);
0206     const double eff_et = et_effs[i_et]; //< integral eff
0207     // Scale to |eta| shape, following the ~85% efficient high-ET bin from Run 1
0208     const double eff = eff_et * (eta_refs[i_eta]/0.85) * ELECTRON_RECOEFF_ATLAS_RUN2(e);
0209     //return ELECTRON_IDEFF_ATLAS_RUN1_TIGHT(e);
0210     return eff;
0211   }
0212 
0213 
0214 
0215   /// ATLAS Run 1 electron reco smearing
0216   inline Particle ELECTRON_SMEAR_ATLAS_RUN1(const Particle& e) {
0217     static const vector<double> edges_eta = {0., 2.5, 3.};
0218     static const vector<double> edges_pt = {0., 0.1, 25.};
0219     static const vector<double> e2s = {0.000, 0.015, 0.005,
0220                                        0.005, 0.005, 0.005,
0221                                        0.107, 0.107, 0.107};
0222     static const vector<double> es = {0.00, 0.00, 0.05,
0223                                       0.05, 0.05, 0.05,
0224                                       2.08, 2.08, 2.08};
0225     static const vector<double> cs = {0.00, 0.00, 0.25,
0226                                       0.25, 0.25, 0.25,
0227                                       0.00, 0.00, 0.00};
0228 
0229     const int i_eta = binIndex(e.abseta(), edges_eta, true);
0230     const int i_pt = binIndex(e.pT()/GeV, edges_pt, true);
0231     const int i = i_eta*edges_pt.size() + i_pt;
0232 
0233     // Calculate absolute resolution in GeV
0234     const double c1 = sqr(e2s[i]), c2 = sqr(es[i]), c3 = sqr(cs[i]);
0235     const double resolution = sqrt(c1*e.E2() + c2*e.E() + c3) * GeV;
0236 
0237     // normal_distribution<> d(e.E(), resolution);
0238     // const double mass = e.mass2() > 0 ? e.mass() : 0; //< numerical carefulness...
0239     // const double smeared_E = max(d(gen), mass); //< can't let the energy go below the mass!
0240     // return Particle(e.pid(), FourMomentum::mkEtaPhiME(e.eta(), e.phi(), mass, smeared_E));
0241     return Particle(e.pid(), P4_SMEAR_E_GAUSS(e, resolution));
0242   }
0243 
0244 
0245   /// ATLAS Run 2 electron reco smearing
0246   ///
0247   /// @todo Currently just a copy of the Run 1 version: fix!
0248   inline Particle ELECTRON_SMEAR_ATLAS_RUN2(const Particle& e) {
0249     return ELECTRON_SMEAR_ATLAS_RUN1(e);
0250   }
0251 
0252 
0253   /// @todo Add charge flip efficiency?
0254 
0255 
0256 
0257   /// CMS Run 1 electron reconstruction efficiency
0258   inline double ELECTRON_EFF_CMS_RUN1(const Particle& e) {
0259     if (e.abspid() != PID::ELECTRON) return 0;
0260     if (e.abseta() > 2.5) return 0;
0261     if (e.pT() < 10*GeV) return 0;
0262     return (e.abseta() < 1.5) ? 0.95 : 0.85;
0263   }
0264   /// CMS Run 1 loose electron reconstruction efficiency
0265   /// @todo Just an alias to generic: improve!
0266   inline double ELECTRON_EFF_CMS_RUN1_LOOSE(const Particle& e) { return ELECTRON_EFF_CMS_RUN1(e); }
0267   /// CMS Run 1 medium electron reconstruction efficiency
0268   /// @todo Just an alias to generic: improve!
0269   inline double ELECTRON_EFF_CMS_RUN1_MEDIUM(const Particle& e) { return ELECTRON_EFF_CMS_RUN1(e); }
0270   /// CMS Run 1 tight electron reconstruction efficiency
0271   /// @todo Just an alias to generic: improve!
0272   inline double ELECTRON_EFF_CMS_RUN1_TIGHT(const Particle& e) { return ELECTRON_EFF_CMS_RUN1(e); }
0273 
0274 
0275   /// CMS Run 2 electron reco efficiency
0276   ///
0277   /// @todo Currently just a copy of Run 1: fix!
0278   inline double ELECTRON_EFF_CMS_RUN2(const Particle& e) {
0279     if (e.abspid() != PID::ELECTRON) return 0;
0280     return ELECTRON_EFF_CMS_RUN1(e);
0281   }
0282   /// CMS Run 2 loose electron reconstruction efficiency
0283   /// @todo Just an alias to generic: improve!
0284   inline double ELECTRON_EFF_CMS_RUN2_LOOSE(const Particle& e) { return ELECTRON_EFF_CMS_RUN2(e); }
0285   /// CMS Run 2 medium electron reconstruction efficiency
0286   /// @todo Just an alias to generic: improve!
0287   inline double ELECTRON_EFF_CMS_RUN2_MEDIUM(const Particle& e) { return ELECTRON_EFF_CMS_RUN2(e); }
0288   /// CMS Run 2 tight electron reconstruction efficiency
0289   /// @todo Just an alias to generic: improve!
0290   inline double ELECTRON_EFF_CMS_RUN2_TIGHT(const Particle& e) { return ELECTRON_EFF_CMS_RUN2(e); }
0291 
0292 
0293   /// @brief CMS electron energy smearing, preserving direction
0294   ///
0295   /// Calculate resolution
0296   /// for pT > 0.1 GeV, E resolution = |eta| < 0.5 -> sqrt(0.06^2 + pt^2 * 1.3e-3^2)
0297   ///                                  |eta| < 1.5 -> sqrt(0.10^2 + pt^2 * 1.7e-3^2)
0298   ///                                  |eta| < 2.5 -> sqrt(0.25^2 + pt^2 * 3.1e-3^2)
0299   inline Particle ELECTRON_SMEAR_CMS_RUN1(const Particle& e) {
0300     // Calculate absolute resolution in GeV from functional form
0301     double resolution = 0;
0302     const double abseta = e.abseta();
0303     if (e.pT() > 0.1*GeV && abseta < 2.5) { //< should be a given from efficiencies
0304       if (abseta < 0.5) {
0305         resolution = add_quad(0.06, 1.3e-3 * e.pT()/GeV) * GeV;
0306       } else if (abseta < 1.5) {
0307         resolution = add_quad(0.10, 1.7e-3 * e.pT()/GeV) * GeV;
0308       } else { // still |eta| < 2.5
0309         resolution = add_quad(0.25, 3.1e-3 * e.pT()/GeV) * GeV;
0310       }
0311     }
0312 
0313     // normal_distribution<> d(e.E(), resolution);
0314     // const double mass = e.mass2() > 0 ? e.mass() : 0; //< numerical carefulness...
0315     // const double smeared_E = max(d(gen), mass); //< can't let the energy go below the mass!
0316     // return Particle(e.pid(), FourMomentum::mkEtaPhiME(e.eta(), e.phi(), mass, smeared_E));
0317     return Particle(e.pid(), P4_SMEAR_E_GAUSS(e, resolution));
0318   }
0319 
0320   /// CMS Run 2 electron reco smearing
0321   ///
0322   /// @todo Currently just a copy of the Run 1 version: fix!
0323   inline Particle ELECTRON_SMEAR_CMS_RUN2(const Particle& e) {
0324     return ELECTRON_SMEAR_CMS_RUN1(e);
0325   }
0326 
0327   /// @}
0328 
0329 
0330 
0331   /// @defgroup smearing_photon Experiment-specific photon efficiency and smearing functions
0332   /// @{
0333 
0334   /// @brief ATLAS Run 2 photon reco efficiency
0335   ///
0336   /// Taken from converted photons, Fig 8, in arXiv:1606.01813
0337   inline double PHOTON_EFF_ATLAS_RUN1(const Particle& y) {
0338     if (y.abspid() != PID::PHOTON) return 0; ///< @todo Allow electron misID? What about jet misID?
0339 
0340     if (y.pT() < 10*GeV) return 0;
0341     if (inRange(y.abseta(), 1.37, 1.52) || y.abseta() > 2.37) return 0;
0342 
0343     static const vector<double> edges_eta = {0., 0.6, 1.37, 1.52, 1.81, 2.37};
0344     static const vector<double> edges_pt = {10., 15., 20., 25., 30., 35., 40., 45.,
0345                                             50., 60., 80., 100., 125., 150., 175., 250.};
0346     static const vector<double> effs = {0.53, 0.65, 0.73, 0.83, 0.86, 0.93, 0.94, 0.96,
0347                                         0.97, 0.98, 0.98, 0.98, 0.98, 0.98, 0.98, 0.98,//
0348                                         0.45, 0.57, 0.67, 0.74, 0.84, 0.87, 0.93, 0.94,
0349                                         0.95, 0.96, 0.97, 0.98, 0.98, 0.99, 0.99, 0.99,//
0350                                         0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00,
0351                                         0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00,//
0352                                         0.48, 0.56, 0.68, 0.76, 0.86, 0.90, 0.93, 0.95,
0353                                         0.96, 0.97, 0.98, 0.99, 0.99, 1.00, 1.00, 1.00,//
0354                                         0.50, 0.61, 0.74, 0.82, 0.88, 0.92, 0.94, 0.95,
0355                                         0.96, 0.97, 0.98, 0.98, 0.98, 0.98, 0.99, 0.99};
0356 
0357     const int i_eta = binIndex(y.abseta(), edges_eta);
0358     const int i_pt = binIndex(y.pT()/GeV, edges_pt, true);
0359     const int i = i_eta*edges_pt.size() + i_pt;
0360     const double eff = effs[i];
0361     return eff;
0362   }
0363   /// ATLAS Run 1 loose photon reconstruction efficiency
0364   /// @todo Just an alias to generic: improve!
0365   inline double PHOTON_EFF_ATLAS_RUN1_LOOSE(const Particle& y) { return PHOTON_EFF_ATLAS_RUN1(y); }
0366   /// ATLAS Run 1 medium photon reconstruction efficiency
0367   /// @todo Just an alias to generic: improve!
0368   inline double PHOTON_EFF_ATLAS_RUN1_MEDIUM(const Particle& y) { return PHOTON_EFF_ATLAS_RUN1(y); }
0369   /// ATLAS Run 1 tight photon reconstruction efficiency
0370   /// @todo Just an alias to generic: improve!
0371   inline double PHOTON_EFF_ATLAS_RUN1_TIGHT(const Particle& y) { return PHOTON_EFF_ATLAS_RUN1(y); }
0372 
0373 
0374   /// @brief ATLAS Run 2 photon reco efficiency
0375   ///
0376   /// Taken from converted photons, Fig 6, in ATL-PHYS-PUB-2016-014
0377   inline double PHOTON_EFF_ATLAS_RUN2(const Particle& y) {
0378     if (y.abspid() != PID::PHOTON) return 0; ///< @todo Allow electron misID? What about jet misID?
0379 
0380     if (y.pT() < 10*GeV) return 0;
0381     if (inRange(y.abseta(), 1.37, 1.52) || y.abseta() > 2.37) return 0;
0382 
0383     static const vector<double> edges_eta = {0., 0.6, 1.37, 1.52, 1.81, 2.37};
0384     static const vector<double> edges_pt = {10., 15., 20., 25., 30., 35., 40., 45.,
0385                                             50., 60., 80., 100., 125., 150., 175., 250.};
0386     static const vector<double> effs = {0.55, 0.70, 0.85, 0.89, 0.93, 0.95, 0.96, 0.96,
0387                                         0.97, 0.97, 0.98, 0.97, 0.97, 0.97, 0.97, 0.97,//
0388                                         0.47, 0.66, 0.79, 0.86, 0.89, 0.94, 0.96, 0.97,
0389                                         0.97, 0.98, 0.97, 0.98, 0.98, 0.98, 0.98, 0.98,//
0390                                         0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00,
0391                                         0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00,//
0392                                         0.54, 0.71, 0.84, 0.88, 0.92, 0.93, 0.94, 0.95,
0393                                         0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96,//
0394                                         0.61, 0.74, 0.83, 0.88, 0.91, 0.94, 0.95, 0.96,
0395                                         0.97, 0.98, 0.98, 0.98, 0.98, 0.98, 0.98, 0.98};
0396 
0397     const int i_eta = binIndex(y.abseta(), edges_eta);
0398     const int i_pt = binIndex(y.pT()/GeV, edges_pt, true);
0399     const int i = i_eta*edges_pt.size() + i_pt;
0400     const double eff = effs[i];
0401     return eff;
0402   }
0403   /// ATLAS Run 2 loose photon reconstruction efficiency
0404   /// @todo Just an alias to generic: improve!
0405   inline double PHOTON_EFF_ATLAS_RUN2_LOOSE(const Particle& y) { return PHOTON_EFF_ATLAS_RUN2(y); }
0406   /// ATLAS Run 2 medium photon reconstruction efficiency
0407   /// @todo Just an alias to generic: improve!
0408   inline double PHOTON_EFF_ATLAS_RUN2_MEDIUM(const Particle& y) { return PHOTON_EFF_ATLAS_RUN2(y); }
0409   /// ATLAS Run 2 tight photon reconstruction efficiency
0410   /// @todo Just an alias to generic: improve!
0411   inline double PHOTON_EFF_ATLAS_RUN2_TIGHT(const Particle& y) { return PHOTON_EFF_ATLAS_RUN2(y); }
0412 
0413 
0414   /// CMS Run 1 photon reco efficiency
0415   ///
0416   /// @todo Currently from Delphes
0417   inline double PHOTON_EFF_CMS_RUN1(const Particle& y) {
0418     if (y.abspid() != PID::PHOTON) return 0; ///< @todo Allow electron misID? What about jet misID?
0419     if (y.pT() < 10*GeV || y.abseta() > 2.5) return 0;
0420     return (y.abseta() < 1.5) ? 0.95 : 0.85;
0421   }
0422   /// CMS Run 1 loose photon reconstruction efficiency
0423   /// @todo Just an alias to generic: improve!
0424   inline double PHOTON_EFF_CMS_RUN1_LOOSE(const Particle& y) { return PHOTON_EFF_CMS_RUN1(y); }
0425   /// CMS Run 1 medium photon reconstruction efficiency
0426   /// @todo Just an alias to generic: improve!
0427   inline double PHOTON_EFF_CMS_RUN1_MEDIUM(const Particle& y) { return PHOTON_EFF_CMS_RUN1(y); }
0428   /// CMS Run 1 tight photon reconstruction efficiency
0429   /// @todo Just an alias to generic: improve!
0430   inline double PHOTON_EFF_CMS_RUN1_TIGHT(const Particle& y) { return PHOTON_EFF_CMS_RUN1(y); }
0431 
0432 
0433   /// CMS Run 2 photon reco efficiency
0434   ///
0435   /// @todo Currently just a copy of Run 1: fix!
0436   inline double PHOTON_EFF_CMS_RUN2(const Particle& y) {
0437     if (y.abspid() != PID::PHOTON) return 0; ///< @todo Allow electron misID? What about jet misID?
0438     return PHOTON_EFF_CMS_RUN1(y);
0439   }
0440   /// CMS Run 2 loose photon reconstruction efficiency
0441   /// @todo Just an alias to generic: improve!
0442   inline double PHOTON_EFF_CMS_RUN2_LOOSE(const Particle& y) { return PHOTON_EFF_CMS_RUN2(y); }
0443   /// CMS Run 2 medium photon reconstruction efficiency
0444   /// @todo Just an alias to generic: improve!
0445   inline double PHOTON_EFF_CMS_RUN2_MEDIUM(const Particle& y) { return PHOTON_EFF_CMS_RUN2(y); }
0446   /// CMS Run 2 tight photon reconstruction efficiency
0447   /// @todo Just an alias to generic: improve!
0448   inline double PHOTON_EFF_CMS_RUN2_TIGHT(const Particle& y) { return PHOTON_EFF_CMS_RUN2(y); }
0449 
0450 
0451   /// @todo Use real photon smearing
0452   inline Particle PHOTON_SMEAR_ATLAS_RUN1(const Particle& y) { return y; }
0453   inline Particle PHOTON_SMEAR_ATLAS_RUN2(const Particle& y) { return y; }
0454   inline Particle PHOTON_SMEAR_CMS_RUN1(const Particle& y) { return y; }
0455   inline Particle PHOTON_SMEAR_CMS_RUN2(const Particle& y) { return y; }
0456 
0457   /// @}
0458 
0459 
0460 
0461   /// @defgroup smearing_muon Experiment-specific muon efficiency and smearing functions
0462   /// @{
0463 
0464   /// ATLAS Run 1 muon reco efficiency
0465   inline double MUON_EFF_ATLAS_RUN1(const Particle& m) {
0466     if (m.abspid() != PID::MUON) return 0;
0467     if (m.abseta() > 2.7) return 0;
0468     if (m.pT() < 10*GeV) return 0;
0469     return (m.abseta() < 1.5) ? 0.95 : 0.85;
0470   }
0471   /// ATLAS Run 1 loose muon reconstruction efficiency
0472   /// @todo Just an alias to generic: improve!
0473   inline double MUON_EFF_ATLAS_RUN1_LOOSE(const Particle& m) { return MUON_EFF_ATLAS_RUN1(m); }
0474   /// ATLAS Run 1 medium muon reconstruction efficiency
0475   /// @todo Just an alias to generic: improve!
0476   inline double MUON_EFF_ATLAS_RUN1_MEDIUM(const Particle& m) { return MUON_EFF_ATLAS_RUN1(m); }
0477   /// ATLAS Run 1 tight muon reconstruction efficiency
0478   /// @todo Just an alias to generic: improve!
0479   inline double MUON_EFF_ATLAS_RUN1_TIGHT(const Particle& m) { return MUON_EFF_ATLAS_RUN1(m); }
0480 
0481 
0482   /// ATLAS Run 2 muon reco efficiency
0483   ///
0484   /// From https://arxiv.org/pdf/1603.05598.pdf , Fig3 top
0485   inline double MUON_RECOEFF_ATLAS_RUN2(const Particle& m) {
0486     if (m.abspid() != PID::MUON) return 0;
0487     if (m.abseta() > 2.5) return 0;
0488     if (m.abseta() < 0.1) return 0.61;
0489     // if (m.pT() < 10*GeV) return 0;
0490     return (m.abseta() < 1) ? 0.98 : 0.99;
0491   }
0492 
0493   /// @brief ATLAS Run 2 muon reco+ID efficiency
0494   ///
0495   /// For medium ID, from Fig 3 of
0496   /// https://cds.cern.ch/record/2047831/files/ATL-PHYS-PUB-2015-037.pdf
0497   inline double MUON_EFF_ATLAS_RUN2(const Particle& m) {
0498     if (m.abspid() != PID::MUON) return 0;
0499     if (m.abseta() > 2.7) return 0;
0500     static const vector<double> edges_pt = {0., 3.5, 4., 5., 6., 7., 8., 10.};
0501     static const vector<double> effs = {0.00, 0.76, 0.94, 0.97, 0.98, 0.98, 0.98, 0.99};
0502     const int i_pt = binIndex(m.pT()/GeV, edges_pt, true);
0503     const double eff = effs[i_pt] * MUON_RECOEFF_ATLAS_RUN2(m);
0504     return eff;
0505   }
0506   /// ATLAS Run 2 loose muon reconstruction efficiency
0507   /// @todo Add muon loose/medium/tight ID efficiencies? All around 95-98%... ignore?
0508   /// @todo Just an alias to generic: improve!
0509   inline double MUON_EFF_ATLAS_RUN2_LOOSE(const Particle& m) { return MUON_EFF_ATLAS_RUN2(m); }
0510   /// ATLAS Run 2 medium muon reconstruction efficiency
0511   /// @todo Add muon loose/medium/tight ID efficiencies? All around 95-98%... ignore?
0512   /// @todo Just an alias to generic: improve!
0513   inline double MUON_EFF_ATLAS_RUN2_MEDIUM(const Particle& m) { return MUON_EFF_ATLAS_RUN2(m); }
0514   /// ATLAS Run 2 tight muon reconstruction efficiency
0515   /// @todo Add muon loose/medium/tight ID efficiencies? All around 95-98%... ignore?
0516   /// @todo Just an alias to generic: improve!
0517   inline double MUON_EFF_ATLAS_RUN2_TIGHT(const Particle& m) { return MUON_EFF_ATLAS_RUN2(m); }
0518 
0519 
0520   /// ATLAS Run 1 muon reco smearing
0521   inline Particle MUON_SMEAR_ATLAS_RUN1(const Particle& m) {
0522     static const vector<double> edges_eta = {0, 1.5, 2.5};
0523     static const vector<double> edges_pt = {0, 0.1, 1.0, 10., 200.};
0524     static const vector<double> res = {0., 0.03, 0.02, 0.03, 0.05,
0525                                        0., 0.04, 0.03, 0.04, 0.05};
0526 
0527     const int i_eta = binIndex(m.abseta(), edges_eta);
0528     const int i_pt = binIndex(m.pT()/GeV, edges_pt, true);
0529     const int i = i_eta*edges_pt.size() + i_pt;
0530 
0531     const double resolution = res[i];
0532 
0533     // Smear by a Gaussian centered on the current pT, with width given by the resolution
0534     // normal_distribution<> d(m.pT(), resolution*m.pT());
0535     // const double smeared_pt = max(d(gen), 0.);
0536     // const double mass = m.mass2() > 0 ? m.mass() : 0; //< numerical carefulness...
0537     // return Particle(m.pid(), FourMomentum::mkEtaPhiMPt(m.eta(), m.phi(), mass, smeared_pt));
0538     return Particle(m.pid(), P4_SMEAR_PT_GAUSS(m, resolution*m.pT()));
0539   }
0540 
0541   /// ATLAS Run 2 muon reco smearing
0542   ///
0543   /// From https://arxiv.org/abs/1603.05598 , eq (10) and Fig 12
0544   inline Particle MUON_SMEAR_ATLAS_RUN2(const Particle& m) {
0545     double mres_pt = 0.015;
0546     if (m.pT() > 50*GeV) mres_pt = 0.014 + 0.01*(m.pT()/GeV-50)/50;
0547     if (m.pT() > 100*GeV) mres_pt = 0.025;
0548     const double ptres_pt = SQRT2 * mres_pt; //< from Eq (10)
0549     const double resolution = (m.abseta() < 1.5 ? 1.0 : 1.25) * ptres_pt;
0550     return Particle(m.pid(), P4_SMEAR_PT_GAUSS(m, resolution*m.pT()));
0551   }
0552 
0553 
0554 
0555   /// CMS Run 1 muon reco efficiency
0556   inline double MUON_EFF_CMS_RUN1(const Particle& m) {
0557     if (m.abspid() != PID::MUON) return 0;
0558     if (m.abseta() > 2.4) return 0;
0559     if (m.pT() < 10*GeV) return 0;
0560     return 0.95 * (m.abseta() < 1.5 ? 1 : exp(0.5 - 5e-4*m.pT()/GeV));
0561   }
0562   /// CMS Run 1 loose muon reconstruction efficiency
0563   /// @todo Add muon loose/medium/tight ID efficiencies? All around 95-98%... ignore?
0564   /// @todo Just an alias to generic: improve!
0565   inline double MUON_EFF_CMS_RUN1_LOOSE(const Particle& m) { return MUON_EFF_CMS_RUN1(m); }
0566   /// CMS Run 1 medium muon reconstruction efficiency
0567   /// @todo Add muon loose/medium/tight ID efficiencies? All around 95-98%... ignore?
0568   /// @todo Just an alias to generic: improve!
0569   inline double MUON_EFF_CMS_RUN1_MEDIUM(const Particle& m) { return MUON_EFF_CMS_RUN1(m); }
0570   /// CMS Run 1 tight muon reconstruction efficiency
0571   /// @todo Add muon loose/medium/tight ID efficiencies? All around 95-98%... ignore?
0572   /// @todo Just an alias to generic: improve!
0573   inline double MUON_EFF_CMS_RUN1_TIGHT(const Particle& m) { return MUON_EFF_CMS_RUN1(m); }
0574 
0575 
0576   /// CMS Run 2 muon reco efficiency
0577   ///
0578   /// @todo Currently just a copy of Run 1: fix!
0579   inline double MUON_EFF_CMS_RUN2(const Particle& m) {
0580     return MUON_EFF_CMS_RUN1(m);
0581   }
0582   /// CMS Run 2 loose muon reconstruction efficiency
0583   /// @todo Add muon loose/medium/tight ID efficiencies? All around 95-98%... ignore?
0584   /// @todo Just an alias to generic: improve!
0585   inline double MUON_EFF_CMS_RUN2_LOOSE(const Particle& m) { return MUON_EFF_CMS_RUN2(m); }
0586   /// CMS Run 2 medium muon reconstruction efficiency
0587   /// @todo Add muon loose/medium/tight ID efficiencies? All around 95-98%... ignore?
0588   /// @todo Just an alias to generic: improve!
0589   inline double MUON_EFF_CMS_RUN2_MEDIUM(const Particle& m) { return MUON_EFF_CMS_RUN2(m); }
0590   /// CMS Run 2 tight muon reconstruction efficiency
0591   /// @todo Add muon loose/medium/tight ID efficiencies? All around 95-98%... ignore?
0592   /// @todo Just an alias to generic: improve!
0593   inline double MUON_EFF_CMS_RUN2_TIGHT(const Particle& m) { return MUON_EFF_CMS_RUN2(m); }
0594 
0595 
0596 
0597   /// CMS Run 1 muon reco smearing
0598   inline Particle MUON_SMEAR_CMS_RUN1(const Particle& m) {
0599     // Calculate fractional resolution
0600     // for pT > 0.1 GeV, mom resolution = |eta| < 0.5 -> sqrt(0.01^2 + pt^2 * 2.0e-4^2)
0601     //                                    |eta| < 1.5 -> sqrt(0.02^2 + pt^2 * 3.0e-4^2)
0602     //                                    |eta| < 2.5 -> sqrt(0.05^2 + pt^2 * 2.6e-4^2)
0603     double resolution = 0;
0604     const double abseta = m.abseta();
0605     if (m.pT() > 0.1*GeV && abseta < 2.5) {
0606       if (abseta < 0.5) {
0607         resolution = add_quad(0.01, 2.0e-4 * m.pT()/GeV);
0608       } else if (abseta < 1.5) {
0609         resolution = add_quad(0.02, 3.0e-4 * m.pT()/GeV);
0610       } else { // still |eta| < 2.5... but isn't CMS' mu acceptance < 2.4?
0611         resolution = add_quad(0.05, 2.6e-4 * m.pT()/GeV);
0612       }
0613     }
0614 
0615     // Smear by a Gaussian centered on the current pT, with width given by the resolution
0616     // normal_distribution<> d(m.pT(), resolution*m.pT());
0617     // const double smeared_pt = max(d(gen), 0.);
0618     // const double mass = m.mass2() > 0 ? m.mass() : 0; //< numerical carefulness...
0619     // return Particle(m.pid(), FourMomentum::mkEtaPhiMPt(m.eta(), m.phi(), mass, smeared_pt));
0620     return Particle(m.pid(), P4_SMEAR_PT_GAUSS(m, resolution*m.pT()));
0621   }
0622 
0623   /// CMS Run 2 muon reco smearing
0624   ///
0625   /// @todo Currently just a copy of the Run 1 version: fix!
0626   inline Particle MUON_SMEAR_CMS_RUN2(const Particle& m) {
0627     return MUON_SMEAR_CMS_RUN1(m);
0628   }
0629 
0630   /// @}
0631 
0632 
0633 
0634   /// @defgroup smearing_tau Experiment-specific tau efficiency and smearing functions
0635   /// @{
0636 
0637   /// @brief ATLAS Run 1 8 TeV tau efficiencies (medium working point)
0638   ///
0639   /// Taken from http://arxiv.org/pdf/1412.7086.pdf
0640   ///   20-40 GeV 1-prong LMT eff|mis = 0.66|1/10, 0.56|1/20, 0.36|1/80
0641   ///   20-40 GeV 3-prong LMT eff|mis = 0.45|1/60, 0.38|1/100, 0.27|1/300
0642   ///   > 40 GeV 1-prong LMT eff|mis = 0.66|1/15, 0.56|1/25, 0.36|1/80
0643   ///   > 40 GeV 3-prong LMT eff|mis = 0.45|1/250, 0.38|1/400, 0.27|1/1300
0644   inline double TAU_EFF_ATLAS_RUN1_MEDIUM(const Particle& t) {
0645     if (t.abseta() > 2.5) return 0; //< hmm... mostly
0646     if (inRange(t.abseta(), 1.37, 1.52)) return 0; //< crack region
0647     double pThadvis = 0;
0648     Particles chargedhadrons;
0649     for (const Particle& p : t.children()) {
0650       if (p.isHadron()) {
0651         pThadvis += p.pT(); //< right definition? Paper is unclear
0652         if (p.charge3() != 0 && p.abseta() < 2.5 && p.pT() > 1*GeV) chargedhadrons += p;
0653       }
0654     }
0655     if (chargedhadrons.empty()) return 0; //< leptonic tau
0656     if (pThadvis < 20*GeV) return 0; //< below threshold
0657     if (pThadvis < 40*GeV) {
0658       if (chargedhadrons.size() == 1) return (t.abspid() == PID::TAU) ? 0.56 : 0; //1/20.;
0659       if (chargedhadrons.size() == 3) return (t.abspid() == PID::TAU) ? 0.38 : 0; //1/100.;
0660     } else {
0661       if (chargedhadrons.size() == 1) return (t.abspid() == PID::TAU) ? 0.56 : 0; //1/25.;
0662       if (chargedhadrons.size() == 3) return (t.abspid() == PID::TAU) ? 0.38 : 0; //1/400.;
0663     }
0664     return 0;
0665   }
0666   /// ATLAS Run 1 medium muon reconstruction efficiency
0667   /// @todo Just an alias to generic: improve!
0668   inline double TAU_EFF_ATLAS_RUN1(const Particle& t) { return TAU_EFF_ATLAS_RUN1_MEDIUM(t); }
0669   /// ATLAS Run 1 loose tau reconstruction efficiency
0670   /// @todo Just an alias to generic: improve!
0671   inline double TAU_EFF_ATLAS_RUN1_LOOSE(const Particle& t) { return TAU_EFF_ATLAS_RUN1(t); }
0672   /// ATLAS Run 1 tight muon reconstruction efficiency
0673   /// @todo Just an alias to generic: improve!
0674   inline double TAU_EFF_ATLAS_RUN1_TIGHT(const Particle& t) { return TAU_EFF_ATLAS_RUN1(t); }
0675 
0676 
0677   /// @brief ATLAS Run 1 8 TeV tau misID rates (medium working point)
0678   ///
0679   /// Taken from http://arxiv.org/pdf/1412.7086.pdf
0680   ///   20-40 GeV 1-prong LMT eff|mis = 0.66|1/10, 0.56|1/20, 0.36|1/80
0681   ///   20-40 GeV 3-prong LMT eff|mis = 0.45|1/60, 0.38|1/100, 0.27|1/300
0682   ///   > 40 GeV 1-prong LMT eff|mis = 0.66|1/15, 0.56|1/25, 0.36|1/80
0683   ///   > 40 GeV 3-prong LMT eff|mis = 0.45|1/250, 0.38|1/400, 0.27|1/1300
0684   inline double TAUJET_EFF_ATLAS_RUN1(const Jet& j) {
0685     if (j.abseta() > 2.5) return 0; //< hmm... mostly
0686     if (inRange(j.abseta(), 1.37, 1.52)) return 0; //< crack region
0687     double pThadvis = 0;
0688     Particles chargedhadrons;
0689     for (const Particle& p : j.particles()) {
0690       if (p.isHadron()) {
0691         pThadvis += p.pT(); //< right definition? Paper is unclear
0692         if (p.charge3() != 0 && p.abseta() < 2.5 && p.pT() > 1*GeV) chargedhadrons += p;
0693       }
0694     }
0695 
0696     if (chargedhadrons.empty()) return 0; //< leptonic tau?
0697     if (pThadvis < 20*GeV) return 0; //< below threshold
0698 
0699     const Particles ttags = j.tauTags(Cuts::pT > 10*GeV);
0700     // MisID rates for jets without truth tau label
0701     if (ttags.empty()) {
0702       if (pThadvis < 40*GeV)
0703         return chargedhadrons.size() == 1 ? 1/20. : chargedhadrons.size() == 3 ? 1/100. : 0; //< fake rates
0704       else
0705         return chargedhadrons.size() == 1 ? 1/25. : chargedhadrons.size() == 3 ? 1/400. : 0; //< fake rates
0706     }
0707     // Efficiencies for jets with a truth tau label
0708     const Particles prongs = ttags[0].stableDescendants(Cuts::charge3 > 0 && Cuts::pT > 1*GeV && Cuts::abseta < 2.5);
0709     return prongs.size() == 1 ? 0.56 : 0.38;
0710   }
0711 
0712 
0713   /// @brief ATLAS Run 2 13 TeV tau efficiencies (medium working point)
0714   ///
0715   /// From https://atlas.web.cern.ch/Atlas/GROUPS/PHYSICS/PUBNOTES/ATL-PHYS-PUB-2015-045/ATL-PHYS-PUB-2015-045.pdf
0716   ///   LMT 1 prong efficiency/mistag = 0.6|1/30, 0.55|1/50, 0.45|1/120
0717   ///   LMT 3 prong efficiency/mistag = 0.5|1/30, 0.4|1/110, 0.3|1/300
0718   inline double TAU_EFF_ATLAS_RUN2_MEDIUM(const Particle& t) {
0719     if (t.abspid() != PID::TAU) return 0;
0720     if (t.abseta() > 2.5) return 0; //< hmm... mostly
0721     if (inRange(t.abseta(), 1.37, 1.52)) return 0; //< crack region
0722     double pThadvis = 0;
0723     Particles chargedhadrons;
0724     for (const Particle& p : t.children()) {
0725       if (p.isHadron()) {
0726         pThadvis += p.pT(); //< right definition? Paper is unclear
0727         if (p.charge3() != 0 && p.abseta() < 2.5 && p.pT() > 1*GeV) chargedhadrons += p;
0728       }
0729     }
0730     if (chargedhadrons.empty()) return 0; //< leptonic tau
0731     if (pThadvis < 20*GeV) return 0; //< below threshold
0732     if (chargedhadrons.size() == 1) return (t.abspid() == PID::TAU) ? 0.55 : 0; //1/50.;
0733     if (chargedhadrons.size() == 3) return (t.abspid() == PID::TAU) ? 0.40 : 0; //1/110.;
0734     return 0;
0735   }
0736   /// ATLAS Run 2 medium muon reconstruction efficiency
0737   /// @todo Just an alias to generic: improve!
0738   inline double TAU_EFF_ATLAS_RUN2(const Particle& t) { return TAU_EFF_ATLAS_RUN2_MEDIUM(t); }
0739   /// ATLAS Run 2 loose tau reconstruction efficiency
0740   /// @todo Just an alias to generic: improve!
0741   inline double TAU_EFF_ATLAS_RUN2_LOOSE(const Particle& t) { return TAU_EFF_ATLAS_RUN2(t); }
0742   /// ATLAS Run 2 tight muon reconstruction efficiency
0743   /// @todo Just an alias to generic: improve!
0744   inline double TAU_EFF_ATLAS_RUN2_TIGHT(const Particle& t) { return TAU_EFF_ATLAS_RUN2(t); }
0745 
0746 
0747   /// @brief ATLAS Run 2 13 TeV tau misID rate (medium working point)
0748   ///
0749   /// From https://atlas.web.cern.ch/Atlas/GROUPS/PHYSICS/PUBNOTES/ATL-PHYS-PUB-2015-045/ATL-PHYS-PUB-2015-045.pdf
0750   ///   LMT 1 prong efficiency/mistag = 0.6|1/30, 0.55|1/50, 0.45|1/120
0751   ///   LMT 3 prong efficiency/mistag = 0.5|1/30, 0.4|1/110, 0.3|1/300
0752   inline double TAUJET_EFF_ATLAS_RUN2(const Jet& j) {
0753     if (j.abseta() > 2.5) return 0; //< hmm... mostly
0754     if (inRange(j.abseta(), 1.37, 1.52)) return 0; //< crack region
0755     double pThadvis = 0;
0756     Particles chargedhadrons;
0757     for (const Particle& p : j.particles()) {
0758       if (p.isHadron()) {
0759         pThadvis += p.pT(); //< right definition? Paper is unclear
0760         if (p.charge3() != 0 && p.abseta() < 2.5 && p.pT() > 1*GeV) chargedhadrons += p;
0761       }
0762     }
0763     if (chargedhadrons.empty()) return 0;
0764     if (pThadvis < 20*GeV) return 0; //< below threshold
0765     const Particles ttags = j.tauTags(Cuts::pT > 10*GeV);
0766     // if (ttags.empty()) {
0767     //   if (pThadvis < 40*GeV)
0768     //     return chargedhadrons.size() == 1 ? 1/50. : 1/110.; //< fake rates
0769     //   else
0770     //     return chargedhadrons.size() == 1 ? 1/25. : 1/400.; //< fake rates
0771     // }
0772     if (ttags.empty()) return chargedhadrons.size() == 1 ? 1/50. : chargedhadrons.size() == 3 ? 1/110. : 0; //< fake rates
0773     const Particles prongs = ttags[0].stableDescendants(Cuts::charge3 > 0 && Cuts::pT > 1*GeV && Cuts::abseta < 2.5);
0774     return prongs.size() == 1 ? 0.55 : 0.40;
0775   }
0776 
0777 
0778   /// ATLAS Run 1 tau smearing
0779   ///
0780   /// @todo Currently a copy of the jet smearing
0781   inline Particle TAU_SMEAR_ATLAS_RUN1(const Particle& t) {
0782     // // Const fractional resolution for now
0783     // static const double resolution = 0.03;
0784 
0785     // // Smear by a Gaussian centered on 1 with width given by the (fractional) resolution
0786     // /// @todo Is this the best way to smear? Should we preserve the energy, or pT, or direction?
0787     // const double fsmear = max(randnorm(1., resolution), 0.);
0788     // const double mass = t.mass2() > 0 ? t.mass() : 0; //< numerical carefulness...
0789     // return Particle(t.pid(), FourMomentum::mkXYZM(t.px()*fsmear, t.py()*fsmear, t.pz()*fsmear, mass));
0790 
0791     // Jet energy resolution lookup
0792     //   Implemented by Matthias Danninger for GAMBIT, based roughly on
0793     //   https://atlas.web.cern.ch/Atlas/GROUPS/PHYSICS/CONFNOTES/ATLAS-CONF-2015-017/
0794     //   Parameterisation can be still improved, but eta dependence is minimal
0795     /// @todo Also need a JES uncertainty component?
0796     static const vector<double> binedges_pt = {0., 50., 70., 100., 150., 200., 1000., 10000.};
0797     static const vector<double> jer = {0.145, 0.115, 0.095, 0.075, 0.07, 0.05, 0.04, 0.04}; //< note overflow value
0798     const int ipt = binIndex(t.pT()/GeV, binedges_pt, true);
0799     if (ipt < 0) return t;
0800     const double resolution = jer.at(ipt);
0801 
0802     // Smear by a Gaussian centered on 1 with width given by the (fractional) resolution
0803     /// @todo Is this the best way to smear? Should we preserve the energy, or pT, or direction?
0804     const double fsmear = max(randnorm(1., resolution), 0.);
0805     const double mass = t.mass2() > 0 ? t.mass() : 0; //< numerical carefulness...
0806     Particle rtn(PID::TAU, FourMomentum::mkXYZM(t.px()*fsmear, t.py()*fsmear, t.pz()*fsmear, mass));
0807     //if (deltaPhi(t, rtn) > 0.01) cout << "jdphi: " << deltaPhi(t, rtn) << endl;
0808     return rtn;
0809   }
0810 
0811 
0812   /// ATLAS Run 2 tau smearing
0813   ///
0814   /// @todo Currently a copy of the Run 1 version
0815   inline Particle TAU_SMEAR_ATLAS_RUN2(const Particle& t) {
0816     return TAU_SMEAR_ATLAS_RUN1(t);
0817   }
0818 
0819 
0820   /// CMS Run 1 tau efficiency
0821   ///
0822   /// @todo Needs work; this is just a copy of the Run 2 version in Delphes 3.3.2
0823   inline double TAU_EFF_CMS_RUN1(const Particle& t) {
0824     if (t.abspid() != PID::TAU) return 0;
0825     return (t.abspid() == PID::TAU) ? 0.6 : 0;
0826   }
0827   /// CMS Run 1 loose tau reconstruction efficiency
0828   /// @todo Just an alias to generic: improve!
0829   inline double TAU_EFF_CMS_RUN1_LOOSE(const Particle& t) { return TAU_EFF_CMS_RUN1(t); }
0830   /// CMS Run 1 medium muon reconstruction efficiency
0831   /// @todo Just an alias to generic: improve!
0832   inline double TAU_EFF_CMS_RUN1_MEDIUM(const Particle& t) { return TAU_EFF_CMS_RUN1(t); }
0833   /// CMS Run 1 tight muon reconstruction efficiency
0834   /// @todo Just an alias to generic: improve!
0835   inline double TAU_EFF_CMS_RUN1_TIGHT(const Particle& t) { return TAU_EFF_CMS_RUN1(t); }
0836 
0837   /// CMS Run 2 tau efficiency
0838   ///
0839   /// @todo Needs work; this is the dumb version from Delphes 3.3.2
0840   inline double TAU_EFF_CMS_RUN2(const Particle& t) {
0841     if (t.abspid() != PID::TAU) return 0;
0842     return (t.abspid() == PID::TAU) ? 0.6 : 0;
0843   }
0844   /// CMS Run 2 loose tau reconstruction efficiency
0845   /// @todo Just an alias to generic: improve!
0846   inline double TAU_EFF_CMS_RUN2_LOOSE(const Particle& t) { return TAU_EFF_CMS_RUN2(t); }
0847   /// CMS Run 2 medium muon reconstruction efficiency
0848   /// @todo Just an alias to generic: improve!
0849   inline double TAU_EFF_CMS_RUN2_MEDIUM(const Particle& t) { return TAU_EFF_CMS_RUN2(t); }
0850   /// CMS Run 2 tight muon reconstruction efficiency
0851   /// @todo Just an alias to generic: improve!
0852   inline double TAU_EFF_CMS_RUN2_TIGHT(const Particle& t) { return TAU_EFF_CMS_RUN2(t); }
0853 
0854 
0855   /// CMS Run 1 tau smearing
0856   ///
0857   /// @todo Currently a copy of the crappy ATLAS one
0858   inline Particle TAU_SMEAR_CMS_RUN1(const Particle& t) {
0859     return TAU_SMEAR_ATLAS_RUN1(t);
0860   }
0861 
0862 
0863   /// CMS Run 2 tau smearing
0864   ///
0865   /// @todo Currently a copy of the Run 1 version
0866   inline Particle TAU_SMEAR_CMS_RUN2(const Particle& t) {
0867     return TAU_SMEAR_CMS_RUN1(t);
0868   }
0869 
0870   /// @}
0871 
0872 
0873 
0874   /// @defgroup smearing_jet Experiment-specific jet efficiency and smearing functions
0875   /// @{
0876 
0877   /// Return the ATLAS Run 1 jet flavour tagging efficiency for the given Jet, from Delphes
0878   inline double JET_BTAG_ATLAS_RUN1(const Jet& j) {
0879     /// @todo This form drops past ~100 GeV, asymptotically to zero efficiency... really?!
0880     if (j.abseta() > 2.5) return 0;
0881     const auto ftagsel = [&](const Particle& p){ return p.pT() > 5*GeV && deltaR(p,j) < 0.3; };
0882     if (j.bTagged(ftagsel)) return 0.80*tanh(0.003*j.pT()/GeV)*(30/(1+0.0860*j.pT()/GeV));
0883     if (j.cTagged(ftagsel)) return 0.20*tanh(0.020*j.pT()/GeV)*( 1/(1+0.0034*j.pT()/GeV));
0884     return 0.002 + 7.3e-6*j.pT()/GeV;
0885   }
0886   /// Alias for naming scheme
0887   inline double JET_BTAG_ATLAS_RUN1_XXX(const Jet& j) { return JET_BTAG_ATLAS_RUN1(j); }
0888 
0889   /// Return the ATLAS Run 2 MC2c20 77% WP jet flavour tagging efficiency for the given Jet
0890   inline double JET_BTAG_ATLAS_RUN2_MV2C20(const Jet& j) {
0891     if (j.abseta() > 2.5) return 0;
0892     if (j.bTagged(Cuts::pT > 5*GeV)) return 0.77;
0893     if (j.cTagged(Cuts::pT > 5*GeV)) return 1/4.5;
0894     return 1/140.;
0895   }
0896 
0897   /// Return the ATLAS Run 2 MC2c10 77% WP jet flavour tagging efficiency for the given Jet
0898   inline double JET_BTAG_ATLAS_RUN2_MV2C10(const Jet& j) {
0899     if (j.abseta() > 2.5) return 0;
0900     if (j.bTagged(Cuts::pT > 5*GeV)) return 0.77;
0901     if (j.cTagged(Cuts::pT > 5*GeV)) return 1/6.0;
0902     return 1/134.;
0903   }
0904 
0905 
0906   /// @brief ATLAS Run 1 jet smearing
0907   ///
0908   /// Implemented by Matthias Danninger for GAMBIT, based roughly on
0909   /// https://atlas.web.cern.ch/Atlas/GROUPS/PHYSICS/CONFNOTES/ATLAS-CONF-2015-017/
0910   inline Jet JET_SMEAR_ATLAS_RUN1(const Jet& j) {
0911     // Jet energy resolution lookup
0912     //   Parameterisation can be still improved, but eta dependence is minimal
0913     /// @todo Also need a JES uncertainty component?
0914     static const vector<double> binedges_pt = {0., 50., 70., 100., 150., 200., 1000., 10000.};
0915     static const vector<double> jer = {0.145, 0.115, 0.095, 0.075, 0.07, 0.05, 0.04, 0.04}; //< note overflow value
0916     const int ipt = binIndex(j.pT()/GeV, binedges_pt, true);
0917     if (ipt < 0) return j;
0918     const double resolution = jer.at(ipt);
0919 
0920     // Smear by a Gaussian centered on 1 with width given by the (fractional) resolution
0921     /// @todo Is this the best way to smear? Should we preserve the energy, or pT, or direction?
0922     const double fsmear = max(randnorm(1., resolution), 0.);
0923     const double mass = j.mass2() > 0 ? j.mass() : 0; //< numerical carefulness...
0924     Jet rtn(FourMomentum::mkXYZM(j.px()*fsmear, j.py()*fsmear, j.pz()*fsmear, mass));
0925     //if (deltaPhi(j, rtn) > 0.01) cout << "jdphi: " << deltaPhi(j, rtn) << endl;
0926     return rtn;
0927   }
0928 
0929   /// ATLAS Run 2 jet smearing
0930   ///
0931   /// @todo Just a copy of the Run 1 one: improve!!
0932   inline Jet JET_SMEAR_ATLAS_RUN2(const Jet& j) { return JET_SMEAR_ATLAS_RUN1(j); }
0933 
0934   /// CMS Run 2 jet smearing
0935   ///
0936   /// @todo Just a copy of the suboptimal ATLAS one: improve!!
0937   inline Jet JET_SMEAR_CMS_RUN1(const Jet& j) { return JET_SMEAR_ATLAS_RUN1(j); }
0938 
0939   /// CMS Run 2 jet smearing
0940   ///
0941   /// @todo Just a copy of the suboptimal ATLAS one: improve!!
0942   inline Jet JET_SMEAR_CMS_RUN2(const Jet& j) { return JET_SMEAR_CMS_RUN1(j); }
0943 
0944   /// @}
0945 
0946 
0947   /// @defgroup smearing_met Experiment-specific missing-ET smearing functions
0948   /// @{
0949 
0950   /// @brief ATLAS Run 1 ETmiss resolution
0951   ///
0952   /// Based on https://arxiv.org/pdf/1108.5602v2.pdf, Figs 14 and 15
0953   inline METSmearParams MET_SMEARPARAMS_ATLAS_RUN1(const Vector3& met, double set) {
0954     // Linearity offset (Fig 14)
0955     Vector3 smeared_met = met;
0956     if (met.mod() < 25*GeV) smeared_met *= 1.05;
0957     else if (met.mod() < 40*GeV) smeared_met *= (1.05 - (0.04/15)*(met.mod()/GeV - 25)); //< linear decrease
0958     else smeared_met *= 1.01;
0959 
0960     return {smeared_met, 0.45*sqrt(set/GeV)*GeV, 0.0};
0961   }
0962 
0963   /// @brief ATLAS Run 1 ETmiss smearing
0964   inline Vector3 MET_SMEAR_ATLAS_RUN1(const Vector3& met, double set) {
0965     return MET_SMEAR_NORM(MET_SMEARPARAMS_ATLAS_RUN1(met, set));
0966   }
0967 
0968 
0969   /// ATLAS Run 2 ETmiss resolution
0970   ///
0971   /// Based on https://arxiv.org/pdf/1802.08168.pdf, Figs 6-9
0972   inline METSmearParams MET_SMEARPARAMS_ATLAS_RUN2(const Vector3& met, double set) {
0973     // Linearity offset (Fig 6)
0974     Vector3 smeared_met = met;
0975     if (met.mod() < 25*GeV) smeared_met *= 1.5;
0976     else smeared_met *= (1 + exp(-(met.mod() - 25*GeV)/(10*GeV)) - 0.02); //< exp approx to Fig 6 curve, approaching -0.02
0977 
0978     // Resolution(sumEt) ~ 0.45 sqrt(sumEt) GeV
0979     // above SET = 180 GeV, and with linear trend from SET = 180 -> 0  to resolution = 0 (Fig 7)
0980     const double resolution1 = (set < 180*GeV ? set/180. : 1) * 0.45 * sqrt(max(set/GeV, 180)) * GeV;
0981     /// @todo Allow smearing function to access the whole event, since Njet also affects? Or assume encoded in SET?
0982 
0983     // Resolution(MET_true) (Fig 9)
0984     const double resolution2 = 15*GeV + 0.5*sqrt(met.mod()/GeV)*GeV;
0985 
0986     // Smearing width is given by the minimum resolution estimator (should mean low-MET events
0987     // with high SET do not get a large smearing, and will be dominated by the linearity effect).
0988     /// @todo Arguably should somehow take the max of this and the linearity... but this is re-used
0989     const double sigma = min(resolution1, resolution2);
0990 
0991     return {smeared_met, sigma, 0.0};
0992   }
0993 
0994   /// ATLAS Run 2 ETmiss smearing
0995   inline Vector3 MET_SMEAR_ATLAS_RUN2(const Vector3& met, double set) {
0996     return MET_SMEAR_NORM(MET_SMEARPARAMS_ATLAS_RUN2(met, set));
0997   }
0998 
0999 
1000   /// CMS Run 1 ETmiss smearing params
1001   ///
1002   /// From https://arxiv.org/pdf/1411.0511.pdf Table 2, p16 (Z channels)
1003   inline METSmearParams MET_SMEARPARAMS_CMS_RUN1(const Vector3& met, double set) {
1004     // Calculate parallel and perpendicular resolutions and combine in quadrature (?)
1005     const double resolution_x = (1.1 + 0.6*sqrt(set/GeV)) * GeV;
1006     const double resolution_y = (1.4 + 0.6*sqrt(set/GeV)) * GeV;
1007     const double resolution = sqrt(sqr(resolution_x) + sqr(resolution_y));
1008     return {met, resolution,0.0};
1009   }
1010 
1011   /// CMS Run 1 ETmiss smearing
1012   inline Vector3 MET_SMEAR_CMS_RUN1(const Vector3& met, double set) {
1013     return MET_SMEAR_NORM(MET_SMEARPARAMS_CMS_RUN1(met, set));
1014   }
1015 
1016 
1017   /// @brief CMS Run 2 ETmiss smearing
1018   ///
1019   /// From http://inspirehep.net/record/1681214/files/JME-17-001-pas.pdf Table 3, p20
1020   inline METSmearParams MET_SMEARPARAMS_CMS_RUN2(const Vector3& met, double set) {
1021     // Calculate parallel and perpendicular resolutions and combine in quadrature (?)
1022     const double resolution_para = ( 2.0 + 0.64*sqrt(set/GeV)) * GeV;
1023     const double resolution_perp = (-1.5 + 0.68*sqrt(set/GeV)) * GeV;
1024     const double resolution = sqrt(sqr(resolution_para) + sqr(resolution_perp));
1025     return {met, resolution,0.0};
1026   }
1027 
1028   /// @brief CMS Run 2 ETmiss smearing
1029   inline Vector3 MET_SMEAR_CMS_RUN2(const Vector3& met, double set) {
1030     return MET_SMEAR_NORM(MET_SMEARPARAMS_CMS_RUN2(met, set));
1031   }
1032 
1033   /// @}
1034 
1035 
1036   /// @defgroup smearing_trk Experiment-specific tracking efficiency and smearing functions
1037   /// @{
1038 
1039   /// ATLAS Run 1 tracking efficiency
1040   inline double TRK_EFF_ATLAS_RUN1(const Particle& p) {
1041     if (p.charge3() == 0) return 0;
1042     if (p.abseta() > 2.5) return 0;
1043     if (p.pT() < 0.1*GeV) return 0;
1044 
1045     if (p.abspid() == PID::ELECTRON) {
1046       if (p.abseta() < 1.5) {
1047         if (p.pT() < 1*GeV) return 0.73;
1048         if (p.pT() < 100*GeV) return 0.95;
1049         return 0.99;
1050       } else {
1051         if (p.pT() < 1*GeV) return 0.50;
1052         if (p.pT() < 100*GeV) return 0.83;
1053         else return 0.90;
1054       }
1055     } else if (p.abspid() == PID::MUON) {
1056       if (p.abseta() < 1.5) {
1057         return (p.pT() < 1*GeV) ? 0.75 : 0.99;
1058       } else {
1059         return (p.pT() < 1*GeV) ? 0.70 : 0.98;
1060       }
1061     } else { // charged hadrons
1062       if (p.abseta() < 1.5) {
1063         return (p.pT() < 1*GeV) ? 0.70 : 0.95;
1064       } else {
1065         return (p.pT() < 1*GeV) ? 0.60 : 0.85;
1066       }
1067     }
1068   }
1069 
1070   /// ATLAS Run 2 tracking efficiency
1071   ///
1072   /// @todo Currently just a copy of Run 1: fix!
1073   inline double TRK_EFF_ATLAS_RUN2(const Particle& p) {
1074     return TRK_EFF_ATLAS_RUN1(p);
1075   }
1076 
1077 
1078   /// ATLAS Run 1 track smearing
1079   ///
1080   /// @todo Currently identity: fix!
1081   inline Particle TRK_SMEAR_ATLAS_RUN1(const Particle& t) {
1082     return PARTICLE_SMEAR_IDENTITY(t);
1083   }
1084   /// ATLAS Run 2 track smearing
1085   ///
1086   /// @todo Currently identity: fix!
1087   inline Particle TRK_SMEAR_ATLAS_RUN2(const Particle& t) {
1088     return PARTICLE_SMEAR_IDENTITY(t);
1089   }
1090 
1091 
1092   /// CMS Run 1 tracking efficiency
1093   inline double TRK_EFF_CMS_RUN1(const Particle& p) {
1094     if (p.charge3() == 0) return 0;
1095     if (p.abseta() > 2.5) return 0;
1096     if (p.pT() < 0.1*GeV) return 0;
1097 
1098     if (p.abspid() == PID::ELECTRON) {
1099       if (p.abseta() < 1.5) {
1100         if (p.pT() < 1*GeV) return 0.73;
1101         if (p.pT() < 100*GeV) return 0.95;
1102         return 0.99;
1103       } else {
1104         if (p.pT() < 1*GeV) return 0.50;
1105         if (p.pT() < 100*GeV) return 0.83;
1106         else return 0.90;
1107       }
1108     } else if (p.abspid() == PID::MUON) {
1109       if (p.abseta() < 1.5) {
1110         return (p.pT() < 1*GeV) ? 0.75 : 0.99;
1111       } else {
1112         return (p.pT() < 1*GeV) ? 0.70 : 0.98;
1113       }
1114     } else { // charged hadrons
1115       if (p.abseta() < 1.5) {
1116         return (p.pT() < 1*GeV) ? 0.70 : 0.95;
1117       } else {
1118         return (p.pT() < 1*GeV) ? 0.60 : 0.85;
1119       }
1120     }
1121   }
1122 
1123   /// CMS Run 2 tracking efficiency
1124   ///
1125   /// @todo Currently just a copy of Run 1: fix!
1126   inline double TRK_EFF_CMS_RUN2(const Particle& p) {
1127     return TRK_EFF_CMS_RUN1(p);
1128   }
1129 
1130 
1131   /// CMS Run 1 track smearing
1132   ///
1133   /// @todo Currently identity: fix!
1134   inline Particle TRK_SMEAR_CMS_RUN1(const Particle& t) {
1135     return PARTICLE_SMEAR_IDENTITY(t);
1136   }
1137   /// CMS Run 2 track smearing
1138   ///
1139   /// @todo Currently identity: fix!
1140   inline Particle TRK_SMEAR_CMS_RUN2(const Particle& t) {
1141     return PARTICLE_SMEAR_IDENTITY(t);
1142   }
1143 
1144   /// @brief Return the efficiency of the ATLAS JVT tagger at > 0.2 W.P.
1145   ///
1146   /// Values taken from fig 17, https://arxiv.org/pdf/1510.03823.pdf.
1147   /// Shockingly, its run1, but it's still the most recent ATLAS JVT public plots
1148   /// And is still cited by all the Run 2 papers
1149   /// Note! The exact translation of weak/medium/tight varies across runs/years
1150   /// If possible double check the score cut off used.
1151   inline double ATLAS_JVT_EFF_LOOSE(const Jet & j){
1152     if (j.abseta() > 2.4) return 1.0;
1153     // No data for jets pt < 20, hopefully not relevant:
1154     if (j.pt() <= 20*GeV || j.pt() >= 60*GeV) return 1.0;
1155 
1156     const static vector<double> binedges_pt = {20*GeV,25*GeV,30*GeV,35*GeV,40*GeV,45*GeV,50*GeV,55*GeV,60*GeV};
1157     const static vector<double> binvals = {0.9, 0.93, 0.94, 0.95, 0.96, 0.96, 0.97, 0.97};
1158 
1159     const size_t bini = binIndex(j.pt(), binedges_pt);
1160     return binvals[bini];
1161   }
1162 
1163   /// @brief Return the efficiency of the ATLAS JVT tagger at > 0.4 W.P.
1164   ///
1165   /// Values taken from fig 17, https://arxiv.org/pdf/1510.03823.pdf.
1166   /// Shockingly, its run1, but it's still the most recent ATLAS JVT public plots
1167   /// And is still cited by all the Run 2 papers
1168   /// Note! The exact translation of weak/medium/tight varies across runs/years
1169   /// If possible double check the score cut off used.
1170   inline double ATLAS_JVT_EFF_MEDIUM(const Jet & j){
1171     if (j.abseta() > 2.4) return 1.0;
1172     // No data for jets pt < 20, hopefully not relevant:
1173     if (j.pt() <= 20*GeV || j.pt() >= 60*GeV) return 1.0;
1174 
1175     const static vector<double> binedges_pt = {20*GeV,25*GeV,30*GeV,35*GeV,40*GeV,45*GeV,50*GeV,55*GeV,60*GeV};
1176     const static vector<double> binvals = {0.86, 0.9, 0.92, 0.93, 0.94, 0.95, 0.95, 0.96};
1177 
1178     const size_t bini = binIndex(j.pt(), binedges_pt);
1179     return binvals[bini];
1180   }
1181 
1182   /// @brief Return the efficiency of the ATLAS JVT tagger at > 0.7 W.P.
1183   ///
1184   /// Values taken from fig 17, https://arxiv.org/pdf/1510.03823.pdf.
1185   /// Shockingly, its run1, but it's still the most recent ATLAS JVT public plots
1186   /// And is still cited by all the Run 2 papers
1187   /// Note! The exact translation of weak/medium/tight varies across runs/years
1188   /// If possible double check the score cut off used.
1189   inline double ATLAS_JVT_EFF_TIGHT(const Jet & j){
1190     if (j.abseta() > 2.4) return 1.0;
1191     // No data for jets pt < 20, hopefully not relevant:
1192     if (j.pt() <= 20*GeV || j.pt() >= 60*GeV) return 1.0;
1193 
1194     const static vector<double> binedges_pt = {20*GeV,25*GeV,30*GeV,35*GeV,40*GeV,45*GeV,50*GeV,55*GeV,60*GeV};
1195     const static vector<double> binvals = {0.81, 0.84, 0.87, 0.90, 0.91, 0.92, 0.93, 0.94};
1196 
1197     const size_t bini = binIndex(j.pt(), binedges_pt);
1198     return binvals[bini];
1199   }
1200 
1201   /// @}
1202 
1203   /// @name A full family of identity smearing functions
1204   /// @{
1205 
1206   inline double ELECTRON_EFF_IDENTITY_LOOSE(const Particle& e) { return PARTICLE_EFF_ONE(e); }
1207   inline double ELECTRON_EFF_IDENTITY_MEDIUM(const Particle& e) { return PARTICLE_EFF_ONE(e); }
1208   inline double ELECTRON_EFF_IDENTITY_TIGHT(const Particle& e) { return PARTICLE_EFF_ONE(e); }
1209   inline double MUON_EFF_IDENTITY_LOOSE(const Particle& m) { return PARTICLE_EFF_ONE(m); }
1210   inline double MUON_EFF_IDENTITY_MEDIUM(const Particle& m) { return PARTICLE_EFF_ONE(m); }
1211   inline double MUON_EFF_IDENTITY_TIGHT(const Particle& m) { return PARTICLE_EFF_ONE(m); }
1212   inline double PHOTON_EFF_IDENTITY_LOOSE(const Particle& y) { return PARTICLE_EFF_ONE(y); }
1213   inline double PHOTON_EFF_IDENTITY_MEDIUM(const Particle& y) { return PARTICLE_EFF_ONE(y); }
1214   inline double PHOTON_EFF_IDENTITY_TIGHT(const Particle& y) { return PARTICLE_EFF_ONE(y); }
1215   inline double TAU_EFF_IDENTITY_LOOSE(const Particle& t) { return PARTICLE_EFF_ONE(t); }
1216   inline double TAU_EFF_IDENTITY_MEDIUM(const Particle& t) { return PARTICLE_EFF_ONE(t); }
1217   inline double TAU_EFF_IDENTITY_TIGHT(const Particle& t) { return PARTICLE_EFF_ONE(t); }
1218 
1219   inline Particle ELECTRON_SMEAR_IDENTITY(const Particle& e) { return PARTICLE_SMEAR_IDENTITY(e); }
1220   inline Particle MUON_SMEAR_IDENTITY(const Particle& m) { return PARTICLE_SMEAR_IDENTITY(m); }
1221   inline Particle PHOTON_SMEAR_IDENTITY(const Particle& y) { return PARTICLE_SMEAR_IDENTITY(y); }
1222   inline Particle TAU_SMEAR_IDENTITY(const Particle& t) { return PARTICLE_SMEAR_IDENTITY(t); }
1223 
1224   inline double TRK_EFF_IDENTITY_TIGHT(const Particle& trk) { return PARTICLE_EFF_ONE(trk); }
1225   inline Particle TRK_SMEAR_IDENTITY(const Particle& trk) { return PARTICLE_SMEAR_IDENTITY(trk); }
1226 
1227   inline double JET_BTAG_IDENTITY_IDENTITY(const Jet& j) { return JET_BTAG_IDENTITY(j); }
1228 
1229   // Already defined: JET_SMEAR_IDENTITY and MET_SMEAR_IDENTITY
1230 
1231   /// @}
1232 
1233   /// @}
1234 
1235 
1236 }
1237 
1238 #endif