testProbability.cc

Go to the documentation of this file.

#include <iostream>

#include <tst/Verify.h>

#include <utl/CoordinateSystem.h>
#include <utl/AugerCoordinateSystem.h>
#include <utl/Probability.h>
#include <utl/CoordinateSystem.h>
#include <utl/CoordinateSystemPtr.h>
#include <utl/Vector.h>
#include <utl/AugerUnits.h>

#include <fwk/CentralConfig.h>
#include <fwk/CoordinateSystemRegistry.h>

#include <cppunit/extensions/HelperMacros.h>

#include <TH1D.h>
#include <TGraph.h>
#include <TF1.h>
#include <TCanvas.h>
#include <TStyle.h>
// #include <Python/geometry.i>

using namespace std;
using namespace utl;
using namespace tst;


class TestProbability : public CppUnit::TestFixture {

  CPPUNIT_TEST_SUITE(TestProbability);
  CPPUNIT_TEST(TestFisher);
  CPPUNIT_TEST(TestVonMises);
  //CPPUNIT_TEST(PlotFisher1);
  //CPPUNIT_TEST(PlotFisher2);

  CPPUNIT_TEST_SUITE_END();

private:
  CoordinateSystemPtr fCS;

public:
  void
  setUp()
  {
    fCS = utl::CoordinateTransformer::GetRootCoordinateSystem();
  }

  void tearDown() { }

  void
  TestVonMises()
  {
    const double kappa = 708.99;
    for (int j = 0; j < 100; ++j) {
      const double x = (0.5*j - 25) * utl::deg;
      const double prob = Probability::GetInstance().GetVonMisesPDF(x, 0, kappa);
      const double probGaus = Probability::GetInstance().GetNormalPDF(x, 0, sqrt(1./kappa));
      CPPUNIT_ASSERT(Verify<CloseTo>(prob, probGaus, 2e-3));
    }
  }

  void
  TestFisher()
  {
    const int nTheta = 100;
    const int nPhi = 100;
    const int nAlpha = 100;
    for (int iTheta = 0; iTheta < nTheta; ++iTheta) {
      const double theta = 0.5 * iTheta / nTheta * M_PI;
      for (int iPhi = 0; iPhi < nPhi; ++iPhi) {
        const double phi = 2. * iPhi / nPhi * M_PI;
        const Vector showeraxis(1, theta, phi, fCS, Vector::kSpherical);
        for (int iAlpha = 0; iAlpha < nAlpha; ++iAlpha) {
            double alpha = 5. * iAlpha/nAlpha * utl::deg;
            const Vector smearedVector = Probability::GetInstance().GetFisher(showeraxis, alpha);
            const double angle = Angle(showeraxis, smearedVector);
            CPPUNIT_ASSERT(Verify<CloseTo>(alpha, angle, 1e-8));
        }
      }
    }
  }

  void
  PlotFisher1()
  {
    const int n = 10000000;
    TH1D h = TH1D("h", "", 90, 0, 0.5);
    TGraph g;
    for (int i = 0; i < n; ++i) {
      h.Fill(Probability::GetInstance().GetRandomFisher(250));
      const double x = 1. * i/n * M_PI;
      g.SetPoint(i, x, Probability::GetFisherPDF(x, 250));
    }
    h.Scale(1/h.Integral("width"));
    // check compatibility
    for (int i = 1; i <= 90; ++i) {
      const double x = h.GetBinCenter(i);
      CPPUNIT_ASSERT(Verify<CloseTo>(h.GetBinContent(i), Probability::GetFisherPDF(x, 250), 1e-2));
    }

    gStyle->SetOptStat("merou");
    TCanvas* const c = new TCanvas();
    c->cd();
    h.GetXaxis()->SetTitle("angle [deg]");
    h.Draw("");
    g.Draw("Psame");
    c->SaveAs("hFisher.png");
    delete c;
  }

  void
  PlotFisher2()
  {
    const int n = 10000000;
    TH1D h = TH1D("h", "", 90, 0, 0.5);
    TGraph g;
    const double sigma = 1 * utl::deg;
    for (int i = 0; i < n; ++i) {
      h.Fill(Probability::GetInstance().GetRandomFisher(1/(sigma*sigma)));
      const double x = 1. * i/n * M_PI;
      g.SetPoint(i, x, Probability::GetRayleighPDF(x, sigma));
    }
    h.Scale(1/h.Integral("width"));

    gStyle->SetOptStat("merou");
    TCanvas* const c = new TCanvas();
    c->cd();
    h.GetXaxis()->SetTitle("angle [deg]");
    h.Draw("");
    g.Draw("Psame");
    c->SaveAs("hFisher_Rayleigh.png");
    delete c;
  }

};


CPPUNIT_TEST_SUITE_REGISTRATION(TestProbability);


// Configure (x)emacs for this file ...
// Local Variables:
// mode: c++
// End: