ObjectBasedMET.hh

// -*- C++ -*-
#ifndef RIVET_TOOLS_OBJECTBASEDMET_HH
#define RIVET_TOOLS_OBJECTBASEDMET_HH
 
#include "Rivet/Tools/ExptResolutionFunctions.hh"
 
namespace Rivet {

  class ObjectBasedMET {
  public:

    ObjectBasedMET() { }

    ObjectBasedMET(ResolutionFunctorPtr<Jet> jet_pt_res,
                   ResolutionFunctorPtr<Jet> jet_phi_res,
 
                   ResolutionFunctorPtr<Particle> el_et_res,
                   ResolutionFunctorPtr<Particle> el_phi_res,
 
                   ResolutionFunctorPtr<Particle> mu_pt_res,
                   ResolutionFunctorPtr<Particle> mu_phi_res,
 
                   ResolutionFunctorPtr<Particle> photon_pt_res = nullptr,
                   ResolutionFunctorPtr<Particle> photon_phi_res = nullptr,
 
                   ResolutionFunctorPtr<Jet> hadtau_pt_res = nullptr,
                   ResolutionFunctorPtr<Jet> hadtau_phi_res = nullptr)
        : _jet_pt_res(std::move(jet_pt_res)),
          _jet_phi_res(std::move(jet_phi_res)),
          _el_et_res(std::move(el_et_res)),
          _el_phi_res(std::move(el_phi_res)),
          _mu_pt_res(std::move(mu_pt_res)),
          _mu_phi_res(std::move(mu_phi_res)),
          _photon_pt_res(std::move(photon_pt_res)),
          _photon_phi_res(std::move(photon_phi_res)),
          _hadtau_pt_res(std::move(hadtau_pt_res)),
          _hadtau_phi_res(std::move(hadtau_phi_res)) { }

    double sig(const ThreeMomentum& MET,
               const Jets& js,
               const Particles& es,
               const Particles& mus,
               const Particles& gammas = {},
               const Jets& hadtaus = {}) {
      // Work out soft term contribution by summing out known contributions.
      ThreeMomentum softVec = MET;
 
      // TODO: maybe Rivet should just have a Matrix2 class?
      Matrix<2> cov_sum = Matrix<2>::mkZero();
      for (const Jet& j : js) {
        softVec += j.p3();
        // TODO Working with rivet matrices is not fun
        const double jptres = _jet_pt_res->resolution(j);
        const double jphires = _jet_phi_res->resolution(j);
        const double metangle = directionalDeltaPhi(MET, j.p3());
 
        update_matrix(cov_sum, jptres, jphires, j.pT2(), metangle);
      }
 
      for (const Particle& e : es) {
        softVec += e.p3();
        const double etres = _el_et_res->resolution(e);
        const double phires = _el_phi_res->resolution(e);
        const double metangle = directionalDeltaPhi(MET, e.p3());
        update_matrix(cov_sum, etres, phires, e.pt2(), metangle);
      }
      for (const Particle& mu : mus) {
        softVec += mu.p3();
        const double ptres = _mu_pt_res->resolution(mu);
        const double phires = _mu_phi_res->resolution(mu);
        const double metangle = directionalDeltaPhi(MET, mu.p3());
        update_matrix(cov_sum, ptres, phires, mu.pt2(), metangle);
      }
 
      if (_photon_phi_res != nullptr && _photon_pt_res != nullptr) {
        for (const Particle& gamma : gammas) {
          softVec += gamma.p3();
          const double ptres = _photon_pt_res->resolution(gamma);
          const double phires = _photon_phi_res->resolution(gamma);
          const double metangle = directionalDeltaPhi(MET, gamma.p3());
          update_matrix(cov_sum, ptres, phires, gamma.pt2(), metangle);
        }
      }
      // If user provides photons, but not photon resolution functors,
      // something unintended may be happening.
      else if (gammas.size() > 0) {
        // TODO: do we have a way of triggering warnings only once? This will be hit on every event
        MSG_WARNING(
            "You have passed photons to a ObjectBasedMET calculator without a photon resolution functor."
            "\n\tThis is probably NOT what you intended to do. Photons will NOT count towards MET "
            "significance.");
      }
 
 
      if (_hadtau_phi_res != nullptr && _hadtau_pt_res != nullptr) {
        for (const Jet& hadtau : hadtaus) {
          softVec += hadtau.p3();
          const double ptres = _hadtau_pt_res->resolution(hadtau);
          const double phires = _hadtau_phi_res->resolution(hadtau);
          const double metangle = directionalDeltaPhi(MET, hadtau.p3());
          update_matrix(cov_sum, ptres, phires, hadtau.pt2(), metangle);
        }
      }
      // If user provides taus, but not tau resolution functors,
      // something unintended may be happening.
      else if (gammas.size() > 0) {
        // TODO: do we have a way of triggering warnings only once? This will be hit on every event
        MSG_WARNING("You have passed taus to a ObjectBasedMET calculator without a tau resolution functor."
                    "\n\tThis is probably NOT what you intended to do. Taus will NOT count towards MET "
                    "significance.");
      }
 
 
      // Assume 10 GeV resolution from soft term
      // TODO: Why can't I initialise rivet vectors/matrices with initialiser lists?
      Matrix<2> softUncert = Matrix<2>::mkIdentity() * 100;
      rotateMatrix2(softUncert, directionalDeltaPhi(MET, softVec));
      cov_sum += softUncert;
 
 
      if (cov_sum.get(0, 0) == 0) return 0;
      // Simple analysis includes a check to exclude unphysically? large rho.
      double rho = cov_sum.get(0, 1) / sqrt(cov_sum.get(0, 0) * cov_sum.get(1, 1));
      rho = rho < 0.9 ? rho : 0.;
      // TODO: The SA code uses a variable called Et, but from my attempts to look inside, Et = pt?
      const double significance = MET.pT() / sqrt(cov_sum.get(0, 0) * (1 - rho * rho));
 
      // TODO: it would be nice to be able to also return the actual (2D) met resolution to use for smearing?
      // However, to get a "consistent" smear, we would need both the "true" objects and the "smeared ones"
      // -- which gets pretty messy pretty fast
 
      return significance;
    }
 
