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

#include <OPALFitterGSL.h>

Inheritance diagram for OPALFitterGSL:

Public Member Functions
virtual double	fit () override

virtual int	getError () const override
	Return error code.

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 Typedef Documentation

◆ ConstraintContainer

typedef std::vector<BaseHardConstraint*> ConstraintContainer

protectedinherited

Definition at line 92 of file BaseFitter.h.

◆ ConstraintIterator

typedef ConstraintContainer::iterator ConstraintIterator

protectedinherited

Definition at line 96 of file BaseFitter.h.

◆ FitObjectContainer

typedef std::vector<BaseFitObject*> FitObjectContainer

protectedinherited

Definition at line 91 of file BaseFitter.h.

◆ FitObjectIterator

typedef FitObjectContainer::iterator FitObjectIterator

protectedinherited

Definition at line 95 of file BaseFitter.h.

◆ SoftConstraintContainer

typedef std::vector<BaseSoftConstraint*> SoftConstraintContainer

protectedinherited

Definition at line 93 of file BaseFitter.h.

◆ SoftConstraintIterator

typedef SoftConstraintContainer::iterator SoftConstraintIterator

protectedinherited

Definition at line 97 of file BaseFitter.h.

Member Enumeration Documentation

◆ anonymous enum

anonymous enum

protected

Definition at line 136 of file OPALFitterGSL.h.

136{NPARMAX = 50, NCONMAX = 20, NUNMMAX = 20};

Constructor & Destructor Documentation

◆ OPALFitterGSL()

OPALFitterGSL ( )

Definition at line 50 of file OPALFitterGSL.cc.

      : npar(0), nmea(0), nunm(0), ncon(0), ierr(0), nit(0),
        fitprob(0), chi2(0),
        f(nullptr), r(nullptr), Fetaxi(nullptr), S(nullptr), Sinv(nullptr), SinvFxi(nullptr), SinvFeta(nullptr),
        W1(nullptr), G(nullptr), H(nullptr), HU(nullptr), IGV(nullptr), V(nullptr), VLU(nullptr), Vinv(nullptr), Vnew(nullptr),
        Minv(nullptr), dxdt(nullptr), Vdxdt(nullptr),
        dxi(nullptr), Fxidxi(nullptr), lambda(nullptr), FetaTlambda(nullptr),
        etaxi(nullptr), etasv(nullptr), y(nullptr), y_eta(nullptr), Vinvy_eta(nullptr), FetaV(nullptr),
        permS(nullptr), permU(nullptr), permV(nullptr), debug(0)
    {}

◆ ~OPALFitterGSL()

~OPALFitterGSL ( )

virtual

Definition at line 62 of file OPALFitterGSL.cc.

    {
      if (f) gsl_vector_free(f);
      if (r) gsl_vector_free(r);
      if (Fetaxi) gsl_matrix_free(Fetaxi);
      if (S) gsl_matrix_free(S);
      if (Sinv) gsl_matrix_free(Sinv);
      if (SinvFxi) gsl_matrix_free(SinvFxi);
      if (SinvFeta) gsl_matrix_free(SinvFeta);
      if (W1) gsl_matrix_free(W1);
      if (G) gsl_matrix_free(G);
      if (H) gsl_matrix_free(H);
      if (HU) gsl_matrix_free(HU);
      if (IGV) gsl_matrix_free(IGV);
      if (V) gsl_matrix_free(V);
      if (VLU) gsl_matrix_free(VLU);
      if (Vinv) gsl_matrix_free(Vinv);
      if (Vnew) gsl_matrix_free(Vnew);
      if (Minv) gsl_matrix_free(Minv);
      if (dxdt) gsl_matrix_free(dxdt);
      if (Vdxdt) gsl_matrix_free(Vdxdt);
      if (dxi) gsl_vector_free(dxi);
      if (Fxidxi) gsl_vector_free(Fxidxi);
      if (lambda) gsl_vector_free(lambda);
      if (FetaTlambda) gsl_vector_free(FetaTlambda);
      if (etaxi) gsl_vector_free(etaxi);
      if (etasv) gsl_vector_free(etasv);
      if (y) gsl_vector_free(y);
      if (y_eta) gsl_vector_free(y_eta);
      if (Vinvy_eta) gsl_vector_free(Vinvy_eta);
      if (FetaV) gsl_matrix_free(FetaV);
      if (permS) gsl_permutation_free(permS);
      if (permU) gsl_permutation_free(permU);
      if (permV) gsl_permutation_free(permV);
 
      f = nullptr;
      r = nullptr;
      Fetaxi = nullptr;
      S = nullptr;
      Sinv = nullptr;
      SinvFxi = nullptr;
      SinvFeta = nullptr;
      W1 = nullptr;
      G = nullptr;
      H = nullptr;
      HU = nullptr;
      IGV = nullptr;
      V = nullptr;
      VLU = nullptr;
      Vinv = nullptr;
      Vnew = nullptr;
      Minv = nullptr;
      dxdt = nullptr;
      Vdxdt = nullptr;
      dxi = nullptr;
      Fxidxi = nullptr;
      lambda = nullptr;
      FetaTlambda = nullptr;
      etaxi = nullptr;
      etasv = nullptr;
      y = nullptr;
      y_eta = nullptr;
      Vinvy_eta = nullptr;
      FetaV = nullptr;
      permS = nullptr;
      permU = nullptr;
      permV = nullptr;
    }

