MomentumSmearingFunctions.hh

// -*- C++ -*-
#ifndef RIVET_MomentumSmearingFunctions_HH
#define RIVET_MomentumSmearingFunctions_HH
 
#include "Rivet/Math/Vector4.hh"
#include "Rivet/Tools/Random.hh"
 
namespace Rivet {
 



  struct P4SmearParams {
    P4SmearParams(const FourMomentum& p4nom,
                  double pres = 0.0,
                  double eres = 0.0,
                  double phires = 0.0,
                  double etares = 0.0)
        : p4Nominal(p4nom),
          pResolution(pres),
          eResolution(pres),
          phiResolution(phires),
          etaResolution(etares) { }
    FourMomentum p4Nominal;
    double pResolution{0.0}, eResolution{0.0}, phiResolution{0.0}, etaResolution{0.0};
  };

  typedef std::function<FourMomentum(const FourMomentum&)> P4SmearFn;

  typedef std::function<double(const FourMomentum&)> P4EffFn;
 

  inline double P4_EFF_ZERO(const FourMomentum&) {
    return 0;
  }

  inline double P4_EFF_ONE(const FourMomentum&) {
    return 1;
  }

  struct P4_EFF_CONST {
    P4_EFF_CONST(double x)
        : _x(x) { }
    double operator()(const FourMomentum&) const {
      return _x;
    }
    double _x;
  };
 

  inline FourMomentum P4_SMEAR_IDENTITY(const FourMomentum& p) {
    return p;
  }

  inline FourMomentum P4_SMEAR_PERFECT(const FourMomentum& p) {
    return p;
  }

  inline FourMomentum P4_SMEAR_E_GAUSS(const FourMomentum& p, double resolution) {
    const double mass = p.mass2() > 0 ? p.mass() : 0; //< numerical carefulness...
    const double smeared_E = max(randnorm(p.E(), resolution),
                                 mass); //< can't let the energy go below the mass!
    return FourMomentum::mkEtaPhiME(p.eta(), p.phi(), mass, smeared_E);
  }

  inline FourMomentum P4_SMEAR_PT_GAUSS(const FourMomentum& p, double resolution) {
    const double smeared_pt = max(randnorm(p.pT(), resolution), 0.);
    const double mass = p.mass2() > 0 ? p.mass() : 0; //< numerical carefulness...
    return FourMomentum::mkEtaPhiMPt(p.eta(), p.phi(), mass, smeared_pt);
  }

  inline FourMomentum P4_SMEAR_MASS_GAUSS(const FourMomentum& p, double resolution) {
    const double smeared_mass = max(randnorm(p.mass(), resolution), 0.);
    return FourMomentum::mkEtaPhiMPt(p.eta(), p.phi(), smeared_mass, p.pT());
  }

 


  inline double P3_EFF_ZERO(const Vector3&) {
    return 0;
  }

  inline double P3_EFF_ONE(const Vector3&) {
    return 1;
  }

  struct P3_EFF_CONST {
    P3_EFF_CONST(double x)
        : _x(x) { }
    double operator()(const Vector3&) const {
      return _x;
    }
    double _x;
  };
 

  inline Vector3 P3_SMEAR_IDENTITY(const Vector3& p) {
    return p;
  }

  inline Vector3 P3_SMEAR_PERFECT(const Vector3& p) {
    return p;
  }

  inline Vector3 P3_SMEAR_LEN_GAUSS(const Vector3& p, double resolution) {
    const double smeared_mod = max(randnorm(p.mod(), resolution),
                                   0.); //< can't let the energy go below the mass!
    return smeared_mod * p.unit();
  }

 


  struct METSmearParams {
    METSmearParams(const Vector3& vmetnom, double pres = 0.0, double phires = 0.0)
        : vmetNominal(vmetnom), pResolution(pres), phiResolution(phires) { }
    Vector3 vmetNominal;
    double pResolution{0.0}, phiResolution{0.0};
  };

  typedef function<METSmearParams(const Vector3&, double)> METSmearParamsFn;

  typedef function<Vector3(const Vector3&, double)> METSmearFn;
 

  inline Vector3 MET_SMEAR_NORM(const METSmearParams& msps) {
    const Vector3& vmet = msps.vmetNominal;
    // MET magnitude smear
    // const double metsmear = max(randnorm(vmet.mod(), resolution), 0.); //< make peak at 0
    const double metsmear = fabs(randnorm(vmet.mod(), msps.pResolution)); //< "reflect" at 0
    // MET phi-angle smearing matrix
    const Matrix3 phismear = Matrix3::mkZRotation(randnorm(0.0, msps.phiResolution));
    // Apply smearings
    return metsmear * phismear * vmet.unit();
  }
 

  inline METSmearParams MET_SMEARPARAMS_IDENTITY(const Vector3& met, double) {
    return METSmearParams(met, 0.0, 0.0);
  }

  inline Vector3 MET_SMEAR_IDENTITY(const Vector3& met, double) {
    return met;
  }


 
}
 
#endif