GeometryChi2.cc

Go to the documentation of this file.

#include "GeometryChi2.h"
#include "SphericalShowerFront.h"
#include "GeometryFitter.h"

#include <vector>
#include <iostream>

#include <utl/ErrorLogger.h>
#include <utl/PhysicalConstants.h>
#include <utl/Math.h>

using std::cout;
using std::endl;

namespace MdGeometryFitterAG {

  double
  GeometryChi2::operator()(const std::vector<double>& par)
    const
  {

    DEBUGLOG("Starting GeometryChi2()");

    const double u = par[0];
    const double v = par[1];
    const double ct0 = par[2];
    const double radius = par[3];
    const double xcore = par[4];
    const double ycore = par[5];
    const double zcore = par[6];

//    cout << "u = " << u << ", v = " << v << ", ct0 = " << ct0 << endl;

    const double w2 = 1 - utl::Sqr(u) - utl::Sqr(v);
    if (w2 < 0) {
      ERROR("complex w");
      double value = 1e30;
      return value;
    }
    const double w = sqrt(w2);

    const TVector3 axis(u,v,w);
    const TVector3 core(xcore, ycore, zcore);

    DEBUGLOG("Calculating the chi2");
//    cout << "detector   observedValue    expectedValue   deviance   variance   chi2" << endl;

    const double t0 = ct0 / utl::kSpeedOfLight;
    SphericalShowerFront showerFront(axis, t0, radius, core);

    double chi2 = 0;
    for (std::vector<TimeData>::const_iterator it=timeData.begin(); it!=timeData.end(); ++it) { 
      const double deviance =  it->time - showerFront(it->pos);  // difference between observed and expected time
      const double variance = it->sigma * it->sigma;
      chi2 += deviance*deviance / variance;
    }
    
    DEBUGLOG("Finishing GeometryChi2()");

    return chi2;
  }

}