Belle II Software  release-06-00-14
OPALFitterGSL.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 #include<iostream>
18 #include<cmath>
19 #include<cassert>
20 
21 #include "analysis/OrcaKinFit/OPALFitterGSL.h"
22 
23 #include "analysis/OrcaKinFit/BaseFitObject.h"
24 #include "analysis/OrcaKinFit/BaseHardConstraint.h"
25 #include "analysis/OrcaKinFit/BaseTracer.h"
26 #include <framework/logging/Logger.h>
27 
28 #include <gsl/gsl_block.h>
29 #include <gsl/gsl_vector.h>
30 #include <gsl/gsl_matrix.h>
31 #include <gsl/gsl_permutation.h>
32 #include <gsl/gsl_linalg.h>
33 #include <gsl/gsl_blas.h>
34 #include <gsl/gsl_cdf.h>
35 
36 using std::endl;
37 using std::abs;
38 
39 //using namespace Belle2 {
40 //using namespace OrcaKinFit {
41 
42 namespace Belle2 {
47  namespace OrcaKinFit {
48 
49 // constructor
50  OPALFitterGSL::OPALFitterGSL()
51  : npar(0), nmea(0), nunm(0), ncon(0), ierr(0), nit(0),
52  fitprob(0), chi2(0),
53  f(nullptr), r(nullptr), Fetaxi(nullptr), S(nullptr), Sinv(nullptr), SinvFxi(nullptr), SinvFeta(nullptr),
54  W1(nullptr), G(nullptr), H(nullptr), HU(nullptr), IGV(nullptr), V(nullptr), VLU(nullptr), Vinv(nullptr), Vnew(nullptr),
55  Minv(nullptr), dxdt(nullptr), Vdxdt(nullptr),
56  dxi(nullptr), Fxidxi(nullptr), lambda(nullptr), FetaTlambda(nullptr),
57  etaxi(nullptr), etasv(nullptr), y(nullptr), y_eta(nullptr), Vinvy_eta(nullptr), FetaV(nullptr),
58  permS(nullptr), permU(nullptr), permV(nullptr), debug(0)
59  {}
60 
61 // destructor
62  OPALFitterGSL::~OPALFitterGSL()
63  {
64  if (f) gsl_vector_free(f);
65  if (r) gsl_vector_free(r);
66  if (Fetaxi) gsl_matrix_free(Fetaxi);
67  if (S) gsl_matrix_free(S);
68  if (Sinv) gsl_matrix_free(Sinv);
69  if (SinvFxi) gsl_matrix_free(SinvFxi);
70  if (SinvFeta) gsl_matrix_free(SinvFeta);
71  if (W1) gsl_matrix_free(W1);
72  if (G) gsl_matrix_free(G);
73  if (H) gsl_matrix_free(H);
74  if (HU) gsl_matrix_free(HU);
75  if (IGV) gsl_matrix_free(IGV);
76  if (V) gsl_matrix_free(V);
77  if (VLU) gsl_matrix_free(VLU);
78  if (Vinv) gsl_matrix_free(Vinv);
79  if (Vnew) gsl_matrix_free(Vnew);
80  if (Minv) gsl_matrix_free(Minv);
81  if (dxdt) gsl_matrix_free(dxdt);
82  if (Vdxdt) gsl_matrix_free(Vdxdt);
83  if (dxi) gsl_vector_free(dxi);
84  if (Fxidxi) gsl_vector_free(Fxidxi);
85  if (lambda) gsl_vector_free(lambda);
86  if (FetaTlambda) gsl_vector_free(FetaTlambda);
87  if (etaxi) gsl_vector_free(etaxi);
88  if (etasv) gsl_vector_free(etasv);
89  if (y) gsl_vector_free(y);
90  if (y_eta) gsl_vector_free(y_eta);
91  if (Vinvy_eta) gsl_vector_free(Vinvy_eta);
92  if (FetaV) gsl_matrix_free(FetaV);
93  if (permS) gsl_permutation_free(permS);
94  if (permU) gsl_permutation_free(permU);
95  if (permV) gsl_permutation_free(permV);
96 
97  f = nullptr;
98  r = nullptr;
99  Fetaxi = nullptr;
100  S = nullptr;
101  Sinv = nullptr;
102  SinvFxi = nullptr;
103  SinvFeta = nullptr;
104  W1 = nullptr;
105  G = nullptr;
106  H = nullptr;
107  HU = nullptr;
108  IGV = nullptr;
109  V = nullptr;
110  VLU = nullptr;
111  Vinv = nullptr;
112  Vnew = nullptr;
113  Minv = nullptr;
114  dxdt = nullptr;
115  Vdxdt = nullptr;
116  dxi = nullptr;
117  Fxidxi = nullptr;
118  lambda = nullptr;
119  FetaTlambda = nullptr;
120  etaxi = nullptr;
121  etasv = nullptr;
122  y = nullptr;
123  y_eta = nullptr;
124  Vinvy_eta = nullptr;
125  FetaV = nullptr;
126  permS = nullptr;
127  permU = nullptr;
128  permV = nullptr;
129  }
130 
131 // do it (~ transcription of WWFGO as of ww113)
132  double OPALFitterGSL::fit()
133  {
134 
135 
136 
137  //
138  // ( ) ^ ^
139  // ( eta ) nmea |
140  // ( ) v |
141  // etaxi = (-----) --- npar
142  // ( ) ^ |
143  // ( xi ) nunm |
144  // ( ) v v
145 
146  // <- ncon ->
147  // ( ) ^ ^
148  // ( Feta^T ) nmea |
149  // ( ) v |
150  // Fetaxi^T = ( -------- ) --- npar
151  // ( ) ^ |
152  // ( Fxi^T ) nunm |
153  // ( ) v v
154  //
155  // Fetaxi are the derivatives of the constraints wrt the fitted parameters, thus A_theta/xi in book
156 
157 
158  //
159  // <- nmea ->|<- nunm ->
160  // ( | ) ^ ^
161  // ( Vetaeta | Vetaxi ) nmea |
162  // ( | ) v |
163  // V = (----------+----------) --- npar
164  // ( | ) ^ |
165  // ( Vxieta | Vxixi ) nunm |
166  // ( | ) v v
167  //
168  // V is the covariance matrix. Before the fit only Vetaeta is non-zero (covariance of measured parameters)
169 
170  //
171  // <- nmea ->|<- nunm ->
172  // ( | ) ^
173  // dxdt^T = ( detadt | dxidt ) nmea
174  // ( | ) v
175  //
176  // dxdt are the partial derivates of the fitted parameters wrt the measured parameters,
177  //
178  // <- nmea ->|<- nunm ->
179  // ( | ) ^
180  // Vdxdt = ( Vdetadt | Vdxidt ) nmea
181  // ( | ) v
182  //
183  // Vdxdt is Vetaeta * dxdt^T, thus Vdxdt[nmea][npar]
184  // both needed for calculation of the full covariance matrix of the fitted parameters
185 
186  // order parameters etc
187  initialize();
188 
189  assert(f && (int)f->size == ncon);
190  assert(r && (int)r->size == ncon);
191  assert(Fetaxi && (int)Fetaxi->size1 == ncon && (int)Fetaxi->size2 == npar);
192  assert(S && (int)S->size1 == ncon && (int)S->size2 == ncon);
193  assert(Sinv && (int)Sinv->size1 == ncon && (int)Sinv->size2 == ncon);
194  assert(nunm == 0 || (SinvFxi && (int)SinvFxi->size1 == ncon && (int)SinvFxi->size2 == nunm));
195  assert(SinvFeta && (int)SinvFeta->size1 == ncon && (int)SinvFeta->size2 == nmea);
196  assert(nunm == 0 || (W1 && (int)W1->size1 == nunm && (int)W1->size2 == nunm));
197  assert(G && (int)G->size1 == nmea && (int)G->size2 == nmea);
198  assert(nunm == 0 || (H && (int)H->size1 == nmea && (int)H->size2 == nunm));
199  assert(nunm == 0 || (HU && (int)HU->size1 == nmea && (int)HU->size2 == nunm));
200  assert(IGV && (int)IGV->size1 == nmea && (int)IGV->size2 == nmea);
201  assert(V && (int)V->size1 == npar && (int)V->size2 == npar);
202  assert(VLU && (int)VLU->size1 == nmea && (int)VLU->size2 == nmea);
203  assert(Vinv && (int)Vinv->size1 == nmea && (int)Vinv->size2 == nmea);
204  assert(Vnew && (int)Vnew->size1 == npar && (int)Vnew->size2 == npar);
205  assert(Minv && (int)Minv->size1 == npar && (int)Minv->size2 == npar);
206  assert(dxdt && (int)dxdt->size1 == npar && (int)dxdt->size2 == nmea);
207  assert(Vdxdt && (int)Vdxdt->size1 == nmea && (int)Vdxdt->size2 == npar);
208  assert(nunm == 0 || (dxi && (int)dxi->size == nunm));
209  assert(nunm == 0 || (Fxidxi && (int)Fxidxi->size == ncon));
210  assert(lambda && (int)lambda->size == ncon);
211  assert(FetaTlambda && (int)FetaTlambda->size == nmea);
212  assert(etaxi && (int)etaxi->size == npar);
213  assert(etasv && (int)etasv->size == npar);
214  assert(y && (int)y->size == nmea);
215  assert(y_eta && (int)y_eta->size == nmea);
216  assert(Vinvy_eta && (int)Vinvy_eta->size == nmea);
217  assert(FetaV && (int)FetaV->size1 == ncon && (int)FetaV->size2 == nmea);
218  assert(permS && (int)permS->size == ncon);
219  assert(nunm == 0 || (permU && (int)permU->size == nunm));
220  assert(permV && (int)permV->size == nmea);
221 
222  // eta is the part of etaxi containing the measured quantities
223  gsl_vector_view eta = gsl_vector_subvector(etaxi, 0, nmea);
224 
225  // Feta is the part of Fetaxi containing the measured quantities
226 
227  B2DEBUG(11, "==== " << ncon << " " << nmea);
228 
229  gsl_matrix_view Feta = gsl_matrix_submatrix(Fetaxi, 0, 0, ncon, nmea);
230 
231  gsl_matrix_view Vetaeta = gsl_matrix_submatrix(V, 0, 0, nmea, nmea);
232 
233  for (unsigned int i = 0; i < fitobjects.size(); ++i) {
234  for (int ilocal = 0; ilocal < fitobjects[i]->getNPar(); ++ilocal) {
235  if (!fitobjects[i]->isParamFixed(ilocal)) {
236  int iglobal = fitobjects[i]->getGlobalParNum(ilocal);
237  assert(iglobal >= 0 && iglobal < npar);
238  gsl_vector_set(etaxi, iglobal, fitobjects[i]->getParam(ilocal));
239  if (fitobjects[i]->isParamMeasured(ilocal)) {
240  assert(iglobal < nmea);
241  gsl_vector_set(y, iglobal, fitobjects[i]->getMParam(ilocal));
242  }
243  B2DEBUG(10, "etaxi[" << iglobal << "] = " << gsl_vector_get(etaxi, iglobal)
244  << " for jet " << i << " and ilocal = " << ilocal);
245  }
246  }
247  }
248 
250  gsl_matrix_set_zero(Fetaxi);
251  for (int k = 0; k < ncon; k++) {
252  constraints[k]->getDerivatives(Fetaxi->size2, Fetaxi->block->data + k * Fetaxi->tda);
253  if (debug > 1) for (int j = 0; j < npar; j++)
254  if (gsl_matrix_get(Fetaxi, k, j) != 0)
255  B2INFO("1: Fetaxi[" << k << "][" << j << "] = " << gsl_matrix_get(Fetaxi, k, j));
256  }
257 
258  // chi2's, step size, # iterations
259  double chinew = 0, chit = 0, chik = 0;
260  double alph = 1.;
261  nit = 0;
262  // convergence criteria (as in WWINIT)
263  int nitmax = 200;
264  double chik0 = 100.;
265  double chit0 = 100.;
266  double dchikc = 1.0E-3;
267  double dchitc = 1.0E-4;
268  double dchikt = 1.0E-2;
269  double dchik = 1.05;
270  double chimxw = 10000.;
271  double almin = 0.05;
272 
273  // repeat with or with out smaller steps size
274  bool repeat = true;
275  bool scut = false;
276  bool calcerr = true;
277 
278 #ifndef FIT_TRACEOFF
279  if (tracer) tracer->initialize(*this);
280 #endif
281 
282 
283  // start of iterations
284  while (repeat) {
285  bool updatesuccess = true;
286 
287 //*-- If necessary, retry smaller step, same direction
288  if (scut) {
289  gsl_vector_memcpy(etaxi, etasv);
290  updatesuccess = updateFitObjects(etaxi->block->data);
291  if (!updatesuccess) {
292  B2ERROR("OPALFitterGSL::fit: old parameters are garbage!");
293  return -1;
294  }
295 
296 
297  gsl_matrix_set_zero(Fetaxi);
298  for (int k = 0; k < ncon; ++k) {
299  constraints[k]->getDerivatives(Fetaxi->size2, Fetaxi->block->data + k * Fetaxi->tda);
300  }
301  if (debug > 1) debug_print(Fetaxi, "1: Fetaxi");
302  } else {
303  gsl_vector_memcpy(etasv, etaxi);
304  chik0 = chik;
305  chit0 = chit;
306  }
307 
308  // Get covariance matrix
309  gsl_matrix_set_zero(V);
310 
311  for (auto& fitobject : fitobjects) {
312  fitobject->addToGlobCov(V->block->data, V->tda);
313  }
314  if (debug > 1) debug_print(V, "V");
315 
316  gsl_matrix_memcpy(VLU, &Vetaeta.matrix);
317 
318  // invert covariance matrix (needed for chi2 calculation later)
319 
320  int signum;
321  int result;
322  result = gsl_linalg_LU_decomp(VLU, permV, &signum);
323  B2DEBUG(11, "gsl_linalg_LU_decomp result=" << result);
324  if (debug > 3) debug_print(VLU, "VLU");
325 
326  result = gsl_linalg_LU_invert(VLU, permV, Vinv);
327  B2DEBUG(11, "gsl_linalg_LU_invert result=" << result);
328 
329  if (debug > 2) debug_print(Vinv, "Vinv");
330 
331 // *-- Evaluate f and S.
332  for (int k = 0; k < ncon; ++k) {
333  gsl_vector_set(f, k, constraints[k]->getValue());
334  }
335  if (debug > 1) debug_print(f, "f");
336 
337  // y_eta = y - eta
338  gsl_vector_memcpy(y_eta, y);
339  gsl_vector_sub(y_eta, &eta.vector);
340  // r=f
341  gsl_vector_memcpy(r, f);
342  // r = 1*Feta*y_eta + 1*r
343  gsl_blas_dgemv(CblasNoTrans, 1, &Feta.matrix, y_eta, 1, r);
344 
345  if (debug > 1) debug_print(r, "r");
346 
347  // S = Feta * V * Feta^T
348 
349  //FetaV = 1*Feta*V + 0*FetaV
350  B2DEBUG(12, "Creating FetaV");
351  gsl_blas_dsymm(CblasRight, CblasUpper, 1, &Vetaeta.matrix, &Feta.matrix, 0, FetaV);
352  // S = 1 * FetaV * Feta^T + 0*S
353  B2DEBUG(12, "Creating S");
354  gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1, FetaV, &Feta.matrix, 0, S);
355 
356  if (nunm > 0) {
357  // New invention by B. List, 6.12.04:
358  // add F_xi * F_xi^T to S, to make method work when some
359  // constraints do not depend on any measured parameter
360 
361  // Fxi is the part of Fetaxi containing the unmeasured quantities, if any
362  gsl_matrix_view Fxi = gsl_matrix_submatrix(Fetaxi, 0, nmea, ncon, nunm);
363 
364  //S = 1*Fxi*Fxi^T + 1*S
365  gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1, &Fxi.matrix, &Fxi.matrix, 1, S);
366  }
367 
368  if (debug > 1) debug_print(S, "S");
369 
370 // *-- Invert S to Sinv; S is destroyed here!
371 // S is symmetric and positive definite
372 
373  gsl_linalg_LU_decomp(S, permS, &signum);
374  int inverr = gsl_linalg_LU_invert(S, permS, Sinv);
375 
376  if (inverr != 0) {
377  B2ERROR("S: gsl_linalg_LU_invert error " << inverr);
378  ierr = 7;
379  calcerr = false;
380  break;
381  }
382 
383  // Calculate S^1*r here, we will need it
384  // Store it in lambda!
385  // lambda = 1*Sinv*r + 0*lambda; Sinv is symmetric
386  gsl_blas_dsymv(CblasUpper, 1, Sinv, r, 0, lambda);
387 
388 // *-- Calculate new unmeasured quantities, if any
389 
390  if (nunm > 0) {
391  // Fxi is the part of Fetaxi containing the unmeasured quantities, if any
392  gsl_matrix_view Fxi = gsl_matrix_submatrix(Fetaxi, 0, nmea, ncon, nunm);
393  // W1 = Fxi^T * Sinv * Fxi
394  // SinvFxi = 1*Sinv*Fxi + 0*SinvFxi
395  gsl_blas_dsymm(CblasLeft, CblasUpper, 1, Sinv, &Fxi.matrix, 0, SinvFxi);
396  // W1 = 1*Fxi^T*SinvFxi + 0*W1
397  gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1, &Fxi.matrix, SinvFxi, 0, W1);
398 
399  if (debug > 1) {
400  debug_print(W1, "W1");
401  // Check symmetry of W1
402  for (int i = 0; i < nunm; ++i) {
403  for (int j = 0; j < nunm; ++j) {
404  if (abs(gsl_matrix_get(W1, i, j) - gsl_matrix_get(W1, j, i)) > 1E-3 * abs(gsl_matrix_get(W1, i, j) + gsl_matrix_get(W1, j, i)))
405  B2INFO("W1[" << i << "][" << j << "] = " << gsl_matrix_get(W1, i, j)
406  << " W1[" << j << "][" << i << "] = " << gsl_matrix_get(W1, j, i)
407  << " => diff=" << abs(gsl_matrix_get(W1, i, j) - gsl_matrix_get(W1, j, i))
408  << " => tol=" << 1E-3 * abs(gsl_matrix_get(W1, i, j) + gsl_matrix_get(W1, j, i)));
409  }
410  }
411  }
412 
413  // calculate shift of unmeasured parameters
414 
415  // dxi = -alph*W1^-1 * Fxi^T*Sinv*r
416  // => dxi is solution of W1*dxi = -alph*Fxi^T*Sinv*r
417  // Compute rights hand side first
418  // Sinv*r was already calculated and is stored in lambda
419  // dxi = -alph*Fxi^T*lambda + 0*dxi
420 
421  B2DEBUG(11, "alph = " << alph);
422  if (debug > 1) {
423  debug_print(lambda, "lambda");
424  debug_print(&(Fxi.matrix), "Fxi");
425  }
426 
427  gsl_blas_dgemv(CblasTrans, -alph, &Fxi.matrix, lambda, 0, dxi);
428 
429  if (debug > 1) {
430  debug_print(dxi, "dxi0");
431  debug_print(W1, "W1");
432  }
433 
434  // now solve the system
435  // Note added 23.12.04: W1 is symmetric and positive definite,
436  // so we can use the Cholesky instead of LU decomposition
437  gsl_linalg_cholesky_decomp(W1);
438  inverr = gsl_linalg_cholesky_svx(W1, dxi);
439 
440  if (debug > 1) debug_print(dxi, "dxi1");
441 
442 
443  if (inverr != 0) {
444  B2ERROR("W1: gsl_linalg_cholesky_svx error " << inverr);
445  ierr = 8;
446  calcerr = false;
447  break;
448  }
449 
450  if (debug > 1) debug_print(dxi, "dxi");
451 
452 // *-- And now update unmeasured parameters; xi is a view of etaxi
453  // xi is the part of etaxi containing the unmeasured quantities, if any
454  gsl_vector_view xi = gsl_vector_subvector(etaxi, nmea, nunm);
455 
456  // xi += dxi
457  gsl_vector_add(&xi.vector, dxi);
458 
459  }
460 
461 // *-- Calculate new Lagrange multipliers.
462  // lambda = Sinv*r + Sinv*Fxi*dxi
463  // lambda is already set to Sinv*r, we just need to add Sinv*Fxi*dxi
464 
465  if (nunm > 0) {
466  // Fxi is the part of Fetaxi containing the unmeasured quantities, if any
467  gsl_matrix_view Fxi = gsl_matrix_submatrix(Fetaxi, 0, nmea, ncon, nunm);
468  // calculate Fxidxi = 1*Fxi*dxi + 0*Fxidxi
469  gsl_blas_dgemv(CblasNoTrans, 1, &Fxi.matrix, dxi, 0, Fxidxi);
470  // add to existing lambda: lambda = 1*Sinv*Fxidxi + 1*lambda; Sinv is symmetric
471  gsl_blas_dsymv(CblasUpper, 1, Sinv, Fxidxi, 1, lambda);
472 
473  }
474 
475  if (debug > 1) debug_print(lambda, "lambda");
476 
477 // *-- Calculate new fitted parameters:
478  // eta = y - V*Feta^T*lambda
479  gsl_vector_memcpy(&eta.vector, y);
480  // FetaTlambda = 1*Feta^T*lambda + 0*FetaTlambda
481  gsl_blas_dgemv(CblasTrans, 1, &Feta.matrix, lambda, 0, FetaTlambda);
482  // eta = -1*V*FetaTlambda + 1*eta; V is symmetric
483  gsl_blas_dsymv(CblasUpper, -1, &Vetaeta.matrix, FetaTlambda, 1, &eta.vector);
484 
485 
486  if (debug > 1) debug_print(&eta.vector, "updated eta");
487 
488 // *-- Calculate constraints and their derivatives.
489  // since the constraints ask the fitobjects for their parameters,
490  // we need to update the fitobjects first!
491  // COULD BE DONE: update also ERRORS! (now only in the very end!)
492  updatesuccess = updateFitObjects(etaxi->block->data);
493 
494  B2DEBUG(10, "After adjustment of all parameters:\n");
495  for (int k = 0; k < ncon; ++k) {
496  B2DEBUG(10, "Value of constraint " << k << " = " << constraints[k]->getValue());
497  }
498  gsl_matrix_set_zero(Fetaxi);
499  for (int k = 0; k < ncon; k++) {
500  constraints[k]->getDerivatives(Fetaxi->size2, Fetaxi->block->data + k * Fetaxi->tda);
501  }
502  if (debug > 1) debug_print(Fetaxi, "2: Fetaxi");
503 
504 
505 // *-- Calculate new chisq.
506  // First, calculate new y-eta
507  // y_eta = y - eta
508  gsl_vector_memcpy(y_eta, y);
509  gsl_vector_sub(y_eta, &eta.vector);
510  // Now calculate Vinv*y_eta [ as solution to V* Vinvy_eta = y_eta]
511  // Vinvy_eta = 1*Vinv*y_eta + 0*Vinvy_eta; Vinv is symmetric
512  gsl_blas_dsymv(CblasUpper, 1, Vinv, y_eta, 0, Vinvy_eta);
513  // Now calculate y_eta *Vinvy_eta
514  gsl_blas_ddot(y_eta, Vinvy_eta, &chit);
515 
516  for (int i = 0; i < nmea; ++i)
517  for (int j = 0; j < nmea; ++j) {
518  double dchit = (gsl_vector_get(y_eta, i)) *
519  gsl_matrix_get(Vinv, i, j) *
520  (gsl_vector_get(y_eta, j));
521  if (dchit != 0)
522  B2DEBUG(11, "chit for i,j = " << i << " , " << j << " = "
523  << dchit);
524  }
525 
526  chik = 0;
527  for (int k = 0; k < ncon; ++k) chik += std::abs(2 * gsl_vector_get(lambda, k) * constraints[k]->getValue());
528  chinew = chit + chik;
529 
530 //*-- Calculate change in chisq, and check constraints are satisfied.
531  nit++;
532 
533  bool sconv = (std::abs(chik - chik0) < dchikc)
534  && (std::abs(chit - chit0) < dchitc * chit)
535  && (chik < dchikt * chit);
536  // Second convergence criterium:
537  // If all constraints are fulfilled to better than 1E-8,
538  // and all parameters have changed by less than 1E-8,
539  // assume convergence
540  // This criterium assumes that all constraints and all parameters
541  // are "of order 1", i.e. their natural values are around 1 to 100,
542  // as for GeV or radians
543  double eps = 1E-6;
544  bool sconv2 = true;
545  for (int k = 0; sconv2 && (k < ncon); ++k)
546  sconv2 &= (std::abs(gsl_vector_get(f, k)) < eps);
547  if (sconv2)
548  B2DEBUG(10, "All constraints fulfilled to better than " << eps);
549 
550  for (int j = 0; sconv2 && (j < npar); ++j)
551  sconv2 &= (std::abs(gsl_vector_get(etaxi, j) - gsl_vector_get(etasv, j)) < eps);
552  if (sconv2)
553  B2DEBUG(10, "All parameters stable to better than " << eps);
554  sconv |= sconv2;
555 
556  bool sbad = (chik > dchik * chik0)
557  && (chik > dchikt * chit)
558  && (chik > chik0 + 1.E-10);
559 
560  scut = false;
561 
562  if (nit > nitmax) {
563 // *-- Out of iterations
564  repeat = false;
565  calcerr = false;
566  ierr = 1;
567  } else if (sconv && updatesuccess) {
568 // *-- Converged
569  repeat = false;
570  calcerr = true;
571  ierr = 0;
572  } else if (nit > 2 && chinew > chimxw && updatesuccess) {
573 // *-- Chi2 crazy ?
574  repeat = false;
575  calcerr = false;
576  ierr = 2;
577  } else if ((sbad && nit > 1) || !updatesuccess) {
578 // *-- ChiK increased, try smaller step
579  if (alph == almin) {
580  repeat = true; // false;
581  ierr = 3;
582  } else {
583  alph = std::max(almin, 0.5 * alph);
584  scut = true;
585  repeat = true;
586  ierr = 4;
587  }
588  } else {
589 // *-- Keep going..
590  alph = std::min(alph + 0.1, 1.);
591  repeat = true;
592  ierr = 5;
593  }
594 
595  B2DEBUG(10, "======== NIT = " << nit << ", CHI2 = " << chinew
596  << ", ierr = " << ierr << ", alph=" << alph);
597 
598  for (unsigned int i = 0; i < fitobjects.size(); ++i)
599  B2DEBUG(10, "fitobject " << i << ": " << *fitobjects[i]);
600 
601 #ifndef FIT_TRACEOFF
602  if (tracer) tracer->step(*this);
603 #endif
604 
605  } // end of while (repeat)
606 
607 // *-- Turn chisq into probability.
608  fitprob = (ncon - nunm > 0) ? gsl_cdf_chisq_Q(chinew, ncon - nunm) : 0.5;
609  chi2 = chinew;
610 
611 // *-- End of iterations - calculate errors.
612 
613 // The result will (ultimately) be stored in Vnew
614 
615  gsl_matrix_set_zero(Vnew);
616  gsl_matrix_set_zero(Minv);
617  if (debug > 2) {
618  for (int i = 0; i < npar; ++i) {
619  for (int j = 0; j < npar; ++j) {
620  B2INFO("Minv[" << i << "," << j << "]=" << gsl_matrix_get(Minv, i, j));
621  }
622  }
623  }
624 
625  B2DEBUG(10, "OPALFitterGSL: calcerr = " << calcerr);
626 
627  if (calcerr) {
628 
629 // *-- As a first step, calculate Minv as in 9.4.2 of Benno's book chapter
630 // (in O.Behnke et al "Data Analysis in High Energy Physics")
631 
632 
633 // *-- Evaluate S and invert.
634 
635  if (debug > 2) {
636  debug_print(&Vetaeta.matrix, "V");
637  debug_print(&Feta.matrix, "Feta");
638  }
639 
640  // CblasRight means C = alpha B A + beta C with symmetric matrix A
641  //FetaV[ncon][nmea] = 1*Feta[ncon][nmea]*V[nmea][nmea] + 0*FetaV
642  gsl_blas_dsymm(CblasRight, CblasUpper, 1, &Vetaeta.matrix, &Feta.matrix, 0, FetaV);
643  // S[ncon][ncon] = 1 * FetaV[ncon][nmea] * Feta^T[nmea][ncon] + 0*S
644  gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1, FetaV, &Feta.matrix, 0, S);
645 
646 
647  if (nunm > 0) {
648  // New invention by B. List, 6.12.04:
649  // add F_xi * F_xi^T to S, to make method work when some
650  // constraints do not depend on any measured parameter
651 
652  // Fxi is the part of Fetaxi containing the unmeasured quantities, if any
653  gsl_matrix_view Fxi = gsl_matrix_submatrix(Fetaxi, 0, nmea, ncon, nunm);
654  //S[ncon][ncon] = 1*Fxi[ncon][nunm]*Fxi^T[nunm][ncon] + 1*S[ncon][ncon]
655  gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1, &Fxi.matrix, &Fxi.matrix, 1, S);
656  }
657 
658  if (debug > 2) debug_print(S, "S");
659 
660 // *-- Invert S, testing for singularity first.
661 // S is symmetric and positive definite
662 
663  int signum;
664  gsl_linalg_LU_decomp(S, permS, &signum);
665  int inverr = gsl_linalg_LU_invert(S, permS, Sinv);
666 
667  if (inverr != 0) {
668  B2ERROR("S: gsl_linalg_LU_invert error " << inverr << " in error calculation");
669  ierr = -1;
670  return -1;
671  }
672 
673 
674 // *-- Calculate G.
675 // (same as W1, but for measured parameters)
676 // G = Feta^T * Sinv * Feta
677 
678  // SinvFeta[ncon][nmea] = 1*Sinv[ncon][ncon]*Feta[ncon][nmea] + 0*SinvFeta
679  gsl_blas_dsymm(CblasLeft, CblasUpper, 1, Sinv, &Feta.matrix, 0, SinvFeta);
680  // G[nmea][nmea] = 1*Feta^T[nmea][ncon]*SinvFeta[ncon][nmea] + 0*G
681  gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1, &Feta.matrix, SinvFeta, 0, G);
682 
683  if (debug > 2) debug_print(G, "G(1)");
684 
685 
686  if (nunm > 0) {
687 
688 // *-- Calculate H
689 
690  // Fxi is the part of Fetaxi containing the unmeasured quantities, if any
691  gsl_matrix_view Fxi = gsl_matrix_submatrix(Fetaxi, 0, nmea, ncon, nunm);
692  // H = Feta^T * Sinv * Fxi
693  // SinvFxi[ncon][nunm] = 1*Sinv[ncon][ncon]*Fxi[ncon][nunm] + 0*SinvFxi
694  gsl_blas_dsymm(CblasLeft, CblasUpper, 1, Sinv, &Fxi.matrix, 0, SinvFxi);
695  // H[nmea][nunm] = 1*Feta^T[nmea][ncon]*SinvFxi[ncon][nunm] + 0*H
696  gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1, &Feta.matrix, SinvFxi, 0, H);
697 
698  if (debug > 2) debug_print(H, "H");
699 
700 // *-- Calculate U**-1 and invert.
701 // (same as W1)
702 // U is a part of Minv
703 
704  gsl_matrix* Uinv = W1;
705  gsl_matrix_view U = gsl_matrix_submatrix(Minv, nmea, nmea, nunm, nunm);
706  // Uinv = Fxi^T * Sinv * Fxi
707  // Uinv[nunm][nunm] = 1*Fxi^T[nunm][ncon]*SinvFxi[ncon][nunm] + 0*W1
708  gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1, &Fxi.matrix, SinvFxi, 0, Uinv);
709 
710  gsl_linalg_LU_decomp(Uinv, permU, &signum);
711  inverr = gsl_linalg_LU_invert(Uinv, permU, &U.matrix);
712 
713  if (debug > 2) debug_print(&U.matrix, "U");
714  for (int i = 0; i < npar; ++i) {
715  for (int j = 0; j < npar; ++j) {
716  B2DEBUG(12, "after U Minv[" << i << "," << j << "]=" << gsl_matrix_get(Minv, i, j));
717  }
718  }
719 
720  if (inverr != 0) {
721  B2ERROR("U: gsl_linalg_LU_invert error " << inverr << " in error calculation ");
722  ierr = -1;
723  return -1;
724  }
725 
726 // *-- Covariance matrix between measured and unmeasured parameters.
727 
728 // HU[nmea][nunm] = 1*H[nmea][nunm]*U[nunm][nunm] + 0*HU
729  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, H, &U.matrix, 0, HU);
730 // Vnewetaxi is a view of Vnew
731  gsl_matrix_view Minvetaxi = gsl_matrix_submatrix(Minv, 0, nmea, nmea, nunm);
732  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, -1, &Vetaeta.matrix, HU, 0, &Minvetaxi.matrix);
733  for (int i = 0; i < npar; ++i) {
734  for (int j = 0; j < npar; ++j) {
735  B2DEBUG(12, "after etaxi Minv[" << i << "," << j << "]=" << gsl_matrix_get(Minv, i, j));
736  }
737  }
738 
739 
740 // *-- Fill in symmetric part:
741 // Vnewxieta is a view of Vnew
742  gsl_matrix_view Minvxieta = gsl_matrix_submatrix(Minv, nmea, 0, nunm, nmea);
743  gsl_matrix_transpose_memcpy(&Minvxieta.matrix, &Minvetaxi.matrix);
744  for (int i = 0; i < npar; ++i) {
745  for (int j = 0; j < npar; ++j) {
746  B2DEBUG(12, "after symmetric: Minv[" << i << "," << j << "]=" << gsl_matrix_get(Minv, i, j));
747  }
748  }
749 
750 // *-- Calculate G-HUH^T:
751 // G = -1*HU*H^T +1*G
752  gsl_blas_dgemm(CblasNoTrans, CblasTrans, -1, HU, H, +1, G);
753 
754  } // endif nunm > 0
755 
756 // *-- Calculate I-GV.
757  // IGV = 1
758  gsl_matrix_set_identity(IGV);
759  // IGV = -1*G*Vetaeta + 1*IGV
760  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, -1, G, &Vetaeta.matrix, 1, IGV);
761 
762 // *-- And finally error matrix on fitted parameters.
763  gsl_matrix_view Minvetaeta = gsl_matrix_submatrix(Minv, 0, 0, nmea, nmea);
764 
765  // Vnewetaeta = 1*Vetaeta*IGV + 0*Vnewetaeta
766  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, &Vetaeta.matrix, IGV, 0, &Minvetaeta.matrix);
767 
768  for (int i = 0; i < npar; ++i) {
769  for (int j = 0; j < npar; ++j) {
770  B2DEBUG(12, "complete Minv[" << i << "," << j << "]=" << gsl_matrix_get(Minv, i, j));
771  }
772  }
773 
774 // *-- now Minv should be complete, can calculate new covariance matrix Vnew
775 // Vnew = (dx / dt)^T * V_t * dx / dt
776 // x are the fitted parameters, t are the measured parameters,
777 // V_t is the covariance matrix of the measured parameters
778 // thus in OPALFitter variables: V_t = Vetaeta, x = (eta, xi) = etaxi
779 // d eta / dt and d xi / dt are given by eqns 9.54 and 9.55, respectively,
780 // as deta/dt = - Minvetaeta * Fetat and dxi/dt = - Minvxieta * Fetat
781 // Fetat is d^2 chi^2 / deta dt = - V^-1, even in presence of soft constraints, since
782 // the soft constraints don't depend on the _measured_ parameters t (but on the fitted ones)
783 // HERE, V (and Vinv) do not include the second derivatives of the soft constraints
784 // so we can use them right away
785 // Note that "F" is used for completely different things:
786 // in OPALFitter it is the derivatives of the constraints wrt to the parameters (A in book)
787 // in the book, F are the second derivatives of the objective function wrt the parameters
788 
789 
790  gsl_matrix_view detadt = gsl_matrix_submatrix(dxdt, 0, 0, nmea, nmea);
791  B2DEBUG(13, "after detadt");
792  // Vdxdt is Vetaeta * dxdt^T, thus Vdxdt[nmea][npar]
793  gsl_matrix_view Vdetadt = gsl_matrix_submatrix(Vdxdt, 0, 0, nmea, nmea);
794  B2DEBUG(13, "after Vdetadt");
795 
796  // detadt = - Minvetaeta * Fetat = -1 * Minvetaeta * (-1) * Vinv + 0 * detadt // replace by symm?
797  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, &Minvetaeta.matrix, Vinv, 0, &detadt.matrix);
798  if (debug > 2) debug_print(&detadt.matrix, "deta/dt");
799 
800  // Vdetadt = 1 * Vetaeta * detadt^T + 0* Vdetadt
801  gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1, &Vetaeta.matrix, &detadt.matrix, 0, &Vdetadt.matrix); // ok
802  if (debug > 2) debug_print(&Vdetadt.matrix, "Vetata * deta/dt");
803 
804  gsl_matrix_view Vnewetaeta = gsl_matrix_submatrix(Vnew, 0, 0, nmea, nmea); //[nmea],[nmea]
805  // Vnewetaeta = 1 * detadt * Vdetadt + 0* Vnewetaeta
806  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, &detadt.matrix, &Vdetadt.matrix, 0, &Vnewetaeta.matrix);
807 
808  if (debug > 2) debug_print(Vnew, "Vnew after part for measured parameters");
809 
810 
811  if (nunm > 0) {
812 
813  gsl_matrix_view Minvxieta = gsl_matrix_submatrix(Minv, nmea, 0, nunm, nmea); //[nunm][nmea]
814  B2DEBUG(13, "after Minvxieta");
815  if (debug > 2) debug_print(&Minvxieta.matrix, "Minvxieta");
816 
817  gsl_matrix_view dxidt = gsl_matrix_submatrix(dxdt, nmea, 0, nunm, nmea); //[nunm][nmea]
818  B2DEBUG(13, "after dxidt");
819  // dxidt[nunm][nmea] = - Minvxieta * Fetat = -1 * Minvxieta[nunm][nmea] * Vinv[nmea][nmea] + 0 * dxidt
820  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, &Minvxieta.matrix, Vinv, 0, &dxidt.matrix); //ok
821  if (debug > 2) debug_print(&dxidt.matrix, "dxi/dt");
822 
823  // Vdxdt = V * dxdt^T => Vdxdt[nmea][npar]
824  gsl_matrix_view Vdxidt = gsl_matrix_submatrix(Vdxdt, 0, nmea, nmea, nunm); //[nmea][nunm]
825  B2DEBUG(13, "after Vdxidt");
826  // Vdxidt = 1 * Vetaeta[nmea][nmea] * dxidt^T[nmea][nunm] + 0* Vdxidt => Vdxidt[nmea][nunm]
827  gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1, &Vetaeta.matrix, &dxidt.matrix, 0, &Vdxidt.matrix); // ok
828  if (debug > 2) debug_print(&Vdxidt.matrix, "Vetaeta * dxi/dt^T");
829 
830  gsl_matrix_view Vnewetaxi = gsl_matrix_submatrix(Vnew, 0, nmea, nmea, nunm); //[nmea][nunm]
831  gsl_matrix_view Vnewxieta = gsl_matrix_submatrix(Vnew, nmea, 0, nunm, nmea); //[nunm][nmea]
832  gsl_matrix_view Vnewxixi = gsl_matrix_submatrix(Vnew, nmea, nmea, nunm, nunm); //[nunm][nunm]
833 
834  // Vnewxieta[nunm][nmea] = 1 * dxidt[nunm][nmea] * Vdetadt[nmea][nmea] + 0* Vnewxieta
835  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, &dxidt.matrix, &Vdetadt.matrix, 0, &Vnewxieta.matrix); // ok
836  if (debug > 2) debug_print(Vnew, "Vnew after xieta part");
837  // Vnewetaxi[nmea][nunm] = 1 * detadt[nmea][nmea] * Vdxidt[nmea][nunm] + 0* Vnewetaxi
838  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, &detadt.matrix, &Vdxidt.matrix, 0, &Vnewetaxi.matrix); // ok
839  if (debug > 2) debug_print(Vnew, "Vnew after etaxi part");
840  // Vnewxixi[nunm][nunm] = 1 * dxidt[nunm][nmea] * Vdxidt[nmea][nunm] + 0* Vnewxixi
841  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, &dxidt.matrix, &Vdxidt.matrix, 0, &Vnewxixi.matrix);
842  if (debug > 2) debug_print(Vnew, "Vnew after xixi part");
843  }
844 
845 // *-- now we finally have Vnew, fill into fitobjects
846 
847  // update errors in fitobjects
848  for (auto& fitobject : fitobjects) {
849  for (int ilocal = 0; ilocal < fitobject->getNPar(); ++ilocal) {
850  int iglobal = fitobject->getGlobalParNum(ilocal);
851  for (int jlocal = ilocal; jlocal < fitobject->getNPar(); ++jlocal) {
852  int jglobal = fitobject->getGlobalParNum(jlocal);
853  if (iglobal >= 0 && jglobal >= 0)
854  fitobject->setCov(ilocal, jlocal, gsl_matrix_get(Vnew, iglobal, jglobal));
855  }
856  }
857  }
858 
859  // Finally, copy covariance matrix
860  if (cov && covDim != nmea + nunm) {
861  delete[] cov;
862  cov = nullptr;
863  }
864  covDim = nmea + nunm;
865  if (!cov) cov = new double[covDim * covDim];
866  for (int i = 0; i < covDim; ++i) {
867  for (int j = 0; j < covDim; ++j) {
868  cov[i * covDim + j] = gsl_matrix_get(Vnew, i, j);
869  }
870  }
871  covValid = true;
872 
873  } // endif calcerr == true
874 
875 #ifndef FIT_TRACEOFF
876  if (tracer) tracer->finish(*this);
877 #endif
878 
879  return fitprob;
880 
881  }
882 
883  bool OPALFitterGSL::initialize()
884  {
885  covValid = false;
886  // tell fitobjects the global ordering of their parameters:
887  int iglobal = 0;
888  // measured parameters first!
889  for (auto& fitobject : fitobjects) {
890  for (int ilocal = 0; ilocal < fitobject->getNPar(); ++ilocal) {
891  if (fitobject->isParamMeasured(ilocal) &&
892  !fitobject->isParamFixed(ilocal)) {
893  fitobject->setGlobalParNum(ilocal, iglobal);
894  B2DEBUG(10, "Object " << fitobject->getName()
895  << " Parameter " << fitobject->getParamName(ilocal)
896  << " is measured, global number " << iglobal);
897  ++iglobal;
898  }
899  }
900  }
901  nmea = iglobal;
902  // now unmeasured parameters!
903  for (auto& fitobject : fitobjects) {
904  for (int ilocal = 0; ilocal < fitobject->getNPar(); ++ilocal) {
905  if (!fitobject->isParamMeasured(ilocal) &&
906  !fitobject->isParamFixed(ilocal)) {
907  fitobject->setGlobalParNum(ilocal, iglobal);
908  B2DEBUG(10, "Object " << fitobject->getName()
909  << " Parameter " << fitobject->getParamName(ilocal)
910  << " is unmeasured, global number " << iglobal);
911  ++iglobal;
912  }
913  }
914  }
915  npar = iglobal;
916  assert(npar <= NPARMAX);
917  nunm = npar - nmea;
918  assert(nunm <= NUNMMAX);
919 
920  // set number of constraints
921  ncon = constraints.size();
922  assert(ncon <= NCONMAX);
923 
924  ini_gsl_vector(f, ncon);
925  ini_gsl_vector(r, ncon);
926 
927  ini_gsl_matrix(Fetaxi, ncon, npar);
928  ini_gsl_matrix(S, ncon, ncon);
929  ini_gsl_matrix(Sinv, ncon, ncon);
930  ini_gsl_matrix(SinvFxi, ncon, nunm);
931  ini_gsl_matrix(SinvFeta, ncon, nmea);
932  ini_gsl_matrix(W1, nunm, nunm);
933  ini_gsl_matrix(G, nmea, nmea);
934  ini_gsl_matrix(H, nmea, nunm);
935  ini_gsl_matrix(HU, nmea, nunm);
936  ini_gsl_matrix(IGV, nmea, nmea);
937  ini_gsl_matrix(V, npar, npar);
938  ini_gsl_matrix(VLU, nmea, nmea);
939  ini_gsl_matrix(Vinv, nmea, nmea);
940  ini_gsl_matrix(Vnew, npar, npar);
941  ini_gsl_matrix(Minv, npar, npar);
942  ini_gsl_matrix(dxdt, npar, nmea);
943  ini_gsl_matrix(Vdxdt, nmea, npar);
944 
945  ini_gsl_vector(dxi, nunm);
946  ini_gsl_vector(Fxidxi, ncon);
947  ini_gsl_vector(lambda, ncon);
948  ini_gsl_vector(FetaTlambda, nmea);
949 
950  ini_gsl_vector(etaxi, npar);
951  ini_gsl_vector(etasv, npar);
952  ini_gsl_vector(y, nmea);
953  ini_gsl_vector(y_eta, nmea);
954  ini_gsl_vector(Vinvy_eta, nmea);
955 
956  ini_gsl_matrix(FetaV, ncon, nmea);
957 
958  ini_gsl_permutation(permS, ncon);
959  ini_gsl_permutation(permU, nunm);
960  ini_gsl_permutation(permV, nmea);
961 
962  assert(f && (int)f->size == ncon);
963  assert(r && (int)r->size == ncon);
964  assert(Fetaxi && (int)Fetaxi->size1 == ncon && (int)Fetaxi->size2 == npar);
965  assert(S && (int)S->size1 == ncon && (int)S->size2 == ncon);
966  assert(Sinv && (int)Sinv->size1 == ncon && (int)Sinv->size2 == ncon);
967  assert(nunm == 0 || (SinvFxi && (int)SinvFxi->size1 == ncon && (int)SinvFxi->size2 == nunm));
968  assert(SinvFeta && (int)SinvFeta->size1 == ncon && (int)SinvFeta->size2 == nmea);
969  assert(nunm == 0 || (W1 && (int)W1->size1 == nunm && (int)W1->size2 == nunm));
970  assert(G && (int)G->size1 == nmea && (int)G->size2 == nmea);
971  assert(nunm == 0 || (H && (int)H->size1 == nmea && (int)H->size2 == nunm));
972  assert(nunm == 0 || (HU && (int)HU->size1 == nmea && (int)HU->size2 == nunm));
973  assert(IGV && (int)IGV->size1 == nmea && (int)IGV->size2 == nmea);
974  assert(V && (int)V->size1 == npar && (int)V->size2 == npar);
975  assert(VLU && (int)VLU->size1 == nmea && (int)VLU->size2 == nmea);
976  assert(Vinv && (int)Vinv->size1 == nmea && (int)Vinv->size2 == nmea);
977  assert(Vnew && (int)Vnew->size1 == npar && (int)Vnew->size2 == npar);
978  assert(nunm == 0 || (dxi && (int)dxi->size == nunm));
979  assert(nunm == 0 || (Fxidxi && (int)Fxidxi->size == ncon));
980  assert(lambda && (int)lambda->size == ncon);
981  assert(FetaTlambda && (int)FetaTlambda->size == nmea);
982  assert(etaxi && (int)etaxi->size == npar);
983  assert(etasv && (int)etasv->size == npar);
984  assert(y && (int)y->size == nmea);
985  assert(y_eta && (int)y_eta->size == nmea);
986  assert(Vinvy_eta && (int)Vinvy_eta->size == nmea);
987  assert(FetaV && (int)FetaV->size1 == ncon && (int)FetaV->size2 == nmea);
988  assert(permS && (int)permS->size == ncon);
989  assert(permS && (int)permS->size == ncon);
990  assert(nunm == 0 || (permU && (int)permU->size == nunm));
991  assert(permV && (int)permV->size == nmea);
992 
993  return true;
994 
995  }
996 
997  bool OPALFitterGSL::updateFitObjects(double eetaxi[])
998  {
999  // changed etaxi -> eetaxi to avoid clash with class member DJeans
1000  //bool debug = false;
1001  bool result = true;
1002  for (auto& fitobject : fitobjects) {
1003  for (int ilocal = 0; ilocal < fitobject->getNPar(); ++ilocal) {
1004  fitobject->updateParams(eetaxi, npar);
1005 // }
1006  }
1007  }
1008  return result;
1009  }
1010 
1011  int OPALFitterGSL::getError() const {return ierr;}
1012  double OPALFitterGSL::getProbability() const {return fitprob;}
1013  double OPALFitterGSL::getChi2() const {return chi2;}
1014  int OPALFitterGSL::getDoF() const {return ncon - nunm;}
1015  int OPALFitterGSL::getIterations() const {return nit;}
1016 
1017  void OPALFitterGSL::ini_gsl_permutation(gsl_permutation*& p, unsigned int size)
1018  {
1019  if (p) {
1020  if (p->size != size) {
1021  gsl_permutation_free(p);
1022  if (size > 0) p = gsl_permutation_alloc(size);
1023  else p = nullptr;
1024  }
1025  } else if (size > 0) p = gsl_permutation_alloc(size);
1026  }
1027 
1028  void OPALFitterGSL::ini_gsl_vector(gsl_vector*& v, unsigned int size)
1029  {
1030 
1031  if (v) {
1032  if (v->size != size) {
1033  gsl_vector_free(v);
1034  if (size > 0) v = gsl_vector_alloc(size);
1035  else v = nullptr;
1036  }
1037  } else if (size > 0) v = gsl_vector_alloc(size);
1038  }
1039 
1040  void OPALFitterGSL::ini_gsl_matrix(gsl_matrix*& m, unsigned int size1, unsigned int size2)
1041  {
1042  if (m) {
1043  if (m->size1 != size1 || m->size2 != size2) {
1044  gsl_matrix_free(m);
1045  if (size1 * size2 > 0) m = gsl_matrix_alloc(size1, size2);
1046  else m = nullptr;
1047  }
1048  } else if (size1 * size2 > 0) m = gsl_matrix_alloc(size1, size2);
1049  }
1050 
1051  void OPALFitterGSL::setDebug(int debuglevel)
1052  {
1053  debug = debuglevel;
1054  }
1055 
1056  void OPALFitterGSL::debug_print(gsl_matrix* m, const char* name)
1057  {
1058  for (unsigned int i = 0; i < m->size1; ++i)
1059  for (unsigned int j = 0; j < m->size2; ++j)
1060  if (gsl_matrix_get(m, i, j) != 0)
1061  B2INFO(name << "[" << i << "][" << j << "]=" << gsl_matrix_get(m, i, j));
1062  }
1063 
1064  void OPALFitterGSL::debug_print(gsl_vector* v, const char* name)
1065  {
1066  for (unsigned int i = 0; i < v->size; ++i)
1067  if (gsl_vector_get(v, i) != 0)
1068  B2INFO(name << "[" << i << "]=" << gsl_vector_get(v, i));
1069  }
1070 
1071  int OPALFitterGSL::getNcon() const {return ncon;}
1072  int OPALFitterGSL::getNsoft() const {return 0;}
1073  int OPALFitterGSL::getNunm() const {return nunm;}
1074  int OPALFitterGSL::getNpar() const {return npar;}
1075 
1076  }// end OrcaKinFit namespace
1078 } // end Belle2 namespace
1079 
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17