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