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: 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, JetAlg::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.hasAncestorWith(Cuts::pid == 15) || mu.hasAncestorWith(Cuts::pid == -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.hasAncestorWith(Cuts::pid == 15) || e.hasAncestorWith(Cuts::pid == -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 = apply<FastJets>(event, "AntiKtJets04").jetsByPt(Cuts::pT > 25*GeV && Cuts::abseta < 4.9);
 191
 192      // ETmiss
 193      Particles vfs_particles = apply<VisibleFinalState>(event, "VFS").particles();
 194      FourMomentum pTmiss;
 195      for (const Particle& p : vfs_particles) pTmiss -= p.momentum();
 196      double eTmiss = pTmiss.pT()/GeV;
 197
 198
 199      //------------------
 200      // Overlap removal
 201
 202      // electron - electron
 203      Particles electron_candidates_2;
 204      for (size_t ie = 0; ie < electron_candidates.size(); ++ie) {
 205        const Particle & e = electron_candidates[ie];
 206        bool away = true;
 207        // If electron pair within dR < 0.1: remove electron with lower pT
 208        for (size_t ie2=0; ie2 < electron_candidates_2.size(); ++ie2) {
 209          if ( deltaR( e.momentum(), electron_candidates_2[ie2].momentum()) < 0.1 ) {
 210            away = false;
 211            break;
 212          }
 213        }
 214        // If isolated keep it
 215        if ( away )
 216          electron_candidates_2.push_back( e );
 217      }
 218      // jet - electron
 219      Jets recon_jets;
 220      for (const Jet& jet : jet_candidates) {
 221        bool away = true;
 222        // if jet within dR < 0.2 of electron: remove jet
 223        for (const Particle& e : electron_candidates_2) {
 224          if (deltaR(e.momentum(), jet.momentum()) < 0.2) {
 225            away = false;
 226            break;
 227          }
 228        }
 229        // jet - tau
 230        if (away)  {
 231          // If jet within dR < 0.2 of tau: remove jet
 232          for (const Particle& tau : tau_candidates) {
 233            if (deltaR(tau.momentum(), jet.momentum()) < 0.2) {
 234              away = false;
 235              break;
 236            }
 237          }
 238        }
 239        // If isolated keep it
 240        if ( away )
 241          recon_jets.push_back( jet );
 242      }
 243
 244
 245      // electron - jet
 246      Particles recon_leptons, recon_e;
 247      for (size_t ie = 0; ie < electron_candidates_2.size(); ++ie) {
 248        const Particle& e = electron_candidates_2[ie];
 249        // If electron within 0.2 < dR < 0.4 from any jets: remove electron
 250        bool away = true;
 251        for (const Jet& jet : recon_jets) {
 252          if (deltaR(e.momentum(), jet.momentum()) < 0.4) {
 253            away = false;
 254            break;
 255          }
 256        }
 257        // electron - muon
 258        // if electron within dR < 0.1 of a muon: remove electron
 259        if (away) {
 260          for (const Particle& mu : muon_candidates) {
 261            if (deltaR(mu.momentum(), e.momentum()) < 0.1) {
 262              away = false;
 263              break;
 264            }
 265          }
 266        }
 267        // If isolated keep it
 268        if (away)  {
 269          recon_e += e;
 270          recon_leptons += e;
 271        }
 272      }
 273
 274
 275      // tau - electron
 276      Particles recon_tau;
 277      for ( const Particle& tau : tau_candidates ) {
 278        bool away = true;
 279        // If tau within dR < 0.2 of an electron: remove tau
 280        for ( const Particle& e : recon_e ) {
 281          if (deltaR( tau.momentum(), e.momentum()) < 0.2) {
 282            away = false;
 283            break;
 284          }
 285        }
 286        // tau - muon
 287        // If tau within dR < 0.2 of a muon: remove tau
 288        if (away)  {
 289          for (const Particle& mu : muon_candidates) {
 290            if (deltaR(tau.momentum(), mu.momentum()) < 0.2) {
 291              away = false;
 292              break;
 293            }
 294          }
 295        }
 296        // If isolated keep it
 297        if (away) recon_tau.push_back( tau );
 298      }
 299
 300      // Muon - jet isolation
 301      Particles recon_mu, trigger_mu;
 302      // If muon within dR < 0.4 of a jet, remove muon
 303      for (const Particle& mu : muon_candidates) {
 304        bool away = true;
 305        for (const Jet& jet : recon_jets) {
 306          if ( deltaR( mu.momentum(), jet.momentum()) < 0.4 ) {
 307            away = false;
 308            break;
 309          }
 310        }
 311        if (away) {
 312          recon_mu.push_back( mu );
 313          recon_leptons.push_back( mu );
 314          if (mu.abseta() < 2.4) trigger_mu.push_back( mu );
 315        }
 316      }
 317
 318      // End overlap removal
 319      //------------------
 320
 321
 322      // Jet cleaning
 323      if (rand()/static_cast<double>(RAND_MAX) <= 0.42) {
 324        for (const Jet& jet : recon_jets) {
 325          const double eta = jet.rapidity();
 326          const double phi = jet.azimuthalAngle(MINUSPI_PLUSPI);
 327          if (jet.pT() > 25*GeV && inRange(eta, -0.1, 1.5) && inRange(phi, -0.9, -0.5)) vetoEvent;
 328        }
 329      }
 330
 331
 332      // Post-isolation event cuts
 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() > 25*GeV)  &&
 338          !( !trigger_mu.empty() && trigger_mu[0].pT() > 25*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].perp()/GeV);
 355        _h_pt_2_3l->fill(recon_leptons[1].perp()/GeV);
 356        _h_pt_3_3l->fill(recon_leptons[2].perp()/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].perp()/GeV);
 364        _h_pt_2_2ltau->fill(recon_leptons[1].perp()/GeV);
 365        _h_pt_3_2ltau->fill(recon_tau[0].perp()/GeV);
 366        HTlep = (recon_leptons[0].pT() + recon_leptons[1].pT() + 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      // Number of prompt e/mu and had taus
 373      _h_e_n  ->fill(recon_e.size());
 374      _h_mu_n ->fill(recon_mu.size());
 375      _h_tau_n->fill(recon_tau.size());
 376
 377      // Calculate HTjets
 378      double HTjets = 0.;
 379      for ( const Jet & jet : recon_jets )
 380        HTjets += jet.perp()/GeV;
 381
 382      // Calculate meff
 383      double meff = eTmiss + HTjets;
 384      Particles all_leptons;
 385      for ( const Particle & e  : recon_e  )  {
 386        meff += e.perp()/GeV;
 387        all_leptons.push_back( e );
 388      }
 389      for ( const Particle & mu : recon_mu )  {
 390        meff += mu.perp()/GeV;
 391        all_leptons.push_back( mu );
 392      }
 393      for ( const Particle & tau : recon_tau )  {
 394        meff += tau.perp()/GeV;
 395        all_leptons.push_back( tau );
 396      }
 397
 398      // Fill histogram of kinematic variables
 399      _h_HTlep_all ->fill(HTlep);
 400      _h_HTjets_all->fill(HTjets);
 401      _h_MET_all   ->fill(eTmiss);
 402      _h_Meff_all  ->fill(meff);
 403
 404      // Determine signal region (3l/2ltau, onZ/offZ)
 405      string basic_signal_region;
 406      if ( recon_mu.size() + recon_e.size() > 2 )
 407        basic_signal_region += "3l_";
 408      else if ( (recon_mu.size() + recon_e.size() == 2) && (recon_tau.size() > 0))
 409        basic_signal_region += "2ltau_";
 410      // Is there an OSSF pair or a three lepton combination with an invariant mass close to the Z mass
 411      int onZ = isonZ(chosen_leptons);
 412      if      (onZ == 1)   basic_signal_region += "onZ";
 413      else if (onZ == 0)   basic_signal_region += "offZ";
 414      // Check in which signal regions this event falls and adjust event counters
 415      fillEventCountsPerSR(basic_signal_region, onZ, HTlep, eTmiss, HTjets, meff);
 416    }
 417
 418
 419    /// Normalise histograms etc., after the run
 420    void finalize() {
 421
 422      // Normalize to an integrated luminosity of 1 fb-1
 423      double norm = crossSection()/femtobarn/sumOfWeights();
 424      string best_signal_region = "";
 425      double ratio_best_SR = 0.;
 426
 427      // Loop over all signal regions and find signal region with best sensitivity (ratio signal events/visible cross-section)
 428      for (size_t i = 0; i < _signal_regions.size(); i++)  {
 429        double signal_events = _eventCountsPerSR[_signal_regions[i]]->val() * norm;
 430        // Use expected upper limits to find best signal region
 431        double UL95  = getUpperLimit(_signal_regions[i], false);
 432        double ratio = signal_events / UL95;
 433        if (ratio > ratio_best_SR)  {
 434          best_signal_region = _signal_regions[i];
 435          ratio_best_SR = ratio;
 436        }
 437      }
 438
 439      double signal_events_best_SR = _eventCountsPerSR[best_signal_region]->val() * norm;
 440      double exp_UL_best_SR = getUpperLimit(best_signal_region, false);
 441      double obs_UL_best_SR = getUpperLimit(best_signal_region, true);
 442
 443      // Print out result
 444      cout << "----------------------------------------------------------------------------------------" << endl;
 445      cout << "Best signal region: " << best_signal_region << endl;
 446      cout << "Normalized number of signal events in this best signal region (per fb-1): " << signal_events_best_SR << endl;
 447      cout << "Efficiency*Acceptance: " <<  _eventCountsPerSR[best_signal_region]->val()/sumOfWeights() << endl;
 448      cout << "Cross-section [fb]: " << crossSection()/femtobarn << endl;
 449      cout << "Expected visible cross-section (per fb-1): " << exp_UL_best_SR << endl;
 450      cout << "Ratio (signal events / expected visible cross-section): " << ratio_best_SR << endl;
 451      cout << "Observed visible cross-section (per fb-1): " << obs_UL_best_SR << endl;
 452      cout << "Ratio (signal events / observed visible cross-section): " <<  signal_events_best_SR/obs_UL_best_SR << endl;
 453      cout << "----------------------------------------------------------------------------------------" << endl;
 454
 455      cout << "Using the EXPECTED limits (visible cross-section) of the analysis: " << endl;
 456      if (signal_events_best_SR > exp_UL_best_SR)  {
 457        cout << "Since the number of signal events > the visible cross-section, this model/grid point is EXCLUDED with 95% CL." << endl;
 458        _h_excluded->fill(1);
 459      }
 460      else  {
 461        cout << "Since the number of signal events < the visible cross-section, this model/grid point is NOT EXCLUDED." << endl;
 462        _h_excluded->fill(0);
 463      }
 464      cout << "----------------------------------------------------------------------------------------" << endl;
 465
 466      cout << "Using the OBSERVED limits (visible cross-section) of the analysis: " << endl;
 467      if (signal_events_best_SR > obs_UL_best_SR)  {
 468        cout << "Since the number of signal events > the visible cross-section, this model/grid point is EXCLUDED with 95% CL." << endl;
 469        _h_excluded->fill(1);
 470      }
 471      else  {
 472        cout << "Since the number of signal events < the visible cross-section, this model/grid point is NOT EXCLUDED." << endl;
 473        _h_excluded->fill(0);
 474      }
 475      cout << "----------------------------------------------------------------------------------------" << endl;
 476
 477
 478      // Normalize to cross section
 479      if (norm != 0)  {
 480        scale(_h_HTlep_all,  norm);
 481        scale(_h_HTjets_all, norm);
 482        scale(_h_MET_all,    norm);
 483        scale(_h_Meff_all,   norm);
 484
 485        scale(_h_pt_1_3l,    norm);
 486        scale(_h_pt_2_3l,    norm);
 487        scale(_h_pt_3_3l,    norm);
 488        scale(_h_pt_1_2ltau, norm);
 489        scale(_h_pt_2_2ltau, norm);
 490        scale(_h_pt_3_2ltau, norm);
 491
 492        scale(_h_e_n,        norm);
 493        scale(_h_mu_n,       norm);
 494        scale(_h_tau_n,      norm);
 495
 496        scale(_h_excluded,   signal_events_best_SR);
 497      }
 498    }
 499
 500
 501    /// Helper functions
 502    /// @{
 503
 504    /// Function giving a list of all signal regions
 505    vector<string> getSignalRegions()  {
 506
 507      // List of basic signal regions
 508      vector<string> basic_signal_regions;
 509      basic_signal_regions.push_back("3l_offZ");
 510      basic_signal_regions.push_back("3l_onZ");
 511      basic_signal_regions.push_back("2ltau_offZ");
 512      basic_signal_regions.push_back("2ltau_onZ");
 513
 514      // List of kinematic variables
 515      vector<string> kinematic_variables;
 516      kinematic_variables.push_back("HTlep");
 517      kinematic_variables.push_back("METStrong");
 518      kinematic_variables.push_back("METWeak");
 519      kinematic_variables.push_back("Meff");
 520      kinematic_variables.push_back("MeffStrong");
 521
 522      vector<string> signal_regions;
 523      // Loop over all kinematic variables and basic signal regions
 524      for (size_t i0 = 0; i0 < kinematic_variables.size(); i0++)  {
 525        for (size_t i1 = 0; i1 < basic_signal_regions.size(); i1++)  {
 526          // Is signal region onZ?
 527          int onZ = (basic_signal_regions[i1].find("onZ") != string::npos) ? 1 : 0;
 528          // Get cut values for this kinematic variable
 529          vector<int> cut_values = getCutsPerSignalRegion(kinematic_variables[i0], onZ);
 530          // Loop over all cut values
 531          for (size_t i2 = 0; i2 < cut_values.size(); i2++)  {
 532            // push signal region into vector
 533            signal_regions.push_back( (kinematic_variables[i0] + "_" + basic_signal_regions[i1] + "_cut_" + toString(i2)) );
 534          }
 535        }
 536      }
 537      return signal_regions;
 538    }
 539
 540
 541    /// Function giving all cut vales per kinematic variable (taking onZ for MET into account)
 542    vector<int> getCutsPerSignalRegion(const string& signal_region, int onZ=0)  {
 543      vector<int> cutValues;
 544
 545      // Cut values for HTlep
 546      if (signal_region.compare("HTlep") == 0)  {
 547        cutValues.push_back(0);
 548        cutValues.push_back(100);
 549        cutValues.push_back(150);
 550        cutValues.push_back(200);
 551        cutValues.push_back(300);
 552      }
 553      // Cut values for METStrong (HTjets > 100 GeV) and METWeak (HTjets < 100 GeV)
 554      else if (signal_region.compare("METStrong") == 0 || signal_region.compare("METWeak") == 0)  {
 555        if      (onZ == 0) cutValues.push_back(0);
 556        else if (onZ == 1) cutValues.push_back(20);
 557        cutValues.push_back(50);
 558        cutValues.push_back(75);
 559      }
 560      // Cut values for Meff and MeffStrong (MET > 75 GeV)
 561      if (signal_region.compare("Meff") == 0 || signal_region.compare("MeffStrong") == 0)  {
 562        cutValues.push_back(0);
 563        cutValues.push_back(150);
 564        cutValues.push_back(300);
 565        cutValues.push_back(500);
 566      }
 567
 568      return cutValues;
 569    }
 570
 571
 572    /// function fills map EventCountsPerSR by looping over all signal regions
 573    /// and looking if the event falls into this signal region
 574    void fillEventCountsPerSR(const string& basic_signal_region, int onZ,
 575                              double HTlep, double eTmiss,
 576                              double HTjets, double meff)  {
 577
 578      // Get cut values for HTlep, loop over them and add event if cut is passed
 579      vector<int> cut_values = getCutsPerSignalRegion("HTlep", onZ);
 580      for (size_t i = 0; i < cut_values.size(); i++)  {
 581        if (HTlep > cut_values[i])
 582          _eventCountsPerSR[("HTlep_" + basic_signal_region + "_cut_" + toString(cut_values[i]))]->fill();
 583      }
 584
 585      // Get cut values for METStrong, loop over them and add event if cut is passed
 586      cut_values = getCutsPerSignalRegion("METStrong", onZ);
 587      for (size_t i = 0; i < cut_values.size(); i++)  {
 588        if (eTmiss > cut_values[i] && HTjets > 100.)
 589          _eventCountsPerSR[("METStrong_" + basic_signal_region + "_cut_" + toString(cut_values[i]))]->fill();
 590      }
 591
 592      // Get cut values for METWeak, loop over them and add event if cut is passed
 593      cut_values = getCutsPerSignalRegion("METWeak", onZ);
 594      for (size_t i = 0; i < cut_values.size(); i++)  {
 595        if (eTmiss > cut_values[i] && HTjets <= 100.)
 596          _eventCountsPerSR[("METWeak_" + basic_signal_region + "_cut_" + toString(cut_values[i]))]->fill();
 597      }
 598
 599      // Get cut values for Meff, loop over them and add event if cut is passed
 600      cut_values = getCutsPerSignalRegion("Meff", onZ);
 601      for (size_t i = 0; i < cut_values.size(); i++)  {
 602        if (meff > cut_values[i])
 603          _eventCountsPerSR[("Meff_" + basic_signal_region + "_cut_" + toString(cut_values[i]))]->fill();
 604      }
 605
 606      // Get cut values for MeffStrong, loop over them and add event if cut is passed
 607      cut_values = getCutsPerSignalRegion("MeffStrong", onZ);
 608      for (size_t i = 0; i < cut_values.size(); i++)  {
 609        if (meff > cut_values[i] && eTmiss > 75.)
 610          _eventCountsPerSR[("MeffStrong_" + basic_signal_region + "_cut_" + toString(cut_values[i]))]->fill();
 611      }
 612    }
 613
 614
 615    /// Function returning 4-vector of daughter-particle if it is a tau neutrino
 616    /// @todo Move to TauFinder and make less HepMC-ish
 617    FourMomentum get_tau_neutrino_mom(const Particle& p)  {
 618      assert(p.abspid() == PID::TAU);
 619      ConstGenVertexPtr dv = p.genParticle()->end_vertex();
 620      assert(dv != nullptr);
 621      for(ConstGenParticlePtr pp: HepMCUtils::particles(dv, Relatives::CHILDREN)){
 622        if (abs(pp->pdg_id()) == PID::NU_TAU) return FourMomentum(pp->momentum());
 623      }
 624      return FourMomentum();
 625    }
 626
 627
 628    /// Function calculating the prong number of taus
 629    /// @todo Move to TauFinder and make less HepMC-ish
 630    void get_prong_number(ConstGenParticlePtr p, unsigned int& nprong, bool& lep_decaying_tau) {
 631      assert(p != nullptr);
 632      //const int tau_barcode = p->barcode();
 633      ConstGenVertexPtr dv = p->end_vertex();
 634      assert(dv != nullptr);
 635      for(ConstGenParticlePtr pp: HepMCUtils::particles(dv, Relatives::CHILDREN)){
 636        // If they have status 1 and are charged they will produce a track and the prong number is +1
 637        if (pp->status() == 1 )  {
 638          const int id = pp->pdg_id();
 639          if (Rivet::PID::charge(id) != 0 ) ++nprong;
 640          // Check if tau decays leptonically
 641          // @todo Can a tau decay include a tau in its decay daughters?!
 642          if ((abs(id) == PID::ELECTRON || abs(id) == PID::MUON || abs(id) == PID::TAU) && abs(p->pdg_id()) == PID::TAU) lep_decaying_tau = true;
 643        }
 644        // If the status of the daughter particle is 2 it is unstable and the further decays are checked
 645        else if (pp->status() == 2 )  {
 646          get_prong_number(pp, nprong, lep_decaying_tau);
 647        }
 648      }
 649    }
 650
 651
 652    /// Function giving fiducial lepton efficiency
 653    double apply_reco_eff(int flavor, const Particle& p) {
 654      float pt = p.pT()/GeV;
 655      float eta = p.eta();
 656
 657      double eff = 0.;
 658      //double err = 0.;
 659
 660      if (flavor == 11) { // weight prompt electron -- now including data/MC ID SF in eff.
 661        //float rho = 0.820;
 662        float p0 = 7.34;  float p1 = 0.8977;
 663        //float ep0= 0.5 ;  float ep1= 0.0087;
 664        eff = p1 - p0/pt;
 665
 666        //double err0 = ep0/pt; // d(eff)/dp0
 667        //double err1 = ep1;    // d(eff)/dp1
 668        //err = sqrt(err0*err0 + err1*err1 - 2*rho*err0*err1);
 669
 670        double avgrate = 0.6867;
 671        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,};
 672        //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,};
 673        int ibin = 3;
 674
 675        if (eta >= -2.5 && eta < -2.0) ibin = 0;
 676        if (eta >= -2.0 && eta < -1.5) ibin = 1;
 677        if (eta >= -1.5 && eta < -1.0) ibin = 2;
 678        if (eta >= -1.0 && eta < -0.5) ibin = 3;
 679        if (eta >= -0.5 && eta < -0.1) ibin = 4;
 680        if (eta >= -0.1 && eta <  0.1) ibin = 5;
 681        if (eta >=  0.1 && eta <  0.5) ibin = 6;
 682        if (eta >=  0.5 && eta <  1.0) ibin = 7;
 683        if (eta >=  1.0 && eta <  1.5) ibin = 8;
 684        if (eta >=  1.5 && eta <  2.0) ibin = 9;
 685        if (eta >=  2.0 && eta <  2.5) ibin = 10;
 686
 687        double eff_eta = wz_ele_eta[ibin];
 688        //double err_eta = ewz_ele_eta[ibin];
 689
 690        eff = (eff*eff_eta)/avgrate;
 691      }
 692
 693      if (flavor == 12)  { // weight electron from tau
 694        //float rho = 0.884;
 695        float p0 = 6.799;  float p1 = 0.842;
 696        //float ep0= 0.664;  float ep1= 0.016;
 697        eff = p1 - p0/pt;
 698
 699        //double err0 = ep0/pt; // d(eff)/dp0
 700        //double err1 = ep1;    // d(eff)/dp1
 701        //err = sqrt(err0*err0 + err1*err1 - 2*rho*err0*err1);
 702
 703        double avgrate = 0.5319;
 704        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,};
 705        //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,};
 706        int ibin = 3;
 707
 708        if (eta >= -2.5 && eta < -2.0) ibin = 0;
 709        if (eta >= -2.0 && eta < -1.5) ibin = 1;
 710        if (eta >= -1.5 && eta < -1.0) ibin = 2;
 711        if (eta >= -1.0 && eta < -0.5) ibin = 3;
 712        if (eta >= -0.5 && eta < -0.1) ibin = 4;
 713        if (eta >= -0.1 && eta <  0.1) ibin = 5;
 714        if (eta >=  0.1 && eta <  0.5) ibin = 6;
 715        if (eta >=  0.5 && eta <  1.0) ibin = 7;
 716        if (eta >=  1.0 && eta <  1.5) ibin = 8;
 717        if (eta >=  1.5 && eta <  2.0) ibin = 9;
 718        if (eta >=  2.0 && eta <  2.5) ibin = 10;
 719
 720        double eff_eta = wz_elet_eta[ibin];
 721        //double err_eta = ewz_elet_eta[ibin];
 722
 723        eff = (eff*eff_eta)/avgrate;
 724
 725      }
 726
 727      if (flavor == 13)  {// weight prompt muon
 728
 729        //if eta>0.1
 730        float p0 = -18.21;  float p1 = 14.83;  float p2 = 0.9312;
 731        //float ep0= 5.06;    float ep1= 1.9;    float ep2=0.00069;
 732
 733        if ( fabs(eta) < 0.1)  {
 734          p0  = 7.459; p1 = 2.615; p2  = 0.5138;
 735          //ep0 = 10.4; ep1 = 4.934; ep2 = 0.0034;
 736        }
 737
 738        double arg = ( pt-p0 )/( 2.*p1 ) ;
 739        eff = 0.5 * p2 * (1.+erf(arg));
 740        //err = 0.1*eff;
 741      }
 742
 743      if (flavor == 14)  {// weight muon from tau
 744
 745        if (fabs(eta) < 0.1) {
 746          float p0 = -1.756;  float p1 = 12.38;  float p2 = 0.4441;
 747          //float ep0= 10.39;   float ep1= 7.9;  float ep2=0.022;
 748          double arg = ( pt-p0 )/( 2.*p1 ) ;
 749          eff = 0.5 * p2 * (1.+erf(arg));
 750          //err = 0.1*eff;
 751        }
 752        else {
 753          float p0 = 2.102;  float p1 = 0.8293;
 754          //float ep0= 0.271;  float ep1= 0.0083;
 755          eff = p1 - p0/pt;
 756          //double err0 = ep0/pt; // d(eff)/dp0
 757          //double err1 = ep1;    // d(eff)/dp1
 758          //err = sqrt(err0*err0 + err1*err1 - 2*rho*err0*err1);
 759        }
 760      }
 761
 762      if (flavor == 15)  {// weight hadronic tau 1p
 763
 764        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,};
 765        //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,};
 766        int ibin = 0;
 767        if (pt > 15)  ibin = 1;
 768        if (pt > 20)  ibin = 2;
 769        if (pt > 25)  ibin = 3;
 770        if (pt > 30)  ibin = 4;
 771        if (pt > 40)  ibin = 5;
 772        if (pt > 50)  ibin = 6;
 773        if (pt > 60)  ibin = 7;
 774        if (pt > 80)  ibin = 8;
 775        if (pt > 100) ibin = 9;
 776        if (pt > 200) ibin = 10;
 777
 778        eff = wz_tau1p[ibin];
 779        //err = ewz_tau1p[ibin];
 780
 781
 782        double avgrate = 0.1718;
 783        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,};
 784        //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,};
 785
 786        ibin = 3;
 787        if (eta >= -2.5 && eta < -2.0) ibin = 0;
 788        if (eta >= -2.0 && eta < -1.5) ibin = 1;
 789        if (eta >= -1.5 && eta < -1.0) ibin = 2;
 790        if (eta >= -1.0 && eta < -0.5) ibin = 3;
 791        if (eta >= -0.5 && eta < -0.1) ibin = 4;
 792        if (eta >= -0.1 && eta <  0.1) ibin = 5;
 793        if (eta >=  0.1 && eta <  0.5) ibin = 6;
 794        if (eta >=  0.5 && eta <  1.0) ibin = 7;
 795        if (eta >=  1.0 && eta <  1.5) ibin = 8;
 796        if (eta >=  1.5 && eta <  2.0) ibin = 9;
 797        if (eta >=  2.0 && eta <  2.5) ibin = 10;
 798
 799        double eff_eta = wz_tau1p_eta[ibin];
 800        //double err_eta = ewz_tau1p_eta[ibin];
 801
 802        eff = (eff*eff_eta)/avgrate;
 803      }
 804
 805      if (flavor == 16)  { //weight hadronic tau 3p
 806
 807        float wz_tau3p[] = {0.000587199,0.00247181,0.0013031,0.00280112,};
 808        //float ewz_tau3p[] ={0.000415091,0.000617187,0.000582385,0.00197792,};
 809
 810        int ibin = 0;
 811        if (pt > 15) ibin = 1;
 812        if (pt > 20) ibin = 2;
 813        if (pt > 40) ibin = 3;
 814        if (pt > 80) ibin = 4;
 815
 816        eff = wz_tau3p[ibin];
 817        //err = ewz_tau3p[ibin];
 818      }
 819
 820      return eff;
 821    }
 822
 823
 824    /// Function giving observed upper limit (visible cross-section)
 825    double getUpperLimit(const string& signal_region, bool observed)  {
 826      map<string,double> upperLimitsObserved;
 827      upperLimitsObserved["HTlep_3l_offZ_cut_0"] = 11.;
 828      upperLimitsObserved["HTlep_3l_offZ_cut_100"] = 8.7;
 829      upperLimitsObserved["HTlep_3l_offZ_cut_150"] = 4.0;
 830      upperLimitsObserved["HTlep_3l_offZ_cut_200"] = 4.4;
 831      upperLimitsObserved["HTlep_3l_offZ_cut_300"] = 1.6;
 832      upperLimitsObserved["HTlep_2ltau_offZ_cut_0"] = 25.;
 833      upperLimitsObserved["HTlep_2ltau_offZ_cut_100"] = 14.;
 834      upperLimitsObserved["HTlep_2ltau_offZ_cut_150"] = 6.1;
 835      upperLimitsObserved["HTlep_2ltau_offZ_cut_200"] = 3.3;
 836      upperLimitsObserved["HTlep_2ltau_offZ_cut_300"] = 1.2;
 837      upperLimitsObserved["HTlep_3l_onZ_cut_0"] = 48.;
 838      upperLimitsObserved["HTlep_3l_onZ_cut_100"] = 38.;
 839      upperLimitsObserved["HTlep_3l_onZ_cut_150"] = 14.;
 840      upperLimitsObserved["HTlep_3l_onZ_cut_200"] = 7.2;
 841      upperLimitsObserved["HTlep_3l_onZ_cut_300"] = 4.5;
 842      upperLimitsObserved["HTlep_2ltau_onZ_cut_0"] = 85.;
 843      upperLimitsObserved["HTlep_2ltau_onZ_cut_100"] = 53.;
 844      upperLimitsObserved["HTlep_2ltau_onZ_cut_150"] = 11.0;
 845      upperLimitsObserved["HTlep_2ltau_onZ_cut_200"] = 5.2;
 846      upperLimitsObserved["HTlep_2ltau_onZ_cut_300"] = 3.0;
 847      upperLimitsObserved["METStrong_3l_offZ_cut_0"] = 2.6;
 848      upperLimitsObserved["METStrong_3l_offZ_cut_50"] = 2.1;
 849      upperLimitsObserved["METStrong_3l_offZ_cut_75"] = 2.1;
 850      upperLimitsObserved["METStrong_2ltau_offZ_cut_0"] = 4.2;
 851      upperLimitsObserved["METStrong_2ltau_offZ_cut_50"] = 3.1;
 852      upperLimitsObserved["METStrong_2ltau_offZ_cut_75"] = 2.6;
 853      upperLimitsObserved["METStrong_3l_onZ_cut_20"] = 11.0;
 854      upperLimitsObserved["METStrong_3l_onZ_cut_50"] = 6.4;
 855      upperLimitsObserved["METStrong_3l_onZ_cut_75"] = 5.1;
 856      upperLimitsObserved["METStrong_2ltau_onZ_cut_20"] = 5.9;
 857      upperLimitsObserved["METStrong_2ltau_onZ_cut_50"] = 3.4;
 858      upperLimitsObserved["METStrong_2ltau_onZ_cut_75"] = 1.2;
 859      upperLimitsObserved["METWeak_3l_offZ_cut_0"] = 11.;
 860      upperLimitsObserved["METWeak_3l_offZ_cut_50"] = 5.3;
 861      upperLimitsObserved["METWeak_3l_offZ_cut_75"] = 3.1;
 862      upperLimitsObserved["METWeak_2ltau_offZ_cut_0"] = 23.;
 863      upperLimitsObserved["METWeak_2ltau_offZ_cut_50"] = 4.3;
 864      upperLimitsObserved["METWeak_2ltau_offZ_cut_75"] = 3.1;
 865      upperLimitsObserved["METWeak_3l_onZ_cut_20"] = 41.;
 866      upperLimitsObserved["METWeak_3l_onZ_cut_50"] = 16.;
 867      upperLimitsObserved["METWeak_3l_onZ_cut_75"] = 8.0;
 868      upperLimitsObserved["METWeak_2ltau_onZ_cut_20"] = 80.;
 869      upperLimitsObserved["METWeak_2ltau_onZ_cut_50"] = 4.4;
 870      upperLimitsObserved["METWeak_2ltau_onZ_cut_75"] = 1.8;
 871      upperLimitsObserved["Meff_3l_offZ_cut_0"] = 11.;
 872      upperLimitsObserved["Meff_3l_offZ_cut_150"] = 8.1;
 873      upperLimitsObserved["Meff_3l_offZ_cut_300"] = 3.1;
 874      upperLimitsObserved["Meff_3l_offZ_cut_500"] = 2.1;
 875      upperLimitsObserved["Meff_2ltau_offZ_cut_0"] = 25.;
 876      upperLimitsObserved["Meff_2ltau_offZ_cut_150"] = 12.;
 877      upperLimitsObserved["Meff_2ltau_offZ_cut_300"] = 3.9;
 878      upperLimitsObserved["Meff_2ltau_offZ_cut_500"] = 2.2;
 879      upperLimitsObserved["Meff_3l_onZ_cut_0"] = 48.;
 880      upperLimitsObserved["Meff_3l_onZ_cut_150"] = 37.;
 881      upperLimitsObserved["Meff_3l_onZ_cut_300"] = 11.;
 882      upperLimitsObserved["Meff_3l_onZ_cut_500"] = 4.8;
 883      upperLimitsObserved["Meff_2ltau_onZ_cut_0"] = 85.;
 884      upperLimitsObserved["Meff_2ltau_onZ_cut_150"] = 28.;
 885      upperLimitsObserved["Meff_2ltau_onZ_cut_300"] = 5.9;
 886      upperLimitsObserved["Meff_2ltau_onZ_cut_500"] = 1.9;
 887      upperLimitsObserved["MeffStrong_3l_offZ_cut_0"] = 3.8;
 888      upperLimitsObserved["MeffStrong_3l_offZ_cut_150"] = 3.8;
 889      upperLimitsObserved["MeffStrong_3l_offZ_cut_300"] = 2.8;
 890      upperLimitsObserved["MeffStrong_3l_offZ_cut_500"] = 2.1;
 891      upperLimitsObserved["MeffStrong_2ltau_offZ_cut_0"] = 3.9;
 892      upperLimitsObserved["MeffStrong_2ltau_offZ_cut_150"] = 4.0;
 893      upperLimitsObserved["MeffStrong_2ltau_offZ_cut_300"] = 2.9;
 894      upperLimitsObserved["MeffStrong_2ltau_offZ_cut_500"] = 1.5;
 895      upperLimitsObserved["MeffStrong_3l_onZ_cut_0"] = 10.0;
 896      upperLimitsObserved["MeffStrong_3l_onZ_cut_150"] = 10.0;
 897      upperLimitsObserved["MeffStrong_3l_onZ_cut_300"] = 6.8;
 898      upperLimitsObserved["MeffStrong_3l_onZ_cut_500"] = 3.9;
 899      upperLimitsObserved["MeffStrong_2ltau_onZ_cut_0"] = 1.6;
 900      upperLimitsObserved["MeffStrong_2ltau_onZ_cut_150"] = 1.4;
 901      upperLimitsObserved["MeffStrong_2ltau_onZ_cut_300"] = 1.5;
 902      upperLimitsObserved["MeffStrong_2ltau_onZ_cut_500"] = 0.9;
 903
 904      // Expected upper limits are also given but not used in this analysis
 905      map<string,double> upperLimitsExpected;
 906      upperLimitsExpected["HTlep_3l_offZ_cut_0"] = 11.;
 907      upperLimitsExpected["HTlep_3l_offZ_cut_100"] = 8.5;
 908      upperLimitsExpected["HTlep_3l_offZ_cut_150"] = 4.6;
 909      upperLimitsExpected["HTlep_3l_offZ_cut_200"] = 3.6;
 910      upperLimitsExpected["HTlep_3l_offZ_cut_300"] = 1.9;
 911      upperLimitsExpected["HTlep_2ltau_offZ_cut_0"] = 23.;
 912      upperLimitsExpected["HTlep_2ltau_offZ_cut_100"] = 14.;
 913      upperLimitsExpected["HTlep_2ltau_offZ_cut_150"] = 6.4;
 914      upperLimitsExpected["HTlep_2ltau_offZ_cut_200"] = 3.6;
 915      upperLimitsExpected["HTlep_2ltau_offZ_cut_300"] = 1.5;
 916      upperLimitsExpected["HTlep_3l_onZ_cut_0"] = 33.;
 917      upperLimitsExpected["HTlep_3l_onZ_cut_100"] = 25.;
 918      upperLimitsExpected["HTlep_3l_onZ_cut_150"] = 12.;
 919      upperLimitsExpected["HTlep_3l_onZ_cut_200"] = 6.5;
 920      upperLimitsExpected["HTlep_3l_onZ_cut_300"] = 3.1;
 921      upperLimitsExpected["HTlep_2ltau_onZ_cut_0"] = 94.;
 922      upperLimitsExpected["HTlep_2ltau_onZ_cut_100"] = 61.;
 923      upperLimitsExpected["HTlep_2ltau_onZ_cut_150"] = 9.9;
 924      upperLimitsExpected["HTlep_2ltau_onZ_cut_200"] = 4.5;
 925      upperLimitsExpected["HTlep_2ltau_onZ_cut_300"] = 1.9;
 926      upperLimitsExpected["METStrong_3l_offZ_cut_0"] = 3.1;
 927      upperLimitsExpected["METStrong_3l_offZ_cut_50"] = 2.4;
 928      upperLimitsExpected["METStrong_3l_offZ_cut_75"] = 2.3;
 929      upperLimitsExpected["METStrong_2ltau_offZ_cut_0"] = 4.8;
 930      upperLimitsExpected["METStrong_2ltau_offZ_cut_50"] = 3.3;
 931      upperLimitsExpected["METStrong_2ltau_offZ_cut_75"] = 2.1;
 932      upperLimitsExpected["METStrong_3l_onZ_cut_20"] = 8.7;
 933      upperLimitsExpected["METStrong_3l_onZ_cut_50"] = 4.9;
 934      upperLimitsExpected["METStrong_3l_onZ_cut_75"] = 3.8;
 935      upperLimitsExpected["METStrong_2ltau_onZ_cut_20"] = 7.3;
 936      upperLimitsExpected["METStrong_2ltau_onZ_cut_50"] = 2.8;
 937      upperLimitsExpected["METStrong_2ltau_onZ_cut_75"] = 1.5;
 938      upperLimitsExpected["METWeak_3l_offZ_cut_0"] = 10.;
 939      upperLimitsExpected["METWeak_3l_offZ_cut_50"] = 4.7;
 940      upperLimitsExpected["METWeak_3l_offZ_cut_75"] = 3.0;
 941      upperLimitsExpected["METWeak_2ltau_offZ_cut_0"] = 21.;
 942      upperLimitsExpected["METWeak_2ltau_offZ_cut_50"] = 4.0;
 943      upperLimitsExpected["METWeak_2ltau_offZ_cut_75"] = 2.6;
 944      upperLimitsExpected["METWeak_3l_onZ_cut_20"] = 30.;
 945      upperLimitsExpected["METWeak_3l_onZ_cut_50"] = 10.;
 946      upperLimitsExpected["METWeak_3l_onZ_cut_75"] = 5.4;
 947      upperLimitsExpected["METWeak_2ltau_onZ_cut_20"] = 88.;
 948      upperLimitsExpected["METWeak_2ltau_onZ_cut_50"] = 5.5;
 949      upperLimitsExpected["METWeak_2ltau_onZ_cut_75"] = 2.2;
 950      upperLimitsExpected["Meff_3l_offZ_cut_0"] = 11.;
 951      upperLimitsExpected["Meff_3l_offZ_cut_150"] = 8.8;
 952      upperLimitsExpected["Meff_3l_offZ_cut_300"] = 3.7;
 953      upperLimitsExpected["Meff_3l_offZ_cut_500"] = 2.1;
 954      upperLimitsExpected["Meff_2ltau_offZ_cut_0"] = 23.;
 955      upperLimitsExpected["Meff_2ltau_offZ_cut_150"] = 13.;
 956      upperLimitsExpected["Meff_2ltau_offZ_cut_300"] = 4.9;
 957      upperLimitsExpected["Meff_2ltau_offZ_cut_500"] = 2.4;
 958      upperLimitsExpected["Meff_3l_onZ_cut_0"] = 33.;
 959      upperLimitsExpected["Meff_3l_onZ_cut_150"] = 25.;
 960      upperLimitsExpected["Meff_3l_onZ_cut_300"] = 9.;
 961      upperLimitsExpected["Meff_3l_onZ_cut_500"] = 3.9;
 962      upperLimitsExpected["Meff_2ltau_onZ_cut_0"] = 94.;
 963      upperLimitsExpected["Meff_2ltau_onZ_cut_150"] = 35.;
 964      upperLimitsExpected["Meff_2ltau_onZ_cut_300"] = 6.8;
 965      upperLimitsExpected["Meff_2ltau_onZ_cut_500"] = 2.5;
 966      upperLimitsExpected["MeffStrong_3l_offZ_cut_0"] = 3.9;
 967      upperLimitsExpected["MeffStrong_3l_offZ_cut_150"] = 3.9;
 968      upperLimitsExpected["MeffStrong_3l_offZ_cut_300"] = 3.0;
 969      upperLimitsExpected["MeffStrong_3l_offZ_cut_500"] = 2.0;
 970      upperLimitsExpected["MeffStrong_2ltau_offZ_cut_0"] = 3.8;
 971      upperLimitsExpected["MeffStrong_2ltau_offZ_cut_150"] = 3.9;
 972      upperLimitsExpected["MeffStrong_2ltau_offZ_cut_300"] = 3.1;
 973      upperLimitsExpected["MeffStrong_2ltau_offZ_cut_500"] = 1.6;
 974      upperLimitsExpected["MeffStrong_3l_onZ_cut_0"] = 6.9;
 975      upperLimitsExpected["MeffStrong_3l_onZ_cut_150"] = 7.1;
 976      upperLimitsExpected["MeffStrong_3l_onZ_cut_300"] = 4.9;
 977      upperLimitsExpected["MeffStrong_3l_onZ_cut_500"] = 3.0;
 978      upperLimitsExpected["MeffStrong_2ltau_onZ_cut_0"] = 2.4;
 979      upperLimitsExpected["MeffStrong_2ltau_onZ_cut_150"] = 2.5;
 980      upperLimitsExpected["MeffStrong_2ltau_onZ_cut_300"] = 2.0;
 981      upperLimitsExpected["MeffStrong_2ltau_onZ_cut_500"] = 1.1;
 982
 983      if (observed) return upperLimitsObserved[signal_region];
 984      else          return upperLimitsExpected[signal_region];
 985    }
 986
 987
 988    /// Function checking if there is an OSSF lepton pair or a combination of 3 leptons with an invariant mass close to the Z mass
 989    /// @todo Should the reference Z mass be 91.2?
 990    int isonZ (const Particles& particles)  {
 991      int onZ = 0;
 992      double best_mass_2 = 999.;
 993      double best_mass_3 = 999.;
 994
 995      // Loop over all 2 particle combinations to find invariant mass of OSSF pair closest to Z mass
 996      for ( const Particle& p1 : particles  )  {
 997        for ( const Particle& p2 : particles  )  {
 998          double mass_difference_2_old = fabs(91.0 - best_mass_2);
 999          double mass_difference_2_new = fabs(91.0 - (p1.momentum() + p2.momentum()).mass()/GeV);
1000
1001          // If particle combination is OSSF pair calculate mass difference to Z mass
1002          if ( (p1.pid()*p2.pid() == -121 || p1.pid()*p2.pid() == -169) )  {
1003
1004            // Get invariant mass closest to Z mass
1005            if (mass_difference_2_new < mass_difference_2_old)
1006              best_mass_2 = (p1.momentum() + p2.momentum()).mass()/GeV;
1007
1008            // In case there is an OSSF pair take also 3rd lepton into account (e.g. from FSR and photon to electron conversion)
1009            for ( const Particle & p3 : particles  )  {
1010              double mass_difference_3_old = fabs(91.0 - best_mass_3);
1011              double mass_difference_3_new = fabs(91.0 - (p1.momentum() + p2.momentum() + p3.momentum()).mass()/GeV);
1012              if (mass_difference_3_new < mass_difference_3_old)
1013                best_mass_3 = (p1.momentum() + p2.momentum() + p3.momentum()).mass()/GeV;
1014            }
1015          }
1016        }
1017      }
1018
1019      // Pick the minimum invariant mass of the best OSSF pair combination and the best 3 lepton combination
1020      // If this mass is in a 20 GeV window around the Z mass, the event is classified as onZ
1021      double best_mass = min(best_mass_2, best_mass_3);
1022      if (fabs(91.0 - best_mass) < 20) onZ = 1;
1023      return onZ;
1024    }
1025
1026    /// @}
1027
1028
1029  private:
1030
1031    /// Histograms
1032    /// @{
1033    Histo1DPtr _h_HTlep_all, _h_HTjets_all, _h_MET_all, _h_Meff_all;
1034    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;
1035    Histo1DPtr _h_e_n, _h_mu_n, _h_tau_n;
1036    Histo1DPtr _h_excluded;
1037    /// @}
1038
1039    /// Fiducial efficiencies to model the effects of the ATLAS detector
1040    bool _use_fiducial_lepton_efficiency;
1041
1042    /// List of signal regions and event counts per signal region
1043    vector<string> _signal_regions;
1044
1045    map<string, CounterPtr> _eventCountsPerSR;
1046
1047  };
1048
1049
1050  RIVET_DECLARE_PLUGIN(ATLAS_2012_I1204447);
1051
1052}