Member Function Documentation

◆ addConstraint() [1/2]

void addConstraint ( BaseConstraint & constraint_ )

virtualinherited

Definition at line 74 of file BaseFitter.cc.

    {
      covValid = false;
      if (auto* hc = dynamic_cast<BaseHardConstraint*>(&constraint_))
        constraints.push_back(hc);
      else if (auto* sc = dynamic_cast<BaseSoftConstraint*>(&constraint_))
        softconstraints.push_back(sc);
    }

◆ addConstraint() [2/2]

void addConstraint ( BaseConstraint * constraint_ )

virtualinherited

Definition at line 60 of file BaseFitter.cc.

    {
      covValid = false;
 
      if (auto* hc = dynamic_cast<BaseHardConstraint*>(constraint_))
        constraints.push_back(hc);
      else if (auto* sc = dynamic_cast<BaseSoftConstraint*>(constraint_))
        softconstraints.push_back(sc);
      else {
        // illegal constraint
        assert(0);
      }
    }

◆ addFitObject() [1/2]

void addFitObject ( BaseFitObject & fitobject_ )

virtualinherited

Definition at line 54 of file BaseFitter.cc.

    {
      covValid = false;
      fitobjects.push_back(&fitobject_);
    }

◆ addFitObject() [2/2]

void addFitObject ( BaseFitObject * fitobject_ )

virtualinherited

Definition at line 48 of file BaseFitter.cc.

    {
      covValid = false;
      fitobjects.push_back(fitobject_);
    }

◆ addHardConstraint() [1/2]

void addHardConstraint ( BaseHardConstraint & constraint_ )

virtualinherited

Definition at line 89 of file BaseFitter.cc.

    {
      covValid = false;
      constraints.push_back(&constraint_);
    }

◆ addHardConstraint() [2/2]

void addHardConstraint ( BaseHardConstraint * constraint_ )

virtualinherited

Definition at line 83 of file BaseFitter.cc.

    {
      covValid = false;
      constraints.push_back(constraint_);
    }

◆ addSoftConstraint() [1/2]

void addSoftConstraint ( BaseSoftConstraint & constraint_ )

virtualinherited

Definition at line 101 of file BaseFitter.cc.

    {
      covValid = false;
      softconstraints.push_back(&constraint_);
    }

◆ addSoftConstraint() [2/2]

void addSoftConstraint ( BaseSoftConstraint * constraint_ )

virtualinherited

Definition at line 95 of file BaseFitter.cc.

    {
      covValid = false;
      softconstraints.push_back(constraint_);
    }

◆ debug_print() [1/2]

void debug_print	(	gsl_matrix *	m,
		const char *	name
	)

staticprotected

Definition at line 1057 of file OPALFitterGSL.cc.

    {
      for (unsigned int  i = 0; i < m->size1; ++i)
        for (unsigned int j = 0; j < m->size2; ++j)
          if (gsl_matrix_get(m, i, j) != 0)
            B2INFO(name << "[" << i << "][" << j << "]=" << gsl_matrix_get(m, i, j));
    }

◆ debug_print() [2/2]

void debug_print	(	gsl_vector *	v,
		const char *	name
	)

staticprotected

