Belle II Software development
Public Member Functions | Public Attributes | Protected Types | Protected Member Functions | Static Protected Member Functions | Protected Attributes | Private Attributes | List of all members
OPALFitterGSL Class Reference

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

#include <OPALFitterGSL.h>

Inheritance diagram for OPALFitterGSL:
BaseFitter

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.

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 {}

◆ ~OPALFitterGSL()

~OPALFitterGSL ( )
overridevirtual

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 }

◆ 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 1056 of file OPALFitterGSL.cc.

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

◆ debug_print() [2/2]

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

Definition at line 1064 of file OPALFitterGSL.cc.

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

◆ 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 if (nunm > 0) {
466 // Fxi is the part of Fetaxi containing the unmeasured quantities, if any
467 gsl_matrix_view Fxi = gsl_matrix_submatrix(Fetaxi, 0, nmea, ncon, nunm);
468 // calculate Fxidxi = 1*Fxi*dxi + 0*Fxidxi
469 gsl_blas_dgemv(CblasNoTrans, 1, &Fxi.matrix, dxi, 0, Fxidxi);
470 // add to existing lambda: lambda = 1*Sinv*Fxidxi + 1*lambda; Sinv is symmetric
471 gsl_blas_dsymv(CblasUpper, 1, Sinv, Fxidxi, 1, lambda);
472
473 }
474
475 if (debug > 1) debug_print(lambda, "lambda");
476
477// *-- Calculate new fitted parameters:
478 // eta = y - V*Feta^T*lambda
479 gsl_vector_memcpy(&eta.vector, y);
480 // FetaTlambda = 1*Feta^T*lambda + 0*FetaTlambda
481 gsl_blas_dgemv(CblasTrans, 1, &Feta.matrix, lambda, 0, FetaTlambda);
482 // eta = -1*V*FetaTlambda + 1*eta; V is symmetric
483 gsl_blas_dsymv(CblasUpper, -1, &Vetaeta.matrix, FetaTlambda, 1, &eta.vector);
484
485
486 if (debug > 1) debug_print(&eta.vector, "updated eta");
487
488// *-- Calculate constraints and their derivatives.
489 // since the constraints ask the fitobjects for their parameters,
490 // we need to update the fitobjects first!
491 // COULD BE DONE: update also ERRORS! (now only in the very end!)
492 updatesuccess = updateFitObjects(etaxi->block->data);
493
494 B2DEBUG(10, "After adjustment of all parameters:\n");
495 for (int k = 0; k < ncon; ++k) {
496 B2DEBUG(10, "Value of constraint " << k << " = " << constraints[k]->getValue());
497 }
498 gsl_matrix_set_zero(Fetaxi);
499 for (int k = 0; k < ncon; k++) {
500 constraints[k]->getDerivatives(Fetaxi->size2, Fetaxi->block->data + k * Fetaxi->tda);
501 }
502 if (debug > 1) debug_print(Fetaxi, "2: Fetaxi");
503
504
505// *-- Calculate new chisq.
506 // First, calculate new y-eta
507 // y_eta = y - eta
508 gsl_vector_memcpy(y_eta, y);
509 gsl_vector_sub(y_eta, &eta.vector);
510 // Now calculate Vinv*y_eta [ as solution to V* Vinvy_eta = y_eta]
511 // Vinvy_eta = 1*Vinv*y_eta + 0*Vinvy_eta; Vinv is symmetric
512 gsl_blas_dsymv(CblasUpper, 1, Vinv, y_eta, 0, Vinvy_eta);
513 // Now calculate y_eta *Vinvy_eta
514 gsl_blas_ddot(y_eta, Vinvy_eta, &chit);
515
516 for (int i = 0; i < nmea; ++i)
517 for (int j = 0; j < nmea; ++j) {
518 double dchit = (gsl_vector_get(y_eta, i)) *
519 gsl_matrix_get(Vinv, i, j) *
520 (gsl_vector_get(y_eta, j));
521 if (dchit != 0)
522 B2DEBUG(11, "chit for i,j = " << i << " , " << j << " = "
523 << dchit);
524 }
525
526 chik = 0;
527 for (int k = 0; k < ncon; ++k) chik += std::abs(2 * gsl_vector_get(lambda, k) * constraints[k]->getValue());
528 chinew = chit + chik;
529
530//*-- Calculate change in chisq, and check constraints are satisfied.
531 nit++;
532
533 bool sconv = (std::abs(chik - chik0) < dchikc)
534 && (std::abs(chit - chit0) < dchitc * chit)
535 && (chik < dchikt * chit);
536 // Second convergence criterium:
537 // If all constraints are fulfilled to better than 1E-8,
538 // and all parameters have changed by less than 1E-8,
539 // assume convergence
540 // This criterium assumes that all constraints and all parameters
541 // are "of order 1", i.e. their natural values are around 1 to 100,
542 // as for GeV or radians
543 double eps = 1E-6;
544 bool sconv2 = true;
545 for (int k = 0; sconv2 && (k < ncon); ++k)
546 sconv2 &= (std::abs(gsl_vector_get(f, k)) < eps);
547 if (sconv2)
548 B2DEBUG(10, "All constraints fulfilled to better than " << eps);
549
550 for (int j = 0; sconv2 && (j < npar); ++j)
551 sconv2 &= (std::abs(gsl_vector_get(etaxi, j) - gsl_vector_get(etasv, j)) < eps);
552 if (sconv2)
553 B2DEBUG(10, "All parameters stable to better than " << eps);
554 sconv |= sconv2;
555
556 bool sbad = (chik > dchik * chik0)
557 && (chik > dchikt * chit)
558 && (chik > chik0 + 1.E-10);
559
560 scut = false;
561
562 if (nit > nitmax) {
563// *-- Out of iterations
564 repeat = false;
565 calcerr = false;
566 ierr = 1;
567 } else if (sconv && updatesuccess) {
568// *-- Converged
569 repeat = false;
570 calcerr = true;
571 ierr = 0;
572 } else if (nit > 2 && chinew > chimxw && updatesuccess) {
573// *-- Chi2 crazy ?
574 repeat = false;
575 calcerr = false;
576 ierr = 2;
577 } else if ((sbad && nit > 1) || !updatesuccess) {
578// *-- ChiK increased, try smaller step
579 if (alph == almin) {
580 repeat = true; // false;
581 ierr = 3;
582 } else {
583 alph = std::max(almin, 0.5 * alph);
584 scut = true;
585 repeat = true;
586 ierr = 4;
587 }
588 } else {
589// *-- Keep going..
590 alph = std::min(alph + 0.1, 1.);
591 repeat = true;
592 ierr = 5;
593 }
594
595 B2DEBUG(10, "======== NIT = " << nit << ", CHI2 = " << chinew
596 << ", ierr = " << ierr << ", alph=" << alph);
597
598 for (unsigned int i = 0; i < fitobjects.size(); ++i)
599 B2DEBUG(10, "fitobject " << i << ": " << *fitobjects[i]);
600
601#ifndef FIT_TRACEOFF
602 if (tracer) tracer->step(*this);
603#endif
604
605 } // end of while (repeat)
606
607// *-- Turn chisq into probability.
608 fitprob = (ncon - nunm > 0) ? gsl_cdf_chisq_Q(chinew, ncon - nunm) : 0.5;
609 chi2 = chinew;
610
611// *-- End of iterations - calculate errors.
612
613// The result will (ultimately) be stored in Vnew
614
615 gsl_matrix_set_zero(Vnew);
616 gsl_matrix_set_zero(Minv);
617 if (debug > 2) {
618 for (int i = 0; i < npar; ++i) {
619 for (int j = 0; j < npar; ++j) {
620 B2INFO("Minv[" << i << "," << j << "]=" << gsl_matrix_get(Minv, i, j));
621 }
622 }
623 }
624
625 B2DEBUG(10, "OPALFitterGSL: calcerr = " << calcerr);
626
627 if (calcerr) {
628
629// *-- As a first step, calculate Minv as in 9.4.2 of Benno's book chapter
630// (in O.Behnke et al "Data Analysis in High Energy Physics")
631
632
633// *-- Evaluate S and invert.
634
635 if (debug > 2) {
636 debug_print(&Vetaeta.matrix, "V");
637 debug_print(&Feta.matrix, "Feta");
638 }
639
640 // CblasRight means C = alpha B A + beta C with symmetric matrix A
641 //FetaV[ncon][nmea] = 1*Feta[ncon][nmea]*V[nmea][nmea] + 0*FetaV
642 gsl_blas_dsymm(CblasRight, CblasUpper, 1, &Vetaeta.matrix, &Feta.matrix, 0, FetaV);
643 // S[ncon][ncon] = 1 * FetaV[ncon][nmea] * Feta^T[nmea][ncon] + 0*S
644 gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1, FetaV, &Feta.matrix, 0, S);
645
646
647 if (nunm > 0) {
648 // New invention by B. List, 6.12.04:
649 // add F_xi * F_xi^T to S, to make method work when some
650 // constraints do not depend on any measured parameter
651
652 // Fxi is the part of Fetaxi containing the unmeasured quantities, if any
653 gsl_matrix_view Fxi = gsl_matrix_submatrix(Fetaxi, 0, nmea, ncon, nunm);
654 //S[ncon][ncon] = 1*Fxi[ncon][nunm]*Fxi^T[nunm][ncon] + 1*S[ncon][ncon]
655 gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1, &Fxi.matrix, &Fxi.matrix, 1, S);
656 }
657
658 if (debug > 2) debug_print(S, "S");
659
660// *-- Invert S, testing for singularity first.
661// S is symmetric and positive definite
662
663 int signum;
664 gsl_linalg_LU_decomp(S, permS, &signum);
665 int inverr = gsl_linalg_LU_invert(S, permS, Sinv);
666
667 if (inverr != 0) {
668 B2ERROR("S: gsl_linalg_LU_invert error " << inverr << " in error calculation");
669 ierr = -1;
670 return -1;
671 }
672
673
674// *-- Calculate G.
675// (same as W1, but for measured parameters)
676// G = Feta^T * Sinv * Feta
677
678 // SinvFeta[ncon][nmea] = 1*Sinv[ncon][ncon]*Feta[ncon][nmea] + 0*SinvFeta
679 gsl_blas_dsymm(CblasLeft, CblasUpper, 1, Sinv, &Feta.matrix, 0, SinvFeta);
680 // G[nmea][nmea] = 1*Feta^T[nmea][ncon]*SinvFeta[ncon][nmea] + 0*G
681 gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1, &Feta.matrix, SinvFeta, 0, G);
682
683 if (debug > 2) debug_print(G, "G(1)");
684
685
686 if (nunm > 0) {
687
688// *-- Calculate H
689
690 // Fxi is the part of Fetaxi containing the unmeasured quantities, if any
691 gsl_matrix_view Fxi = gsl_matrix_submatrix(Fetaxi, 0, nmea, ncon, nunm);
692 // H = Feta^T * Sinv * Fxi
693 // SinvFxi[ncon][nunm] = 1*Sinv[ncon][ncon]*Fxi[ncon][nunm] + 0*SinvFxi
694 gsl_blas_dsymm(CblasLeft, CblasUpper, 1, Sinv, &Fxi.matrix, 0, SinvFxi);
695 // H[nmea][nunm] = 1*Feta^T[nmea][ncon]*SinvFxi[ncon][nunm] + 0*H
696 gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1, &Feta.matrix, SinvFxi, 0, H);
697
698 if (debug > 2) debug_print(H, "H");
699
700// *-- Calculate U**-1 and invert.
701// (same as W1)
702// U is a part of Minv
703
704 gsl_matrix* Uinv = W1;
705 gsl_matrix_view U = gsl_matrix_submatrix(Minv, nmea, nmea, nunm, nunm);
706 // Uinv = Fxi^T * Sinv * Fxi
707 // Uinv[nunm][nunm] = 1*Fxi^T[nunm][ncon]*SinvFxi[ncon][nunm] + 0*W1
708 gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1, &Fxi.matrix, SinvFxi, 0, Uinv);
709
710 gsl_linalg_LU_decomp(Uinv, permU, &signum);
711 inverr = gsl_linalg_LU_invert(Uinv, permU, &U.matrix);
712
713 if (debug > 2) debug_print(&U.matrix, "U");
714 for (int i = 0; i < npar; ++i) {
715 for (int j = 0; j < npar; ++j) {
716 B2DEBUG(12, "after U Minv[" << i << "," << j << "]=" << gsl_matrix_get(Minv, i, j));
717 }
718 }
719
720 if (inverr != 0) {
721 B2ERROR("U: gsl_linalg_LU_invert error " << inverr << " in error calculation ");
722 ierr = -1;
723 return -1;
724 }
725
726// *-- Covariance matrix between measured and unmeasured parameters.
727
728// HU[nmea][nunm] = 1*H[nmea][nunm]*U[nunm][nunm] + 0*HU
729 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, H, &U.matrix, 0, HU);
730// Vnewetaxi is a view of Vnew
731 gsl_matrix_view Minvetaxi = gsl_matrix_submatrix(Minv, 0, nmea, nmea, nunm);
732 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, -1, &Vetaeta.matrix, HU, 0, &Minvetaxi.matrix);
733 for (int i = 0; i < npar; ++i) {
734 for (int j = 0; j < npar; ++j) {
735 B2DEBUG(12, "after etaxi Minv[" << i << "," << j << "]=" << gsl_matrix_get(Minv, i, j));
736 }
737 }
738
739
740// *-- Fill in symmetric part:
741// Vnewxieta is a view of Vnew
742 gsl_matrix_view Minvxieta = gsl_matrix_submatrix(Minv, nmea, 0, nunm, nmea);
743 gsl_matrix_transpose_memcpy(&Minvxieta.matrix, &Minvetaxi.matrix);
744 for (int i = 0; i < npar; ++i) {
745 for (int j = 0; j < npar; ++j) {
746 B2DEBUG(12, "after symmetric: Minv[" << i << "," << j << "]=" << gsl_matrix_get(Minv, i, j));
747 }
748 }
749
750// *-- Calculate G-HUH^T:
751// G = -1*HU*H^T +1*G
752 gsl_blas_dgemm(CblasNoTrans, CblasTrans, -1, HU, H, +1, G);
753
754 } // endif nunm > 0
755
756// *-- Calculate I-GV.
757 // IGV = 1
758 gsl_matrix_set_identity(IGV);
759 // IGV = -1*G*Vetaeta + 1*IGV
760 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, -1, G, &Vetaeta.matrix, 1, IGV);
761
762// *-- And finally error matrix on fitted parameters.
763 gsl_matrix_view Minvetaeta = gsl_matrix_submatrix(Minv, 0, 0, nmea, nmea);
764
765 // Vnewetaeta = 1*Vetaeta*IGV + 0*Vnewetaeta
766 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, &Vetaeta.matrix, IGV, 0, &Minvetaeta.matrix);
767
768 for (int i = 0; i < npar; ++i) {
769 for (int j = 0; j < npar; ++j) {
770 B2DEBUG(12, "complete Minv[" << i << "," << j << "]=" << gsl_matrix_get(Minv, i, j));
771 }
772 }
773
774// *-- now Minv should be complete, can calculate new covariance matrix Vnew
775// Vnew = (dx / dt)^T * V_t * dx / dt
776// x are the fitted parameters, t are the measured parameters,
777// V_t is the covariance matrix of the measured parameters
778// thus in OPALFitter variables: V_t = Vetaeta, x = (eta, xi) = etaxi
779// d eta / dt and d xi / dt are given by eqns 9.54 and 9.55, respectively,
780// as deta/dt = - Minvetaeta * Fetat and dxi/dt = - Minvxieta * Fetat
781// Fetat is d^2 chi^2 / deta dt = - V^-1, even in presence of soft constraints, since
782// the soft constraints don't depend on the _measured_ parameters t (but on the fitted ones)
783// HERE, V (and Vinv) do not include the second derivatives of the soft constraints
784// so we can use them right away
785// Note that "F" is used for completely different things:
786// in OPALFitter it is the derivatives of the constraints wrt to the parameters (A in book)
787// in the book, F are the second derivatives of the objective function wrt the parameters
788
789
790 gsl_matrix_view detadt = gsl_matrix_submatrix(dxdt, 0, 0, nmea, nmea);
791 B2DEBUG(13, "after detadt");
792 // Vdxdt is Vetaeta * dxdt^T, thus Vdxdt[nmea][npar]
793 gsl_matrix_view Vdetadt = gsl_matrix_submatrix(Vdxdt, 0, 0, nmea, nmea);
794 B2DEBUG(13, "after Vdetadt");
795
796 // detadt = - Minvetaeta * Fetat = -1 * Minvetaeta * (-1) * Vinv + 0 * detadt // replace by symm?
797 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, &Minvetaeta.matrix, Vinv, 0, &detadt.matrix);
798 if (debug > 2) debug_print(&detadt.matrix, "deta/dt");
799
800 // Vdetadt = 1 * Vetaeta * detadt^T + 0* Vdetadt
801 gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1, &Vetaeta.matrix, &detadt.matrix, 0, &Vdetadt.matrix); // ok
802 if (debug > 2) debug_print(&Vdetadt.matrix, "Vetata * deta/dt");
803
804 gsl_matrix_view Vnewetaeta = gsl_matrix_submatrix(Vnew, 0, 0, nmea, nmea); //[nmea],[nmea]
805 // Vnewetaeta = 1 * detadt * Vdetadt + 0* Vnewetaeta
806 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, &detadt.matrix, &Vdetadt.matrix, 0, &Vnewetaeta.matrix);
807
808 if (debug > 2) debug_print(Vnew, "Vnew after part for measured parameters");
809
810
811 if (nunm > 0) {
812
813 gsl_matrix_view Minvxieta = gsl_matrix_submatrix(Minv, nmea, 0, nunm, nmea); //[nunm][nmea]
814 B2DEBUG(13, "after Minvxieta");
815 if (debug > 2) debug_print(&Minvxieta.matrix, "Minvxieta");
816
817 gsl_matrix_view dxidt = gsl_matrix_submatrix(dxdt, nmea, 0, nunm, nmea); //[nunm][nmea]
818 B2DEBUG(13, "after dxidt");
819 // dxidt[nunm][nmea] = - Minvxieta * Fetat = -1 * Minvxieta[nunm][nmea] * Vinv[nmea][nmea] + 0 * dxidt
820 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, &Minvxieta.matrix, Vinv, 0, &dxidt.matrix); //ok
821 if (debug > 2) debug_print(&dxidt.matrix, "dxi/dt");
822
823 // Vdxdt = V * dxdt^T => Vdxdt[nmea][npar]
824 gsl_matrix_view Vdxidt = gsl_matrix_submatrix(Vdxdt, 0, nmea, nmea, nunm); //[nmea][nunm]
825 B2DEBUG(13, "after Vdxidt");
826 // Vdxidt = 1 * Vetaeta[nmea][nmea] * dxidt^T[nmea][nunm] + 0* Vdxidt => Vdxidt[nmea][nunm]
827 gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1, &Vetaeta.matrix, &dxidt.matrix, 0, &Vdxidt.matrix); // ok
828 if (debug > 2) debug_print(&Vdxidt.matrix, "Vetaeta * dxi/dt^T");
829
830 gsl_matrix_view Vnewetaxi = gsl_matrix_submatrix(Vnew, 0, nmea, nmea, nunm); //[nmea][nunm]
831 gsl_matrix_view Vnewxieta = gsl_matrix_submatrix(Vnew, nmea, 0, nunm, nmea); //[nunm][nmea]
832 gsl_matrix_view Vnewxixi = gsl_matrix_submatrix(Vnew, nmea, nmea, nunm, nunm); //[nunm][nunm]
833
834 // Vnewxieta[nunm][nmea] = 1 * dxidt[nunm][nmea] * Vdetadt[nmea][nmea] + 0* Vnewxieta
835 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, &dxidt.matrix, &Vdetadt.matrix, 0, &Vnewxieta.matrix); // ok
836 if (debug > 2) debug_print(Vnew, "Vnew after xieta part");
837 // Vnewetaxi[nmea][nunm] = 1 * detadt[nmea][nmea] * Vdxidt[nmea][nunm] + 0* Vnewetaxi
838 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, &detadt.matrix, &Vdxidt.matrix, 0, &Vnewetaxi.matrix); // ok
839 if (debug > 2) debug_print(Vnew, "Vnew after etaxi part");
840 // Vnewxixi[nunm][nunm] = 1 * dxidt[nunm][nmea] * Vdxidt[nmea][nunm] + 0* Vnewxixi
841 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, &dxidt.matrix, &Vdxidt.matrix, 0, &Vnewxixi.matrix);
842 if (debug > 2) debug_print(Vnew, "Vnew after xixi part");
843 }
844
845// *-- now we finally have Vnew, fill into fitobjects
846
847 // update errors in fitobjects
848 for (auto& fitobject : fitobjects) {
849 for (int ilocal = 0; ilocal < fitobject->getNPar(); ++ilocal) {
850 int iglobal = fitobject->getGlobalParNum(ilocal);
851 for (int jlocal = ilocal; jlocal < fitobject->getNPar(); ++jlocal) {
852 int jglobal = fitobject->getGlobalParNum(jlocal);
853 if (iglobal >= 0 && jglobal >= 0)
854 fitobject->setCov(ilocal, jlocal, gsl_matrix_get(Vnew, iglobal, jglobal));
855 }
856 }
857 }
858
859 // Finally, copy covariance matrix
860 if (cov && covDim != nmea + nunm) {
861 delete[] cov;
862 cov = nullptr;
863 }
864 covDim = nmea + nunm;
865 if (!cov) cov = new double[covDim * covDim];
866 for (int i = 0; i < covDim; ++i) {
867 for (int j = 0; j < covDim; ++j) {
868 cov[i * covDim + j] = gsl_matrix_get(Vnew, i, j);
869 }
870 }
871 covValid = true;
872
873 } // endif calcerr == true
874
875#ifndef FIT_TRACEOFF
876 if (tracer) tracer->finish(*this);
877#endif
878
879 return fitprob;
880
881 }
E
R E
internal precision of FFTW codelets
Definition DiscreteCosineTransform_31points.cc:36

