NewtonFitterGSL.h

/**************************************************************************
 * BASF2 (Belle Analysis Framework 2)                                     *
 * See https://github.com/tferber/OrcaKinfit, forked from                 *
 * https://github.com/iLCSoft/MarlinKinfit                                *
 *                                                                        *
 * Further information about the fit engine and the user interface        *
 * provided in MarlinKinfit can be found at                               *
 * https://www.desy.de/~blist/kinfit/doc/html/                            *
 * and in the LCNotes LC-TOOL-2009-001 and LC-TOOL-2009-004 available     *
 * from http://www-flc.desy.de/lcnotes/                                   *
 *                                                                        *
 * Adopted by: Torben Ferber (torben.ferber@desy.de) (TF)                 *
 *                                                                        *
 * This software is provided "as is" without any warranty.                *
 **************************************************************************/
 
#ifndef __NEWTONFITTERGSL_H
#define __NEWTONFITTERGSL_H
 
#include "analysis/OrcaKinFit/BaseFitter.h"
 
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_eigen.h>
 
// Class NewtonFitterGSL
 
namespace Belle2 {
  namespace OrcaKinFit {
 
    class NewtonFitterGSL : public BaseFitter {
    public:
      NewtonFitterGSL();
      virtual ~NewtonFitterGSL();
 
      virtual double fit() override;
 
      virtual int getError() const override;
 
      virtual double getProbability() const override;
      virtual double getChi2() const override;
      virtual int    getDoF() const override;
      virtual int  getIterations() const override;
 
      virtual int    getNcon() const;
 
      virtual int    getNsoft() const;
 
      virtual int    getNpar() const;
 
      virtual int    getNunm() const;
 
      virtual bool initialize() override;
 
      virtual void setDebug(int debuglevel);
 
    protected:
      virtual double calcChi2();
 
      int calcDx();
 
      int calcDxSVD();
 
      void printMy(double M[], double y[], int idim);
 
      bool updateParams(gsl_vector* xnew);
      void fillxold();
      void fillperr();
      int calcM(bool errorpropagation = false);
      int calcy();
      int optimizeScale();
      int invertM();
 
      void calcCovMatrix();
 
      double meritFunction(double mu);
      double meritFunctionDeriv();
 
      enum {NPARMAX = 50, NCONMAX = 10, NUNMMAX = 10};
 
      int npar;      
      int ncon;      
      int nsoft;     
      int nunm;      
      int ierr;      
      int nit;       
 
      double fitprob;   
      double chi2;      
 
      static void ini_gsl_permutation(gsl_permutation*& p, unsigned int size);
      static void ini_gsl_vector(gsl_vector*& v, int unsigned size);
      static void ini_gsl_matrix(gsl_matrix*& m, int unsigned size1, unsigned int size2);
 
      static void debug_print(gsl_matrix* m, const char* name);
      static void debug_print(gsl_vector* v, const char* name);
 
    private:
      unsigned int idim;
      gsl_vector* x;
      gsl_vector* xold;
      gsl_vector* xbest;
      gsl_vector* dx;
      gsl_vector* dxscal;
      gsl_vector* grad;
      gsl_vector* y;
      gsl_vector* yscal;
      gsl_vector* perr;
      gsl_vector* v1;
      gsl_vector* v2;
      gsl_vector* Meval;
 
      gsl_matrix* M;
      gsl_matrix* Mscal;
      gsl_matrix* M1;
      gsl_matrix* M2;
      gsl_matrix* M3;
      gsl_matrix* M4;
      gsl_matrix* M5;
      gsl_matrix* Mevec;
      gsl_matrix* CC;
      gsl_matrix* CC1;
      gsl_matrix* CCinv;
 
      gsl_permutation* permM;
      gsl_eigen_symmv_workspace* ws;
      unsigned int wsdim;
 
      double chi2best;
      double chi2new;
      double chi2old;
      double fvalbest;
      double scale;
      double scalebest;
      double stepsize;
      double stepbest;
      enum {NITMAX = 100};
      double scalevals[NITMAX];
      double fvals[NITMAX];
 
      int imerit;
 
      int debug;
    };
 
  }// end OrcaKinFit namespace
} // end Belle2 namespace
 
#endif // __NEWTONFITTERGSL_H