  protected:


    void update_matrix(Matrix<2>& mat,
                       const double ptres,
                       const double phires,
                       const double pt2,
                       const double angle) const {
      Matrix<2> particle_uncert = Matrix<2>::mkZero();
      particle_uncert.set(0, 0, ptres * ptres * pt2);
      particle_uncert.set(1, 1, phires * phires * pt2);
      rotateMatrix2(particle_uncert, angle);
      mat += particle_uncert;
    }

    void rotateMatrix2(Matrix<2>& in, const double phi) const {
      const double c = cos(phi);
      const double s = sin(phi);
      const double cc = c * c;
      const double ss = s * s;
      const double cs = c * s;
 
      const double _00 = in.get(0, 0) * cc + in.get(1, 1) * ss - cs * (in.get(1, 0) + in.get(0, 1));
      const double _01 = in.get(0, 1) * cc - in.get(1, 0) * ss + cs * (in.get(0, 0) - in.get(1, 1));
      const double _10 = in.get(1, 0) * cc - in.get(0, 1) * ss + cs * (in.get(0, 0) - in.get(1, 1));
      const double _11 = in.get(0, 0) * ss + in.get(1, 1) * cc - cs * (in.get(1, 0) + in.get(0, 1));
 
      in.set(0, 0, _00);
      in.set(1, 0, _10);
      in.set(0, 1, _01);
      in.set(1, 1, _11);
    }

    double directionalDeltaPhi(const ThreeMomentum& first, const ThreeMomentum& second) const {
      double dPhi = first.phi() - second.phi();
      if (dPhi > PI / 2.) return dPhi - PI;
      if (dPhi < -PI / 2.) return dPhi + PI;
      return dPhi;
    }

 

    ResolutionFunctorPtr<Jet> _jet_pt_res;
    ResolutionFunctorPtr<Jet> _jet_phi_res;
 
    ResolutionFunctorPtr<Particle> _el_et_res;
    ResolutionFunctorPtr<Particle> _el_phi_res;
 
    ResolutionFunctorPtr<Particle> _mu_pt_res;
    ResolutionFunctorPtr<Particle> _mu_phi_res;
 
    ResolutionFunctorPtr<Particle> _photon_pt_res;
    ResolutionFunctorPtr<Particle> _photon_phi_res;
 
    ResolutionFunctorPtr<Jet> _hadtau_pt_res;
    ResolutionFunctorPtr<Jet> _hadtau_phi_res;


    Rivet::Log& getLog() const {
      return Rivet::Log::getLog("Rivet.ObjectBaseMET");
    }
  };

  template <typename jet_pt_res,
            typename jet_phi_res,
            typename el_pt_res,
            typename el_phi_res,
            typename mu_pt_res,
            typename mu_phi_res,
            typename photon_pt_res = void,
            typename photon_phi_res = void,
            typename tau_pt_res = void,
            typename tau_phi_res = void>
  ObjectBasedMET makeObjectBasedMET() {
 
    // Deal with possible default arguments
    ResolutionFunctorPtr<Particle> photon_pt_ptr = nullptr;
    ResolutionFunctorPtr<Particle> photon_phi_ptr = nullptr;
    ResolutionFunctorPtr<Jet> tau_pt_ptr = nullptr;
    ResolutionFunctorPtr<Jet> tau_phi_ptr = nullptr;
 
    if constexpr (!std::is_void_v<photon_pt_res>) {
      photon_pt_ptr = std::make_unique<photon_pt_res>();
      if constexpr (!std::is_void_v<photon_phi_res>) {
        photon_phi_ptr = std::make_unique<photon_phi_res>();
      }
    }
    if constexpr (!std::is_void_v<tau_pt_res>) {
      tau_pt_ptr = std::make_unique<tau_pt_res>();
      if constexpr (!std::is_void_v<tau_phi_res>) {
        tau_phi_ptr = std::make_unique<tau_phi_res>();
      }
    }
 
    return ObjectBasedMET(std::make_unique<jet_pt_res>(), std::make_unique<jet_phi_res>(),
                          std::make_unique<el_pt_res>(), std::make_unique<el_phi_res>(),
                          std::make_unique<mu_pt_res>(), std::make_unique<mu_phi_res>(),
                          std::move(photon_pt_ptr), std::move(photon_phi_ptr), std::move(tau_pt_ptr),
                          std::move(tau_phi_ptr));
  }

  ObjectBasedMET ATLAS_RUN2_DEFAULT_METSIG() {
    return makeObjectBasedMET<ATLAS_RUN2_EMTOPO_PT_RES, ATLAS_RUN2_JET_PHI_RES,
 
                              ATLAS_RUN2_ELECTRON_ET_RES, ATLAS_RUN2_ELECTRON_PHI_RES,
 
                              ATLAS_RUN2_MUON_PT_RES, ATLAS_RUN2_MUON_PHI_RES>();
    // Photons and Taus deliberately excluded because unvalidated.
    // If you need either of them, you can add them by manually
    // constructing the ObjectBasedMET object.
    // Check them carefully (and let us know if they work!)
  }
 
}
 
#endif