Belle II Software  release-08-01-10
NewtonFitterGSL.cc
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 #undef NDEBUG
18 
19 #include "analysis/OrcaKinFit/NewtonFitterGSL.h"
20 
21 #include<iostream>
22 #include<cmath>
23 #include<cassert>
24 
25 #include "analysis/OrcaKinFit/BaseFitObject.h"
26 #include "analysis/OrcaKinFit/BaseHardConstraint.h"
27 #include "analysis/OrcaKinFit/BaseSoftConstraint.h"
28 #include "analysis/OrcaKinFit/BaseTracer.h"
29 #include <framework/logging/Logger.h>
30 
31 #include <gsl/gsl_block.h>
32 #include <gsl/gsl_vector.h>
33 #include <gsl/gsl_matrix.h>
34 #include <gsl/gsl_permutation.h>
35 #include <gsl/gsl_linalg.h>
36 #include <gsl/gsl_blas.h>
37 #include <gsl/gsl_cdf.h>
38 
39 using std::abs;
40 
41 static int nitdebug = 100;
42 static int nitcalc = 0;
43 static int nitsvd = 0;
44 
45 namespace Belle2 {
50  namespace OrcaKinFit {
51 
52 // constructor
53  NewtonFitterGSL::NewtonFitterGSL()
54  : npar(0), ncon(0), nsoft(0), nunm(0), ierr(0), nit(0), fitprob(0), chi2(0),
55  idim(0),
56  x(nullptr), xold(nullptr), xbest(nullptr), dx(nullptr), dxscal(nullptr), grad(nullptr), y(nullptr), yscal(nullptr),
57  perr(nullptr), v1(nullptr), v2(nullptr), Meval(nullptr),
58  M(nullptr), Mscal(nullptr), M1(nullptr), M2(nullptr), M3(nullptr), M4(nullptr), M5(nullptr), Mevec(nullptr),
59  CC(nullptr), CC1(nullptr), CCinv(nullptr), permM(nullptr), ws(nullptr),
60  wsdim(0), chi2best(0),
61  chi2new(0),
62  chi2old(0),
63  fvalbest(0), scale(0), scalebest(0), stepsize(0), stepbest(0),
64  scalevals{}, fvals{},
65  imerit(1),
66  debug(0)
67  {}
68 
69 // destructor
70  NewtonFitterGSL::~NewtonFitterGSL()
71  {
72 
73  if (x) gsl_vector_free(x);
74  if (xold) gsl_vector_free(xold);
75  if (xbest) gsl_vector_free(xbest);
76  if (dx) gsl_vector_free(dx);
77  if (dxscal) gsl_vector_free(dxscal);
78  if (grad) gsl_vector_free(grad);
79  if (y) gsl_vector_free(y);
80  if (yscal) gsl_vector_free(yscal);
81  if (perr) gsl_vector_free(perr);
82  if (v1) gsl_vector_free(v1);
83  if (v2) gsl_vector_free(v2);
84  if (Meval) gsl_vector_free(Meval);
85  if (M) gsl_matrix_free(M);
86  if (Mscal) gsl_matrix_free(Mscal);
87  if (M1) gsl_matrix_free(M1);
88  if (M2) gsl_matrix_free(M2);
89  if (M3) gsl_matrix_free(M3);
90  if (M4) gsl_matrix_free(M4);
91  if (M5) gsl_matrix_free(M5);
92  if (Mevec) gsl_matrix_free(Mevec);
93  if (CC) gsl_matrix_free(CC);
94  if (CC1) gsl_matrix_free(CC1);
95  if (CCinv) gsl_matrix_free(CCinv);
96  if (permM) gsl_permutation_free(permM);
97  if (ws) gsl_eigen_symmv_free(ws);
98 
99  x = nullptr;
100  xold = nullptr;
101  xbest = nullptr;
102  dx = nullptr;
103  dxscal = nullptr;
104  grad = nullptr;
105  y = nullptr;
106  yscal = nullptr;
107  perr = nullptr;
108  v1 = nullptr;
109  v2 = nullptr;
110  Meval = nullptr;
111  M = nullptr;
112  Mscal = nullptr;
113  M1 = nullptr;
114  M2 = nullptr;
115  M3 = nullptr;
116  M4 = nullptr;
117  M5 = nullptr;
118  Mevec = nullptr;
119  CC = nullptr;
120  CC1 = nullptr;
121  CCinv = nullptr;
122  permM = nullptr;
123  ws = nullptr; wsdim = 0;
124  }
125 
126 
127 
128  double NewtonFitterGSL::fit()
129  {
130 
131  // order parameters etc
132  initialize();
133 
134  // initialize eta, etasv, y
135  assert(x && (unsigned int)x->size == idim);
136  assert(dx && (unsigned int)dx->size == idim);
137  assert(y && (unsigned int)y->size == idim);
138  assert(perr && (unsigned int)perr->size == idim);
139  assert(v1 && (unsigned int)v1->size == idim);
140  assert(v2 && (unsigned int)v2->size == idim);
141  assert(Meval && (unsigned int)Meval->size == idim);
142  assert(M && (unsigned int)M->size1 == idim);
143  assert(M1 && (unsigned int)M1->size1 == idim);
144  assert(Mevec && (unsigned int)Mevec->size1 == idim);
145  assert(permM && (unsigned int)permM->size == idim);
146 
147  gsl_vector_set_zero(x);
148  gsl_vector_set_zero(y);
149  gsl_vector_set_all(perr, 1);
150 
151  // Store old x values in xold
152  fillxold();
153  // make sure parameters are consistent
154  updateParams(xold);
155  fillxold();
156  // Get starting values into x
157  gsl_vector_memcpy(x, xold);
158 
159 
160  // LET THE GAMES BEGIN
161 
162 #ifndef FIT_TRACEOFF
163  if (tracer) tracer->initialize(*this);
164 #endif
165 
166  bool converged = false;
167  ierr = 0;
168 
169  chi2new = calcChi2();
170  nit = 0;
171  if (debug > 1) {
172  B2INFO("Fit objects:\n");
173  for (auto fo : fitobjects) {
174  assert(fo);
175  B2INFO(fo->getName() << ": " << *fo << ", chi2=" << fo->getChi2());
176  for (int ilocal = 0; ilocal < fo->getNPar(); ++ilocal) {
177  if (!fo->isParamFixed(ilocal)) {
178  int iglobal = fo->getGlobalParNum(ilocal);
179  B2INFO(" par " << fo->getParamName(ilocal) << ": local: " << ilocal << ": global: " << iglobal
180  << " value=" << fo->getParam(ilocal) << " +- " << fo->getError(ilocal));
181  if (fo->isParamMeasured(ilocal))
182  B2INFO(" measured: " << fo->getMParam(ilocal));
183  } else {
184  B2INFO(" par " << fo->getParamName(ilocal) << ": local: " << ilocal << " -- fixed -- "
185  << " value=" << fo->getParam(ilocal) << " +- " << fo->getError(ilocal));
186 
187  }
188  }
189  }
190  B2INFO("constraints:\n");
191  for (auto i = constraints.begin(); i != constraints.end(); ++i) {
192  BaseHardConstraint* c = *i;
193  assert(c);
194  B2INFO(i - constraints.begin() << " " << c->getName() << ": " << c->getValue() << "+-" << c->getError());
195  int kglobal = c->getGlobalNum();
196  B2INFO(" global number: " << kglobal);
197  }
198  B2INFO("soft constraints:\n");
199  for (auto i = softconstraints.begin(); i != softconstraints.end(); ++i) {
200  BaseSoftConstraint* c = *i;
201  assert(c);
202  B2INFO(i - softconstraints.begin() << " " << c->getName() << ": " << c->getValue() << "+-" << c->getError()
203  << ", chi2: " << c->getChi2());
204  }
205  }
206 
207  do {
208 
209  chi2old = chi2new;
210 
211  if (nit == 0 || nit < nitdebug) {
212  B2DEBUG(11, "===================\nStarting iteration " << nit);
213  std::string printout = "\n Fit objects:\n";
214  for (auto fo : fitobjects) {
215  assert(fo);
216  printout += fo->getName();
217  printout += ", chi2=" ;
218  printout += fo->getChi2();
219  printout += "\n";
220  }
221  printout += "constraints:\n";
222  for (auto c : constraints) {
223  assert(c);
224  printout += c->getName() ;
225  printout += ": ";
226  printout += c->getValue();
227  printout += "+-" ;
228  printout += c->getError() ;
229  printout += "\n";
230  }
231  printout += "soft constraints:\n";
232  for (auto c : softconstraints) {
233  assert(c);
234  printout += c->getName() ;
235  printout += ": ";
236  printout += c->getValue() ;
237  printout += "+-" ;
238  printout += c->getError();
239  printout += ", chi2: ";
240  printout += c->getChi2() ;
241  printout += "\n";
242  }
243  B2DEBUG(12, printout);
244  }
245 
246  // Store old x values in xold
247  fillxold();
248  // Fill errors into perr
249  fillperr();
250 
251  // Compose M:
252  calcM();
253 
254  // Now, calculate the result vector y with the values of the derivatives
255  // d chi^2/d x
256  calcy();
257 
258  if (nit == 0 || nit < nitdebug) {
259  B2DEBUG(13, "After setting up equations: \n");
260  debug_print(M, "M");
261  debug_print(Mscal, "Mscal");
262  debug_print(y, "y");
263  debug_print(yscal, "yscal");
264  debug_print(perr, "perr");
265  debug_print(x, "x");
266  debug_print(xold, "xold");
267  }
268 
269  scalevals[0] = 0;
270  fvals[0] = 0.5 * pow(gsl_blas_dnrm2(yscal), 2);
271  fvalbest = fvals[0];
272  stepsize = 0;
273  scalebest = 0;
274  stepbest = 0;
275 
276  int ifail = calcDx();
277  if (ifail != 0) {
278  B2INFO("NewtonFitterGSL::fit: calcDx: ifail=" << ifail);
279  ierr = 2;
280  break;
281  }
282  // Update values in Fitobjects
283  updateParams(xbest);
284 
285  if (debug > 1) {
286  debug_print(xbest, "new parameters");
287  }
288 
289  calcy();
290  B2DEBUG(13, "New fval: " << 0.5 * pow(gsl_blas_dnrm2(yscal), 2));
291  chi2new = calcChi2();
292  B2DEBUG(13, "chi2: " << chi2old << " -> " << chi2new);
293 
294  if (nit == 0 || nit < nitdebug) {
295  B2DEBUG(13, "After solving equations: \n");
296  debug_print(xbest, "xbest");
297  }
298 
299 
300 // *-- Convergence criteria
301 
302  if (nit < nitdebug) {
303  B2DEBUG(11, "old chi2: " << chi2old << ", new chi2: " << chi2new << ", diff=" << chi2old - chi2new);
304  }
305  ++nit;
306  if (nit > 200) ierr = 1;
307 
308  converged = (abs(chi2new - chi2old) < 0.001 && fvalbest < 1E-3 &&
309  (fvalbest < 1E-6 || abs(fvals[0] - fvalbest) < 0.2 * fvalbest));
310 
311 // if (abs (chi2new - chi2old) >= 0.001)
312 // B2INFO( "abs (chi2new - chi2old)=" << abs (chi2new - chi2old) << " -> try again\n");
313 // if (fvalbest >= 1E-3)
314 // B2INFO("fvalbest=" << fvalbest << " -> try again\n");
315 // if (fvalbest >= 1E-6 && abs(fvals[0]-fvalbest) >= 0.2*fvalbest )
316 // B2INFO("fvalbest=" << fvalbest
317 // << ", abs(fvals[0]-fvalbest)=" << abs(fvals[0]-fvalbest)<< " -> try again\n");
318 // if (stepbest >= 1E-3)
319 // B2INFO("stepbest=" << stepbest << " -> try again\n");
320 // B2INFO( "converged=" << converged);
321  if (converged) {
322  B2DEBUG(10, "abs (chi2new - chi2old)=" << abs(chi2new - chi2old) << "\n"
323  << "fvalbest=" << fvalbest << "\n"
324  << "abs(fvals[0]-fvalbest)=" << abs(fvals[0] - fvalbest) << "\n");
325  }
326 
327 #ifndef FIT_TRACEOFF
328  if (tracer) tracer->step(*this);
329 #endif
330 
331  } while (!(converged || ierr));
332 
333 // *-- End of iterations - calculate errors.
334 
335 // ERROR CALCULATION
336 
337  if (!ierr) {
338 
339  calcCovMatrix();
340 
341  // update errors in fitobjects
342  for (auto& fitobject : fitobjects) {
343  for (int ilocal = 0; ilocal < fitobject->getNPar(); ++ilocal) {
344  int iglobal = fitobject->getGlobalParNum(ilocal);
345  for (int jlocal = ilocal; jlocal < fitobject->getNPar(); ++jlocal) {
346  int jglobal = fitobject->getGlobalParNum(jlocal);
347  if (iglobal >= 0 && jglobal >= 0)
348  fitobject->setCov(ilocal, jlocal, gsl_matrix_get(CCinv, iglobal, jglobal));
349  }//CCinv
350  }
351  }
352  }
353 
354 
355  std::string printout = "\n ========= END =========\n";
356  printout += "Fit objects:\n";
357  for (auto fo : fitobjects) {
358  assert(fo);
359  printout += fo->getName() ;
360  printout += ": ";
361  printout += ", chi2=" ;
362  printout += fo->getChi2() ;
363  printout += "\n";
364  }
365  printout += "constraints:\n";
366  for (auto c : constraints) {
367  assert(c);
368  printout += c->getName();
369  printout += ": ";
370  printout += c->getValue();
371  printout += "+-" ;
372  printout += c->getError();
373  printout += "\n";
374  }
375  printout += "=============================================\n";
376  B2DEBUG(11, printout);
377 
378 // *-- Turn chisq into probability.
379  fitprob = (chi2new >= 0 && ncon + nsoft - nunm > 0) ? gsl_cdf_chisq_Q(chi2new, ncon + nsoft - nunm) : -1;
380 
381 #ifndef FIT_TRACEOFF
382  if (tracer) tracer->finish(*this);
383 #endif
384 
385  B2DEBUG(10, "NewtonFitterGSL::fit: converged=" << converged
386  << ", nit=" << nit << ", fitprob=" << fitprob);
387 
388  if (ierr > 0) fitprob = -1;
389 
390  return fitprob;
391 
392  }
393 
394  bool NewtonFitterGSL::initialize()
395  {
396  covValid = false;
397 
398  // tell fitobjects the global ordering of their parameters:
399  npar = 0;
400  nunm = 0;
401 
402  for (auto& fitobject : fitobjects) {
403  for (int ilocal = 0; ilocal < fitobject->getNPar(); ++ilocal) {
404  if (!fitobject->isParamFixed(ilocal)) {
405  B2DEBUG(13, "NewtonFitterGSL::initialize: parameter " << ilocal
406  << " of fitobject " << fitobject->getName()
407  << " gets global number " << npar);
408  fitobject->setGlobalParNum(ilocal, npar);
409  ++npar;
410  if (!fitobject->isParamMeasured(ilocal)) ++nunm;
411  }
412  }
413  }
414 
415  // set number of constraints
416  ncon = constraints.size();
417  // Tell the constraints their numbers
418  for (int icon = 0; icon < ncon; ++icon) {
419  BaseHardConstraint* c = constraints[icon];
420  assert(c);
421  B2DEBUG(13, "NewtonFitterGSL::initialize: constraint " << c->getName()
422  << " gets global number " << npar + icon);
423  c->setGlobalNum(npar + icon);
424  B2DEBUG(14, "Constraint " << icon << " -> global " << c->getGlobalNum());
425  }
426 
427  nsoft = softconstraints.size();
428 
429  if (nunm > ncon + nsoft) {
430  B2ERROR("NewtonFitterGSL::initialize: nunm=" << nunm << " > ncon+nsoft="
431  << ncon << "+" << nsoft);
432  }
433 
434  idim = npar + ncon;
435 
436  ini_gsl_vector(x, idim);
437  ini_gsl_vector(xold, idim);
438  ini_gsl_vector(xbest, idim);
439  ini_gsl_vector(dx, idim);
440  ini_gsl_vector(dxscal, idim);
441  ini_gsl_vector(grad, idim);
442  ini_gsl_vector(y, idim);
443  ini_gsl_vector(yscal, idim);
444  ini_gsl_vector(perr, idim);
445  ini_gsl_vector(v1, idim);
446  ini_gsl_vector(v2, idim);
447  ini_gsl_vector(Meval, idim);
448 
449  ini_gsl_matrix(M, idim, idim);
450  ini_gsl_matrix(Mscal, idim, idim);
451  ini_gsl_matrix(M1, idim, idim);
452  ini_gsl_matrix(M2, idim, idim);
453  ini_gsl_matrix(M3, idim, idim);
454  ini_gsl_matrix(M4, idim, idim);
455  ini_gsl_matrix(M5, idim, idim);
456  ini_gsl_matrix(Mevec, idim, idim);
457  ini_gsl_matrix(CC, idim, idim);
458  ini_gsl_matrix(CC1, idim, idim);
459  ini_gsl_matrix(CCinv, idim, idim);
460 
461  ini_gsl_permutation(permM, idim);
462 
463  if (ws && wsdim != idim) {
464  gsl_eigen_symmv_free(ws);
465  ws = nullptr;
466  }
467  if (ws == nullptr) ws = gsl_eigen_symmv_alloc(idim);
468  wsdim = idim;
469 
470  return true;
471 
472  }
473 
474  double NewtonFitterGSL::calcChi2()
475  {
476  chi2 = 0;
477  for (auto fo : fitobjects) {
478  assert(fo);
479  chi2 += fo->getChi2();
480  }
481  for (auto bsc : softconstraints) {
482  assert(bsc);
483  chi2 += bsc->getChi2();
484  }
485  return chi2;
486  }
487 
488  void NewtonFitterGSL::printMy(const double Ma[], const double yo[], int idime)
489  {
490  for (int i = 0; i < idime; ++i) {
491  B2INFO(i << " [ " << Ma[idime * i + 0]);
492  for (int j = 1; j < idime; ++j) B2INFO(", " << Ma[idime * i + j]);
493  B2INFO("] [" << yo[i] << "]\n");
494  }
495  }
496 
497  int NewtonFitterGSL::getError() const {return ierr;}
498  double NewtonFitterGSL::getProbability() const {return fitprob;}
499  double NewtonFitterGSL::getChi2() const {return chi2;}
500  int NewtonFitterGSL::getDoF() const {return ncon + nsoft - nunm;}
501  int NewtonFitterGSL::getIterations() const {return nit;}
502 
503  int NewtonFitterGSL::calcDx()
504  {
505  B2DEBUG(11, "entering calcDx");
506  nitcalc++;
507  // from x_(n+1) = x_n - y/y' = x_n - M^(-1)*y we have M*(x_n-x_(n+1)) = y,
508  // which we solve for dx = x_n-x_(n+1) and hence x_(n+1) = x_n-dx
509 
510  gsl_matrix_memcpy(M1, Mscal);
511 
512 
513  int ifail = 0;
514 
515  int signum;
516  int result = gsl_linalg_LU_decomp(M1, permM, &signum);
517  B2DEBUG(11, "calcDx: gsl_linalg_LU_decomp result=" << result);
518  // Solve M1*dx = y
519  ifail = gsl_linalg_LU_solve(M1, permM, yscal, dxscal);
520  B2DEBUG(11, "calcDx: gsl_linalg_LU_solve result=" << ifail);
521 
522  if (ifail != 0) {
523  B2ERROR("NewtonFitter::calcDx: ifail from gsl_linalg_LU_solve=" << ifail);
524 // return calcDxSVD();
525  return -1;
526  }
527  stepsize = std::abs(gsl_vector_get(dxscal, gsl_blas_idamax(dxscal)));
528 
529  // dx = dxscal*perr (component wise)
530  gsl_vector_memcpy(dx, dxscal);
531  gsl_vector_mul(dx, perr);
532 
533  // optimize scale
534 
535  gsl_vector_memcpy(xbest, xold);
536  chi2best = chi2old;
537 
538  optimizeScale();
539 
540  if (scalebest < 0.01) {
541  B2DEBUG(11, "NewtonFitter::calcDx: reverting to calcDxSVD\n");
542  return calcDxSVD();
543  }
544 
545  return 0;
546  }
547 
548  int NewtonFitterGSL::calcDxSVD()
549  {
550 
551  nitsvd++;
552  // from x_(n+1) = x_n - y/y' = x_n - M^(-1)*y we have M*(x_n-x_(n+1)) = y,
553  // which we solve for dx = x_n-x_(n+1) and hence x_(n+1) = x_n-dx
554 
555  for (unsigned int i = 0; i < idim; ++i) assert(gsl_vector_get(perr, i) > 0);
556 
557  // Get eigenvalues and eigenvectors of Mscal
558  int ierrc = 0;
559  gsl_matrix_memcpy(M1, Mscal);
560  B2DEBUG(13, "NewtonFitterGSL::calcDxSVD: Calling gsl_eigen_symmv");
561  ierrc = gsl_eigen_symmv(M1, Meval, Mevec, ws);
562  B2DEBUG(13, "NewtonFitterGSL::calcDxSVD: result of gsl_eigen_symmv: ");
563  if (ierrc != 0) {
564  B2ERROR("NewtonFitter::calcDxSVD: ierr=" << ierrc << "from gsl_eigen_symmv!\n");
565  }
566  // Sort the eigenvalues and eigenvectors in descending order in magnitude
567  ierrc = gsl_eigen_symmv_sort(Meval, Mevec, GSL_EIGEN_SORT_ABS_DESC);
568  if (ierrc != 0) {
569  B2ERROR("NewtonFitter::calcDxSVD: ierr=" << ierrc << "from gsl_eigen_symmv_sort!\n");
570  }
571 
572 
573  // The eigenvectors are stored in the columns of Mevec;
574  // the eigenvectors are orthonormal, therefore Mevec^T = Mevec^-1
575  // Therefore, Mscal = Mevec * diag(Meval)* Mevec^T, and
576  // Mscal^-1 = (Mevec^T)^-1 * diag(Meval)^-1 * Mevec^-1
577  // = Mevec * diag(Meval)^-1 * Mevec^T
578  // So, the solution of M*dx = y is given by
579  // dx = M^-1 * y = Mevec * diag(Meval)^-1 * Mevec^-1 *y
580  // = Mevec * v2
581  // For the pseudoinverse, the last elements of Meveal^-1 are set
582  // to 0, therefore the last elements of v2 will be 0,
583  // therefore we can restrict the calculation of Mevec * v2
584  // to the first ndim rows.
585  // So, we calculate v2 only once, with only the inverse of zero eigenvalues
586  // set to 0, and then calculate Mevec * v2 for fewer and fewer rows
587 
588 
589  // Now M = U * s * V^T
590  // We want to solve M*dx = y, hence dx = V * s^-1 * U^T * y
591  // Calculate UTy first; we need only the first ndim entries
592  unsigned int ndim = 0;
593  while (ndim < idim && gsl_vector_get(Meval, ndim) != 0) ++ndim;
594 
595  if (ndim < idim) {
596  B2INFO("calcDxSVD: idim = " << idim << " > ndim = " << ndim);
597  }
598 
599 
600  // Calculate v2 = 1*Mevec^T*y + 0*v2
601  gsl_blas_dgemv(CblasTrans, 1, Mevec, yscal, 0, v2);
602 
603  // Divide by nonzero eigenvalues
604  for (unsigned int i = 0; i < idim; ++i) {
605  if (double e = gsl_vector_get(Meval, i)) gsl_vector_set(v2, i, gsl_vector_get(v2, i) / e);
606  else gsl_vector_set(v2, i, 0);
607  }
608 
609  stepsize = 0;
610 
611  do {
612  gsl_vector_view v2part = gsl_vector_subvector(v2, 0, ndim);
613  gsl_matrix_view Mevecpart = gsl_matrix_submatrix(Mevec, 0, 0, idim, ndim);
614 
615  // Calculate dx = 1*Mevecpart^T*v2 + 0*dx
616  gsl_blas_dgemv(CblasNoTrans, 1, &Mevecpart.matrix, &v2part.vector, 0, dxscal);
617  // get maximum element
618 // for (unsigned int i = 0; i < idim; ++i) {
619 // if(std::abs(gsl_vector_get (dxscal, i))>stepsize)
620 // stepsize=std::abs(gsl_vector_get (dxscal, i));
621 // }
622  stepsize = std::abs(gsl_vector_get(dxscal, gsl_blas_idamax(dxscal)));
623 
624  // dx = dxscal*perr (component wise)
625  gsl_vector_memcpy(dx, dxscal);
626  gsl_vector_mul(dx, perr);
627 
628  if (debug > 1) {
629  B2INFO("calcDxSVD: Optimizing scale for ndim=" << ndim);
630  debug_print(dxscal, "dxscal");
631  }
632 
633  optimizeScale();
634 
635  --ndim;
636 
637  if (scalebest < 0.01 || ndim < idim - 1) {
638  B2DEBUG(11, "ndim=" << ndim << ", scalebest=" << scalebest);
639  }
640  } while (ndim > 0 && scalebest < 0.01);
641 
642 
643  return 0;
644  }
645 
646  void NewtonFitterGSL::ini_gsl_permutation(gsl_permutation*& p, unsigned int size)
647  {
648  if (p) {
649  if (p->size != size) {
650  gsl_permutation_free(p);
651  if (size > 0) p = gsl_permutation_alloc(size);
652  else p = nullptr;
653  }
654  } else if (size > 0) p = gsl_permutation_alloc(size);
655  }
656 
657  void NewtonFitterGSL::ini_gsl_vector(gsl_vector*& v, unsigned int size)
658  {
659 
660  if (v) {
661  if (v->size != size) {
662  gsl_vector_free(v);
663  if (size > 0) v = gsl_vector_alloc(size);
664  else v = nullptr;
665  }
666  } else if (size > 0) v = gsl_vector_alloc(size);
667  }
668 
669  void NewtonFitterGSL::ini_gsl_matrix(gsl_matrix*& m, unsigned int size1, unsigned int size2)
670  {
671  if (m) {
672  if (m->size1 != size1 || m->size2 != size2) {
673  gsl_matrix_free(m);
674  if (size1 * size2 > 0) m = gsl_matrix_alloc(size1, size2);
675  else m = nullptr;
676  }
677  } else if (size1 * size2 > 0) m = gsl_matrix_alloc(size1, size2);
678  }
679 
680  void NewtonFitterGSL::debug_print(gsl_matrix* m, const char* name)
681  {
682  for (unsigned int i = 0; i < m->size1; ++i)
683  for (unsigned int j = 0; j < m->size2; ++j)
684  if (gsl_matrix_get(m, i, j) != 0)
685  B2INFO(name << "[" << i << "][" << j << "]=" << gsl_matrix_get(m, i, j));
686  }
687 
688  void NewtonFitterGSL::debug_print(gsl_vector* v, const char* name)
689  {
690  for (unsigned int i = 0; i < v->size; ++i)
691  if (gsl_vector_get(v, i) != 0)
692  B2INFO(name << "[" << i << "]=" << gsl_vector_get(v, i));
693  }
694 
695  int NewtonFitterGSL::getNcon() const {return ncon;}
696  int NewtonFitterGSL::getNsoft() const {return nsoft;}
697  int NewtonFitterGSL::getNunm() const {return nunm;}
698  int NewtonFitterGSL::getNpar() const {return npar;}
699 
700  bool NewtonFitterGSL::updateParams(gsl_vector* xnew)
701  {
702  assert(xnew);
703  assert(xnew->size == idim);
704  bool significant = false;
705  for (auto i = fitobjects.begin(); i != fitobjects.end(); ++i) {
706  BaseFitObject* fo = *i;
707  assert(fo);
708  bool s = fo->updateParams(xnew->block->data, xnew->size);
709  significant |= s;
710  if (nit < nitdebug && s) {
711  B2DEBUG(11, "Significant update for FO " << i - fitobjects.begin() << " ("
712  << fo->getName() << ")\n");
713  }
714  }
715  return significant;
716  }
717 
718 
719  void NewtonFitterGSL::fillxold()
720  {
721  assert(xold);
722  assert(xold->size == idim);
723  for (auto fo : fitobjects) {
724  assert(fo);
725  for (int ilocal = 0; ilocal < fo->getNPar(); ++ilocal) {
726  if (!fo->isParamFixed(ilocal)) {
727  int iglobal = fo->getGlobalParNum(ilocal);
728  assert(iglobal >= 0 && iglobal < npar);
729  gsl_vector_set(xold, iglobal, fo->getParam(ilocal));
730  }
731  }
732  }
733  for (auto c : constraints) {
734  assert(c);
735  int iglobal = c->getGlobalNum();
736  assert(iglobal >= 0 && iglobal < (int)idim);
737  gsl_vector_set(xold, iglobal, gsl_vector_get(x, iglobal));
738  }
739  }
740 
741  void NewtonFitterGSL::fillperr()
742  {
743  assert(perr);
744  assert(perr->size == idim);
745  gsl_vector_set_all(perr, 1);
746  for (auto fo : fitobjects) {
747  assert(fo);
748  for (int ilocal = 0; ilocal < fo->getNPar(); ++ilocal) {
749  if (!fo->isParamFixed(ilocal)) {
750  int iglobal = fo->getGlobalParNum(ilocal);
751  assert(iglobal >= 0 && iglobal < npar);
752  double e = std::abs(fo->getError(ilocal));
753  gsl_vector_set(perr, iglobal, e ? e : 1);
754  }
755  }
756  }
757  for (auto c : constraints) {
758  assert(c);
759  int iglobal = c->getGlobalNum();
760  assert(iglobal >= 0 && iglobal < (int)idim);
761  double e = c->getError();
762  gsl_vector_set(perr, iglobal, e ? 1 / e : 1);
763  }
764  }
765 
766  int NewtonFitterGSL::calcM(bool errorpropagation)
767  {
768  assert(M);
769  assert(M->size1 == idim && M->size2 == idim);
770 
771  gsl_matrix_set_zero(M);
772 
773  // First, all terms d^2 chi^2/dx1 dx2
774  for (auto fo : fitobjects) {
775  assert(fo);
776  fo->addToGlobalChi2DerMatrix(M->block->data, M->tda);
777  }
778  if (debug > 3) {
779  B2INFO("After adding covariances from fit objects:\n");
780  debug_print(M, "M");
781  }
782 
783  // Second, all terms d^2 chi^2/dlambda dx,
784  // i.e. the first derivatives of the constraints,
785  // plus the second derivatives times the lambda values
786  for (auto c : constraints) {
787  assert(c);
788  int kglobal = c->getGlobalNum();
789  assert(kglobal >= 0 && kglobal < (int)idim);
790  c->add1stDerivativesToMatrix(M->block->data, M->tda);
791  if (debug > 3) {
792  B2INFO("After adding first derivatives of constraint " << c->getName());
793  debug_print(M, "M");
794  B2INFO("errorpropagation = " << errorpropagation);
795  }
796  // for error propagation after fit,
797  //2nd derivatives of constraints times lambda should _not_ be included!
798  if (!errorpropagation) c->add2ndDerivativesToMatrix(M->block->data, M->tda, gsl_vector_get(x, kglobal));
799  }
800  if (debug > 3) {
801  B2INFO("After adding derivatives of constraints::\n");
802  debug_print(M, "M");
803  B2INFO("===========================================::\n");
804  }
805 
806 
807  // Finally, treat the soft constraints
808 
809  for (auto bsc : softconstraints) {
810  assert(bsc);
811  bsc->add2ndDerivativesToMatrix(M->block->data, M->tda);
812  }
813  if (debug > 3) {
814  B2INFO("After adding soft constraints::\n");
815  debug_print(M, "M");
816  B2INFO("===========================================::\n");
817  }
818 
819  // Rescale columns and rows by perr
820  for (unsigned int i = 0; i < idim; ++i)
821  for (unsigned int j = 0; j < idim; ++j)
822  gsl_matrix_set(Mscal, i, j,
823  gsl_vector_get(perr, i)*gsl_vector_get(perr, j)*gsl_matrix_get(M, i, j));
824 
825  return 0;
826  }
827 
828 
829  int NewtonFitterGSL::calcy()
830  {
831  assert(y);
832  assert(y->size == idim);
833  gsl_vector_set_zero(y);
834  // First, for the parameters
835  for (auto fo : fitobjects) {
836  assert(fo);
837  fo->addToGlobalChi2DerVector(y->block->data, y->size);
838  }
839 
840  // Now add lambda*derivatives of constraints,
841  // And finally, the derivatives w.r.t. to the constraints, i.e. the constraints themselves
842  for (auto c : constraints) {
843  assert(c);
844  int kglobal = c->getGlobalNum();
845  assert(kglobal >= 0 && kglobal < (int)idim);
846  c->addToGlobalChi2DerVector(y->block->data, y->size, gsl_vector_get(x, kglobal));
847  gsl_vector_set(y, kglobal, c->getValue());
848  }
849 
850  // Finally, treat the soft constraints
851 
852  for (auto bsc : softconstraints) {
853  assert(bsc);
854  bsc->addToGlobalChi2DerVector(y->block->data, y->size);
855  }
856  gsl_vector_memcpy(yscal, y);
857  gsl_vector_mul(yscal, perr);
858  return 0;
859  }
860 
861  int NewtonFitterGSL::optimizeScale()
862  {
863  updateParams(xold);
864  calcy();
865  if (debug > 1) {
866  B2INFO("NewtonFitterGSL::optimizeScale");
867  debug_print(xold, "xold");
868  debug_print(yscal, "yscal");
869  debug_print(dx, "dx");
870  debug_print(dxscal, "dxscal");
871  }
872  scalevals[0] = 0;
873  fvals[0] = 0.5 * pow(gsl_blas_dnrm2(yscal), 2);
874  B2DEBUG(11, "NewtonFitterGSL::optimizeScale: fvals[0] = " << fvals[0]);
875  // -dx is a direction
876  // we want to minimize f=0.5*|y|^2 in that direction with y=grad chi2
877  // The gradient grad f is given by grad f = d^chi^2/dxi dxj * dchi^2/dxj
878  // = Mscal*yscal
879 
880  // Calculate grad = 1*Mscal*yscal + 0*grad
881  gsl_blas_dgemv(CblasNoTrans, 1, Mscal, yscal, 0, grad);
882  if (debug > 1) {
883  debug_print(grad, "grad");
884  }
885 
886  // Code adapted from Numerical Recipes (3rd ed), page 479
887  // routine lnsrch
888 
889  int nite = 0;
890 
891  static const double ALF = 1E-4;
892 
893  stepsize = std::abs(gsl_vector_get(dxscal, gsl_blas_idamax(dxscal)));
894  static const double maxstepsize = 5;
895  double scalefactor = maxstepsize / stepsize;
896  if (stepsize > maxstepsize) {
897  gsl_vector_scale(dxscal, maxstepsize / stepsize);
898  B2DEBUG(12, "NewtonFitterGSL::optimizeScale: Rescaling dxscal by factor " << scalefactor);
899  stepsize = std::abs(gsl_vector_get(dxscal, gsl_blas_idamax(dxscal)));
900  if (debug > 1) {
901  debug_print(dxscal, "dxscal");
902  }
903  }
904 
905  double slope;
906  gsl_blas_ddot(dxscal, grad, &slope);
907  slope *= -1;
908  B2DEBUG(12, "NewtonFitterGSL::optimizeScale: slope=" << slope
909  << ", 2*fvals[0]*factor=" << 2 * fvals[0]*scalefactor);
910 
911  scale = 1;
912  double tmpscale, scaleold = 1;
913  do {
914  // sets x = xold
915  gsl_vector_memcpy(x, xold);
916  // x = scale*dx + x
917  if (debug > 1) {
918  debug_print(x, "x(1)");
919  }
920  gsl_blas_daxpy(-scale, dx, x);
921  if (debug > 1) {
922  debug_print(x, "x(2)");
923  }
924 
925  updateParams(x);
926  calcy();
927  if (debug > 1) {
928  debug_print(x, "x(3)");
929  debug_print(yscal, "yscal");
930  }
931  ++nite;
932  scalevals[nite] = scale;
933  fvals[nite] = 0.5 * pow(gsl_blas_dnrm2(yscal), 2);
934 
935  chi2new = calcChi2();
936 
937 
938 // if (chi2new <= chi2best && fvals[nite] <= fvalbest) {
939 // if ((fvals[nite] < fvalbest && chi2new <= chi2best) ||
940 // (fvals[nite] < 1E-4 && chi2new < chi2best)) {
941  if ((fvals[nite] < fvalbest)) {
942  B2DEBUG(13, "new best value: "
943  << " scale " << scalevals[nite] << " -> |y|^2 = " << fvals[nite]
944  << ", chi2=" << chi2new << ", old best chi2: " << chi2best);
945  gsl_vector_memcpy(xbest, x);
946  chi2best = chi2new;
947  fvalbest = fvals[nite];
948  scalebest = scale;
949  stepbest = scale * stepsize;
950  }
951 
952  if (fvals[nite] < fvals[0] + ALF * scale * slope) break;
953  if (nite == 1) {
954  tmpscale = -slope / (2 * (fvals[nite] - fvals[0] - slope));
955  B2DEBUG(13, "quadratic estimate for best scale: " << tmpscale);
956  } else {
957  double rhs1 = fvals[nite] - fvals[0] - scale * slope;
958  double rhs2 = fvals[nite - 1] - fvals[0] - scaleold * slope;
959  double a = (rhs1 / (scale * scale) - rhs2 / (scaleold * scaleold)) / (scale - scaleold);
960  double b = (-scaleold * rhs1 / (scale * scale) + scale * rhs2 / (scaleold * scaleold)) / (scale - scaleold);
961  if (a == 0) tmpscale = -slope / (2 * b);
962  else {
963  double disc = b * b - 3 * a * slope;
964  if (disc < 0) tmpscale = 0.5 * scale;
965  else if (b <= 0) tmpscale = (-b + sqrt(disc)) / (3 * a);
966  else tmpscale = -slope / (b + sqrt(disc));
967  }
968  B2DEBUG(13, "cubic estimate for best scale: " << tmpscale);
969  if (tmpscale > 0.5 * scale) tmpscale = 0.5 * scale;
970  }
971  scaleold = scale;
972  scale = (tmpscale < 0.1 * scale) ? 0.1 * scale : tmpscale;
973  B2DEBUG(11, "New scale: " << scale);
974 
975  } while (nite < NITMAX && scale > 0.0001);
976 
977  if (debug > 1) {
978  for (int it = 0; it <= nite; ++it) {
979  B2INFO(" scale " << scalevals[it] << " -> |y|^2 = " << fvals[it]
980  << " should be " << fvals[0] + ALF * scale * slope);
981  }
982  }
983  return chi2best < chi2old;
984  }
985 
986  int NewtonFitterGSL::invertM()
987  {
988  gsl_matrix_memcpy(M1, M);
989 
990  int ifail = 0;
991 
992  int signum;
993  // Calculate LU decomposition of M into M1
994  int result = gsl_linalg_LU_decomp(M1, permM, &signum);
995  B2DEBUG(11, "invertM: gsl_linalg_LU_decomp result=" << result);
996  // Calculate inverse of M
997  ifail = gsl_linalg_LU_invert(M1, permM, M);
998  B2DEBUG(11, "invertM: gsl_linalg_LU_solve result=" << ifail);
999 
1000  if (ifail != 0) {
1001  B2ERROR("NewtonFitter::invert: ifail from gsl_linalg_LU_invert=" << ifail);
1002  }
1003  return ifail;
1004  }
1005 
1006  void NewtonFitterGSL::setDebug(int debuglevel)
1007  {
1008  debug = debuglevel;
1009  }
1010 
1011 
1012  void NewtonFitterGSL::calcCovMatrix()
1013  {
1014  // Set up equation system M*dadeta + dydeta = 0
1015  // here, dadeta is d a / d eta, i.e. the derivatives of the fitted
1016  // parameters a w.r.t. to the measured parameters eta,
1017  // and dydeta is the derivative of the set of equations
1018  // w.r.t eta, i.e. simply d^2 chi^2 / d a d eta.
1019  // Now, if chi2 = (a-eta)^T*Vinv((a-eta), we have simply
1020  // d^2 chi^2 / d a d eta = - d^2 chi^2 / d a d a
1021  // and can use the method addToGlobalChi2DerMatrix.
1022 
1023  gsl_matrix_set_zero(M1);
1024  gsl_matrix_set_zero(M2);
1025  // First, all terms d^2 chi^2/dx1 dx2
1026  for (auto fo : fitobjects) {
1027  assert(fo);
1028  fo->addToGlobalChi2DerMatrix(M1->block->data, M1->tda);
1029  fo->addToGlobCov(M2->block->data, M2->tda);
1030  }
1031 
1032  // multiply by -1
1033  gsl_matrix_scale(M1, -1);
1034 
1035  // JL: dy/eta are the derivatives of the "objective function" with respect to the MEASURED parameters.
1036  // Since the soft constraints do not depend on the measured, but only on the fitted (!) parameters,
1037  // dy/deta stays -1*M1 also in the presence of soft constraints!
1038  gsl_matrix_view dydeta = gsl_matrix_submatrix(M1, 0, 0, idim, npar);
1039  gsl_matrix_view Cov_eta = gsl_matrix_submatrix(M2, 0, 0, npar, npar);
1040 
1041  if (debug > 3) {
1042  B2INFO("NewtonFitterGSL::calcCovMatrix\n");
1043  debug_print(&dydeta.matrix, "dydeta");
1044  debug_print(&Cov_eta.matrix, "Cov_eta"); \
1045  }
1046 
1047  // JL: calculates d^2 chi^2 / dx1 dx2 + first (!) derivatives of hard & soft constraints, and the
1048  // second derivatives of the soft constraints times the values of the fitted parameters
1049  // - all of the with respect to the FITTED parameters, therefore with soft constraints like in the fit itself.
1050  calcM(true);
1051 
1052  if (debug > 3) {
1053  debug_print(M, "M");
1054  }
1055 
1056  // Now, solve M*dadeta = dydeta
1057 
1058  // Calculate LU decomposition of M into M3
1059  int signum;
1060  int result = gsl_linalg_LU_decomp(M, permM, &signum);
1061 
1062  if (debug > 3) {
1063  B2INFO("invertM: gsl_linalg_LU_decomp result=" << result);
1064  debug_print(M, "M_LU");
1065  }
1066 
1067  // Calculate inverse of M, store in M3
1068  int ifail = gsl_linalg_LU_invert(M, permM, M3);
1069 
1070  if (debug > 3) {
1071  B2INFO("invertM: gsl_linalg_LU_invert ifail=" << ifail);
1072  debug_print(M3, "Minv");
1073  }
1074 
1075  // Calculate dadeta = M3*dydeta
1076  gsl_matrix_set_zero(M4);
1077  gsl_matrix_view dadeta = gsl_matrix_submatrix(M4, 0, 0, idim, npar);
1078 
1079  if (debug > 3) {
1080  debug_print(&dadeta.matrix, "dadeta");
1081  }
1082 
1083  // dadeta = 1*M*dydeta + 0*dadeta
1084  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, M3, &dydeta.matrix, 0, &dadeta.matrix);
1085 
1086 
1087  // Now calculate Cov_a = dadeta*Cov_eta*dadeta^T
1088 
1089  // First, calculate M3 = Cov_eta*dadeta^T as
1090  gsl_matrix_view M3part = gsl_matrix_submatrix(M3, 0, 0, npar, idim);
1091  gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1, &Cov_eta.matrix, &dadeta.matrix, 0, &M3part.matrix);
1092  // Now Cov_a = dadeta*M3part
1093  gsl_matrix_set_zero(M5);
1094  gsl_matrix_view Cov_a = gsl_matrix_submatrix(M5, 0, 0, npar, npar);
1095  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, &dadeta.matrix, &M3part.matrix, 0, M5);
1096  gsl_matrix_memcpy(CCinv, M5);
1097 
1098  if (debug > 3) {
1099  debug_print(&Cov_a.matrix, "Cov_a");
1100  debug_print(CCinv, "full Cov from err prop");
1101  debug_print(M1, "uncorr Cov from err prop");
1102  }
1103 
1104  // Finally, copy covariance matrix
1105  if (cov && covDim != npar) {
1106  delete[] cov;
1107  cov = nullptr;
1108  }
1109  covDim = npar;
1110  if (!cov) cov = new double[covDim * covDim];
1111  for (int i = 0; i < covDim; ++i) {
1112  for (int j = 0; j < covDim; ++j) {
1113  cov[i * covDim + j] = gsl_matrix_get(&Cov_a.matrix, i, j);
1114  }
1115  }
1116  covValid = true;
1117  }
1118 
1119 
1120  double NewtonFitterGSL::meritFunction(double mu)
1121  {
1122  double result = 0;
1123  switch (imerit) {
1124  case 1: // l1 penalty function, Nocedal&Wright Eq. (15.24)
1125  result = chi2;
1126  for (auto c : constraints) {
1127  assert(c);
1128  int kglobal = c->getGlobalNum();
1129  assert(kglobal >= 0 && kglobal < (int)idim);
1130  result += mu * std::fabs(c->getValue());
1131  }
1132  break;
1133  default: assert(0);
1134  }
1135  return result;
1136  }
1137 
1138  double NewtonFitterGSL::meritFunctionDeriv()
1139  {
1140  double result = 0;
1141  switch (imerit) {
1142  case 1:
1143  break;
1144  default: assert(0);
1145  }
1146  return result;
1147  }
1148 
1149  }// end OrcaKinFit namespace
1151 } // end Belle2 namespace
1152 
E
R E
internal precision of FFTW codelets
Definition: DiscreteCosineTransform_31points.cc:36
Belle2::OrcaKinFit::BaseFitObject::updateParams
virtual bool updateParams(double p[], int idim)
Read values from global vector, readjust vector; return: significant change.
Definition: BaseFitObject.cc:196
Belle2::OrcaKinFit::BaseFitObject::getNPar
virtual int getNPar() const =0
Get total number of parameters of this FitObject.
Belle2::OrcaKinFit::BaseFitObject::addToGlobalChi2DerMatrix
virtual void addToGlobalChi2DerMatrix(double *M, int idim) const
Add 2nd derivatives of chi squared to global derivative matrix.
Definition: BaseFitObject.cc:462
Belle2::OrcaKinFit::BaseFitObject::getError
virtual double getError(int ilocal) const
Get error of parameter ilocal.
Definition: BaseFitObject.cc:383
Belle2::OrcaKinFit::BaseFitObject::getGlobalParNum
virtual int getGlobalParNum(int ilocal) const
Get global parameter number of parameter ilocal.
Definition: BaseFitObject.cc:365
Belle2::OrcaKinFit::BaseFitObject::getName
virtual const char * getName() const
Get object's name.
Definition: BaseFitObject.cc:107
Belle2::OrcaKinFit::BaseFitObject::isParamFixed
virtual bool isParamFixed(int ilocal) const
Returns whether parameter is fixed.
Definition: BaseFitObject.cc:408
Belle2::OrcaKinFit::BaseFitObject::getParam
virtual double getParam(int ilocal) const
Get current value of parameter ilocal.
Definition: BaseFitObject.cc:371
Belle2::OrcaKinFit::BaseFitObject::addToGlobalChi2DerVector
virtual int addToGlobalChi2DerVector(double *y, int idim) const
Add derivatives of chi squared to global derivative vector.
Definition: BaseFitObject.cc:482
Belle2::OrcaKinFit::BaseFitObject::addToGlobCov
virtual void addToGlobCov(double *glcov, int idim) const
Add covariance matrix elements to global covariance matrix of size idim x idim.
Definition: BaseFitObject.cc:217
Belle2::OrcaKinFit::BaseFitter::cov
double * cov
global covariance matrix of last fit problem
Definition: BaseFitter.h:104
Belle2::OrcaKinFit::BaseFitter::covValid
bool covValid
Flag whether global covariance is valid.
Definition: BaseFitter.h:105
Belle2::OrcaKinFit::BaseFitter::covDim
int covDim
dimension of global covariance matrix
Definition: BaseFitter.h:103
Belle2::OrcaKinFit::BaseHardConstraint
Abstract base class for constraints of kinematic fits.
Definition: BaseHardConstraint.h:75
Belle2::OrcaKinFit::BaseSoftConstraint
Abstract base class for soft constraints of kinematic fits.
Definition: BaseSoftConstraint.h:71
Belle2::OrcaKinFit::BaseTracer::step
virtual void step(BaseFitter &fitter)
Called at the end of each step.
Definition: BaseTracer.cc:35
Belle2::OrcaKinFit::BaseTracer::initialize
virtual void initialize(BaseFitter &fitter)
Called at the start of a new fit (during initialization)
Definition: BaseTracer.cc:30
Belle2::OrcaKinFit::BaseTracer::finish
virtual void finish(BaseFitter &fitter)
Called at the end of a fit.
Definition: BaseTracer.cc:45
Belle2::OrcaKinFit::NewtonFitterGSL::ncon
int ncon
total number of hard constraints
Definition: NewtonFitterGSL.h:107
Belle2::OrcaKinFit::NewtonFitterGSL::calcDxSVD
int calcDxSVD()
Calculate the vector dx to update the parameters; returns fail code, 0=OK.
Definition: NewtonFitterGSL.cc:548
Belle2::OrcaKinFit::NewtonFitterGSL::npar
int npar
total number of parameters
Definition: NewtonFitterGSL.h:106
Belle2::OrcaKinFit::NewtonFitterGSL::calcChi2
virtual double calcChi2()
Calculate the chi2.
Definition: NewtonFitterGSL.cc:474
Belle2::OrcaKinFit::NewtonFitterGSL::printMy
void printMy(const double M[], const double y[], int idim)
Print a Matrix M and a vector y of dimension idim.
Definition: NewtonFitterGSL.cc:488
Belle2::OrcaKinFit::NewtonFitterGSL::getNsoft
virtual int getNsoft() const
Get the number of soft constraints of the last fit.
Definition: NewtonFitterGSL.cc:696
Belle2::OrcaKinFit::NewtonFitterGSL::fit
virtual double fit() override
The fit method, returns the fit probability.
Definition: NewtonFitterGSL.cc:128
Belle2::OrcaKinFit::NewtonFitterGSL::initialize
virtual bool initialize() override
Initialize the fitter.
Definition: NewtonFitterGSL.cc:394
Belle2::OrcaKinFit::NewtonFitterGSL::getNcon
virtual int getNcon() const
Get the number of hard constraints of the last fit.
Definition: NewtonFitterGSL.cc:695
Belle2::OrcaKinFit::NewtonFitterGSL::nsoft
int nsoft
total number of soft constraints
Definition: NewtonFitterGSL.h:108
Belle2::OrcaKinFit::NewtonFitterGSL::getError
virtual int getError() const override
Get the error code of the last fit: 0=OK, 1=failed.
Definition: NewtonFitterGSL.cc:497
Belle2::OrcaKinFit::NewtonFitterGSL::nunm
int nunm
total number of unmeasured parameters
Definition: NewtonFitterGSL.h:109
Belle2::OrcaKinFit::NewtonFitterGSL::calcDx
int calcDx()
Calculate the vector dx to update the parameters; returns fail code, 0=OK.
Definition: NewtonFitterGSL.cc:503
Belle2::OrcaKinFit::NewtonFitterGSL::getNunm
virtual int getNunm() const
Get the number of unmeasured parameters of the last fit.
Definition: NewtonFitterGSL.cc:697
Belle2::OrcaKinFit::NewtonFitterGSL::getIterations
virtual int getIterations() const override
Get the number of iterations of the last fit.
Definition: NewtonFitterGSL.cc:501
Belle2::OrcaKinFit::NewtonFitterGSL::getProbability
virtual double getProbability() const override
Get the fit probability of the last fit.
Definition: NewtonFitterGSL.cc:498
Belle2::OrcaKinFit::NewtonFitterGSL::setDebug
virtual void setDebug(int debuglevel)
Set the Debug Level.
Definition: NewtonFitterGSL.cc:1006
Belle2::OrcaKinFit::NewtonFitterGSL::getDoF
virtual int getDoF() const override
Get the number of degrees of freedom of the last fit.
Definition: NewtonFitterGSL.cc:500
Belle2::OrcaKinFit::NewtonFitterGSL::getChi2
virtual double getChi2() const override
Get the chi**2 of the last fit.
Definition: NewtonFitterGSL.cc:499
Belle2::OrcaKinFit::NewtonFitterGSL::getNpar
virtual int getNpar() const
Get the number of all parameters of the last fit.
Definition: NewtonFitterGSL.cc:698
Belle2::OrcaKinFit::NewtonFitterGSL::~NewtonFitterGSL
virtual ~NewtonFitterGSL()
Virtual destructor.
Definition: NewtonFitterGSL.cc:70
Belle2::sqrt
double sqrt(double a)
sqrt for double
Definition: beamHelpers.h:28
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17