From 92929a7efc25939a34473830bb69cecbc5371427 Mon Sep 17 00:00:00 2001 From: cfuselli Date: Sun, 26 Mar 2023 20:23:31 +0000 Subject: [PATCH 1/4] linposfit works --- straxen/plugins/peaks/__init__.py | 3 + straxen/plugins/peaks/peak_positions_lpf.py | 814 ++++++++++++++++++++ 2 files changed, 817 insertions(+) create mode 100644 straxen/plugins/peaks/peak_positions_lpf.py diff --git a/straxen/plugins/peaks/__init__.py b/straxen/plugins/peaks/__init__.py index 38de82403..94e197ed9 100644 --- a/straxen/plugins/peaks/__init__.py +++ b/straxen/plugins/peaks/__init__.py @@ -22,6 +22,9 @@ from . import peak_positions_mlp from .peak_positions_mlp import * +from . import peak_positions_lpf +from .peak_positions_lpf import * + from . import peak_proximity from .peak_proximity import * diff --git a/straxen/plugins/peaks/peak_positions_lpf.py b/straxen/plugins/peaks/peak_positions_lpf.py new file mode 100644 index 000000000..437e3f08d --- /dev/null +++ b/straxen/plugins/peaks/peak_positions_lpf.py @@ -0,0 +1,814 @@ +import strax +import straxen +from straxen.plugins.peaks._peak_positions_base import PeakPositionsBaseNT + + +export, __all__ = strax.exporter() + +import strax +import numpy as np +import numba +import straxen + + +@export +class PeakPositionsLPF(strax.Plugin): + """ + Computes S2 position with linearized log-likelihood fit. + Returns xy position and the parameters of the fit. + See LinPosFit https://github.com/XAMS-nikhef/LinPosFit/ from A.P. Colijn + """ + depends_on=('peaks','peak_basics') + provides = 'peak_positions_lpf' + rechunk_on_save=False + dtype= [('xml',np.float), + ('yml',np.float), + ('r0',np.float), + ('gamma', np.float), + ('logl',np.float), + ('n',np.int)] + dtype += strax.time_fields + + n_top_pmts = straxen.URLConfig(default=straxen.n_top_pmts, infer_type=False, + help="Number of top PMTs") + + pmt_to_lxe = straxen.URLConfig(default=7, infer_type=False, + help="Distance between the PMTs and the liquid interface") + + def setup(self,): + inch = 2.54 # cm + pmt_pos = straxen.pmt_positions() + self.pmt_pos = list(zip(pmt_pos['x'].values,pmt_pos['y'].values,np.repeat(self.pmt_to_lxe, self.n_top_pmts))) + self.pmt_surface=(3*inch)**2*np.pi/4. + + def compute(self,peaks): + res = np.zeros(len(peaks),self.dtype) + + for ix, p in enumerate(peaks): + if p['type']!=2: + #Only reconstruct s2 peaks. We do need to set the time of the peaks + res[ix]['time'] = p['time'] + res[ix]['endtime'] = p['endtime'] + continue + try: + #Some really small single electron s2s fail the minimization + fit_result,_,_ = lpf_execute(self.pmt_pos[:self.n_top_pmts],p['area_per_channel'][:self.n_top_pmts],self.pmt_surface) + + res[ix]['time'] = p['time'] + res[ix]['endtime'] = p['endtime'] + res[ix]['xml'] = fit_result[0] + res[ix]['yml'] = fit_result[1] + res[ix]['r0'] = fit_result[2] + res[ix]['logl'] = fit_result[3] + res[ix]['n'] = fit_result[4] + res[ix]['gamma'] = fit_result[5] + + except Exception as e: + print(e) + res[ix]['time'] = p['time'] + res[ix]['endtime'] = p['endtime'] + continue + + return res + + + + +import numpy as np +from numba import njit + +# from matplotlib.colors import BoundaryNorm +import matplotlib.pyplot as plt +from matplotlib.gridspec import GridSpec +from matplotlib.ticker import MaxNLocator +from IPython.display import clear_output + +# lpf = linearized position fit +# +# A.P. Colijn - Feb. 2021 +lpf_iter_max = 100 +lpf_min_position_cor = 0.01 # convergence criterium for position minimizer +lpf_npar = 4 # number of fit parameters +# x y r0 gamma +lpf_fix_par = np.array([False, False, False, False]) # which parameters to fix +# +debug = False + + +def lpf_execute(xhit, nhit, area_sensor): + #%snakeviz + """ + Execute the fitter + :param xhit: + :param nhit: + :param area_sensor: + :param kwargs: + :return: + """ + xhit = np.array(xhit) + nhit = np.array(nhit) + # minimize the likelihood + fit_result, xiter, hit_is_used = lpf_minimize(xhit, nhit, area_sensor) + return fit_result, xiter, hit_is_used + +@njit +def lpf_factorial(x): + n = 1 + for i in range(2, x+1): + n *= i + return n + +@njit +def lpf_lnlike(xhit, nhit, xf, nuv, gamma, area_sensor): + r0 = nuv * area_sensor * xhit[0][2] / 4 / np.pi + logl = 0 + + # print('next') + for ih in range(len(nhit)): + nexp = lpf_nexp(xhit[ih], xf, r0, gamma) + # print(xhit[ih], nexp, xf, r0) + + logl = logl + nexp - nhit[ih] * np.log(nexp) + + return logl + +@njit +def lpf_lnlike_plot(xhit, nhit, hit_is_used, xf, r0, gamma): + logl = 0 + + for ih in range(len(nhit)): + nexp = lpf_nexp(xhit[ih], xf, r0, gamma) + d = lpf_dist(xhit[ih],xf) + if hit_is_used[ih] == 1: + logl = logl + nexp - nhit[ih] * np.log(nexp) + + return logl + +@njit +def lpf_lnlike_r(xhit, nhit, xf, r0, gamma): + logl = 0 + + for ih in range(len(nhit)): + nexp = lpf_nexp(xhit[ih], xf, r0, gamma) + logl = logl + nexp - nhit[ih] * np.log(nexp) + + return logl + +@njit +def lpf_initialize(xhit, nhit): + """ + Initialize the lpf_fitter. Estimate the initial position and uv photons + + :param xhit: + :param nhit: + :return: + """ + nmax = -1 + + xhit_max = np.zeros(3) + # sort hits in descending order. helps when looping over the hits + # + indices = np.argsort(-nhit) # descending order + + # find a cluster of hits + # + # 1. we start with the highest hit in the event + # 2. if it is surrrounded by other hits then we found a cluster and we are done + # 3. if not: continue at 1. with the second highest hit in the event + # 4. repeat until we checked all possibilities + # 5. if no cluster was found we use the highest hit as the initial fit position + # + + # only consider hits within distance dist_min + # + dist_min = 25. + # - maximum fraction of the cluster signal allowed to be inside the leading hit + # - below this fraction we assume that we are dealing with a hit cluster + hit_frac_cut = 0.75 + found_cluster = False + + # loop starting from the highest hit + # + xcluster = np.zeros(3) + for i in range(len(indices)-1): + index = indices[i] + xhit_max = xhit[index] + # the seed hit has signal xmax + # + nmax = nhit[index] + + # we have reached the zero-signal PMTs..... no cluster was found, so leave this loop + # + if nmax == 0: + break + + # n_around indicates the amount of signal that is found around the seed hit + # + n_around = 0 + + # check the other hits.... + # + xcluster = nhit[index]*xhit[index] + for j in range(i+1, len(indices)): + jindex = indices[j] + xtest = xhit[jindex] + dist = lpf_dist(xhit_max, xtest) + # consider the hit if it is near the seed hit + # + if dist1000: + # print('WARNING:: result = ',result) + r_new = np.sqrt((xstart[0]-(xfit[0] - result[0]))**2 + (xstart[1] - (xfit[1] - result[1]))**2) + + weight = 0.5*(1 + np.math.erf((r_max-r_new)/1.0)) + # weight = 1 + # update fit result + #if weight<1e-3: + # print('ping ping ping.... weight =',weight) + ###weight = 1.0 + xfit[0] = xfit[0] - result[0]*weight + xfit[1] = xfit[1] - result[1]*weight + + #weight = np.exp(-lpf_dist(xstart,xfit)/25.) + r0 = r0 - result[2]*weight + + gamma = gamma - result[3]*weight + + # and store the intermediate results + # + xiter[lpf_iter + 1][0] = xfit[0] + xiter[lpf_iter + 1][1] = xfit[1] + xiter[lpf_iter + 1][2] = r0 + xiter[lpf_iter + 1][3] = lpf_lnlike_r(xhit, nhit, xfit, r0, gamma) + xiter[lpf_iter + 1][4] = n_in_fit + xiter[lpf_iter + 1][5] = gamma + + + # if (lpf_iter>=5) and (abs(result[0]) < 0.01) and (abs(result[1]) < 0.01) or (r0<0): # if position no longer changes -> terminate loop + if ((lpf_iter>1) and abs(result[0]) < lpf_min_position_cor) and (abs(result[1]) < lpf_min_position_cor) or (r0<0): # if position no longer changes -> terminate loop + break + + # calculate the number of uv photons + # + nuv = 4 * np.pi * r0 / area_sensor / xhit[0][2] + + # store the fit results + # + fit_result = np.zeros(7) + fit_result[0] = xfit[0] + fit_result[1] = xfit[1] + fit_result[2] = nuv + + logl = 0 + for ih in range(len(nhit)): + if nhit[ih] > 0: + logprob = lpf_hit_prob(xhit[ih], nhit[ih], xfit, r0, gamma) + if logprob>log_pmin: + logl = logl + lpf_hit_prob(xhit[ih], nhit[ih], xfit, r0, gamma) + + fit_result[3] = logl # lpf_lnlike_r(xhit, nhit, xfit, r0) + fit_result[4] = n_in_fit + fit_result[5] = gamma + fit_result[6] = r0 + + + return fit_result, xiter, hit_is_used + + +@njit +def lpf_dist(x0, x1): + d2 = (x0[0] - x1[0]) ** 2 + (x0[1] - x1[1]) ** 2 + (x0[2] - x1[2]) ** 2 + return np.sqrt(d2) + + +@njit +def lpf_nexp(xi, xf, r0, gamma): + """ + Calculate the expected number of p.e. for a lightsensor + + :param xi: sensor position + :param xf: assumed hit position + :param r0: assumed nuumber of UV photons (re-normalized) + :param gamma: assumed constant background + :return nexp: number of expected photons + + A.P. Colijn + """ + delta = lpf_dist(xi, xf) + nexp = r0 / delta ** 3 + gamma + + return nexp + +@njit +def lpf_f(i, xi, ni, xf, r0, gamma): + """ + Calculate the minimizer functions + :param i: + 0=F0 (x) + 1=F1 (y) + 2=F2 (r0) + 3=F3 (gamma) + :param xi: sensor position + :param ni: hits for the sensor + :param xf: assumed fit position + :param r0: assumed number of UV photons (re-normalized) + :param gamma: assumed constant background + :return f: function value + """ + + f = 0 + if i < 2: + f = -3 * (xf[i] - xi[i]) * (lpf_nexp(xi, xf, r0, gamma) - ni) / lpf_dist(xi, xf) ** 2 + elif i == 2: + f = (lpf_nexp(xi, xf, r0, gamma) - ni) / r0 + elif i == 3: + f = 1. - ni / lpf_nexp(xi, xf, r0, gamma) + + return f + +@njit +def lpf_deriv_f(i, j, xi, ni, xf, r0, gamma): + """ + Derivatives of the minimizer functions + + :param i: + 0=F0 + 1=F1 + 2=F2 + :param j: + 0=x + 1=y + z=r0 + :param xi: hit position + :param ni: number of hits + :param xf: fit position + :param r0: number of photons + :return: + """ + + d = lpf_dist(xi, xf) + n0 = lpf_nexp(xi, xf, r0, gamma) + + deriv = 0 + if i == 0: + dx = xf[0] - xi[0] + if j == 0: # dF0/dx + deriv = -3 * (n0 - ni) / d ** 2 + deriv = deriv - 3 * dx * (n0 - ni) * lpf_deriv_dist_min2(0, xi, xf) + deriv = deriv - 3 * dx * lpf_deriv_n(0, xi, xf, r0, gamma) / d ** 2 + elif j == 1: # dF0/dy + deriv = - 3 * dx * (n0 - ni) * lpf_deriv_dist_min2(1, xi, xf) + deriv = deriv - 3 * dx * lpf_deriv_n(1, xi, xf, r0, gamma) / d ** 2 + elif j == 2: # dF0/dr0 + deriv = -3 * dx * lpf_deriv_n(2, xi, xf, r0, gamma) / d ** 2 + elif j == 3: #dF0/dgamma + deriv = -3 * dx * lpf_deriv_n(3, xi, xf, r0, gamma) / d ** 2 + elif i == 1: + dy = xf[1] - xi[1] + if j == 0: # dF1/dx + deriv = - 3 * dy * (n0 - ni) * lpf_deriv_dist_min2(0, xi, xf) + deriv = deriv - 3 * dy * lpf_deriv_n(0, xi, xf, r0, gamma) / d ** 2 + elif j == 1: # dF1/dy + deriv = -3 * (n0 - ni) / d ** 2 + deriv = deriv - 3 * dy * (n0 - ni) * lpf_deriv_dist_min2(1, xi, xf) + deriv = deriv - 3 * dy * lpf_deriv_n(1, xi, xf, r0, gamma) / d ** 2 + elif j == 2: # dF1/dr0 + deriv = -3 * dy * lpf_deriv_n(2, xi, xf, r0, gamma) / d ** 2 + elif j == 3: #dF0/dgamma + deriv = -3 * dy * lpf_deriv_n(3, xi, xf, r0, gamma) / d ** 2 + elif i == 2: + #if j == 0: + # deriv = lpf_deriv_n(0, xi, xf, r0) + #elif j == 1: + # deriv = lpf_deriv_n(1, xi, xf, r0) + #elif j == 2: + # deriv = lpf_deriv_n(2, xi, xf, r0) + if j == 0: + deriv = lpf_deriv_n(0, xi, xf, r0, gamma) / r0 + elif j == 1: + deriv = lpf_deriv_n(1, xi, xf, r0, gamma) / r0 + elif j == 2: + deriv = lpf_deriv_n(2, xi, xf, r0, gamma) / r0 - (n0 - ni) / r0 ** 2 + elif j == 3: + deriv = lpf_deriv_n(3, xi, xf, r0, gamma) / r0 + elif i == 3: + deriv = ni * lpf_deriv_n(j, xi, xf, r0, gamma) / n0 ** 2 + else: + deriv = 0. + + + return deriv + + +@njit +def lpf_deriv_n(i, xi, xf, r0, gamma): + """ + Derivative of n wrt to fit parameters + + :param i: + 0=x + 1=y + 2=r0 + 3=gamma + :param xi: hit position + :param xf: fit position + :param r0: number of photons + :param gamma: constant background term + :return: dn/di + """ + + if i < 2: + deriv = -3 * lpf_nexp(xi, xf, r0, gamma) * (xf[i] - xi[i]) / lpf_dist(xi, xf) ** 2 + elif i == 2: + deriv = lpf_nexp(xi, xf, r0, gamma) / r0 + elif i == 3: + deriv = 1. + else: + deriv = 0. + + return deriv + + +@njit +def lpf_deriv_dist(i, xi, xf): + """ + Derivative of distance wrt fit parameters + + :param i: 0=x + 1=y + 2=r0 + :param xi: hit position + :param xf: fit position + :return: dDist/di + """ + + if i < 2: + deriv = (xf[i] - xi[i]) / lpf_dist(xi, xf) + else: + deriv = 0.0 + + return deriv + + +@njit +def lpf_deriv_dist_min2(i, xi, xf): + """ + Derivative of 1/dist**2 + + :param i: 0=x + 1=y + 2=r0 + :param xi: hit position + :param xf: fit position + :return: d(1/Dist**2)/di + """ + + deriv = 0.0 + if i < 2: + d = lpf_dist(xi, xf) + deriv = -(2 / d ** 3) * lpf_deriv_dist(i, xi, xf) + + return deriv + + +# --------------------------------------------------------------------------------------- # +def lpf_event_display(xhit, nhit, fit_result, hit_is_used, xiter, **kwargs): + """ + Event display + + :param xhit: + :param nhit: + :param fit_result: + :param hit-is_used: + :param xiter: + :param kwargs: + :return: + """ + + xhit = np.array(xhit) + nhit = np.array(nhit) + + plot_range = kwargs.pop('range', None) + zoom = kwargs.pop('zoom',-1) + nbins = kwargs.pop('nbins', 15) + + if plot_range == 'None': + plot_range = ((0, 100), (0, 100)) + + if zoom > 0: + plot_range = ((fit_result[0]-zoom/2,fit_result[0]+zoom/2),(fit_result[1]-zoom/2,fit_result[1]+zoom/2)) + + print("Reconstruction::lpf_event_display() ") + fig = plt.figure(figsize=(16, 6)) + + gs = GridSpec(2, 2, figure=fig) + ax0 = plt.subplot(gs.new_subplotspec((0, 0), rowspan=2)) + ax1 = plt.subplot(gs.new_subplotspec((0, 1), rowspan=1)) + ax2 = plt.subplot(gs.new_subplotspec((1, 1), rowspan=1)) + + # ax1 = plt.subplot(222) + # ax2 = plt.subplot(224) + + # fig, ax0 = plt.subplots(nrows=1) + + # make a list of the intermediate steps + xp = [] + yp = [] + nn = [] + iiter = [] + for i in range(len(xiter)): + if xiter[i][3] != 0: + xp.append(xiter[i][0]) + yp.append(xiter[i][1]) + nn.append(xiter[i][2]) + iiter.append(i) + xp = np.array(xp) + yp = np.array(yp) + nn = np.array(nn) + iiter = np.array(iiter) + + niter = len(xp) + ax1.plot(iiter, xp-fit_result[0]) + ax1.plot(iiter, yp-fit_result[1]) + ax1.set_xlabel('iteration') + ax1.set_ylabel('difference between position and fit result') + ax2.plot(iiter, nn) + ax2.set_xlabel('iteration') + ax2.set_ylabel('number of photons') + + # make these smaller to increase the resolution + dx, dy = (plot_range[0][1]-plot_range[0][0])/400,(plot_range[1][1]-plot_range[1][0])/400 + + # generate 2 2d grids for the x & y bounds + x = np.arange(plot_range[0][0], plot_range[0][1], dx) + y = np.arange(plot_range[1][0], plot_range[1][1], dy) + z = np.zeros((len(y), len(x))) + + #print(x,y) + for i in range(len(x)): + for j in range(len(y)): + xx = x[i] + yy = y[j] + xff = np.array([xx,yy,0]) + #print('xfit =',xff,' r0 =',xiter[niter-1][2]) + z[j][i] = lpf_lnlike_plot(xhit, nhit, hit_is_used, xff, xiter[niter-1][2], xiter[niter-1][5]) + + # z = z[:-1, :-1] + levels = MaxNLocator(nbins=nbins).tick_values(z.min(), z.max()) + + # cmap = plt.get_cmap('afmhot') + cmap = plt.get_cmap('PiYG') + + # norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True) + + # ax0 = fig.gca() + + cf = ax0.contourf(x + dx / 2., y + dy / 2., z, levels=levels, cmap=cmap) + fig.colorbar(cf, ax=ax0) + title_string = 'x = {:8.2f} y = {:8.2f} r0= {:8.2f}'.format(fit_result[0], fit_result[1], fit_result[2]) + ax0.set_title(title_string) + + # add the light detectors + mx_eff = -1 + for ih in range(len(nhit)): + if nhit[ih] > mx_eff: + mx_eff = nhit[ih] + + for ih in range(len(nhit)): + # draw location of SiPM + xs = xhit[ih] + + # plot sensor only if in range + if (xs[0] > plot_range[0][0]) & (xs[0] < plot_range[0][1]) & \ + (xs[1] > plot_range[1][0]) & (xs[1] < plot_range[1][1]): + #dx = nhit[ih] / mx_eff *10 + rr = (nhit[ih]+0.25) / mx_eff *10 + #sq = plt.Rectangle(xy=(xs[0] - dx / 2, xs[1] - dx / 2), + # height=dx, + # width=dx, + # fill=False, color='red') + if nhit[ih]>0: + if hit_is_used[ih]: + color = 'red' + else: + color = 'black' + sq = plt.Circle(xy=(xs[0], xs[1]), + radius=rr, + fill=False, color=color) + else: + sq = plt.Circle(xy=(xs[0], xs[1]), + radius=rr, + fill=False, color='black') + + ax0.add_artist(sq) + # write number of detected photons + ### txs = str(nhit[ih]) + ### ax0.text(xs[0] + dx / 2 + 2.5, xs[1], txs, color='red') + + # initial position + ax0.plot(xiter[0][0], xiter[0][1], 'o', markersize=10, color='cyan') + ax0.plot(xp, yp, 'w-o', markersize=5) + + # true position + # plt.plot(self.sim.get_x0()[0], self.sim.get_x0()[1], 'x', markersize=14, color='cyan') + # reconstructed position + # if abs(self.fdata['xr']) < 100: + ax0.plot(fit_result[0], fit_result[1], 'wo', markersize=10) + ax0.set_xlabel('x (cm)', fontsize=18) + ax0.set_ylabel('y (cm)', fontsize=18) + ax0.set_xlim([plot_range[0][0],plot_range[0][1]]) + ax0.set_ylim([plot_range[1][0],plot_range[1][1]]) + plt.show() + + istat = int(input("Type: 0 to continue, 1 to make pdf, 2 to quit....")) + + if istat == 1: + fname = 'event.pdf' + fig.savefig(fname) + fname = 'event.png' + fig.savefig(fname) + + clear_output() + return istat + + +#------------------------------ +def difference(pmt_pos,nhit,xml,yml,r0,gamma): + """ + Give a plot of the expected number of hits and the detected number of hits + Also shows the impact of the gamma factor on the expected number of hits + """ + + #different lists used + difference_n = [] + x_lijst = [] + y_lijst = [] + n_exp_lijst = [] + n_exp_gamma_lijst = [] + x_nexp_lijst = [] + + #loops through all the PMTs + for iii in range(len(pmt_pos)): + x_nexp_lijst.append(iii) + + # calculate expected number of hits + n_exp = lpf_nexp(pmt_pos[iii],(xml,yml,0),r0,gamma) + difference = (((n_exp-nhit[iii])))/(n_exp**0.5) + difference_n.append(difference) + n_exp_lijst.append(n_exp) + n_exp_gamma_lijst.append(n_exp - gamma) + x_lijst.append(pmt_pos[iii][0]) + y_lijst.append(pmt_pos[iii][1]) + + # bins the detected hits into clear datapoints with error values + nbins5 = 100 + n15, _____ = np.histogram(x_nexp_lijst, bins=nbins5) + sy15, _____ = np.histogram(x_nexp_lijst, bins=nbins5, weights=(nhit)) + sy215, _____ = np.histogram(x_nexp_lijst, bins=nbins5, weights=(nhit)*(nhit)) + mean15 = sy15 / n15 + std15 = np.sqrt(sy215/n15 - mean15*mean15) + + # plots the expected number of photons with and without gamma + plt.errorbar((_____[1:] + _____[:-1])/2, mean15, yerr=std15, fmt='ko',zorder=10) + plt.plot(n_exp_lijst,label= "Fitted hits with gamma",linewidth = 2,color = 'red',zorder=5) + plt.plot(n_exp_gamma_lijst, '--',label= "Fitted hits without gamma",linewidth = 2,color = 'red',zorder=0) + plt.xlabel("PMT number",fontsize=14) + plt.ylabel("number of hits",fontsize=14) + plt.legend(loc="upper right") + plt.yscale('log') + plt.xlim(0,250) + plt.ylim(0,10000) + plt.show() From b550b2491333bf8d2f9134380b281abfd91e83dd Mon Sep 17 00:00:00 2001 From: cfuselli Date: Sun, 26 Mar 2023 23:02:28 +0000 Subject: [PATCH 2/4] restyle and add ezra errorcalc --- straxen/plugins/peaks/peak_positions_lpf.py | 360 +++++--------------- 1 file changed, 90 insertions(+), 270 deletions(-) diff --git a/straxen/plugins/peaks/peak_positions_lpf.py b/straxen/plugins/peaks/peak_positions_lpf.py index 437e3f08d..be7f90911 100644 --- a/straxen/plugins/peaks/peak_positions_lpf.py +++ b/straxen/plugins/peaks/peak_positions_lpf.py @@ -1,33 +1,31 @@ -import strax -import straxen -from straxen.plugins.peaks._peak_positions_base import PeakPositionsBaseNT - - -export, __all__ = strax.exporter() - import strax import numpy as np import numba import straxen +from numba import njit +export, __all__ = strax.exporter() + +# A.P. Colijn - Feb. 2021 +# Parameters for the fit +lpf_iter_max = 100 +lpf_min_position_cor = 0.01 # convergence criterium for position minimizer +lpf_npar = 4 # number of fit parameters +# x y r0 gamma +lpf_fix_par = np.array([False, False, False, False]) # which parameters to fix @export class PeakPositionsLPF(strax.Plugin): """ + LPF = Linearized Position Fit Computes S2 position with linearized log-likelihood fit. Returns xy position and the parameters of the fit. - See LinPosFit https://github.com/XAMS-nikhef/LinPosFit/ from A.P. Colijn + See LinPosFit https://github.com/XAMS-nikhef/LinPosFit/ + from A.P. Colijn (Feb 2021) """ - depends_on=('peaks','peak_basics') + depends_on =('peaks','peak_basics') provides = 'peak_positions_lpf' - rechunk_on_save=False - dtype= [('xml',np.float), - ('yml',np.float), - ('r0',np.float), - ('gamma', np.float), - ('logl',np.float), - ('n',np.int)] - dtype += strax.time_fields + rechunk_on_save = False n_top_pmts = straxen.URLConfig(default=straxen.n_top_pmts, infer_type=False, help="Number of top PMTs") @@ -35,6 +33,21 @@ class PeakPositionsLPF(strax.Plugin): pmt_to_lxe = straxen.URLConfig(default=7, infer_type=False, help="Distance between the PMTs and the liquid interface") + def infer_dtype(self): + + dtype = [('x_lpf',np.float32, 'Reconstructed x (cm) position, LinPosFit'), + ('y_lpf',np.float32, 'Reconstructed y (cm) position, LinPosFit'), + ('lpf_logl',np.float32, 'LinPosFit LogLikelihood value'), + ('lpf_r0',np.float32, 'LinPosFit r0 parameter, reduced rate'), + ('lpf_gamma', np.float32, 'LinPosFit gamma parameter, reflection'), + ('lpf_n',np.float32, 'LinPosFit parameter n, number of photons generated'), + ('lpf_err_x',np.float32, 'LinPosFit error on x (cm), 95p conficence interval'), + ('lpf_err_y',np.float32, 'LinPosFit error on y (cm), 95p conficence interval') + ] + + dtype += strax.time_fields + return dtype + def setup(self,): inch = 2.54 # cm pmt_pos = straxen.pmt_positions() @@ -42,57 +55,51 @@ def setup(self,): self.pmt_surface=(3*inch)**2*np.pi/4. def compute(self,peaks): - res = np.zeros(len(peaks),self.dtype) - for ix, p in enumerate(peaks): - if p['type']!=2: - #Only reconstruct s2 peaks. We do need to set the time of the peaks - res[ix]['time'] = p['time'] - res[ix]['endtime'] = p['endtime'] - continue - try: - #Some really small single electron s2s fail the minimization - fit_result,_,_ = lpf_execute(self.pmt_pos[:self.n_top_pmts],p['area_per_channel'][:self.n_top_pmts],self.pmt_surface) - - res[ix]['time'] = p['time'] - res[ix]['endtime'] = p['endtime'] - res[ix]['xml'] = fit_result[0] - res[ix]['yml'] = fit_result[1] - res[ix]['r0'] = fit_result[2] - res[ix]['logl'] = fit_result[3] - res[ix]['n'] = fit_result[4] - res[ix]['gamma'] = fit_result[5] - - except Exception as e: - print(e) - res[ix]['time'] = p['time'] - res[ix]['endtime'] = p['endtime'] - continue - - return res - - + result = np.ones(len(peaks), dtype=self.dtype) + + result['time'], result['endtime'] = peaks['time'], strax.endtime(peaks) + result['x_lpf'] *= float('nan') + result['y_lpf'] *= float('nan') + result['lpf_r0'] *= float('nan') + result['lpf_logl'] *= float('nan') + result['lpf_n'] *= float('nan') + result['lpf_gamma'] *= float('nan') + result['lpf_err_x'] *= float('nan') + result['lpf_err_y'] *= float('nan') + for ix, p in enumerate(peaks): + if p['type']==2: + # Only reconstruct s2 peaks. We do need to set the time of the peaks + try: + # Some really small single electron s2s fail the minimization + fit_result, xiter, user_hits = lpf_execute(self.pmt_pos[:self.n_top_pmts], + p['area_per_channel'][:self.n_top_pmts], + self.pmt_surface) + + err_x, err_y = lpf_deriv_95(self.pmt_pos[:self.n_top_pmts], + p['area_per_channel'][:self.n_top_pmts], + fit_result, + used_hits) + + result[ix]['time'] = p['time'] + result[ix]['endtime'] = p['endtime'] + result[ix]['x_lpf'] = fit_result[0] + result[ix]['y_lpf'] = fit_result[1] + result[ix]['lpf_r0'] = fit_result[2] + result[ix]['lpf_logl'] = fit_result[3] + result[ix]['lpf_n'] = fit_result[4] + result[ix]['lpf_gamma'] = fit_result[5] + result[ix]['lpf_err_x'] = err_x + result[ix]['lpf_err_y'] = err_y + + except Exception as e: + # Sometimes inverting the matrix fails.. + # Let's just leave it as a NaN + pass + + return result -import numpy as np -from numba import njit - -# from matplotlib.colors import BoundaryNorm -import matplotlib.pyplot as plt -from matplotlib.gridspec import GridSpec -from matplotlib.ticker import MaxNLocator -from IPython.display import clear_output - -# lpf = linearized position fit -# -# A.P. Colijn - Feb. 2021 -lpf_iter_max = 100 -lpf_min_position_cor = 0.01 # convergence criterium for position minimizer -lpf_npar = 4 # number of fit parameters -# x y r0 gamma -lpf_fix_par = np.array([False, False, False, False]) # which parameters to fix -# -debug = False def lpf_execute(xhit, nhit, area_sensor): @@ -235,10 +242,7 @@ def lpf_initialize(xhit, nhit): else: # we calculate the center-of-gravity as initial position xhit_max = xcluster / (nmax+n_around) - if debug: - print('--- lpf_initialize --- hit fraction = ',hit_frac) - print('--- lpf_initialize --- found_cluster = ',found_cluster) - + # refine the determination of the initial position # logl_min = 1e12 # xtemp = np.zeros(3) @@ -603,212 +607,28 @@ def lpf_deriv_dist_min2(i, xi, xf): return deriv - -# --------------------------------------------------------------------------------------- # -def lpf_event_display(xhit, nhit, fit_result, hit_is_used, xiter, **kwargs): +@njit +def lpf_deriv_95(xhit, nhit, fit_result, hit_is_used, **kwargs): """ - Event display - - :param xhit: - :param nhit: - :param fit_result: - :param hit-is_used: - :param xiter: - :param kwargs: - :return: + Error estimation for LinPosFit + Bachelor thesis: Ezra de Cleen, Nikhef, June 2022 """ + + deriv_x = 0 + deriv_y = 0 - xhit = np.array(xhit) - nhit = np.array(nhit) + for isensor in range(len(nhit)): - plot_range = kwargs.pop('range', None) - zoom = kwargs.pop('zoom',-1) - nbins = kwargs.pop('nbins', 15) + xff = np.array([fit_result[0],fit_result[1],0]) + + deriv_x += lpf_deriv_f(0, 0, xhit[isensor], nhit[isensor], xff, fit_result[2], fit_result[5]) - if plot_range == 'None': - plot_range = ((0, 100), (0, 100)) + xff = np.array([fit_result[0],fit_result[1],0]) - if zoom > 0: - plot_range = ((fit_result[0]-zoom/2,fit_result[0]+zoom/2),(fit_result[1]-zoom/2,fit_result[1]+zoom/2)) - - print("Reconstruction::lpf_event_display() ") - fig = plt.figure(figsize=(16, 6)) - - gs = GridSpec(2, 2, figure=fig) - ax0 = plt.subplot(gs.new_subplotspec((0, 0), rowspan=2)) - ax1 = plt.subplot(gs.new_subplotspec((0, 1), rowspan=1)) - ax2 = plt.subplot(gs.new_subplotspec((1, 1), rowspan=1)) - - # ax1 = plt.subplot(222) - # ax2 = plt.subplot(224) - - # fig, ax0 = plt.subplots(nrows=1) - - # make a list of the intermediate steps - xp = [] - yp = [] - nn = [] - iiter = [] - for i in range(len(xiter)): - if xiter[i][3] != 0: - xp.append(xiter[i][0]) - yp.append(xiter[i][1]) - nn.append(xiter[i][2]) - iiter.append(i) - xp = np.array(xp) - yp = np.array(yp) - nn = np.array(nn) - iiter = np.array(iiter) - - niter = len(xp) - ax1.plot(iiter, xp-fit_result[0]) - ax1.plot(iiter, yp-fit_result[1]) - ax1.set_xlabel('iteration') - ax1.set_ylabel('difference between position and fit result') - ax2.plot(iiter, nn) - ax2.set_xlabel('iteration') - ax2.set_ylabel('number of photons') - - # make these smaller to increase the resolution - dx, dy = (plot_range[0][1]-plot_range[0][0])/400,(plot_range[1][1]-plot_range[1][0])/400 - - # generate 2 2d grids for the x & y bounds - x = np.arange(plot_range[0][0], plot_range[0][1], dx) - y = np.arange(plot_range[1][0], plot_range[1][1], dy) - z = np.zeros((len(y), len(x))) - - #print(x,y) - for i in range(len(x)): - for j in range(len(y)): - xx = x[i] - yy = y[j] - xff = np.array([xx,yy,0]) - #print('xfit =',xff,' r0 =',xiter[niter-1][2]) - z[j][i] = lpf_lnlike_plot(xhit, nhit, hit_is_used, xff, xiter[niter-1][2], xiter[niter-1][5]) - - # z = z[:-1, :-1] - levels = MaxNLocator(nbins=nbins).tick_values(z.min(), z.max()) - - # cmap = plt.get_cmap('afmhot') - cmap = plt.get_cmap('PiYG') - - # norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True) - - # ax0 = fig.gca() - - cf = ax0.contourf(x + dx / 2., y + dy / 2., z, levels=levels, cmap=cmap) - fig.colorbar(cf, ax=ax0) - title_string = 'x = {:8.2f} y = {:8.2f} r0= {:8.2f}'.format(fit_result[0], fit_result[1], fit_result[2]) - ax0.set_title(title_string) - - # add the light detectors - mx_eff = -1 - for ih in range(len(nhit)): - if nhit[ih] > mx_eff: - mx_eff = nhit[ih] + deriv_y += lpf_deriv_f(1, 1, xhit[isensor], nhit[isensor], xff, fit_result[2], fit_result[5]) + + x_95 = 1*(1/np.sqrt(deriv_x)) + y_95 = 1*(1/np.sqrt(deriv_y)) - for ih in range(len(nhit)): - # draw location of SiPM - xs = xhit[ih] - - # plot sensor only if in range - if (xs[0] > plot_range[0][0]) & (xs[0] < plot_range[0][1]) & \ - (xs[1] > plot_range[1][0]) & (xs[1] < plot_range[1][1]): - #dx = nhit[ih] / mx_eff *10 - rr = (nhit[ih]+0.25) / mx_eff *10 - #sq = plt.Rectangle(xy=(xs[0] - dx / 2, xs[1] - dx / 2), - # height=dx, - # width=dx, - # fill=False, color='red') - if nhit[ih]>0: - if hit_is_used[ih]: - color = 'red' - else: - color = 'black' - sq = plt.Circle(xy=(xs[0], xs[1]), - radius=rr, - fill=False, color=color) - else: - sq = plt.Circle(xy=(xs[0], xs[1]), - radius=rr, - fill=False, color='black') - - ax0.add_artist(sq) - # write number of detected photons - ### txs = str(nhit[ih]) - ### ax0.text(xs[0] + dx / 2 + 2.5, xs[1], txs, color='red') - - # initial position - ax0.plot(xiter[0][0], xiter[0][1], 'o', markersize=10, color='cyan') - ax0.plot(xp, yp, 'w-o', markersize=5) - - # true position - # plt.plot(self.sim.get_x0()[0], self.sim.get_x0()[1], 'x', markersize=14, color='cyan') - # reconstructed position - # if abs(self.fdata['xr']) < 100: - ax0.plot(fit_result[0], fit_result[1], 'wo', markersize=10) - ax0.set_xlabel('x (cm)', fontsize=18) - ax0.set_ylabel('y (cm)', fontsize=18) - ax0.set_xlim([plot_range[0][0],plot_range[0][1]]) - ax0.set_ylim([plot_range[1][0],plot_range[1][1]]) - plt.show() - - istat = int(input("Type: 0 to continue, 1 to make pdf, 2 to quit....")) - - if istat == 1: - fname = 'event.pdf' - fig.savefig(fname) - fname = 'event.png' - fig.savefig(fname) - - clear_output() - return istat - - -#------------------------------ -def difference(pmt_pos,nhit,xml,yml,r0,gamma): - """ - Give a plot of the expected number of hits and the detected number of hits - Also shows the impact of the gamma factor on the expected number of hits - """ - - #different lists used - difference_n = [] - x_lijst = [] - y_lijst = [] - n_exp_lijst = [] - n_exp_gamma_lijst = [] - x_nexp_lijst = [] - - #loops through all the PMTs - for iii in range(len(pmt_pos)): - x_nexp_lijst.append(iii) - # calculate expected number of hits - n_exp = lpf_nexp(pmt_pos[iii],(xml,yml,0),r0,gamma) - difference = (((n_exp-nhit[iii])))/(n_exp**0.5) - difference_n.append(difference) - n_exp_lijst.append(n_exp) - n_exp_gamma_lijst.append(n_exp - gamma) - x_lijst.append(pmt_pos[iii][0]) - y_lijst.append(pmt_pos[iii][1]) - - # bins the detected hits into clear datapoints with error values - nbins5 = 100 - n15, _____ = np.histogram(x_nexp_lijst, bins=nbins5) - sy15, _____ = np.histogram(x_nexp_lijst, bins=nbins5, weights=(nhit)) - sy215, _____ = np.histogram(x_nexp_lijst, bins=nbins5, weights=(nhit)*(nhit)) - mean15 = sy15 / n15 - std15 = np.sqrt(sy215/n15 - mean15*mean15) - - # plots the expected number of photons with and without gamma - plt.errorbar((_____[1:] + _____[:-1])/2, mean15, yerr=std15, fmt='ko',zorder=10) - plt.plot(n_exp_lijst,label= "Fitted hits with gamma",linewidth = 2,color = 'red',zorder=5) - plt.plot(n_exp_gamma_lijst, '--',label= "Fitted hits without gamma",linewidth = 2,color = 'red',zorder=0) - plt.xlabel("PMT number",fontsize=14) - plt.ylabel("number of hits",fontsize=14) - plt.legend(loc="upper right") - plt.yscale('log') - plt.xlim(0,250) - plt.ylim(0,10000) - plt.show() + return x_95, y_95 \ No newline at end of file From 5f38fc0ab4645326d8dbea4ba9c48c5a56e08b48 Mon Sep 17 00:00:00 2001 From: cfuselli Date: Sun, 26 Mar 2023 23:05:54 +0000 Subject: [PATCH 3/4] restyle and add ezra errorcalc --- straxen/plugins/peaks/peak_positions_lpf.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/straxen/plugins/peaks/peak_positions_lpf.py b/straxen/plugins/peaks/peak_positions_lpf.py index be7f90911..e399c6e77 100644 --- a/straxen/plugins/peaks/peak_positions_lpf.py +++ b/straxen/plugins/peaks/peak_positions_lpf.py @@ -72,18 +72,17 @@ def compute(self,peaks): if p['type']==2: # Only reconstruct s2 peaks. We do need to set the time of the peaks try: - # Some really small single electron s2s fail the minimization + # Execute linearised position reconstruction fit fit_result, xiter, user_hits = lpf_execute(self.pmt_pos[:self.n_top_pmts], p['area_per_channel'][:self.n_top_pmts], self.pmt_surface) + # Error estimation function err_x, err_y = lpf_deriv_95(self.pmt_pos[:self.n_top_pmts], p['area_per_channel'][:self.n_top_pmts], fit_result, used_hits) - result[ix]['time'] = p['time'] - result[ix]['endtime'] = p['endtime'] result[ix]['x_lpf'] = fit_result[0] result[ix]['y_lpf'] = fit_result[1] result[ix]['lpf_r0'] = fit_result[2] From b4b14ab79ae8d6a962398af44e2e5f91a38f08fc Mon Sep 17 00:00:00 2001 From: cfuselli Date: Thu, 13 Apr 2023 13:56:40 +0000 Subject: [PATCH 4/4] it works --- straxen/plugins/events/__init__.py | 3 + .../plugins/events/event_s2_position_lpf.py | 109 ++++++++++++++++++ straxen/plugins/peaks/peak_positions_lpf.py | 69 ++++++----- 3 files changed, 153 insertions(+), 28 deletions(-) create mode 100644 straxen/plugins/events/event_s2_position_lpf.py diff --git a/straxen/plugins/events/__init__.py b/straxen/plugins/events/__init__.py index 0079976a5..44501f43f 100644 --- a/straxen/plugins/events/__init__.py +++ b/straxen/plugins/events/__init__.py @@ -55,6 +55,9 @@ from . import event_s2_position_gcn from .event_s2_position_gcn import * +from . import event_s2_position_lpf +from .event_s2_position_lpf import * + from . import local_minimum_info from .local_minimum_info import * diff --git a/straxen/plugins/events/event_s2_position_lpf.py b/straxen/plugins/events/event_s2_position_lpf.py new file mode 100644 index 000000000..db0a67104 --- /dev/null +++ b/straxen/plugins/events/event_s2_position_lpf.py @@ -0,0 +1,109 @@ +import strax +import straxen +import numpy as np +from warnings import warn +export, __all__ = strax.exporter() + +from straxen.plugins.peaks.peak_positions_lpf import lpf_execute + +@export +class EventS2PositionLPF(strax.Plugin): + """ + LinPosFit pluging for S2 position reconstruction at event level + """ + __version__ = '0.0.0' + depends_on = ('event_area_per_channel', 'event_basics') + provides = "event_s2_position_lpf" + + algorithm = 'lpf' + compressor = 'zstd' + parallel = True # can set to "process" after #82 + + min_reconstruction_area = straxen.URLConfig( + help='Skip reconstruction if area (PE) is less than this', + default=10, infer_type=False, ) + n_top_pmts = straxen.URLConfig( + default=straxen.n_top_pmts, infer_type=False, + help="Number of top PMTs") + pmt_to_lxe = straxen.URLConfig(default=7, infer_type=False, + help="Distance between the PMTs and the liquid interface") + + + def infer_dtype(self): + dtype = [] + for ptype in ['', 'alt_']: + for xy in ['x', 'y']: + dtype += [(f'event_{ptype}s2_{xy}_' + self.algorithm, np.float32, + f'Reconstructed {self.algorithm} {ptype}S2 {xy} position (cm), uncorrected'), + (f'event_{ptype}s2_err_{xy}_' + self.algorithm, np.float32, + f'Error on {self.algorithm} {ptype}S2 {xy} position (cm), 95p confidence interval'), + ] + + dtype += [(f'event_{ptype}s2_lpf_logl',np.float32, f'LinPosFit {ptype}S2 LogLikelihood value'), + (f'event_{ptype}s2_lpf_nuv',np.float32, f'LinPosFit {ptype}S2 parameter n, number of photons generated'), + (f'event_{ptype}s2_lpf_gamma', np.float32, f'LinPosFit {ptype}S2 gamma parameter, reflection'), + (f'event_{ptype}s2_lpf_n_in_fit',np.float32, f'LinPosFit {ptype}S2 number of PMTs included in the fit'), + (f'event_{ptype}s2_lpf_r0',np.float32, f'LinPosFit {ptype}S2 r0 parameter, reduced rate'), + ] + + dtype += strax.time_fields + return dtype + + def setup(self,): + inch = 2.54 # cm + pmt_pos = straxen.pmt_positions() + self.pmt_pos = list(zip(pmt_pos['x'].values,pmt_pos['y'].values,np.repeat(self.pmt_to_lxe, self.n_top_pmts))) + self.pmt_surface=(3*inch)**2*np.pi/4. + + + def compute(self, events): + + result = np.ones(len(events), dtype=self.dtype) + result['time'], result['endtime'] = events['time'], strax.endtime(events) + + for p_type in ['s2', 'alt_s2']: + result[f'event_{p_type}_x_lpf'] *= float('nan') + result[f'event_{p_type}_err_x_lpf'] *= float('nan') + result[f'event_{p_type}_y_lpf'] *= float('nan') + result[f'event_{p_type}_err_y_lpf'] *= float('nan') + result[f'event_{p_type}_lpf_logl'] *= float('nan') + result[f'event_{p_type}_lpf_nuv'] *= float('nan') + result[f'event_{p_type}_lpf_gamma'] *= float('nan') + result[f'event_{p_type}_lpf_n_in_fit'] *= float('nan') + result[f'event_{p_type}_lpf_r0'] *= float('nan') + + + for ix, ev in enumerate(events): + for p_type in ['s2', 'alt_s2']: + if ev[p_type+'_area'] > self.min_reconstruction_area: + + _top_pattern = ev[p_type + '_area_per_channel'][0:self.n_top_pmts] + + try: + + # Execute linearised position reconstruction fit + fit_result, xiter, used_hits = lpf_execute(self.pmt_pos[:self.n_top_pmts], + _top_pattern, + self.pmt_surface) + + result[ix][f'event_{p_type}_x_lpf'] = fit_result[0] + result[ix][f'event_{p_type}_y_lpf'] = fit_result[1] + result[ix][f'event_{p_type}_lpf_nuv'] = fit_result[2] + result[ix][f'event_{p_type}_lpf_logl'] = fit_result[3] + result[ix][f'event_{p_type}_lpf_n_in_fit'] = fit_result[4] + result[ix][f'event_{p_type}_lpf_gamma'] = fit_result[5] + result[ix][f'event_{p_type}_lpf_r0'] = fit_result[6] + result[ix][f'event_{p_type}_err_x_lpf'] = fit_result[7] + result[ix][f'event_{p_type}_err_y_lpf'] = fit_result[8] + + except Exception as e: + # Sometimes inverting the matrix fails.. + # Let's just leave it as a NaN + print(e) + pass + + return result + + + + diff --git a/straxen/plugins/peaks/peak_positions_lpf.py b/straxen/plugins/peaks/peak_positions_lpf.py index e399c6e77..cd78f2ccd 100644 --- a/straxen/plugins/peaks/peak_positions_lpf.py +++ b/straxen/plugins/peaks/peak_positions_lpf.py @@ -38,9 +38,10 @@ def infer_dtype(self): dtype = [('x_lpf',np.float32, 'Reconstructed x (cm) position, LinPosFit'), ('y_lpf',np.float32, 'Reconstructed y (cm) position, LinPosFit'), ('lpf_logl',np.float32, 'LinPosFit LogLikelihood value'), - ('lpf_r0',np.float32, 'LinPosFit r0 parameter, reduced rate'), + ('lpf_nuv',np.float32, 'LinPosFit parameter n, number of photons generated'), ('lpf_gamma', np.float32, 'LinPosFit gamma parameter, reflection'), - ('lpf_n',np.float32, 'LinPosFit parameter n, number of photons generated'), + ('lpf_n_in_fit',np.float32, 'LinPosFit number of PMTs included in the fit'), + ('lpf_r0',np.float32, 'LinPosFit r0 parameter, reduced rate'), ('lpf_err_x',np.float32, 'LinPosFit error on x (cm), 95p conficence interval'), ('lpf_err_y',np.float32, 'LinPosFit error on y (cm), 95p conficence interval') ] @@ -61,10 +62,11 @@ def compute(self,peaks): result['time'], result['endtime'] = peaks['time'], strax.endtime(peaks) result['x_lpf'] *= float('nan') result['y_lpf'] *= float('nan') - result['lpf_r0'] *= float('nan') + result['lpf_nuv'] *= float('nan') result['lpf_logl'] *= float('nan') - result['lpf_n'] *= float('nan') + result['lpf_n_in_fit'] *= float('nan') result['lpf_gamma'] *= float('nan') + result['lpf_r0'] *= float('nan') result['lpf_err_x'] *= float('nan') result['lpf_err_y'] *= float('nan') @@ -73,24 +75,19 @@ def compute(self,peaks): # Only reconstruct s2 peaks. We do need to set the time of the peaks try: # Execute linearised position reconstruction fit - fit_result, xiter, user_hits = lpf_execute(self.pmt_pos[:self.n_top_pmts], + fit_result, xiter, used_hits = lpf_execute(self.pmt_pos[:self.n_top_pmts], p['area_per_channel'][:self.n_top_pmts], self.pmt_surface) - # Error estimation function - err_x, err_y = lpf_deriv_95(self.pmt_pos[:self.n_top_pmts], - p['area_per_channel'][:self.n_top_pmts], - fit_result, - used_hits) - result[ix]['x_lpf'] = fit_result[0] result[ix]['y_lpf'] = fit_result[1] - result[ix]['lpf_r0'] = fit_result[2] + result[ix]['lpf_nuv'] = fit_result[2] result[ix]['lpf_logl'] = fit_result[3] - result[ix]['lpf_n'] = fit_result[4] + result[ix]['lpf_n_in_fit'] = fit_result[4] result[ix]['lpf_gamma'] = fit_result[5] - result[ix]['lpf_err_x'] = err_x - result[ix]['lpf_err_y'] = err_y + result[ix]['lpf_r0'] = fit_result[6] + result[ix]['lpf_err_x'] = fit_result[7] + result[ix]['lpf_err_y'] = fit_result[8] except Exception as e: # Sometimes inverting the matrix fails.. @@ -115,6 +112,7 @@ def lpf_execute(xhit, nhit, area_sensor): nhit = np.array(nhit) # minimize the likelihood fit_result, xiter, hit_is_used = lpf_minimize(xhit, nhit, area_sensor) + return fit_result, xiter, hit_is_used @njit @@ -390,7 +388,7 @@ def lpf_minimize(xhit, nhit, area_sensor): # store the fit results # - fit_result = np.zeros(7) + fit_result = np.zeros(9) fit_result[0] = xfit[0] fit_result[1] = xfit[1] fit_result[2] = nuv @@ -407,7 +405,11 @@ def lpf_minimize(xhit, nhit, area_sensor): fit_result[5] = gamma fit_result[6] = r0 + errors = lpf_deriv_95(xhit, nhit, hit_is_used, xfit, nuv, gamma) + fit_result[7] = errors[0] + fit_result[8] = errors[1] + return fit_result, xiter, hit_is_used @@ -607,27 +609,38 @@ def lpf_deriv_dist_min2(i, xi, xf): return deriv @njit -def lpf_deriv_95(xhit, nhit, fit_result, hit_is_used, **kwargs): +def lpf_deriv_95(xhit, nhit, hit_is_used, xfit, nuv, gamma): """ Error estimation for LinPosFit Bachelor thesis: Ezra de Cleen, Nikhef, June 2022 """ - deriv_x = 0 - deriv_y = 0 + # deriv_x = 0 + # deriv_y = 0 - for isensor in range(len(nhit)): + # for isensor in np.arange(len(nhit)): - xff = np.array([fit_result[0],fit_result[1],0]) + # xff = np.zeros(3) + # xff[0] = xfit[0] + # xff[1] = xfit[1] + # xff[2] = 0. - deriv_x += lpf_deriv_f(0, 0, xhit[isensor], nhit[isensor], xff, fit_result[2], fit_result[5]) + # # m[i][j] = m[i][j] + lpf_deriv_f(i, j, xhit[isensor], nhit[isensor], xfit, r0, gamma) - xff = np.array([fit_result[0],fit_result[1],0]) + # deriv_x = deriv_x + lpf_deriv_f(float(0), float(0), xhit[isensor], nhit[isensor], xff, nuv, gamma) + # deriv_y = deriv_y + lpf_deriv_f(float(1), float(1), xhit[isensor], nhit[isensor], xff, nuv, gamma) + + + # return x_95, y_95 - deriv_y += lpf_deriv_f(1, 1, xhit[isensor], nhit[isensor], xff, fit_result[2], fit_result[5]) - x_95 = 1*(1/np.sqrt(deriv_x)) - y_95 = 1*(1/np.sqrt(deriv_y)) + err = np.zeros(2) + + for i in np.arange(2): + for isensor in range(len(nhit)): + err[i] = err[i] + lpf_deriv_f(i, i, xhit[isensor], nhit[isensor], xfit, nuv, gamma) + + err[i] = 1/np.sqrt(err[i]) + + return err - - return x_95, y_95 \ No newline at end of file