rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

ATLAS_2014_I1327229

Inclusive multilepton search at 8 TeV
Experiment: ATLAS (LHC)
Inspire ID: 1327229
Status: VALIDATED
Authors:
  • Joern Mahlstedt
References: Beams: p+ p+
Beam energies: (4000.0, 4000.0) GeV
Run details:
  • Any process producing at least 3 leptons (e.g. pair production of doubly-charged Higgs or excited leptons)

A generic search for anomalous production of events with at least three charged leptons is presented. The data sample consists of $pp$ collisions at $\sqrt{s} = 8$\,TeV collected in 2012 by the ATLAS experiment at the CERN Large Hadron Collider, and corresponds to an integrated luminosity of 20.3\,$\text{fb}^{-1}$. Events are required to have at least three selected lepton candidates, at least two of which must be electrons or muons, while the third may be a hadronically decaying tau. Selected events are categorized based on their lepton flavour content and signal regions are constructed using several kinematic variables of interest. No significant deviations from Standard Model predictions are observed. Model-independent upper limits on contributions from beyond the Standard Model phenomena are provided for each signal region, along with prescription to re-interpret the limits for any model. Constraints are also placed on models predicting doubly charged Higgs bosons and excited leptons. For doubly charged Higgs bosons decaying to $e\tau$ or $\muon\tau$, lower limits on the mass are set at 400\,GeV at 95\,\% confidence level. For excited leptons, constraints are provided as functions of both the mass of the excited state and the compositeness scale $\Lambda$, with the strongest mass constraints arising in regions where the mass equals $\Lambda$. In such scenarios, lower mass limits are set at 3.0\,TeV for excited electrons and muons, 2.5\,TeV for excited taus, and 1.6\,TeV for every excited-neutrino flavour.

