1 #include <CLHEP/Matrix/Vector.h>
2 #include <CLHEP/Matrix/Matrix.h>
8 #define _DBG__ cerr<<__FILE__<<":"<<__LINE__<<endl;
12 printf (
"print_state: iter: %3u "
15 gsl_blas_dnrm2 (s->f));
24 cout <<
"number of points less than or equal to number of parameters\n";
28 covar = gsl_matrix_alloc(p, p);
29 T = gsl_multifit_fdfsolver_lmsder;
30 s = gsl_multifit_fdfsolver_alloc (
T, (
size_t)n, (
size_t)p);
35 gsl_multifit_fdfsolver_free(
s);
36 gsl_matrix_free(
covar);
41 cout <<
"setStartParams called\n";
42 cout <<
"start params: " << xIn << endl;
46 if (
debug_level >=2) cout <<
"number of data points = " << n <<
" number of parameters = " << p << endl;
55 xGslPtr = gsl_vector_alloc((
size_t)p);
56 for (
unsigned int i = 0; i < p; i++) {
57 gsl_vector_set(xGslPtr, (
size_t)i, xIn(i + 1));
58 if (
debug_level >=2) cout << i <<
" " << gsl_vector_get(xGslPtr, (
size_t)i) << endl;
60 if (
debug_level >=2) cout <<
"initialize the solver" << endl;
61 gsl_multifit_fdfsolver_set(
s, &
f, xGslPtr);
62 gsl_vector_free(xGslPtr);
67 for (
unsigned int i = 0; i < p; i++) {
68 xOut(i + 1) = gsl_vector_get (
s->x, (
size_t)i);
72 double chi = gsl_blas_dnrm2(
s->f);
83 status = gsl_multifit_fdfsolver_iterate(
s);
85 printf(
"chisqMin::fit: iter = %d, status = %d, text = %s\n", iter,
status, gsl_strerror (
status));
89 if (
debug_level >= 2) cout <<
"chisqMin::fit: continue not taken, do delta test" << endl;
90 status = gsl_multifit_test_delta(
s->dx,
s->x, 1
e-4, 1
e-4);
92 printf(
"chisqMin::fit: from delta test, status = %d, text = %s\n",
status, gsl_strerror (
status));
109 gsl_matrix* gslCovarPtr = gsl_matrix_alloc(nparams, nparams);
111 gsl_multifit_covar(
s->J, 0.0, gslCovarPtr);
112 HepMatrix
covar(nparams, nparams);
114 for (
int i = 0; i < nparams; i++) {
115 for (
int j = 0; j < nparams; j++) {
116 covar(i + 1, j + 1) = gsl_matrix_get(gslCovarPtr, i, j);
119 gsl_matrix_free(gslCovarPtr);
virtual unsigned int getP()=0
virtual unsigned int getN()=0
chisqMin(residFunc *f_in, int level=1)
gsl_multifit_function_fdf f
int fGsl(const gsl_vector *x, void *data, gsl_vector *f)
const gsl_multifit_fdfsolver_type * T
gsl_multifit_fdfsolver * s
void setStartParams(HepVector &x)
int fdfGsl(const gsl_vector *x, void *data, gsl_vector *f, gsl_matrix *J)
printf("string=%s", string)
int dfGsl(const gsl_vector *x, void *data, gsl_matrix *J)
void print_state(size_t iter, gsl_multifit_fdfsolver *s)