add scripts for fitting arrival direction and energy

fitting/fit-energy.py

+import numpy as np
+import matplotlib.pyplot as plt
+
+
+# load data
+with np.load('../RawData/showers-A1-0.npz') as data:
+    vd = data['detector']
+    vc = data['showercore']
+    va = data['showeraxis']
+    logE = data['logE']
+    signal = data['signal'].sum(axis=-1)  # total signal
+    nb_events = len(vc)
+    nb_stations = len(vd)
+
+
+# for each event, calculate distances of detectors to shower core, plane and axis
+xd = np.repeat(np.expand_dims(vd, 0), nb_events, axis=0)
+xc = np.repeat(np.expand_dims(vc, 1), nb_stations, axis=1)
+xa = np.repeat(np.expand_dims(va, 1), nb_stations, axis=1)
+
+d1 = np.linalg.norm(xd - xc, axis=-1)  # distance to shower core
+d2 = np.sum((xd - xc) * xa, axis=-1)  # distance to shower plane
+d3 = np.linalg.norm(xd - xc - xa * d2[..., np.newaxis], axis=-1)  # distance to shower axis
+
+
+# fit lateral distribution to each event, S(r) = S1000 * (r/1000)**a
+logS1000 = np.zeros(nb_events)
+for i in range(nb_events):
+    m = ~np.isnan(signal[i])  # mask for triggered stations
+    m *= d3[i] > 200  # cut stations close to the shower core
+    x = np.log10(d3[i][m] / 1000)  # log10(r/1000m) distance to shower axis
+    y = np.log10(signal[i][m] - 100)  # log10(signal - noise level)
+    p = np.polyfit(x, y, 1)  # fit linear model
+    logS1000[i] = p[1]
+
+    # plt.figure()
+    # plt.scatter(x + 3, y)
+    # x = np.linspace(-1, 1)
+    # y = np.polyval(p, x)
+    # plt.plot(x + 3, y, 'r')
+    # plt.xlabel('$\log_{10}(r)$')
+    # plt.ylabel('$\log_{10}(S)$')
+    # plt.grid()
+
+
+# fit zenith angle dependency: logS1000 = logS38 + p * (cos(zenith)**2 - cos(38)**2)
+x, y, z = va.T
+zenith = np.pi / 2 - np.arctan2(z, (x * x + y * y) ** .5)
+cz = np.cos(zenith)**2 - np.cos(np.deg2rad(38))**2
+
+pcal1 = np.polyfit(cz, logS1000, 1)
+cic = np.polyval(pcal1, cz)
+logS38 = logS1000 / cic  # remove the zenith angle dependency
+
+plt.figure()
+plt.scatter(cz, logS1000)
+plt.plot(cz, cic, 'r')
+plt.xlabel(r'$\cos(\theta)^2 - \cos(38)^2$')
+plt.ylabel(r'$\log_{10}(S_{1000})$')
+plt.savefig('energy-calibration1.png')
+
+
+# fit energy dependency: logE = p0 + p1 * logS
+pcal2 = np.polyfit(logS38, logE, 1)
+logE_rec = np.polyval(pcal2, logS38)
+
+plt.figure()
+plt.scatter(logS38, logE)
+x = np.linspace(min(logS38), max(logS38))
+y = np.polyval(pcal2, x)
+plt.plot(x, y, 'r')
+plt.xlabel('$\log_{10}(S_{38})$')
+plt.ylabel('$\log_{10}(E / \mathrm{eV})$')
+plt.grid()
+plt.savefig('energy-calibration2.png')
+
+
+# evaluate resolution
+r = 10**(logE_rec - logE) - 1  # (E_rec - E_true) / E_true
+plt.figure()
+plt.hist(r, bins=np.linspace(-0.6, 0.6, 20))
+plt.xlabel('$(E_\mathrm{rec} - E_\mathrm{true}) / E_\mathrm{true}$')
+plt.ylabel('\#')
+plt.grid()
+plt.text(
+    0.95, 0.95,
+    'mean = %.3f\n std = %.3f' % (np.mean(r), np.std(r)),
+    transform=plt.gca().transAxes, ha='right', va='top',
+    bbox=dict(boxstyle='round', facecolor='white'))
+plt.savefig('energy-resolution.png')

fitting/fit-planewave.py

+import numpy as np
+import matplotlib.pyplot as plt
+
+
+# load data
+with np.load('../RawData/showers-A1-0.npz') as data:
+    vd = data['detector']
+    va = data['showeraxis']
+    time = data['time']
+    nb_events = len(va)
+
+
+# reconstruct shower directions with plane wave fit
+# cf. ftp://ftp.mi.ingv.it/download/augliera/PER-MARA/article_del_pezzo.pdf
+va_rec = np.zeros((nb_events, 3))
+
+for i in range(nb_events):
+    m = ~np.isnan(time[i])  # mask for triggered stations
+    t = time[i][m]
+    x, y, z = vd[m].T
+
+    # distances between all pairs of stations
+    idx = np.triu_indices(len(t), k=1)
+    dt = np.subtract.outer(t, t)[idx]
+    dx = np.subtract.outer(x, x)[idx]
+    dy = np.subtract.outer(y, y)[idx]
+
+    # all stations on a plane --> neglect z-position
+    dv = np.stack((dx, dy), axis=-1)
+    a1 = np.linalg.inv(np.dot(dv.T, dv))
+
+    # # for using z-position as well
+    # dz = np.subtract.outer(z, z)[idx]
+    # dv = np.stack((dx, dy, dz), axis=-1)
+    # a1 = np.linalg.pinv(np.dot(dv.T, dv))
+
+    a2 = np.einsum('ij,i->j', dv, dt)
+    a3 = np.einsum('ij,j->i', a1, a2)
+
+    # determine z-component using known speed of propagation = c
+    px, py = -a3 * 3E8
+    pz = (1 - px**2 - py**2)**.5
+    va_rec[i] = [px, py, pz]
+
+
+# plot angular separation between true and reconstructed angle
+ang = np.arccos(np.clip(np.sum(va * va_rec, axis=1), -1, 1)) * 180 / np.pi
+r68 = np.percentile(ang, 68)
+
+plt.figure()
+plt.hist(ang, bins=np.linspace(0, 5, 41))
+plt.axvline(r68, c='r')
+plt.xlabel('Angular distance [deg]')
+plt.ylabel('\#')
+plt.grid()
+plt.savefig('fit-planewave.png', bbox_inches='tight')