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