diff options
Diffstat (limited to 'src/fit.c')
-rw-r--r-- | src/fit.c | 37 |
1 files changed, 26 insertions, 11 deletions
@@ -5243,7 +5243,13 @@ int fit_event2(event *ev, double *xopt, double *fmin, int *id, size_t n, double double x[100], ss[100], lb[100], ub[100], fval, n_d2o, x0[100], T0, Tmin, mass, pos[3], t0, dir[3]; struct timeval tv_start, tv_stop; double time_elapsed; + double peak_theta[MAX_NPEAKS]; + double peak_phi[MAX_NPEAKS]; + size_t npeaks; + size_t *result; + size_t nvertices; + /* Create the minimizer object. */ opt = nlopt_create(NLOPT_LN_BOBYQA, 4+3*n); nlopt_set_min_objective(opt,nlopt_nll2,&fpars); @@ -5261,49 +5267,54 @@ int fit_event2(event *ev, double *xopt, double *fmin, int *id, size_t n, double t0 = guess_t0(ev,pos); } + /* Set the initial parameters for the parameters which don't change during + * the fast minimization. */ x0[0] = pos[0]; x0[1] = pos[1]; x0[2] = pos[2]; x0[3] = t0; + /* Set the initial step sizes for the parameters which don't change during + * the fast minimization. */ ss[0] = 10.0; ss[1] = 10.0; ss[2] = 10.0; ss[3] = 1.0; + /* Set the lower bound for the parameters which don't change during the + * fast minimization. */ lb[0] = -1000.0; lb[1] = -1000.0; lb[2] = -1000.0; lb[3] = 0.0; + /* Set the upper bound for the parameters which don't change during the + * fast minimization. */ ub[0] = +1000.0; ub[1] = +1000.0; ub[2] = +1000.0; ub[3] = GTVALID; - double peak_theta[MAX_NPEAKS]; - double peak_phi[MAX_NPEAKS]; - size_t npeaks; - + /* Find the peaks in the Hough transform of the event. */ find_peaks(ev, pos, 100, 100, peak_theta, peak_phi, &npeaks, LEN(peak_theta), 0.5); - size_t *result = malloc(sizeof(size_t)*n*ipow(npeaks,n)); - - size_t nvertices; + result = malloc(sizeof(size_t)*n*ipow(npeaks,n)); + /* Compute the set of all unique combinations of particles and directions. */ unique_vertices(id,n,npeaks,result,&nvertices); + /* Compute the index of refraction in D2O. */ n_d2o = get_index_snoman_d2o(400.0); fpars.ev = ev; + /* Set the global iteration count back to zero. + * + * FIXME: Some way to do this which is thread safe? */ iter = 0; /* First we do a set of "quick" minimizations to try and start the - * minimizer close to the minimum. To make these function evaluations - * faster, we set the absolute tolerance on the likelihood to 1.0, the - * maximum number of function evaluations to 100, and the relative - * tolerance on the numerical integration to 10%. */ + * minimizer close to the minimum. */ fpars.dx = 1.0; fpars.dx_shower = 10.0; fpars.fast = 1; @@ -5318,6 +5329,7 @@ int fit_event2(event *ev, double *xopt, double *fmin, int *id, size_t n, double fpars.id[i] = id[i]; for (i = 0; i < nvertices; i++) { + /* Copy the starting parameters to the `x` array. */ memcpy(x,x0,sizeof(x)); for (j = 0; j < n; j++) { x[5+3*j] = peak_theta[result[i*n+j]]; @@ -5419,6 +5431,7 @@ int fit_event2(event *ev, double *xopt, double *fmin, int *id, size_t n, double ub[2] = +1000.0; ub[3] = GTVALID; + /* Set the energy step size to 2% of the energy. */ for (i = 0; i < nvertices; i++) { for (j = 0; j < n; j++) { ss[4+3*j] = xopt[4+3*j]*0.02; @@ -5476,11 +5489,13 @@ int fit_event2(event *ev, double *xopt, double *fmin, int *id, size_t n, double end: + free(result); nlopt_destroy(opt); return 0; stop: + free(result); nlopt_destroy(opt); return NLOPT_FORCED_STOP; |