Rivet analyses

// -*- 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"

namespace Rivet {


  /// Two-lepton supersymmetry search
  class ATLAS_2011_I894578 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2011_I894578);


    /// @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 > 20 * GeV);
      elecs.acceptIdPair(PID::ELECTRON);
      declare(elecs, "elecs");


      // veto region electrons
      Cut vetocut = Cuts::absetaIn(1.37, 1.52);
      IdentifiedFinalState veto_elecs(vetocut && Cuts::pT > 10 * GeV);
      veto_elecs.acceptIdPair(PID::ELECTRON);
      declare(veto_elecs, "veto_elecs");


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


      // jet finder
      VetoedFinalState vfs;
      vfs.addVetoPairId(PID::MUON);
      declare(FastJets(vfs, JetAlg::ANTIKT, 0.4), "AntiKtJets04");


      // all tracks (to do deltaR with leptons)
      declare(ChargedFinalState(Cuts::abseta < 3 && Cuts::pT > 0.5 * GeV), "cfs");


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


      /// book histograms
      book(_count_OS_e_mu, "count_OS_e+-mu-+", 1, 0., 1.);
      book(_count_OS_e_e, "count_OS_e+e-", 1, 0., 1.);
      book(_count_OS_mu_mu, "count_OS_mu+mu-", 1, 0., 1.);
      book(_count_SS_e_mu, "count_SS_e+-mu+-", 1, 0., 1.);
      book(_count_SS_e_e, "count_SS_e+-e+-", 1, 0., 1.);
      book(_count_SS_mu_mu, "count_SS_mu+-mu+-", 1, 0., 1.);

      book(_hist_eTmiss_OS, "Et_miss_OS", 20, 0., 400.);
      book(_hist_eTmiss_SS, "Et_miss_SS", 20, 0., 400.);
    }


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


      Particles veto_e = apply<IdentifiedFinalState>(event, "veto_elecs").particles();
      if (!veto_e.empty()) {
        MSG_DEBUG("electrons in veto region");
        vetoEvent;
      }

      Jets cand_jets = apply<FastJets>(event, "AntiKtJets04")
                           .jetsByPt(Cuts::pT > 20 * GeV && Cuts::abseta < 2.5);

      Particles cand_e = apply<IdentifiedFinalState>(event, "elecs").particlesByPt();

      // charged particle for isolation
      Particles chg_tracks = apply<ChargedFinalState>(event, "cfs").particles();

      // apply muon isolation
      Particles cand_mu;
      // pTcone around muon track
      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);
      }

      // Discard jets that overlap with electrons
      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);
      }

      // Leptons far from jet
      Particles recon_e;
      for (const Particle& e : cand_e) {
        bool e_near_jet = false;
        for (const Jet& jet : recon_jets) {
          if (deltaR(e.momentum(), jet.momentum()) < 0.4) {
            e_near_jet = true;
            break;
          }
        }
        // Electron isolation criterion
        if (!e_near_jet) {
          double EtinCone = -e.Et();
          for (const Particle& track : chg_tracks) {
            if (deltaR(e.momentum(), track.momentum()) <= 0.2) EtinCone += track.Et();
          }
          if (EtinCone / e.pT() <= 0.15) recon_e.push_back(e);
        }
      }

      Particles recon_mu;
      for (const Particle& mu : cand_mu) {
        bool mu_near_jet = false;
        for (const Jet& jet : recon_jets) {
          if (deltaR(mu.momentum(), jet.momentum()) < 0.4) {
            mu_near_jet = true;
            break;
          }
        }
        if (!mu_near_jet) 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();

      // Exactly two leptons for each event
      if (recon_mu.size() + recon_e.size() != 2) vetoEvent;

      // Lepton pair mass
      FourMomentum p_leptons;
      for (Particle e : recon_e) {
        p_leptons += e.momentum();
      }
      for (Particle mu : recon_mu) {
        p_leptons += mu.momentum();
      }

      if (p_leptons.mass() <= 5.0 * GeV) vetoEvent;

      // ==================== FILL ====================


      // electron, electron
      if (recon_e.size() == 2) {

        // SS ee
        if (recon_e[0].pid() * recon_e[1].pid() > 0) {
          _hist_eTmiss_SS->fill(eTmiss);
          if (eTmiss > 100) {
            MSG_DEBUG("Hits SS e+/-e+/-");
            _count_SS_e_e->fill(0.5);
          }
        }

        // OS ee
        else if (recon_e[0].pid() * recon_e[1].pid() < 0) {
          _hist_eTmiss_OS->fill(eTmiss);
          if (eTmiss > 150) {
            MSG_DEBUG("Hits OS e+e-");
            _count_OS_e_e->fill(0.5);
          }
        }
      }


      // muon, electron
      else if (recon_e.size() == 1) {

        // SS mu_e
        if (recon_e[0].pid() * recon_mu[0].pid() > 0) {
          _hist_eTmiss_SS->fill(eTmiss);
          if (eTmiss > 100) {
            MSG_DEBUG("Hits SS e+/-mu+/-");
            _count_SS_e_mu->fill(0.5);
          }
        }

        // OS mu_e
        else if (recon_e[0].pid() * recon_mu[0].pid() < 0) {
          _hist_eTmiss_OS->fill(eTmiss);
          if (eTmiss > 150) {
            MSG_DEBUG("Hits OS e+mu-");
            _count_OS_e_mu->fill(0.5);
          }
        }
      }


      // muon, muon
      else if (recon_mu.size() == 2) {

        // SS mu_mu
        if (recon_mu[0].pid() * recon_mu[1].pid() > 0) {
          _hist_eTmiss_SS->fill(eTmiss);
          if (eTmiss > 100) {
            MSG_DEBUG("Hits SS mu+/-mu+/-");
            _count_SS_mu_mu->fill(0.5);
          }
        }

        // OS mu_mu
        else if (recon_mu[0].pid() * recon_mu[1].pid() < 0) {
          _hist_eTmiss_OS->fill(eTmiss);
          if (eTmiss > 150) {
            MSG_DEBUG("Hits OS mu+mu-");
            _count_OS_mu_mu->fill(0.5);
          }
        }
      }
    }

    /// @}


    void finalize() {
      double norm = crossSection() / picobarn * 35. / sumOfWeights();
      // event counts
      scale(_count_OS_e_mu, norm);
      scale(_count_OS_e_e, norm);
      scale(_count_OS_mu_mu, norm);
      scale(_count_SS_e_mu, norm);
      scale(_count_SS_e_e, norm);
      scale(_count_SS_mu_mu, norm);
      scale(_hist_eTmiss_OS, 10. * norm);
      scale(_hist_eTmiss_SS, 10. * norm);
    }


  private:

    /// @name Histograms
    /// @{
    Histo1DPtr _count_OS_e_mu;
    Histo1DPtr _count_OS_e_e;
    Histo1DPtr _count_OS_mu_mu;
    Histo1DPtr _count_SS_e_mu;
    Histo1DPtr _count_SS_e_e;
    Histo1DPtr _count_SS_mu_mu;
    Histo1DPtr _hist_eTmiss_OS;
    Histo1DPtr _hist_eTmiss_SS;
    /// @}
  };


  RIVET_DECLARE_ALIASED_PLUGIN(ATLAS_2011_I894578, ATLAS_2011_S9019561);

}