diff options
Diffstat (limited to 'src/fit.c')
| -rw-r--r-- | src/fit.c | 205 |
1 files changed, 205 insertions, 0 deletions
diff --git a/src/fit.c b/src/fit.c new file mode 100644 index 0000000..b6e441a --- /dev/null +++ b/src/fit.c @@ -0,0 +1,205 @@ +#include "likelihood.h" +#include <stdio.h> +#include "zebra.h" +#include "Record_Info.h" +#include "event.h" +#include "zdab_utils.h" +#include <gsl/gsl_vector.h> +#include <gsl/gsl_multimin.h> +#include "scattering.h" +#include "pmt.h" +#include "sno_charge.h" +#include "db.h" +#include "dqxx.h" + +#define EV_RECORD 0x45562020 // 'EV ' (as written to ZDAB file) + +double nll(const gsl_vector *v, void *params) +{ + double T, theta, phi, t0; + double pos[3], dir[3]; + + T = gsl_vector_get(v,0); + + /* FIXME: should use a minimization algorithm with bounds. */ + if (T < 0) return 1e9; + + pos[0] = gsl_vector_get(v,1); + pos[1] = gsl_vector_get(v,2); + pos[2] = gsl_vector_get(v,3); + + theta = gsl_vector_get(v,4); + phi = gsl_vector_get(v,5); + dir[0] = sin(theta)*cos(phi); + dir[1] = sin(theta)*sin(phi); + dir[2] = cos(theta); + + t0 = gsl_vector_get(v,6); + + return nll_muon((event *) params, T, pos, dir, t0); +} + +int main(int argc, char **argv) +{ + int i; + zebraFile *f; + bank b; + int rv; + PMTBank bpmt; + EVBank bev; + event ev = {0}; + int crate, card, channel; + int id; + + if (argc < 2) { + fprintf(stderr, "usage: fit [filename]\n"); + return 1; + } + + f = zebra_open(argv[1]); + + if (!f) { + fprintf(stderr, "%s\n", zebra_err); + return 1; + } + + load_pmt_info(); + + for (i = 0; i < 10000; i++) { + ev.pmt_hits[i].hit = 0; + ev.pmt_hits[i].flags = 0; + } + + ev.run = -1; + + const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex2; + gsl_multimin_fminimizer *s = NULL; + gsl_vector *ss, *x; + gsl_multimin_function minex_func; + + size_t iter = 0; + int status; + double size; + + /* Starting point */ + x = gsl_vector_alloc(7); + + ss = gsl_vector_alloc(7); + + minex_func.n = 7; + minex_func.f = nll; + minex_func.params = &ev; + + init_interpolation(); + init_charge(); + dict *db = db_init(); + + if (load_file(db, "DQXX_0000010000.dat")) { + fprintf(stderr, "failed to load DQXX_0000010000.dat: %s\n", db_err); + exit(1); + } + + if (dqxx_init(db, &ev)) { + fprintf(stderr, "failed to initialize DQXX bank: %s\n", dqxx_err); + exit(1); + } + + while (1) { + rv = next_bank(f, &b); + + if (rv == -1) { + fprintf(stderr, "error getting bank: %s\n", zebra_err); + goto err; + } else if (rv == 1) { + /* EOF */ + break; + } + + if (b.name == PMT_RECORD) { + unpack_pmt(b.data, &bpmt); + card = bpmt.pin/1024; + crate = (bpmt.pin % 1024)/32; + channel = bpmt.pin % 32; + id = crate*512 + card*32 + channel; + ev.pmt_hits[id].hit = 1; + ev.pmt_hits[id].t = bpmt.pt; + ev.pmt_hits[id].qhl = bpmt.phl; + ev.pmt_hits[id].qhs = bpmt.phs; + ev.pmt_hits[id].qlx = bpmt.plx; + //ev.pmt_hits[id].flags |= bpmt.pn & KPF_DIS; + } else if (b.name == EV_RECORD) { + unpack_ev(b.data, &bev); + if (ev.run != -1) { + if (ev.run != bev.run || ev.gtid != bev.gtr_id) { + /* New event, so we need to fit the old event. */ + iter = 0; + + gsl_vector_set(x,0,1000.0); + gsl_vector_set(x,1,100.0); + gsl_vector_set(x,2,100.0); + gsl_vector_set(x,3,100.0); + gsl_vector_set(x,4,1.57); + gsl_vector_set(x,5,0.0); + gsl_vector_set(x,6,120.0); + + gsl_vector_set(ss,0,10.0); + gsl_vector_set(ss,1,10.0); + gsl_vector_set(ss,2,10.0); + gsl_vector_set(ss,3,10.0); + gsl_vector_set(ss,4,0.1); + gsl_vector_set(ss,5,0.1); + gsl_vector_set(ss,6,1.0); + + s = gsl_multimin_fminimizer_alloc(T, 7); + gsl_multimin_fminimizer_set(s, &minex_func, x, ss); + + do { + iter++; + status = gsl_multimin_fminimizer_iterate(s); + + if (status) break; + + size = gsl_multimin_fminimizer_size(s); + status = gsl_multimin_test_size(size,1e-2); + + if (status == GSL_SUCCESS) { + printf("converged to minimum at\n"); + } + + printf("%5zu %.2g %.2g %.2g %.2g %.2g %.2g %.2g f() = %7.3e size = %.3f\n", + iter, + gsl_vector_get(s->x,0), + gsl_vector_get(s->x,1), + gsl_vector_get(s->x,2), + gsl_vector_get(s->x,3), + gsl_vector_get(s->x,4), + gsl_vector_get(s->x,5), + gsl_vector_get(s->x,6), + s->fval, + size); + } while (status == GSL_CONTINUE && iter < 1000); + + gsl_multimin_fminimizer_free(s); + } + } + ev.run = bev.run; + ev.gtid = bev.gtr_id; + } + } + + free_interpolation(); + + db_free(db); + + gsl_vector_free(x); + gsl_vector_free(ss); + + zebra_close(f); + + return 0; + +err: + + zebra_close(f); + return 1; +} |