Chi2ForSphericalWaveFitVarC.cc

Go to the documentation of this file.

#include "Chi2ForSphericalWaveFitVarC.h"

using namespace std;
using namespace utl;

  
Chi2ForSphericalWaveFitVarC::Chi2ForSphericalWaveFitVarC()
{
  Spherical = true;
  sigma_gamma = 1./60.;
  doPenalty = true;
  computeSigmaGamma = true;
}
Chi2ForSphericalWaveFitVarC::~Chi2ForSphericalWaveFitVarC(){}

void Chi2ForSphericalWaveFitVarC::SetPenalty(bool dP){doPenalty = dP;}

void Chi2ForSphericalWaveFitVarC::SetSigmaGamma(double _sigma_gamma){
  sigma_gamma = _sigma_gamma;
  computeSigmaGamma = false;
}

// Set Antenna Times and Positions
void Chi2ForSphericalWaveFitVarC::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;
  }

  if(computeSigmaGamma) {
    sigma_gamma = kSpeedOfLight/125.*utl::m*TMath::Sqrt(2)*AntennaTimesError[0];  // (c/125m)*sqrt(2)*sigma
  }
}



// Objective function
double Chi2ForSphericalWaveFitVarC::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
  
  double gamma = 0;
  
  // Input values
  utl::Vector TestPosition(x[0],x[1],x[2],fgLocalCS,utl::Vector::kSpherical);        // Test Position
  
  gamma = x[3];
  if (gamma==0.){gamma=1.e-12;}          // avoid division by zero
  
  // 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*gamma);
    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);
  }
  if(doPenalty){
    chi += pow(1.-gamma,2)/pow(sigma_gamma,2);  // penalty term
  }

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