Probability.cc

Go to the documentation of this file.

#include <utl/Probability.h>

#include <utl/Vector.h>
#include <utl/CoordinateSystemPtr.h>
#include <utl/CoordinateSystem.h>

#include <TRandom3.h>
#include <TVector3.h>
#include <TRotation.h>
#include <TMath.h>

#include <cmath>


namespace utl {

  Probability::Probability() :
    fRand(0)
  { }


  void
  Probability::SetSeed(const unsigned int seed)
  {
    fRand.SetSeed(seed);
  }


  double
  Probability::GetRayleighPDF(const double x, const double sigma)
  {
    if (x < 0)
      return std::numeric_limits<double>::min();
    return x / (sigma * sigma) * exp(- x * x / (2 * sigma * sigma));
  }


  double
  Probability::GetRayleighCDF(const double x, const double sigma)
  {
    if (x > 0)
      return 1 - exp(- x * x / (2 * sigma * sigma));
    return std::numeric_limits<double>::min();
  }


  double
  Probability::GetNormalPDF(const double x, const double mu, const double sigma)
  {
    return 1. / sqrt(2 * M_PI * sigma * sigma) *
      exp(-1. * (x - mu) * (x - mu) / (2 * sigma * sigma));
  }


  long double
  Probability::GetVonMisesPDF(const double x, const double mu, const double kappa)
  {
    if (kappa > 709)
      return GetNormalPDF(x, mu, sqrt(1./kappa));
    return exp(kappa * cos(x - mu)) / (2 * M_PI * TMath::BesselI0(kappa));
  }


  double
  Probability::GetFisherPDF(const double x, const double kappa)
  {
    const double tmp = kappa * sin(x) / (2 * sinh(kappa)) * exp(kappa * cos(x));
    return (tmp == 0) ? std::numeric_limits<double>::min() : tmp;
  }

  double
  Probability::GetFisherCDF(const double x, const double kappa)
  {
    const double tmp = (exp(kappa) - exp(kappa * cos(x))) / (2 * sinh(kappa));
    return (tmp == 0) ? std::numeric_limits<double>::min() : tmp;
  }


  double
  Probability::GetRandomRayleigh(const double sigma)
  {
    return sqrt(-2 * sigma * sigma * log(1- fRand.Uniform()));
  }


  double
  Probability::GetRandomFisher(const double k)
  {
    return acos(1 + 1 / k * log(1 - fRand.Uniform() * (1 - exp(-2 * k))));
  }


  Vector
  Probability::GetRandomFisher(const Vector& meanDirection, const double k)
  {
    const double alpha = GetRandomFisher(k);
    return GetFisher(meanDirection, alpha);
  }


  Vector
  Probability::GetFisher(const Vector& meanDirection, const double alpha)
  {
    const CoordinateSystemPtr cs = CoordinateSystem::GetRootCoordinateSystem();
    const TVector3 md(meanDirection.GetX(cs), meanDirection.GetY(cs), meanDirection.GetZ(cs));
    const TVector3 rotDirection = GetRandomVectorOnSphere();
    TVector3 rotAxis = md.Cross(rotDirection);
    rotAxis = rotAxis.Unit();

    TRotation rot;
    rot.Rotate(alpha, rotAxis);

    //RotationMatrix M(rotAxis, alpha);
    const TVector3 vnew = rot * md;
    return Vector(vnew.X(), vnew.Y(), vnew.Z(), cs);
  }


  TVector3
  Probability::GetRandomVectorOnSphere()
  {
    const double z = fRand.Uniform(-1, 1);
    const double t = fRand.Uniform(-M_PI, M_PI);
    const double r = sqrt(1 - z * z);
    return TVector3(r * cos(t), r * sin(t), z);
  }

}