Source code: ATLAS_2014_I1327229.cc
   1// -*- C++ -*-
   2#include "Rivet/Analysis.hh"
   3#include "Rivet/Projections/FinalState.hh"
   4#include "Rivet/Projections/ChargedFinalState.hh"
   5#include "Rivet/Projections/VisibleFinalState.hh"
   6#include "Rivet/Projections/VetoedFinalState.hh"
   7#include "Rivet/Projections/IdentifiedFinalState.hh"
   8#include "Rivet/Projections/UnstableParticles.hh"
   9#include "Rivet/Projections/FastJets.hh"
  10#include "Rivet/Tools/Random.hh"
  11
  12namespace Rivet {
  13
  14
  15  /// Inclusive multilepton search at 8 TeV
  16  class ATLAS_2014_I1327229 : public Analysis {
  17  public:
  18
  19    /// Constructor
  20    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2014_I1327229);
  21
  22
  23    /// Book histograms and initialise projections before the run
  24    void init() {
  25
  26      // To calculate the acceptance without having the fiducial lepton efficiencies included, this part can be turned off
  27      _use_fiducial_lepton_efficiency = true;
  28
  29      // Read in all signal regions
  30      _signal_regions = getSignalRegions();
  31
  32      // Set number of events per signal region to 0
  33      for (size_t i = 0; i < _signal_regions.size(); i++)
  34        book(_eventCountsPerSR[_signal_regions[i]], "_eventCountsPerSR_" + _signal_regions[i]);
  35
  36      // Final state including all charged and neutral particles
  37      const FinalState fs((Cuts::etaIn(-5.0, 5.0) && Cuts::pT >=  1*GeV));
  38      declare(fs, "FS");
  39
  40      // Final state including all charged particles
  41      declare(ChargedFinalState((Cuts::etaIn(-2.5, 2.5) && Cuts::pT >=  1*GeV)), "CFS");
  42
  43      // Final state including all visible particles (to calculate MET, Jets etc.)
  44      declare(VisibleFinalState((Cuts::etaIn(-5.0,5.0))),"VFS");
  45
  46      // Final state including all AntiKt 04 Jets
  47      VetoedFinalState vfs;
  48      vfs.addVetoPairId(PID::MUON);
  49      declare(FastJets(vfs, JetAlg::ANTIKT, 0.4), "AntiKtJets04");
  50
  51      // Final state including all unstable particles (including taus)
  52      declare(UnstableParticles(Cuts::abseta < 5.0 && Cuts::pT > 5*GeV),"UFS");
  53
  54      // Final state including all electrons
  55      IdentifiedFinalState elecs(Cuts::abseta < 2.47 && Cuts::pT > 10*GeV);
  56      elecs.acceptIdPair(PID::ELECTRON);
  57      declare(elecs, "elecs");
  58
  59      // Final state including all muons
  60      IdentifiedFinalState muons(Cuts::abseta < 2.5 && Cuts::pT > 10*GeV);
  61      muons.acceptIdPair(PID::MUON);
  62      declare(muons, "muons");
  63
  64
  65
  66      /// Book histograms:
  67      book(_h_HTlep_all ,"HTlep_all", 30,0,3000);
  68      book(_h_HTjets_all ,"HTjets_all", 30,0,3000);
  69      book(_h_MET_all ,"MET_all", 30,0,1500);
  70      book(_h_Meff_all ,"Meff_all", 50,0,5000);
  71      book(_h_min_pT_all ,"min_pT_all", 50, 0, 2000);
  72      book(_h_mT_all ,"mT_all", 50, 0, 2000);
  73
  74      book(_h_e_n ,"e_n", 10, -0.5, 9.5);
  75      book(_h_mu_n ,"mu_n", 10, -0.5, 9.5);
  76      book(_h_tau_n ,"tau_n", 10, -0.5, 9.5);
  77
  78      book(_h_pt_1_3l ,"pt_1_3l", 100, 0, 2000);
  79      book(_h_pt_2_3l ,"pt_2_3l", 100, 0, 2000);
  80      book(_h_pt_3_3l ,"pt_3_3l", 100, 0, 2000);
  81      book(_h_pt_1_2ltau ,"pt_1_2ltau", 100, 0, 2000);
  82      book(_h_pt_2_2ltau ,"pt_2_2ltau", 100, 0, 2000);
  83      book(_h_pt_3_2ltau ,"pt_3_2ltau", 100, 0, 2000);
  84
  85      book(_h_excluded ,"excluded", 2, -0.5, 1.5);
  86    }
  87
  88
  89    /// Perform the per-event analysis
  90    void analyze(const Event& event) {
  91
  92      // Muons
  93      Particles muon_candidates;
  94      const Particles charged_tracks = apply<ChargedFinalState>(event, "CFS").particles();
  95      const Particles visible_particles = apply<VisibleFinalState>(event, "VFS").particles();
  96      for (const Particle& mu : apply<IdentifiedFinalState>(event, "muons").particlesByPt() ) {
  97
  98        // Calculate pTCone30 variable (pT of all tracks within dR<0.3 - pT of muon itself)
  99        double pTinCone = -mu.pT();
 100        for (const Particle& track : charged_tracks ) {
 101          if (deltaR(mu.momentum(),track.momentum()) < 0.3 )
 102            pTinCone += track.pT();
 103        }
 104
 105        // Calculate eTCone30 variable (pT of all visible particles within dR<0.3)
 106        double eTinCone = 0.;
 107        for (const Particle& visible_particle : visible_particles) {
 108          if (visible_particle.abspid() != PID::MUON && inRange(deltaR(mu.momentum(),visible_particle.momentum()), 0.1, 0.3))
 109            eTinCone += visible_particle.pT();
 110        }
 111
 112        // Apply reconstruction efficiency and simulate reconstruction
 113        int muon_id = 13;
 114        if (mu.hasAncestorWith(Cuts::pid == PID::TAU) || mu.hasAncestorWith(Cuts::pid == -PID::TAU)) muon_id = 14;
 115        const double eff = (_use_fiducial_lepton_efficiency) ? apply_reco_eff(muon_id,mu) : 1.0;
 116        const bool keep_muon = rand01()<=eff;
 117
 118        // Keep muon if pTCone30/pT < 0.15 and eTCone30/pT < 0.2 and reconstructed
 119        if (keep_muon && pTinCone/mu.pT() <= 0.1 && eTinCone/mu.pT() < 0.1)
 120          muon_candidates.push_back(mu);
 121      }
 122
 123      // Electrons
 124      Particles electron_candidates;
 125      for (const Particle& e : apply<IdentifiedFinalState>(event, "elecs").particlesByPt() ) {
 126        // Neglect electrons in crack regions
 127        if (inRange(e.abseta(), 1.37, 1.52)) continue;
 128
 129        // Calculate pTCone30 variable (pT of all tracks within dR<0.3 - pT of electron itself)
 130        double pTinCone = -e.pT();
 131        for (const Particle& track : charged_tracks) {
 132          if (deltaR(e.momentum(), track.momentum()) < 0.3 ) pTinCone += track.pT();
 133        }
 134
 135        // Calculate eTCone30 variable (pT of all visible particles (except muons) within dR<0.3)
 136        double eTinCone = 0.;
 137        for (const Particle& visible_particle : visible_particles) {
 138          if (visible_particle.abspid() != PID::MUON && inRange(deltaR(e.momentum(),visible_particle.momentum()), 0.1, 0.3))
 139            eTinCone += visible_particle.pT();
 140        }
 141
 142        // Apply reconstruction efficiency and simulate reconstruction
 143        int elec_id = 11;
 144        if (e.hasAncestorWith(Cuts::pid == 15) || e.hasAncestorWith(Cuts::pid == -15)) elec_id = 12;
 145        const double eff = (_use_fiducial_lepton_efficiency) ? apply_reco_eff(elec_id,e) : 1.0;
 146        const bool keep_elec = rand01()<=eff;
 147
 148        // Keep electron if pTCone30/pT < 0.13 and eTCone30/pT < 0.2 and reconstructed
 149        if (keep_elec && pTinCone/e.pT() <= 0.1  && eTinCone/e.pT() < 0.1)
 150          electron_candidates.push_back(e);
 151      }
 152
 153
 154      // Taus
 155      Particles tau_candidates;
 156      for (const Particle& tau : apply<UnstableParticles>(event, "UFS").particles() ) {
 157        // Only pick taus out of all unstable particles
 158        if ( tau.abspid() != PID::TAU) continue;
 159        // Check that tau has decayed into daughter particles
 160        if (tau.genParticle()->end_vertex() == 0) continue;
 161        // Calculate visible tau momentum using the tau neutrino momentum in the tau decay
 162        FourMomentum daughter_tau_neutrino_momentum = get_tau_neutrino_momentum(tau);
 163        Particle tau_vis = tau;
 164        tau_vis.setMomentum(tau.momentum()-daughter_tau_neutrino_momentum);
 165        // keep only taus in certain eta region and above 15 GeV of visible tau pT
 166        if ( tau_vis.pT()/GeV <= 15.0 || tau_vis.abseta() > 2.5) continue;
 167
 168        // Get prong number (number of tracks) in tau decay and check if tau decays leptonically
 169        unsigned int nprong = 0;
 170        bool lep_decaying_tau = false;
 171        get_prong_number(tau.genParticle(),nprong,lep_decaying_tau);
 172
 173        // Apply reconstruction efficiency and simulate reconstruction
 174        int tau_id = 15;
 175        if (nprong == 1) tau_id = 15;
 176        else if (nprong == 3) tau_id = 16;
 177
 178
 179        const double eff = (_use_fiducial_lepton_efficiency) ? apply_reco_eff(tau_id,tau_vis) : 1.0;
 180        const bool keep_tau = rand01()<=eff;
 181
 182        // Keep tau if nprong = 1, it decays hadronically and it is reconstructed
 183        if ( !lep_decaying_tau && nprong == 1 && keep_tau) tau_candidates.push_back(tau_vis);
 184      }
 185
 186      // Jets (all anti-kt R=0.4 jets with pT > 30 GeV and eta < 4.9
 187      Jets jet_candidates = apply<FastJets>(event, "AntiKtJets04").jetsByPt(Cuts::pT > 30.0*GeV && Cuts::abseta < 4.9);
 188
 189      // ETmiss
 190      Particles vfs_particles = apply<VisibleFinalState>(event, "VFS").particles();
 191      FourMomentum pTmiss;
 192      for (const Particle& p : vfs_particles)
 193        pTmiss -= p.momentum();
 194      double eTmiss = pTmiss.pT()/GeV;
 195
 196
 197      // -------------------------
 198      // Overlap removal
 199
 200      // electron - electron
 201      Particles electron_candidates_2;
 202      for(size_t ie = 0; ie < electron_candidates.size(); ++ie) {
 203        const Particle& e = electron_candidates[ie];
 204        bool away = true;
 205        // If electron pair within dR < 0.1: remove electron with lower pT
 206        for(size_t ie2 = 0; ie2 < electron_candidates_2.size(); ++ie2) {
 207          if (deltaR(e.momentum(),electron_candidates_2[ie2].momentum()) < 0.1 ) {
 208            away = false;
 209            break;
 210          }
 211        }
 212        // If isolated keep it
 213        if ( away )
 214          electron_candidates_2.push_back( e );
 215      }
 216
 217      // jet - electron
 218      Jets recon_jets;
 219      for (const Jet& jet : jet_candidates) {
 220        bool away = true;
 221        // If jet within dR < 0.2 of electron: remove jet
 222        for (const Particle& e : electron_candidates_2) {
 223          if (deltaR(e.momentum(), jet.momentum()) < 0.2 ) {
 224            away = false;
 225            break;
 226          }
 227        }
 228
 229        // jet - tau
 230        if ( away )  {
 231          // If jet within dR < 0.2 of tau: remove jet
 232          for (const Particle& tau : tau_candidates) {
 233            if (deltaR(tau.momentum(), jet.momentum()) < 0.2 ) {
 234              away = false;
 235              break;
 236            }
 237          }
 238        }
 239        // If isolated keep it
 240        if ( away )
 241          recon_jets.push_back( jet );
 242      }
 243
 244      // electron - jet
 245      Particles recon_leptons, recon_e;
 246      for (size_t ie = 0; ie < electron_candidates_2.size(); ++ie) {
 247        const Particle& e = electron_candidates_2[ie];
 248        // If electron within 0.2 < dR < 0.4 from any jets: remove electron
 249        bool away = true;
 250        for (const Jet& jet : recon_jets) {
 251          if (deltaR(e.momentum(), jet.momentum()) < 0.4 ) {
 252            away = false;
 253            break;
 254          }
 255        }
 256        // electron - muon
 257        // If electron within dR < 0.1 of a muon: remove electron
 258        if (away) {
 259          for (const Particle& mu : muon_candidates) {
 260            if (deltaR(mu.momentum(),e.momentum()) < 0.1) {
 261              away = false;
 262              break;
 263            }
 264          }
 265        }
 266        // If isolated keep it
 267        if ( away )  {
 268          recon_e.push_back( e );
 269          recon_leptons.push_back( e );
 270          }
 271      }
 272
 273      // tau - electron
 274      Particles recon_tau;
 275      for (const Particle& tau : tau_candidates) {
 276        bool away = true;
 277        // If tau within dR < 0.2 of an electron: remove tau
 278        for (const Particle & e : recon_e) {
 279          if (deltaR(tau.momentum(),e.momentum()) < 0.2 ) {
 280            away = false;
 281            break;
 282          }
 283        }
 284        // tau - muon
 285        // If tau within dR < 0.2 of a muon: remove tau
 286        if (away)  {
 287          for (const Particle& mu : muon_candidates) {
 288            if (deltaR(tau.momentum(), mu.momentum()) < 0.2 ) {
 289              away = false;
 290              break;
 291            }
 292          }
 293        }
 294        // If isolated keep it
 295        if (away) recon_tau.push_back( tau );
 296      }
 297
 298      // muon - jet
 299      Particles recon_mu, trigger_mu;
 300      // If muon within dR < 0.4 of a jet: remove muon
 301      for (const Particle& mu : muon_candidates ) {
 302        bool away = true;
 303        for (const Jet& jet : recon_jets) {
 304          if (deltaR(mu.momentum(), jet.momentum()) < 0.4 ) {
 305            away = false;
 306            break;
 307          }
 308        }
 309        if (away)  {
 310          recon_mu.push_back( mu );
 311          recon_leptons.push_back( mu );
 312          if (mu.abseta() < 2.4) trigger_mu.push_back( mu );
 313        }
 314      }
 315
 316      // End overlap removal
 317      // ---------------------
 318
 319      // Jet cleaning
 320      if (rand01() <= 0.42) {
 321        for (const Jet& jet : recon_jets ) {
 322          const double eta = jet.rapidity();
 323          const double phi = jet.azimuthalAngle(MINUSPI_PLUSPI);
 324          if(jet.pT() > 25*GeV && inRange(eta,-0.1,1.5) && inRange(phi,-0.9,-0.5)) vetoEvent;
 325        }
 326      }
 327
 328      // Event selection
 329      // Require at least 3 charged tracks in event
 330      if (charged_tracks.size() < 3) vetoEvent;
 331
 332      // And at least one e/mu passing trigger
 333      if( !( !recon_e.empty() && recon_e[0].pT()>26.*GeV)  &&
 334          !( !trigger_mu.empty() && trigger_mu[0].pT()>26.*GeV) ) {
 335        MSG_DEBUG("Hardest lepton fails trigger");
 336        vetoEvent;
 337      }
 338
 339      // And only accept events with at least 2 electrons and muons and at least 3 leptons in total
 340      if (recon_mu.size() + recon_e.size() + recon_tau.size() < 3 || recon_leptons.size() < 2) vetoEvent;
 341
 342      // Sort leptons by decreasing pT
 343      isortByPt(recon_leptons);
 344      isortByPt(recon_tau);
 345
 346      // Calculate HTlep, fill lepton pT histograms & store chosen combination of 3 leptons
 347      double HTlep = 0.;
 348      Particles chosen_leptons;
 349      if (recon_leptons.size() > 2) {
 350        _h_pt_1_3l->fill(recon_leptons[0].pT()/GeV);
 351        _h_pt_2_3l->fill(recon_leptons[1].pT()/GeV);
 352        _h_pt_3_3l->fill(recon_leptons[2].pT()/GeV);
 353        HTlep = (recon_leptons[0].pT() + recon_leptons[1].pT() + recon_leptons[2].pT())/GeV;
 354        chosen_leptons.push_back( recon_leptons[0] );
 355        chosen_leptons.push_back( recon_leptons[1] );
 356        chosen_leptons.push_back( recon_leptons[2] );
 357      }
 358      else {
 359        _h_pt_1_2ltau->fill(recon_leptons[0].pT()/GeV);
 360        _h_pt_2_2ltau->fill(recon_leptons[1].pT()/GeV);
 361        _h_pt_3_2ltau->fill(recon_tau[0].pT()/GeV);
 362        HTlep = recon_leptons[0].pT()/GeV + recon_leptons[1].pT()/GeV + recon_tau[0].pT()/GeV;
 363        chosen_leptons.push_back( recon_leptons[0] );
 364        chosen_leptons.push_back( recon_leptons[1] );
 365        chosen_leptons.push_back( recon_tau[0] );
 366      }
 367
 368      // Calculate mT and mTW variable
 369      Particles mT_leptons;
 370      Particles mTW_leptons;
 371      for (size_t i1 = 0; i1 < 3; i1 ++)  {
 372        for (size_t i2 = i1+1; i2 < 3; i2 ++)  {
 373          double OSSF_inv_mass = isOSSF_mass(chosen_leptons[i1],chosen_leptons[i2]);
 374          if (OSSF_inv_mass != 0.)  {
 375            for (size_t i3 = 0; i3 < 3 ; i3 ++)  {
 376              if (i3 != i2 && i3 != i1)  {
 377                mT_leptons.push_back(chosen_leptons[i3]);
 378                if ( fabs(91.0 - OSSF_inv_mass) < 20. )
 379                  mTW_leptons.push_back(chosen_leptons[i3]);
 380              }
 381            }
 382          }
 383          else  {
 384            mT_leptons.push_back(chosen_leptons[0]);
 385            mTW_leptons.push_back(chosen_leptons[0]);
 386          }
 387        }
 388      }
 389
 390      isortByPt(mT_leptons);
 391      isortByPt(mTW_leptons);
 392
 393      double mT = sqrt(2*pTmiss.pT()/GeV*mT_leptons[0].pT()/GeV*(1-cos(pTmiss.phi()-mT_leptons[0].phi())));
 394      double mTW = sqrt(2*pTmiss.pT()/GeV*mTW_leptons[0].pT()/GeV*(1-cos(pTmiss.phi()-mTW_leptons[0].phi())));
 395
 396      // Calculate Min pT variable
 397      double min_pT = chosen_leptons[2].pT()/GeV;
 398
 399      // Number of prompt e/mu and had taus
 400      _h_e_n->fill(recon_e.size());
 401      _h_mu_n->fill(recon_mu.size());
 402      _h_tau_n->fill(recon_tau.size());
 403
 404      // Calculate HTjets variable
 405      double HTjets = 0.;
 406      for (const Jet& jet : recon_jets)
 407        HTjets += jet.pT()/GeV;
 408
 409      // Calculate meff variable
 410      double meff = eTmiss + HTjets;
 411      Particles all_leptons;
 412      for (const Particle& e : recon_e ) {
 413        meff += e.pT()/GeV;
 414        all_leptons.push_back( e );
 415      }
 416      for (const Particle& mu : recon_mu) {
 417        meff += mu.pT()/GeV;
 418        all_leptons.push_back( mu );
 419      }
 420      for (const Particle& tau : recon_tau) {
 421        meff += tau.pT()/GeV;
 422        all_leptons.push_back( tau );
 423      }
 424
 425      // Fill histograms of kinematic variables
 426      _h_HTlep_all->fill(HTlep);
 427      _h_HTjets_all->fill(HTjets);
 428      _h_MET_all->fill(eTmiss);
 429      _h_Meff_all->fill(meff);
 430      _h_min_pT_all->fill(min_pT);
 431      _h_mT_all->fill(mT);
 432
 433      // Determine signal region (3l / 2ltau , onZ / offZ OSSF / offZ no-OSSF)
 434      // 3l vs. 2ltau
 435      string basic_signal_region;
 436      if (recon_mu.size() + recon_e.size() > 2)
 437        basic_signal_region += "3l_";
 438      else if ( (recon_mu.size() + recon_e.size() == 2) && (recon_tau.size() > 0))
 439        basic_signal_region += "2ltau_";
 440
 441      // Is there an OSSF pair or a three lepton combination with an invariant mass close to the Z mass
 442      int onZ = isonZ(chosen_leptons);
 443      if (onZ == 1) basic_signal_region += "onZ";
 444      else if (onZ == 0)  {
 445        bool OSSF = isOSSF(chosen_leptons);
 446        if (OSSF) basic_signal_region += "offZ_OSSF";
 447        else basic_signal_region += "offZ_noOSSF";
 448        }
 449
 450      // Check in which signal regions this event falls and adjust event counters
 451      // INFO: The b-jet signal regions of the paper are not included in this Rivet implementation
 452      fillEventCountsPerSR(basic_signal_region,onZ,HTlep,eTmiss,HTjets,meff,min_pT,mTW);
 453    }
 454
 455
 456    /// Normalise histograms etc., after the run
 457    void finalize() {
 458
 459      // Normalize to an integrated luminosity of 1 fb-1
 460      double norm = crossSection()/femtobarn/sumOfWeights();
 461
 462      string best_signal_region = "";
 463      double ratio_best_SR = 0.;
 464
 465      // Loop over all signal regions and find signal region with best sensitivity (ratio signal events/visible cross-section)
 466      for (size_t i = 0; i < _signal_regions.size(); i++) {
 467        double signal_events = _eventCountsPerSR[_signal_regions[i]]->val() * norm;
 468        // Use expected upper limits to find best signal region:
 469        double UL95 = getUpperLimit(_signal_regions[i],false);
 470        double ratio = signal_events / UL95;
 471        if (ratio > ratio_best_SR)  {
 472          best_signal_region = _signal_regions.at(i);
 473          ratio_best_SR = ratio;
 474        }
 475      }
 476
 477      double signal_events_best_SR = _eventCountsPerSR[best_signal_region]->val() * norm;
 478      double exp_UL_best_SR = getUpperLimit(best_signal_region, false);
 479      double obs_UL_best_SR = getUpperLimit(best_signal_region, true);
 480
 481
 482      // Print out result
 483      cout << "----------------------------------------------------------------------------------------" << '\n';
 484      cout << "Number of total events: " << sumOfWeights() << '\n';
 485      cout << "Best signal region: " << best_signal_region << '\n';
 486      cout << "Normalized number of signal events in this best signal region (per fb-1): " << signal_events_best_SR << '\n';
 487      cout << "Efficiency*Acceptance: " << _eventCountsPerSR[best_signal_region]->val()/sumOfWeights() << '\n';
 488      cout << "Cross-section [fb]: " << crossSection()/femtobarn << '\n';
 489      cout << "Expected visible cross-section (per fb-1): " << exp_UL_best_SR << '\n';
 490      cout << "Ratio (signal events / expected visible cross-section): " << ratio_best_SR << '\n';
 491      cout << "Observed visible cross-section (per fb-1): " << obs_UL_best_SR << '\n';
 492      cout << "Ratio (signal events / observed visible cross-section): " <<  signal_events_best_SR/obs_UL_best_SR << '\n';
 493      cout << "----------------------------------------------------------------------------------------" << '\n';
 494
 495      cout << "Using the EXPECTED limits (visible cross-section) of the analysis: " << '\n';
 496      if (signal_events_best_SR > exp_UL_best_SR)  {
 497        cout << "Since the number of signal events > the visible cross-section, this model/grid point is EXCLUDED with 95% C.L." << '\n';
 498        _h_excluded->fill(1);
 499      }
 500      else  {
 501        cout << "Since the number of signal events < the visible cross-section, this model/grid point is NOT EXCLUDED." << '\n';
 502        _h_excluded->fill(0);
 503      }
 504      cout << "----------------------------------------------------------------------------------------" << '\n';
 505
 506      cout << "Using the OBSERVED limits (visible cross-section) of the analysis: " << '\n';
 507      if (signal_events_best_SR > obs_UL_best_SR)  {
 508        cout << "Since the number of signal events > the visible cross-section, this model/grid point is EXCLUDED with 95% C.L." << '\n';
 509        _h_excluded->fill(1);
 510      }
 511      else  {
 512        cout << "Since the number of signal events < the visible cross-section, this model/grid point is NOT EXCLUDED." << '\n';
 513        _h_excluded->fill(0);
 514      }
 515      cout << "----------------------------------------------------------------------------------------" << '\n';
 516      cout << "INFO: The b-jet signal regions of the paper are not included in this Rivet implementation." << '\n';
 517      cout << "----------------------------------------------------------------------------------------" << '\n';
 518
 519
 520      /// Normalize to cross section
 521
 522      if (norm != 0)  {
 523        scale(_h_HTlep_all, norm);
 524        scale(_h_HTjets_all, norm);
 525        scale(_h_MET_all, norm);
 526        scale(_h_Meff_all, norm);
 527        scale(_h_min_pT_all, norm);
 528        scale(_h_mT_all, norm);
 529
 530        scale(_h_pt_1_3l, norm);
 531        scale(_h_pt_2_3l, norm);
 532        scale(_h_pt_3_3l, norm);
 533        scale(_h_pt_1_2ltau, norm);
 534        scale(_h_pt_2_2ltau, norm);
 535        scale(_h_pt_3_2ltau, norm);
 536
 537        scale(_h_e_n, norm);
 538        scale(_h_mu_n, norm);
 539        scale(_h_tau_n, norm);
 540
 541        scale(_h_excluded, norm);
 542      }
 543
 544    }
 545
 546
 547    /// Helper functions
 548    /// @{
 549    /// Function giving a list of all signal regions
 550    vector<string> getSignalRegions()  {
 551
 552      // List of basic signal regions
 553      vector<string> basic_signal_regions;
 554      basic_signal_regions.push_back("3l_offZ_OSSF");
 555      basic_signal_regions.push_back("3l_offZ_noOSSF");
 556      basic_signal_regions.push_back("3l_onZ");
 557      basic_signal_regions.push_back("2ltau_offZ_OSSF");
 558      basic_signal_regions.push_back("2ltau_offZ_noOSSF");
 559      basic_signal_regions.push_back("2ltau_onZ");
 560
 561      // List of kinematic variables
 562      vector<string> kinematic_variables;
 563      kinematic_variables.push_back("HTlep");
 564      kinematic_variables.push_back("METStrong");
 565      kinematic_variables.push_back("METWeak");
 566      kinematic_variables.push_back("Meff");
 567      kinematic_variables.push_back("MeffStrong");
 568      kinematic_variables.push_back("MeffMt");
 569      kinematic_variables.push_back("MinPt");
 570
 571      vector<string> signal_regions;
 572      // Loop over all kinematic variables and basic signal regions
 573      for (size_t i0 = 0; i0 < kinematic_variables.size(); i0++)  {
 574        for (size_t i1 = 0; i1 < basic_signal_regions.size(); i1++)  {
 575          // Is signal region onZ?
 576          int onZ = (basic_signal_regions[i1].find("onZ") != string::npos) ? 1 : 0;
 577          // Get cut values for this kinematic variable
 578          vector<int> cut_values = getCutsPerSignalRegion(kinematic_variables[i0], onZ);
 579          // Loop over all cut values
 580          for (size_t i2 = 0; i2 < cut_values.size(); i2++)  {
 581            // Push signal region into vector
 582            signal_regions.push_back( kinematic_variables[i0] + "_" + basic_signal_regions[i1] + "_cut_" + toString(cut_values[i2]) );
 583          }
 584        }
 585      }
 586      return signal_regions;
 587    }
 588
 589
 590
 591    /// Function giving all cut values per kinematic variable
 592    vector<int> getCutsPerSignalRegion(const string& signal_region, int onZ = 0)  {
 593      vector<int> cutValues;
 594
 595      // Cut values for HTlep
 596      if (signal_region.compare("HTlep") == 0)  {
 597        cutValues.push_back(0);
 598        cutValues.push_back(200);
 599        cutValues.push_back(500);
 600        cutValues.push_back(800);
 601        }
 602      // Cut values for MinPt
 603      else if (signal_region.compare("MinPt") == 0)  {
 604        cutValues.push_back(0);
 605        cutValues.push_back(50);
 606        cutValues.push_back(100);
 607        cutValues.push_back(150);
 608        }
 609      // Cut values for METStrong (HTjets > 150 GeV) and METWeak (HTjets < 150 GeV)
 610      else if (signal_region.compare("METStrong") == 0 || signal_region.compare("METWeak") == 0)  {
 611        cutValues.push_back(0);
 612        cutValues.push_back(100);
 613        cutValues.push_back(200);
 614        cutValues.push_back(300);
 615        }
 616      // Cut values for Meff
 617      if (signal_region.compare("Meff") == 0)  {
 618        cutValues.push_back(0);
 619        cutValues.push_back(600);
 620        cutValues.push_back(1000);
 621        cutValues.push_back(1500);
 622        }
 623      // Cut values for MeffStrong (MET > 100 GeV)
 624      if ((signal_region.compare("MeffStrong") == 0 || signal_region.compare("MeffMt") == 0) && onZ ==1)  {
 625        cutValues.push_back(0);
 626        cutValues.push_back(600);
 627        cutValues.push_back(1200);
 628        }
 629
 630      return cutValues;
 631    }
 632
 633    /// function fills map _eventCountsPerSR by looping over all signal regions
 634    /// and looking if the event falls into this signal region
 635    void fillEventCountsPerSR(const string& basic_signal_region, int onZ,
 636                              double HTlep, double eTmiss, double HTjets,
 637                              double meff, double min_pT, double mTW)  {
 638
 639      // Get cut values for HTlep, loop over them and add event if cut is passed
 640      vector<int> cut_values = getCutsPerSignalRegion("HTlep", onZ);
 641      for (size_t i = 0; i < cut_values.size(); i++)  {
 642        if (HTlep > cut_values[i])
 643          _eventCountsPerSR[("HTlep_" + basic_signal_region + "_cut_" + toString(cut_values[i]))]->fill();
 644      }
 645
 646      // Get cut values for MinPt, loop over them and add event if cut is passed
 647      cut_values = getCutsPerSignalRegion("MinPt", onZ);
 648      for (size_t i = 0; i < cut_values.size(); i++)  {
 649        if (min_pT > cut_values[i])
 650          _eventCountsPerSR[("MinPt_" + basic_signal_region + "_cut_" + toString(cut_values[i]))]->fill();
 651      }
 652
 653      // Get cut values for METStrong, loop over them and add event if cut is passed
 654      cut_values = getCutsPerSignalRegion("METStrong", onZ);
 655      for (size_t i = 0; i < cut_values.size(); i++)  {
 656        if (eTmiss > cut_values[i] && HTjets > 150.)
 657          _eventCountsPerSR[("METStrong_" + basic_signal_region + "_cut_" + toString(cut_values[i]))]->fill();
 658      }
 659
 660      // Get cut values for METWeak, loop over them and add event if cut is passed
 661      cut_values = getCutsPerSignalRegion("METWeak", onZ);
 662      for (size_t i = 0; i < cut_values.size(); i++)  {
 663        if (eTmiss > cut_values[i] && HTjets <= 150.)
 664          _eventCountsPerSR[("METWeak_" + basic_signal_region + "_cut_" + toString(cut_values[i]))]->fill();
 665      }
 666
 667      // Get cut values for Meff, loop over them and add event if cut is passed
 668      cut_values = getCutsPerSignalRegion("Meff", onZ);
 669      for (size_t i = 0; i < cut_values.size(); i++)  {
 670        if (meff > cut_values[i])
 671          _eventCountsPerSR[("Meff_" + basic_signal_region + "_cut_" + toString(cut_values[i]))]->fill();
 672      }
 673
 674      // Get cut values for MeffStrong, loop over them and add event if cut is passed
 675      cut_values = getCutsPerSignalRegion("MeffStrong", onZ);
 676      for (size_t i = 0; i < cut_values.size(); i++)  {
 677        if (meff > cut_values[i] && eTmiss > 100.)
 678          _eventCountsPerSR[("MeffStrong_" + basic_signal_region + "_cut_" + toString(cut_values[i]))]->fill();
 679      }
 680
 681      // Get cut values for MeffMt, loop over them and add event if cut is passed
 682      cut_values = getCutsPerSignalRegion("MeffMt", onZ);
 683      for (size_t i = 0; i < cut_values.size(); i++)  {
 684        if (meff > cut_values[i] && mTW > 100. && onZ == 1)
 685          _eventCountsPerSR[("MeffMt_" + basic_signal_region + "_cut_" + toString(cut_values[i]))]->fill();
 686      }
 687
 688    }
 689
 690    /// Function returning 4-momentum of daughter-particle if it is a tau neutrino
 691    FourMomentum get_tau_neutrino_momentum(const Particle& p)  {
 692      assert(p.abspid() == PID::TAU);
 693      ConstGenVertexPtr dv = p.genParticle()->end_vertex();
 694      assert(dv != nullptr);
 695      // Loop over all daughter particles
 696      for(ConstGenParticlePtr pp: HepMCUtils::particles(dv, Relatives::CHILDREN)){
 697        if (abs(pp->pdg_id()) == PID::NU_TAU) return FourMomentum(pp->momentum());
 698      }
 699      return FourMomentum();
 700    }
 701
 702    /// Function calculating the prong number of taus
 703    void get_prong_number(ConstGenParticlePtr p, unsigned int& nprong, bool& lep_decaying_tau)  {
 704      assert(p != nullptr);
 705      ConstGenVertexPtr dv = p->end_vertex();
 706      assert(dv != nullptr);
 707      for(ConstGenParticlePtr pp: HepMCUtils::particles(dv, Relatives::CHILDREN)){
 708        // If they have status 1 and are charged they will produce a track and the prong number is +1
 709        if (pp->status() == 1 )  {
 710          const int id = pp->pdg_id();
 711          if (Rivet::PID::charge(id) != 0 ) ++nprong;
 712          // Check if tau decays leptonically
 713          if (( abs(id) == PID::ELECTRON || abs(id) == PID::MUON || abs(id) == PID::TAU ) && abs(p->pdg_id()) == PID::TAU) lep_decaying_tau = true;
 714        }
 715        // If the status of the daughter particle is 2 it is unstable and the further decays are checked
 716        else if (pp->status() == 2 )  {
 717          get_prong_number(pp,nprong,lep_decaying_tau);
 718        }
 719      }
 720    }
 721
 722    /// Function giving fiducial lepton efficiency
 723    double apply_reco_eff(int flavor, const Particle& p) {
 724      double pt = p.pT()/GeV;
 725      double eta = p.eta();
 726
 727      double eff = 0.;
 728
 729      if (flavor == 11) { // weight prompt electron -- now including data/MC ID SF in eff.
 730        double avgrate = 0.685;
 731        const static double wz_ele[] =  {0.0256,0.522,0.607,0.654,0.708,0.737,0.761,0.784,0.815,0.835,0.851,0.841,0.898};
 732        // double ewz_ele[] = {0.000257,0.00492,0.00524,0.00519,0.00396,0.00449,0.00538,0.00513,0.00773,0.00753,0.0209,0.0964,0.259};
 733        int ibin = 0;
 734        if(pt > 10  && pt < 15) ibin = 0;
 735        if(pt > 15  && pt < 20) ibin = 1;
 736        if(pt > 20  && pt < 25) ibin = 2;
 737        if(pt > 25  && pt < 30) ibin = 3;
 738        if(pt > 30  && pt < 40) ibin = 4;
 739        if(pt > 40  && pt < 50) ibin = 5;
 740        if(pt > 50  && pt < 60) ibin = 6;
 741        if(pt > 60  && pt < 80) ibin = 7;
 742        if(pt > 80  && pt < 100) ibin = 8;
 743        if(pt > 100 && pt < 200) ibin = 9;
 744        if(pt > 200 && pt < 400) ibin = 10;
 745        if(pt > 400 && pt < 600) ibin = 11;
 746        if(pt > 600) ibin = 12;
 747        double eff_pt = 0.;
 748        eff_pt = wz_ele[ibin];
 749
 750        eta = fabs(eta);
 751
 752        const static double wz_ele_eta[] =  {0.65,0.714,0.722,0.689,0.635,0.615};
 753        // double ewz_ele_eta[] = {0.00642,0.00355,0.00335,0.004,0.00368,0.00422};
 754        ibin = 0;
 755        if(eta > 0 && eta < 0.1) ibin = 0;
 756        if(eta > 0.1 && eta < 0.5) ibin = 1;
 757        if(eta > 0.5 && eta < 1.0) ibin = 2;
 758        if(eta > 1.0 && eta < 1.5) ibin = 3;
 759        if(eta > 1.5 && eta < 2.0) ibin = 4;
 760        if(eta > 2.0 && eta < 2.5) ibin = 5;
 761        double eff_eta = 0.;
 762        eff_eta = wz_ele_eta[ibin];
 763
 764        eff = (eff_pt * eff_eta) / avgrate;
 765      }
 766
 767      if (flavor == 12) { // weight electron from tau
 768        double avgrate = 0.476;
 769        const static double wz_ele[] =  {0.00855,0.409,0.442,0.55,0.632,0.616,0.615,0.642,0.72,0.617};
 770        // double ewz_ele[] = {0.000573,0.0291,0.0366,0.0352,0.0363,0.0474,0.0628,0.0709,0.125,0.109};
 771        int ibin = 0;
 772        if(pt > 10  && pt < 15) ibin = 0;
 773        if(pt > 15  && pt < 20) ibin = 1;
 774        if(pt > 20  && pt < 25) ibin = 2;
 775        if(pt > 25  && pt < 30) ibin = 3;
 776        if(pt > 30  && pt < 40) ibin = 4;
 777        if(pt > 40  && pt < 50) ibin = 5;
 778        if(pt > 50  && pt < 60) ibin = 6;
 779        if(pt > 60  && pt < 80) ibin = 7;
 780        if(pt > 80  && pt < 100) ibin = 8;
 781        if(pt > 100)           ibin = 9;
 782        double eff_pt = 0.;
 783        eff_pt = wz_ele[ibin];
 784
 785        eta = fabs(eta);
 786
 787        const static double wz_ele_eta[] =  {0.546,0.5,0.513,0.421,0.47,0.433};
 788        //double ewz_ele_eta[] = {0.0566,0.0257,0.0263,0.0263,0.0303,0.0321};
 789        ibin = 0;
 790        if(eta > 0 && eta < 0.1) ibin = 0;
 791        if(eta > 0.1 && eta < 0.5) ibin = 1;
 792        if(eta > 0.5 && eta < 1.0) ibin = 2;
 793        if(eta > 1.0 && eta < 1.5) ibin = 3;
 794        if(eta > 1.5 && eta < 2.0) ibin = 4;
 795        if(eta > 2.0 && eta < 2.5) ibin = 5;
 796        double eff_eta = 0.;
 797        eff_eta = wz_ele_eta[ibin];
 798
 799        eff = (eff_pt * eff_eta) / avgrate;
 800      }
 801
 802      if (flavor == 13) { // weight prompt muon
 803        int ibin = 0;
 804        if(pt > 10  && pt < 15) ibin = 0;
 805        if(pt > 15  && pt < 20) ibin = 1;
 806        if(pt > 20  && pt < 25) ibin = 2;
 807        if(pt > 25  && pt < 30) ibin = 3;
 808        if(pt > 30  && pt < 40) ibin = 4;
 809        if(pt > 40  && pt < 50) ibin = 5;
 810        if(pt > 50  && pt < 60) ibin = 6;
 811        if(pt > 60  && pt < 80) ibin = 7;
 812        if(pt > 80  && pt < 100) ibin = 8;
 813        if(pt > 100 && pt < 200) ibin = 9;
 814        if(pt > 200 && pt < 400) ibin = 10;
 815        if(pt > 400) ibin = 11;
 816        if(fabs(eta) < 0.1) {
 817          const static double wz_mu[] =  {0.00705,0.402,0.478,0.49,0.492,0.499,0.527,0.512,0.53,0.528,0.465,0.465};
 818          //double ewz_mu[] = {0.000298,0.0154,0.017,0.0158,0.0114,0.0123,0.0155,0.0133,0.0196,0.0182,0.0414,0.0414};
 819          double eff_pt = 0.;
 820          eff_pt = wz_mu[ibin];
 821          eff = eff_pt;
 822        }
 823        if(fabs(eta) > 0.1) {
 824          const static double wz_mu[] =  {0.0224,0.839,0.887,0.91,0.919,0.923,0.925,0.925,0.922,0.918,0.884,0.834};
 825          //double ewz_mu[] = {0.000213,0.00753,0.0074,0.007,0.00496,0.00534,0.00632,0.00583,0.00849,0.00804,0.0224,0.0963};
 826          double eff_pt = 0.;
 827          eff_pt = wz_mu[ibin];
 828          eff = eff_pt;
 829        }
 830      }
 831
 832      if (flavor == 14) { // weight muon from tau
 833        int ibin = 0;
 834        if(pt > 10  && pt < 15) ibin = 0;
 835        if(pt > 15  && pt < 20) ibin = 1;
 836        if(pt > 20  && pt < 25) ibin = 2;
 837        if(pt > 25  && pt < 30) ibin = 3;
 838        if(pt > 30  && pt < 40) ibin = 4;
 839        if(pt > 40  && pt < 50) ibin = 5;
 840        if(pt > 50  && pt < 60) ibin = 6;
 841        if(pt > 60  && pt < 80) ibin = 7;
 842        if(pt > 80  && pt < 100) ibin = 8;
 843        if(pt > 100) ibin = 9;
 844
 845        if(fabs(eta) < 0.1) {
 846          const static double wz_mu[] =  {0.0,0.664,0.124,0.133,0.527,0.283,0.495,0.25,0.5,0.331};
 847          //double ewz_mu[] = {0.0,0.192,0.0437,0.0343,0.128,0.107,0.202,0.125,0.25,0.191};
 848          double eff_pt = 0.;
 849          eff_pt = wz_mu[ibin];
 850          eff = eff_pt;
 851        }
 852        if(fabs(eta) > 0.1) {
 853          const static double wz_mu[] =  {0.0,0.617,0.655,0.676,0.705,0.738,0.712,0.783,0.646,0.745};
 854          //double ewz_mu[] = {0.0,0.043,0.0564,0.0448,0.0405,0.0576,0.065,0.0825,0.102,0.132};
 855          double eff_pt = 0.;
 856          eff_pt = wz_mu[ibin];
 857          eff = eff_pt;
 858        }
 859      }
 860
 861      if (flavor == 15) { // weight hadronic tau 1p
 862        double avgrate = 0.16;
 863        const static double wz_tau1p[] =  {0.0,0.0311,0.148,0.229,0.217,0.292,0.245,0.307,0.227,0.277};
 864        //double ewz_tau1p[] = {0.0,0.00211,0.0117,0.0179,0.0134,0.0248,0.0264,0.0322,0.0331,0.0427};
 865        int ibin = 0;
 866        if(pt > 10  && pt < 15) ibin = 0;
 867        if(pt > 15  && pt < 20) ibin = 1;
 868        if(pt > 20  && pt < 25) ibin = 2;
 869        if(pt > 25  && pt < 30) ibin = 3;
 870        if(pt > 30  && pt < 40) ibin = 4;
 871        if(pt > 40  && pt < 50) ibin = 5;
 872        if(pt > 50  && pt < 60) ibin = 6;
 873        if(pt > 60  && pt < 80) ibin = 7;
 874        if(pt > 80  && pt < 100) ibin = 8;
 875        if(pt > 100) ibin = 9;
 876        double eff_pt = 0.;
 877        eff_pt = wz_tau1p[ibin];
 878
 879        const static double wz_tau1p_eta[] = {0.166,0.15,0.188,0.175,0.142,0.109};
 880        //double ewz_tau1p_eta[] ={0.0166,0.00853,0.0097,0.00985,0.00949,0.00842};
 881        ibin = 0;
 882        if(eta > 0.0 && eta < 0.1) ibin = 0;
 883        if(eta > 0.1 && eta < 0.5) ibin = 1;
 884        if(eta > 0.5 && eta < 1.0) ibin = 2;
 885        if(eta > 1.0 && eta < 1.5) ibin = 3;
 886        if(eta > 1.5 && eta < 2.0) ibin = 4;
 887        if(eta > 2.0 && eta < 2.5) ibin = 5;
 888        double eff_eta = 0.;
 889        eff_eta = wz_tau1p_eta[ibin];
 890
 891        eff = (eff_pt * eff_eta) / avgrate;
 892      }
 893
 894      return eff;
 895    }
 896
 897
 898    /// Function giving observed and expected upper limits (on the visible cross-section)
 899    double getUpperLimit(const string& signal_region, bool observed)  {
 900
 901      map<string,double> upperLimitsObserved;
 902      map<string,double> upperLimitsExpected;
 903
 904      upperLimitsObserved["HTlep_3l_offZ_OSSF_cut_0"] = 2.435;
 905      upperLimitsObserved["HTlep_3l_offZ_OSSF_cut_200"] = 0.704;
 906      upperLimitsObserved["HTlep_3l_offZ_OSSF_cut_500"] = 0.182;
 907      upperLimitsObserved["HTlep_3l_offZ_OSSF_cut_800"] = 0.147;
 908      upperLimitsObserved["HTlep_2ltau_offZ_OSSF_cut_0"] = 13.901;
 909      upperLimitsObserved["HTlep_2ltau_offZ_OSSF_cut_200"] = 1.677;
 910      upperLimitsObserved["HTlep_2ltau_offZ_OSSF_cut_500"] = 0.141;
 911      upperLimitsObserved["HTlep_2ltau_offZ_OSSF_cut_800"] = 0.155;
 912      upperLimitsObserved["HTlep_3l_offZ_noOSSF_cut_0"] = 1.054;
 913      upperLimitsObserved["HTlep_3l_offZ_noOSSF_cut_200"] = 0.341;
 914      upperLimitsObserved["HTlep_3l_offZ_noOSSF_cut_500"] = 0.221;
 915      upperLimitsObserved["HTlep_3l_offZ_noOSSF_cut_800"] = 0.140;
 916      upperLimitsObserved["HTlep_2ltau_offZ_noOSSF_cut_0"] = 4.276;
 917      upperLimitsObserved["HTlep_2ltau_offZ_noOSSF_cut_200"] = 0.413;
 918      upperLimitsObserved["HTlep_2ltau_offZ_noOSSF_cut_500"] = 0.138;
 919      upperLimitsObserved["HTlep_2ltau_offZ_noOSSF_cut_800"] = 0.150;
 920      upperLimitsObserved["HTlep_3l_onZ_cut_0"] = 29.804;
 921      upperLimitsObserved["HTlep_3l_onZ_cut_200"] = 3.579;
 922      upperLimitsObserved["HTlep_3l_onZ_cut_500"] = 0.466;
 923      upperLimitsObserved["HTlep_3l_onZ_cut_800"] = 0.298;
 924      upperLimitsObserved["HTlep_2ltau_onZ_cut_0"] = 205.091;
 925      upperLimitsObserved["HTlep_2ltau_onZ_cut_200"] = 3.141;
 926      upperLimitsObserved["HTlep_2ltau_onZ_cut_500"] = 0.290;
 927      upperLimitsObserved["HTlep_2ltau_onZ_cut_800"] = 0.157;
 928      upperLimitsObserved["METStrong_3l_offZ_OSSF_cut_0"] = 1.111;
 929      upperLimitsObserved["METStrong_3l_offZ_OSSF_cut_100"] = 0.354;
 930      upperLimitsObserved["METStrong_3l_offZ_OSSF_cut_200"] = 0.236;
 931      upperLimitsObserved["METStrong_3l_offZ_OSSF_cut_300"] = 0.150;
 932      upperLimitsObserved["METStrong_2ltau_offZ_OSSF_cut_0"] = 1.881;
 933      upperLimitsObserved["METStrong_2ltau_offZ_OSSF_cut_100"] = 0.406;
 934      upperLimitsObserved["METStrong_2ltau_offZ_OSSF_cut_200"] = 0.194;
 935      upperLimitsObserved["METStrong_2ltau_offZ_OSSF_cut_300"] = 0.134;
 936      upperLimitsObserved["METStrong_3l_offZ_noOSSF_cut_0"] = 0.770;
 937      upperLimitsObserved["METStrong_3l_offZ_noOSSF_cut_100"] = 0.295;
 938      upperLimitsObserved["METStrong_3l_offZ_noOSSF_cut_200"] = 0.149;
 939      upperLimitsObserved["METStrong_3l_offZ_noOSSF_cut_300"] = 0.140;
 940      upperLimitsObserved["METStrong_2ltau_offZ_noOSSF_cut_0"] = 2.003;
 941      upperLimitsObserved["METStrong_2ltau_offZ_noOSSF_cut_100"] = 0.806;
 942      upperLimitsObserved["METStrong_2ltau_offZ_noOSSF_cut_200"] = 0.227;
 943      upperLimitsObserved["METStrong_2ltau_offZ_noOSSF_cut_300"] = 0.138;
 944      upperLimitsObserved["METStrong_3l_onZ_cut_0"] = 6.383;
 945      upperLimitsObserved["METStrong_3l_onZ_cut_100"] = 0.959;
 946      upperLimitsObserved["METStrong_3l_onZ_cut_200"] = 0.549;
 947      upperLimitsObserved["METStrong_3l_onZ_cut_300"] = 0.182;
 948      upperLimitsObserved["METStrong_2ltau_onZ_cut_0"] = 10.658;
 949      upperLimitsObserved["METStrong_2ltau_onZ_cut_100"] = 0.637;
 950      upperLimitsObserved["METStrong_2ltau_onZ_cut_200"] = 0.291;
 951      upperLimitsObserved["METStrong_2ltau_onZ_cut_300"] = 0.227;
 952      upperLimitsObserved["METWeak_3l_offZ_OSSF_cut_0"] = 1.802;
 953      upperLimitsObserved["METWeak_3l_offZ_OSSF_cut_100"] = 0.344;
 954      upperLimitsObserved["METWeak_3l_offZ_OSSF_cut_200"] = 0.189;
 955      upperLimitsObserved["METWeak_3l_offZ_OSSF_cut_300"] = 0.148;
 956      upperLimitsObserved["METWeak_2ltau_offZ_OSSF_cut_0"] = 12.321;
 957      upperLimitsObserved["METWeak_2ltau_offZ_OSSF_cut_100"] = 0.430;
 958      upperLimitsObserved["METWeak_2ltau_offZ_OSSF_cut_200"] = 0.137;
 959      upperLimitsObserved["METWeak_2ltau_offZ_OSSF_cut_300"] = 0.134;
 960      upperLimitsObserved["METWeak_3l_offZ_noOSSF_cut_0"] = 0.562;
 961      upperLimitsObserved["METWeak_3l_offZ_noOSSF_cut_100"] = 0.153;
 962      upperLimitsObserved["METWeak_3l_offZ_noOSSF_cut_200"] = 0.154;
 963      upperLimitsObserved["METWeak_3l_offZ_noOSSF_cut_300"] = 0.141;
 964      upperLimitsObserved["METWeak_2ltau_offZ_noOSSF_cut_0"] = 2.475;
 965      upperLimitsObserved["METWeak_2ltau_offZ_noOSSF_cut_100"] = 0.244;
 966      upperLimitsObserved["METWeak_2ltau_offZ_noOSSF_cut_200"] = 0.141;
 967      upperLimitsObserved["METWeak_2ltau_offZ_noOSSF_cut_300"] = 0.142;
 968      upperLimitsObserved["METWeak_3l_onZ_cut_0"] = 24.769;
 969      upperLimitsObserved["METWeak_3l_onZ_cut_100"] = 0.690;
 970      upperLimitsObserved["METWeak_3l_onZ_cut_200"] = 0.198;
 971      upperLimitsObserved["METWeak_3l_onZ_cut_300"] = 0.138;
 972      upperLimitsObserved["METWeak_2ltau_onZ_cut_0"] = 194.360;
 973      upperLimitsObserved["METWeak_2ltau_onZ_cut_100"] = 0.287;
 974      upperLimitsObserved["METWeak_2ltau_onZ_cut_200"] = 0.144;
 975      upperLimitsObserved["METWeak_2ltau_onZ_cut_300"] = 0.130;
 976      upperLimitsObserved["Meff_3l_offZ_OSSF_cut_0"] = 2.435;
 977      upperLimitsObserved["Meff_3l_offZ_OSSF_cut_600"] = 0.487;
 978      upperLimitsObserved["Meff_3l_offZ_OSSF_cut_1000"] = 0.156;
 979      upperLimitsObserved["Meff_3l_offZ_OSSF_cut_1500"] = 0.140;
 980      upperLimitsObserved["Meff_2ltau_offZ_OSSF_cut_0"] = 13.901;
 981      upperLimitsObserved["Meff_2ltau_offZ_OSSF_cut_600"] = 0.687;
 982      upperLimitsObserved["Meff_2ltau_offZ_OSSF_cut_1000"] = 0.224;
 983      upperLimitsObserved["Meff_2ltau_offZ_OSSF_cut_1500"] = 0.155;
 984      upperLimitsObserved["Meff_3l_offZ_noOSSF_cut_0"] = 1.054;
 985      upperLimitsObserved["Meff_3l_offZ_noOSSF_cut_600"] = 0.249;
 986      upperLimitsObserved["Meff_3l_offZ_noOSSF_cut_1000"] = 0.194;
 987      upperLimitsObserved["Meff_3l_offZ_noOSSF_cut_1500"] = 0.145;
 988      upperLimitsObserved["Meff_2ltau_offZ_noOSSF_cut_0"] = 4.276;
 989      upperLimitsObserved["Meff_2ltau_offZ_noOSSF_cut_600"] = 0.772;
 990      upperLimitsObserved["Meff_2ltau_offZ_noOSSF_cut_1000"] = 0.218;
 991      upperLimitsObserved["Meff_2ltau_offZ_noOSSF_cut_1500"] = 0.204;
 992      upperLimitsObserved["Meff_3l_onZ_cut_0"] = 29.804;
 993      upperLimitsObserved["Meff_3l_onZ_cut_600"] = 2.933;
 994      upperLimitsObserved["Meff_3l_onZ_cut_1000"] = 0.912;
 995      upperLimitsObserved["Meff_3l_onZ_cut_1500"] = 0.225;
 996      upperLimitsObserved["Meff_2ltau_onZ_cut_0"] = 205.091;
 997      upperLimitsObserved["Meff_2ltau_onZ_cut_600"] = 1.486;
 998      upperLimitsObserved["Meff_2ltau_onZ_cut_1000"] = 0.641;
 999      upperLimitsObserved["Meff_2ltau_onZ_cut_1500"] = 0.204;
1000      upperLimitsObserved["MeffStrong_3l_offZ_OSSF_cut_0"] = 0.479;
1001      upperLimitsObserved["MeffStrong_3l_offZ_OSSF_cut_600"] = 0.353;
1002      upperLimitsObserved["MeffStrong_3l_offZ_OSSF_cut_1200"] = 0.187;
1003      upperLimitsObserved["MeffStrong_2ltau_offZ_OSSF_cut_0"] = 0.617;
1004      upperLimitsObserved["MeffStrong_2ltau_offZ_OSSF_cut_600"] = 0.320;
1005      upperLimitsObserved["MeffStrong_2ltau_offZ_OSSF_cut_1200"] = 0.281;
1006      upperLimitsObserved["MeffStrong_3l_offZ_noOSSF_cut_0"] = 0.408;
1007      upperLimitsObserved["MeffStrong_3l_offZ_noOSSF_cut_600"] = 0.240;
1008      upperLimitsObserved["MeffStrong_3l_offZ_noOSSF_cut_1200"] = 0.150;
1009      upperLimitsObserved["MeffStrong_2ltau_offZ_noOSSF_cut_0"] = 0.774;
1010      upperLimitsObserved["MeffStrong_2ltau_offZ_noOSSF_cut_600"] = 0.417;
1011      upperLimitsObserved["MeffStrong_2ltau_offZ_noOSSF_cut_1200"] = 0.266;
1012      upperLimitsObserved["MeffStrong_3l_onZ_cut_0"] = 1.208;
1013      upperLimitsObserved["MeffStrong_3l_onZ_cut_600"] = 0.837;
1014      upperLimitsObserved["MeffStrong_3l_onZ_cut_1200"] = 0.269;
1015      upperLimitsObserved["MeffStrong_2ltau_onZ_cut_0"] = 0.605;
1016      upperLimitsObserved["MeffStrong_2ltau_onZ_cut_600"] = 0.420;
1017      upperLimitsObserved["MeffStrong_2ltau_onZ_cut_1200"] = 0.141;
1018      upperLimitsObserved["MeffMt_3l_onZ_cut_0"] = 1.832;
1019      upperLimitsObserved["MeffMt_3l_onZ_cut_600"] = 0.862;
1020      upperLimitsObserved["MeffMt_3l_onZ_cut_1200"] = 0.222;
1021      upperLimitsObserved["MeffMt_2ltau_onZ_cut_0"] = 1.309;
1022      upperLimitsObserved["MeffMt_2ltau_onZ_cut_600"] = 0.481;
1023      upperLimitsObserved["MeffMt_2ltau_onZ_cut_1200"] = 0.146;
1024      upperLimitsObserved["MinPt_3l_offZ_OSSF_cut_0"] = 2.435;
1025      upperLimitsObserved["MinPt_3l_offZ_OSSF_cut_50"] = 0.500;
1026      upperLimitsObserved["MinPt_3l_offZ_OSSF_cut_100"] = 0.203;
1027      upperLimitsObserved["MinPt_3l_offZ_OSSF_cut_150"] = 0.128;
1028      upperLimitsObserved["MinPt_2ltau_offZ_OSSF_cut_0"] = 13.901;
1029      upperLimitsObserved["MinPt_2ltau_offZ_OSSF_cut_50"] = 0.859;
1030      upperLimitsObserved["MinPt_2ltau_offZ_OSSF_cut_100"] = 0.158;
1031      upperLimitsObserved["MinPt_2ltau_offZ_OSSF_cut_150"] = 0.155;
1032      upperLimitsObserved["MinPt_3l_offZ_noOSSF_cut_0"] = 1.054;
1033      upperLimitsObserved["MinPt_3l_offZ_noOSSF_cut_50"] = 0.295;
1034      upperLimitsObserved["MinPt_3l_offZ_noOSSF_cut_100"] = 0.148;
1035      upperLimitsObserved["MinPt_3l_offZ_noOSSF_cut_150"] = 0.137;
1036      upperLimitsObserved["MinPt_2ltau_offZ_noOSSF_cut_0"] = 4.276;
1037      upperLimitsObserved["MinPt_2ltau_offZ_noOSSF_cut_50"] = 0.314;
1038      upperLimitsObserved["MinPt_2ltau_offZ_noOSSF_cut_100"] = 0.134;
1039      upperLimitsObserved["MinPt_2ltau_offZ_noOSSF_cut_150"] = 0.140;
1040      upperLimitsObserved["MinPt_3l_onZ_cut_0"] = 29.804;
1041      upperLimitsObserved["MinPt_3l_onZ_cut_50"] = 1.767;
1042      upperLimitsObserved["MinPt_3l_onZ_cut_100"] = 0.690;
1043      upperLimitsObserved["MinPt_3l_onZ_cut_150"] = 0.301;
1044      upperLimitsObserved["MinPt_2ltau_onZ_cut_0"] = 205.091;
1045      upperLimitsObserved["MinPt_2ltau_onZ_cut_50"] = 1.050;
1046      upperLimitsObserved["MinPt_2ltau_onZ_cut_100"] = 0.155;
1047      upperLimitsObserved["MinPt_2ltau_onZ_cut_150"] = 0.146;
1048      upperLimitsObserved["nbtag_3l_offZ_OSSF_cut_0"] = 2.435;
1049      upperLimitsObserved["nbtag_3l_offZ_OSSF_cut_1"] = 0.865;
1050      upperLimitsObserved["nbtag_3l_offZ_OSSF_cut_2"] = 0.474;
1051      upperLimitsObserved["nbtag_2ltau_offZ_OSSF_cut_0"] = 13.901;
1052      upperLimitsObserved["nbtag_2ltau_offZ_OSSF_cut_1"] = 1.566;
1053      upperLimitsObserved["nbtag_2ltau_offZ_OSSF_cut_2"] = 0.426;
1054      upperLimitsObserved["nbtag_3l_offZ_noOSSF_cut_0"] = 1.054;
1055      upperLimitsObserved["nbtag_3l_offZ_noOSSF_cut_1"] = 0.643;
1056      upperLimitsObserved["nbtag_3l_offZ_noOSSF_cut_2"] = 0.321;
1057      upperLimitsObserved["nbtag_2ltau_offZ_noOSSF_cut_0"] = 4.276;
1058      upperLimitsObserved["nbtag_2ltau_offZ_noOSSF_cut_1"] = 2.435;
1059      upperLimitsObserved["nbtag_2ltau_offZ_noOSSF_cut_2"] = 1.073;
1060      upperLimitsObserved["nbtag_3l_onZ_cut_0"] = 29.804;
1061      upperLimitsObserved["nbtag_3l_onZ_cut_1"] = 3.908;
1062      upperLimitsObserved["nbtag_3l_onZ_cut_2"] = 0.704;
1063      upperLimitsObserved["nbtag_2ltau_onZ_cut_0"] = 205.091;
1064      upperLimitsObserved["nbtag_2ltau_onZ_cut_1"] = 9.377;
1065      upperLimitsObserved["nbtag_2ltau_onZ_cut_2"] = 0.657;
1066      upperLimitsExpected["HTlep_3l_offZ_OSSF_cut_0"] = 2.893;
1067      upperLimitsExpected["HTlep_3l_offZ_OSSF_cut_200"] = 1.175;
1068      upperLimitsExpected["HTlep_3l_offZ_OSSF_cut_500"] = 0.265;
1069      upperLimitsExpected["HTlep_3l_offZ_OSSF_cut_800"] = 0.155;
1070      upperLimitsExpected["HTlep_2ltau_offZ_OSSF_cut_0"] = 14.293;
1071      upperLimitsExpected["HTlep_2ltau_offZ_OSSF_cut_200"] = 1.803;
1072      upperLimitsExpected["HTlep_2ltau_offZ_OSSF_cut_500"] = 0.159;
1073      upperLimitsExpected["HTlep_2ltau_offZ_OSSF_cut_800"] = 0.155;
1074      upperLimitsExpected["HTlep_3l_offZ_noOSSF_cut_0"] = 0.836;
1075      upperLimitsExpected["HTlep_3l_offZ_noOSSF_cut_200"] = 0.340;
1076      upperLimitsExpected["HTlep_3l_offZ_noOSSF_cut_500"] = 0.218;
1077      upperLimitsExpected["HTlep_3l_offZ_noOSSF_cut_800"] = 0.140;
1078      upperLimitsExpected["HTlep_2ltau_offZ_noOSSF_cut_0"] = 4.132;
1079      upperLimitsExpected["HTlep_2ltau_offZ_noOSSF_cut_200"] = 0.599;
1080      upperLimitsExpected["HTlep_2ltau_offZ_noOSSF_cut_500"] = 0.146;
1081      upperLimitsExpected["HTlep_2ltau_offZ_noOSSF_cut_800"] = 0.148;
1082      upperLimitsExpected["HTlep_3l_onZ_cut_0"] = 32.181;
1083      upperLimitsExpected["HTlep_3l_onZ_cut_200"] = 4.879;
1084      upperLimitsExpected["HTlep_3l_onZ_cut_500"] = 0.473;
1085      upperLimitsExpected["HTlep_3l_onZ_cut_800"] = 0.266;
1086      upperLimitsExpected["HTlep_2ltau_onZ_cut_0"] = 217.801;
1087      upperLimitsExpected["HTlep_2ltau_onZ_cut_200"] = 3.676;
1088      upperLimitsExpected["HTlep_2ltau_onZ_cut_500"] = 0.235;
1089      upperLimitsExpected["HTlep_2ltau_onZ_cut_800"] = 0.150;
1090      upperLimitsExpected["METStrong_3l_offZ_OSSF_cut_0"] = 1.196;
1091      upperLimitsExpected["METStrong_3l_offZ_OSSF_cut_100"] = 0.423;
1092      upperLimitsExpected["METStrong_3l_offZ_OSSF_cut_200"] = 0.208;
1093      upperLimitsExpected["METStrong_3l_offZ_OSSF_cut_300"] = 0.158;
1094      upperLimitsExpected["METStrong_2ltau_offZ_OSSF_cut_0"] = 2.158;
1095      upperLimitsExpected["METStrong_2ltau_offZ_OSSF_cut_100"] = 0.461;
1096      upperLimitsExpected["METStrong_2ltau_offZ_OSSF_cut_200"] = 0.186;
1097      upperLimitsExpected["METStrong_2ltau_offZ_OSSF_cut_300"] = 0.138;
1098      upperLimitsExpected["METStrong_3l_offZ_noOSSF_cut_0"] = 0.495;
1099      upperLimitsExpected["METStrong_3l_offZ_noOSSF_cut_100"] = 0.284;
1100      upperLimitsExpected["METStrong_3l_offZ_noOSSF_cut_200"] = 0.150;
1101      upperLimitsExpected["METStrong_3l_offZ_noOSSF_cut_300"] = 0.146;
1102      upperLimitsExpected["METStrong_2ltau_offZ_noOSSF_cut_0"] = 1.967;
1103      upperLimitsExpected["METStrong_2ltau_offZ_noOSSF_cut_100"] = 0.732;
1104      upperLimitsExpected["METStrong_2ltau_offZ_noOSSF_cut_200"] = 0.225;
1105      upperLimitsExpected["METStrong_2ltau_offZ_noOSSF_cut_300"] = 0.147;
1106      upperLimitsExpected["METStrong_3l_onZ_cut_0"] = 7.157;
1107      upperLimitsExpected["METStrong_3l_onZ_cut_100"] = 1.342;
1108      upperLimitsExpected["METStrong_3l_onZ_cut_200"] = 0.508;
1109      upperLimitsExpected["METStrong_3l_onZ_cut_300"] = 0.228;
1110      upperLimitsExpected["METStrong_2ltau_onZ_cut_0"] = 12.441;
1111      upperLimitsExpected["METStrong_2ltau_onZ_cut_100"] = 0.534;
1112      upperLimitsExpected["METStrong_2ltau_onZ_cut_200"] = 0.243;
1113      upperLimitsExpected["METStrong_2ltau_onZ_cut_300"] = 0.218;
1114      upperLimitsExpected["METWeak_3l_offZ_OSSF_cut_0"] = 2.199;
1115      upperLimitsExpected["METWeak_3l_offZ_OSSF_cut_100"] = 0.391;
1116      upperLimitsExpected["METWeak_3l_offZ_OSSF_cut_200"] = 0.177;
1117      upperLimitsExpected["METWeak_3l_offZ_OSSF_cut_300"] = 0.144;
1118      upperLimitsExpected["METWeak_2ltau_offZ_OSSF_cut_0"] = 12.431;
1119      upperLimitsExpected["METWeak_2ltau_offZ_OSSF_cut_100"] = 0.358;
1120      upperLimitsExpected["METWeak_2ltau_offZ_OSSF_cut_200"] = 0.150;
1121      upperLimitsExpected["METWeak_2ltau_offZ_OSSF_cut_300"] = 0.135;
1122      upperLimitsExpected["METWeak_3l_offZ_noOSSF_cut_0"] = 0.577;
1123      upperLimitsExpected["METWeak_3l_offZ_noOSSF_cut_100"] = 0.214;
1124      upperLimitsExpected["METWeak_3l_offZ_noOSSF_cut_200"] = 0.155;
1125      upperLimitsExpected["METWeak_3l_offZ_noOSSF_cut_300"] = 0.140;
1126      upperLimitsExpected["METWeak_2ltau_offZ_noOSSF_cut_0"] = 2.474;
1127      upperLimitsExpected["METWeak_2ltau_offZ_noOSSF_cut_100"] = 0.382;
1128      upperLimitsExpected["METWeak_2ltau_offZ_noOSSF_cut_200"] = 0.144;
1129      upperLimitsExpected["METWeak_2ltau_offZ_noOSSF_cut_300"] = 0.146;
1130      upperLimitsExpected["METWeak_3l_onZ_cut_0"] = 26.305;
1131      upperLimitsExpected["METWeak_3l_onZ_cut_100"] = 1.227;
1132      upperLimitsExpected["METWeak_3l_onZ_cut_200"] = 0.311;
1133      upperLimitsExpected["METWeak_3l_onZ_cut_300"] = 0.188;
1134      upperLimitsExpected["METWeak_2ltau_onZ_cut_0"] = 205.198;
1135      upperLimitsExpected["METWeak_2ltau_onZ_cut_100"] = 0.399;
1136      upperLimitsExpected["METWeak_2ltau_onZ_cut_200"] = 0.166;
1137      upperLimitsExpected["METWeak_2ltau_onZ_cut_300"] = 0.140;
1138      upperLimitsExpected["Meff_3l_offZ_OSSF_cut_0"] = 2.893;
1139      upperLimitsExpected["Meff_3l_offZ_OSSF_cut_600"] = 0.649;
1140      upperLimitsExpected["Meff_3l_offZ_OSSF_cut_1000"] = 0.252;
1141      upperLimitsExpected["Meff_3l_offZ_OSSF_cut_1500"] = 0.150;
1142      upperLimitsExpected["Meff_2ltau_offZ_OSSF_cut_0"] = 14.293;
1143      upperLimitsExpected["Meff_2ltau_offZ_OSSF_cut_600"] = 0.657;
1144      upperLimitsExpected["Meff_2ltau_offZ_OSSF_cut_1000"] = 0.226;
1145      upperLimitsExpected["Meff_2ltau_offZ_OSSF_cut_1500"] = 0.154;
1146      upperLimitsExpected["Meff_3l_offZ_noOSSF_cut_0"] = 0.836;
1147      upperLimitsExpected["Meff_3l_offZ_noOSSF_cut_600"] = 0.265;
1148      upperLimitsExpected["Meff_3l_offZ_noOSSF_cut_1000"] = 0.176;
1149      upperLimitsExpected["Meff_3l_offZ_noOSSF_cut_1500"] = 0.146;
1150      upperLimitsExpected["Meff_2ltau_offZ_noOSSF_cut_0"] = 4.132;
1151      upperLimitsExpected["Meff_2ltau_offZ_noOSSF_cut_600"] = 0.678;
1152      upperLimitsExpected["Meff_2ltau_offZ_noOSSF_cut_1000"] = 0.243;
1153      upperLimitsExpected["Meff_2ltau_offZ_noOSSF_cut_1500"] = 0.184;
1154      upperLimitsExpected["Meff_3l_onZ_cut_0"] = 32.181;
1155      upperLimitsExpected["Meff_3l_onZ_cut_600"] = 3.219;
1156      upperLimitsExpected["Meff_3l_onZ_cut_1000"] = 0.905;
1157      upperLimitsExpected["Meff_3l_onZ_cut_1500"] = 0.261;
1158      upperLimitsExpected["Meff_2ltau_onZ_cut_0"] = 217.801;
1159      upperLimitsExpected["Meff_2ltau_onZ_cut_600"] = 1.680;
1160      upperLimitsExpected["Meff_2ltau_onZ_cut_1000"] = 0.375;
1161      upperLimitsExpected["Meff_2ltau_onZ_cut_1500"] = 0.178;
1162      upperLimitsExpected["MeffStrong_3l_offZ_OSSF_cut_0"] = 0.571;
1163      upperLimitsExpected["MeffStrong_3l_offZ_OSSF_cut_600"] = 0.386;
1164      upperLimitsExpected["MeffStrong_3l_offZ_OSSF_cut_1200"] = 0.177;
1165      upperLimitsExpected["MeffStrong_2ltau_offZ_OSSF_cut_0"] = 0.605;
1166      upperLimitsExpected["MeffStrong_2ltau_offZ_OSSF_cut_600"] = 0.335;
1167      upperLimitsExpected["MeffStrong_2ltau_offZ_OSSF_cut_1200"] = 0.249;
1168      upperLimitsExpected["MeffStrong_3l_offZ_noOSSF_cut_0"] = 0.373;
1169      upperLimitsExpected["MeffStrong_3l_offZ_noOSSF_cut_600"] = 0.223;
1170      upperLimitsExpected["MeffStrong_3l_offZ_noOSSF_cut_1200"] = 0.150;
1171      upperLimitsExpected["MeffStrong_2ltau_offZ_noOSSF_cut_0"] = 0.873;
1172      upperLimitsExpected["MeffStrong_2ltau_offZ_noOSSF_cut_600"] = 0.428;
1173      upperLimitsExpected["MeffStrong_2ltau_offZ_noOSSF_cut_1200"] = 0.210;
1174      upperLimitsExpected["MeffStrong_3l_onZ_cut_0"] = 2.034;
1175      upperLimitsExpected["MeffStrong_3l_onZ_cut_600"] = 1.093;
1176      upperLimitsExpected["MeffStrong_3l_onZ_cut_1200"] = 0.293;
1177      upperLimitsExpected["MeffStrong_2ltau_onZ_cut_0"] = 0.690;
1178      upperLimitsExpected["MeffStrong_2ltau_onZ_cut_600"] = 0.392;
1179      upperLimitsExpected["MeffStrong_2ltau_onZ_cut_1200"] = 0.156;
1180      upperLimitsExpected["MeffMt_3l_onZ_cut_0"] = 2.483;
1181      upperLimitsExpected["MeffMt_3l_onZ_cut_600"] = 0.845;
1182      upperLimitsExpected["MeffMt_3l_onZ_cut_1200"] = 0.255;
1183      upperLimitsExpected["MeffMt_2ltau_onZ_cut_0"] = 1.448;
1184      upperLimitsExpected["MeffMt_2ltau_onZ_cut_600"] = 0.391;
1185      upperLimitsExpected["MeffMt_2ltau_onZ_cut_1200"] = 0.146;
1186      upperLimitsExpected["MinPt_3l_offZ_OSSF_cut_0"] = 2.893;
1187      upperLimitsExpected["MinPt_3l_offZ_OSSF_cut_50"] = 0.703;
1188      upperLimitsExpected["MinPt_3l_offZ_OSSF_cut_100"] = 0.207;
1189      upperLimitsExpected["MinPt_3l_offZ_OSSF_cut_150"] = 0.143;
1190      upperLimitsExpected["MinPt_2ltau_offZ_OSSF_cut_0"] = 14.293;
1191      upperLimitsExpected["MinPt_2ltau_offZ_OSSF_cut_50"] = 0.705;
1192      upperLimitsExpected["MinPt_2ltau_offZ_OSSF_cut_100"] = 0.149;
1193      upperLimitsExpected["MinPt_2ltau_offZ_OSSF_cut_150"] = 0.155;
1194      upperLimitsExpected["MinPt_3l_offZ_noOSSF_cut_0"] = 0.836;
1195      upperLimitsExpected["MinPt_3l_offZ_noOSSF_cut_50"] = 0.249;
1196      upperLimitsExpected["MinPt_3l_offZ_noOSSF_cut_100"] = 0.135;
1197      upperLimitsExpected["MinPt_3l_offZ_noOSSF_cut_150"] = 0.136;
1198      upperLimitsExpected["MinPt_2ltau_offZ_noOSSF_cut_0"] = 4.132;
1199      upperLimitsExpected["MinPt_2ltau_offZ_noOSSF_cut_50"] = 0.339;
1200      upperLimitsExpected["MinPt_2ltau_offZ_noOSSF_cut_100"] = 0.149;
1201      upperLimitsExpected["MinPt_2ltau_offZ_noOSSF_cut_150"] = 0.145;
1202      upperLimitsExpected["MinPt_3l_onZ_cut_0"] = 32.181;
1203      upperLimitsExpected["MinPt_3l_onZ_cut_50"] = 2.260;
1204      upperLimitsExpected["MinPt_3l_onZ_cut_100"] = 0.438;
1205      upperLimitsExpected["MinPt_3l_onZ_cut_150"] = 0.305;
1206      upperLimitsExpected["MinPt_2ltau_onZ_cut_0"] = 217.801;
1207      upperLimitsExpected["MinPt_2ltau_onZ_cut_50"] = 1.335;
1208      upperLimitsExpected["MinPt_2ltau_onZ_cut_100"] = 0.162;
1209      upperLimitsExpected["MinPt_2ltau_onZ_cut_150"] = 0.149;
1210      upperLimitsExpected["nbtag_3l_offZ_OSSF_cut_0"] = 2.893;
1211      upperLimitsExpected["nbtag_3l_offZ_OSSF_cut_1"] = 0.923;
1212      upperLimitsExpected["nbtag_3l_offZ_OSSF_cut_2"] = 0.452;
1213      upperLimitsExpected["nbtag_2ltau_offZ_OSSF_cut_0"] = 14.293;
1214      upperLimitsExpected["nbtag_2ltau_offZ_OSSF_cut_1"] = 1.774;
1215      upperLimitsExpected["nbtag_2ltau_offZ_OSSF_cut_2"] = 0.549;
1216      upperLimitsExpected["nbtag_3l_offZ_noOSSF_cut_0"] = 0.836;
1217      upperLimitsExpected["nbtag_3l_offZ_noOSSF_cut_1"] = 0.594;
1218      upperLimitsExpected["nbtag_3l_offZ_noOSSF_cut_2"] = 0.298;
1219      upperLimitsExpected["nbtag_2ltau_offZ_noOSSF_cut_0"] = 4.132;
1220      upperLimitsExpected["nbtag_2ltau_offZ_noOSSF_cut_1"] = 2.358;
1221      upperLimitsExpected["nbtag_2ltau_offZ_noOSSF_cut_2"] = 0.958;
1222      upperLimitsExpected["nbtag_3l_onZ_cut_0"] = 32.181;
1223      upperLimitsExpected["nbtag_3l_onZ_cut_1"] = 3.868;
1224      upperLimitsExpected["nbtag_3l_onZ_cut_2"] = 0.887;
1225      upperLimitsExpected["nbtag_2ltau_onZ_cut_0"] = 217.801;
1226      upperLimitsExpected["nbtag_2ltau_onZ_cut_1"] = 9.397;
1227      upperLimitsExpected["nbtag_2ltau_onZ_cut_2"] = 0.787;
1228
1229
1230
1231      if (observed) return upperLimitsObserved[signal_region];
1232      else          return upperLimitsExpected[signal_region];
1233    }
1234
1235
1236    /// Function checking if there is an OSSF lepton pair or a combination of 3 leptons with an invariant mass close to the Z mass
1237    int isonZ (const Particles& particles) {
1238      int onZ = 0;
1239      double best_mass_2 = 999.;
1240      double best_mass_3 = 999.;
1241
1242      // Loop over all 2 particle combinations to find invariant mass of OSSF pair closest to Z mass
1243      for (const Particle& p1 : particles)  {
1244        for (const Particle& p2 : particles)  {
1245          double mass_difference_2_old = fabs(91.0 - best_mass_2);
1246          double mass_difference_2_new = fabs(91.0 - (p1.momentum() + p2.momentum()).mass()/GeV);
1247
1248          // If particle combination is OSSF pair calculate mass difference to Z mass
1249          if ((p1.pid()*p2.pid() == -121 || p1.pid()*p2.pid() == -169))  {
1250
1251            // Get invariant mass closest to Z mass
1252            if (mass_difference_2_new < mass_difference_2_old)
1253              best_mass_2 = (p1.momentum() + p2.momentum()).mass()/GeV;
1254          // In case there is an OSSF pair take also 3rd lepton into account (e.g. from FSR and photon to electron conversion)
1255            for (const Particle& p3 : particles  )  {
1256            double mass_difference_3_old = fabs(91.0 - best_mass_3);
1257            double mass_difference_3_new = fabs(91.0 - (p1.momentum() + p2.momentum() + p3.momentum()).mass()/GeV);
1258            if (mass_difference_3_new < mass_difference_3_old)
1259              best_mass_3 = (p1.momentum() + p2.momentum() + p3.momentum()).mass()/GeV;
1260            }
1261          }
1262        }
1263      }
1264
1265      // Pick the minimum invariant mass of the best OSSF pair combination and the best 3 lepton combination
1266      double best_mass = min(best_mass_2,best_mass_3);
1267      // if this mass is in a 20 GeV window around the Z mass, the event is classified as onZ
1268      if ( fabs(91.0 - best_mass) < 20. ) onZ = 1;
1269      return onZ;
1270    }
1271
1272    /// function checking if two leptons are an OSSF lepton pair and giving out the invariant mass (0 if no OSSF pair)
1273    double isOSSF_mass (const Particle& p1, const Particle& p2) {
1274      double inv_mass = 0.;
1275      // Is particle combination OSSF pair?
1276      if ((p1.pid()*p2.pid() == -121 || p1.pid()*p2.pid() == -169))  {
1277        // Get invariant mass
1278        inv_mass = (p1.momentum() + p2.momentum()).mass()/GeV;
1279      }
1280      return inv_mass;
1281    }
1282
1283    /// Function checking if there is an OSSF lepton pair
1284    bool isOSSF (const Particles& particles)  {
1285      for (size_t i1=0 ; i1 < 3 ; i1 ++)  {
1286        for (size_t i2 = i1+1 ; i2 < 3 ; i2 ++)  {
1287          if ((particles[i1].pid()*particles[i2].pid() == -121 || particles[i1].pid()*particles[i2].pid() == -169))  {
1288            return true;
1289          }
1290        }
1291      }
1292      return false;
1293    }
1294
1295    /// @}
1296
1297  private:
1298
1299    /// Histograms
1300    /// @{
1301    Histo1DPtr _h_HTlep_all, _h_HTjets_all, _h_MET_all, _h_Meff_all, _h_min_pT_all, _h_mT_all;
1302    Histo1DPtr _h_pt_1_3l, _h_pt_2_3l, _h_pt_3_3l, _h_pt_1_2ltau, _h_pt_2_2ltau, _h_pt_3_2ltau;
1303    Histo1DPtr _h_e_n, _h_mu_n, _h_tau_n;
1304    Histo1DPtr _h_excluded;
1305    /// @}
1306
1307    /// Fiducial efficiencies to model the effects of the ATLAS detector
1308    bool _use_fiducial_lepton_efficiency;
1309
1310    /// List of signal regions and event counts per signal region
1311    vector<string> _signal_regions;
1312    map<string, CounterPtr> _eventCountsPerSR;
1313
1314  };
1315
1316
1317  RIVET_DECLARE_PLUGIN(ATLAS_2014_I1327229);
1318
1319}