Deprecated/RdWaveFit/Chi2ForPlaneWaveFit.cc

Go to the documentation of this file.

#include "Chi2ForPlaneWaveFit.h"


using namespace std;
using namespace utl;

Chi2ForPlaneWaveFit::Chi2ForPlaneWaveFit(){}
Chi2ForPlaneWaveFit::~Chi2ForPlaneWaveFit(){}

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

// Objective function
double Chi2ForPlaneWaveFit::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
//   utl::Vector::CoordinateType kSpherical;
  
  utl::Vector ShowerAxis(1.*utl::m,x[0],x[1], fgLocalCS, utl::Vector::kSpherical);      // Unit vector pointing to where the wave is coming from!
//   ShowerAxis.SetMagThetaPhi(1.,x[0],x[1]); 
  
  // Calculate the expected times
  int n = AntennaPositions.size();
  double exp_time;

  for (int i=0;i<n;i++)
  {
    exp_time = (-(ShowerAxis*AntennaPositions[i]))/c; // calculate (negative) projection
    ExpectedTimes.push_back(exp_time);

  }
  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 x
  return chi;

}