Chi2ForSphericalWaveFit.cc

Go to the documentation of this file.

#include "Chi2ForSphericalWaveFit.h"

using namespace std;
using namespace utl;

  
Chi2ForSphericalWaveFit::Chi2ForSphericalWaveFit()
{
  Spherical = true;
}
Chi2ForSphericalWaveFit::~Chi2ForSphericalWaveFit(){}

// Use cartesic coordinate system
void Chi2ForSphericalWaveFit::UseCartesicCoordinates(){Spherical = false;}
// Use spherical coordinate system (default)
void Chi2ForSphericalWaveFit::UseSphericalCoordinates(){Spherical = true;}

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

  double tau0;
  tau0 = TMath::Mean(AntennaTimes.begin(),AntennaTimes.end());

  /* 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
   */
  int n = AntennaTimes.size();
  for (int i=0;i<n;i++)
  {
    AntennaTimes[i] -= tau0;
  }
}

// Objective function
double Chi2ForSphericalWaveFit::DoEval(const double* x) const
{
  double chi = 0.;                      // Chi Square / NDF
  double c = kSpeedOfLight;                     // in AugerUnits
  std::vector<double> ExpectedTimes;     // Expected time for each antenna if the test position is the source position
  

  // Input values
  utl::Vector TestPosition(x[0],x[1],x[2],fgLocalCS,utl::Vector::kSpherical);        // Test Position

  
  // Calculate the expected times
  int n = AntennaPositions.size();
  double exp_time;
  for (int i=0;i<n;i++)
  {
    exp_time = ((AntennaPositions[i]-TestPosition).GetMag())/c;
    ExpectedTimes.push_back(exp_time);
  }

  // get the minimum of ExpectedTimes
  double t0;
  t0 = TMath::Mean(ExpectedTimes.begin(),ExpectedTimes.end());

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

  // Return Chi Square for input Test Position
  return chi;
}