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"
#include "Rivet/Tools/Random.hh"

namespace Rivet {


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

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2012_I1095236);


    /// @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");

      // Projection to find the muons
      IdentifiedFinalState muons(Cuts::abseta < 2.4 && Cuts::pT > 10 * 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.0), "cfs");

      // Used for pTmiss
      declare(VisibleFinalState((Cuts::etaIn(-4.9, 4.9))), "vfs");

      // Book histograms
      book(_count_SR0_A1, "count_SR0_A1", 1, 0., 1.);
      book(_count_SR0_B1, "count_SR0_B1", 1, 0., 1.);
      book(_count_SR0_C1, "count_SR0_C1", 1, 0., 1.);
      book(_count_SR0_A2, "count_SR0_A2", 1, 0., 1.);
      book(_count_SR0_B2, "count_SR0_B2", 1, 0., 1.);
      book(_count_SR0_C2, "count_SR0_C2", 1, 0., 1.);
      book(_count_SR1_D, "count_SR1_D", 1, 0., 1.);
      book(_count_SR1_E, "count_SR1_E", 1, 0., 1.);

      book(_hist_meff_SR0_A1, "hist_m_eff_SR0_A1", 14, 400., 1800.);
      book(_hist_meff_SR0_A2, "hist_m_eff_SR0_A2", 14, 400., 1800.);
      book(_hist_meff_SR1_D_e, "hist_meff_SR1_D_e", 16, 600., 2200.);
      book(_hist_meff_SR1_D_mu, "hist_meff_SR1_D_mu", 16, 600., 2200.);

      book(_hist_met_SR0_A1, "hist_met_SR0_A1", 14, 0., 700.);
      book(_hist_met_SR0_A2, "hist_met_SR0_A2", 14, 0., 700.);
      book(_hist_met_SR0_D_e, "hist_met_SR1_D_e", 15, 0., 600.);
      book(_hist_met_SR0_D_mu, "hist_met_SR1_D_mu", 15, 0., 600.);
    }


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

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

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

      const Particles cand_mu = apply<IdentifiedFinalState>(event, "muons").particlesByPt();
      // Resolve jet-lepton overlap for jets with |eta| < 2.8
      Jets recon_jets;
      for (const Jet& jet : cand_jets) {
        if (fabs(jet.eta()) >= 2.8) continue;
        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);
      }

      // get the loose leptons used to define the 0 lepton channel
      Particles loose_e, loose_mu;
      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;
          }
        }
        if (away) loose_e.push_back(e);
      }
      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;
          }
        }
        if (away) loose_mu.push_back(mu);
      }
      // tight leptons for the 1-lepton channel
      Particles tight_mu;
      Particles chg_tracks = apply<ChargedFinalState>(event, "cfs").particles();
      for (const Particle& mu : loose_mu) {
        if (mu.perp() < 20.) continue;
        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) tight_mu.push_back(mu);
      }
      Particles tight_e;
      for (const Particle& e : loose_e) {
        if (e.perp() < 25.) 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) {
          tight_e.push_back(e);
        }
      }

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

      // get the number of b-tagged jets
      unsigned int ntagged = 0;
      for (const Jet& jet : recon_jets) {
        if (jet.perp() > 50. && abs(jet.eta()) < 2.5 && jet.bTagged() && rand01() <= 0.60) ++ntagged;
      }

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

      // at least 1 b tag
      if (ntagged == 0) vetoEvent;

      // minumum Et miss
      if (eTmiss < 80.) vetoEvent;

      // at least 3 jets pT > 50
      if (recon_jets.size() < 3 || recon_jets[2].perp() < 50.) vetoEvent;

      // m_eff
      double m_eff = eTmiss;
      for (unsigned int ix = 0; ix < 3; ++ix) m_eff += recon_jets[ix].perp();

      // delta Phi
      double min_dPhi = 999.999;
      double pTmiss_phi = pTmiss.phi();
      for (unsigned int ix = 0; ix < 3; ++ix) {
        min_dPhi = min(min_dPhi, deltaPhi(pTmiss_phi, recon_jets[ix].phi()));
      }

      // 0-lepton channels
      if (loose_e.empty() && loose_mu.empty() && recon_jets[0].perp() > 130. && eTmiss > 130.
          && eTmiss / m_eff > 0.25 && min_dPhi > 0.4) {
        // jet charge cut
        bool jetCharge = true;
        for (unsigned int ix = 0; ix < 3; ++ix) {
          if (fabs(recon_jets[ix].eta()) > 2.) continue;
          double trackpT = 0;
          for (const Particle& p : recon_jets[ix].particles()) {
            if (PID::charge3(p.pid()) == 0) continue;
            trackpT += p.perp();
          }
          if (trackpT / recon_jets[ix].perp() < 0.05) jetCharge = false;
        }
        if (jetCharge) {
          // SR0-A region
          if (m_eff > 500.) {
            _count_SR0_A1->fill(0.5);
            _hist_meff_SR0_A1->fill(m_eff);
            _hist_met_SR0_A1->fill(eTmiss);
            if (ntagged >= 2) {
              _count_SR0_A2->fill(0.5);
              _hist_meff_SR0_A2->fill(m_eff);
              _hist_met_SR0_A2->fill(eTmiss);
            }
          }
          // SR0-B
          if (m_eff > 700.) {
            _count_SR0_B1->fill(0.5);
            if (ntagged >= 2) _count_SR0_B2->fill(0.5);
          }
          // SR0-C
          if (m_eff > 900.) {
            _count_SR0_C1->fill(0.5);
            if (ntagged >= 2) _count_SR0_C2->fill(0.5);
          }
        }
      }

      // 1-lepton channels
      if (tight_e.size() + tight_mu.size() == 1 && recon_jets.size() >= 4 && recon_jets[3].perp() > 50.
          && recon_jets[0].perp() > 60.) {
        Particle lepton = tight_e.empty() ? tight_mu[0] : tight_e[0];
        m_eff += lepton.perp() + recon_jets[3].perp();
        // transverse mass cut
        double mT = 2. * (lepton.perp() * eTmiss - lepton.px() * pTmiss.px() - lepton.py() * pTmiss.py());
        mT = sqrt(mT);
        if (mT > 100. && m_eff > 700.) {
          // D region
          _count_SR1_D->fill(0.5);
          if (lepton.abspid() == PID::ELECTRON) {
            _hist_meff_SR1_D_e->fill(m_eff);
            _hist_met_SR0_D_e->fill(eTmiss);
          }
          else {
            _hist_meff_SR1_D_mu->fill(m_eff);
            _hist_met_SR0_D_mu->fill(eTmiss);
          }
          // E region
          if (eTmiss > 200.) {
            _count_SR1_E->fill(0.5);
          }
        }
      }
    }


    void finalize() {

      double norm = crossSection() / femtobarn * 2.05 / sumOfWeights();
      // these are number of events at 2.05fb^-1 per 100 GeV
      scale(_hist_meff_SR0_A1, 100. * norm);
      scale(_hist_meff_SR0_A2, 100. * norm);
      scale(_hist_meff_SR1_D_e, 100. * norm);
      scale(_hist_meff_SR1_D_mu, 100. * norm);
      // these are number of events at 2.05fb^-1 per 50 GeV
      scale(_hist_met_SR0_A1, 50. * norm);
      scale(_hist_met_SR0_A2, 40. * norm);
      // these are number of events at 2.05fb^-1 per 40 GeV
      scale(_hist_met_SR0_D_e, 40. * norm);
      scale(_hist_met_SR0_D_mu, 40. * norm);
      // these are number of events at 2.05fb^-1
      scale(_count_SR0_A1, norm);
      scale(_count_SR0_B1, norm);
      scale(_count_SR0_C1, norm);
      scale(_count_SR0_A2, norm);
      scale(_count_SR0_B2, norm);
      scale(_count_SR0_C2, norm);
      scale(_count_SR1_D, norm);
      scale(_count_SR1_E, norm);
    }

    /// @}


  private:

    Histo1DPtr _count_SR0_A1;
    Histo1DPtr _count_SR0_B1;
    Histo1DPtr _count_SR0_C1;
    Histo1DPtr _count_SR0_A2;
    Histo1DPtr _count_SR0_B2;
    Histo1DPtr _count_SR0_C2;
    Histo1DPtr _count_SR1_D;
    Histo1DPtr _count_SR1_E;

    Histo1DPtr _hist_meff_SR0_A1;
    Histo1DPtr _hist_meff_SR0_A2;
    Histo1DPtr _hist_meff_SR1_D_e;
    Histo1DPtr _hist_meff_SR1_D_mu;
    Histo1DPtr _hist_met_SR0_A1;
    Histo1DPtr _hist_met_SR0_A2;
    Histo1DPtr _hist_met_SR0_D_e;
    Histo1DPtr _hist_met_SR0_D_mu;
  };


  RIVET_DECLARE_PLUGIN(ATLAS_2012_I1095236);

}