Rivet analyses

// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/DileptonFinder.hh"
#include "Rivet/Projections/FastJets.hh"
#include "Rivet/Projections/VetoedFinalState.hh"

namespace Rivet {


  /// Z + jets in pp at 13 TeV
  class ATLAS_2017_I1514251 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2017_I1514251);

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

      // Get options from the new option system
      // default to combined.
      _mode = 2;
      if (getOption("LMODE") == "EL") _mode = 0;
      if (getOption("LMODE") == "MU") _mode = 1;
      if (getOption("LMODE") == "EMU") _mode = 2;

      Cut cuts = Cuts::pT > 25 * GeV && Cuts::abseta < 2.5;
      DileptonFinder zeefinder(91.2 * GeV, 0.1, cuts && Cuts::abspid == PID::ELECTRON,
                               Cuts::massIn(71 * GeV, 111 * GeV));
      DileptonFinder zmumufinder(91.2 * GeV, 0.1, cuts && Cuts::abspid == PID::MUON,
                                 Cuts::massIn(71 * GeV, 111 * GeV));
      declare(zeefinder, "zeefinder");
      declare(zmumufinder, "zmumufinder");

      // Define veto FS in order to prevent Z-decay products entering the jet algorithm
      VetoedFinalState had_fs;
      had_fs.addVetoOnThisFinalState(zeefinder);
      had_fs.addVetoOnThisFinalState(zmumufinder);
      FastJets jets(had_fs, JetAlg::ANTIKT, 0.4, JetMuons::ALL, JetInvisibles::DECAY);
      declare(jets, "jets");

      // individual channels
      book(_h_Njets_excl, _mode + 1, 1, 1);
      book(_h_Njets, _mode + 4, 1, 1);
      book(_h_Njets_Ratio, _mode + 7, 1, 1);

      book(_h_leading_jet_pT_eq1jet, _mode + 10, 1, 1);
      book(_h_leading_jet_pT, _mode + 13, 1, 1);
      book(_h_leading_jet_pT_2jet, _mode + 16, 1, 1);
      book(_h_leading_jet_pT_3jet, _mode + 19, 1, 1);
      book(_h_leading_jet_pT_4jet, _mode + 22, 1, 1);
      book(_h_leading_jet_rap, _mode + 25, 1, 1);
      book(_h_HT, _mode + 28, 1, 1);
      book(_h_jet_dphi, _mode + 31, 1, 1);
      book(_h_jet_mass, _mode + 34, 1, 1);
    }


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

      const DileptonFinder& zeefinder = apply<DileptonFinder>(event, "zeefinder");
      const DileptonFinder& zmumufinder = apply<DileptonFinder>(event, "zmumufinder");

      const Particles& zees = zeefinder.bosons();
      const Particles& zmumus = zmumufinder.bosons();

      //Veto Z->mumu in electron mode, and vice versa:
      if (_mode == 0 && (zees.size() != 1 || zmumus.size())) vetoEvent;

      if (_mode == 1 && (zees.size() || zmumus.size() != 1)) vetoEvent;

      if (zees.size() + zmumus.size() != 1) {
        // Running in combined mode, we did not find exactly one Z. Not good.
        MSG_DEBUG("Did not find exactly one good Z candidate");
        vetoEvent;
      }

      // Find the (dressed!) leptons
      const Particles& leptons = zees.size() ? zeefinder.constituents() : zmumufinder.constituents();
      if (leptons.size() != 2) vetoEvent;

      Jets jets = apply<JetFinder>(event, "jets").jetsByPt(Cuts::pT > 30 * GeV && Cuts::absrap < 2.5);

      bool veto = false;
      for (const Jet& j : jets) {
        for (const Particle& l : leptons) {
          veto |= deltaR(j, l) < 0.4;
        }
      }
      if (veto) vetoEvent;

      double HT = 0;
      for (const Particle& l : leptons) {
        HT += l.pT();
      }

      const size_t Njets = jets.size();
      _h_Njets_excl->fill(Njets);
      for (size_t i = 0; i <= Njets; ++i) {
        _h_Njets->fill(i);
      }

      if (Njets < 1) vetoEvent;


      for (size_t i = 0; i < Njets; ++i) {
        HT += jets[i].pT();
      }
      const double pT = jets[0].pT();
      const double rap = jets[0].rapidity();

      _h_HT->fill(HT);
      _h_leading_jet_rap->fill(fabs(rap));
      _h_leading_jet_pT->fill(pT);
      if (Njets == 1) _h_leading_jet_pT_eq1jet->fill(pT);
      if (Njets > 1) {
        _h_leading_jet_pT_2jet->fill(pT);
        _h_jet_dphi->fill(deltaPhi(jets[0], jets[1]));
        _h_jet_mass->fill((jets[0].momentum() + jets[1].momentum()).mass());
      }

      if (Njets > 2) _h_leading_jet_pT_3jet->fill(pT);
      if (Njets > 3) _h_leading_jet_pT_4jet->fill(pT);
    }

    void finalize() {
      for (size_t i = 1; i < _h_Njets->numBins(); ++i) {
        double n = _h_Njets->bin(i + 1).sumW();
        double dN = _h_Njets->bin(i + 1).sumW2();
        double d = _h_Njets->bin(i).sumW();
        double dD = _h_Njets->bin(i).sumW2();
        double r = safediv(n, d);
        double e = sqrt(safediv(r * (1 - r), d));
        if (_h_Njets->effNumEntries() != _h_Njets->numEntries()) {
          // use F. James's approximation for weighted events:
          e = sqrt(safediv((1 - 2 * r) * dN + r * r * dD, d * d));
        }
        _h_Njets_Ratio->bin(i).set(r, e);
      }

      // when running in combined mode, need to average to get lepton xsec
      double normfac = crossSectionPerEvent();
      if (_mode == 2) normfac = 0.5 * normfac;

      scale(_h_Njets, normfac);
      scale(_h_Njets_excl, normfac);
      scale(_h_HT, normfac);
      scale(_h_leading_jet_rap, normfac);
      scale(_h_leading_jet_pT, normfac);
      scale(_h_leading_jet_pT_eq1jet, normfac);
      scale(_h_leading_jet_pT_2jet, normfac);
      scale(_h_leading_jet_pT_3jet, normfac);
      scale(_h_leading_jet_pT_4jet, normfac);
      scale(_h_jet_dphi, normfac);
      scale(_h_jet_mass, normfac);
    }


  protected:

    size_t _mode;


  private:

    Estimate1DPtr _h_Njets_Ratio;
    Histo1DPtr _h_Njets;
    Estimate1DPtr _h_Njets_excl_Ratio;
    Histo1DPtr _h_Njets_excl;
    Histo1DPtr _h_HT;
    Histo1DPtr _h_leading_jet_rap;
    Histo1DPtr _h_leading_jet_pT;
    Histo1DPtr _h_leading_jet_pT_eq1jet;
    Histo1DPtr _h_leading_jet_pT_2jet;
    Histo1DPtr _h_leading_jet_pT_3jet;
    Histo1DPtr _h_leading_jet_pT_4jet;
    Histo1DPtr _h_jet_dphi;
    Histo1DPtr _h_jet_mass;
  };


  RIVET_DECLARE_PLUGIN(ATLAS_2017_I1514251);

}