diff options
Diffstat (limited to 'src/fit.c')
| -rw-r--r-- | src/fit.c | 85 |
1 files changed, 58 insertions, 27 deletions
@@ -22,6 +22,7 @@ #include <signal.h> /* for signal() */ #include "release.h" #include "id_particles.h" +#include "misc.h" char *GitSHA1(void); char *GitDirty(void); @@ -40,8 +41,14 @@ typedef struct fitParams { int n; int fast; int print; + int id; } fitParams; +int particles[] = { + IDP_E_MINUS, + IDP_MU_MINUS, +}; + /* In order to start the fitter close to the minimum, we first do a series of * "quick" minimizations starting at the following points. We keep track of the * parameters with the best likelihood value and then start the "real" @@ -4996,7 +5003,7 @@ double nll(unsigned int n, const double *x, double *grad, void *params) z2[0] = x[8]; gettimeofday(&tv_start, NULL); - fval = nll_muon(fpars->ev, IDP_MU_MINUS, T, pos, dir, t0, z1, z2, 1, fpars->n, fpars->fast); + fval = nll_muon(fpars->ev, fpars->id, T, pos, dir, t0, z1, z2, 1, fpars->n, fpars->fast); gettimeofday(&tv_stop, NULL); long long elapsed = (tv_stop.tv_sec - tv_start.tv_sec)*1000 + (tv_stop.tv_usec - tv_start.tv_usec)/1000; @@ -5099,11 +5106,11 @@ void guess_direction(event *ev, double *pos, double *theta, double *phi) *phi = atan2(dir[1],dir[0]); } -int fit_event(event *ev, double *xopt, double *fmin) +int fit_event(event *ev, double *xopt, double *fmin, int id) { size_t i; fitParams fpars; - double x[9], ss[9], lb[9], ub[9], fval, n, qhs_sum, x0[9], T0, Tmin; + double x[9], ss[9], lb[9], ub[9], fval, n, qhs_sum, x0[9], T0, Tmin, mass; struct timeval tv_start, tv_stop; opt = nlopt_create(NLOPT_LN_BOBYQA, 9); @@ -5124,12 +5131,29 @@ int fit_event(event *ev, double *xopt, double *fmin) n = get_index_snoman_d2o(400.0); + switch (id) { + case IDP_E_MINUS: + case IDP_E_PLUS: + mass = ELECTRON_MASS; + break; + case IDP_MU_MINUS: + case IDP_MU_PLUS: + mass = MUON_MASS; + break; + case IDP_PROTON: + mass = PROTON_MASS; + break; + default: + fprintf(stderr, "unknown particle id %i\n", id); + exit(1); + } + /* Calculate the Cerenkov threshold for a muon. */ - Tmin = MUON_MASS/sqrt(1.0-1.0/(n*n)) - MUON_MASS; + Tmin = mass/sqrt(1.0-1.0/(n*n)) - mass; /* If our guess is below the Cerenkov threshold, start at the Cerenkov * threshold. */ - double Tmin2 = sqrt(MUON_MASS*MUON_MASS/(1-pow(0.9,2)))-MUON_MASS; + double Tmin2 = sqrt(mass*mass/(1-pow(0.9,2)))-mass; if (T0 < Tmin2) T0 = Tmin2; x0[0] = T0; @@ -5178,6 +5202,7 @@ int fit_event(event *ev, double *xopt, double *fmin) nlopt_set_initial_step(opt, ss); fpars.ev = ev; + fpars.id = id; iter = 0; @@ -5538,32 +5563,38 @@ int main(int argc, char **argv) rv = get_event(f,&ev,&b); - gettimeofday(&tv_start, NULL); - if (fit_event(&ev,xopt,&fmin) == NLOPT_FORCED_STOP) { - printf("ctrl-c caught. quitting...\n"); - goto end; - } - gettimeofday(&tv_stop, NULL); - - long long elapsed = (tv_stop.tv_sec - tv_start.tv_sec)*1000 + (tv_stop.tv_usec - tv_start.tv_usec)/1000; - if (fout) { if (first_ev) fprintf(fout, " ev:\n"); fprintf(fout, " - gtid: %i\n", ev.gtid); fprintf(fout, " fit:\n"); - fprintf(fout, " -\n"); - fprintf(fout, " energy: %.2f\n", xopt[0]); - fprintf(fout, " posx: %.2f\n", xopt[1]); - fprintf(fout, " posy: %.2f\n", xopt[2]); - fprintf(fout, " posz: %.2f\n", xopt[3]); - fprintf(fout, " theta: %.4f\n", xopt[4]); - fprintf(fout, " phi: %.4f\n", xopt[5]); - fprintf(fout, " t0: %.2f\n", xopt[6]); - fprintf(fout, " z1: %.2f\n", xopt[7]); - fprintf(fout, " z2: %.2f\n", xopt[8]); - fprintf(fout, " fmin: %.2f\n", fmin); - fprintf(fout, " time: %lld\n", elapsed); - fflush(fout); + } + + for (i = 0; i < LEN(particles); i++) { + gettimeofday(&tv_start, NULL); + if (fit_event(&ev,xopt,&fmin,particles[i]) == NLOPT_FORCED_STOP) { + printf("ctrl-c caught. quitting...\n"); + goto end; + } + gettimeofday(&tv_stop, NULL); + + long long elapsed = (tv_stop.tv_sec - tv_start.tv_sec)*1000 + (tv_stop.tv_usec - tv_start.tv_usec)/1000; + + if (fout) { + fprintf(fout, " -\n"); + fprintf(fout, " id: %i\n", particles[i]); + fprintf(fout, " energy: %.2f\n", xopt[0]); + fprintf(fout, " posx: %.2f\n", xopt[1]); + fprintf(fout, " posy: %.2f\n", xopt[2]); + fprintf(fout, " posz: %.2f\n", xopt[3]); + fprintf(fout, " theta: %.4f\n", xopt[4]); + fprintf(fout, " phi: %.4f\n", xopt[5]); + fprintf(fout, " t0: %.2f\n", xopt[6]); + fprintf(fout, " z1: %.2f\n", xopt[7]); + fprintf(fout, " z2: %.2f\n", xopt[8]); + fprintf(fout, " fmin: %.2f\n", fmin); + fprintf(fout, " time: %lld\n", elapsed); + fflush(fout); + } } first_ev = 0; |