Chi2ForConicalWaveFit.cc

Go to the documentation of this file.

/*
 * Chi2ForConicalWaveFit.cc
 *
 *  Created on: Oct 30, 2012
 *      Author: leurs
 * modified by Qader Dorosti 
 *          on: Jun 14, 2013 
 */
#include "Chi2ForConicalWaveFit.h"
#include <cmath>

using namespace std;
using namespace utl;


Chi2ForConicalWaveFit::Chi2ForConicalWaveFit()
{
  Cylindrical = true;
}


Chi2ForConicalWaveFit::~Chi2ForConicalWaveFit()
{ }


// Use cartesic coordinate system
void
Chi2ForConicalWaveFit::UseCartesicCoordinates()
{
  Cylindrical = false;
}


// Use cylindrical coordinate system (default)
void
Chi2ForConicalWaveFit::UseCylindricalCoordinates()
{
  Cylindrical = true;
}


// Set Antenna Times and Positions
void
Chi2ForConicalWaveFit::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 chosen to be the mean of all measured/test times
   */
  double tau0 = 0;
  tau0 = TMath::Mean(AntennaTimes.begin(), AntennaTimes.end());
  int n = AntennaTimes.size();
  for (int i = 0; i < n; ++i) {
    AntennaTimes[i] -= tau0;
  }
}


// Objective function
double
Chi2ForConicalWaveFit::DoEval(const double* const x)
  const
{
  double chi = 0; // Chi^2/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 ShowerAxis(1.*utl::m, x[1], x[2], 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) {

    double L = fabs((fabs(ShowerAxis*AntennaPositions[i])*ShowerAxis-AntennaPositions[i]).GetMag());
    exp_time = (-ShowerAxis*AntennaPositions[i] + tan(x[0]) * L) / c;

    ExpectedTimes.push_back(exp_time);

  }

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

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

  //cout << "chi = " << chi << endl;

  // return chi squared
  return chi;
}