Definition at line 1065 of file OPALFitterGSL.cc.

    {
      for (unsigned int  i = 0; i < v->size; ++i)
        if (gsl_vector_get(v, i) != 0)
          B2INFO(name << "[" << i << "]=" << gsl_vector_get(v, i));
    }

◆ 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;
 
    }

◆ getChi2()

double getChi2 ( ) const

overridevirtual

Get the chi**2 of the last fit.

Implements BaseFitter.

Definition at line 1014 of file OPALFitterGSL.cc.

1014{return chi2;}

◆ getConstraints()

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

virtualinherited

Definition at line 112 of file BaseFitter.cc.

    {
      return &constraints;
    }

◆ getDoF()

int getDoF ( ) const

overridevirtual

Get the number of degrees of freedom of the last fit.

Implements BaseFitter.

Definition at line 1015 of file OPALFitterGSL.cc.

1015{return ncon - nunm;}

◆ 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.

1012{return ierr;}

◆ getFitObjects()

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

virtualinherited

Definition at line 107 of file BaseFitter.cc.

    {
      return &fitobjects;
    }

◆ getGlobalCovarianceMatrix() [1/2]

double * getGlobalCovarianceMatrix ( int & idim )

virtualinherited

Parameters

idim	1st dimension of global covariance matrix

Definition at line 158 of file BaseFitter.cc.

    {
      if (covValid && cov) {
        idim = covDim;
        return cov;
      }
      idim = 0;
      return nullptr;
    }

◆ 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.

    {
      if (covValid && cov) {
        idim = covDim;
        return cov;
      }
      idim = 0;
      return nullptr;
    }

◆ getIterations()

int getIterations ( ) const

overridevirtual

Get the number of iterations of the last fit.

Implements BaseFitter.

Definition at line 1016 of file OPALFitterGSL.cc.

1016{return nit;}

◆ getNcon()

int getNcon ( ) const

virtual

Get the number of hard constraints of the last fit.

Definition at line 1072 of file OPALFitterGSL.cc.

1072{return ncon;}

◆ getNpar()

int getNpar ( ) const

virtual

Get the number of all parameters of the last fit.

Definition at line 1075 of file OPALFitterGSL.cc.

1075{return npar;}

◆ getNsoft()

int getNsoft ( ) const

virtual

Get the number of soft constraints of the last fit.

Definition at line 1073 of file OPALFitterGSL.cc.

1073{return 0;}

◆ getNunm()

int getNunm ( ) const

virtual

Get the number of unmeasured parameters of the last fit.

Definition at line 1074 of file OPALFitterGSL.cc.

1074{return nunm;}

◆ getProbability()

double getProbability ( ) const

overridevirtual

Get the fit probability of the last fit.

Implements BaseFitter.

Definition at line 1013 of file OPALFitterGSL.cc.

1013{return fitprob;}

◆ getSoftConstraints()

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

virtualinherited

Definition at line 117 of file BaseFitter.cc.

    {
      return &softconstraints;
    }

◆ getTracer() [1/2]

BaseTracer * getTracer ( )

virtualinherited

Definition at line 130 of file BaseFitter.cc.

    {
      return tracer;
    }

◆ getTracer() [2/2]

const BaseTracer * getTracer ( ) const

virtualinherited

Definition at line 134 of file BaseFitter.cc.

    {
      return tracer;
    }

◆ ini_gsl_matrix()

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

staticprotected

Definition at line 1041 of file OPALFitterGSL.cc.

    {
      if (m) {
        if (m->size1 != size1 || m->size2 != size2) {
          gsl_matrix_free(m);
          if (size1 * size2 > 0) m = gsl_matrix_alloc(size1, size2);
          else m = nullptr;
        }
      } else if (size1 * size2 > 0) m = gsl_matrix_alloc(size1, size2);
    }

◆ ini_gsl_permutation()

void ini_gsl_permutation	(	gsl_permutation *&	p,
		unsigned int	size
	)

staticprotected

Definition at line 1018 of file OPALFitterGSL.cc.

    {
      if (p) {
        if (p->size != size) {
          gsl_permutation_free(p);
          if (size > 0) p = gsl_permutation_alloc(size);
          else p = nullptr;
        }
      } else if (size > 0) p = gsl_permutation_alloc(size);
    }

◆ ini_gsl_vector()

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

staticprotected

