Belle II Software light-2406-ragdoll
OPALFitterGSL Class Reference

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

#include <OPALFitterGSL.h>

Inheritance diagram for OPALFitterGSL:
Collaboration 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 BaseTracergetTracer ()
 
virtual const BaseTracergetTracer () 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.
 
BaseTracertracer
 

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()

Definition at line 50 of file OPALFitterGSL.cc.

51 : npar(0), nmea(0), nunm(0), ncon(0), ierr(0), nit(0),
52 fitprob(0), chi2(0),
53 f(nullptr), r(nullptr), Fetaxi(nullptr), S(nullptr), Sinv(nullptr), SinvFxi(nullptr), SinvFeta(nullptr),
54 W1(nullptr), G(nullptr), H(nullptr), HU(nullptr), IGV(nullptr), V(nullptr), VLU(nullptr), Vinv(nullptr), Vnew(nullptr),
55 Minv(nullptr), dxdt(nullptr), Vdxdt(nullptr),
56 dxi(nullptr), Fxidxi(nullptr), lambda(nullptr), FetaTlambda(nullptr),
57 etaxi(nullptr), etasv(nullptr), y(nullptr), y_eta(nullptr), Vinvy_eta(nullptr), FetaV(nullptr),
58 permS(nullptr), permU(nullptr), permV(nullptr), debug(0)
59 {}
gsl_permutation * permU
Helper permutation vector.
gsl_permutation * permV
Helper permutation vector.
gsl_permutation * permS
Helper permutation vector.

◆ ~OPALFitterGSL()

~OPALFitterGSL ( )
virtual

Definition at line 62 of file OPALFitterGSL.cc.

63 {
64 if (f) gsl_vector_free(f);
65 if (r) gsl_vector_free(r);
66 if (Fetaxi) gsl_matrix_free(Fetaxi);
67 if (S) gsl_matrix_free(S);
68 if (Sinv) gsl_matrix_free(Sinv);
69 if (SinvFxi) gsl_matrix_free(SinvFxi);
70 if (SinvFeta) gsl_matrix_free(SinvFeta);
71 if (W1) gsl_matrix_free(W1);
72 if (G) gsl_matrix_free(G);
73 if (H) gsl_matrix_free(H);
74 if (HU) gsl_matrix_free(HU);
75 if (IGV) gsl_matrix_free(IGV);
76 if (V) gsl_matrix_free(V);
77 if (VLU) gsl_matrix_free(VLU);
78 if (Vinv) gsl_matrix_free(Vinv);
79 if (Vnew) gsl_matrix_free(Vnew);
80 if (Minv) gsl_matrix_free(Minv);
81 if (dxdt) gsl_matrix_free(dxdt);
82 if (Vdxdt) gsl_matrix_free(Vdxdt);
83 if (dxi) gsl_vector_free(dxi);
84 if (Fxidxi) gsl_vector_free(Fxidxi);
85 if (lambda) gsl_vector_free(lambda);
86 if (FetaTlambda) gsl_vector_free(FetaTlambda);
87 if (etaxi) gsl_vector_free(etaxi);
88 if (etasv) gsl_vector_free(etasv);
89 if (y) gsl_vector_free(y);
90 if (y_eta) gsl_vector_free(y_eta);
91 if (Vinvy_eta) gsl_vector_free(Vinvy_eta);
92 if (FetaV) gsl_matrix_free(FetaV);
93 if (permS) gsl_permutation_free(permS);
94 if (permU) gsl_permutation_free(permU);
95 if (permV) gsl_permutation_free(permV);
96
97 f = nullptr;
98 r = nullptr;
99 Fetaxi = nullptr;
100 S = nullptr;
101 Sinv = nullptr;
102 SinvFxi = nullptr;
103 SinvFeta = nullptr;
104 W1 = nullptr;
105 G = nullptr;
106 H = nullptr;
107 HU = nullptr;
108 IGV = nullptr;
109 V = nullptr;
110 VLU = nullptr;
111 Vinv = nullptr;
112 Vnew = nullptr;
113 Minv = nullptr;
114 dxdt = nullptr;
115 Vdxdt = nullptr;
116 dxi = nullptr;
117 Fxidxi = nullptr;
118 lambda = nullptr;
119 FetaTlambda = nullptr;
120 etaxi = nullptr;
121 etasv = nullptr;
122 y = nullptr;
123 y_eta = nullptr;
124 Vinvy_eta = nullptr;
125 FetaV = nullptr;
126 permS = nullptr;
127 permU = nullptr;
128 permV = nullptr;
129 }

Member Function Documentation

◆ addConstraint() [1/2]

void addConstraint ( BaseConstraint constraint_)
virtualinherited

Definition at line 74 of file BaseFitter.cc.

75 {
76 covValid = false;
77 if (auto* hc = dynamic_cast<BaseHardConstraint*>(&constraint_))
78 constraints.push_back(hc);
79 else if (auto* sc = dynamic_cast<BaseSoftConstraint*>(&constraint_))
80 softconstraints.push_back(sc);
81 }
bool covValid
Flag whether global covariance is valid.
Definition: BaseFitter.h:105

◆ addConstraint() [2/2]

void addConstraint ( BaseConstraint constraint_)
virtualinherited

Definition at line 60 of file BaseFitter.cc.

61 {
62 covValid = false;
63
64 if (auto* hc = dynamic_cast<BaseHardConstraint*>(constraint_))
65 constraints.push_back(hc);
66 else if (auto* sc = dynamic_cast<BaseSoftConstraint*>(constraint_))
67 softconstraints.push_back(sc);
68 else {
69 // illegal constraint
70 assert(0);
71 }
72 }

◆ addFitObject() [1/2]

void addFitObject ( BaseFitObject fitobject_)
virtualinherited

Definition at line 54 of file BaseFitter.cc.

55 {
56 covValid = false;
57 fitobjects.push_back(&fitobject_);
58 }

◆ addFitObject() [2/2]

void addFitObject ( BaseFitObject fitobject_)
virtualinherited

Definition at line 48 of file BaseFitter.cc.

49 {
50 covValid = false;
51 fitobjects.push_back(fitobject_);
52 }

◆ addHardConstraint() [1/2]

void addHardConstraint ( BaseHardConstraint constraint_)
virtualinherited

Definition at line 89 of file BaseFitter.cc.

90 {
91 covValid = false;
92 constraints.push_back(&constraint_);
93 }

◆ addHardConstraint() [2/2]

void addHardConstraint ( BaseHardConstraint constraint_)
virtualinherited

Definition at line 83 of file BaseFitter.cc.

84 {
85 covValid = false;
86 constraints.push_back(constraint_);
87 }

◆ addSoftConstraint() [1/2]

void addSoftConstraint ( BaseSoftConstraint constraint_)
virtualinherited

Definition at line 101 of file BaseFitter.cc.

102 {
103 covValid = false;
104 softconstraints.push_back(&constraint_);
105 }

◆ addSoftConstraint() [2/2]

