Rivet analyses

// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FastJets.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/LeptonFinder.hh"
#include "Rivet/Projections/PromptFinalState.hh"
#include "Rivet/Projections/VetoedFinalState.hh"

namespace Rivet {


  /// @brief ZZ production at 13.6 TeV
  class ATLAS_2023_I2723369 : public Analysis {
  public:

    /// Default constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2023_I2723369);

    void init() {

      /// Dressed leptons
      Cut cut_lep = (Cuts::abspid == PID::MUON && Cuts::abseta < 2.5 && Cuts::pT > 5 * GeV)
          || (Cuts::abspid == PID::ELECTRON && Cuts::abseta < 2.47 && Cuts::pT > 7 * GeV);
      PromptFinalState photons(Cuts::abspid == PID::PHOTON);
      PromptFinalState leptons(Cuts::abspid == PID::MUON || Cuts::abspid == PID::ELECTRON);
      LeptonFinder dLeptons(leptons, photons, 0.1, cut_lep);
      declare(dLeptons, "AllLeptons");

      /// Jet inputs
      FinalState fs_jet(Cuts::abseta < 5.0);
      VetoedFinalState jet_input(fs_jet);

      // reject all leptons dressed with only prompt photons from jet input
      FinalState all_leptons(Cuts::abspid == PID::ELECTRON || Cuts::abspid == PID::MUON);
      LeptonFinder reject_leptons(all_leptons, photons, 0.1);
      jet_input.addVetoOnThisFinalState(reject_leptons);

      // reject prompt invisibles, including from tau decays
      VetoedFinalState invis_fs_jet(fs_jet);
      invis_fs_jet.addVetoOnThisFinalState(VisibleFinalState(fs_jet));
      PromptFinalState invis_pfs_jet(invis_fs_jet, TauDecaysAs::PROMPT);
      jet_input.addVetoOnThisFinalState(invis_pfs_jet);

      // declare jets
      FastJets jets(jet_input, JetAlg::ANTIKT, 0.4, JetMuons::NONE, JetInvisibles::DECAY);
      declare(jets, "Jets");

      // Book histograms
      book(_h["m4l"], 1, 1, 1);
      book(_h["pt4l"], 2, 1, 1);
    }

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


      const DressedLeptons& all_leps = apply<LeptonFinder>(e, "AllLeptons").dressedLeptons();
      size_t n_parts = all_leps.size();
      size_t n_OSSF_pairs = 0;

      // Form Z candidate (opposite-sign same-flavour) lepton pairs
      std::vector<Zstate> dileptons;
      for (size_t i = 0; i < n_parts; ++i) {
        for (size_t j = i + 1; j < n_parts; ++j) {
          if (isOSSF(all_leps[i], all_leps[j])) {
            n_OSSF_pairs += 1;
            // Set positive charge first for later ME calculation
            if (all_leps[i].charge() > 0) {
              dileptons.emplace_back(make_pair(all_leps[i], all_leps[j]));
            }
            else {
              dileptons.emplace_back(make_pair(all_leps[j], all_leps[i]));
            }
          }
        }
      }

      // At least two pairs required to select ZZ->llll final state
      if (n_OSSF_pairs < 2) vetoEvent;

      // Form the quadruplet of two lepon pairs passing kinematics cuts
      std::vector<Quadruplet> quadruplets;
      for (size_t i = 0; i < dileptons.size(); ++i) {
        for (size_t j = i + 1; j < dileptons.size(); ++j) {
          // Only use unique leptons
          if (isSame(dileptons[i].first, dileptons[j].first)) continue;
          if (isSame(dileptons[i].first, dileptons[j].second)) continue;
          if (isSame(dileptons[i].second, dileptons[j].first)) continue;
          if (isSame(dileptons[i].second, dileptons[j].second)) continue;

          Particles leptons{dileptons[i].first, dileptons[i].second, dileptons[j].first, dileptons[j].second};
          isortByPt(leptons);

          // Apply kinematic cuts
          if (leptons[0].pT() < 20 * GeV) continue;
          if (leptons[1].pT() < 10 * GeV) continue;

          // Form the quad with pair closest to Z pole first
          if (dileptons[i].Zdist() < dileptons[j].Zdist()) {
            quadruplets.emplace_back(dileptons[i], dileptons[j]);
          }
          else {
            quadruplets.emplace_back(dileptons[j], dileptons[i]);
          }
        }
      }

      // Veto if no quad passes kinematic selection
      size_t n_quads = quadruplets.size();
      if (n_quads == 0) vetoEvent;

