Belle II Software development
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
36using std::endl;
37using std::abs;
38
39//using namespace Belle2 {
40//using namespace OrcaKinFit {
41
42namespace 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
E
R E
internal precision of FFTW codelets
Definition: DiscreteCosineTransform_31points.cc:36
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17