void addSoftConstraint ( BaseSoftConstraint constraint_)
virtualinherited

Definition at line 95 of file BaseFitter.cc.

96 {
97 covValid = false;
98 softconstraints.push_back(constraint_);
99 }

◆ debug_print() [1/2]

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

Definition at line 1057 of file OPALFitterGSL.cc.

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

◆ debug_print() [2/2]

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

Definition at line 1065 of file OPALFitterGSL.cc.

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

◆ fit()

double fit ( )
overridevirtual

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

Implements BaseFitter.

Definition at line 132 of file OPALFitterGSL.cc.

133 {
134
135
136
137 //
138 // ( ) ^ ^
139 // ( eta ) nmea |
140 // ( ) v |
141 // etaxi = (-----) --- npar
142 // ( ) ^ |
143 // ( xi ) nunm |
144 // ( ) v v
145
146 // <- ncon ->
147 // ( ) ^ ^
148 // ( Feta^T ) nmea |
149 // ( ) v |
150 // Fetaxi^T = ( -------- ) --- npar
151 // ( ) ^ |
152 // ( Fxi^T ) nunm |
153 // ( ) v v
154 //
155 // Fetaxi are the derivatives of the constraints wrt the fitted parameters, thus A_theta/xi in book
156
157
158 //
159 // <- nmea ->|<- nunm ->
160 // ( | ) ^ ^
161 // ( Vetaeta | Vetaxi ) nmea |
162 // ( | ) v |
163 // V = (----------+----------) --- npar
164 // ( | ) ^ |
165 // ( Vxieta | Vxixi ) nunm |
166 // ( | ) v v
167 //
168 // V is the covariance matrix. Before the fit only Vetaeta is non-zero (covariance of measured parameters)
169
170 //
171 // <- nmea ->|<- nunm ->
172 // ( | ) ^
173 // dxdt^T = ( detadt | dxidt ) nmea
174 // ( | ) v
175 //
176 // dxdt are the partial derivates of the fitted parameters wrt the measured parameters,
177 //
178 // <- nmea ->|<- nunm ->
179 // ( | ) ^
180 // Vdxdt = ( Vdetadt | Vdxidt ) nmea
181 // ( | ) v
182 //
183 // Vdxdt is Vetaeta * dxdt^T, thus Vdxdt[nmea][npar]
184 // both needed for calculation of the full covariance matrix of the fitted parameters
185
186 // order parameters etc
187 initialize();
188
189 assert(f && (int)f->size == ncon);
190 assert(r && (int)r->size == ncon);
191 assert(Fetaxi && (int)Fetaxi->size1 == ncon && (int)Fetaxi->size2 == npar);
192 assert(S && (int)S->size1 == ncon && (int)S->size2 == ncon);
193 assert(Sinv && (int)Sinv->size1 == ncon && (int)Sinv->size2 == ncon);
194 assert(nunm == 0 || (SinvFxi && (int)SinvFxi->size1 == ncon && (int)SinvFxi->size2 == nunm));
195 assert(SinvFeta && (int)SinvFeta->size1 == ncon && (int)SinvFeta->size2 == nmea);
196 assert(nunm == 0 || (W1 && (int)W1->size1 == nunm && (int)W1->size2 == nunm));
197 assert(G && (int)G->size1 == nmea && (int)G->size2 == nmea);
198 assert(nunm == 0 || (H && (int)H->size1 == nmea && (int)H->size2 == nunm));
199 assert(nunm == 0 || (HU && (int)HU->size1 == nmea && (int)HU->size2 == nunm));
200 assert(IGV && (int)IGV->size1 == nmea && (int)IGV->size2 == nmea);
201 assert(V && (int)V->size1 == npar && (int)V->size2 == npar);
202 assert(VLU && (int)VLU->size1 == nmea && (int)VLU->size2 == nmea);
203 assert(Vinv && (int)Vinv->size1 == nmea && (int)Vinv->size2 == nmea);
204 assert(Vnew && (int)Vnew->size1 == npar && (int)Vnew->size2 == npar);
205 assert(Minv && (int)Minv->size1 == npar && (int)Minv->size2 == npar);
206 assert(dxdt && (int)dxdt->size1 == npar && (int)dxdt->size2 == nmea);
207 assert(Vdxdt && (int)Vdxdt->size1 == nmea && (int)Vdxdt->size2 == npar);
208 assert(nunm == 0 || (dxi && (int)dxi->size == nunm));
209 assert(nunm == 0 || (Fxidxi && (int)Fxidxi->size == ncon));
210 assert(lambda && (int)lambda->size == ncon);
211 assert(FetaTlambda && (int)FetaTlambda->size == nmea);
212 assert(etaxi && (int)etaxi->size == npar);
213 assert(etasv && (int)etasv->size == npar);
214 assert(y && (int)y->size == nmea);
215 assert(y_eta && (int)y_eta->size == nmea);
216 assert(Vinvy_eta && (int)Vinvy_eta->size == nmea);
217 assert(FetaV && (int)FetaV->size1 == ncon && (int)FetaV->size2 == nmea);
218 assert(permS && (int)permS->size == ncon);
219 assert(nunm == 0 || (permU && (int)permU->size == nunm));
220 assert(permV && (int)permV->size == nmea);
221
222 // eta is the part of etaxi containing the measured quantities
223 gsl_vector_view eta = gsl_vector_subvector(etaxi, 0, nmea);
224
225 // Feta is the part of Fetaxi containing the measured quantities
226
227 B2DEBUG(11, "==== " << ncon << " " << nmea);
228
229 gsl_matrix_view Feta = gsl_matrix_submatrix(Fetaxi, 0, 0, ncon, nmea);
230
231 gsl_matrix_view Vetaeta = gsl_matrix_submatrix(V, 0, 0, nmea, nmea);
232
233 for (unsigned int i = 0; i < fitobjects.size(); ++i) {
234 for (int ilocal = 0; ilocal < fitobjects[i]->getNPar(); ++ilocal) {
235 if (!fitobjects[i]->isParamFixed(ilocal)) {
236 int iglobal = fitobjects[i]->getGlobalParNum(ilocal);
237 assert(iglobal >= 0 && iglobal < npar);
238 gsl_vector_set(etaxi, iglobal, fitobjects[i]->getParam(ilocal));
239 if (fitobjects[i]->isParamMeasured(ilocal)) {
240 assert(iglobal < nmea);
241 gsl_vector_set(y, iglobal, fitobjects[i]->getMParam(ilocal));
242 }
243 B2DEBUG(10, "etaxi[" << iglobal << "] = " << gsl_vector_get(etaxi, iglobal)
244 << " for jet " << i << " and ilocal = " << ilocal);
245 }
246 }
247 }
248
250 gsl_matrix_set_zero(Fetaxi);
251 for (int k = 0; k < ncon; k++) {
252 constraints[k]->getDerivatives(Fetaxi->size2, Fetaxi->block->data + k * Fetaxi->tda);
253 if (debug > 1) for (int j = 0; j < npar; j++)
254 if (gsl_matrix_get(Fetaxi, k, j) != 0)
255 B2INFO("1: Fetaxi[" << k << "][" << j << "] = " << gsl_matrix_get(Fetaxi, k, j));
256 }
257
258 // chi2's, step size, # iterations
259 double chinew = 0, chit = 0, chik = 0;
260 double alph = 1.;
261 nit = 0;
262 // convergence criteria (as in WWINIT)
263 int nitmax = 200;
264 double chik0 = 100.;
265 double chit0 = 100.;
266 double dchikc = 1.0E-3;
267 double dchitc = 1.0E-4;
268 double dchikt = 1.0E-2;
269 double dchik = 1.05;
270 double chimxw = 10000.;
271 double almin = 0.05;
272
273 // repeat with or with out smaller steps size
274 bool repeat = true;
275 bool scut = false;
276 bool calcerr = true;
277
278#ifndef FIT_TRACEOFF
279 if (tracer) tracer->initialize(*this);
280#endif
281
282
283 // start of iterations
284 while (repeat) {
285 bool updatesuccess = true;
286
287//*-- If necessary, retry smaller step, same direction
288 if (scut) {
289 gsl_vector_memcpy(etaxi, etasv);
290 updatesuccess = updateFitObjects(etaxi->block->data);
291 if (!updatesuccess) {
292 B2ERROR("OPALFitterGSL::fit: old parameters are garbage!");
293 return -1;
294 }
295
296
297 gsl_matrix_set_zero(Fetaxi);
298 for (int k = 0; k < ncon; ++k) {
299 constraints[k]->getDerivatives(Fetaxi->size2, Fetaxi->block->data + k * Fetaxi->tda);
300 }
301 if (debug > 1) debug_print(Fetaxi, "1: Fetaxi");
302 } else {
303 gsl_vector_memcpy(etasv, etaxi);
304 chik0 = chik;
305 chit0 = chit;
306 }
307
308 // Get covariance matrix
309 gsl_matrix_set_zero(V);
310
311 for (auto& fitobject : fitobjects) {
312 fitobject->addToGlobCov(V->block->data, V->tda);
313 }
314 if (debug > 1) debug_print(V, "V");
315
316 gsl_matrix_memcpy(VLU, &Vetaeta.matrix);
317
318 // invert covariance matrix (needed for chi2 calculation later)
319
320 int signum;
321 int result;
322 result = gsl_linalg_LU_decomp(VLU, permV, &signum);
323 B2DEBUG(11, "gsl_linalg_LU_decomp result=" << result);
324 if (debug > 3) debug_print(VLU, "VLU");
325
326 result = gsl_linalg_LU_invert(VLU, permV, Vinv);
327 B2DEBUG(11, "gsl_linalg_LU_invert result=" << result);
328
329 if (debug > 2) debug_print(Vinv, "Vinv");
330
331// *-- Evaluate f and S.
332 for (int k = 0; k < ncon; ++k) {
333 gsl_vector_set(f, k, constraints[k]->getValue());
334 }
335 if (debug > 1) debug_print(f, "f");
336
337 // y_eta = y - eta
338 gsl_vector_memcpy(y_eta, y);
339 gsl_vector_sub(y_eta, &eta.vector);
340 // r=f
341 gsl_vector_memcpy(r, f);
342 // r = 1*Feta*y_eta + 1*r
343 gsl_blas_dgemv(CblasNoTrans, 1, &Feta.matrix, y_eta, 1, r);
344
345 if (debug > 1) debug_print(r, "r");
346
347 // S = Feta * V * Feta^T
348
349 //FetaV = 1*Feta*V + 0*FetaV
350 B2DEBUG(12, "Creating FetaV");
351 gsl_blas_dsymm(CblasRight, CblasUpper, 1, &Vetaeta.matrix, &Feta.matrix, 0, FetaV);
352 // S = 1 * FetaV * Feta^T + 0*S
353 B2DEBUG(12, "Creating S");
354 gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1, FetaV, &Feta.matrix, 0, S);
355
356 if (nunm > 0) {
357 // New invention by B. List, 6.12.04:
358 // add F_xi * F_xi^T to S, to make method work when some
359 // constraints do not depend on any measured parameter
360
361 // Fxi is the part of Fetaxi containing the unmeasured quantities, if any
362 gsl_matrix_view Fxi = gsl_matrix_submatrix(Fetaxi, 0, nmea, ncon, nunm);
363
364 //S = 1*Fxi*Fxi^T + 1*S
365 gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1, &Fxi.matrix, &Fxi.matrix, 1, S);
366 }
367
368 if (debug > 1) debug_print(S, "S");
369
370// *-- Invert S to Sinv; S is destroyed here!
371// S is symmetric and positive definite
372
373 gsl_linalg_LU_decomp(S, permS, &signum);
374 int inverr = gsl_linalg_LU_invert(S, permS, Sinv);
375
376 if (inverr != 0) {
377 B2ERROR("S: gsl_linalg_LU_invert error " << inverr);
378 ierr = 7;
379 calcerr = false;
380 break;
381 }
382
383 // Calculate S^1*r here, we will need it
384 // Store it in lambda!
385 // lambda = 1*Sinv*r + 0*lambda; Sinv is symmetric
386 gsl_blas_dsymv(CblasUpper, 1, Sinv, r, 0, lambda);
387
388// *-- Calculate new unmeasured quantities, if any
389
390 if (nunm > 0) {
391 // Fxi is the part of Fetaxi containing the unmeasured quantities, if any
392 gsl_matrix_view Fxi = gsl_matrix_submatrix(Fetaxi, 0, nmea, ncon, nunm);
393 // W1 = Fxi^T * Sinv * Fxi
394 // SinvFxi = 1*Sinv*Fxi + 0*SinvFxi
395 gsl_blas_dsymm(CblasLeft, CblasUpper, 1, Sinv, &Fxi.matrix, 0, SinvFxi);
396 // W1 = 1*Fxi^T*SinvFxi + 0*W1
397 gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1, &Fxi.matrix, SinvFxi, 0, W1);
398
399 if (debug > 1) {
400 debug_print(W1, "W1");
401 // Check symmetry of W1
402 for (int i = 0; i < nunm; ++i) {
403 for (int j = 0; j < nunm; ++j) {
404 if (abs(gsl_matrix_get(W1, i, j) - gsl_matrix_get(W1, j, i)) > 1E-3 * abs(gsl_matrix_get(W1, i, j) + gsl_matrix_get(W1, j, i)))
405 B2INFO("W1[" << i << "][" << j << "] = " << gsl_matrix_get(W1, i, j)
406 << " W1[" << j << "][" << i << "] = " << gsl_matrix_get(W1, j, i)
407 << " => diff=" << abs(gsl_matrix_get(W1, i, j) - gsl_matrix_get(W1, j, i))
408 << " => tol=" << 1E-3 * abs(gsl_matrix_get(W1, i, j) + gsl_matrix_get(W1, j, i)));
409 }
410 }
411 }
412
413 // calculate shift of unmeasured parameters
414
415 // dxi = -alph*W1^-1 * Fxi^T*Sinv*r
416 // => dxi is solution of W1*dxi = -alph*Fxi^T*Sinv*r
417 // Compute rights hand side first
418 // Sinv*r was already calculated and is stored in lambda
419 // dxi = -alph*Fxi^T*lambda + 0*dxi
420
421 B2DEBUG(11, "alph = " << alph);
422 if (debug > 1) {
423 debug_print(lambda, "lambda");
424 debug_print(&(Fxi.matrix), "Fxi");
425 }
426
427 gsl_blas_dgemv(CblasTrans, -alph, &Fxi.matrix, lambda, 0, dxi);
428
429 if (debug > 1) {
430 debug_print(dxi, "dxi0");
431 debug_print(W1, "W1");
432 }
433
434 // now solve the system
435 // Note added 23.12.04: W1 is symmetric and positive definite,
436 // so we can use the Cholesky instead of LU decomposition
437 gsl_linalg_cholesky_decomp(W1);
438 inverr = gsl_linalg_cholesky_svx(W1, dxi);
439
440 if (debug > 1) debug_print(dxi, "dxi1");
441
442
443 if (inverr != 0) {
444 B2ERROR("W1: gsl_linalg_cholesky_svx error " << inverr);
445 ierr = 8;
446 calcerr = false;
447 break;
448 }
449
450 if (debug > 1) debug_print(dxi, "dxi");
451
452// *-- And now update unmeasured parameters; xi is a view of etaxi
453 // xi is the part of etaxi containing the unmeasured quantities, if any
454 gsl_vector_view xi = gsl_vector_subvector(etaxi, nmea, nunm);
455
456 // xi += dxi
457 gsl_vector_add(&xi.vector, dxi);
458
459 }
460
461// *-- Calculate new Lagrange multipliers.
462 // lambda = Sinv*r + Sinv*Fxi*dxi
463 // lambda is already set to Sinv*r, we just need to add Sinv*Fxi*dxi
464
465 // cppcheck-suppress duplicateCondition
466 if (nunm > 0) {
467 // Fxi is the part of Fetaxi containing the unmeasured quantities, if any
468 gsl_matrix_view Fxi = gsl_matrix_submatrix(Fetaxi, 0, nmea, ncon, nunm);
469 // calculate Fxidxi = 1*Fxi*dxi + 0*Fxidxi
470 gsl_blas_dgemv(CblasNoTrans, 1, &Fxi.matrix, dxi, 0, Fxidxi);
471 // add to existing lambda: lambda = 1*Sinv*Fxidxi + 1*lambda; Sinv is symmetric
472 gsl_blas_dsymv(CblasUpper, 1, Sinv, Fxidxi, 1, lambda);
473
474 }
475
476 if (debug > 1) debug_print(lambda, "lambda");
477
478// *-- Calculate new fitted parameters:
479 // eta = y - V*Feta^T*lambda
480 gsl_vector_memcpy(&eta.vector, y);
481 // FetaTlambda = 1*Feta^T*lambda + 0*FetaTlambda
482 gsl_blas_dgemv(CblasTrans, 1, &Feta.matrix, lambda, 0, FetaTlambda);
483 // eta = -1*V*FetaTlambda + 1*eta; V is symmetric
484 gsl_blas_dsymv(CblasUpper, -1, &Vetaeta.matrix, FetaTlambda, 1, &eta.vector);
485
486
487 if (debug > 1) debug_print(&eta.vector, "updated eta");
488
489// *-- Calculate constraints and their derivatives.
490 // since the constraints ask the fitobjects for their parameters,
491 // we need to update the fitobjects first!
492 // COULD BE DONE: update also ERRORS! (now only in the very end!)
493 updatesuccess = updateFitObjects(etaxi->block->data);
494
495 B2DEBUG(10, "After adjustment of all parameters:\n");
496 for (int k = 0; k < ncon; ++k) {
497 B2DEBUG(10, "Value of constraint " << k << " = " << constraints[k]->getValue());
498 }
499 gsl_matrix_set_zero(Fetaxi);
500 for (int k = 0; k < ncon; k++) {
501 constraints[k]->getDerivatives(Fetaxi->size2, Fetaxi->block->data + k * Fetaxi->tda);
502 }
503 if (debug > 1) debug_print(Fetaxi, "2: Fetaxi");
504
505
506// *-- Calculate new chisq.
507 // First, calculate new y-eta
508 // y_eta = y - eta
509 gsl_vector_memcpy(y_eta, y);
510 gsl_vector_sub(y_eta, &eta.vector);
511 // Now calculate Vinv*y_eta [ as solution to V* Vinvy_eta = y_eta]
512 // Vinvy_eta = 1*Vinv*y_eta + 0*Vinvy_eta; Vinv is symmetric
513 gsl_blas_dsymv(CblasUpper, 1, Vinv, y_eta, 0, Vinvy_eta);
514 // Now calculate y_eta *Vinvy_eta
515 gsl_blas_ddot(y_eta, Vinvy_eta, &chit);
516
517 for (int i = 0; i < nmea; ++i)
518 for (int j = 0; j < nmea; ++j) {
519 double dchit = (gsl_vector_get(y_eta, i)) *
520 gsl_matrix_get(Vinv, i, j) *
521 (gsl_vector_get(y_eta, j));
522 if (dchit != 0)
523 B2DEBUG(11, "chit for i,j = " << i << " , " << j << " = "
524 << dchit);
525 }
526
527 chik = 0;
528 for (int k = 0; k < ncon; ++k) chik += std::abs(2 * gsl_vector_get(lambda, k) * constraints[k]->getValue());
529 chinew = chit + chik;
530
531//*-- Calculate change in chisq, and check constraints are satisfied.
532 nit++;
533
534 bool sconv = (std::abs(chik - chik0) < dchikc)
535 && (std::abs(chit - chit0) < dchitc * chit)
536 && (chik < dchikt * chit);
537 // Second convergence criterium:
538 // If all constraints are fulfilled to better than 1E-8,
539 // and all parameters have changed by less than 1E-8,
540 // assume convergence
541 // This criterium assumes that all constraints and all parameters
542 // are "of order 1", i.e. their natural values are around 1 to 100,
543 // as for GeV or radians
544 double eps = 1E-6;
545 bool sconv2 = true;
546 for (int k = 0; sconv2 && (k < ncon); ++k)
547 sconv2 &= (std::abs(gsl_vector_get(f, k)) < eps);
548 if (sconv2)
549 B2DEBUG(10, "All constraints fulfilled to better than " << eps);
550
551 for (int j = 0; sconv2 && (j < npar); ++j)
552 sconv2 &= (std::abs(gsl_vector_get(etaxi, j) - gsl_vector_get(etasv, j)) < eps);
553 if (sconv2)
554 B2DEBUG(10, "All parameters stable to better than " << eps);
555 sconv |= sconv2;
556
557 bool sbad = (chik > dchik * chik0)
558 && (chik > dchikt * chit)
559 && (chik > chik0 + 1.E-10);
560
561 scut = false;
562
563 if (nit > nitmax) {
564// *-- Out of iterations
565 repeat = false;
566 calcerr = false;
567 ierr = 1;
568 } else if (sconv && updatesuccess) {
569// *-- Converged
570 repeat = false;
571 calcerr = true;
572 ierr = 0;
573 } else if (nit > 2 && chinew > chimxw && updatesuccess) {
574// *-- Chi2 crazy ?
575 repeat = false;
576 calcerr = false;
577 ierr = 2;
578 } else if ((sbad && nit > 1) || !updatesuccess) {
579// *-- ChiK increased, try smaller step
580 if (alph == almin) {
581 repeat = true; // false;
582 ierr = 3;
583 } else {
584 alph = std::max(almin, 0.5 * alph);
585 scut = true;
586 repeat = true;
587 ierr = 4;
588 }
589 } else {
590// *-- Keep going..
591 alph = std::min(alph + 0.1, 1.);
592 repeat = true;
593 ierr = 5;
594 }
595
596 B2DEBUG(10, "======== NIT = " << nit << ", CHI2 = " << chinew
597 << ", ierr = " << ierr << ", alph=" << alph);
598
599 for (unsigned int i = 0; i < fitobjects.size(); ++i)
600 B2DEBUG(10, "fitobject " << i << ": " << *fitobjects[i]);
601
602#ifndef FIT_TRACEOFF
603 if (tracer) tracer->step(*this);
604#endif
605
606 } // end of while (repeat)
607
608// *-- Turn chisq into probability.
609 fitprob = (ncon - nunm > 0) ? gsl_cdf_chisq_Q(chinew, ncon - nunm) : 0.5;
610 chi2 = chinew;
611
612// *-- End of iterations - calculate errors.
613
614// The result will (ultimately) be stored in Vnew
615
616 gsl_matrix_set_zero(Vnew);
617 gsl_matrix_set_zero(Minv);
618 if (debug > 2) {
619 for (int i = 0; i < npar; ++i) {
620 for (int j = 0; j < npar; ++j) {
621 B2INFO("Minv[" << i << "," << j << "]=" << gsl_matrix_get(Minv, i, j));
622 }
623 }
624 }
625
626 B2DEBUG(10, "OPALFitterGSL: calcerr = " << calcerr);
627
628 if (calcerr) {
629
630// *-- As a first step, calculate Minv as in 9.4.2 of Benno's book chapter
631// (in O.Behnke et al "Data Analysis in High Energy Physics")
632
633
634// *-- Evaluate S and invert.
635
636 if (debug > 2) {
637 debug_print(&Vetaeta.matrix, "V");
638 debug_print(&Feta.matrix, "Feta");
639 }
640
641 // CblasRight means C = alpha B A + beta C with symmetric matrix A
642 //FetaV[ncon][nmea] = 1*Feta[ncon][nmea]*V[nmea][nmea] + 0*FetaV
643 gsl_blas_dsymm(CblasRight, CblasUpper, 1, &Vetaeta.matrix, &Feta.matrix, 0, FetaV);
644 // S[ncon][ncon] = 1 * FetaV[ncon][nmea] * Feta^T[nmea][ncon] + 0*S
645 gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1, FetaV, &Feta.matrix, 0, S);
646
647
648 if (nunm > 0) {
649 // New invention by B. List, 6.12.04:
650 // add F_xi * F_xi^T to S, to make method work when some
651 // constraints do not depend on any measured parameter
652
653 // Fxi is the part of Fetaxi containing the unmeasured quantities, if any
654 gsl_matrix_view Fxi = gsl_matrix_submatrix(Fetaxi, 0, nmea, ncon, nunm);
655 //S[ncon][ncon] = 1*Fxi[ncon][nunm]*Fxi^T[nunm][ncon] + 1*S[ncon][ncon]
656 gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1, &Fxi.matrix, &Fxi.matrix, 1, S);
657 }
658
659 if (debug > 2) debug_print(S, "S");
660
661// *-- Invert S, testing for singularity first.
662// S is symmetric and positive definite
663
664 int signum;
665 gsl_linalg_LU_decomp(S, permS, &signum);
666 int inverr = gsl_linalg_LU_invert(S, permS, Sinv);
667
668 if (inverr != 0) {
669 B2ERROR("S: gsl_linalg_LU_invert error " << inverr << " in error calculation");
670 ierr = -1;
671 return -1;
672 }
673
674
675// *-- Calculate G.
676// (same as W1, but for measured parameters)
677// G = Feta^T * Sinv * Feta
678
679 // SinvFeta[ncon][nmea] = 1*Sinv[ncon][ncon]*Feta[ncon][nmea] + 0*SinvFeta
680 gsl_blas_dsymm(CblasLeft, CblasUpper, 1, Sinv, &Feta.matrix, 0, SinvFeta);
681 // G[nmea][nmea] = 1*Feta^T[nmea][ncon]*SinvFeta[ncon][nmea] + 0*G
682 gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1, &Feta.matrix, SinvFeta, 0, G);
683
684 if (debug > 2) debug_print(G, "G(1)");
685
686
687 if (nunm > 0) {
688
689// *-- Calculate H
690
691 // Fxi is the part of Fetaxi containing the unmeasured quantities, if any
692 gsl_matrix_view Fxi = gsl_matrix_submatrix(Fetaxi, 0, nmea, ncon, nunm);
693 // H = Feta^T * Sinv * Fxi
694 // SinvFxi[ncon][nunm] = 1*Sinv[ncon][ncon]*Fxi[ncon][nunm] + 0*SinvFxi
695 gsl_blas_dsymm(CblasLeft, CblasUpper, 1, Sinv, &Fxi.matrix, 0, SinvFxi);
696 // H[nmea][nunm] = 1*Feta^T[nmea][ncon]*SinvFxi[ncon][nunm] + 0*H
697 gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1, &Feta.matrix, SinvFxi, 0, H);
698
699 if (debug > 2) debug_print(H, "H");
700
701// *-- Calculate U**-1 and invert.
702// (same as W1)
703// U is a part of Minv
704
705 gsl_matrix* Uinv = W1;
706 gsl_matrix_view U = gsl_matrix_submatrix(Minv, nmea, nmea, nunm, nunm);
707 // Uinv = Fxi^T * Sinv * Fxi
708 // Uinv[nunm][nunm] = 1*Fxi^T[nunm][ncon]*SinvFxi[ncon][nunm] + 0*W1
709 gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1, &Fxi.matrix, SinvFxi, 0, Uinv);
710
711 gsl_linalg_LU_decomp(Uinv, permU, &signum);
712 inverr = gsl_linalg_LU_invert(Uinv, permU, &U.matrix);
713
714 if (debug > 2) debug_print(&U.matrix, "U");
715 for (int i = 0; i < npar; ++i) {
716 for (int j = 0; j < npar; ++j) {
717 B2DEBUG(12, "after U Minv[" << i << "," << j << "]=" << gsl_matrix_get(Minv, i, j));
718 }
719 }
720
721 if (inverr != 0) {
722 B2ERROR("U: gsl_linalg_LU_invert error " << inverr << " in error calculation ");
723 ierr = -1;
724 return -1;
725 }
726
727// *-- Covariance matrix between measured and unmeasured parameters.
728
729// HU[nmea][nunm] = 1*H[nmea][nunm]*U[nunm][nunm] + 0*HU
730 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, H, &U.matrix, 0, HU);
731// Vnewetaxi is a view of Vnew
732 gsl_matrix_view Minvetaxi = gsl_matrix_submatrix(Minv, 0, nmea, nmea, nunm);
733 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, -1, &Vetaeta.matrix, HU, 0, &Minvetaxi.matrix);
734 for (int i = 0; i < npar; ++i) {
735 for (int j = 0; j < npar; ++j) {
736 B2DEBUG(12, "after etaxi Minv[" << i << "," << j << "]=" << gsl_matrix_get(Minv, i, j));
737 }
738 }
739
740
741// *-- Fill in symmetric part:
742// Vnewxieta is a view of Vnew
743 gsl_matrix_view Minvxieta = gsl_matrix_submatrix(Minv, nmea, 0, nunm, nmea);
744 gsl_matrix_transpose_memcpy(&Minvxieta.matrix, &Minvetaxi.matrix);
745 for (int i = 0; i < npar; ++i) {
746 for (int j = 0; j < npar; ++j) {
747 B2DEBUG(12, "after symmetric: Minv[" << i << "," << j << "]=" << gsl_matrix_get(Minv, i, j));
748 }
749 }
750
751// *-- Calculate G-HUH^T:
752// G = -1*HU*H^T +1*G
753 gsl_blas_dgemm(CblasNoTrans, CblasTrans, -1, HU, H, +1, G);
754
755 } // endif nunm > 0
756
757// *-- Calculate I-GV.
758 // IGV = 1
759 gsl_matrix_set_identity(IGV);
760 // IGV = -1*G*Vetaeta + 1*IGV
761 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, -1, G, &Vetaeta.matrix, 1, IGV);
762
763// *-- And finally error matrix on fitted parameters.
764 gsl_matrix_view Minvetaeta = gsl_matrix_submatrix(Minv, 0, 0, nmea, nmea);
765
766 // Vnewetaeta = 1*Vetaeta*IGV + 0*Vnewetaeta
767 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, &Vetaeta.matrix, IGV, 0, &Minvetaeta.matrix);
768
769 for (int i = 0; i < npar; ++i) {
770 for (int j = 0; j < npar; ++j) {
771 B2DEBUG(12, "complete Minv[" << i << "," << j << "]=" << gsl_matrix_get(Minv, i, j));
772 }
773 }
774
775// *-- now Minv should be complete, can calculate new covariance matrix Vnew
776// Vnew = (dx / dt)^T * V_t * dx / dt
777// x are the fitted parameters, t are the measured parameters,
778// V_t is the covariance matrix of the measured parameters
779// thus in OPALFitter variables: V_t = Vetaeta, x = (eta, xi) = etaxi
780// d eta / dt and d xi / dt are given by eqns 9.54 and 9.55, respectively,
781// as deta/dt = - Minvetaeta * Fetat and dxi/dt = - Minvxieta * Fetat
782// Fetat is d^2 chi^2 / deta dt = - V^-1, even in presence of soft constraints, since
783// the soft constraints don't depend on the _measured_ parameters t (but on the fitted ones)
784// HERE, V (and Vinv) do not include the second derivatives of the soft constraints
785// so we can use them right away
786// Note that "F" is used for completely different things:
787// in OPALFitter it is the derivatives of the constraints wrt to the parameters (A in book)
788// in the book, F are the second derivatives of the objective function wrt the parameters
789
790
791 gsl_matrix_view detadt = gsl_matrix_submatrix(dxdt, 0, 0, nmea, nmea);
792 B2DEBUG(13, "after detadt");
793 // Vdxdt is Vetaeta * dxdt^T, thus Vdxdt[nmea][npar]
794 gsl_matrix_view Vdetadt = gsl_matrix_submatrix(Vdxdt, 0, 0, nmea, nmea);
795 B2DEBUG(13, "after Vdetadt");
796
797 // detadt = - Minvetaeta * Fetat = -1 * Minvetaeta * (-1) * Vinv + 0 * detadt // replace by symm?
798 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, &Minvetaeta.matrix, Vinv, 0, &detadt.matrix);
799 if (debug > 2) debug_print(&detadt.matrix, "deta/dt");
800
801 // Vdetadt = 1 * Vetaeta * detadt^T + 0* Vdetadt
802 gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1, &Vetaeta.matrix, &detadt.matrix, 0, &Vdetadt.matrix); // ok
803 if (debug > 2) debug_print(&Vdetadt.matrix, "Vetata * deta/dt");
804
805 gsl_matrix_view Vnewetaeta = gsl_matrix_submatrix(Vnew, 0, 0, nmea, nmea); //[nmea],[nmea]
806 // Vnewetaeta = 1 * detadt * Vdetadt + 0* Vnewetaeta
807 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, &detadt.matrix, &Vdetadt.matrix, 0, &Vnewetaeta.matrix);
808
809 if (debug > 2) debug_print(Vnew, "Vnew after part for measured parameters");
810
811
812 if (nunm > 0) {
813
814 gsl_matrix_view Minvxieta = gsl_matrix_submatrix(Minv, nmea, 0, nunm, nmea); //[nunm][nmea]
815 B2DEBUG(13, "after Minvxieta");
816 if (debug > 2) debug_print(&Minvxieta.matrix, "Minvxieta");
817
818 gsl_matrix_view dxidt = gsl_matrix_submatrix(dxdt, nmea, 0, nunm, nmea); //[nunm][nmea]
819 B2DEBUG(13, "after dxidt");
820 // dxidt[nunm][nmea] = - Minvxieta * Fetat = -1 * Minvxieta[nunm][nmea] * Vinv[nmea][nmea] + 0 * dxidt
821 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, &Minvxieta.matrix, Vinv, 0, &dxidt.matrix); //ok
822 if (debug > 2) debug_print(&dxidt.matrix, "dxi/dt");
823
824 // Vdxdt = V * dxdt^T => Vdxdt[nmea][npar]
825 gsl_matrix_view Vdxidt = gsl_matrix_submatrix(Vdxdt, 0, nmea, nmea, nunm); //[nmea][nunm]
826 B2DEBUG(13, "after Vdxidt");
827 // Vdxidt = 1 * Vetaeta[nmea][nmea] * dxidt^T[nmea][nunm] + 0* Vdxidt => Vdxidt[nmea][nunm]
828 gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1, &Vetaeta.matrix, &dxidt.matrix, 0, &Vdxidt.matrix); // ok
829 if (debug > 2) debug_print(&Vdxidt.matrix, "Vetaeta * dxi/dt^T");
830
831 gsl_matrix_view Vnewetaxi = gsl_matrix_submatrix(Vnew, 0, nmea, nmea, nunm); //[nmea][nunm]
832 gsl_matrix_view Vnewxieta = gsl_matrix_submatrix(Vnew, nmea, 0, nunm, nmea); //[nunm][nmea]
833 gsl_matrix_view Vnewxixi = gsl_matrix_submatrix(Vnew, nmea, nmea, nunm, nunm); //[nunm][nunm]
834
835 // Vnewxieta[nunm][nmea] = 1 * dxidt[nunm][nmea] * Vdetadt[nmea][nmea] + 0* Vnewxieta
836 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, &dxidt.matrix, &Vdetadt.matrix, 0, &Vnewxieta.matrix); // ok
837 if (debug > 2) debug_print(Vnew, "Vnew after xieta part");
838 // Vnewetaxi[nmea][nunm] = 1 * detadt[nmea][nmea] * Vdxidt[nmea][nunm] + 0* Vnewetaxi
839 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, &detadt.matrix, &Vdxidt.matrix, 0, &Vnewetaxi.matrix); // ok
840 if (debug > 2) debug_print(Vnew, "Vnew after etaxi part");
841 // Vnewxixi[nunm][nunm] = 1 * dxidt[nunm][nmea] * Vdxidt[nmea][nunm] + 0* Vnewxixi
842 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, &dxidt.matrix, &Vdxidt.matrix, 0, &Vnewxixi.matrix);
843 if (debug > 2) debug_print(Vnew, "Vnew after xixi part");
844 }
845
846// *-- now we finally have Vnew, fill into fitobjects
847
848 // update errors in fitobjects
849 for (auto& fitobject : fitobjects) {
850 for (int ilocal = 0; ilocal < fitobject->getNPar(); ++ilocal) {
851 int iglobal = fitobject->getGlobalParNum(ilocal);
852 for (int jlocal = ilocal; jlocal < fitobject->getNPar(); ++jlocal) {
853 int jglobal = fitobject->getGlobalParNum(jlocal);
854 if (iglobal >= 0 && jglobal >= 0)
855 fitobject->setCov(ilocal, jlocal, gsl_matrix_get(Vnew, iglobal, jglobal));
856 }
857 }
858 }
859
860 // Finally, copy covariance matrix
861 if (cov && covDim != nmea + nunm) {
862 delete[] cov;
863 cov = nullptr;
864 }
865 covDim = nmea + nunm;
866 if (!cov) cov = new double[covDim * covDim];
867 for (int i = 0; i < covDim; ++i) {
868 for (int j = 0; j < covDim; ++j) {
869 cov[i * covDim + j] = gsl_matrix_get(Vnew, i, j);
870 }
871 }
872 covValid = true;
873
874 } // endif calcerr == true
875
876#ifndef FIT_TRACEOFF
877 if (tracer) tracer->finish(*this);
878#endif
879
880 return fitprob;
881
882 }
double * cov
global covariance matrix of last fit problem
Definition: BaseFitter.h:104
int covDim
dimension of global covariance matrix
Definition: BaseFitter.h:103
virtual void step(BaseFitter &fitter)
Called at the end of each step.
Definition: BaseTracer.cc:35
virtual void initialize(BaseFitter &fitter)
Called at the start of a new fit (during initialization)
Definition: BaseTracer.cc:30
virtual void finish(BaseFitter &fitter)
Called at the end of a fit.
Definition: BaseTracer.cc:45

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

