Belle II Software  release-08-01-10
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  // cppcheck-suppress duplicateCondition
466  if (nunm > 0) {
467  // Fxi is the part of Fetaxi containing the unmeasured quantities, if any
468  gsl_matrix_view Fxi = gsl_matrix_submatrix(Fetaxi, 0, nmea, ncon, nunm);
469  // calculate Fxidxi = 1*Fxi*dxi + 0*Fxidxi
470  gsl_blas_dgemv(CblasNoTrans, 1, &Fxi.matrix, dxi, 0, Fxidxi);
471  // add to existing lambda: lambda = 1*Sinv*Fxidxi + 1*lambda; Sinv is symmetric
472  gsl_blas_dsymv(CblasUpper, 1, Sinv, Fxidxi, 1, lambda);
473 
474  }
475 
476  if (debug > 1) debug_print(lambda, "lambda");
477 
478 // *-- Calculate new fitted parameters:
479  // eta = y - V*Feta^T*lambda
480  gsl_vector_memcpy(&eta.vector, y);
481  // FetaTlambda = 1*Feta^T*lambda + 0*FetaTlambda
482  gsl_blas_dgemv(CblasTrans, 1, &Feta.matrix, lambda, 0, FetaTlambda);
483  // eta = -1*V*FetaTlambda + 1*eta; V is symmetric
484  gsl_blas_dsymv(CblasUpper, -1, &Vetaeta.matrix, FetaTlambda, 1, &eta.vector);
485 
486 
487  if (debug > 1) debug_print(&eta.vector, "updated eta");
488 
489 // *-- Calculate constraints and their derivatives.
490  // since the constraints ask the fitobjects for their parameters,
491  // we need to update the fitobjects first!
492  // COULD BE DONE: update also ERRORS! (now only in the very end!)
493  updatesuccess = updateFitObjects(etaxi->block->data);
494 
495  B2DEBUG(10, "After adjustment of all parameters:\n");
496  for (int k = 0; k < ncon; ++k) {
497  B2DEBUG(10, "Value of constraint " << k << " = " << constraints[k]->getValue());
498  }
499  gsl_matrix_set_zero(Fetaxi);
500  for (int k = 0; k < ncon; k++) {
501  constraints[k]->getDerivatives(Fetaxi->size2, Fetaxi->block->data + k * Fetaxi->tda);
502  }
503  if (debug > 1) debug_print(Fetaxi, "2: Fetaxi");
504 
505 
506 // *-- Calculate new chisq.
507  // First, calculate new y-eta
508  // y_eta = y - eta
509  gsl_vector_memcpy(y_eta, y);
510  gsl_vector_sub(y_eta, &eta.vector);
511  // Now calculate Vinv*y_eta [ as solution to V* Vinvy_eta = y_eta]
512  // Vinvy_eta = 1*Vinv*y_eta + 0*Vinvy_eta; Vinv is symmetric
513  gsl_blas_dsymv(CblasUpper, 1, Vinv, y_eta, 0, Vinvy_eta);
514  // Now calculate y_eta *Vinvy_eta
515  gsl_blas_ddot(y_eta, Vinvy_eta, &chit);
516 
517  for (int i = 0; i < nmea; ++i)
518  for (int j = 0; j < nmea; ++j) {
519  double dchit = (gsl_vector_get(y_eta, i)) *
520  gsl_matrix_get(Vinv, i, j) *
521  (gsl_vector_get(y_eta, j));
522  if (dchit != 0)
523  B2DEBUG(11, "chit for i,j = " << i << " , " << j << " = "
524  << dchit);
525  }
526 
527  chik = 0;
528  for (int k = 0; k < ncon; ++k) chik += std::abs(2 * gsl_vector_get(lambda, k) * constraints[k]->getValue());
529  chinew = chit + chik;
530 
531 //*-- Calculate change in chisq, and check constraints are satisfied.
532  nit++;
533 
534  bool sconv = (std::abs(chik - chik0) < dchikc)
535  && (std::abs(chit - chit0) < dchitc * chit)
536  && (chik < dchikt * chit);
537  // Second convergence criterium:
538  // If all constraints are fulfilled to better than 1E-8,
539  // and all parameters have changed by less than 1E-8,
540  // assume convergence
541  // This criterium assumes that all constraints and all parameters
542  // are "of order 1", i.e. their natural values are around 1 to 100,
543  // as for GeV or radians
544  double eps = 1E-6;
545  bool sconv2 = true;
546  for (int k = 0; sconv2 && (k < ncon); ++k)
547  sconv2 &= (std::abs(gsl_vector_get(f, k)) < eps);
548  if (sconv2)
549  B2DEBUG(10, "All constraints fulfilled to better than " << eps);
550 
551  for (int j = 0; sconv2 && (j < npar); ++j)
552  sconv2 &= (std::abs(gsl_vector_get(etaxi, j) - gsl_vector_get(etasv, j)) < eps);
553  if (sconv2)
554  B2DEBUG(10, "All parameters stable to better than " << eps);
555  sconv |= sconv2;
556 
557  bool sbad = (chik > dchik * chik0)
558  && (chik > dchikt * chit)
559  && (chik > chik0 + 1.E-10);
560 
561  scut = false;
562 
563  if (nit > nitmax) {
564 // *-- Out of iterations
565  repeat = false;
566  calcerr = false;
567  ierr = 1;
568  } else if (sconv && updatesuccess) {
569 // *-- Converged
570  repeat = false;
571  calcerr = true;
572  ierr = 0;
573  } else if (nit > 2 && chinew > chimxw && updatesuccess) {
574 // *-- Chi2 crazy ?
575  repeat = false;
576  calcerr = false;
577  ierr = 2;
578  } else if ((sbad && nit > 1) || !updatesuccess) {
579 // *-- ChiK increased, try smaller step
580  if (alph == almin) {
581  repeat = true; // false;
582  ierr = 3;
583  } else {
584  alph = std::max(almin, 0.5 * alph);
585  scut = true;
586  repeat = true;
587  ierr = 4;
588  }
589  } else {
590 // *-- Keep going..
591  alph = std::min(alph + 0.1, 1.);
592  repeat = true;
593  ierr = 5;
594  }
595 
596  B2DEBUG(10, "======== NIT = " << nit << ", CHI2 = " << chinew
597  << ", ierr = " << ierr << ", alph=" << alph);
598 
599  for (unsigned int i = 0; i < fitobjects.size(); ++i)
600  B2DEBUG(10, "fitobject " << i << ": " << *fitobjects[i]);
601 
602 #ifndef FIT_TRACEOFF
603  if (tracer) tracer->step(*this);
604 #endif
605 
606  } // end of while (repeat)
607 
608 // *-- Turn chisq into probability.
609  fitprob = (ncon - nunm > 0) ? gsl_cdf_chisq_Q(chinew, ncon - nunm) : 0.5;
610  chi2 = chinew;
611 
612 // *-- End of iterations - calculate errors.
613 
614 // The result will (ultimately) be stored in Vnew
615 
616  gsl_matrix_set_zero(Vnew);
617  gsl_matrix_set_zero(Minv);
618  if (debug > 2) {
619  for (int i = 0; i < npar; ++i) {
620  for (int j = 0; j < npar; ++j) {
621  B2INFO("Minv[" << i << "," << j << "]=" << gsl_matrix_get(Minv, i, j));
622  }
623  }
624  }
625 
626  B2DEBUG(10, "OPALFitterGSL: calcerr = " << calcerr);
627 
628  if (calcerr) {
629 
630 // *-- As a first step, calculate Minv as in 9.4.2 of Benno's book chapter
631 // (in O.Behnke et al "Data Analysis in High Energy Physics")
632 
633 
634 // *-- Evaluate S and invert.
635 
636  if (debug > 2) {
637  debug_print(&Vetaeta.matrix, "V");
638  debug_print(&Feta.matrix, "Feta");
639  }
640 
641  // CblasRight means C = alpha B A + beta C with symmetric matrix A
642  //FetaV[ncon][nmea] = 1*Feta[ncon][nmea]*V[nmea][nmea] + 0*FetaV
643  gsl_blas_dsymm(CblasRight, CblasUpper, 1, &Vetaeta.matrix, &Feta.matrix, 0, FetaV);
644  // S[ncon][ncon] = 1 * FetaV[ncon][nmea] * Feta^T[nmea][ncon] + 0*S
645  gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1, FetaV, &Feta.matrix, 0, S);
646 
647 
648  if (nunm > 0) {
649  // New invention by B. List, 6.12.04:
650  // add F_xi * F_xi^T to S, to make method work when some
651  // constraints do not depend on any measured parameter
652 
653  // Fxi is the part of Fetaxi containing the unmeasured quantities, if any
654  gsl_matrix_view Fxi = gsl_matrix_submatrix(Fetaxi, 0, nmea, ncon, nunm);
655  //S[ncon][ncon] = 1*Fxi[ncon][nunm]*Fxi^T[nunm][ncon] + 1*S[ncon][ncon]
656  gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1, &Fxi.matrix, &Fxi.matrix, 1, S);
657  }
658 
659  if (debug > 2) debug_print(S, "S");
660 
661 // *-- Invert S, testing for singularity first.
662 // S is symmetric and positive definite
663 
664  int signum;
665  gsl_linalg_LU_decomp(S, permS, &signum);
666  int inverr = gsl_linalg_LU_invert(S, permS, Sinv);
667 
668  if (inverr != 0) {
669  B2ERROR("S: gsl_linalg_LU_invert error " << inverr << " in error calculation");
670  ierr = -1;
671  return -1;
672  }
673 
674 
675 // *-- Calculate G.
676 // (same as W1, but for measured parameters)
677 // G = Feta^T * Sinv * Feta
678 
679  // SinvFeta[ncon][nmea] = 1*Sinv[ncon][ncon]*Feta[ncon][nmea] + 0*SinvFeta
680  gsl_blas_dsymm(CblasLeft, CblasUpper, 1, Sinv, &Feta.matrix, 0, SinvFeta);
681  // G[nmea][nmea] = 1*Feta^T[nmea][ncon]*SinvFeta[ncon][nmea] + 0*G
682  gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1, &Feta.matrix, SinvFeta, 0, G);
683 
684  if (debug > 2) debug_print(G, "G(1)");
685 
686 
687  if (nunm > 0) {
688 
689 // *-- Calculate H
690 
691  // Fxi is the part of Fetaxi containing the unmeasured quantities, if any
692  gsl_matrix_view Fxi = gsl_matrix_submatrix(Fetaxi, 0, nmea, ncon, nunm);
693  // H = Feta^T * Sinv * Fxi
694  // SinvFxi[ncon][nunm] = 1*Sinv[ncon][ncon]*Fxi[ncon][nunm] + 0*SinvFxi
695  gsl_blas_dsymm(CblasLeft, CblasUpper, 1, Sinv, &Fxi.matrix, 0, SinvFxi);
696  // H[nmea][nunm] = 1*Feta^T[nmea][ncon]*SinvFxi[ncon][nunm] + 0*H
697  gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1, &Feta.matrix, SinvFxi, 0, H);
698 
699  if (debug > 2) debug_print(H, "H");
700 
701 // *-- Calculate U**-1 and invert.
702 // (same as W1)
703 // U is a part of Minv
704 
705  gsl_matrix* Uinv = W1;
706  gsl_matrix_view U = gsl_matrix_submatrix(Minv, nmea, nmea, nunm, nunm);
707  // Uinv = Fxi^T * Sinv * Fxi
708  // Uinv[nunm][nunm] = 1*Fxi^T[nunm][ncon]*SinvFxi[ncon][nunm] + 0*W1
709  gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1, &Fxi.matrix, SinvFxi, 0, Uinv);
710 
711  gsl_linalg_LU_decomp(Uinv, permU, &signum);
712  inverr = gsl_linalg_LU_invert(Uinv, permU, &U.matrix);
713 
714  if (debug > 2) debug_print(&U.matrix, "U");
715  for (int i = 0; i < npar; ++i) {
716  for (int j = 0; j < npar; ++j) {
717  B2DEBUG(12, "after U Minv[" << i << "," << j << "]=" << gsl_matrix_get(Minv, i, j));
718  }
719  }
720 
721  if (inverr != 0) {
722  B2ERROR("U: gsl_linalg_LU_invert error " << inverr << " in error calculation ");
723  ierr = -1;
724  return -1;
725  }
726 
727 // *-- Covariance matrix between measured and unmeasured parameters.
728 
729 // HU[nmea][nunm] = 1*H[nmea][nunm]*U[nunm][nunm] + 0*HU
730  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, H, &U.matrix, 0, HU);
731 // Vnewetaxi is a view of Vnew
732  gsl_matrix_view Minvetaxi = gsl_matrix_submatrix(Minv, 0, nmea, nmea, nunm);
733  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, -1, &Vetaeta.matrix, HU, 0, &Minvetaxi.matrix);
734  for (int i = 0; i < npar; ++i) {
735  for (int j = 0; j < npar; ++j) {
736  B2DEBUG(12, "after etaxi Minv[" << i << "," << j << "]=" << gsl_matrix_get(Minv, i, j));
737  }
738  }
739 
740 
741 // *-- Fill in symmetric part:
742 // Vnewxieta is a view of Vnew
743  gsl_matrix_view Minvxieta = gsl_matrix_submatrix(Minv, nmea, 0, nunm, nmea);
744  gsl_matrix_transpose_memcpy(&Minvxieta.matrix, &Minvetaxi.matrix);
745  for (int i = 0; i < npar; ++i) {
746  for (int j = 0; j < npar; ++j) {
747  B2DEBUG(12, "after symmetric: Minv[" << i << "," << j << "]=" << gsl_matrix_get(Minv, i, j));
748  }
749  }
750 
751 // *-- Calculate G-HUH^T:
752 // G = -1*HU*H^T +1*G
753  gsl_blas_dgemm(CblasNoTrans, CblasTrans, -1, HU, H, +1, G);
754 
755  } // endif nunm > 0
756 
757 // *-- Calculate I-GV.
758  // IGV = 1
759  gsl_matrix_set_identity(IGV);
760  // IGV = -1*G*Vetaeta + 1*IGV
761  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, -1, G, &Vetaeta.matrix, 1, IGV);
762 
763 // *-- And finally error matrix on fitted parameters.
764  gsl_matrix_view Minvetaeta = gsl_matrix_submatrix(Minv, 0, 0, nmea, nmea);
765 
766  // Vnewetaeta = 1*Vetaeta*IGV + 0*Vnewetaeta
767  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, &Vetaeta.matrix, IGV, 0, &Minvetaeta.matrix);
768 
769  for (int i = 0; i < npar; ++i) {
770  for (int j = 0; j < npar; ++j) {
771  B2DEBUG(12, "complete Minv[" << i << "," << j << "]=" << gsl_matrix_get(Minv, i, j));
772  }
773  }
774 
775 // *-- now Minv should be complete, can calculate new covariance matrix Vnew
776 // Vnew = (dx / dt)^T * V_t * dx / dt
777 // x are the fitted parameters, t are the measured parameters,
778 // V_t is the covariance matrix of the measured parameters
779 // thus in OPALFitter variables: V_t = Vetaeta, x = (eta, xi) = etaxi
780 // d eta / dt and d xi / dt are given by eqns 9.54 and 9.55, respectively,
781 // as deta/dt = - Minvetaeta * Fetat and dxi/dt = - Minvxieta * Fetat
782 // Fetat is d^2 chi^2 / deta dt = - V^-1, even in presence of soft constraints, since
783 // the soft constraints don't depend on the _measured_ parameters t (but on the fitted ones)
784 // HERE, V (and Vinv) do not include the second derivatives of the soft constraints
785 // so we can use them right away
786 // Note that "F" is used for completely different things:
787 // in OPALFitter it is the derivatives of the constraints wrt to the parameters (A in book)
788 // in the book, F are the second derivatives of the objective function wrt the parameters
789 
790 
791  gsl_matrix_view detadt = gsl_matrix_submatrix(dxdt, 0, 0, nmea, nmea);
792  B2DEBUG(13, "after detadt");
793  // Vdxdt is Vetaeta * dxdt^T, thus Vdxdt[nmea][npar]
794  gsl_matrix_view Vdetadt = gsl_matrix_submatrix(Vdxdt, 0, 0, nmea, nmea);
795  B2DEBUG(13, "after Vdetadt");
796 
797  // detadt = - Minvetaeta * Fetat = -1 * Minvetaeta * (-1) * Vinv + 0 * detadt // replace by symm?
798  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, &Minvetaeta.matrix, Vinv, 0, &detadt.matrix);
799  if (debug > 2) debug_print(&detadt.matrix, "deta/dt");
800 
801  // Vdetadt = 1 * Vetaeta * detadt^T + 0* Vdetadt
802  gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1, &Vetaeta.matrix, &detadt.matrix, 0, &Vdetadt.matrix); // ok
803  if (debug > 2) debug_print(&Vdetadt.matrix, "Vetata * deta/dt");
804 
805  gsl_matrix_view Vnewetaeta = gsl_matrix_submatrix(Vnew, 0, 0, nmea, nmea); //[nmea],[nmea]
806  // Vnewetaeta = 1 * detadt * Vdetadt + 0* Vnewetaeta
807  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, &detadt.matrix, &Vdetadt.matrix, 0, &Vnewetaeta.matrix);
808 
809  if (debug > 2) debug_print(Vnew, "Vnew after part for measured parameters");
810 
811 
812  if (nunm > 0) {
813 
814  gsl_matrix_view Minvxieta = gsl_matrix_submatrix(Minv, nmea, 0, nunm, nmea); //[nunm][nmea]
815  B2DEBUG(13, "after Minvxieta");
816  if (debug > 2) debug_print(&Minvxieta.matrix, "Minvxieta");
817 
818  gsl_matrix_view dxidt = gsl_matrix_submatrix(dxdt, nmea, 0, nunm, nmea); //[nunm][nmea]
819  B2DEBUG(13, "after dxidt");
820  // dxidt[nunm][nmea] = - Minvxieta * Fetat = -1 * Minvxieta[nunm][nmea] * Vinv[nmea][nmea] + 0 * dxidt
821  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, &Minvxieta.matrix, Vinv, 0, &dxidt.matrix); //ok
822  if (debug > 2) debug_print(&dxidt.matrix, "dxi/dt");
823 
824  // Vdxdt = V * dxdt^T => Vdxdt[nmea][npar]
825  gsl_matrix_view Vdxidt = gsl_matrix_submatrix(Vdxdt, 0, nmea, nmea, nunm); //[nmea][nunm]
826  B2DEBUG(13, "after Vdxidt");
827  // Vdxidt = 1 * Vetaeta[nmea][nmea] * dxidt^T[nmea][nunm] + 0* Vdxidt => Vdxidt[nmea][nunm]
828  gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1, &Vetaeta.matrix, &dxidt.matrix, 0, &Vdxidt.matrix); // ok
829  if (debug > 2) debug_print(&Vdxidt.matrix, "Vetaeta * dxi/dt^T");
830 
831  gsl_matrix_view Vnewetaxi = gsl_matrix_submatrix(Vnew, 0, nmea, nmea, nunm); //[nmea][nunm]
832  gsl_matrix_view Vnewxieta = gsl_matrix_submatrix(Vnew, nmea, 0, nunm, nmea); //[nunm][nmea]
833  gsl_matrix_view Vnewxixi = gsl_matrix_submatrix(Vnew, nmea, nmea, nunm, nunm); //[nunm][nunm]
834 
835  // Vnewxieta[nunm][nmea] = 1 * dxidt[nunm][nmea] * Vdetadt[nmea][nmea] + 0* Vnewxieta
836  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, &dxidt.matrix, &Vdetadt.matrix, 0, &Vnewxieta.matrix); // ok
837  if (debug > 2) debug_print(Vnew, "Vnew after xieta part");
838  // Vnewetaxi[nmea][nunm] = 1 * detadt[nmea][nmea] * Vdxidt[nmea][nunm] + 0* Vnewetaxi
839  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, &detadt.matrix, &Vdxidt.matrix, 0, &Vnewetaxi.matrix); // ok
840  if (debug > 2) debug_print(Vnew, "Vnew after etaxi part");
841  // Vnewxixi[nunm][nunm] = 1 * dxidt[nunm][nmea] * Vdxidt[nmea][nunm] + 0* Vnewxixi
842  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, &dxidt.matrix, &Vdxidt.matrix, 0, &Vnewxixi.matrix);
843  if (debug > 2) debug_print(Vnew, "Vnew after xixi part");
844  }
845 
846 // *-- now we finally have Vnew, fill into fitobjects
847 
848  // update errors in fitobjects
849  for (auto& fitobject : fitobjects) {
850  for (int ilocal = 0; ilocal < fitobject->getNPar(); ++ilocal) {
851  int iglobal = fitobject->getGlobalParNum(ilocal);
852  for (int jlocal = ilocal; jlocal < fitobject->getNPar(); ++jlocal) {
853  int jglobal = fitobject->getGlobalParNum(jlocal);
854  if (iglobal >= 0 && jglobal >= 0)
855  fitobject->setCov(ilocal, jlocal, gsl_matrix_get(Vnew, iglobal, jglobal));
856  }
857  }
858  }
859 
860  // Finally, copy covariance matrix
861  if (cov && covDim != nmea + nunm) {
862  delete[] cov;
863  cov = nullptr;
864  }
865  covDim = nmea + nunm;
866  if (!cov) cov = new double[covDim * covDim];
867  for (int i = 0; i < covDim; ++i) {
868  for (int j = 0; j < covDim; ++j) {
869  cov[i * covDim + j] = gsl_matrix_get(Vnew, i, j);
870  }
871  }
872  covValid = true;
873 
874  } // endif calcerr == true
875 
876 #ifndef FIT_TRACEOFF
877  if (tracer) tracer->finish(*this);
878 #endif
879 
880  return fitprob;
881 
882  }
883 
884  bool OPALFitterGSL::initialize()
885  {
886  covValid = false;
887  // tell fitobjects the global ordering of their parameters:
888  int iglobal = 0;
889  // measured parameters first!
890  for (auto& fitobject : fitobjects) {
891  for (int ilocal = 0; ilocal < fitobject->getNPar(); ++ilocal) {
892  if (fitobject->isParamMeasured(ilocal) &&
893  !fitobject->isParamFixed(ilocal)) {
894  fitobject->setGlobalParNum(ilocal, iglobal);
895  B2DEBUG(10, "Object " << fitobject->getName()
896  << " Parameter " << fitobject->getParamName(ilocal)
897  << " is measured, global number " << iglobal);
898  ++iglobal;
899  }
900  }
901  }
902  nmea = iglobal;
903  // now unmeasured parameters!
904  for (auto& fitobject : fitobjects) {
905  for (int ilocal = 0; ilocal < fitobject->getNPar(); ++ilocal) {
906  if (!fitobject->isParamMeasured(ilocal) &&
907  !fitobject->isParamFixed(ilocal)) {
908  fitobject->setGlobalParNum(ilocal, iglobal);
909  B2DEBUG(10, "Object " << fitobject->getName()
910  << " Parameter " << fitobject->getParamName(ilocal)
911  << " is unmeasured, global number " << iglobal);
912  ++iglobal;
913  }
914  }
915  }
916  npar = iglobal;
917  assert(npar <= NPARMAX);
918  nunm = npar - nmea;
919  assert(nunm <= NUNMMAX);
920 
921  // set number of constraints
922  ncon = constraints.size();
923  assert(ncon <= NCONMAX);
924 
925  ini_gsl_vector(f, ncon);
926  ini_gsl_vector(r, ncon);
927 
928  ini_gsl_matrix(Fetaxi, ncon, npar);
929  ini_gsl_matrix(S, ncon, ncon);
930  ini_gsl_matrix(Sinv, ncon, ncon);
931  ini_gsl_matrix(SinvFxi, ncon, nunm);
932  ini_gsl_matrix(SinvFeta, ncon, nmea);
933  ini_gsl_matrix(W1, nunm, nunm);
934  ini_gsl_matrix(G, nmea, nmea);
935  ini_gsl_matrix(H, nmea, nunm);
936  ini_gsl_matrix(HU, nmea, nunm);
937  ini_gsl_matrix(IGV, nmea, nmea);
938  ini_gsl_matrix(V, npar, npar);
939  ini_gsl_matrix(VLU, nmea, nmea);
940  ini_gsl_matrix(Vinv, nmea, nmea);
941  ini_gsl_matrix(Vnew, npar, npar);
942  ini_gsl_matrix(Minv, npar, npar);
943  ini_gsl_matrix(dxdt, npar, nmea);
944  ini_gsl_matrix(Vdxdt, nmea, npar);
945 
946  ini_gsl_vector(dxi, nunm);
947  ini_gsl_vector(Fxidxi, ncon);
948  ini_gsl_vector(lambda, ncon);
949  ini_gsl_vector(FetaTlambda, nmea);
950 
951  ini_gsl_vector(etaxi, npar);
952  ini_gsl_vector(etasv, npar);
953  ini_gsl_vector(y, nmea);
954  ini_gsl_vector(y_eta, nmea);
955  ini_gsl_vector(Vinvy_eta, nmea);
956 
957  ini_gsl_matrix(FetaV, ncon, nmea);
958 
959  ini_gsl_permutation(permS, ncon);
960  ini_gsl_permutation(permU, nunm);
961  ini_gsl_permutation(permV, nmea);
962 
963  assert(f && (int)f->size == ncon);
964  assert(r && (int)r->size == ncon);
965  assert(Fetaxi && (int)Fetaxi->size1 == ncon && (int)Fetaxi->size2 == npar);
966  assert(S && (int)S->size1 == ncon && (int)S->size2 == ncon);
967  assert(Sinv && (int)Sinv->size1 == ncon && (int)Sinv->size2 == ncon);
968  assert(nunm == 0 || (SinvFxi && (int)SinvFxi->size1 == ncon && (int)SinvFxi->size2 == nunm));
969  assert(SinvFeta && (int)SinvFeta->size1 == ncon && (int)SinvFeta->size2 == nmea);
970  assert(nunm == 0 || (W1 && (int)W1->size1 == nunm && (int)W1->size2 == nunm));
971  assert(G && (int)G->size1 == nmea && (int)G->size2 == nmea);
972  assert(nunm == 0 || (H && (int)H->size1 == nmea && (int)H->size2 == nunm));
973  assert(nunm == 0 || (HU && (int)HU->size1 == nmea && (int)HU->size2 == nunm));
974  assert(IGV && (int)IGV->size1 == nmea && (int)IGV->size2 == nmea);
975  assert(V && (int)V->size1 == npar && (int)V->size2 == npar);
976  assert(VLU && (int)VLU->size1 == nmea && (int)VLU->size2 == nmea);
977  assert(Vinv && (int)Vinv->size1 == nmea && (int)Vinv->size2 == nmea);
978  assert(Vnew && (int)Vnew->size1 == npar && (int)Vnew->size2 == npar);
979  assert(nunm == 0 || (dxi && (int)dxi->size == nunm));
980  assert(nunm == 0 || (Fxidxi && (int)Fxidxi->size == ncon));
981  assert(lambda && (int)lambda->size == ncon);
982  assert(FetaTlambda && (int)FetaTlambda->size == nmea);
983  assert(etaxi && (int)etaxi->size == npar);
984  assert(etasv && (int)etasv->size == npar);
985  assert(y && (int)y->size == nmea);
986  assert(y_eta && (int)y_eta->size == nmea);
987  assert(Vinvy_eta && (int)Vinvy_eta->size == nmea);
988  assert(FetaV && (int)FetaV->size1 == ncon && (int)FetaV->size2 == nmea);
989  assert(permS && (int)permS->size == ncon);
990  assert(permS && (int)permS->size == ncon);
991  assert(nunm == 0 || (permU && (int)permU->size == nunm));
992  assert(permV && (int)permV->size == nmea);
993 
994  return true;
995 
996  }
997 
998  bool OPALFitterGSL::updateFitObjects(double eetaxi[])
999  {
1000  // changed etaxi -> eetaxi to avoid clash with class member DJeans
1001  //bool debug = false;
1002  bool result = true;
1003  for (auto& fitobject : fitobjects) {
1004  for (int ilocal = 0; ilocal < fitobject->getNPar(); ++ilocal) {
1005  fitobject->updateParams(eetaxi, npar);
1006 // }
1007  }
1008  }
1009  return result;
1010  }
1011 
1012  int OPALFitterGSL::getError() const {return ierr;}
1013  double OPALFitterGSL::getProbability() const {return fitprob;}
1014  double OPALFitterGSL::getChi2() const {return chi2;}
1015  int OPALFitterGSL::getDoF() const {return ncon - nunm;}
1016  int OPALFitterGSL::getIterations() const {return nit;}
1017 
1018  void OPALFitterGSL::ini_gsl_permutation(gsl_permutation*& p, unsigned int size)
1019  {
1020  if (p) {
1021  if (p->size != size) {
1022  gsl_permutation_free(p);
1023  if (size > 0) p = gsl_permutation_alloc(size);
1024  else p = nullptr;
1025  }
1026  } else if (size > 0) p = gsl_permutation_alloc(size);
1027  }
1028 
1029  void OPALFitterGSL::ini_gsl_vector(gsl_vector*& v, unsigned int size)
1030  {
1031 
1032  if (v) {
1033  if (v->size != size) {
1034  gsl_vector_free(v);
1035  if (size > 0) v = gsl_vector_alloc(size);
1036  else v = nullptr;
1037  }
1038  } else if (size > 0) v = gsl_vector_alloc(size);
1039  }
1040 
1041  void OPALFitterGSL::ini_gsl_matrix(gsl_matrix*& m, unsigned int size1, unsigned int size2)
1042  {
1043  if (m) {
1044  if (m->size1 != size1 || m->size2 != size2) {
1045  gsl_matrix_free(m);
1046  if (size1 * size2 > 0) m = gsl_matrix_alloc(size1, size2);
1047  else m = nullptr;
1048  }
1049  } else if (size1 * size2 > 0) m = gsl_matrix_alloc(size1, size2);
1050  }
1051 
1052  void OPALFitterGSL::setDebug(int debuglevel)
1053  {
1054  debug = debuglevel;
1055  }
1056 
1057  void OPALFitterGSL::debug_print(gsl_matrix* m, const char* name)
1058  {
1059  for (unsigned int i = 0; i < m->size1; ++i)
1060  for (unsigned int j = 0; j < m->size2; ++j)
1061  if (gsl_matrix_get(m, i, j) != 0)
1062  B2INFO(name << "[" << i << "][" << j << "]=" << gsl_matrix_get(m, i, j));
1063  }
1064 
1065  void OPALFitterGSL::debug_print(gsl_vector* v, const char* name)
1066  {
1067  for (unsigned int i = 0; i < v->size; ++i)
1068  if (gsl_vector_get(v, i) != 0)
1069  B2INFO(name << "[" << i << "]=" << gsl_vector_get(v, i));
1070  }
1071 
1072  int OPALFitterGSL::getNcon() const {return ncon;}
1073  int OPALFitterGSL::getNsoft() const {return 0;}
1074  int OPALFitterGSL::getNunm() const {return nunm;}
1075  int OPALFitterGSL::getNpar() const {return npar;}
1076 
1077  }// end OrcaKinFit namespace
1079 } // end Belle2 namespace
1080 
R E
internal precision of FFTW codelets
Abstract base class for different kinds of events.