GeometryFitter.h

Go to the documentation of this file.

#ifndef _MdGeometryFitterAG_GeometryFitter_h_
#define _MdGeometryFitterAG_GeometryFitter_h_

#include <vector>
#include <memory>
#include <bitset> 

#include <TVector3.h>
#include <TMatrixDSym.h>
#include <Minuit2/MnUserParameters.h>
#include <Minuit2/MnUserCovariance.h>

#include "TimeData.h"

namespace MdGeometryFitterAG {

  class GeometryFitter {

    public:

      GeometryFitter();

      //  Add a detector to the data used in the fitting  
      void AddDetector(int id, TVector3 pos, double time, double sigma); 

      void RemoveDetectors() {timeData.clear();}
      int GetNDetectors() const {return timeData.size();}
      bool HasDetector(int detectorId) const;

      // Parameter configuration
      void SetRadius(double radius) {parSeed.SetValue("radius", radius);}
      void SetCurvatureFree() {parSeed.Release("radius");}
      void SetCurvatureFix() {parSeed.Fix("radius");}
      bool IsCurvatureFix() const;

      void SetCore(const TVector3 &core);
      void SetCoreFree();
      void SetCoreFix();
      bool IsCoreFix() const;

      int Fit();
      void SetMinuitOutput(bool flag=true) {minuitOutput=flag;} 

      // After the minimisation the fit parameters are saved in the parFitted data member of this class 
      // Fitter parameters and their errors can be retrieved afterwards using the following getters
      double GetU() const {return parFitted->Value("u");}   
      double GetV() const {return parFitted->Value("v");}
      double GetW() const {return sqrt(1-GetU()*GetU()-GetV()*GetV());}

      TVector3 GetAxis() const;
      double GetTheta() const {return acos(GetW());}
      double GetThetaError() const;
      double GetPhi() const {return atan2(GetV(),GetU());}
      double GetPhiError() const;
      double GetThetaPhiCorrelation() const;

      double GetCt0() const {return parFitted->Value("ct0");}
      double GetCt0Error() const {return parFitted->Error("ct0");}

      double GetRadius() const {return parFitted->Value("radius");}
      double GetRadiusError() const;

      TVector3 GetCore() const;

      // Residual stuff
      double GetTimeResidual(int detectorId) const;
      double GetTimeError(int detectorId) const;
      double GetMeanTimeResidual() const;
      double GetChi2() const;
      double GetNdof() const;
      double GetTimeResidualSpread() const;

      // Delay in nanoseconds of a detector signal wrt to a plane shower front passing through the core position 
      // and oriented with the reconstructed axis. 
      double GetPlaneFrontDelay(int detectorId) const;   

    private:

      // List of detector data used in the geometrical reconstruction: id, position, time observable value, and time observable error
      std::vector<TimeData> timeData;

      bool minuitOutput;

      // Fit parameters and their errors
      int nParameters;
      typedef std::unique_ptr<ROOT::Minuit2::MnUserParameters> ParameterPtr;  
      ParameterPtr parFitted;   // seed and fitted parameters
      ROOT::Minuit2::MnUserParameters parSeed;

      typedef std::unique_ptr<ROOT::Minuit2::MnUserCovariance> CovariancePtr;
      CovariancePtr covariance;  

  // analytic linear fit in local CS; approximation by neglecting
  // z-coordinate, since z ~ 0 in local bary coordinate system,
  // therefore loosing the w component of the axis. 
      int FitPlane(double &u, double &v, double &ct0) const;   
      double GetTimeResidual(const TimeData& data) const;

      // Studentised residual (i.e. residual over sigma)
      double GetTimeResidualSigmas(const TimeData& data) const
      { return GetTimeResidual(data)/data.sigma; }

      double GetPlaneFrontDelay(const TimeData& data) const;   

  }; // end GeometryFitter class


} // end namespace MdGeometryFitterAG

#endif