113 {
114 return &constraints;
115 }

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

108 {
109 return &fitobjects;
110 }

◆ getGlobalCovarianceMatrix() [1/2]

double * getGlobalCovarianceMatrix ( int &  idim)
virtualinherited
Parameters
idim1st dimension of global covariance matrix

Definition at line 158 of file BaseFitter.cc.

159 {
160 if (covValid && cov) {
161 idim = covDim;
162 return cov;
163 }
164 idim = 0;
165 return nullptr;
166 }

◆ getGlobalCovarianceMatrix() [2/2]

const double * getGlobalCovarianceMatrix ( int &  idim) const
virtualinherited
Parameters
idim1st dimension of global covariance matrix

Definition at line 148 of file BaseFitter.cc.

149 {
150 if (covValid && cov) {
151 idim = covDim;
152 return cov;
153 }
154 idim = 0;
155 return nullptr;
156 }

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

118 {
119 return &softconstraints;
120 }

◆ getTracer() [1/2]

BaseTracer * getTracer ( )
virtualinherited

Definition at line 130 of file BaseFitter.cc.

131 {
132 return tracer;
133 }

◆ getTracer() [2/2]

const BaseTracer * getTracer ( ) const
virtualinherited

Definition at line 134 of file BaseFitter.cc.

135 {
136 return tracer;
137 }

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