      // To resolve ambiguities in lepton pairing order quads by channel priority first, then m12 - mz and m34 - mz
      // The first in every channel is considered nominal
      std::sort(quadruplets.begin(), quadruplets.end(), [](const Quadruplet& q1, const Quadruplet& q2) {
        if (q1.ch_priority == q2.ch_priority) {
          // if rarely, Z1 the same distance from the Z pole, compare Z2
          if (fabs(q1.Z1().Zdist() - q2.Z1().Zdist()) < 1.e-5) {
            return q1.Z2().Zdist() < q2.Z2().Zdist();
          }
          return q1.Z1().Zdist() < q2.Z1().Zdist();
        }
        return q1.ch_priority < q2.ch_priority;
      });


      // Select the best quad
      Particles leptons_sel4l;
      Quadruplet quadSel;
      bool atleastonequadpassed = false;

      for (size_t iquad = 0; iquad < n_quads; ++iquad) {

        // Veto event if nominal quad was not selected in 4 lepton case
        if (n_parts == 4 && iquad > 0) vetoEvent;

        const Quadruplet& quad = quadruplets[iquad];

        // Z invariant mass requirements
        if (!(inRange(quad.Z1().mom().mass(), 66 * GeV, 116 * GeV))) continue;
        if (!(inRange(quad.Z2().mom().mass(), 66 * GeV, 116 * GeV))) continue;

        // Lepton separation and J/Psi veto
        bool b_pass_leptonseparation = true;
        bool b_pass_jpsi = true;
        leptons_sel4l.clear();
        leptons_sel4l.push_back(quad.Z1().first);
        leptons_sel4l.push_back(quad.Z1().second);
        leptons_sel4l.push_back(quad.Z2().first);
        leptons_sel4l.push_back(quad.Z2().second);

        for (size_t i = 0; i < 4; ++i) {
          for (size_t j = i + 1; j < 4; ++j) {
            if (deltaR(leptons_sel4l[i], leptons_sel4l[j]) < 0.1) b_pass_leptonseparation = false;
            if (isOSSF(leptons_sel4l[i], leptons_sel4l[j])
                && (leptons_sel4l[i].mom() + leptons_sel4l[j].mom()).mass() <= 5. * GeV)
              b_pass_jpsi = false;
          }
        }
        if (b_pass_leptonseparation == false || b_pass_jpsi == false) continue;

        isortByPt(leptons_sel4l);
        atleastonequadpassed = true;
        quadSel = quad;
        break;
      }

      if (!atleastonequadpassed) vetoEvent;

      // Veto if quad not in Higgs mass window
      const FourMomentum ZZ = quadSel.mom();
      _h["pt4l"]->fill(ZZ.pT() / GeV);
      _h["m4l"]->fill(ZZ.mass() / GeV);
    }


    void finalize() {
      scale(_h, crossSection() / femtobarn / sumOfWeights());
    }

  private:

    map<string, Histo1DPtr> _h;

    /// Generic Z candidate
    struct Zstate : public ParticlePair {

      Zstate() { }

      Zstate(ParticlePair&& _particlepair)
          : ParticlePair(std::move(_particlepair)) { }

      FourMomentum mom() const {
        return first.momentum() + second.momentum();
      }

      double Zdist() const {
        return fabs(mom().mass() - 91.1876 * GeV);
      }

      int flavour() const {
        return first.abspid();
      }
    };

    /// Generic quadruplet
    struct Quadruplet {

      // find out which type it is: 4mu = 0, 4e = 1, 2mu2e = 2, 2e2mu = 3 (mm, ee, me, em)
      // channel priority is 4m, 2e2m, 2m2e, 4e
      enum class FlavCombi { mm = 0, ee, me, em, undefined };

      Zstate _z1, _z2;

      FlavCombi _type;

      int ch_priority;

      Quadruplet() { }

      Quadruplet(const Zstate& z1, const Zstate& z2)
          : _z1(z1), _z2(z2) {
        if (_z1.flavour() == 13 && _z2.flavour() == 13) {
          _type = FlavCombi::mm;
          ch_priority = 0;
        }
        else if (_z1.flavour() == 11 && _z2.flavour() == 11) {
          _type = FlavCombi::ee;
          ch_priority = 3;
        }
        else if (_z1.flavour() == 13 && _z2.flavour() == 11) {
          _type = FlavCombi::me;
          ch_priority = 2;
        }
        else if (_z1.flavour() == 11 && _z2.flavour() == 13) {
          _type = FlavCombi::em;
          ch_priority = 1;
        }
        else {
          _type = FlavCombi::undefined;
        }
      }

      const Zstate& Z1() const {
        return _z1;
      }

      const Zstate& Z2() const {
        return _z2;
      }

      FourMomentum mom() const {
        return _z1.mom() + _z2.mom();
      }

      const FlavCombi& type() const {
        return _type;
      }
    };
  };

  RIVET_DECLARE_PLUGIN(ATLAS_2023_I2723369);
}