Definition at line 1029 of file OPALFitterGSL.cc.

    {
 
      if (v) {
        if (v->size != size) {
          gsl_vector_free(v);
          if (size > 0) v = gsl_vector_alloc(size);
          else v = nullptr;
        }
      } else if (size > 0) v = gsl_vector_alloc(size);
    }

◆ initialize()

bool initialize ( )

overridevirtual

Implements BaseFitter.

Definition at line 884 of file OPALFitterGSL.cc.

    {
      covValid = false;
      // tell fitobjects the global ordering of their parameters:
      int iglobal = 0;
      // measured parameters first!
      for (auto& fitobject : fitobjects) {
        for (int ilocal = 0; ilocal < fitobject->getNPar(); ++ilocal) {
          if (fitobject->isParamMeasured(ilocal) &&
              !fitobject->isParamFixed(ilocal)) {
            fitobject->setGlobalParNum(ilocal, iglobal);
            B2DEBUG(10, "Object " << fitobject->getName()
                    << " Parameter " << fitobject->getParamName(ilocal)
                    << " is measured, global number " << iglobal);
            ++iglobal;
          }
        }
      }
      nmea = iglobal;
      // now  unmeasured parameters!
      for (auto& fitobject : fitobjects) {
        for (int ilocal = 0; ilocal < fitobject->getNPar(); ++ilocal) {
          if (!fitobject->isParamMeasured(ilocal) &&
              !fitobject->isParamFixed(ilocal)) {
            fitobject->setGlobalParNum(ilocal, iglobal);
            B2DEBUG(10, "Object " << fitobject->getName()
                    << " Parameter " << fitobject->getParamName(ilocal)
                    << " is unmeasured, global number " << iglobal);
            ++iglobal;
          }
        }
      }
      npar = iglobal;
      assert(npar <= NPARMAX);
      nunm = npar - nmea;
      assert(nunm <= NUNMMAX);
 
      // set number of constraints
      ncon = constraints.size();
      assert(ncon <= NCONMAX);
 
      ini_gsl_vector(f, ncon);
      ini_gsl_vector(r, ncon);
 
      ini_gsl_matrix(Fetaxi, ncon, npar);
      ini_gsl_matrix(S, ncon, ncon);
      ini_gsl_matrix(Sinv, ncon, ncon);
      ini_gsl_matrix(SinvFxi, ncon, nunm);
      ini_gsl_matrix(SinvFeta, ncon, nmea);
      ini_gsl_matrix(W1, nunm, nunm);
      ini_gsl_matrix(G, nmea, nmea);
      ini_gsl_matrix(H, nmea, nunm);
      ini_gsl_matrix(HU, nmea, nunm);
      ini_gsl_matrix(IGV, nmea, nmea);
      ini_gsl_matrix(V, npar, npar);
      ini_gsl_matrix(VLU, nmea, nmea);
      ini_gsl_matrix(Vinv, nmea, nmea);
      ini_gsl_matrix(Vnew, npar, npar);
      ini_gsl_matrix(Minv, npar, npar);
      ini_gsl_matrix(dxdt, npar, nmea);
      ini_gsl_matrix(Vdxdt, nmea, npar);
 
      ini_gsl_vector(dxi, nunm);
      ini_gsl_vector(Fxidxi, ncon);
      ini_gsl_vector(lambda, ncon);
      ini_gsl_vector(FetaTlambda, nmea);
 
      ini_gsl_vector(etaxi, npar);
      ini_gsl_vector(etasv, npar);
      ini_gsl_vector(y, nmea);
      ini_gsl_vector(y_eta, nmea);
      ini_gsl_vector(Vinvy_eta, nmea);
 
      ini_gsl_matrix(FetaV, ncon, nmea);
 
      ini_gsl_permutation(permS, ncon);
      ini_gsl_permutation(permU, nunm);
      ini_gsl_permutation(permV, nmea);
 
      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(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(permS && (int)permS->size == ncon);
      assert(nunm == 0 || (permU && (int)permU->size == nunm));
      assert(permV && (int)permV->size == nmea);
 
      return true;
 
    }

◆ reset()

void reset ( )

virtualinherited

Definition at line 122 of file BaseFitter.cc.

    {
      fitobjects.resize(0);
      constraints.resize(0);
      softconstraints.resize(0);
      covValid = false;
    }

◆ setDebug()

void setDebug ( int debuglevel )

virtual

Set the Debug Level.

Definition at line 1052 of file OPALFitterGSL.cc.

    {
      debug = debuglevel;
    }

◆ setTracer() [1/2]

void setTracer ( BaseTracer & newTracer )

virtualinherited

Definition at line 143 of file BaseFitter.cc.

    {
      tracer = &newTracer;
    }

◆ setTracer() [2/2]

void setTracer ( BaseTracer * newTracer )

virtualinherited

Definition at line 138 of file BaseFitter.cc.

    {
      tracer = newTracer;
    }

◆ updateFitObjects()

bool updateFitObjects ( double etaxi[] )

protectedvirtual

Definition at line 998 of file OPALFitterGSL.cc.

    {
      // changed etaxi -> eetaxi to avoid clash with class member DJeans
      //bool debug = false;
      bool result = true;
      for (auto& fitobject : fitobjects) {
        for (int ilocal = 0; ilocal < fitobject->getNPar(); ++ilocal) {
          fitobject->updateParams(eetaxi, npar);
//       }
        }
      }
      return result;
    }

Member Data Documentation

◆ chi2

double chi2

protected

Definition at line 146 of file OPALFitterGSL.h.

◆ constraints

ConstraintContainer constraints

protectedinherited

Definition at line 100 of file BaseFitter.h.

◆ cov

double* cov

protectedinherited

global covariance matrix of last fit problem

Definition at line 104 of file BaseFitter.h.

◆ covDim

int covDim

protectedinherited

dimension of global covariance matrix

Definition at line 103 of file BaseFitter.h.

◆ covValid

bool covValid

protectedinherited

Flag whether global covariance is valid.

Definition at line 105 of file BaseFitter.h.

◆ debug

int debug

private

Definition at line 191 of file OPALFitterGSL.h.

◆ dxdt

gsl_matrix* dxdt

private

Definition at line 173 of file OPALFitterGSL.h.

◆ dxi

gsl_vector* dxi

private

Definition at line 175 of file OPALFitterGSL.h.

◆ etasv

gsl_vector* etasv

private

Definition at line 180 of file OPALFitterGSL.h.

◆ etaxi

gsl_vector* etaxi

private

Definition at line 179 of file OPALFitterGSL.h.

◆ f

gsl_vector* f

private

Definition at line 156 of file OPALFitterGSL.h.

◆ FetaTlambda

gsl_vector* FetaTlambda

private

Definition at line 178 of file OPALFitterGSL.h.

◆ FetaV

gsl_matrix* FetaV

private

Definition at line 185 of file OPALFitterGSL.h.

◆ Fetaxi

gsl_matrix* Fetaxi

private

Definition at line 158 of file OPALFitterGSL.h.

◆ fitobjects

FitObjectContainer fitobjects

protectedinherited

Definition at line 99 of file BaseFitter.h.

◆ fitprob

double fitprob

protected

Definition at line 145 of file OPALFitterGSL.h.

◆ Fxidxi

gsl_vector* Fxidxi

private

Definition at line 176 of file OPALFitterGSL.h.

◆ G

gsl_matrix* G

private

Definition at line 164 of file OPALFitterGSL.h.

◆ H

gsl_matrix* H

private

Definition at line 165 of file OPALFitterGSL.h.

◆ HU

gsl_matrix* HU

private

Definition at line 166 of file OPALFitterGSL.h.

◆ ierr

int ierr

protected

Definition at line 142 of file OPALFitterGSL.h.

◆ IGV

gsl_matrix* IGV

private

Definition at line 167 of file OPALFitterGSL.h.

◆ lambda

gsl_vector* lambda

private

Definition at line 177 of file OPALFitterGSL.h.

◆ Minv

gsl_matrix* Minv

private

Definition at line 172 of file OPALFitterGSL.h.

◆ ncon

int ncon

protected

Definition at line 141 of file OPALFitterGSL.h.

◆ nit

int nit

protected

Definition at line 143 of file OPALFitterGSL.h.

◆ nmea

int nmea

protected

Definition at line 139 of file OPALFitterGSL.h.

◆ npar

int npar

protected

Definition at line 138 of file OPALFitterGSL.h.

◆ nunm

int nunm

protected

Definition at line 140 of file OPALFitterGSL.h.

◆ permS

gsl_permutation* permS

private

Helper permutation vector.

Definition at line 187 of file OPALFitterGSL.h.

◆ permU

gsl_permutation* permU

private

Helper permutation vector.

Definition at line 188 of file OPALFitterGSL.h.

◆ permV

gsl_permutation* permV

private

Helper permutation vector.

Definition at line 189 of file OPALFitterGSL.h.

◆ r

gsl_vector* r

private

Definition at line 157 of file OPALFitterGSL.h.

◆ S

gsl_matrix* S

private

Definition at line 159 of file OPALFitterGSL.h.

◆ Sinv

gsl_matrix* Sinv

private

Definition at line 160 of file OPALFitterGSL.h.

◆ SinvFeta

gsl_matrix* SinvFeta

private

Definition at line 162 of file OPALFitterGSL.h.

◆ SinvFxi

gsl_matrix* SinvFxi

private

Definition at line 161 of file OPALFitterGSL.h.

◆ softconstraints

SoftConstraintContainer softconstraints

protectedinherited

Definition at line 101 of file BaseFitter.h.

◆ tracer

BaseTracer* tracer

protectedinherited

Definition at line 108 of file BaseFitter.h.

◆ traceValues

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

inherited

Definition at line 112 of file BaseFitter.h.

◆ V

gsl_matrix* V

private

Definition at line 168 of file OPALFitterGSL.h.

◆ Vdxdt

gsl_matrix* Vdxdt

private

Definition at line 174 of file OPALFitterGSL.h.

◆ Vinv

gsl_matrix* Vinv

private

Definition at line 170 of file OPALFitterGSL.h.

◆ Vinvy_eta

gsl_vector* Vinvy_eta

private

Definition at line 183 of file OPALFitterGSL.h.

◆ VLU

gsl_matrix* VLU

private

Definition at line 169 of file OPALFitterGSL.h.

◆ Vnew

gsl_matrix* Vnew

private

Definition at line 171 of file OPALFitterGSL.h.

◆ W1

gsl_matrix* W1

private

Definition at line 163 of file OPALFitterGSL.h.

◆ y

gsl_vector* y

private

Definition at line 181 of file OPALFitterGSL.h.

◆ y_eta

gsl_vector* y_eta

private

Definition at line 182 of file OPALFitterGSL.h.

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

Protected Attributes

Private Attributes

Detailed Description

Member Typedef Documentation

◆ ConstraintContainer

◆ ConstraintIterator

◆ FitObjectContainer

◆ FitObjectIterator

◆ SoftConstraintContainer

◆ SoftConstraintIterator

Member Enumeration Documentation

◆ anonymous enum

Constructor & Destructor Documentation

◆ OPALFitterGSL()

◆ ~OPALFitterGSL()

Member Function Documentation

◆ addConstraint() [1/2]

◆ addConstraint() [2/2]

◆ addFitObject() [1/2]

◆ addFitObject() [2/2]

◆ addHardConstraint() [1/2]

◆ addHardConstraint() [2/2]

◆ addSoftConstraint() [1/2]

◆ addSoftConstraint() [2/2]

◆ debug_print() [1/2]

◆ debug_print() [2/2]

◆ fit()

◆ getChi2()

◆ getConstraints()

◆ getDoF()

◆ getError()

◆ getFitObjects()

◆ getGlobalCovarianceMatrix() [1/2]

◆ getGlobalCovarianceMatrix() [2/2]

◆ getIterations()

◆ getNcon()

◆ getNpar()

◆ getNsoft()

◆ getNunm()

◆ getProbability()

◆ getSoftConstraints()

◆ getTracer() [1/2]

◆ getTracer() [2/2]

◆ ini_gsl_matrix()

◆ ini_gsl_permutation()

◆ ini_gsl_vector()

◆ initialize()

◆ reset()

◆ setDebug()

◆ setTracer() [1/2]

◆ setTracer() [2/2]

◆ updateFitObjects()

Member Data Documentation

◆ chi2

◆ constraints

◆ cov

◆ covDim

◆ covValid

◆ debug

◆ dxdt

◆ dxi

◆ etasv

◆ etaxi

◆ f

◆ FetaTlambda

◆ FetaV

◆ Fetaxi

◆ fitobjects

◆ fitprob

◆ Fxidxi

◆ G

◆ H

◆ HU

◆ ierr

◆ IGV