1042 {
1043 if (m) {
1044 if (m->size1 != size1 || m->size2 != size2) {
1045 gsl_matrix_free(m);
1046 if (size1 * size2 > 0) m = gsl_matrix_alloc(size1, size2);
1047 else m = nullptr;
1048 }
1049 } else if (size1 * size2 > 0) m = gsl_matrix_alloc(size1, size2);
1050 }

◆ ini_gsl_permutation()

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

Definition at line 1018 of file OPALFitterGSL.cc.

1019 {
1020 if (p) {
1021 if (p->size != size) {
1022 gsl_permutation_free(p);
1023 if (size > 0) p = gsl_permutation_alloc(size);
1024 else p = nullptr;
1025 }
1026 } else if (size > 0) p = gsl_permutation_alloc(size);
1027 }

◆ ini_gsl_vector()

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

Definition at line 1029 of file OPALFitterGSL.cc.

1030 {
1031
1032 if (v) {
1033 if (v->size != size) {
1034 gsl_vector_free(v);
1035 if (size > 0) v = gsl_vector_alloc(size);
1036 else v = nullptr;
1037 }
1038 } else if (size > 0) v = gsl_vector_alloc(size);
1039 }

◆ initialize()

bool initialize ( )
overridevirtual

Implements BaseFitter.

