Rivet/Projections/SmearedMET.hh

0001 // -*- C++ -*-
0002 #ifndef RIVET_SmearedMET_HH
0003 #define RIVET_SmearedMET_HH
0004
0005 #include "Rivet/Projection.hh"
0006 #include "Rivet/Projections/METFinder.hh"
0007 #include "Rivet/Projections/MissingMomentum.hh"
0008 #include "Rivet/Tools/SmearingFunctions.hh"
0009 #include <functional>
0010
0011 namespace Rivet {
0012
0013
0014   /// Wrapper projection for smearing missing (transverse) energy/momentum with detector resolutions
0015   class SmearedMET : public METFinder {
0016   public:
0017
0018     /// @name Constructors etc.
0019     /// @{
0020
0021     /// @brief Constructor from a MissingMomentum projection and a smearing-params function
0022     ///
0023     /// Smearing-params function maps a 3-vector MET and scalar SET to a new METSmearParams struct for sampling
0024     template<typename SMEARPARAMSFN, typename std::enable_if_t<is_same_v<invoke_result_t<SMEARPARAMSFN, Vector3, double>, METSmearParams>, int> = 0>
0025     SmearedMET(const MissingMomentum& mm, const SMEARPARAMSFN& metSmearParamsFn)
0026       : _metSmearParamsFn(metSmearParamsFn), _metSmearFn(nullptr)
0027         // _metSmearFn([&](const Vector3& met, double set) -> Vector3 {
0028         //   const METSmearParams msps = metSmearParamsFn(met, set);
0029         //   return MET_SMEAR_NORM(msps);
0030         // })
0031     {
0032       setName("SmearedMET");
0033       declare(mm, "TruthMET");
0034       _noSmear = getEnvParam<bool>("RIVET_DISABLE_SMEARING", false);
0035     }
0036
0037     /// @brief Constructor from a Cut (on the particles used to determine missing momentum) and a smearing-params function
0038     ///
0039     /// Smearing-params function maps a 3-vector MET and scalar SET to a new METSmearParams struct for sampling
0040     template<typename SMEARPARAMSFN, typename std::enable_if_t<is_same_v<invoke_result_t<SMEARPARAMSFN, Vector3, double>, METSmearParams>, int> = 0>
0041     SmearedMET(const SMEARPARAMSFN& metSmearParamsFn, const Cut& cut=Cuts::OPEN)
0042       : SmearedMET(MissingMomentum(cut), metSmearParamsFn)
0043     {  }
0044
0045     /// @brief Constructor from a MissingMomentum projection and a smearing function
0046     ///
0047     /// Smearing function maps a 3-vector MET and scalar SET to a new MET 3-vector: f(V3, double) -> V3
0048     template<typename SMEARFN, typename std::enable_if_t<is_same_v<invoke_result_t<SMEARFN, Vector3, double>, Vector3>, int> = 0>
0049     SmearedMET(const MissingMomentum& mm, const SMEARFN& metSmearFn)
0050        : _metSmearParamsFn(nullptr), _metSmearFn(metSmearFn)
0051     {
0052       setName("SmearedMET");
0053       declare(mm, "TruthMET");
0054       _noSmear = getEnvParam<bool>("RIVET_DISABLE_SMEARING", false);
0055     }
0056
0057     /// @brief Constructor from a Cut (on the particles used to determine missing momentum) and a smearing function
0058     ///
0059     /// Smearing function maps a 3-vector MET and scalar SET to a new MET 3-vector: f(V3, double) -> V3
0060     template<typename SMEARFN, typename std::enable_if_t<is_same_v<invoke_result_t<SMEARFN, Vector3, double>, Vector3>, int> = 0>
0061     SmearedMET(const SMEARFN& metSmearFn, const Cut& cut=Cuts::OPEN)
0062      : SmearedMET(MissingMomentum(cut), metSmearFn)
0063     {  }
0064
0065     /// @brief Constructor from a MissingMomentum projection and a pair of smearing-params and smearing functions
0066     ///
0067     /// Smearing-params function maps a 3-vector MET and scalar SET to a new METSmearParams struct for sampling
0068     ///
0069     /// Smearing function maps a 3-vector MET and scalar SET to a new MET 3-vector: f(V3, double) -> V3
0070     template <typename SMEARPARAMSFN, typename SMEARFN,
0071               typename std::enable_if_t<is_same_v<invoke_result_t<SMEARFN, Vector3, double>, Vector3> &&
0072                                         is_same_v<invoke_result_t<SMEARPARAMSFN, Vector3, double>, METSmearParams>, int> = 0>
0073     SmearedMET(const MissingMomentum& mm, const SMEARPARAMSFN& metSmearParamsFn, const SMEARFN& metSmearFn)
0074       : _metSmearParamsFn(metSmearParamsFn), _metSmearFn(metSmearFn)
0075     {
0076       setName("SmearedMET");
0077       declare(mm, "TruthMET");
0078       _noSmear = getEnvParam<bool>("RIVET_DISABLE_SMEARING", false);
0079     }
0080
0081     /// @brief Constructor from a Cut (on the particles used to determine missing momentum) and a pair of smearing (params) functions
0082     ///
0083     /// Smearing-params function maps a 3-vector MET and scalar SET to a new METSmearParams struct for sampling
0084     ///
0085     /// Smearing function maps a 3-vector MET and scalar SET to a new MET 3-vector: f(V3, double) -> V3
0086     template <typename SMEARPARAMSFN, typename SMEARFN,
0087               typename std::enable_if_t<is_same_v<invoke_result_t<SMEARFN, Vector3, double>, Vector3> &&
0088                                         is_same_v<invoke_result_t<SMEARPARAMSFN, Vector3, double>, METSmearParams>, int> = 0>
0089     SmearedMET(const SMEARPARAMSFN& metSmearParamsFn, const SMEARFN& metSmearFn, const Cut& cut=Cuts::OPEN)
0090       : SmearedMET(MissingMomentum(cut), metSmearParamsFn, metSmearFn)
0091     {  }
0092
0093
0094     /// Clone on the heap.
0095     RIVET_DEFAULT_PROJ_CLONE(SmearedMET);
0096
0097     /// @}
0098
0099     /// Import to avoid warnings about overload-hiding
0100     using Projection::operator =;
0101
0102
0103     /// Compare to another SmearedMET
0104     CmpState compare(const Projection& p) const {
0105       // const SmearedMET& other = dynamic_cast<const SmearedMET&>(p);
0106       // if (get_address(_metSmearParamsFn) == 0) return cmp((size_t)this, (size_t)&p);
0107       // if (get_address(_metSmearFn) == 0) return cmp((size_t)this, (size_t)&p);
0108       // MSG_TRACE("Smear hashes (params) = " << get_address(_metSmearParamsFn) << "," << get_address(other._metSmearParamsFn));
0109       // MSG_TRACE("Smear hashes (smear) = " << get_address(_metSmearFn) << "," << get_address(other._metSmearFn));
0110       // return mkPCmp(other, "TruthMET") ||
0111       //   cmp(get_address(_metSmearParamsFn), get_address(other._metSmearParamsFn));
0112       //   cmp(get_address(_metSmearFn), get_address(other._metSmearFn));
0113       if (_noSmear) return CmpState::EQ;
0114       return CmpState::UNDEF;
0115     }
0116
0117
0118     /// Perform the MET finding & smearing calculation
0119     void project(const Event& e) {
0120       const METFinder& mm = apply<MissingMomentum>(e, "TruthMET");
0121       _spt = mm.scalarPt();
0122       _set = mm.scalarEt();
0123       _vpt = mm.vectorPt();
0124       _vet = mm.vectorEt();
0125
0126       // Short-circuit if smearing is disabled
0127       if (_noSmear) return;
0128
0129       if (_metSmearFn) {
0130         _vpt = _metSmearFn(_vpt, _spt); //< custom smearing
0131         _vet = _metSmearFn(_vet, _set); //< custom smearing
0132       } else if (_metSmearParamsFn) {
0133         const METSmearParams msps_p = _metSmearParamsFn(_vpt, _set); //< custom smear params: note still calibrated wrt SET, not SPT
0134         _vpt = MET_SMEAR_NORM(msps_p); //< normal-distribution smearing from custom params
0135         const METSmearParams msps_e = _metSmearParamsFn(_vet, _set); //< custom smear params
0136         _vet = MET_SMEAR_NORM(msps_e); //< normal-distribution smearing from custom params
0137       } else {
0138         throw SmearError("Attempt to smear MPT and MET with neither a smearing function nor a smearing-params function set");
0139       }
0140     }
0141
0142
0143     /// @name Transverse-momentum functions
0144     ///
0145     /// @note This may be what you want, even if the paper calls it "missing Et"!
0146     ///
0147     /// NOTE: all other MPT functions are implemented in METFinder via the first two here.
0148     ///
0149     /// @{
0150
0151     /// The vector-summed visible transverse momentum in the event, as a 3-vector with z=0
0152     const Vector3& vectorSumPt() const { return _vpt; }
0153
0154     /// The scalar-summed visible transverse momentum in the event
0155     double scalarSumPt() const { return _spt; }
0156
0157     /// Obtain an approximation to the MPT resolution for this event
0158     double missingPtResolution() const {
0159       if (!_metSmearParamsFn)
0160         throw UserError("Trying to compute MPT significance without a registered significance function");
0161       METSmearParams msps = _metSmearParamsFn(vectorPt(), scalarPt());
0162       return msps.pResolution;
0163     }
0164
0165     /// Obtain an approximation to the MPT significance (value/resolution) for this event
0166     double missingPtSignf() const {
0167       return missingPt() / missingPtResolution();
0168     }
0169
0170     /// @}
0171
0172
0173     /// @name Transverse-energy functions
0174     ///
0175     /// @warning Despite the common names "MET" and "SET", what's often meant is the pT functions above!
0176     ///
0177     /// NOTE: all other MET functions are implemented in METFinder via the first two here.
0178     ///
0179     /// @{
0180
0181     /// The vector-summed visible transverse energy in the event, as a 3-vector with z=0
0182     ///
0183     /// @note Reverse this vector with operator- to get the missing ET vector.
0184     const Vector3& vectorSumEt() const { return _vet; }
0185
0186     /// The scalar-summed visible transverse energy in the event
0187     double scalarSumEt() const { return _set; }
0188
0189     /// Obtain an approximation to the MET resolution for this event
0190     double missingEtResolution() const {
0191       if (!_metSmearParamsFn)
0192         throw UserError("Trying to compute MET significance without a registered significance function");
0193       METSmearParams msps = _metSmearParamsFn(vectorEt(), scalarEt());
0194       //return max(msps.eResolution, msps.pResolution);
0195       return msps.pResolution; //< @todo Check?
0196     }
0197
0198     /// Obtain an approximation to the MET significance (value/resolution) for this event
0199     double missingEtSignf() const {
0200       return missingEt() / missingEtResolution();
0201     }
0202
0203     /// @}
0204
0205
0206     /// Reset the projection. Smearing functions will be unchanged.
0207     void reset() {  }
0208
0209
0210   protected:
0211
0212     Vector3 _vpt, _vet;
0213     double _spt, _set;
0214
0215     /// Stored smearing-params function
0216     METSmearParamsFn _metSmearParamsFn;
0217
0218     /// Stored smearing function
0219     METSmearFn _metSmearFn;
0220
0221     /// Flag to disable smearing
0222     ///
0223     /// @todo Could/should be handled statically and centrally
0224     bool _noSmear{false};
0225
0226   };
0227
0228
0229 }
0230
0231 #endif