Belle II Software light-2406-ragdoll
NewFitterGSL.h
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * Forked from https://github.com/iLCSoft/MarlinKinfit *
6 * *
7 * Further information about the fit engine and the user interface *
8 * provided in MarlinKinfit can be found at *
9 * https://www.desy.de/~blist/kinfit/doc/html/ *
10 * and in the LCNotes LC-TOOL-2009-001 and LC-TOOL-2009-004 available *
11 * from http://www-flc.desy.de/lcnotes/ *
12 * *
13 * See git log for contributors and copyright holders. *
14 * This file is licensed under LGPL-3.0, see LICENSE.md. *
15 **************************************************************************/
16
17#ifndef __NEWFITTERGSL_H
18#define __NEWFITTERGSL_H
19
20#include "analysis/OrcaKinFit/BaseFitter.h"
21
22#include <gsl/gsl_vector.h>
23#include <gsl/gsl_matrix.h>
24#include <gsl/gsl_permutation.h>
25#include <gsl/gsl_eigen.h>
26
27namespace Belle2 {
33 namespace OrcaKinFit {
34
35// Class NewFitterGSL
37
51 class NewFitterGSL : public BaseFitter {
52 public:
56 virtual ~NewFitterGSL();
57
59 virtual double fit() override;
60
62 virtual int getError() const override;
63
65 virtual double getProbability() const override;
67 virtual double getChi2() const override;
69 virtual int getDoF() const override;
71 virtual int getIterations() const override;
72
74 virtual int getNcon() const;
75
77 virtual int getNsoft() const;
78
80 virtual int getNpar() const;
81
83 virtual int getNunm() const;
84
86 virtual bool initialize() override;
87
89 virtual void setDebug(int debuglevel);
90
92 virtual int determineLambdas(gsl_vector* vecxnew,
93 const gsl_matrix* MatM,
94 const gsl_vector* vecx,
95 gsl_matrix* MatW,
96 gsl_vector* vecw,
97 double eps = 0
98 );
99
101 virtual void calc2ndOrderCorr(gsl_vector* vecdxhat,
102 const gsl_vector* vecxnew,
103 const gsl_matrix* MatM,
104 gsl_matrix* MatW,
105 gsl_vector* vecw,
106 double eps = 0
107 );
108
110 virtual gsl_matrix_view calcZ(int& rankA,
111 gsl_matrix* MatW1,
112 gsl_matrix* MatW2,
113 gsl_vector* vecw1,
114 gsl_vector* vecw2,
115 gsl_permutation* permW,
116 double eps = 0
117 );
118
120 virtual gsl_matrix_view calcReducedHessian(int& rankH,
121 gsl_matrix* MatW1,
122 const gsl_vector* vecx,
123 gsl_matrix* MatW2,
124 gsl_matrix* MatW3,
125 gsl_vector* vecw1,
126 gsl_vector* vecw2,
127 gsl_permutation* permW,
128 double eps = 0
129 );
130
131 virtual gsl_vector_view calcReducedHessianEigenvalues(int& rankH,
132 gsl_matrix* MatW1,
133 const gsl_vector* vecx,
134 gsl_matrix* MatW2,
135 gsl_matrix* MatW3,
136 gsl_vector* vecw1,
137 gsl_vector* vecw2,
138 gsl_permutation* permW,
139 gsl_eigen_symm_workspace* eigenws,
140 double eps = 0
141 );
142 public:
144 virtual double calcChi2();
145
146 void fillx(gsl_vector* vecx);
147 void fillperr(gsl_vector* vece);
148
149 // Fill matrix MatM, using lambdas from vecx
150 void assembleM(gsl_matrix* MatM, const gsl_vector* vecx, bool errorpropagation = false);
151
152 // Fill matrix MatM with 2nd derivative of Lagrangian, using lambdas from vecx
153 void assembleG(gsl_matrix* MatM, const gsl_vector* vecx);
154
155 // Scale matrix MatM to MatMscal using errors from vece
156 void scaleM(gsl_matrix* MatMscal, const gsl_matrix* MatM, const gsl_vector* vece);
157
158 // Fill vector y, using lambdas from vecx
159 void assembley(gsl_vector* vecy, const gsl_vector* vecx);
160
161 // Scale vector vecy to vecyscal using errors from vece
162 void scaley(gsl_vector* vecyscal, const gsl_vector* vecy, const gsl_vector* vece);
163
164 // Fill chi2 derivatives into vector y
165 int assembleChi2Der(gsl_vector* vecy);
166
167 // Fill vector y with values of constraints
168 void addConstraints(gsl_vector* vecy);
169
170 // Fill constraint derivatives into Matrix M
171 void assembleConstDer(gsl_matrix* MatM);
172
173 // Calculate Newton step update vector vecdx from current point vecx and errors vece
174 int calcNewtonDx(gsl_vector* vecdx,
175 gsl_vector* vecdxscal,
176 gsl_vector* vecx,
177 const gsl_vector* vece,
178 gsl_matrix* MatM,
179 gsl_matrix* MatMscal,
180 gsl_vector* vecy,
181 gsl_vector* vecyscal,
182 gsl_matrix* MatW,
183 gsl_matrix* MatW2,
184 gsl_permutation* permW,
185 gsl_vector* vecw
186 );
187
188 // Calculate limited step after linesearch
189 int calcLimitedDx(double& alpha,
190 double& mu,
191 gsl_vector* vecxnew,
192 int imode,
193 gsl_vector* vecx,
194 gsl_vector* vecdxhat,
195 const gsl_vector* vecdx,
196 const gsl_vector* vecdxscal,
197 const gsl_vector* vece,
198 const gsl_matrix* MatM,
199 const gsl_matrix* MatMscal,
200 gsl_matrix* MatW,
201 gsl_vector* vecw
202 );
203
204 // Perform a line search
205 int doLineSearch(double& alpha,
206 gsl_vector* vecxnew,
207 int imode,
208 double phi0,
209 double dphi0,
210// double phiR, ///< Merit function for given alpha
211 double eta,
212 double zeta,
213 double mu,
214 const gsl_vector* vecx,
215 const gsl_vector* vecdx,
216 const gsl_vector* vece,
217 gsl_vector* vecw
218 );
219
220 // Calculate mu for the merit function
221 double calcMu(const gsl_vector* vecx,
222 const gsl_vector* vece,
223 const gsl_vector* vecdx,
224 const gsl_vector* vecdxscal,
225 const gsl_vector* xnew,
226 const gsl_matrix* MatM,
227 const gsl_matrix* MatMscal,
228 gsl_vector* vecw
229 );
230
231 // Calculate the merit function at x
232 double meritFunction(double mu,
233 const gsl_vector* vecx,
234 const gsl_vector* vece
235 );
236 // Calculate the directional derivative of the merit function at x
237 double meritFunctionDeriv(double mu,
238 const gsl_vector* vecx,
239 const gsl_vector* vece,
240 const gsl_vector* vecdx,
241 gsl_vector* vecw
242 );
243
244 // Transfer values from vecx to FitObjects
245 bool updateParams(gsl_vector* vecx);
246
247
248 int invertM();
249
250 int calcCovMatrix(gsl_matrix* MatW, gsl_permutation* permW, gsl_vector* vecx);
251
252 enum {NPARMAX = 50, NCONMAX = 10, NUNMMAX = 10};
253
254 int npar;
255 int ncon;
256 int nsoft;
257 int nunm;
258 int ierr;
259 int nit;
260
261 double fitprob;
262 double chi2;
263
264 static void ini_gsl_permutation(gsl_permutation*& p, unsigned int size);
265 static void ini_gsl_vector(gsl_vector*& v, int unsigned size);
266 static void ini_gsl_matrix(gsl_matrix*& m, int unsigned size1, unsigned int size2);
267
268 static void debug_print(const gsl_matrix* m, const char* name);
269 static void debug_print(const gsl_vector* v, const char* name);
270
271 // z = x + a y
272 static void add(gsl_vector* vecz, const gsl_vector* vecx, double a, const gsl_vector* vecy);
273
274 // Check whether all elements are finite
275 static bool isfinite(const gsl_vector* vec);
276
277 // Check whether all elements are finite
278 static bool isfinite(const gsl_matrix* mat);
279
281 static void MoorePenroseInverse(gsl_matrix* Ainv,
282 gsl_matrix* A,
283 gsl_matrix* W,
284 gsl_vector* w,
285 double eps = 0
286 );
287 double calcpTLp(const gsl_vector* vecdx,
288 const gsl_matrix* MatM,
289 gsl_vector* vecw
290 );
291
293 int solveSystem(gsl_vector* vecdxscal,
294 double& detW,
295 const gsl_vector* vecyscal,
296 const gsl_matrix* MatMscal,
297 gsl_matrix* MatW,
298 gsl_matrix* MatW2,
299 gsl_vector* vecw,
300 double epsLU,
301 double epsSV
302 );
303
305 int solveSystemLU(gsl_vector* vecdxscal,
306 double& detW,
307 const gsl_vector* vecyscal,
308 const gsl_matrix* MatMscal,
309 gsl_matrix* MatW,
310 gsl_vector* vecw,
311 double eps
312 );
313
315 int solveSystemSVD(gsl_vector* vecdxscal,
316// double& detW,
317 const gsl_vector* vecyscal,
318 const gsl_matrix* MatMscal,
319 gsl_matrix* MatW,
320 gsl_matrix* MatW2,
321 gsl_vector* vecw,
322 double eps
323 );
324
325 public:
326 unsigned int idim;
327 gsl_vector* x;
328 gsl_vector* xold;
329 gsl_vector* xnew;
330// gsl_vector *xbest;
331 gsl_vector* dx;
332 gsl_vector* dxscal;
333// gsl_vector *grad;
334 gsl_vector* y;
335 gsl_vector* yscal;
336 gsl_vector* perr;
337 gsl_vector* v1;
338 gsl_vector* v2;
339// gsl_vector *Meval;
340
341 gsl_matrix* M;
342 gsl_matrix* Mscal;
343 gsl_matrix* W;
344 gsl_matrix* W2;
345 gsl_matrix* W3;
346 // these are only used locally in calcCovMatrix
347 gsl_matrix* M1;
348 gsl_matrix* M2;
349 gsl_matrix* M3;
350 gsl_matrix* M4;
351 gsl_matrix* M5;
352// gsl_matrix *Mevec;
353 gsl_matrix* CC;
354 gsl_matrix* CC1;
355 gsl_matrix* CCinv;
356
357 gsl_permutation* permW;
358 gsl_eigen_symm_workspace* eigenws;
359// gsl_eigen_symmv_workspace *eigenwsv;
360 unsigned int eigenwsdim;
361
362 double chi2best;
363 double chi2new;
364 double chi2old;
365 double fvalbest;
366 double scale;
367 double scalebest;
368 double stepsize;
369 double stepbest;
370 enum {NITMAX = 100};
371 double scalevals[NITMAX];
372 double fvals[NITMAX];
373
374 int imerit;
375 bool try2ndOrderCorr;
376
377 int debug;
378 };
379
380 }// end OrcaKinFit namespace
382} // end Belle2 namespace
383
384#endif // __NEWFITTERGSL_H
Abstract base class for fitting engines of kinematic fits.
Definition: BaseFitter.h:47
A kinematic fitter using the Newton-Raphson method to solve the equations.
Definition: NewFitterGSL.h:51
int ncon
total number of hard constraints
Definition: NewFitterGSL.h:255
int npar
total number of parameters
Definition: NewFitterGSL.h:254
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)
virtual double calcChi2()
Calculate the chi2.
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.
virtual int getNsoft() const
Get the number of soft constraints of the last fit.
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
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.
virtual int 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.
virtual double fit() override
The fit method, returns the fit probability.
virtual bool initialize() override
Initialize the fitter.
int solveSystemSVD(gsl_vector *vecdxscal, 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
virtual int getNcon() const
Get the number of hard constraints of the last fit.
virtual ~NewFitterGSL()
Virtual destructor.
Definition: NewFitterGSL.cc:78
int doLineSearch(double &alpha, gsl_vector *vecxnew, int imode, double phi0, double dphi0, double eta, double zeta, double mu, const gsl_vector *vecx, const gsl_vector *vecdx, const gsl_vector *vece, gsl_vector *vecw)
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)
double fitprob
fit probability
Definition: NewFitterGSL.h:261
int nsoft
total number of soft constraints
Definition: NewFitterGSL.h:256
virtual int getError() const override
Get the error code of the last fit: 0=OK, 1=failed.
int nunm
total number of unmeasured parameters
Definition: NewFitterGSL.h:257
int nit
Number of iterations.
Definition: NewFitterGSL.h:259
virtual int getNunm() const
Get the number of unmeasured parameters of the last fit.
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
double meritFunction(double mu, const gsl_vector *vecx, const gsl_vector *vece)
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)
virtual int getIterations() const override
Get the number of iterations of the last fit.
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...
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.
double meritFunctionDeriv(double mu, const gsl_vector *vecx, const gsl_vector *vece, const gsl_vector *vecdx, gsl_vector *vecw)
virtual double getProbability() const override
Get the fit probability of the last fit.
virtual void setDebug(int debuglevel)
Set the Debug Level.
virtual int getDoF() const override
Get the number of degrees of freedom of the last fit.
virtual double getChi2() const override
Get the chi**2 of the last fit.
virtual int getNpar() const
Get the number of all parameters of the last fit.
double calcpTLp(const gsl_vector *vecdx, const gsl_matrix *MatM, gsl_vector *vecw)
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)
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24