Description of the fit algorithm and interface: More...

#include <OPALFitterGSL.h>

Inheritance diagram for OPALFitterGSL:

Collaboration diagram for OPALFitterGSL:

[legend]

Public Member Functions
virtual double	fit () override

virtual int	getError () const override
	Return error code. More...

virtual double	getProbability () const override
	Get the fit probability of the last fit.

virtual double	getChi2 () const override
	Get the chi**2 of the last fit.

virtual int	getDoF () const override
	Get the number of degrees of freedom of the last fit.

virtual int	getIterations () const override
	Get the number of iterations of the last fit.

virtual int	getNcon () const
	Get the number of hard constraints of the last fit.

virtual int	getNsoft () const
	Get the number of soft constraints of the last fit.

virtual int	getNpar () const
	Get the number of all parameters of the last fit.

virtual int	getNunm () const
	Get the number of unmeasured parameters of the last fit.

virtual bool	initialize () override

virtual void	setDebug (int debuglevel)
	Set the Debug Level.

virtual void	addFitObject (BaseFitObject *fitobject_)

virtual void	addFitObject (BaseFitObject &fitobject_)

virtual void	addConstraint (BaseConstraint *constraint_)

virtual void	addConstraint (BaseConstraint &constraint_)

virtual void	addHardConstraint (BaseHardConstraint *constraint_)

virtual void	addHardConstraint (BaseHardConstraint &constraint_)

virtual void	addSoftConstraint (BaseSoftConstraint *constraint_)

virtual void	addSoftConstraint (BaseSoftConstraint &constraint_)

virtual std::vector< BaseFitObject * > *	getFitObjects ()

virtual std::vector< BaseHardConstraint * > *	getConstraints ()

virtual std::vector< BaseSoftConstraint * > *	getSoftConstraints ()

virtual void	reset ()

virtual BaseTracer *	getTracer ()

virtual const BaseTracer *	getTracer () const

virtual void	setTracer (BaseTracer *newTracer)

virtual void	setTracer (BaseTracer &newTracer)

virtual const double *	getGlobalCovarianceMatrix (int &idim) const

virtual double *	getGlobalCovarianceMatrix (int &idim)

Public Attributes
std::map< std::string, double >	traceValues

Protected Types
enum	{ NPARMAX = 50 , NCONMAX = 20 , NUNMMAX = 20 }

typedef std::vector< BaseFitObject * >	FitObjectContainer

typedef std::vector< BaseHardConstraint * >	ConstraintContainer

typedef std::vector< BaseSoftConstraint * >	SoftConstraintContainer

typedef FitObjectContainer::iterator	FitObjectIterator

typedef ConstraintContainer::iterator	ConstraintIterator

typedef SoftConstraintContainer::iterator	SoftConstraintIterator

Protected Member Functions
virtual bool	updateFitObjects (double etaxi[])

Static Protected Member Functions
static void	ini_gsl_permutation (gsl_permutation *&p, unsigned int size)

static void	ini_gsl_vector (gsl_vector *&v, int unsigned size)

static void	ini_gsl_matrix (gsl_matrix *&m, int unsigned size1, unsigned int size2)

static void	debug_print (gsl_matrix m, const char name)

static void	debug_print (gsl_vector v, const char name)

Protected Attributes
int	npar

int	nmea

int	nunm

int	ncon

int	ierr

int	nit

double	fitprob

double	chi2

FitObjectContainer	fitobjects

ConstraintContainer	constraints

SoftConstraintContainer	softconstraints

int	covDim
	dimension of global covariance matrix

double *	cov
	global covariance matrix of last fit problem

bool	covValid
	Flag whether global covariance is valid.

BaseTracer *	tracer

Private Attributes
gsl_vector *	f

gsl_vector *	r

gsl_matrix *	Fetaxi

gsl_matrix *	S

gsl_matrix *	Sinv

gsl_matrix *	SinvFxi

gsl_matrix *	SinvFeta

gsl_matrix *	W1

gsl_matrix *	G

gsl_matrix *	H

gsl_matrix *	HU

gsl_matrix *	IGV

gsl_matrix *	V

gsl_matrix *	VLU

gsl_matrix *	Vinv

gsl_matrix *	Vnew

gsl_matrix *	Minv

gsl_matrix *	dxdt

gsl_matrix *	Vdxdt

gsl_vector *	dxi

gsl_vector *	Fxidxi

gsl_vector *	lambda

gsl_vector *	FetaTlambda

gsl_vector *	etaxi

gsl_vector *	etasv

gsl_vector *	y

gsl_vector *	y_eta

gsl_vector *	Vinvy_eta

gsl_matrix *	FetaV

gsl_permutation *	permS
	Helper permutation vector.

gsl_permutation *	permU
	Helper permutation vector.

gsl_permutation *	permV
	Helper permutation vector.

int	debug

Detailed Description

Description of the fit algorithm and interface:

The OPALFitter object holds a set (fitobjects) of BaseFitObject objects, which represent the objects (particles, jets) whose fourvectors shall be fitted. It also holds a set (constraints) of BaseConstraint objects that represent the constraints (momentum or mass constraints).

OPALFitter::initialize goes over the list of fit objects and determines the number of parameters (measured and unmeasured), and assigns global numbers to them.

OPALFitter::fit first initializes the global parameter numbering using OPALFitter::initialize.

Then it initializes vectors etaxi ( $\eta\xi$ ) and y ( $\vec y$ ) with the current parameter values.

