Rivet analyses

title: ATLAS_2012_I1112263

3 lepton plus missing transverse energy SUSY search

Experiment: ATLAS (LHC)

Inspire ID: 1112263

Status: VALIDATED

Authors: - Peter Richardson

References: - Expt page: ATLAS-SUSY-2011-25 - Expt page: ATLAS-CONF-2012-023 - arXiv: 1204.5638

Beams: p+ p+

Beam energies: (3500.0, 3500.0)GeV

Run details: - BSM signal events at 7000 GeV.

Search for SUSY using events with 3 leptons in association with missing transverse energy in proton-proton collisions at a centre-of-mass energy of 7 TeV. The data sample has a total integrated luminosity of 2.06 fb$^{-1}$. There is no reference data and in addition to the control plots from the paper the number of events in the two signal regions, correctly normalized to an integrated luminosity 2.06 fb$^{-1}$, are calculated.

Source code:ATLAS_2012_I1112263.cc

```c++ // -- C++ --

include "Rivet/Analysis.hh"

include "Rivet/Projections/ChargedFinalState.hh"

include "Rivet/Projections/FastJets.hh"

include "Rivet/Projections/FinalState.hh"

include "Rivet/Projections/IdentifiedFinalState.hh"

include "Rivet/Projections/VetoedFinalState.hh"

include "Rivet/Projections/VisibleFinalState.hh"

include "Rivet/Tools/Random.hh"

include "Rivet/Tools/RivetMT2.hh"

namespace Rivet {

/// @author Peter Richardson class ATLAS_2012_I1112263 : public Analysis { public:

/// Constructor
RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2012_I1112263);


/// @name Analysis methods
/// @{

/// Book histograms and initialise projections before the run
void init() {

  // projection to find the electrons
  IdentifiedFinalState elecs(Cuts::abseta < 2.47 && Cuts::pT > 10 * GeV);
  elecs.acceptIdPair(PID::ELECTRON);
  declare(elecs, "elecs");

  // projection to find the muons
  IdentifiedFinalState muons(Cuts::abseta < 2.4 && Cuts::pT > 10 * GeV);
  muons.acceptIdPair(PID::MUON);
  declare(muons, "muons");

  // for pTmiss
  declare(VisibleFinalState(Cuts::abseta < 4.9), "vfs");

  VetoedFinalState vfs;
  vfs.addVetoPairId(PID::MUON);

  /// Jet finder
  declare(FastJets(vfs, JetAlg::ANTIKT, 0.4), "AntiKtJets04");

  // all tracks (to do deltaR with leptons)
  declare(ChargedFinalState(Cuts::abseta < 3.0), "cfs");

  // Book histograms
  {
    Histo1DPtr tmp;
    _hist_leptonpT_SR1.push_back(book(tmp, "hist_lepton_pT_1_SR1", 11, 0., 220.));
  }
  {
    Histo1DPtr tmp;
    _hist_leptonpT_SR1.push_back(book(tmp, "hist_lepton_pT_2_SR1", 7, 0., 140.));
  }
  {
    Histo1DPtr tmp;
    _hist_leptonpT_SR1.push_back(book(tmp, "hist_lepton_pT_3_SR1", 8, 0., 160.));
  }
  {
    Histo1DPtr tmp;
    _hist_leptonpT_SR2.push_back(book(tmp, "hist_lepton_pT_1_SR2", 11, 0., 220.));
  }
  {
    Histo1DPtr tmp;
    _hist_leptonpT_SR2.push_back(book(tmp, "hist_lepton_pT_2_SR2", 7, 0., 140.));
  }
  {
    Histo1DPtr tmp;
    _hist_leptonpT_SR2.push_back(book(tmp, "hist_lepton_pT_3_SR2", 8, 0., 160.));
  }
  book(_hist_etmiss_SR1_A, "hist_etmiss_SR1_A", 15, 10., 310.);
  book(_hist_etmiss_SR1_B, "hist_etmiss_SR1_B", 9, 10., 190.);
  book(_hist_etmiss_SR2_A, "hist_etmiss_SR2_A", 15, 10., 310.);
  book(_hist_etmiss_SR2_B, "hist_etmiss_SR2_B", 9, 10., 190.);
  book(_hist_mSFOS, "hist_mSFOF", 9, 0., 180.);

  book(_count_SR1, "count_SR1", 1, 0., 1.);
  book(_count_SR2, "count_SR2", 1, 0., 1.);
}


/// Perform the per-event analysis
void analyze(const Event& event) {

  // Get the jet candidates
  Jets cand_jets = apply<FastJets>(event, "AntiKtJets04")
                       .jetsByPt(Cuts::pT > 20 * GeV && Cuts::abseta < 2.8);

  // Candidate muons
  Particles cand_mu;
  Particles chg_tracks = apply<ChargedFinalState>(event, "cfs").particles();
  for (const Particle& mu : apply<IdentifiedFinalState>(event, "muons").particlesByPt()) {
    double pTinCone = -mu.pT();
    for (const Particle& track : chg_tracks) {
      if (deltaR(mu.momentum(), track.momentum()) <= 0.2) pTinCone += track.pT();
    }
    if (pTinCone < 1.8 * GeV) cand_mu.push_back(mu);
  }

  // Candidate electrons
  Particles cand_e;
  for (const Particle& e : apply<IdentifiedFinalState>(event, "elecs").particlesByPt()) {
    double eta = e.eta();
    // Remove electrons with pT<15 in old veto region
    // (NOT EXPLICIT IN THIS PAPER BUT IN SIMILAR 4 LEPTON PAPER and THIS DESCRPITION
    //  IS MUCH WORSE SO ASSUME THIS IS DONE)
    if (fabs(eta) > 1.37 && fabs(eta) < 1.52 && e.perp() < 15. * GeV) continue;
    double pTinCone = -e.perp();
    for (const Particle& track : chg_tracks) {
      if (deltaR(e.momentum(), track.momentum()) <= 0.2) pTinCone += track.pT();
    }
    if (pTinCone / e.perp() < 0.1) {
      cand_e.push_back(e);
    }
  }

  // Resolve jet/lepton ambiguity
  // (NOT EXPLICIT IN THIS PAPER BUT IN SIMILAR 4 LEPTON PAPER and THIS DESCRPITION
  //  IS MUCH WORSE SO ASSUME THIS IS DONE)
  Jets recon_jets;
  for (const Jet& jet : cand_jets) {
    bool away_from_e = true;
    for (const Particle& e : cand_e) {
      if (deltaR(e.momentum(), jet.momentum()) <= 0.2) {
        away_from_e = false;
        break;
      }
    }
    if (away_from_e) recon_jets.push_back(jet);
  }

  // Only keep electrons more than R=0.4 from jets
  Particles recon_e;
  for (const Particle& e : cand_e) {
    bool away = true;
    for (const Jet& jet : recon_jets) {
      if (deltaR(e.momentum(), jet.momentum()) < 0.4) {
        away = false;
        break;
      }
    }
    // ... and 0.1 from any muons
    if (!away) {
      for (const Particle& mu : cand_e) {
        if (deltaR(mu.momentum(), e.momentum()) < 0.1) {
          away = false;
          break;
        }
      }
    }
    if (away) recon_e.push_back(e);
  }
  // Only keep muons more than R=0.4 from jets
  Particles recon_mu;
  for (const Particle& mu : cand_mu) {
    bool away = true;
    for (const Jet& jet : recon_jets) {
      if (deltaR(mu.momentum(), jet.momentum()) < 0.4) {
        away = false;
        break;
      }
    }
    // ... and 0.1 from any electrona
    if (!away) {
      for (const Particle& e : cand_e) {
        if (deltaR(mu.momentum(), e.momentum()) < 0.1) {
          away = false;
          break;
        }
      }
    }
    if (away) recon_mu.push_back(mu);
  }

  // pTmiss
  Particles vfs_particles = apply<VisibleFinalState>(event, "vfs").particles();
  FourMomentum pTmiss;
  for (const Particle& p : vfs_particles) {
    pTmiss -= p.momentum();
  }
  double eTmiss = pTmiss.pT();

  // Now only use recon_jets, recon_mu, recon_e

  // Reject events with wrong number of leptons
  if (recon_mu.size() + recon_e.size() != 3) {
    MSG_DEBUG("To few charged leptons left after selection");
    vetoEvent;
  }

  // ATLAS calo problem
  if (rand01() <= 0.42) {
    for (const Particle& e : recon_e) {
      double eta = e.eta();
      double phi = e.azimuthalAngle(MINUSPI_PLUSPI);
      if (inRange(eta, -0.1, 1.5) && inRange(phi, -0.9, -0.5)) vetoEvent;
    }
    for (const Jet& jet : recon_jets) {
      const double eta = jet.rapidity();
      const double phi = jet.azimuthalAngle(MINUSPI_PLUSPI);
      if (jet.perp() > 40 * GeV && inRange(eta, -0.1, 1.5) && inRange(phi, -0.9, -0.5)) vetoEvent;
    }
  }

  if (!(!recon_e.empty() && recon_e[0].perp() > 25 * GeV)
      && !(!recon_mu.empty() && recon_mu[0].perp() > 20 * GeV)) {
    MSG_DEBUG("Hardest lepton fails trigger");
    vetoEvent;
  }

  // eTmiss cut
  if (eTmiss < 50 * GeV) vetoEvent;

  // Check at least 1 SFOS pair
  double mSFOS = 1e30, mdiff = 1e30 * GeV;
  size_t nSFOS = 0;
  for (size_t ix = 0; ix < recon_e.size(); ++ix) {
    for (size_t iy = ix + 1; iy < recon_e.size(); ++iy) {
      if (recon_e[ix].pid() * recon_e[iy].pid() > 0) continue;
      ++nSFOS;
      double mtest = (recon_e[ix].momentum() + recon_e[iy].momentum()).mass();
      // Veto is mass<20
      if (mtest < 20 * GeV) vetoEvent;
      if (fabs(mtest - 90 * GeV) < mdiff) {
        mSFOS = mtest;
        mdiff = fabs(mtest - 90 * GeV);
      }
    }
  }
  for (size_t ix = 0; ix < recon_mu.size(); ++ix) {
    for (size_t iy = ix + 1; iy < recon_mu.size(); ++iy) {
      if (recon_mu[ix].pid() * recon_mu[iy].pid() > 0) continue;
      ++nSFOS;
      double mtest = (recon_mu[ix].momentum() + recon_mu[iy].momentum()).mass();
      // Veto is mass < 20*GeV
      if (mtest < 20 * GeV) vetoEvent;
      if (fabs(mtest - 90 * GeV) < mdiff) {
        mSFOS = mtest;
        mdiff = fabs(mtest - 90 * GeV);
      }
    }
  }
  // Require at least 1 SFOS pair
  if (nSFOS == 0) vetoEvent;
  // b-jet veto in SR!
  if (mdiff > 10 * GeV) {
    for (const Jet& jet : recon_jets) {
      if (jet.bTagged() && rand01() <= 0.60) vetoEvent;
    }
  }

  // Histogram filling

  // Region SR1, Z depleted
  if (mdiff > 10 * GeV) {
    _count_SR1->fill(0.5);
    _hist_etmiss_SR1_A->fill(eTmiss);
    _hist_etmiss_SR1_B->fill(eTmiss);
    _hist_mSFOS->fill(mSFOS);
  }
  // Region SR2, Z enriched
  else {
    _count_SR2->fill(0.5);
    _hist_etmiss_SR2_A->fill(eTmiss);
    _hist_etmiss_SR2_B->fill(eTmiss);
  }
  // Make the control plots
  // lepton pT
  size_t ie = 0, imu = 0;
  for (size_t ix = 0; ix < 3; ++ix) {
    Histo1DPtr hist = (mdiff > 10 * GeV) ? _hist_leptonpT_SR1[ix] : _hist_leptonpT_SR2[ix];
    double pTe = (ie < recon_e.size()) ? recon_e[ie].perp() : -1 * GeV;
    double pTmu = (imu < recon_mu.size()) ? recon_mu[imu].perp() : -1 * GeV;
    if (pTe > pTmu) {
      hist->fill(pTe);
      ++ie;
    }
    else {
      hist->fill(pTmu);
      ++imu;
    }
  }
}


void finalize() {
  const double norm = crossSection() / femtobarn * 2.06 / sumOfWeights();
  // These are number of events at 2.06fb^-1 per 20 GeV
  for (size_t ix = 0; ix < 3; ++ix) {
    scale(_hist_leptonpT_SR1[ix], norm * 20.);
    scale(_hist_leptonpT_SR2[ix], norm * 20.);
  }
  scale(_hist_etmiss_SR1_A, norm * 20.);
  scale(_hist_etmiss_SR1_B, norm * 20.);
  scale(_hist_etmiss_SR2_A, norm * 20.);
  scale(_hist_etmiss_SR2_B, norm * 20.);
  scale(_hist_mSFOS, norm * 20.);
  // These are number of events at 2.06fb^-1
  scale(_count_SR1, norm);
  scale(_count_SR2, norm);
}

/// @}

private:

/// @name Histograms
/// @{
vector<Histo1DPtr> _hist_leptonpT_SR1, _hist_leptonpT_SR2;
Histo1DPtr _hist_etmiss_SR1_A, _hist_etmiss_SR1_B, _hist_etmiss_SR2_A, _hist_etmiss_SR2_B;
Histo1DPtr _hist_mSFOS, _count_SR1, _count_SR2;
/// @}

};

RIVET_DECLARE_PLUGIN(ATLAS_2012_I1112263);

} ```