GeometryFitter.cc

Go to the documentation of this file.

#include "GeometryFitter.h"
#include "GeometryChi2.h"
#include "SphericalShowerFront.h"
#include "PlaneShowerFront.h"
 
#include <utl/ErrorLogger.h>
#include <utl/PhysicalConstants.h>
#include <utl/Math.h>

#include <TMatrixD.h>

#include <Minuit2/MnUserCovariance.h>
#include <Minuit2/MnMigrad.h>
#include <Minuit2/FunctionMinimum.h>
#include <Minuit2/MnUserParameters.h>
#include <Minuit2/MnPrint.h>


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

namespace MdGeometryFitterAG {

  GeometryFitter::GeometryFitter() 
  { 

    // Here goes the definition of the fit parameters
    // They are always referenced by the names defined here 
    parSeed.Add("u", 0, .01);
    parSeed.Add("v", 0, .01);
    parSeed.Add("ct0", 0, 10);
    parSeed.Add("radius", 7500, 500);
    parSeed.Add("xcore", 0, 100);
    parSeed.Add("ycore", 0, 100);
    parSeed.Add("zcore", 0, 0);

    nParameters = 7; // must match the number of parameters defined above

    // Only the axis and core time are free by default
    parSeed.Fix("radius");
    parSeed.Fix("xcore");
    parSeed.Fix("ycore");
    parSeed.Fix("zcore");
    
  }

  void 
  GeometryFitter::AddDetector(int id, TVector3 pos, double time, double sigma) 
  {
      TimeData detector = { id, pos, time, sigma };
      timeData.push_back( detector );
  } 

  void 
  GeometryFitter::SetCore(const TVector3 &core) 
  {
    parSeed.SetValue("xcore", core.x());
    parSeed.SetValue("ycore", core.y());
    parSeed.SetValue("zcore", core.z());
  }

  void 
  GeometryFitter::SetCoreFree() 
  {
    parSeed.Release("xcore"); 
    parSeed.Release("ycore");   // zcore is always fixed
  }

  void 
  GeometryFitter::SetCoreFix() 
  {
    parSeed.Fix("xcore"); 
    parSeed.Fix("ycore");
    parSeed.Fix("zcore");
  }

  int
  GeometryFitter::FitPlane(double &u, double &v, double &ct0)
    const
  {

    double s1 = 0;
    double sx = 0;
    double sy = 0;
    double st = 0;
    double sxx = 0;
    double sxy = 0;
    double syy = 0;
    double sxt = 0;
    double syt = 0;
    //double stt = 0;

    for (std::vector<TimeData>::const_iterator it=timeData.begin(); it!=timeData.end(); ++it) {
      const double sigma = it->sigma;
      const TVector3& pos = it->pos;
      const double invSigma2 = 1./(sigma*sigma);
      s1 += invSigma2;
      const double x_s = pos.x() * invSigma2;
      sx += x_s;
      sxx += pos.x() * x_s;
      sxy += pos.y() * x_s;
      sxt += (it->time) * x_s;
      const double y_s = pos.y() * invSigma2;
      sy += y_s;
      syy += pos.y() * y_s;
      syt += (it->time) * y_s;
      st += (it->time) * invSigma2;
      //stt += t*t;
    }

    st *= utl::kSpeedOfLight;
    sxt *= utl::kSpeedOfLight;
    syt *= utl::kSpeedOfLight;
    //stt *= utl::kSpeedOfLight*utl::kSpeedOfLight;

    TMatrixDSym mat(3);
    mat(0,0) = sxx; 
    mat(0,1) = sxy; 
    mat(0,2) = -sx;
    mat(1,1) = syy; 
    mat(1,2) = -sy;
    mat(2,2) = s1;

    double det;
    mat.InvertFast(&det);

    if (det==0) {
      INFO("Matrix inversion failed.");
      return EXIT_FAILURE;
    }

    TMatrixD k(3,1);
    k(0,0) = -sxt;
    k(1,0) = -syt;
    k(2,0) = st;

    const TMatrixD p = TMatrixD(mat,TMatrixD::kMult,k);

    const double w2 = 1 - utl::Sqr(p(0,0)) - utl::Sqr(p(1,0));

    if (w2 < 0) {
      INFO("Unphysical direction cosines");
      return EXIT_FAILURE;
    }

    // Return fit results
    u = p(0,0);
    v = p(1,0);
    ct0 = p(2,0);

    return EXIT_SUCCESS;
  }

  int
  GeometryFitter::Fit()
  {

      double uplane, vplane;  // Seeds of the axis direction cosines
      double ct0plane;     // Seed of the 'core time' times the speed of light

      // Fit a plane shower front to set the seeds for u, v, and ct0
      if (FitPlane(uplane, vplane, ct0plane)!=EXIT_SUCCESS) {
        return EXIT_FAILURE;  
      }
      
      // Initializing the cosine directors and the core time with the output from the plane fit
      parSeed.SetValue("u", uplane); 
      parSeed.SetValue("v", vplane); 
      parSeed.SetValue("ct0", ct0plane);      

      GeometryChi2 weightFunction(timeData);  // Create the minimization function
      
      ROOT::Minuit2::MnMigrad migrad(weightFunction, parSeed); // Create the minimizer

      ROOT::Minuit2::FunctionMinimum min = migrad(); // minimize

      if (minuitOutput==true) {
        cout << min;        
      }
            
      // check that the minimization converged
      if ( !min.IsValid() || !min.HasValidCovariance() || min.HesseFailed() || min.HasReachedCallLimit() ) {
        ERROR("Minuit minimize error");
        return EXIT_FAILURE;
      }

      // Retrieve fit parameters
            const ROOT::Minuit2::MnUserParameterState &parameters0 = min.UserParameters();

      // Check that the director cosines are normalised
      const double utemp = parameters0.Parameter(0).Value();
      const double vtemp = parameters0.Parameter(1).Value();
      const double w2 = 1 - utl::Sqr(utemp) - utl::Sqr(vtemp);
      if (w2 < 0) {
        ERROR("  Fit failed due to unphysical angle cosines.");
        return EXIT_FAILURE;
      }

      // Save the fitted parameters and their covariance matrix
            parFitted = ParameterPtr( new ROOT::Minuit2::MnUserParameters( min.UserParameters() ) );
      covariance = CovariancePtr( new ROOT::Minuit2::MnUserCovariance( min.UserCovariance() ) );

      return EXIT_SUCCESS;
  }


  double 
  GeometryFitter::GetTimeResidual(int detectorId) 
    const
  {

    for (std::vector<TimeData>::const_iterator it=timeData.begin(); it!=timeData.end(); ++it) {
      if (it->id==detectorId) {
        return GetTimeResidual(*it);
      }
    }

    ERROR("Detector not found");
    exit(-1); // stop execution
    
  }

  double 
  GeometryFitter::GetPlaneFrontDelay(int detectorId) 
    const
  {

    for (std::vector<TimeData>::const_iterator it=timeData.begin(); it!=timeData.end(); ++it) {
      if (it->id==detectorId) {
        return GetPlaneFrontDelay(*it);
      }
    }

    ERROR("Detector not found");
    exit(-1); // stop execution
    
  }

  double 
  GeometryFitter::GetTimeError(int detectorId) 
    const
  {

    for (std::vector<TimeData>::const_iterator it=timeData.begin(); it!=timeData.end(); ++it) {
      if (it->id==detectorId) {
        return it->sigma;
      }
    }

    ERROR("Detector not found");
    exit(-1); // stop execution
    
  }

  bool 
  GeometryFitter::HasDetector(int detectorId) 
    const
  {

    for (std::vector<TimeData>::const_iterator it=timeData.begin(); it!=timeData.end(); ++it) {
      if (it->id==detectorId) {
        return true;
      }
    }

    return false;
  }

  double 
  GeometryFitter::GetTimeResidual(const TimeData& data) 
    const
  {

    const double t0 = GetCt0() / utl::kSpeedOfLight;

    SphericalShowerFront showerFront(GetAxis(), t0, GetRadius(), GetCore());

    const double residual = data.time - showerFront(data.pos);

    return residual;
  }

  double 
  GeometryFitter::GetPlaneFrontDelay(const TimeData& data) 
    const
  {

    const double t0 = GetCt0() / utl::kSpeedOfLight;

    PlaneShowerFront showerFront(GetAxis(), t0, GetCore());

    const double deviance = data.time - showerFront(data.pos);

    return deviance;
  }

  double
  GeometryFitter::GetMeanTimeResidual()
    const
  {

    double sum = 0;
    for (std::vector<TimeData>::const_iterator it=timeData.begin(); it!=timeData.end(); ++it) {
      sum += GetTimeResidual(*it);
    }

    int n = timeData.size();
    double meanResidual = sum / n;
    
    return meanResidual;
  }

  double
  GeometryFitter::GetChi2()
    const
  {

    double chi2 = 0;
    for (std::vector<TimeData>::const_iterator it=timeData.begin(); it!=timeData.end(); ++it) {
      const double residual = GetTimeResidualSigmas(*it);
      chi2 += utl::Sqr(residual);
    }

    return chi2;
  }

  double
  GeometryFitter::GetNdof()
    const
  {

    int nDetectors = GetNDetectors();

    int nFreeParameters(0); 
    for(int i=0; i<nParameters; ++i) {
      if (parSeed.Parameter(i).IsFixed()==false) {
        nFreeParameters++;  
      }
    }

    int ndof = nDetectors - nFreeParameters;

    return ndof;
  }

  double 
  GeometryFitter::GetTimeResidualSpread() 
    const
  {
    double chi2 = GetChi2();
    int n = timeData.size(); 
    double spread = sqrt(chi2/(n - 1));
    return spread;
  }

  TVector3 
  GeometryFitter::GetAxis()
    const
  {
    return TVector3 (GetU(), GetV(), GetW());
  }

  double 
  GeometryFitter::GetThetaError() 
    const
  {

      double thetaError, phiError, thetaPhiCorrelation;
      const utl::Vector::Triple axis (GetU(), GetV(), GetW());
      const utl::Vector::Triple cova ((*covariance)(0,0), (*covariance)(0,1), (*covariance)(1,1));
      
      utl::PropagateAxisErrors(axis, cova, thetaError, phiError, thetaPhiCorrelation);

      return thetaError;
  }

  double 
  GeometryFitter::GetPhiError() 
    const
  {

      double thetaError, phiError, thetaPhiCorrelation;
      const utl::Vector::Triple axis (GetU(), GetV(), GetW());
      const utl::Vector::Triple cova ((*covariance)(0,0), (*covariance)(0,1), (*covariance)(1,1));
      
      utl::PropagateAxisErrors(axis, cova, thetaError, phiError, thetaPhiCorrelation);

      return phiError;
  }

  double 
  GeometryFitter::GetThetaPhiCorrelation() 
    const
  {

      double thetaError, phiError, thetaPhiCorrelation;
      const utl::Vector::Triple axis (GetU(), GetV(), GetW());
      const utl::Vector::Triple cova ((*covariance)(0,0), (*covariance)(0,1), (*covariance)(1,1));
      
      utl::PropagateAxisErrors(axis, cova, thetaError, phiError, thetaPhiCorrelation);

      return thetaPhiCorrelation;
  }

  TVector3 
  GeometryFitter::GetCore() 
    const 
  {
    double x = parFitted->Value("xcore");
    double y = parFitted->Value("ycore");
    double z = parFitted->Value("zcore");
    
    TVector3 core(x, y, z);

    return core;
  }

  double GeometryFitter::GetRadiusError() 
    const 
  {
    if (IsCurvatureFix()) {
      return 0;
    } else { 
      return parFitted->Error("radius");
    }
  }

  bool GeometryFitter::IsCurvatureFix() 
    const 
  {
    unsigned int i = parSeed.Index("radius"); 
    return parSeed.Parameter(i).IsFixed();
  }

  bool GeometryFitter::IsCoreFix() 
    const 
  {
    unsigned int i = parSeed.Index("xcore"); 
    return parSeed.Parameter(i).IsFixed();
  }

} // end namespace MdGeometryFitterAG