Definition at line 884 of file OPALFitterGSL.cc.

885 {
886 covValid = false;
887 // tell fitobjects the global ordering of their parameters:
888 int iglobal = 0;
889 // measured parameters first!
890 for (auto& fitobject : fitobjects) {
891 for (int ilocal = 0; ilocal < fitobject->getNPar(); ++ilocal) {
892 if (fitobject->isParamMeasured(ilocal) &&
893 !fitobject->isParamFixed(ilocal)) {
894 fitobject->setGlobalParNum(ilocal, iglobal);
895 B2DEBUG(10, "Object " << fitobject->getName()
896 << " Parameter " << fitobject->getParamName(ilocal)
897 << " is measured, global number " << iglobal);
898 ++iglobal;
899 }
900 }
901 }
902 nmea = iglobal;
903 // now unmeasured parameters!
904 for (auto& fitobject : fitobjects) {
905 for (int ilocal = 0; ilocal < fitobject->getNPar(); ++ilocal) {
906 if (!fitobject->isParamMeasured(ilocal) &&
907 !fitobject->isParamFixed(ilocal)) {
908 fitobject->setGlobalParNum(ilocal, iglobal);
909 B2DEBUG(10, "Object " << fitobject->getName()
910 << " Parameter " << fitobject->getParamName(ilocal)
911 << " is unmeasured, global number " << iglobal);
912 ++iglobal;
913 }
914 }
915 }
916 npar = iglobal;
917 assert(npar <= NPARMAX);
918 nunm = npar - nmea;
919 assert(nunm <= NUNMMAX);
920
921 // set number of constraints
922 ncon = constraints.size();
923 assert(ncon <= NCONMAX);
924
925 ini_gsl_vector(f, ncon);
926 ini_gsl_vector(r, ncon);
927
928 ini_gsl_matrix(Fetaxi, ncon, npar);
929 ini_gsl_matrix(S, ncon, ncon);
930 ini_gsl_matrix(Sinv, ncon, ncon);
931 ini_gsl_matrix(SinvFxi, ncon, nunm);
932 ini_gsl_matrix(SinvFeta, ncon, nmea);
933 ini_gsl_matrix(W1, nunm, nunm);
934 ini_gsl_matrix(G, nmea, nmea);
935 ini_gsl_matrix(H, nmea, nunm);
936 ini_gsl_matrix(HU, nmea, nunm);
937 ini_gsl_matrix(IGV, nmea, nmea);
938 ini_gsl_matrix(V, npar, npar);
939 ini_gsl_matrix(VLU, nmea, nmea);
940 ini_gsl_matrix(Vinv, nmea, nmea);
941 ini_gsl_matrix(Vnew, npar, npar);
942 ini_gsl_matrix(Minv, npar, npar);
943 ini_gsl_matrix(dxdt, npar, nmea);
944 ini_gsl_matrix(Vdxdt, nmea, npar);
945
946 ini_gsl_vector(dxi, nunm);
947 ini_gsl_vector(Fxidxi, ncon);
948 ini_gsl_vector(lambda, ncon);
949 ini_gsl_vector(FetaTlambda, nmea);
950
951 ini_gsl_vector(etaxi, npar);
952 ini_gsl_vector(etasv, npar);
953 ini_gsl_vector(y, nmea);
954 ini_gsl_vector(y_eta, nmea);
955 ini_gsl_vector(Vinvy_eta, nmea);
956
957 ini_gsl_matrix(FetaV, ncon, nmea);
958
959 ini_gsl_permutation(permS, ncon);
960 ini_gsl_permutation(permU, nunm);
961 ini_gsl_permutation(permV, nmea);
962
963 assert(f && (int)f->size == ncon);
964 assert(r && (int)r->size == ncon);
965 assert(Fetaxi && (int)Fetaxi->size1 == ncon && (int)Fetaxi->size2 == npar);
966 assert(S && (int)S->size1 == ncon && (int)S->size2 == ncon);
967 assert(Sinv && (int)Sinv->size1 == ncon && (int)Sinv->size2 == ncon);
968 assert(nunm == 0 || (SinvFxi && (int)SinvFxi->size1 == ncon && (int)SinvFxi->size2 == nunm));
969 assert(SinvFeta && (int)SinvFeta->size1 == ncon && (int)SinvFeta->size2 == nmea);
970 assert(nunm == 0 || (W1 && (int)W1->size1 == nunm && (int)W1->size2 == nunm));
971 assert(G && (int)G->size1 == nmea && (int)G->size2 == nmea);
972 assert(nunm == 0 || (H && (int)H->size1 == nmea && (int)H->size2 == nunm));
973 assert(nunm == 0 || (HU && (int)HU->size1 == nmea && (int)HU->size2 == nunm));
974 assert(IGV && (int)IGV->size1 == nmea && (int)IGV->size2 == nmea);
975 assert(V && (int)V->size1 == npar && (int)V->size2 == npar);
976 assert(VLU && (int)VLU->size1 == nmea && (int)VLU->size2 == nmea);
977 assert(Vinv && (int)Vinv->size1 == nmea && (int)Vinv->size2 == nmea);
978 assert(Vnew && (int)Vnew->size1 == npar && (int)Vnew->size2 == npar);
979 assert(nunm == 0 || (dxi && (int)dxi->size == nunm));
980 assert(nunm == 0 || (Fxidxi && (int)Fxidxi->size == ncon));
981 assert(lambda && (int)lambda->size == ncon);
982 assert(FetaTlambda && (int)FetaTlambda->size == nmea);
983 assert(etaxi && (int)etaxi->size == npar);
984 assert(etasv && (int)etasv->size == npar);
985 assert(y && (int)y->size == nmea);
986 assert(y_eta && (int)y_eta->size == nmea);
987 assert(Vinvy_eta && (int)Vinvy_eta->size == nmea);
988 assert(FetaV && (int)FetaV->size1 == ncon && (int)FetaV->size2 == nmea);
989 assert(permS && (int)permS->size == ncon);
990 assert(permS && (int)permS->size == ncon);
991 assert(nunm == 0 || (permU && (int)permU->size == nunm));
992 assert(permV && (int)permV->size == nmea);
993
994 return true;
995
996 }

