rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

ATLAS_2012_I1204447

Inclusive multi-lepton search
Experiment: ATLAS (LHC)
Inspire ID: 1204447
Status: VALIDATED
Authors:
  • Joern Mahlstedt
References:
  • arXiv: 1211.6312
  • Phys. Rev. D 87, 052002 (2013)
Beams: p+ p+
Beam energies: (3500.0, 3500.0) GeV
Run details:
  • Any process producing at least 3 leptons (e.g. pair production of doubly-charged Higgs)

A generic search for anomalous production of events with at least three charged leptons is presented. The search uses a pp-collision data sample at a center-of-mass energy of $\sqrt{s}$ = 7 TeV corresponding to 4.6/fb of integrated luminosity collected in 2011 by the ATLAS detector at the CERN Large Hadron Collider. Events are required to contain at least two electrons or muons, while the third lepton may either be an additional electron or muon, or a hadronically decaying tau lepton. Events are categorized by the presence or absence of a reconstructed tau-lepton or Z-boson candidate decaying to leptons. No significant excess above backgrounds expected from Standard Model processes is observed. Results are presented as upper limits on event yields from non-Standard-Model processes producing at least three prompt, isolated leptons, given as functions of lower bounds on several kinematic variables. Fiducial efficiencies for model testing are also provided. This Rivet module implements the event selection and the fiducial efficiencies to test various models for their exclusion based on observed/excluded limits.

Source code: ATLAS_2012_I1204447.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  class ATLAS_2012_I1204447 : public Analysis {
  15  public:
  16
  17    /// Constructor
  18    ATLAS_2012_I1204447()
  19      : Analysis("ATLAS_2012_I1204447") {  }
  20
  21
  22    /// Book histograms and initialise projections before the run
  23    void init() {
  24
  25      // To calculate the acceptance without having the fiducial lepton efficiencies included, this part can be turned off
  26      _use_fiducial_lepton_efficiency = true;
  27
  28      // Random numbers for simulation of ATLAS detector reconstruction efficiency
  29      srand(160385);
  30
  31      // Read in all signal regions
  32      _signal_regions = getSignalRegions();
  33
  34      // Set number of events per signal region to 0
  35      for (size_t i = 0; i < _signal_regions.size(); i++)
  36        book(_eventCountsPerSR[_signal_regions[i]], "_eventCountsPerSR_" + _signal_regions[i]);
  37
  38      // Final state including all charged and neutral particles
  39      const FinalState fs((Cuts::etaIn(-5.0, 5.0) && Cuts::pT >=  1*GeV));
  40      declare(fs, "FS");
  41
  42      // Final state including all charged particles
  43      declare(ChargedFinalState(Cuts::abseta < 2.5 && Cuts::pT > 1*GeV), "CFS");
  44
  45      // Final state including all visible particles (to calculate MET, Jets etc.)
  46      declare(VisibleFinalState(Cuts::abseta < 5.0), "VFS");
  47
  48      // Final state including all AntiKt 04 Jets
  49      VetoedFinalState vfs;
  50      vfs.addVetoPairId(PID::MUON);
  51      declare(FastJets(vfs, FastJets::ANTIKT, 0.4), "AntiKtJets04");
  52
  53      // Final state including all unstable particles (including taus)
  54      declare(UnstableParticles(Cuts::abseta < 5.0 && Cuts::pT > 5*GeV), "UFS");
  55
  56      // Final state including all electrons
  57      IdentifiedFinalState elecs(Cuts::abseta < 2.47 && Cuts::pT > 10*GeV);
  58      elecs.acceptIdPair(PID::ELECTRON);
  59      declare(elecs, "elecs");
  60
  61      // Final state including all muons
  62      IdentifiedFinalState muons(Cuts::abseta < 2.5 && Cuts::pT > 10*GeV);
  63      muons.acceptIdPair(PID::MUON);
  64      declare(muons, "muons");
  65
  66      // Book histograms
  67      book(_h_HTlep_all  ,"HTlep_all" , 30, 0, 1500);
  68      book(_h_HTjets_all ,"HTjets_all", 30, 0, 1500);
  69      book(_h_MET_all    ,"MET_all"   , 20, 0, 1000);
  70      book(_h_Meff_all   ,"Meff_all"  , 30, 0, 3000);
  71
  72      book(_h_e_n        ,"e_n"  , 10, -0.5, 9.5);
  73      book(_h_mu_n       ,"mu_n" , 10, -0.5, 9.5);
  74      book(_h_tau_n      ,"tau_n", 10, -0.5, 9.5);
  75
  76      book(_h_pt_1_3l    ,"pt_1_3l", 100, 0, 2000);
  77      book(_h_pt_2_3l    ,"pt_2_3l", 100, 0, 2000);
  78      book(_h_pt_3_3l    ,"pt_3_3l", 100, 0, 2000);
  79      book(_h_pt_1_2ltau ,"pt_1_2ltau", 100, 0, 2000);
  80      book(_h_pt_2_2ltau ,"pt_2_2ltau", 100, 0, 2000);
  81      book(_h_pt_3_2ltau ,"pt_3_2ltau", 100, 0, 2000);
  82
  83      book(_h_excluded   ,"excluded", 2, -0.5, 1.5);
  84    }
  85
  86
  87    /// Perform the per-event analysis
  88    void analyze(const Event& event) {
  89
  90      // Muons
  91      Particles muon_candidates;
  92      const Particles charged_tracks    = apply<ChargedFinalState>(event, "CFS").particles();
  93      const Particles visible_particles = apply<VisibleFinalState>(event, "VFS").particles();
  94      for (const Particle& mu : apply<IdentifiedFinalState>(event, "muons").particlesByPt()) {
  95        // Calculate pTCone30 variable (pT of all tracks within dR<0.3 - pT of muon itself)
  96        double pTinCone = -mu.pT();
  97        for (const Particle& track : charged_tracks) {
  98          if (deltaR(mu.momentum(), track.momentum()) < 0.3)
  99            pTinCone += track.pT();
 100        }
 101
 102        // Calculate eTCone30 variable (pT of all visible particles within dR<0.3)
 103        double eTinCone = 0.;
 104        for (const Particle& visible_particle : visible_particles) {
 105          if (visible_particle.abspid() != PID::MUON && inRange(deltaR(mu.momentum(), visible_particle.momentum()), 0.1, 0.3))
 106            eTinCone += visible_particle.pT();
 107        }
 108
 109        // Apply reconstruction efficiency and simulate reco
 110        int muon_id = 13;
 111        if ( mu.hasAncestor(15) || mu.hasAncestor(-15)) muon_id = 14;
 112        const double eff = (_use_fiducial_lepton_efficiency) ? apply_reco_eff(muon_id, mu) : 1.0;
 113        const bool keep_muon = rand()/static_cast<double>(RAND_MAX) <= eff;
 114
 115        // Keep muon if pTCone30/pT < 0.15 and eTCone30/pT < 0.2 and reconstructed
 116        if (keep_muon && pTinCone/mu.pT() <= 0.15 && eTinCone/mu.pT() < 0.2)
 117          muon_candidates.push_back(mu);
 118      }
 119
 120
 121      // Electrons
 122      Particles electron_candidates;
 123      for (const Particle& e : apply<IdentifiedFinalState>(event, "elecs").particlesByPt()) {
 124        // Neglect electrons in crack regions
 125        if (inRange(e.abseta(), 1.37, 1.52)) continue;
 126
 127        // Calculate pTCone30 variable (pT of all tracks within dR<0.3 - pT of electron itself)
 128        double pTinCone = -e.pT();
 129        for (const Particle& track : charged_tracks) {
 130          if (deltaR(e.momentum(), track.momentum()) < 0.3) pTinCone += track.pT();
 131        }
 132
 133        // Calculate eTCone30 variable (pT of all visible particles (except muons) within dR<0.3)
 134        double eTinCone = 0.;
 135        for (const Particle& visible_particle : visible_particles) {
 136          if (visible_particle.abspid() != PID::MUON && inRange(deltaR(e.momentum(), visible_particle.momentum()), 0.1, 0.3))
 137            eTinCone += visible_particle.pT();
 138        }
 139
 140        // Apply reconstruction efficiency and simulate reco
 141        int elec_id = 11;
 142        if (e.hasAncestor(15) || e.hasAncestor(-15)) elec_id = 12;
 143        const double eff = (_use_fiducial_lepton_efficiency) ? apply_reco_eff(elec_id, e) : 1.0;
 144        const bool keep_elec = rand()/static_cast<double>(RAND_MAX) <= eff;
 145
 146        // Keep electron if pTCone30/pT < 0.13 and eTCone30/pT < 0.2 and reconstructed
 147        if (keep_elec && pTinCone/e.pT() <= 0.13 && eTinCone/e.pT() < 0.2)
 148          electron_candidates.push_back(e);
 149      }
 150
 151
 152      // Taus
 153      /// @todo This could benefit from a tau finder projection
 154      Particles tau_candidates;
 155      for (const Particle& tau : apply<UnstableParticles>(event, "UFS").particlesByPt()) {
 156        // Only pick taus out of all unstable particles
 157        if (tau.abspid() != PID::TAU) continue;
 158
 159        // Check that tau has decayed into daughter particles
 160        /// @todo Huh? Unstable taus with no decay vtx? Can use Particle.isStable()? But why in this situation?
 161        if (tau.genParticle()->end_vertex() == 0) continue;
 162
 163        // Calculate visible tau pT from pT of tau neutrino in tau decay for pT and |eta| cuts
 164        FourMomentum daughter_tau_neutrino_momentum = get_tau_neutrino_mom(tau);
 165        Particle tau_vis = tau;
 166        tau_vis.setMomentum(tau.momentum()-daughter_tau_neutrino_momentum);
 167        // keep only taus in certain eta region and above 15 GeV of visible tau pT
 168        if ( tau_vis.pT() <= 15.0*GeV || tau_vis.abseta() > 2.5) continue;
 169
 170        // Get prong number (number of tracks) in tau decay and check if tau decays leptonically
 171        unsigned int nprong = 0;
 172        bool lep_decaying_tau = false;
 173        get_prong_number(tau.genParticle(), nprong, lep_decaying_tau);
 174
 175        // Apply reconstruction efficiency
 176        int tau_id = 15;
 177        if (nprong == 1) tau_id = 15;
 178        else if (nprong == 3) tau_id = 16;
 179
 180        // Get fiducial lepton efficiency simulate reco efficiency
 181        const double eff = (_use_fiducial_lepton_efficiency) ? apply_reco_eff(tau_id, tau_vis) : 1.0;
 182        const bool keep_tau = rand()/static_cast<double>(RAND_MAX) <= eff;
 183
 184        // Keep tau if nprong = 1, it decays hadronically, and it's reconstructed by the detector
 185        if ( !lep_decaying_tau && nprong == 1 && keep_tau) tau_candidates.push_back(tau_vis);
 186      }
 187
 188
 189      // Jets (all anti-kt R=0.4 jets with pT > 25 GeV and eta < 4.9)
 190      Jets jet_candidates;
 191      for (const Jet& jet : apply<FastJets>(event, "AntiKtJets04").jetsByPt(25*GeV)) {
 192        if (jet.abseta() < 4.9) jet_candidates.push_back(jet);
 193      }
 194
 195
 196      // ETmiss
 197      Particles vfs_particles = apply<VisibleFinalState>(event, "VFS").particles();
 198      FourMomentum pTmiss;
 199      for (const Particle& p : vfs_particles) pTmiss -= p.momentum();
 200      double eTmiss = pTmiss.pT()/GeV;
 201
 202
 203      //------------------
 204      // Overlap removal
 205
 206      // electron - electron
 207      Particles electron_candidates_2;
 208      for (size_t ie = 0; ie < electron_candidates.size(); ++ie) {
 209        const Particle & e = electron_candidates[ie];
 210        bool away = true;
 211        // If electron pair within dR < 0.1: remove electron with lower pT
 212        for (size_t ie2=0; ie2 < electron_candidates_2.size(); ++ie2) {
 213          if ( deltaR( e.momentum(), electron_candidates_2[ie2].momentum()) < 0.1 ) {
 214            away = false;
 215            break;
 216          }
 217        }
 218        // If isolated keep it
 219        if ( away )
 220          electron_candidates_2.push_back( e );
 221      }
 222      // jet - electron
 223      Jets recon_jets;
 224      for (const Jet& jet : jet_candidates) {
 225        bool away = true;
 226        // if jet within dR < 0.2 of electron: remove jet
 227        for (const Particle& e : electron_candidates_2) {
 228          if (deltaR(e.momentum(), jet.momentum()) < 0.2) {
 229            away = false;
 230            break;
 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
 249      // electron - jet
 250      Particles recon_leptons, recon_e;
 251      for (size_t ie = 0; ie < electron_candidates_2.size(); ++ie) {
 252        const Particle& e = electron_candidates_2[ie];
 253        // If electron within 0.2 < dR < 0.4 from any jets: remove electron
 254        bool away = true;
 255        for (const Jet& jet : recon_jets) {
 256          if (deltaR(e.momentum(), jet.momentum()) < 0.4) {
 257            away = false;
 258            break;
 259          }
 260        }
 261        // electron - muon
 262        // if electron within dR < 0.1 of a muon: remove electron
 263        if (away) {
 264          for (const Particle& mu : muon_candidates) {
 265            if (deltaR(mu.momentum(), e.momentum()) < 0.1) {
 266              away = false;
 267              break;
 268            }
 269          }
 270        }
 271        // If isolated keep it
 272        if (away)  {
 273          recon_e += e;
 274          recon_leptons += e;
 275        }
 276      }
 277
 278
 279      // tau - electron
 280      Particles recon_tau;
 281      for ( const Particle& tau : tau_candidates ) {
 282        bool away = true;
 283        // If tau within dR < 0.2 of an electron: remove tau
 284        for ( const Particle& e : recon_e ) {
 285          if (deltaR( tau.momentum(), e.momentum()) < 0.2) {
 286            away = false;
 287            break;
 288          }
 289        }
 290        // tau - muon
 291        // If tau within dR < 0.2 of a muon: remove tau
 292        if (away)  {
 293          for (const Particle& mu : muon_candidates) {
 294            if (deltaR(tau.momentum(), mu.momentum()) < 0.2) {
 295              away = false;
 296              break;
 297            }
 298          }
 299        }
 300        // If isolated keep it
 301        if (away) recon_tau.push_back( tau );
 302      }
 303
 304      // Muon - jet isolation
 305      Particles recon_mu, trigger_mu;
 306      // If muon within dR < 0.4 of a jet, remove muon
 307      for (const Particle& mu : muon_candidates) {
 308        bool away = true;
 309        for (const Jet& jet : recon_jets) {
 310          if ( deltaR( mu.momentum(), jet.momentum()) < 0.4 ) {
 311            away = false;
 312            break;
 313          }
 314        }
 315        if (away) {
 316          recon_mu.push_back( mu );
 317          recon_leptons.push_back( mu );
 318          if (mu.abseta() < 2.4) trigger_mu.push_back( mu );
 319        }
 320      }
 321
 322      // End overlap removal
 323      //------------------
 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
 336      // Post-isolation event cuts
 337      // Require at least 3 charged tracks in event
 338      if (charged_tracks.size() < 3) vetoEvent;
 339
 340      // And at least one e/mu passing trigger
 341      if (!( !recon_e   .empty() && recon_e[0]   .pT() > 25*GeV)  &&
 342          !( !trigger_mu.empty() && trigger_mu[0].pT() > 25*GeV) ) {
 343        MSG_DEBUG("Hardest lepton fails trigger");
 344        vetoEvent;
 345      }
 346
 347      // And only accept events with at least 2 electrons and muons and at least 3 leptons in total
 348      if (recon_mu.size() + recon_e.size() + recon_tau.size() < 3 || recon_leptons.size() < 2) vetoEvent;
 349
 350      // Sort leptons by decreasing pT
 351      sortByPt(recon_leptons);
 352      sortByPt(recon_tau);
 353
 354      // Calculate HTlep, fill lepton pT histograms & store chosen combination of 3 leptons
 355      double HTlep = 0.;
 356      Particles chosen_leptons;
 357      if ( recon_leptons.size() > 2 ) {
 358        _h_pt_1_3l->fill(recon_leptons[0].perp()/GeV);
 359        _h_pt_2_3l->fill(recon_leptons[1].perp()/GeV);
 360        _h_pt_3_3l->fill(recon_leptons[2].perp()/GeV);
 361        HTlep = (recon_leptons[0].pT() + recon_leptons[1].pT() + recon_leptons[2].pT())/GeV;
 362        chosen_leptons.push_back( recon_leptons[0] );
 363        chosen_leptons.push_back( recon_leptons[1] );
 364        chosen_leptons.push_back( recon_leptons[2] );
 365      }
 366      else {
 367        _h_pt_1_2ltau->fill(recon_leptons[0].perp()/GeV);
 368        _h_pt_2_2ltau->fill(recon_leptons[1].perp()/GeV);
 369        _h_pt_3_2ltau->fill(recon_tau[0].perp()/GeV);
 370        HTlep = (recon_leptons[0].pT() + recon_leptons[1].pT() + recon_tau[0].pT())/GeV ;
 371        chosen_leptons.push_back( recon_leptons[0] );
 372        chosen_leptons.push_back( recon_leptons[1] );
 373        chosen_leptons.push_back( recon_tau[0] );
 374      }
 375
 376      // Number of prompt e/mu and had taus
 377      _h_e_n  ->fill(recon_e.size());
 378      _h_mu_n ->fill(recon_mu.size());
 379      _h_tau_n->fill(recon_tau.size());
 380
 381      // Calculate HTjets
 382      double HTjets = 0.;
 383      for ( const Jet & jet : recon_jets )
 384        HTjets += jet.perp()/GeV;
 385
 386      // Calculate meff
 387      double meff = eTmiss + HTjets;
 388      Particles all_leptons;
 389      for ( const Particle & e  : recon_e  )  {
 390        meff += e.perp()/GeV;
 391        all_leptons.push_back( e );
 392      }
 393      for ( const Particle & mu : recon_mu )  {
 394        meff += mu.perp()/GeV;
 395        all_leptons.push_back( mu );
 396      }
 397      for ( const Particle & tau : recon_tau )  {
 398        meff += tau.perp()/GeV;
 399        all_leptons.push_back( tau );
 400      }
 401
 402      // Fill histogram of kinematic variables
 403      _h_HTlep_all ->fill(HTlep);
 404      _h_HTjets_all->fill(HTjets);
 405      _h_MET_all   ->fill(eTmiss);
 406      _h_Meff_all  ->fill(meff);
 407
 408      // Determine signal region (3l/2ltau, onZ/offZ)
 409      string basic_signal_region;
 410      if ( recon_mu.size() + recon_e.size() > 2 )
 411        basic_signal_region += "3l_";
 412      else if ( (recon_mu.size() + recon_e.size() == 2) && (recon_tau.size() > 0))
 413        basic_signal_region += "2ltau_";
 414      // Is there an OSSF pair or a three lepton combination with an invariant mass close to the Z mass
 415      int onZ = isonZ(chosen_leptons);
 416      if      (onZ == 1)   basic_signal_region += "onZ";
 417      else if (onZ == 0)   basic_signal_region += "offZ";
 418      // Check in which signal regions this event falls and adjust event counters
 419      fillEventCountsPerSR(basic_signal_region, onZ, HTlep, eTmiss, HTjets, meff);
 420    }
 421
 422
 423    /// Normalise histograms etc., after the run
 424    void finalize() {
 425
 426      // Normalize to an integrated luminosity of 1 fb-1
 427      double norm = crossSection()/femtobarn/sumOfWeights();
 428      string best_signal_region = "";
 429      double ratio_best_SR = 0.;
 430
 431      // Loop over all signal regions and find signal region with best sensitivity (ratio signal events/visible cross-section)
 432      for (size_t i = 0; i < _signal_regions.size(); i++)  {
 433        double signal_events = _eventCountsPerSR[_signal_regions[i]]->val() * norm;
 434        // Use expected upper limits to find best signal region
 435        double UL95  = getUpperLimit(_signal_regions[i], false);
 436        double ratio = signal_events / UL95;
 437        if (ratio > ratio_best_SR)  {
 438          best_signal_region = _signal_regions[i];
 439          ratio_best_SR = ratio;
 440        }
 441      }
 442
 443      double signal_events_best_SR = _eventCountsPerSR[best_signal_region]->val() * norm;
 444      double exp_UL_best_SR = getUpperLimit(best_signal_region, false);
 445      double obs_UL_best_SR = getUpperLimit(best_signal_region, true);
 446
 447      // Print out result
 448      cout << "----------------------------------------------------------------------------------------" << endl;
 449      cout << "Best signal region: " << best_signal_region << endl;
 450      cout << "Normalized number of signal events in this best signal region (per fb-1): " << signal_events_best_SR << endl;
 451      cout << "Efficiency*Acceptance: " <<  _eventCountsPerSR[best_signal_region]->val()/sumOfWeights() << endl;
 452      cout << "Cross-section [fb]: " << crossSection()/femtobarn << endl;
 453      cout << "Expected visible cross-section (per fb-1): " << exp_UL_best_SR << endl;
 454      cout << "Ratio (signal events / expected visible cross-section): " << ratio_best_SR << endl;
 455      cout << "Observed visible cross-section (per fb-1): " << obs_UL_best_SR << endl;
 456      cout << "Ratio (signal events / observed visible cross-section): " <<  signal_events_best_SR/obs_UL_best_SR << endl;
 457      cout << "----------------------------------------------------------------------------------------" << endl;
 458
 459      cout << "Using the EXPECTED limits (visible cross-section) of the analysis: " << endl;
 460      if (signal_events_best_SR > exp_UL_best_SR)  {
 461        cout << "Since the number of signal events > the visible cross-section, this model/grid point is EXCLUDED with 95% CL." << endl;
 462        _h_excluded->fill(1);
 463      }
 464      else  {
 465        cout << "Since the number of signal events < the visible cross-section, this model/grid point is NOT EXCLUDED." << endl;
 466        _h_excluded->fill(0);
 467      }
 468      cout << "----------------------------------------------------------------------------------------" << endl;
 469
 470      cout << "Using the OBSERVED limits (visible cross-section) of the analysis: " << endl;
 471      if (signal_events_best_SR > obs_UL_best_SR)  {
 472        cout << "Since the number of signal events > the visible cross-section, this model/grid point is EXCLUDED with 95% CL." << endl;
 473        _h_excluded->fill(1);
 474      }
 475      else  {
 476        cout << "Since the number of signal events < the visible cross-section, this model/grid point is NOT EXCLUDED." << endl;
 477        _h_excluded->fill(0);
 478      }
 479      cout << "----------------------------------------------------------------------------------------" << endl;
 480
 481
 482      // Normalize to cross section
 483      if (norm != 0)  {
 484        scale(_h_HTlep_all,  norm);
 485        scale(_h_HTjets_all, norm);
 486        scale(_h_MET_all,    norm);
 487        scale(_h_Meff_all,   norm);
 488
 489        scale(_h_pt_1_3l,    norm);
 490        scale(_h_pt_2_3l,    norm);
 491        scale(_h_pt_3_3l,    norm);
 492        scale(_h_pt_1_2ltau, norm);
 493        scale(_h_pt_2_2ltau, norm);
 494        scale(_h_pt_3_2ltau, norm);
 495
 496        scale(_h_e_n,        norm);
 497        scale(_h_mu_n,       norm);
 498        scale(_h_tau_n,      norm);
 499
 500        scale(_h_excluded,   signal_events_best_SR);
 501      }
 502    }
 503
 504
 505    /// Helper functions
 506    //@{
 507
 508    /// Function giving a list of all signal regions
 509    vector<string> getSignalRegions()  {
 510
 511      // List of basic signal regions
 512      vector<string> basic_signal_regions;
 513      basic_signal_regions.push_back("3l_offZ");
 514      basic_signal_regions.push_back("3l_onZ");
 515      basic_signal_regions.push_back("2ltau_offZ");
 516      basic_signal_regions.push_back("2ltau_onZ");
 517
 518      // List of kinematic variables
 519      vector<string> kinematic_variables;
 520      kinematic_variables.push_back("HTlep");
 521      kinematic_variables.push_back("METStrong");
 522      kinematic_variables.push_back("METWeak");
 523      kinematic_variables.push_back("Meff");
 524      kinematic_variables.push_back("MeffStrong");
 525
 526      vector<string> signal_regions;
 527      // Loop over all kinematic variables and basic signal regions
 528      for (size_t i0 = 0; i0 < kinematic_variables.size(); i0++)  {
 529        for (size_t i1 = 0; i1 < basic_signal_regions.size(); i1++)  {
 530          // Is signal region onZ?
 531          int onZ = (basic_signal_regions[i1].find("onZ") != string::npos) ? 1 : 0;
 532          // Get cut values for this kinematic variable
 533          vector<int> cut_values = getCutsPerSignalRegion(kinematic_variables[i0], onZ);
 534          // Loop over all cut values
 535          for (size_t i2 = 0; i2 < cut_values.size(); i2++)  {
 536            // push signal region into vector
 537            signal_regions.push_back( (kinematic_variables[i0] + "_" + basic_signal_regions[i1] + "_cut_" + toString(i2)) );
 538          }
 539        }
 540      }
 541      return signal_regions;
 542    }
 543
 544
 545    /// Function giving all cut vales per kinematic variable (taking onZ for MET into account)
 546    vector<int> getCutsPerSignalRegion(const string& signal_region, int onZ=0)  {
 547      vector<int> cutValues;
 548
 549      // Cut values for HTlep
 550      if (signal_region.compare("HTlep") == 0)  {
 551        cutValues.push_back(0);
 552        cutValues.push_back(100);
 553        cutValues.push_back(150);
 554        cutValues.push_back(200);
 555        cutValues.push_back(300);
 556      }
 557      // Cut values for METStrong (HTjets > 100 GeV) and METWeak (HTjets < 100 GeV)
 558      else if (signal_region.compare("METStrong") == 0 || signal_region.compare("METWeak") == 0)  {
 559        if      (onZ == 0) cutValues.push_back(0);
 560        else if (onZ == 1) cutValues.push_back(20);
 561        cutValues.push_back(50);
 562        cutValues.push_back(75);
 563      }
 564      // Cut values for Meff and MeffStrong (MET > 75 GeV)
 565      if (signal_region.compare("Meff") == 0 || signal_region.compare("MeffStrong") == 0)  {
 566        cutValues.push_back(0);
 567        cutValues.push_back(150);
 568        cutValues.push_back(300);
 569        cutValues.push_back(500);
 570      }
 571
 572      return cutValues;
 573    }
 574
 575
 576    /// function fills map EventCountsPerSR by looping over all signal regions
 577    /// and looking if the event falls into this signal region
 578    void fillEventCountsPerSR(const string& basic_signal_region, int onZ,
 579                              double HTlep, double eTmiss,
 580                              double HTjets, double meff)  {
 581
 582      // Get cut values for HTlep, loop over them and add event if cut is passed
 583      vector<int> cut_values = getCutsPerSignalRegion("HTlep", onZ);
 584      for (size_t i = 0; i < cut_values.size(); i++)  {
 585        if (HTlep > cut_values[i])
 586          _eventCountsPerSR[("HTlep_" + basic_signal_region + "_cut_" + toString(cut_values[i]))]->fill();
 587      }
 588
 589      // Get cut values for METStrong, loop over them and add event if cut is passed
 590      cut_values = getCutsPerSignalRegion("METStrong", onZ);
 591      for (size_t i = 0; i < cut_values.size(); i++)  {
 592        if (eTmiss > cut_values[i] && HTjets > 100.)
 593          _eventCountsPerSR[("METStrong_" + basic_signal_region + "_cut_" + toString(cut_values[i]))]->fill();
 594      }
 595
 596      // Get cut values for METWeak, loop over them and add event if cut is passed
 597      cut_values = getCutsPerSignalRegion("METWeak", onZ);
 598      for (size_t i = 0; i < cut_values.size(); i++)  {
 599        if (eTmiss > cut_values[i] && HTjets <= 100.)
 600          _eventCountsPerSR[("METWeak_" + basic_signal_region + "_cut_" + toString(cut_values[i]))]->fill();
 601      }
 602
 603      // Get cut values for Meff, loop over them and add event if cut is passed
 604      cut_values = getCutsPerSignalRegion("Meff", onZ);
 605      for (size_t i = 0; i < cut_values.size(); i++)  {
 606        if (meff > cut_values[i])
 607          _eventCountsPerSR[("Meff_" + basic_signal_region + "_cut_" + toString(cut_values[i]))]->fill();
 608      }
 609
 610      // Get cut values for MeffStrong, loop over them and add event if cut is passed
 611      cut_values = getCutsPerSignalRegion("MeffStrong", onZ);
 612      for (size_t i = 0; i < cut_values.size(); i++)  {
 613        if (meff > cut_values[i] && eTmiss > 75.)
 614          _eventCountsPerSR[("MeffStrong_" + basic_signal_region + "_cut_" + toString(cut_values[i]))]->fill();
 615      }
 616    }
 617
 618
 619    /// Function returning 4-vector of daughter-particle if it is a tau neutrino
 620    /// @todo Move to TauFinder and make less HepMC-ish
 621    FourMomentum get_tau_neutrino_mom(const Particle& p)  {
 622      assert(p.abspid() == PID::TAU);
 623      ConstGenVertexPtr dv = p.genParticle()->end_vertex();
 624      assert(dv != nullptr);
 625      for(ConstGenParticlePtr pp: HepMCUtils::particles(dv, Relatives::CHILDREN)){
 626        if (abs(pp->pdg_id()) == PID::NU_TAU) return FourMomentum(pp->momentum());
 627      }
 628      return FourMomentum();
 629    }
 630
 631
 632    /// Function calculating the prong number of taus
 633    /// @todo Move to TauFinder and make less HepMC-ish
 634    void get_prong_number(ConstGenParticlePtr p, unsigned int& nprong, bool& lep_decaying_tau) {
 635      assert(p != nullptr);
 636      //const int tau_barcode = p->barcode();
 637      ConstGenVertexPtr dv = p->end_vertex();
 638      assert(dv != nullptr);
 639      for(ConstGenParticlePtr pp: HepMCUtils::particles(dv, Relatives::CHILDREN)){
 640        // If they have status 1 and are charged they will produce a track and the prong number is +1
 641        if (pp->status() == 1 )  {
 642          const int id = pp->pdg_id();
 643          if (Rivet::PID::charge(id) != 0 ) ++nprong;
 644          // Check if tau decays leptonically
 645          // @todo Can a tau decay include a tau in its decay daughters?!
 646          if ((abs(id) == PID::ELECTRON || abs(id) == PID::MUON || abs(id) == PID::TAU) && abs(p->pdg_id()) == PID::TAU) lep_decaying_tau = true;
 647        }
 648        // If the status of the daughter particle is 2 it is unstable and the further decays are checked
 649        else if (pp->status() == 2 )  {
 650          get_prong_number(pp, nprong, lep_decaying_tau);
 651        }
 652      }
 653    }
 654
 655
 656    /// Function giving fiducial lepton efficiency
 657    double apply_reco_eff(int flavor, const Particle& p) {
 658      float pt = p.pT()/GeV;
 659      float eta = p.eta();
 660
 661      double eff = 0.;
 662      //double err = 0.;
 663
 664      if (flavor == 11) { // weight prompt electron -- now including data/MC ID SF in eff.
 665        //float rho = 0.820;
 666        float p0 = 7.34;  float p1 = 0.8977;
 667        //float ep0= 0.5 ;  float ep1= 0.0087;
 668        eff = p1 - p0/pt;
 669
 670        //double err0 = ep0/pt; // d(eff)/dp0
 671        //double err1 = ep1;    // d(eff)/dp1
 672        //err = sqrt(err0*err0 + err1*err1 - 2*rho*err0*err1);
 673
 674        double avgrate = 0.6867;
 675        float wz_ele_eta[] = {0.588717,0.603674,0.666135,0.747493,0.762202,0.675051,0.751606,0.745569,0.665333,0.610432,0.592693,};
 676        //float ewz_ele_eta[] ={0.00292902,0.002476,0.00241209,0.00182319,0.00194339,0.00299785,0.00197339,0.00182004,0.00241793,0.00245997,0.00290394,};
 677        int ibin = 3;
 678
 679        if (eta >= -2.5 && eta < -2.0) ibin = 0;
 680        if (eta >= -2.0 && eta < -1.5) ibin = 1;
 681        if (eta >= -1.5 && eta < -1.0) ibin = 2;
 682        if (eta >= -1.0 && eta < -0.5) ibin = 3;
 683        if (eta >= -0.5 && eta < -0.1) ibin = 4;
 684        if (eta >= -0.1 && eta <  0.1) ibin = 5;
 685        if (eta >=  0.1 && eta <  0.5) ibin = 6;
 686        if (eta >=  0.5 && eta <  1.0) ibin = 7;
 687        if (eta >=  1.0 && eta <  1.5) ibin = 8;
 688        if (eta >=  1.5 && eta <  2.0) ibin = 9;
 689        if (eta >=  2.0 && eta <  2.5) ibin = 10;
 690
 691        double eff_eta = wz_ele_eta[ibin];
 692        //double err_eta = ewz_ele_eta[ibin];
 693
 694        eff = (eff*eff_eta)/avgrate;
 695      }
 696
 697      if (flavor == 12)  { // weight electron from tau
 698        //float rho = 0.884;
 699        float p0 = 6.799;  float p1 = 0.842;
 700        //float ep0= 0.664;  float ep1= 0.016;
 701        eff = p1 - p0/pt;
 702
 703        //double err0 = ep0/pt; // d(eff)/dp0
 704        //double err1 = ep1;    // d(eff)/dp1
 705        //err = sqrt(err0*err0 + err1*err1 - 2*rho*err0*err1);
 706
 707        double avgrate = 0.5319;
 708        float wz_elet_eta[] = {0.468945,0.465953,0.489545,0.58709,0.59669,0.515829,0.59284,0.575828,0.498181,0.463536,0.481738,};
 709        //float ewz_elet_eta[] ={0.00933795,0.00780868,0.00792679,0.00642083,0.00692652,0.0101568,0.00698452,0.00643524,0.0080002,0.00776238,0.0094699,};
 710        int ibin = 3;
 711
 712        if (eta >= -2.5 && eta < -2.0) ibin = 0;
 713        if (eta >= -2.0 && eta < -1.5) ibin = 1;
 714        if (eta >= -1.5 && eta < -1.0) ibin = 2;
 715        if (eta >= -1.0 && eta < -0.5) ibin = 3;
 716        if (eta >= -0.5 && eta < -0.1) ibin = 4;
 717        if (eta >= -0.1 && eta <  0.1) ibin = 5;
 718        if (eta >=  0.1 && eta <  0.5) ibin = 6;
 719        if (eta >=  0.5 && eta <  1.0) ibin = 7;
 720        if (eta >=  1.0 && eta <  1.5) ibin = 8;
 721        if (eta >=  1.5 && eta <  2.0) ibin = 9;
 722        if (eta >=  2.0 && eta <  2.5) ibin = 10;
 723
 724        double eff_eta = wz_elet_eta[ibin];
 725        //double err_eta = ewz_elet_eta[ibin];
 726
 727        eff = (eff*eff_eta)/avgrate;
 728
 729      }
 730
 731      if (flavor == 13)  {// weight prompt muon
 732
 733        //if eta>0.1
 734        float p0 = -18.21;  float p1 = 14.83;  float p2 = 0.9312;
 735        //float ep0= 5.06;    float ep1= 1.9;    float ep2=0.00069;
 736
 737        if ( fabs(eta) < 0.1)  {
 738          p0  = 7.459; p1 = 2.615; p2  = 0.5138;
 739          //ep0 = 10.4; ep1 = 4.934; ep2 = 0.0034;
 740        }
 741
 742        double arg = ( pt-p0 )/( 2.*p1 ) ;
 743        eff = 0.5 * p2 * (1.+erf(arg));
 744        //err = 0.1*eff;
 745      }
 746
 747      if (flavor == 14)  {// weight muon from tau
 748
 749        if (fabs(eta) < 0.1) {
 750          float p0 = -1.756;  float p1 = 12.38;  float p2 = 0.4441;
 751          //float ep0= 10.39;   float ep1= 7.9;  float ep2=0.022;
 752          double arg = ( pt-p0 )/( 2.*p1 ) ;
 753          eff = 0.5 * p2 * (1.+erf(arg));
 754          //err = 0.1*eff;
 755        }
 756        else {
 757          float p0 = 2.102;  float p1 = 0.8293;
 758          //float ep0= 0.271;  float ep1= 0.0083;
 759          eff = p1 - p0/pt;
 760          //double err0 = ep0/pt; // d(eff)/dp0
 761          //double err1 = ep1;    // d(eff)/dp1
 762          //err = sqrt(err0*err0 + err1*err1 - 2*rho*err0*err1);
 763        }
 764      }
 765
 766      if (flavor == 15)  {// weight hadronic tau 1p
 767
 768        float wz_tau1p[] = {0.0249278,0.146978,0.225049,0.229212,0.21519,0.206152,0.201559,0.197917,0.209249,0.228336,0.193548,};
 769        //float ewz_tau1p[] ={0.00178577,0.00425252,0.00535052,0.00592126,0.00484684,0.00612941,0.00792099,0.0083006,0.0138307,0.015568,0.0501751,};
 770        int ibin = 0;
 771        if (pt > 15)  ibin = 1;
 772        if (pt > 20)  ibin = 2;
 773        if (pt > 25)  ibin = 3;
 774        if (pt > 30)  ibin = 4;
 775        if (pt > 40)  ibin = 5;
 776        if (pt > 50)  ibin = 6;
 777        if (pt > 60)  ibin = 7;
 778        if (pt > 80)  ibin = 8;
 779        if (pt > 100) ibin = 9;
 780        if (pt > 200) ibin = 10;
 781
 782        eff = wz_tau1p[ibin];
 783        //err = ewz_tau1p[ibin];
 784
 785
 786        double avgrate = 0.1718;
 787        float wz_tau1p_eta[] = {0.162132,0.176393,0.139619,0.178813,0.185144,0.210027,0.203937,0.178688,0.137034,0.164216,0.163713,};
 788        //float ewz_tau1p_eta[] ={0.00706705,0.00617989,0.00506798,0.00525172,0.00581865,0.00865675,0.00599245,0.00529877,0.00506368,0.00617025,0.00726219,};
 789
 790        ibin = 3;
 791        if (eta >= -2.5 && eta < -2.0) ibin = 0;
 792        if (eta >= -2.0 && eta < -1.5) ibin = 1;
 793        if (eta >= -1.5 && eta < -1.0) ibin = 2;
 794        if (eta >= -1.0 && eta < -0.5) ibin = 3;
 795        if (eta >= -0.5 && eta < -0.1) ibin = 4;
 796        if (eta >= -0.1 && eta <  0.1) ibin = 5;
 797        if (eta >=  0.1 && eta <  0.5) ibin = 6;
 798        if (eta >=  0.5 && eta <  1.0) ibin = 7;
 799        if (eta >=  1.0 && eta <  1.5) ibin = 8;
 800        if (eta >=  1.5 && eta <  2.0) ibin = 9;
 801        if (eta >=  2.0 && eta <  2.5) ibin = 10;
 802
 803        double eff_eta = wz_tau1p_eta[ibin];
 804        //double err_eta = ewz_tau1p_eta[ibin];
 805
 806        eff = (eff*eff_eta)/avgrate;
 807      }
 808
 809      if (flavor == 16)  { //weight hadronic tau 3p
 810
 811        float wz_tau3p[] = {0.000587199,0.00247181,0.0013031,0.00280112,};
 812        //float ewz_tau3p[] ={0.000415091,0.000617187,0.000582385,0.00197792,};
 813
 814        int ibin = 0;
 815        if (pt > 15) ibin = 1;
 816        if (pt > 20) ibin = 2;
 817        if (pt > 40) ibin = 3;
 818        if (pt > 80) ibin = 4;
 819
 820        eff = wz_tau3p[ibin];
 821        //err = ewz_tau3p[ibin];
 822      }
 823
 824      return eff;
 825    }
 826
 827
 828    /// Function giving observed upper limit (visible cross-section)
 829    double getUpperLimit(const string& signal_region, bool observed)  {
 830      map<string,double> upperLimitsObserved;
 831      upperLimitsObserved["HTlep_3l_offZ_cut_0"] = 11.;
 832      upperLimitsObserved["HTlep_3l_offZ_cut_100"] = 8.7;
 833      upperLimitsObserved["HTlep_3l_offZ_cut_150"] = 4.0;
 834      upperLimitsObserved["HTlep_3l_offZ_cut_200"] = 4.4;
 835      upperLimitsObserved["HTlep_3l_offZ_cut_300"] = 1.6;
 836      upperLimitsObserved["HTlep_2ltau_offZ_cut_0"] = 25.;
 837      upperLimitsObserved["HTlep_2ltau_offZ_cut_100"] = 14.;
 838      upperLimitsObserved["HTlep_2ltau_offZ_cut_150"] = 6.1;
 839      upperLimitsObserved["HTlep_2ltau_offZ_cut_200"] = 3.3;
 840      upperLimitsObserved["HTlep_2ltau_offZ_cut_300"] = 1.2;
 841      upperLimitsObserved["HTlep_3l_onZ_cut_0"] = 48.;
 842      upperLimitsObserved["HTlep_3l_onZ_cut_100"] = 38.;
 843      upperLimitsObserved["HTlep_3l_onZ_cut_150"] = 14.;
 844      upperLimitsObserved["HTlep_3l_onZ_cut_200"] = 7.2;
 845      upperLimitsObserved["HTlep_3l_onZ_cut_300"] = 4.5;
 846      upperLimitsObserved["HTlep_2ltau_onZ_cut_0"] = 85.;
 847      upperLimitsObserved["HTlep_2ltau_onZ_cut_100"] = 53.;
 848      upperLimitsObserved["HTlep_2ltau_onZ_cut_150"] = 11.0;
 849      upperLimitsObserved["HTlep_2ltau_onZ_cut_200"] = 5.2;
 850      upperLimitsObserved["HTlep_2ltau_onZ_cut_300"] = 3.0;
 851      upperLimitsObserved["METStrong_3l_offZ_cut_0"] = 2.6;
 852      upperLimitsObserved["METStrong_3l_offZ_cut_50"] = 2.1;
 853      upperLimitsObserved["METStrong_3l_offZ_cut_75"] = 2.1;
 854      upperLimitsObserved["METStrong_2ltau_offZ_cut_0"] = 4.2;
 855      upperLimitsObserved["METStrong_2ltau_offZ_cut_50"] = 3.1;
 856      upperLimitsObserved["METStrong_2ltau_offZ_cut_75"] = 2.6;
 857      upperLimitsObserved["METStrong_3l_onZ_cut_20"] = 11.0;
 858      upperLimitsObserved["METStrong_3l_onZ_cut_50"] = 6.4;
 859      upperLimitsObserved["METStrong_3l_onZ_cut_75"] = 5.1;
 860      upperLimitsObserved["METStrong_2ltau_onZ_cut_20"] = 5.9;
 861      upperLimitsObserved["METStrong_2ltau_onZ_cut_50"] = 3.4;
 862      upperLimitsObserved["METStrong_2ltau_onZ_cut_75"] = 1.2;
 863      upperLimitsObserved["METWeak_3l_offZ_cut_0"] = 11.;
 864      upperLimitsObserved["METWeak_3l_offZ_cut_50"] = 5.3;
 865      upperLimitsObserved["METWeak_3l_offZ_cut_75"] = 3.1;
 866      upperLimitsObserved["METWeak_2ltau_offZ_cut_0"] = 23.;
 867      upperLimitsObserved["METWeak_2ltau_offZ_cut_50"] = 4.3;
 868      upperLimitsObserved["METWeak_2ltau_offZ_cut_75"] = 3.1;
 869      upperLimitsObserved["METWeak_3l_onZ_cut_20"] = 41.;
 870      upperLimitsObserved["METWeak_3l_onZ_cut_50"] = 16.;
 871      upperLimitsObserved["METWeak_3l_onZ_cut_75"] = 8.0;
 872      upperLimitsObserved["METWeak_2ltau_onZ_cut_20"] = 80.;
 873      upperLimitsObserved["METWeak_2ltau_onZ_cut_50"] = 4.4;
 874      upperLimitsObserved["METWeak_2ltau_onZ_cut_75"] = 1.8;
 875      upperLimitsObserved["Meff_3l_offZ_cut_0"] = 11.;
 876      upperLimitsObserved["Meff_3l_offZ_cut_150"] = 8.1;
 877      upperLimitsObserved["Meff_3l_offZ_cut_300"] = 3.1;
 878      upperLimitsObserved["Meff_3l_offZ_cut_500"] = 2.1;
 879      upperLimitsObserved["Meff_2ltau_offZ_cut_0"] = 25.;
 880      upperLimitsObserved["Meff_2ltau_offZ_cut_150"] = 12.;
 881      upperLimitsObserved["Meff_2ltau_offZ_cut_300"] = 3.9;
 882      upperLimitsObserved["Meff_2ltau_offZ_cut_500"] = 2.2;
 883      upperLimitsObserved["Meff_3l_onZ_cut_0"] = 48.;
 884      upperLimitsObserved["Meff_3l_onZ_cut_150"] = 37.;
 885      upperLimitsObserved["Meff_3l_onZ_cut_300"] = 11.;
 886      upperLimitsObserved["Meff_3l_onZ_cut_500"] = 4.8;
 887      upperLimitsObserved["Meff_2ltau_onZ_cut_0"] = 85.;
 888      upperLimitsObserved["Meff_2ltau_onZ_cut_150"] = 28.;
 889      upperLimitsObserved["Meff_2ltau_onZ_cut_300"] = 5.9;
 890      upperLimitsObserved["Meff_2ltau_onZ_cut_500"] = 1.9;
 891      upperLimitsObserved["MeffStrong_3l_offZ_cut_0"] = 3.8;
 892      upperLimitsObserved["MeffStrong_3l_offZ_cut_150"] = 3.8;
 893      upperLimitsObserved["MeffStrong_3l_offZ_cut_300"] = 2.8;
 894      upperLimitsObserved["MeffStrong_3l_offZ_cut_500"] = 2.1;
 895      upperLimitsObserved["MeffStrong_2ltau_offZ_cut_0"] = 3.9;
 896      upperLimitsObserved["MeffStrong_2ltau_offZ_cut_150"] = 4.0;
 897      upperLimitsObserved["MeffStrong_2ltau_offZ_cut_300"] = 2.9;
 898      upperLimitsObserved["MeffStrong_2ltau_offZ_cut_500"] = 1.5;
 899      upperLimitsObserved["MeffStrong_3l_onZ_cut_0"] = 10.0;
 900      upperLimitsObserved["MeffStrong_3l_onZ_cut_150"] = 10.0;
 901      upperLimitsObserved["MeffStrong_3l_onZ_cut_300"] = 6.8;
 902      upperLimitsObserved["MeffStrong_3l_onZ_cut_500"] = 3.9;
 903      upperLimitsObserved["MeffStrong_2ltau_onZ_cut_0"] = 1.6;
 904      upperLimitsObserved["MeffStrong_2ltau_onZ_cut_150"] = 1.4;
 905      upperLimitsObserved["MeffStrong_2ltau_onZ_cut_300"] = 1.5;
 906      upperLimitsObserved["MeffStrong_2ltau_onZ_cut_500"] = 0.9;
 907
 908      // Expected upper limits are also given but not used in this analysis
 909      map<string,double> upperLimitsExpected;
 910      upperLimitsExpected["HTlep_3l_offZ_cut_0"] = 11.;
 911      upperLimitsExpected["HTlep_3l_offZ_cut_100"] = 8.5;
 912      upperLimitsExpected["HTlep_3l_offZ_cut_150"] = 4.6;
 913      upperLimitsExpected["HTlep_3l_offZ_cut_200"] = 3.6;
 914      upperLimitsExpected["HTlep_3l_offZ_cut_300"] = 1.9;
 915      upperLimitsExpected["HTlep_2ltau_offZ_cut_0"] = 23.;
 916      upperLimitsExpected["HTlep_2ltau_offZ_cut_100"] = 14.;
 917      upperLimitsExpected["HTlep_2ltau_offZ_cut_150"] = 6.4;
 918      upperLimitsExpected["HTlep_2ltau_offZ_cut_200"] = 3.6;
 919      upperLimitsExpected["HTlep_2ltau_offZ_cut_300"] = 1.5;
 920      upperLimitsExpected["HTlep_3l_onZ_cut_0"] = 33.;
 921      upperLimitsExpected["HTlep_3l_onZ_cut_100"] = 25.;
 922      upperLimitsExpected["HTlep_3l_onZ_cut_150"] = 12.;
 923      upperLimitsExpected["HTlep_3l_onZ_cut_200"] = 6.5;
 924      upperLimitsExpected["HTlep_3l_onZ_cut_300"] = 3.1;
 925      upperLimitsExpected["HTlep_2ltau_onZ_cut_0"] = 94.;
 926      upperLimitsExpected["HTlep_2ltau_onZ_cut_100"] = 61.;
 927      upperLimitsExpected["HTlep_2ltau_onZ_cut_150"] = 9.9;
 928      upperLimitsExpected["HTlep_2ltau_onZ_cut_200"] = 4.5;
 929      upperLimitsExpected["HTlep_2ltau_onZ_cut_300"] = 1.9;
 930      upperLimitsExpected["METStrong_3l_offZ_cut_0"] = 3.1;
 931      upperLimitsExpected["METStrong_3l_offZ_cut_50"] = 2.4;
 932      upperLimitsExpected["METStrong_3l_offZ_cut_75"] = 2.3;
 933      upperLimitsExpected["METStrong_2ltau_offZ_cut_0"] = 4.8;
 934      upperLimitsExpected["METStrong_2ltau_offZ_cut_50"] = 3.3;
 935      upperLimitsExpected["METStrong_2ltau_offZ_cut_75"] = 2.1;
 936      upperLimitsExpected["METStrong_3l_onZ_cut_20"] = 8.7;
 937      upperLimitsExpected["METStrong_3l_onZ_cut_50"] = 4.9;
 938      upperLimitsExpected["METStrong_3l_onZ_cut_75"] = 3.8;
 939      upperLimitsExpected["METStrong_2ltau_onZ_cut_20"] = 7.3;
 940      upperLimitsExpected["METStrong_2ltau_onZ_cut_50"] = 2.8;
 941      upperLimitsExpected["METStrong_2ltau_onZ_cut_75"] = 1.5;
 942      upperLimitsExpected["METWeak_3l_offZ_cut_0"] = 10.;
 943      upperLimitsExpected["METWeak_3l_offZ_cut_50"] = 4.7;
 944      upperLimitsExpected["METWeak_3l_offZ_cut_75"] = 3.0;
 945      upperLimitsExpected["METWeak_2ltau_offZ_cut_0"] = 21.;
 946      upperLimitsExpected["METWeak_2ltau_offZ_cut_50"] = 4.0;
 947      upperLimitsExpected["METWeak_2ltau_offZ_cut_75"] = 2.6;
 948      upperLimitsExpected["METWeak_3l_onZ_cut_20"] = 30.;
 949      upperLimitsExpected["METWeak_3l_onZ_cut_50"] = 10.;
 950      upperLimitsExpected["METWeak_3l_onZ_cut_75"] = 5.4;
 951      upperLimitsExpected["METWeak_2ltau_onZ_cut_20"] = 88.;
 952      upperLimitsExpected["METWeak_2ltau_onZ_cut_50"] = 5.5;
 953      upperLimitsExpected["METWeak_2ltau_onZ_cut_75"] = 2.2;
 954      upperLimitsExpected["Meff_3l_offZ_cut_0"] = 11.;
 955      upperLimitsExpected["Meff_3l_offZ_cut_150"] = 8.8;
 956      upperLimitsExpected["Meff_3l_offZ_cut_300"] = 3.7;
 957      upperLimitsExpected["Meff_3l_offZ_cut_500"] = 2.1;
 958      upperLimitsExpected["Meff_2ltau_offZ_cut_0"] = 23.;
 959      upperLimitsExpected["Meff_2ltau_offZ_cut_150"] = 13.;
 960      upperLimitsExpected["Meff_2ltau_offZ_cut_300"] = 4.9;
 961      upperLimitsExpected["Meff_2ltau_offZ_cut_500"] = 2.4;
 962      upperLimitsExpected["Meff_3l_onZ_cut_0"] = 33.;
 963      upperLimitsExpected["Meff_3l_onZ_cut_150"] = 25.;
 964      upperLimitsExpected["Meff_3l_onZ_cut_300"] = 9.;
 965      upperLimitsExpected["Meff_3l_onZ_cut_500"] = 3.9;
 966      upperLimitsExpected["Meff_2ltau_onZ_cut_0"] = 94.;
 967      upperLimitsExpected["Meff_2ltau_onZ_cut_150"] = 35.;
 968      upperLimitsExpected["Meff_2ltau_onZ_cut_300"] = 6.8;
 969      upperLimitsExpected["Meff_2ltau_onZ_cut_500"] = 2.5;
 970      upperLimitsExpected["MeffStrong_3l_offZ_cut_0"] = 3.9;
 971      upperLimitsExpected["MeffStrong_3l_offZ_cut_150"] = 3.9;
 972      upperLimitsExpected["MeffStrong_3l_offZ_cut_300"] = 3.0;
 973      upperLimitsExpected["MeffStrong_3l_offZ_cut_500"] = 2.0;
 974      upperLimitsExpected["MeffStrong_2ltau_offZ_cut_0"] = 3.8;
 975      upperLimitsExpected["MeffStrong_2ltau_offZ_cut_150"] = 3.9;
 976      upperLimitsExpected["MeffStrong_2ltau_offZ_cut_300"] = 3.1;
 977      upperLimitsExpected["MeffStrong_2ltau_offZ_cut_500"] = 1.6;
 978      upperLimitsExpected["MeffStrong_3l_onZ_cut_0"] = 6.9;
 979      upperLimitsExpected["MeffStrong_3l_onZ_cut_150"] = 7.1;
 980      upperLimitsExpected["MeffStrong_3l_onZ_cut_300"] = 4.9;
 981      upperLimitsExpected["MeffStrong_3l_onZ_cut_500"] = 3.0;
 982      upperLimitsExpected["MeffStrong_2ltau_onZ_cut_0"] = 2.4;
 983      upperLimitsExpected["MeffStrong_2ltau_onZ_cut_150"] = 2.5;
 984      upperLimitsExpected["MeffStrong_2ltau_onZ_cut_300"] = 2.0;
 985      upperLimitsExpected["MeffStrong_2ltau_onZ_cut_500"] = 1.1;
 986
 987      if (observed) return upperLimitsObserved[signal_region];
 988      else          return upperLimitsExpected[signal_region];
 989    }
 990
 991
 992    /// Function checking if there is an OSSF lepton pair or a combination of 3 leptons with an invariant mass close to the Z mass
 993    /// @todo Should the reference Z mass be 91.2?
 994    int isonZ (const Particles& particles)  {
 995      int onZ = 0;
 996      double best_mass_2 = 999.;
 997      double best_mass_3 = 999.;
 998
 999      // Loop over all 2 particle combinations to find invariant mass of OSSF pair closest to Z mass
1000      for ( const Particle& p1 : particles  )  {
1001        for ( const Particle& p2 : particles  )  {
1002          double mass_difference_2_old = fabs(91.0 - best_mass_2);
1003          double mass_difference_2_new = fabs(91.0 - (p1.momentum() + p2.momentum()).mass()/GeV);
1004
1005          // If particle combination is OSSF pair calculate mass difference to Z mass
1006          if ( (p1.pid()*p2.pid() == -121 || p1.pid()*p2.pid() == -169) )  {
1007
1008            // Get invariant mass closest to Z mass
1009            if (mass_difference_2_new < mass_difference_2_old)
1010              best_mass_2 = (p1.momentum() + p2.momentum()).mass()/GeV;
1011
1012            // In case there is an OSSF pair take also 3rd lepton into account (e.g. from FSR and photon to electron conversion)
1013            for ( const Particle & p3 : particles  )  {
1014              double mass_difference_3_old = fabs(91.0 - best_mass_3);
1015              double mass_difference_3_new = fabs(91.0 - (p1.momentum() + p2.momentum() + p3.momentum()).mass()/GeV);
1016              if (mass_difference_3_new < mass_difference_3_old)
1017                best_mass_3 = (p1.momentum() + p2.momentum() + p3.momentum()).mass()/GeV;
1018            }
1019          }
1020        }
1021      }
1022
1023      // Pick the minimum invariant mass of the best OSSF pair combination and the best 3 lepton combination
1024      // If this mass is in a 20 GeV window around the Z mass, the event is classified as onZ
1025      double best_mass = min(best_mass_2, best_mass_3);
1026      if (fabs(91.0 - best_mass) < 20) onZ = 1;
1027      return onZ;
1028    }
1029
1030    //@}
1031
1032
1033  private:
1034
1035    /// Histograms
1036    //@{
1037    Histo1DPtr _h_HTlep_all, _h_HTjets_all, _h_MET_all, _h_Meff_all;
1038    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;
1039    Histo1DPtr _h_e_n, _h_mu_n, _h_tau_n;
1040    Histo1DPtr _h_excluded;
1041    //@}
1042
1043    /// Fiducial efficiencies to model the effects of the ATLAS detector
1044    bool _use_fiducial_lepton_efficiency;
1045
1046    /// List of signal regions and event counts per signal region
1047    vector<string> _signal_regions;
1048
1049    map<string, CounterPtr> _eventCountsPerSR;
1050
1051  };
1052
1053
1054
1055  RIVET_DECLARE_PLUGIN(ATLAS_2012_I1204447);
1056
1057}