27 #ifndef __NEWFITTERGSL_H
28 #define __NEWFITTERGSL_H
32 #include <gsl/gsl_vector.h>
33 #include <gsl/gsl_matrix.h>
34 #include <gsl/gsl_permutation.h>
35 #include <gsl/gsl_eigen.h>
71 virtual int getDoF()
const;
91 virtual void setDebug (
int debuglevel);
95 const gsl_matrix *MatM,
96 const gsl_vector *vecx,
104 const gsl_vector *vecxnew,
105 const gsl_matrix *MatM,
112 virtual gsl_matrix_view
calcZ (
int& rankA,
117 gsl_permutation *permW,
124 const gsl_vector *vecx,
129 gsl_permutation *permW,
135 const gsl_vector *vecx,
140 gsl_permutation *permW,
141 gsl_eigen_symm_workspace *eigenws,
148 void fillx (gsl_vector *vecx);
149 void fillperr(gsl_vector *vece);
152 void assembleM (gsl_matrix *MatM,
const gsl_vector *vecx,
bool errorpropagation =
false);
155 void assembleG (gsl_matrix *MatM,
const gsl_vector *vecx);
158 void scaleM (gsl_matrix *MatMscal,
const gsl_matrix *MatM,
const gsl_vector *vece);
161 void assembley(gsl_vector *vecy,
const gsl_vector *vecx);
164 void scaley (gsl_vector *vecyscal,
const gsl_vector *vecy,
const gsl_vector *vece);
167 void assembleChi2Der(gsl_vector *vecy);
170 void addConstraints(gsl_vector *vecy);
173 void assembleConstDer(gsl_matrix *MatM);
177 gsl_vector *vecdxscal,
179 const gsl_vector *vece,
181 gsl_matrix *MatMscal,
183 gsl_vector *vecyscal,
186 gsl_permutation *permW,
196 gsl_vector *vecdxhat,
197 const gsl_vector *vecdx,
198 const gsl_vector *vecdxscal,
199 const gsl_vector *vece,
200 const gsl_matrix *MatM,
201 const gsl_matrix *MatMscal,
216 const gsl_vector *vecx,
217 const gsl_vector *vecdx,
218 const gsl_vector *vece,
223 double calcMu (
const gsl_vector *vecx,
224 const gsl_vector *vece,
225 const gsl_vector *vecdx,
226 const gsl_vector *vecdxscal,
227 const gsl_vector *xnew,
228 const gsl_matrix *MatM,
229 const gsl_matrix *MatMscal,
235 const gsl_vector *vecx,
236 const gsl_vector *vece
240 const gsl_vector *vecx,
241 const gsl_vector *vece,
242 const gsl_vector *vecdx,
247 bool updateParams (gsl_vector *vecx);
252 void calcCovMatrix(gsl_matrix *MatW, gsl_permutation *permW, gsl_vector *vecx);
254 enum {NPARMAX=50, NCONMAX=10, NUNMMAX=10};
266 static void ini_gsl_permutation (gsl_permutation *&p,
unsigned int size);
267 static void ini_gsl_vector (gsl_vector *&v,
int unsigned size);
268 static void ini_gsl_matrix (gsl_matrix *&m,
int unsigned size1,
unsigned int size2);
270 static void debug_print (
const gsl_matrix *m,
const char *name);
271 static void debug_print (
const gsl_vector *v,
const char *name);
274 static void add (gsl_vector *vecz,
const gsl_vector *vecx,
double a,
const gsl_vector *vecy);
277 static bool isfinite (
const gsl_vector *vec);
280 static bool isfinite (
const gsl_matrix *mat);
289 double calcpTLp (
const gsl_vector *vecdx,
290 const gsl_matrix *MatM,
297 const gsl_vector *vecyscal,
298 const gsl_matrix *MatMscal,
309 const gsl_vector *vecyscal,
310 const gsl_matrix *MatMscal,
319 const gsl_vector *vecyscal,
320 const gsl_matrix *MatMscal,
359 gsl_permutation *permW;
360 gsl_eigen_symm_workspace *eigenws;
362 unsigned int eigenwsdim;
373 double scalevals[NITMAX];
374 double fvals[NITMAX];
377 bool try2ndOrderCorr;
382 #endif // __NEWFITTERGSL_H
virtual int getNcon() const
Get the number of hard constraints of the last fit.
Definition: NewFitterGSL.cc:463
int nunm
total number of unmeasured parameters
Definition: NewFitterGSL.h:259
virtual gsl_vector_view calcReducedHessianEigenvalues(int &rankH, gsl_matrix *MatW1, const gsl_vector *vecx, gsl_matrix *MatW2, gsl_matrix *MatW3, gsl_vector *vecw1, gsl_vector *vecw2, gsl_permutation *permW, gsl_eigen_symm_workspace *eigenws, double eps=0)
Definition: NewFitterGSL.cc:1885
double chi2
final chi2
Definition: NewFitterGSL.h:264
virtual void determineLambdas(gsl_vector *vecxnew, const gsl_matrix *MatM, const gsl_vector *vecx, gsl_matrix *MatW, gsl_vector *vecw, double eps=0)
Determine best lambda values.
Definition: NewFitterGSL.cc:1456
Abstract base class for fitting engines of kinematic fits.
Definition: BaseFitter.h:63
virtual gsl_matrix_view calcZ(int &rankA, gsl_matrix *MatW1, gsl_matrix *MatW2, gsl_vector *vecw1, gsl_vector *vecw2, gsl_permutation *permW, double eps=0)
Calculate null space of constraints; return value is a matrix view of MatW1, with column vectors span...
Definition: NewFitterGSL.cc:1807
int ncon
total number of hard constraints
Definition: NewFitterGSL.h:257
virtual double fit()
The fit method, returns the fit probability.
Definition: NewFitterGSL.cc:114
int solveSystem(gsl_vector *vecdxscal, double &detW, const gsl_vector *vecyscal, const gsl_matrix *MatMscal, gsl_matrix *MatW, gsl_matrix *MatW2, gsl_vector *vecw, double epsLU, double epsSV)
solve system of equations Mscal*dxscal = yscal
Definition: NewFitterGSL.cc:1680
int solveSystemLU(gsl_vector *vecdxscal, double &detW, const gsl_vector *vecyscal, const gsl_matrix *MatMscal, gsl_matrix *MatW, gsl_vector *vecw, double eps)
solve system of equations Mscal*dxscal = yscal using LU decomposition
Definition: NewFitterGSL.cc:1714
int calcNewtonDx(gsl_vector *vecdx, gsl_vector *vecdxscal, gsl_vector *vecx, const gsl_vector *vece, gsl_matrix *MatM, gsl_matrix *MatMscal, gsl_vector *vecy, gsl_vector *vecyscal, gsl_matrix *MatW, gsl_matrix *MatW2, gsl_permutation *permW, gsl_vector *vecw)
Definition: NewFitterGSL.cc:767
double meritFunction(double mu, const gsl_vector *vecx, const gsl_vector *vece)
Definition: NewFitterGSL.cc:1228
double calcpTLp(const gsl_vector *vecdx, const gsl_matrix *MatM, gsl_vector *vecw)
Definition: NewFitterGSL.cc:1590
Declares class BaseFitter.
virtual double getChi2() const
Get the chi**2 of the last fit.
Definition: NewFitterGSL.cc:398
A kinematic fitter using the Newton-Raphson method to solve the equations.
Definition: NewFitterGSL.h:53
virtual int getNpar() const
Get the number of all parameters of the last fit.
Definition: NewFitterGSL.cc:466
int nsoft
total number of soft constraints
Definition: NewFitterGSL.h:258
virtual bool initialize()
Initialize the fitter.
Definition: NewFitterGSL.cc:300
int ierr
Error status.
Definition: NewFitterGSL.h:260
double fitprob
fit probability
Definition: NewFitterGSL.h:263
virtual int getNunm() const
Get the number of unmeasured parameters of the last fit.
Definition: NewFitterGSL.cc:465
int calcLimitedDx(double &alpha, double &mu, gsl_vector *vecxnew, int imode, gsl_vector *vecx, gsl_vector *vecdxhat, const gsl_vector *vecdx, const gsl_vector *vecdxscal, const gsl_vector *vece, const gsl_matrix *MatM, const gsl_matrix *MatMscal, gsl_matrix *MatW, gsl_vector *vecw)
Definition: NewFitterGSL.cc:880
int nit
Number of iterations.
Definition: NewFitterGSL.h:261
NewFitterGSL()
Constructor.
Definition: NewFitterGSL.cc:58
virtual int getDoF() const
Get the number of degrees of freedom of the last fit.
Definition: NewFitterGSL.cc:399
virtual ~NewFitterGSL()
Virtual destructor.
Definition: NewFitterGSL.cc:79
virtual double calcChi2()
Calculate the chi2.
Definition: NewFitterGSL.cc:381
double meritFunctionDeriv(double mu, const gsl_vector *vecx, const gsl_vector *vece, const gsl_vector *vecdx, gsl_vector *vecw)
Definition: NewFitterGSL.cc:1264
int npar
total number of parameters
Definition: NewFitterGSL.h:256
virtual void calc2ndOrderCorr(gsl_vector *vecdxhat, const gsl_vector *vecxnew, const gsl_matrix *MatM, gsl_matrix *MatW, gsl_vector *vecw, double eps=0)
Calculate 2nd order correction step.
Definition: NewFitterGSL.cc:1609
virtual void setDebug(int debuglevel)
Set the Debug Level.
Definition: NewFitterGSL.cc:1342
int doLineSearch(double &alpha, gsl_vector *vecxnew, int imode, double phi0, double dphi0, double phiR, double eta, double zeta, double mu, const gsl_vector *vecx, const gsl_vector *vecdx, const gsl_vector *vece, gsl_vector *vecw)
Definition: NewFitterGSL.cc:1018
static void MoorePenroseInverse(gsl_matrix *Ainv, gsl_matrix *A, gsl_matrix *W, gsl_vector *w, double eps=0)
Compute the Moore-Penrose pseudo-inverse A+ of A, using SVD.
Definition: NewFitterGSL.cc:1547
virtual int getError() const
Get the error code of the last fit: 0=OK, 1=failed.
Definition: NewFitterGSL.cc:396
virtual gsl_matrix_view calcReducedHessian(int &rankH, gsl_matrix *MatW1, const gsl_vector *vecx, gsl_matrix *MatW2, gsl_matrix *MatW3, gsl_vector *vecw1, gsl_vector *vecw2, gsl_permutation *permW, double eps=0)
Calculate reduced Hessian; return value is a matrix view of MatW1 giving the reduced Hessian...
Definition: NewFitterGSL.cc:1842
virtual double getProbability() const
Get the fit probability of the last fit.
Definition: NewFitterGSL.cc:397
virtual int getIterations() const
Get the number of iterations of the last fit.
Definition: NewFitterGSL.cc:400
virtual int getNsoft() const
Get the number of soft constraints of the last fit.
Definition: NewFitterGSL.cc:464
int solveSystemSVD(gsl_vector *vecdxscal, double &detW, const gsl_vector *vecyscal, const gsl_matrix *MatMscal, gsl_matrix *MatW, gsl_matrix *MatW2, gsl_vector *vecw, double eps)
solve system of equations Mscal*dxscal = yscal using SVD decomposition
Definition: NewFitterGSL.cc:1767
double calcMu(const gsl_vector *vecx, const gsl_vector *vece, const gsl_vector *vecdx, const gsl_vector *vecdxscal, const gsl_vector *xnew, const gsl_matrix *MatM, const gsl_matrix *MatMscal, gsl_vector *vecw)
Definition: NewFitterGSL.cc:1136