SmearedMET.hh

// -*- C++ -*-
#ifndef RIVET_SmearedMET_HH
#define RIVET_SmearedMET_HH
 
#include "Rivet/Projection.hh"
#include "Rivet/Projections/METFinder.hh"
#include "Rivet/Projections/MissingMomentum.hh"
#include "Rivet/Tools/SmearingFunctions.hh"
#include <functional>
 
namespace Rivet {
 

  class SmearedMET : public METFinder {
  public:


    template <
        typename SMEARPARAMSFN,
        typename std::enable_if_t<is_same_v<invoke_result_t<SMEARPARAMSFN, Vector3, double>, METSmearParams>,
                                  int> = 0>
    SmearedMET(const MissingMomentum& mm, const SMEARPARAMSFN& metSmearParamsFn)
        : _metSmearParamsFn(metSmearParamsFn),
          _metSmearFn(nullptr)
    // _metSmearFn([&](const Vector3& met, double set) -> Vector3 {
    //   const METSmearParams msps = metSmearParamsFn(met, set);
    //   return MET_SMEAR_NORM(msps);
    // })
    {
      setName("SmearedMET");
      declare(mm, "TruthMET");
      _noSmear = getEnvParam<bool>("RIVET_DISABLE_SMEARING", false);
    }

    template <
        typename SMEARPARAMSFN,
        typename std::enable_if_t<is_same_v<invoke_result_t<SMEARPARAMSFN, Vector3, double>, METSmearParams>,
                                  int> = 0>
    SmearedMET(const SMEARPARAMSFN& metSmearParamsFn, const Cut& cut = Cuts::OPEN)
        : SmearedMET(MissingMomentum(cut), metSmearParamsFn) { }

    template <
        typename SMEARFN,
        typename std::enable_if_t<is_same_v<invoke_result_t<SMEARFN, Vector3, double>, Vector3>, int> = 0>
    SmearedMET(const MissingMomentum& mm, const SMEARFN& metSmearFn)
        : _metSmearParamsFn(nullptr), _metSmearFn(metSmearFn) {
      setName("SmearedMET");
      declare(mm, "TruthMET");
      _noSmear = getEnvParam<bool>("RIVET_DISABLE_SMEARING", false);
    }

    template <
        typename SMEARFN,
        typename std::enable_if_t<is_same_v<invoke_result_t<SMEARFN, Vector3, double>, Vector3>, int> = 0>
    SmearedMET(const SMEARFN& metSmearFn, const Cut& cut = Cuts::OPEN)
        : SmearedMET(MissingMomentum(cut), metSmearFn) { }

    template <typename SMEARPARAMSFN,
              typename SMEARFN,
              typename std::enable_if_t<
                  is_same_v<invoke_result_t<SMEARFN, Vector3, double>, Vector3>
                      && is_same_v<invoke_result_t<SMEARPARAMSFN, Vector3, double>, METSmearParams>,
                  int> = 0>
    SmearedMET(const MissingMomentum& mm, const SMEARPARAMSFN& metSmearParamsFn, const SMEARFN& metSmearFn)
        : _metSmearParamsFn(metSmearParamsFn), _metSmearFn(metSmearFn) {
      setName("SmearedMET");
      declare(mm, "TruthMET");
      _noSmear = getEnvParam<bool>("RIVET_DISABLE_SMEARING", false);
    }

    template <typename SMEARPARAMSFN,
              typename SMEARFN,
              typename std::enable_if_t<
                  is_same_v<invoke_result_t<SMEARFN, Vector3, double>, Vector3>
                      && is_same_v<invoke_result_t<SMEARPARAMSFN, Vector3, double>, METSmearParams>,
                  int> = 0>
    SmearedMET(const SMEARPARAMSFN& metSmearParamsFn, const SMEARFN& metSmearFn, const Cut& cut = Cuts::OPEN)
        : SmearedMET(MissingMomentum(cut), metSmearParamsFn, metSmearFn) { }
 

    RIVET_DEFAULT_PROJ_CLONE(SmearedMET);


    using Projection::operator=;
 

    CmpState compare(const Projection& p) const {
      // const SmearedMET& other = dynamic_cast<const SmearedMET&>(p);
      // if (get_address(_metSmearParamsFn) == 0) return cmp((size_t)this, (size_t)&p);
      // if (get_address(_metSmearFn) == 0) return cmp((size_t)this, (size_t)&p);
      // MSG_TRACE("Smear hashes (params) = " << get_address(_metSmearParamsFn) << "," << get_address(other._metSmearParamsFn));
      // MSG_TRACE("Smear hashes (smear) = " << get_address(_metSmearFn) << "," << get_address(other._metSmearFn));
      // return mkPCmp(other, "TruthMET") ||
      //   cmp(get_address(_metSmearParamsFn), get_address(other._metSmearParamsFn));
      //   cmp(get_address(_metSmearFn), get_address(other._metSmearFn));
      if (_noSmear) return CmpState::EQ;
      return CmpState::UNDEF;
    }
 

    void project(const Event& e) {
      const METFinder& mm = apply<MissingMomentum>(e, "TruthMET");
      _spt = mm.scalarPt();
      _set = mm.scalarEt();
      _vpt = mm.vectorPt();
      _vet = mm.vectorEt();
 
      // Short-circuit if smearing is disabled
      if (_noSmear) return;
 
      if (_metSmearFn) {
        _vpt = _metSmearFn(_vpt, _spt); //< custom smearing
        _vet = _metSmearFn(_vet, _set); //< custom smearing
      }
      else if (_metSmearParamsFn) {
        const METSmearParams msps_p = _metSmearParamsFn(
            _vpt,
            _set);                     //< custom smear params: note still calibrated wrt SET, not SPT
        _vpt = MET_SMEAR_NORM(msps_p); //< normal-distribution smearing from custom params
        const METSmearParams msps_e = _metSmearParamsFn(_vet, _set); //< custom smear params
        _vet = MET_SMEAR_NORM(msps_e); //< normal-distribution smearing from custom params
      }
      else {
        throw SmearError("Attempt to smear MPT and MET with neither a smearing function nor a "
                         "smearing-params function set");
      }
    }
 


    const Vector3& vectorSumPt() const {
      return _vpt;
    }

    double scalarSumPt() const {
      return _spt;
    }

    double missingPtResolution() const {
      if (!_metSmearParamsFn)
        throw UserError("Trying to compute MPT significance without a registered significance function");
      METSmearParams msps = _metSmearParamsFn(vectorPt(), scalarPt());
      return msps.pResolution;
    }

    double missingPtSignf() const {
      return missingPt() / missingPtResolution();
    }

 


    const Vector3& vectorSumEt() const {
      return _vet;
    }

    double scalarSumEt() const {
      return _set;
    }

    double missingEtResolution() const {
      if (!_metSmearParamsFn)
        throw UserError("Trying to compute MET significance without a registered significance function");
      METSmearParams msps = _metSmearParamsFn(vectorEt(), scalarEt());
      //return max(msps.eResolution, msps.pResolution);
      return msps.pResolution; //< @todo Check?
    }

    double missingEtSignf() const {
      return missingEt() / missingEtResolution();
    }

 

    void reset() { }
 
 
  protected:
 
    Vector3 _vpt, _vet;
    double _spt, _set;

    METSmearParamsFn _metSmearParamsFn;

    METSmearFn _metSmearFn;

    bool _noSmear{false};
  };
 
 
}
 
#endif