◆ getChi2()

double getChi2 ( ) const
overridevirtual

Get the chi**2 of the last fit.

Implements BaseFitter.

Definition at line 1013 of file OPALFitterGSL.cc.

1013{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 1014 of file OPALFitterGSL.cc.

1014{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 1011 of file OPALFitterGSL.cc.

1011{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 1015 of file OPALFitterGSL.cc.

1015{return nit;}

◆ getNcon()

int getNcon ( ) const
virtual

Get the number of hard constraints of the last fit.

Definition at line 1071 of file OPALFitterGSL.cc.

1071{return ncon;}

◆ getNpar()

int getNpar ( ) const
virtual

Get the number of all parameters of the last fit.

Definition at line 1074 of file OPALFitterGSL.cc.

1074{return npar;}

◆ getNsoft()

int getNsoft ( ) const
virtual

Get the number of soft constraints of the last fit.

Definition at line 1072 of file OPALFitterGSL.cc.

1072{return 0;}

◆ getNunm()

int getNunm ( ) const
virtual

Get the number of unmeasured parameters of the last fit.

Definition at line 1073 of file OPALFitterGSL.cc.

1073{return nunm;}

◆ getProbability()

double getProbability ( ) const
overridevirtual

Get the fit probability of the last fit.

Implements BaseFitter.

Definition at line 1012 of file OPALFitterGSL.cc.

1012{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 1040 of file OPALFitterGSL.cc.

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

◆ ini_gsl_permutation()

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

Definition at line 1017 of file OPALFitterGSL.cc.

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

◆ ini_gsl_vector()

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

Definition at line 1028 of file OPALFitterGSL.cc.

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

◆ initialize()

bool initialize ( )
overridevirtual

Implements BaseFitter.

Definition at line 883 of file OPALFitterGSL.cc.

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

◆ 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 1051 of file OPALFitterGSL.cc.

1052 {
1053 debug = debuglevel;
1054 }

◆ 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 997 of file OPALFitterGSL.cc.

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

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: