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