Next it initializes the matrix Feta that represents $d \vec F / d \vec \eta\xi$

Used methods in OPALFitter::initialize:

Used methods in OPALFitter::updateFitObjects:

Used methods in OPALFitter::fit:

Interaction between Constraints and Particles

The Fitter does not care how the constraints and BaseFitObject objects interact. The constraints are expressed in terms of four-vector components of the BaseFitObjects. A constraint object keeps a set of BaseFitObject objects, and, using the chain rule, calculates its derivatives w.r.t. the global parameters.

Definition at line 88 of file OPALFitterGSL.h.

Member Function Documentation

◆ fit()

double fit ( )

overridevirtual

initialize Fetaxi ( = d F / d eta,xi)

Implements BaseFitter.

Definition at line 132 of file OPALFitterGSL.cc.

     {
  
  
  
       //
       //           (     )   ^     ^
       //           ( eta )  nmea   |
       //           (     )   v     |
       //   etaxi = (-----)  ---   npar
       //           (     )   ^     |
       //           ( xi  )  nunm   |
       //           (     )   v     v
  
       //              <- ncon ->
       //             (          )   ^     ^
       //             (   Feta^T )  nmea   |
       //             (          )   v     |
       //  Fetaxi^T = ( -------- )  ---   npar
       //             (          )   ^     |
       //             (   Fxi^T  )  nunm   |
       //             (          )   v     v
       //
       //  Fetaxi are the derivatives of the constraints wrt the fitted parameters, thus A_theta/xi in book
  
  
       //
       //            <- nmea ->|<- nunm ->
       //           (          |          )   ^     ^
       //           ( Vetaeta  |  Vetaxi  )  nmea   |
       //           (          |          )   v     |
       //  V =      (----------+----------)  ---   npar
       //           (          |          )   ^     |
       //           (  Vxieta  |  Vxixi   )  nunm   |
       //           (          |          )   v     v
       //
       //  V is the covariance matrix. Before the fit only Vetaeta is non-zero (covariance of measured parameters)
  
       //
       //                 <- nmea ->|<- nunm ->
       //                (          |         )   ^
       //     dxdt^T =   (   detadt |   dxidt )  nmea
       //                (          |         )   v
       //
       //  dxdt are the partial derivates of the fitted parameters wrt the measured parameters,
       //
       //                 <- nmea ->|<- nunm ->
       //                (          |          )   ^
       //     Vdxdt  =   ( Vdetadt  |  Vdxidt  )  nmea
       //                (          |          )   v
       //
       //  Vdxdt is Vetaeta * dxdt^T, thus Vdxdt[nmea][npar]
       //  both needed for calculation of the full covariance matrix of the fitted parameters
  
       // order parameters etc
       initialize();
  
       assert(f && (int)f->size == ncon);
       assert(r && (int)r->size == ncon);
       assert(Fetaxi && (int)Fetaxi->size1 == ncon && (int)Fetaxi->size2 == npar);
       assert(S && (int)S->size1 == ncon && (int)S->size2 == ncon);
       assert(Sinv && (int)Sinv->size1 == ncon && (int)Sinv->size2 == ncon);
       assert(nunm == 0 || (SinvFxi && (int)SinvFxi->size1 == ncon && (int)SinvFxi->size2 == nunm));
       assert(SinvFeta && (int)SinvFeta->size1 == ncon && (int)SinvFeta->size2 == nmea);
       assert(nunm == 0 || (W1 && (int)W1->size1 == nunm && (int)W1->size2 == nunm));
       assert(G && (int)G->size1 == nmea && (int)G->size2 == nmea);
       assert(nunm == 0 || (H && (int)H->size1 == nmea && (int)H->size2 == nunm));
       assert(nunm == 0 || (HU && (int)HU->size1 == nmea && (int)HU->size2 == nunm));
       assert(IGV && (int)IGV->size1 == nmea && (int)IGV->size2 == nmea);
       assert(V && (int)V->size1 == npar && (int)V->size2 == npar);
       assert(VLU && (int)VLU->size1 == nmea && (int)VLU->size2 == nmea);
       assert(Vinv && (int)Vinv->size1 == nmea && (int)Vinv->size2 == nmea);
       assert(Vnew && (int)Vnew->size1 == npar && (int)Vnew->size2 == npar);
       assert(Minv && (int)Minv->size1 == npar && (int)Minv->size2 == npar);
       assert(dxdt && (int)dxdt->size1 == npar && (int)dxdt->size2 == nmea);
       assert(Vdxdt  && (int)Vdxdt->size1  == nmea && (int)Vdxdt->size2  == npar);
       assert(nunm == 0 || (dxi && (int)dxi->size == nunm));
       assert(nunm == 0 || (Fxidxi && (int)Fxidxi->size == ncon));
       assert(lambda && (int)lambda->size == ncon);
       assert(FetaTlambda && (int)FetaTlambda->size == nmea);
       assert(etaxi && (int)etaxi->size == npar);
       assert(etasv && (int)etasv->size == npar);
       assert(y && (int)y->size == nmea);
       assert(y_eta && (int)y_eta->size == nmea);
       assert(Vinvy_eta && (int)Vinvy_eta->size == nmea);
       assert(FetaV && (int)FetaV->size1 == ncon && (int)FetaV->size2 == nmea);
       assert(permS && (int)permS->size == ncon);
       assert(nunm == 0 || (permU && (int)permU->size == nunm));
       assert(permV && (int)permV->size == nmea);
  
       // eta is the part of etaxi containing the measured quantities
       gsl_vector_view eta = gsl_vector_subvector(etaxi, 0, nmea);
  
       // Feta is the part of Fetaxi containing the measured quantities
  
       B2DEBUG(11, "==== " << ncon << " " << nmea);
  
       gsl_matrix_view Feta = gsl_matrix_submatrix(Fetaxi, 0, 0, ncon, nmea);
  
       gsl_matrix_view Vetaeta   = gsl_matrix_submatrix(V, 0, 0, nmea, nmea);
  
       for (unsigned int i = 0; i < fitobjects.size(); ++i) {
         for (int ilocal = 0; ilocal < fitobjects[i]->getNPar(); ++ilocal) {
           if (!fitobjects[i]->isParamFixed(ilocal)) {
             int iglobal = fitobjects[i]->getGlobalParNum(ilocal);
             assert(iglobal >= 0 && iglobal < npar);
             gsl_vector_set(etaxi, iglobal, fitobjects[i]->getParam(ilocal));
             if (fitobjects[i]->isParamMeasured(ilocal)) {
               assert(iglobal < nmea);
               gsl_vector_set(y, iglobal, fitobjects[i]->getMParam(ilocal));
             }
             B2DEBUG(10, "etaxi[" << iglobal << "] = " << gsl_vector_get(etaxi, iglobal)
                     << " for jet " << i << " and ilocal = " << ilocal);
           }
         }
       }
  
       gsl_matrix_set_zero(Fetaxi);
       for (int k = 0; k < ncon; k++) {
         constraints[k]->getDerivatives(Fetaxi->size2, Fetaxi->block->data + k * Fetaxi->tda);
         if (debug > 1) for (int j = 0; j < npar; j++)
             if (gsl_matrix_get(Fetaxi, k, j) != 0)
               B2INFO("1: Fetaxi[" << k << "][" << j << "] = " << gsl_matrix_get(Fetaxi, k, j));
       }
  
       // chi2's, step size, # iterations
       double chinew = 0, chit = 0, chik = 0;
       double alph = 1.;
       nit = 0;
       // convergence criteria (as in WWINIT)
       int nitmax = 200;
       double chik0 = 100.;
       double chit0 = 100.;
       double dchikc = 1.0E-3;
       double dchitc = 1.0E-4;
       double dchikt = 1.0E-2;
       double dchik  = 1.05;
       double chimxw = 10000.;
       double almin = 0.05;
  
       // repeat with or with out smaller steps size
       bool repeat = true;
       bool scut = false;
       bool calcerr = true;
  
 #ifndef FIT_TRACEOFF
       if (tracer) tracer->initialize(*this);
 #endif
  
  
       // start of iterations
       while (repeat) {
         bool updatesuccess = true;
  
 //*-- If necessary, retry smaller step, same direction
         if (scut) {
           gsl_vector_memcpy(etaxi, etasv);
           updatesuccess = updateFitObjects(etaxi->block->data);
           if (!updatesuccess) {
             B2ERROR("OPALFitterGSL::fit: old parameters are garbage!");
             return -1;
           }
  
  
           gsl_matrix_set_zero(Fetaxi);
           for (int k = 0; k < ncon; ++k)  {
             constraints[k]->getDerivatives(Fetaxi->size2, Fetaxi->block->data + k * Fetaxi->tda);
           }
           if (debug > 1) debug_print(Fetaxi, "1: Fetaxi");
         } else {
           gsl_vector_memcpy(etasv, etaxi);
           chik0 = chik;
           chit0 = chit;
         }
  
         // Get covariance matrix
         gsl_matrix_set_zero(V);
  
         for (auto& fitobject : fitobjects) {
           fitobject->addToGlobCov(V->block->data, V->tda);
         }
         if (debug > 1)  debug_print(V, "V");
  
         gsl_matrix_memcpy(VLU, &Vetaeta.matrix);
  
         // invert covariance matrix (needed for chi2 calculation later)
  
         int signum;
         int result;
         result = gsl_linalg_LU_decomp(VLU, permV, &signum);
         B2DEBUG(11, "gsl_linalg_LU_decomp result=" << result);
         if (debug > 3)  debug_print(VLU, "VLU");
  
         result = gsl_linalg_LU_invert(VLU, permV, Vinv);
         B2DEBUG(11, "gsl_linalg_LU_invert result=" << result);
  
         if (debug > 2) debug_print(Vinv, "Vinv");
  
 // *-- Evaluate f and S.
         for (int k = 0; k < ncon; ++k) {
           gsl_vector_set(f, k, constraints[k]->getValue());
         }
         if (debug > 1)  debug_print(f, "f");
  
         // y_eta = y - eta
         gsl_vector_memcpy(y_eta, y);
         gsl_vector_sub(y_eta, &eta.vector);
         // r=f
         gsl_vector_memcpy(r, f);
         // r = 1*Feta*y_eta + 1*r
         gsl_blas_dgemv(CblasNoTrans, 1, &Feta.matrix, y_eta, 1, r);
  
         if (debug > 1) debug_print(r, "r");
  
         // S = Feta * V * Feta^T
  
         //FetaV = 1*Feta*V + 0*FetaV
         B2DEBUG(12, "Creating FetaV");
         gsl_blas_dsymm(CblasRight, CblasUpper, 1, &Vetaeta.matrix, &Feta.matrix, 0,  FetaV);
         // S = 1 * FetaV * Feta^T + 0*S
         B2DEBUG(12, "Creating S");
         gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1, FetaV, &Feta.matrix, 0, S);
  
         if (nunm > 0) {
           // New invention by B. List, 6.12.04:
           // add F_xi * F_xi^T to S, to make method work when some
           // constraints do not depend on any measured parameter
  
           // Fxi is the part of Fetaxi containing the unmeasured quantities, if any
           gsl_matrix_view Fxi = gsl_matrix_submatrix(Fetaxi,  0, nmea, ncon, nunm);
  
           //S = 1*Fxi*Fxi^T + 1*S
           gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1, &Fxi.matrix, &Fxi.matrix, 1, S);
         }
  
         if (debug > 1) debug_print(S, "S");
  
 // *-- Invert S to Sinv; S is destroyed here!
 // S is symmetric and positive definite
  
         gsl_linalg_LU_decomp(S, permS, &signum);
         int inverr = gsl_linalg_LU_invert(S, permS, Sinv);
  
         if (inverr != 0) {
           B2ERROR("S: gsl_linalg_LU_invert error " << inverr);
           ierr = 7;
           calcerr = false;
           break;
         }
  
         // Calculate S^1*r here, we will need it
         // Store it in lambda!
         // lambda = 1*Sinv*r + 0*lambda; Sinv is symmetric
         gsl_blas_dsymv(CblasUpper, 1, Sinv, r, 0, lambda);
  
 // *-- Calculate new unmeasured quantities, if any
  
         if (nunm > 0) {
           // Fxi is the part of Fetaxi containing the unmeasured quantities, if any
           gsl_matrix_view Fxi = gsl_matrix_submatrix(Fetaxi,  0, nmea, ncon, nunm);
           // W1 = Fxi^T * Sinv * Fxi
           // SinvFxi = 1*Sinv*Fxi + 0*SinvFxi
           gsl_blas_dsymm(CblasLeft, CblasUpper, 1, Sinv, &Fxi.matrix, 0,  SinvFxi);
           // W1 = 1*Fxi^T*SinvFxi + 0*W1
           gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1, &Fxi.matrix, SinvFxi, 0, W1);
  
           if (debug > 1) {
             debug_print(W1, "W1");
             // Check symmetry of W1
             for (int i = 0; i < nunm; ++i) {
               for (int j = 0; j < nunm; ++j) {
                 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)))
                   B2INFO("W1[" << i << "][" << j << "] = " << gsl_matrix_get(W1, i, j)
                          << "   W1[" << j << "][" << i << "] = " << gsl_matrix_get(W1, j, i)
                          << "   => diff=" << abs(gsl_matrix_get(W1, i, j) - gsl_matrix_get(W1, j, i))
                          << "   => tol=" << 1E-3 * abs(gsl_matrix_get(W1, i, j) + gsl_matrix_get(W1, j, i)));
               }
             }
           }
  
           // calculate shift of unmeasured parameters
  
           // dxi = -alph*W1^-1 * Fxi^T*Sinv*r
           // => dxi is solution of W1*dxi = -alph*Fxi^T*Sinv*r
           // Compute rights hand side first
           // Sinv*r was already calculated and is stored in lambda
           // dxi = -alph*Fxi^T*lambda + 0*dxi
  
           B2DEBUG(11, "alph = " << alph);
           if (debug > 1) {
             debug_print(lambda, "lambda");
             debug_print(&(Fxi.matrix), "Fxi");
           }
  
           gsl_blas_dgemv(CblasTrans, -alph, &Fxi.matrix, lambda, 0, dxi);
  
           if (debug > 1) {
             debug_print(dxi, "dxi0");
             debug_print(W1, "W1");
           }
  
           // now solve the system
           // Note added 23.12.04: W1 is symmetric and positive definite,
           // so we can use the Cholesky instead of LU decomposition
           gsl_linalg_cholesky_decomp(W1);
           inverr = gsl_linalg_cholesky_svx(W1, dxi);
  
           if (debug > 1) debug_print(dxi, "dxi1");
  
  
           if (inverr != 0) {
             B2ERROR("W1: gsl_linalg_cholesky_svx error " << inverr);
             ierr = 8;
             calcerr = false;
             break;
           }
  
           if (debug > 1) debug_print(dxi, "dxi");
  
 //    *-- And now update unmeasured parameters; xi is a view of etaxi
           // xi is the part of etaxi containing the unmeasured quantities, if any
           gsl_vector_view xi = gsl_vector_subvector(etaxi, nmea, nunm);
  
           // xi += dxi
           gsl_vector_add(&xi.vector, dxi);
  
         }
  
 // *-- Calculate new Lagrange multipliers.
         // lambda = Sinv*r + Sinv*Fxi*dxi
         // lambda is already set to Sinv*r, we just need to add Sinv*Fxi*dxi
  
         // cppcheck-suppress duplicateCondition
         if (nunm > 0) {
           // Fxi is the part of Fetaxi containing the unmeasured quantities, if any
           gsl_matrix_view Fxi = gsl_matrix_submatrix(Fetaxi,  0, nmea, ncon, nunm);
           // calculate Fxidxi = 1*Fxi*dxi + 0*Fxidxi
           gsl_blas_dgemv(CblasNoTrans, 1, &Fxi.matrix, dxi, 0, Fxidxi);
           // add to existing lambda: lambda = 1*Sinv*Fxidxi + 1*lambda; Sinv is symmetric
           gsl_blas_dsymv(CblasUpper, 1, Sinv, Fxidxi, 1, lambda);
  
         }
  
         if (debug > 1) debug_print(lambda, "lambda");
  
 // *-- Calculate new fitted parameters:
         // eta = y - V*Feta^T*lambda
         gsl_vector_memcpy(&eta.vector, y);
         // FetaTlambda = 1*Feta^T*lambda + 0*FetaTlambda
         gsl_blas_dgemv(CblasTrans, 1, &Feta.matrix, lambda, 0, FetaTlambda);
         // eta = -1*V*FetaTlambda + 1*eta; V is symmetric
         gsl_blas_dsymv(CblasUpper, -1, &Vetaeta.matrix, FetaTlambda, 1, &eta.vector);
  
  
         if (debug > 1) debug_print(&eta.vector, "updated eta");
  
 // *-- Calculate constraints and their derivatives.
         // since the constraints ask the fitobjects for their parameters,
         // we need to update the fitobjects first!
         // COULD BE DONE: update also ERRORS! (now only in the very end!)
         updatesuccess = updateFitObjects(etaxi->block->data);
  
         B2DEBUG(10, "After adjustment of all parameters:\n");
         for (int k = 0; k < ncon; ++k) {
           B2DEBUG(10, "Value of constraint " << k << " = " << constraints[k]->getValue());
         }
         gsl_matrix_set_zero(Fetaxi);
         for (int k = 0; k < ncon; k++) {
           constraints[k]->getDerivatives(Fetaxi->size2, Fetaxi->block->data + k * Fetaxi->tda);
         }
         if (debug > 1)  debug_print(Fetaxi, "2: Fetaxi");
  
  
 // *-- Calculate new chisq.
         // First, calculate new y-eta
         // y_eta = y - eta
         gsl_vector_memcpy(y_eta, y);
         gsl_vector_sub(y_eta, &eta.vector);
         // Now calculate Vinv*y_eta [ as solution to V* Vinvy_eta = y_eta]
         // Vinvy_eta = 1*Vinv*y_eta + 0*Vinvy_eta; Vinv is symmetric
         gsl_blas_dsymv(CblasUpper, 1, Vinv, y_eta, 0, Vinvy_eta);
         // Now calculate y_eta *Vinvy_eta
         gsl_blas_ddot(y_eta, Vinvy_eta, &chit);
  
         for (int i = 0; i < nmea; ++i)
           for (int j = 0; j < nmea; ++j) {
             double dchit = (gsl_vector_get(y_eta, i)) *
                            gsl_matrix_get(Vinv, i, j) *
                            (gsl_vector_get(y_eta, j));
             if (dchit != 0)
               B2DEBUG(11, "chit for i,j = " << i << " , " << j << " = "
                       << dchit);
           }
  
         chik = 0;
         for (int k = 0; k < ncon; ++k) chik += std::abs(2 * gsl_vector_get(lambda, k) * constraints[k]->getValue());
         chinew = chit + chik;
  
 //*-- Calculate change in chisq, and check constraints are satisfied.
         nit++;
  
         bool sconv = (std::abs(chik - chik0) < dchikc)
                      && (std::abs(chit - chit0) < dchitc * chit)
                      && (chik < dchikt * chit);
         // Second convergence criterium:
         // If all constraints are fulfilled to better than 1E-8,
         // and all parameters have changed by less than 1E-8,
         // assume convergence
         // This criterium assumes that all constraints and all parameters
         // are "of order 1", i.e. their natural values are around 1 to 100,
         // as for GeV or radians
         double eps = 1E-6;
         bool sconv2 = true;
         for (int k = 0; sconv2 && (k < ncon); ++k)
           sconv2 &= (std::abs(gsl_vector_get(f, k)) < eps);
         if (sconv2)
           B2DEBUG(10, "All constraints fulfilled to better than " << eps);
  
         for (int j = 0; sconv2 && (j < npar); ++j)
           sconv2 &= (std::abs(gsl_vector_get(etaxi, j) - gsl_vector_get(etasv, j)) < eps);
         if (sconv2)
           B2DEBUG(10, "All parameters stable to better than " << eps);
         sconv |= sconv2;
  
         bool sbad  = (chik > dchik * chik0)
                      && (chik > dchikt * chit)
                      && (chik > chik0 + 1.E-10);
  
         scut = false;
  
         if (nit > nitmax) {
 // *-- Out of iterations
           repeat = false;
           calcerr = false;
           ierr = 1;
         } else if (sconv && updatesuccess) {
 // *-- Converged
           repeat = false;
           calcerr = true;
           ierr = 0;
         } else if (nit > 2 && chinew > chimxw  && updatesuccess) {
 // *-- Chi2  crazy ?
           repeat = false;
           calcerr = false;
           ierr = 2;
         } else if ((sbad && nit > 1) || !updatesuccess) {
 // *-- ChiK increased, try smaller step
           if (alph == almin) {
             repeat = true;   // false;
             ierr = 3;
           } else {
             alph  =  std::max(almin, 0.5 * alph);
             scut  =  true;
             repeat = true;
             ierr = 4;
           }
         } else {
 // *-- Keep going..
           alph  =  std::min(alph + 0.1, 1.);
           repeat = true;
           ierr = 5;
         }
  
         B2DEBUG(10, "======== NIT = " << nit << ",  CHI2 = " << chinew
                 << ",  ierr = " << ierr << ", alph=" << alph);
  
         for (unsigned int i = 0; i < fitobjects.size(); ++i)
           B2DEBUG(10, "fitobject " << i << ": " << *fitobjects[i]);
  
 #ifndef FIT_TRACEOFF
         if (tracer) tracer->step(*this);
 #endif
  
       }   // end of while (repeat)
  
 // *-- Turn chisq into probability.
       fitprob = (ncon - nunm > 0) ? gsl_cdf_chisq_Q(chinew, ncon - nunm) : 0.5;
       chi2 = chinew;
  
 // *-- End of iterations - calculate errors.
  
 // The result will (ultimately) be stored in Vnew
  
       gsl_matrix_set_zero(Vnew);
       gsl_matrix_set_zero(Minv);
       if (debug > 2) {
         for (int i = 0; i < npar; ++i) {
           for (int j = 0; j < npar; ++j) {
             B2INFO("Minv[" << i << "," << j << "]=" << gsl_matrix_get(Minv, i, j));
           }
         }
       }
  
       B2DEBUG(10, "OPALFitterGSL: calcerr = " << calcerr);
  
       if (calcerr) {
  
 // *-- As a first step, calculate Minv as in 9.4.2 of Benno's book chapter
 //                    (in O.Behnke et al "Data Analysis in High Energy Physics")
  
  
 // *-- Evaluate S and invert.
  
         if (debug > 2) {
           debug_print(&Vetaeta.matrix, "V");
           debug_print(&Feta.matrix, "Feta");
         }
  
         // CblasRight means C = alpha B A + beta C with symmetric matrix A
         //FetaV[ncon][nmea] = 1*Feta[ncon][nmea]*V[nmea][nmea] + 0*FetaV
         gsl_blas_dsymm(CblasRight, CblasUpper, 1, &Vetaeta.matrix, &Feta.matrix, 0,  FetaV);
         // S[ncon][ncon] = 1 * FetaV[ncon][nmea] * Feta^T[nmea][ncon] + 0*S
         gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1, FetaV, &Feta.matrix, 0, S);
  
  
         if (nunm > 0) {
           // New invention by B. List, 6.12.04:
           // add F_xi * F_xi^T to S, to make method work when some
           // constraints do not depend on any measured parameter
  
           // Fxi is the part of Fetaxi containing the unmeasured quantities, if any
           gsl_matrix_view Fxi = gsl_matrix_submatrix(Fetaxi,  0, nmea, ncon, nunm);
           //S[ncon][ncon] = 1*Fxi[ncon][nunm]*Fxi^T[nunm][ncon] + 1*S[ncon][ncon]
           gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1, &Fxi.matrix, &Fxi.matrix, 1, S);
         }
  
         if (debug > 2) debug_print(S, "S");
  
 // *-- Invert S, testing for singularity first.
 // S is symmetric and positive definite
  
         int signum;
         gsl_linalg_LU_decomp(S, permS, &signum);
         int inverr = gsl_linalg_LU_invert(S, permS, Sinv);
  
         if (inverr != 0) {
           B2ERROR("S: gsl_linalg_LU_invert error " << inverr << " in error calculation");
           ierr = -1;
           return -1;
         }
  
  
 // *-- Calculate G.
 //  (same as W1, but for measured parameters)
 // G = Feta^T * Sinv * Feta
  
         // SinvFeta[ncon][nmea] = 1*Sinv[ncon][ncon]*Feta[ncon][nmea] + 0*SinvFeta
         gsl_blas_dsymm(CblasLeft, CblasUpper, 1, Sinv, &Feta.matrix, 0,  SinvFeta);
         // G[nmea][nmea] = 1*Feta^T[nmea][ncon]*SinvFeta[ncon][nmea] + 0*G
         gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1, &Feta.matrix, SinvFeta, 0, G);
  
         if (debug > 2) debug_print(G, "G(1)");
  
  
         if (nunm > 0) {
  
 // *-- Calculate H
  
           // Fxi is the part of Fetaxi containing the unmeasured quantities, if any
           gsl_matrix_view Fxi = gsl_matrix_submatrix(Fetaxi,  0, nmea, ncon, nunm);
           // H = Feta^T * Sinv * Fxi
           // SinvFxi[ncon][nunm] = 1*Sinv[ncon][ncon]*Fxi[ncon][nunm] + 0*SinvFxi
           gsl_blas_dsymm(CblasLeft, CblasUpper, 1, Sinv, &Fxi.matrix, 0,  SinvFxi);
           // H[nmea][nunm] = 1*Feta^T[nmea][ncon]*SinvFxi[ncon][nunm] + 0*H
           gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1, &Feta.matrix, SinvFxi, 0, H);
  
           if (debug > 2) debug_print(H, "H");
  
 // *-- Calculate U**-1 and invert.
 //   (same as W1)
 //   U is a part of Minv
  
           gsl_matrix* Uinv = W1;
           gsl_matrix_view U = gsl_matrix_submatrix(Minv, nmea, nmea, nunm, nunm);
           // Uinv = Fxi^T * Sinv * Fxi
           // Uinv[nunm][nunm] = 1*Fxi^T[nunm][ncon]*SinvFxi[ncon][nunm] + 0*W1
           gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1, &Fxi.matrix, SinvFxi, 0, Uinv);
  
           gsl_linalg_LU_decomp(Uinv, permU, &signum);
           inverr = gsl_linalg_LU_invert(Uinv, permU, &U.matrix);
  
           if (debug > 2) debug_print(&U.matrix, "U");
           for (int i = 0; i < npar; ++i) {
             for (int j = 0; j < npar; ++j) {
               B2DEBUG(12, "after U Minv[" << i << "," << j << "]=" << gsl_matrix_get(Minv, i, j));
             }
           }
  
           if (inverr != 0) {
             B2ERROR("U: gsl_linalg_LU_invert error " << inverr << " in error calculation ");
             ierr = -1;
             return -1;
           }
  
 // *-- Covariance matrix between measured and unmeasured parameters.
  
 //    HU[nmea][nunm] = 1*H[nmea][nunm]*U[nunm][nunm] + 0*HU
           gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, H, &U.matrix, 0, HU);
 //    Vnewetaxi is a view of Vnew
           gsl_matrix_view Minvetaxi = gsl_matrix_submatrix(Minv, 0, nmea, nmea, nunm);
           gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, -1, &Vetaeta.matrix, HU, 0, &Minvetaxi.matrix);
           for (int i = 0; i < npar; ++i) {
             for (int j = 0; j < npar; ++j) {
               B2DEBUG(12, "after etaxi Minv[" << i << "," << j << "]=" << gsl_matrix_get(Minv, i, j));
             }
           }
  
  
 // *-- Fill in symmetric part:
 //    Vnewxieta is a view of Vnew
           gsl_matrix_view Minvxieta = gsl_matrix_submatrix(Minv, nmea, 0, nunm, nmea);
           gsl_matrix_transpose_memcpy(&Minvxieta.matrix, &Minvetaxi.matrix);
           for (int i = 0; i < npar; ++i) {
             for (int j = 0; j < npar; ++j) {
               B2DEBUG(12, "after symmetric: Minv[" << i << "," << j << "]=" << gsl_matrix_get(Minv, i, j));
             }
           }
  
 // *-- Calculate G-HUH^T:
 //    G = -1*HU*H^T +1*G
           gsl_blas_dgemm(CblasNoTrans, CblasTrans, -1, HU, H, +1, G);
  
         }  // endif nunm > 0
  
 // *-- Calculate I-GV.
         // IGV = 1
         gsl_matrix_set_identity(IGV);
         // IGV = -1*G*Vetaeta + 1*IGV
         gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, -1, G, &Vetaeta.matrix, 1, IGV);
  
 // *-- And finally error matrix on fitted parameters.
         gsl_matrix_view Minvetaeta = gsl_matrix_submatrix(Minv, 0, 0, nmea, nmea);
  
         // Vnewetaeta = 1*Vetaeta*IGV + 0*Vnewetaeta
         gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, &Vetaeta.matrix, IGV, 0, &Minvetaeta.matrix);
  
         for (int i = 0; i < npar; ++i) {
           for (int j = 0; j < npar; ++j) {
             B2DEBUG(12, "complete Minv[" << i << "," << j << "]=" << gsl_matrix_get(Minv, i, j));
           }
         }
  
 // *-- now Minv should be complete, can calculate new covariance matrix Vnew
 //       Vnew = (dx / dt)^T  * V_t * dx / dt
 //       x are the fitted parameters, t are the measured parameters,
 //       V_t is the covariance matrix of the measured parameters
 //       thus in OPALFitter variables:  V_t = Vetaeta,  x = (eta, xi) = etaxi
 //       d eta / dt and d xi / dt are given by eqns 9.54 and 9.55, respectively,
 //       as deta/dt = - Minvetaeta * Fetat and dxi/dt = - Minvxieta * Fetat
 //       Fetat is d^2 chi^2 / deta dt = - V^-1, even in presence of soft constraints, since
 //       the soft constraints don't depend on the _measured_ parameters t (but on the fitted ones)
 //       HERE,  V (and Vinv) do not include the second derivatives of the soft constraints
 //       so we can use them right away
 //       Note that "F" is used for completely different things:
 //       in OPALFitter it is the derivatives of the constraints wrt to the parameters (A in book)
 //       in the book, F are the second derivatives of the objective function wrt the parameters
  
  
         gsl_matrix_view detadt = gsl_matrix_submatrix(dxdt, 0, 0, nmea, nmea);
         B2DEBUG(13, "after detadt");
         //  Vdxdt is Vetaeta * dxdt^T, thus Vdxdt[nmea][npar]
         gsl_matrix_view Vdetadt = gsl_matrix_submatrix(Vdxdt, 0, 0, nmea, nmea);
         B2DEBUG(13, "after Vdetadt");
  
         // detadt = - Minvetaeta * Fetat = -1 * Minvetaeta * (-1) * Vinv + 0 * detadt   // replace by symm?
         gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, &Minvetaeta.matrix, Vinv, 0, &detadt.matrix);
         if (debug > 2) debug_print(&detadt.matrix, "deta/dt");
  
         // Vdetadt = 1 * Vetaeta * detadt^T + 0* Vdetadt
         gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1, &Vetaeta.matrix, &detadt.matrix, 0, &Vdetadt.matrix);   // ok
         if (debug > 2) debug_print(&Vdetadt.matrix, "Vetata * deta/dt");
  
         gsl_matrix_view Vnewetaeta = gsl_matrix_submatrix(Vnew, 0, 0, nmea, nmea);    //[nmea],[nmea]
         // Vnewetaeta = 1 * detadt * Vdetadt + 0* Vnewetaeta
         gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, &detadt.matrix, &Vdetadt.matrix, 0, &Vnewetaeta.matrix);
  
         if (debug > 2) debug_print(Vnew, "Vnew after part for measured parameters");
  
  
         if (nunm > 0) {
  
           gsl_matrix_view Minvxieta = gsl_matrix_submatrix(Minv, nmea, 0, nunm, nmea);   //[nunm][nmea]
           B2DEBUG(13, "after Minvxieta");
           if (debug > 2) debug_print(&Minvxieta.matrix, "Minvxieta");
  
           gsl_matrix_view dxidt = gsl_matrix_submatrix(dxdt, nmea, 0, nunm, nmea);       //[nunm][nmea]
           B2DEBUG(13, "after dxidt");
           // dxidt[nunm][nmea] = - Minvxieta * Fetat = -1 * Minvxieta[nunm][nmea] * Vinv[nmea][nmea] + 0 * dxidt
           gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, &Minvxieta.matrix, Vinv, 0, &dxidt.matrix);    //ok
           if (debug > 2) debug_print(&dxidt.matrix, "dxi/dt");
  
           // Vdxdt = V * dxdt^T => Vdxdt[nmea][npar]
           gsl_matrix_view Vdxidt = gsl_matrix_submatrix(Vdxdt, 0, nmea, nmea, nunm);     //[nmea][nunm]
           B2DEBUG(13, "after Vdxidt");
           // Vdxidt = 1 * Vetaeta[nmea][nmea] * dxidt^T[nmea][nunm] + 0* Vdxidt => Vdxidt[nmea][nunm]
           gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1, &Vetaeta.matrix, &dxidt.matrix, 0, &Vdxidt.matrix);   // ok
           if (debug > 2) debug_print(&Vdxidt.matrix, "Vetaeta * dxi/dt^T");
  
           gsl_matrix_view Vnewetaxi = gsl_matrix_submatrix(Vnew, 0, nmea, nmea, nunm);     //[nmea][nunm]
           gsl_matrix_view Vnewxieta = gsl_matrix_submatrix(Vnew, nmea, 0, nunm, nmea);     //[nunm][nmea]
           gsl_matrix_view Vnewxixi = gsl_matrix_submatrix(Vnew, nmea, nmea, nunm, nunm);   //[nunm][nunm]
  
           // Vnewxieta[nunm][nmea] = 1 * dxidt[nunm][nmea] * Vdetadt[nmea][nmea] + 0* Vnewxieta
           gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, &dxidt.matrix, &Vdetadt.matrix, 0, &Vnewxieta.matrix);   // ok
           if (debug > 2) debug_print(Vnew, "Vnew after xieta part");
           // Vnewetaxi[nmea][nunm] = 1 * detadt[nmea][nmea] * Vdxidt[nmea][nunm] + 0* Vnewetaxi
           gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, &detadt.matrix, &Vdxidt.matrix, 0, &Vnewetaxi.matrix);   // ok
           if (debug > 2) debug_print(Vnew, "Vnew after etaxi part");
           // Vnewxixi[nunm][nunm] = 1 * dxidt[nunm][nmea] * Vdxidt[nmea][nunm] + 0* Vnewxixi
           gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, &dxidt.matrix, &Vdxidt.matrix, 0, &Vnewxixi.matrix);
           if (debug > 2) debug_print(Vnew, "Vnew after xixi part");
         }
  
 // *-- now we finally have Vnew, fill into fitobjects
  
         // update errors in fitobjects
         for (auto& fitobject : fitobjects) {
           for (int ilocal = 0; ilocal < fitobject->getNPar(); ++ilocal) {
             int iglobal = fitobject->getGlobalParNum(ilocal);
             for (int jlocal = ilocal; jlocal < fitobject->getNPar(); ++jlocal) {
               int jglobal = fitobject->getGlobalParNum(jlocal);
               if (iglobal >= 0 && jglobal >= 0)
                 fitobject->setCov(ilocal, jlocal, gsl_matrix_get(Vnew, iglobal, jglobal));
             }
           }
         }
  
         // Finally, copy covariance matrix
         if (cov && covDim != nmea + nunm) {
           delete[] cov;
           cov = nullptr;
         }
         covDim = nmea + nunm;
         if (!cov) cov = new double[covDim * covDim];
         for (int i = 0; i < covDim; ++i) {
           for (int j = 0; j < covDim; ++j) {
             cov[i * covDim + j] = gsl_matrix_get(Vnew, i, j);
           }
         }
         covValid = true;
  
       } // endif calcerr == true
  
 #ifndef FIT_TRACEOFF
       if (tracer) tracer->finish(*this);
 #endif
  
       return fitprob;
  
     }

◆ getError()

int getError ( ) const

overridevirtual

Return error code.

Error code meanings:

0: No error
1: out of iterations
2: crazy chi^2
3: minimum step size reached, chiK still increasing
4: (step size decreased)
5: (keep going) Get the error code of the last fit: 0=OK, 1=failed

Implements BaseFitter.

Definition at line 1012 of file OPALFitterGSL.cc.

◆ getGlobalCovarianceMatrix() [1/2]

double * getGlobalCovarianceMatrix ( int & idim )

virtualinherited

Parameters

idim	1st dimension of global covariance matrix

Definition at line 158 of file BaseFitter.cc.

◆ getGlobalCovarianceMatrix() [2/2]

const double * getGlobalCovarianceMatrix ( int & idim ) const

virtualinherited

Parameters

idim	1st dimension of global covariance matrix

Definition at line 148 of file BaseFitter.cc.

The documentation for this class was generated from the following files:

analysis/OrcaKinFit/include/OPALFitterGSL.h
analysis/OrcaKinFit/src/OPALFitterGSL.cc

Public Member Functions

Public Attributes

Protected Types

Protected Member Functions

Static Protected Member Functions