RdPreWaveFitter/Chi2ForPlaneWaveFit.cc

Go to the documentation of this file.

#include "Chi2ForPlaneWaveFit.h"

#include <utl/PhysicalConstants.h>

using namespace std;
using namespace utl;


namespace RdPreWaveFitter {

  // Set Antenna Times and Positions
  void
  Chi2ForPlaneWaveFit::Set(const std::vector< utl::Vector >& _AntennaPositions, const std::vector< double >& _AntennaTimes,
                           const std::vector< double >& _AntennaTimesError, const utl::CoordinateSystemPtr& _fgLocalCS)
  {
    AntennaPositions = _AntennaPositions;
    AntennaTimes = _AntennaTimes;
    AntennaTimesError = _AntennaTimesError;
    fgLocalCS = _fgLocalCS;

    /* an absolute event time for both the test and measured times. The absolute event time
       is choosen to be the mean of all measured/test times */
    const double tau0 = TMath::Mean(AntennaTimes.begin(), AntennaTimes.end());
    for (auto& time : AntennaTimes)
    {
      time -= tau0;
    }
  }


  // Objective function
  double
  Chi2ForPlaneWaveFit::DoEval(const double* x) const
  {
    // Unit vector pointing to where the wave is coming from!
    utl::Vector ShowerAxis(1.*utl::m, x[0], x[1], fgLocalCS, utl::Vector::kSpherical);

    // Calculate the expected times
    std::vector<double> ExpectedTimes; // Expected time for each antenna if the test position is the source position
    double exp_time;
    for (auto pos: AntennaPositions)
    {
      exp_time = (-ShowerAxis * pos) / kSpeedOfLight; // calculate (negative) projection
      ExpectedTimes.push_back(exp_time);

    }
    const double t0 = TMath::Mean(ExpectedTimes.begin(), ExpectedTimes.end());

    // Calculate chi square
    double chi = 0.;
    const int n = AntennaTimes.size();
    for (int i = 0; i < n; i++)
    {
      chi += pow(AntennaTimes[i] - (ExpectedTimes[i]-t0), 2) / pow(AntennaTimesError[i], 2);
    }

    // Return Chi Square for input x
    return chi;
  }
}