#ifndef MATRIXCALCULATOR_H
#define MATRIXCALCULATOR_H
//_____________________________________________________________________________
/// 
/// MatrixCalculator - performs matrix calculations for CandFitTrackSA 
/// track fitter. Calculates the multiple Coulomb scattering covariance 
/// matrix and the propagator matrix; performs a fit iteration by
/// solving the matrix equation that minimizes chi-square function.
///
/// \author Sergei Avvakumov avva@fnal.gov
///

#include "TMatrixD.h"

#include "DataFT.h"
#include "ConstFT.h"

class FitResult;

class MatrixCalculator {

public:
    MatrixCalculator();
    MatrixCalculator(const AlgConfig& ac, const TrackContext& tc);
    
    ~MatrixCalculator();
    
    ///
    /// The main method called by AlgFitTrackSA - based on the
    /// track data from DataFT object, calculate covariance and
    /// propagator matrices and solve the matrix equation; return
    /// error status - 0 is succes, otherwise failure  
    ///
    Int_t Solve(const DataFT& data);

    ///
    /// get the error matrix of fit parameters
    /// 
    TMatrixD GetFitErrorMatrix() const { return fFitErrM; };
    
    ///
    /// get an [i,j] element of the error matrix
    ///
    double GetFitErrorMatrix(int i, int j) const { return fFitErrM(i,j); };
    
    ///
    /// get the result vector of fit parameters
    ///
    TVectorD GetTrackOut() const { return fTrackOut; };
    
    ///
    /// get the i-th element of the result vector
    ///
    double GetTrackOut(int i) const { return fTrackOut(i); };
    
    ///
    /// get the fit chi-square
    ///
    double   GetChi2() const { return fChi2; };
    
    ///
    /// get the "delta chi-square" between the input
    /// and output vectors of the fit parameters
    /// 
    double   GetDChi2() const { return fDChi2; };
    
    ///
    /// get fit results
    ///
    FitResult GetFitResult() const;
    
private:
    
    ///
    /// calculate the propagator matrix
    ///
    Int_t MakePropagatorMatrix(const DataFT& data);

    ///
    /// calculate the covariance matrix
    ///
    Int_t MakeCovarianceMatrix(const DataFT& data);
    
    ///
    /// calcuate a diagonal element of the covariance matrix
    ///
    Double_t DiagonalElement(Int_t i, Int_t n, 
                        const TVectorD& theta_i, const DataFT& data) const ;
    
    ///
    /// calcuate a non-diagonal element of the covariance matrix
    ///
    Double_t NonDiagonalElement(Int_t i, Int_t k, Int_t n, 
                        const TVectorD& theta_i, const DataFT& data) const;
    
    ///
    /// vector of hit coordinates (single column TMatrixD)
    ///
    TMatrixD               fC;
    
    ///
    /// vector of hit residuals
    ///
    TVectorD               fRes;
    
    ///
    /// propagation matrix
    ///
    TMatrixD               fA;
    
    ///
    /// covariance matrix
    ///
    TMatrixD               fV;
    
    ///
    /// weight matrix (inverse of covariance matrix)
    ///
    TMatrixDSym            fW;

    ///        
    /// fit parameters covariance matrix
    ///
    TMatrixD               fFitCovM;
    
    ///
    /// fit parameters error matrix
    ///
    TMatrixD               fFitErrM;         

    ///
    /// fit chi-square
    ///
    double                 fChi2;
    
    ///
    /// "delta chi-square" between the input and output vectors 
    /// of the fit parameters
    /// 
    double                 fDChi2;           

    ///        
    /// initial vector of the fit parameters (u, du/dz, v, dv/dz, q/p)
    ///
    TVectorD               fTrackIn;
    
    ///         
    /// solution vector of the fit parameters (u, du/dz, v, dv/dz, q/p)
    ///
    TVectorD               fTrackOut;        

    ///
    /// number of planes used in the fit
    ///
    Int_t                  fNPlanesUsed;
};

#endif // MATRIXCALCULATOR_H