◆ reset()

void reset ( )
virtualinherited

Definition at line 122 of file BaseFitter.cc.

123 {
124 fitobjects.resize(0);
125 constraints.resize(0);
126 softconstraints.resize(0);
127 covValid = false;
128 }

◆ setDebug()

void setDebug ( int  debuglevel)
virtual

Set the Debug Level.

Definition at line 1052 of file OPALFitterGSL.cc.

1053 {
1054 debug = debuglevel;
1055 }

◆ setTracer() [1/2]

void setTracer ( BaseTracer newTracer)
virtualinherited

Definition at line 143 of file BaseFitter.cc.

144 {
145 tracer = &newTracer;
146 }

◆ setTracer() [2/2]

void setTracer ( BaseTracer newTracer)
virtualinherited

Definition at line 138 of file BaseFitter.cc.

139 {
140 tracer = newTracer;
141 }

◆ updateFitObjects()

bool updateFitObjects ( double  etaxi[])
protectedvirtual

Definition at line 998 of file OPALFitterGSL.cc.

999 {
1000 // changed etaxi -> eetaxi to avoid clash with class member DJeans
1001 //bool debug = false;
1002 bool result = true;
1003 for (auto& fitobject : fitobjects) {
1004 for (int ilocal = 0; ilocal < fitobject->getNPar(); ++ilocal) {
1005 fitobject->updateParams(eetaxi, npar);
1006// }
1007 }
1008 }
1009 return result;
1010 }

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: