refactor and tuning

atmosphere.py

@@ -433,12 +433,10 @@ class Atmosphere():
         mask_flat, mask_taylor, mask_numeric = self.__get_method_mask(zenith)
         tmp = np.zeros_like(zenith)
         if np.sum(mask_numeric):
-            print("get vertical height numeric")
             tmp[mask_numeric] = self._get_vertical_height_numeric(*self.__get_arguments(mask_numeric, zenith, X))
         if np.sum(mask_taylor):
             tmp[mask_taylor] = self._get_vertical_height_numeric_taylor(*self.__get_arguments(mask_taylor, zenith, X))
         if np.sum(mask_flat):
-            print("get vertical height flat")
             tmp[mask_flat] = self._get_vertical_height_flat(*self.__get_arguments(mask_flat, zenith, X))
         return tmp
 

plotting.py

@@ -8,6 +8,7 @@ from matplotlib.colors import Normalize
 from utils import distance2showeraxis
 from pylab import *
 
+
 def plot_array(v_stations, values, label='', vmin=None, vmax=None):
     """Plot a map *values* for an detector array specified by *v_stations*.
     """
@@ -29,37 +30,32 @@ def plot_array(v_stations, values, label='', vmin=None, vmax=None):
     ax.set_ylabel('y [km]')
 
 
-def plot_traces_of_array_for_one_event(Smu, Sem, v_axis, v_stations, arraysize = 5):
-    '''plot time traces of tanks around the tank with the smallest distance to shower core
-    Sem = EM-SigmalTrace, Smu = Muon-SignalTrace, arraysize = length of array'''
-    t = np.arange(0, 2001, 25)
-    fig, axes = subplots(arraysize, arraysize, sharex=True, figsize=(29, 16), dpi=80, facecolor='w', edgecolor='k')
+def plot_traces_of_array_for_one_event(Smu, Sem, v_axis, v_stations, n=5):
+    """ Plot time traces of the n^2 central tanks
+        Sem = EM-SigmalTrace, Smu = Muon-SignalTrace
+    """
+    n0 = int(len(v_stations)**.5)
+    i0 = (n0 - n) // 2
+
+    S1 = Smu.reshape(n0, n0, -1)
+    S2 = Sem.reshape(n0, n0, -1)
+
+    fig, axes = subplots(n, n, sharex=True, figsize=(29, 16), facecolor='w')
+    plt.tight_layout()
     t = np.arange(12.5, 2001, 25)
-    tight_layout()
-    coordVcore = np.argmin(distance2showeraxis(v_stations, v_axis))
-    coordX = int(coordVcore/11)
-    coordY = coordVcore%11
-    coords = np.array(0)
-    for j in range(arraysize):
-        for k in range(arraysize):
-            coords = np.append(coords, (coordX-2+j)*11 + (coordY-2+k))
-    coords = coords[1:arraysize*arraysize+1]
-    for i, ax in enumerate(axes.flat):
-        try:
-            h1, h2 = Smu[coords[i]], Sem[coords[i]]
-        except TypeError:
-            h2 = np.histogram(0, bins=t)[0].astype(float)
-        ax.step(t, h1 + h2, c='k', where='mid')
-        ax.step(t, h1, label='$\mu$', where='mid')
-        ax.step(t, h2, label='$e\gamma$', where='mid')
-        ax.legend(title='r=%i' % coords[i], fontsize='x-small')
-        ax.set_xlim(0, 1500)
-        ax.grid(True)
-    for k in range(arraysize):
-        for l in range(arraysize):
-            axes[k, l].set_xlabel('$t$ / ns')
-            axes[k, l].set_ylabel('$S$ / VEM')
-    savefig('trace_VEM.png')#, bbox_inches='tight')
+    for ix in range(n):
+        for iy in range(n):
+            ax = axes[ix, iy]
+            h1 = S1[ix + i0, iy + i0]
+            h2 = S2[ix + i0, iy + i0]
+            ax.step(t, h1 + h2, c='k', where='mid')
+            ax.step(t, h1, label='$\mu$', where='mid')
+            ax.step(t, h2, label='$e\gamma$', where='mid')
+            ax.legend(title='%i, %i' % (ix + i0, iy + i0), fontsize='x-small')
+            ax.set_xlim(0, 1500)
+            ax.grid(True)
+            ax.set_xlabel('$t$ / ns', fontsize='x-small')
+            ax.set_ylabel('$S$ / VEM', fontsize='x-small')
 
 
 if __name__ == '__main__':

shower.py

 from __future__ import division
 
 import numpy as np
-from utils import *
+import utils
 import atmosphere
 import gumbel
-from pylab import *
-atm = atmosphere.Atmosphere()
-import gumbel
 import plotting
+
+
+atm = atmosphere.Atmosphere()
+HEIGHT = 1400
+SPACING = 1500
+
+
 def station_response(r, dX, logE=19, A=1):
     """ Simulate time trace f(t) for a given SD station and shower.
     Args:
@@ -20,12 +24,11 @@ def station_response(r, dX, logE=19, A=1):
     """
     # strength of muon and em signal
     # scaling with energy
-    S0 = 1200 * 10**(logE - 19)
-    # scaling with mass number
+    S0 = 2000 * 10**(logE - 19)
+    # scaling with mass number and relative scaling mu/em
     a = A**0.15
-    b = 2 ## EM scaling
-    S1 = S0 * a / (a + 1) * 1/(b+1)
-    S2 = S0 * 1 / (a + 1) * b/(b+1)
+    S1 = S0 * a / (a + 1) * 1
+    S2 = S0 * 1 / (a + 1) * 3
     # TODO: add scaling with zenith angle (CIC)
     # ...
     # scaling with distance of station to shower axis
@@ -33,11 +36,11 @@ def station_response(r, dX, logE=19, A=1):
     S2 *= np.minimum((r / 1000)**-5.5, 1000)
     # scaling with traversed atmosphere to station
     S1 *= np.minimum((dX / 100)**-0.1, 10)
-    S2 *= np.minimum((dX / 100)**-0.8, 10)
+    S2 *= np.minimum((dX / 100)**-0.4, 10)
 
     # limit total signal, otherwise we get memory problems when drawing that many samples from distribution
     Stot = S1 + S2
-    Smax = 1E6
+    Smax = 1E7
     if Stot > Smax:
         S1 *= Smax / Stot
         S2 *= Smax / Stot
@@ -47,26 +50,26 @@ def station_response(r, dX, logE=19, A=1):
     N2 = np.maximum(int(S2 + S2**.5 * np.random.randn()), 0)  # number of em particles
 
     # parameters of log-normal distributions
-    mean1 = np.log(200 * (r / 750)**1.2 * (1 - 0.2 * dX / 1000.))
-    mean2 = np.log(320 * (r / 750)**1.2 * (1 - 0.1 * dX / 1000.))
+    mean1 = np.log(50 + 140 * (r / 750)**1.4 * (1 - 0.2 * dX / 1000.))
+    mean2 = np.log(80 + 200 * (r / 750)**1.4 * (1 - 0.1 * dX / 1000.))
     sigma1 = 0.7
     sigma2 = 0.7
 
     # draw samples from distributions and create histograms
-    shift = np.exp(mean2) - np.exp(mean1)
+    shift = (np.exp(mean2) - np.exp(mean1)) / 1.5
     bins = np.arange(0, 2001, 25)  # time bins [ns]
     h1 = np.histogram(np.random.lognormal(mean1, sigma1, size=N1), bins=bins)[0]
-    h2 = np.histogram(np.random.lognormal(mean2, sigma2, size=N2) + 0.5*shift, bins=bins)[0]
+    h2 = np.histogram(np.random.lognormal(mean2, sigma2, size=N2) + shift, bins=bins)[0]
 
     # total signal (simplify: 1 particle = 1 VEM)
     return h1, h2
 
 
-def detector_response(logE, mass, v_axis, v_max, v_stations):
+def detector_response(logE, mass, v_axis, v_core, v_max, v_stations):
     """ Simulate the detector response for all SD stations and one event
     """
-    r = distance2showeraxis(v_stations, v_axis)  # radial distance to shower axis
-    phi, zen = vec2ang(v_max - v_stations)  # direction of shower maximum relative to stations
+    r = utils.distance2showeraxis(v_stations, v_core, v_axis)  # radial distance to shower axis
+    phi, zen = utils.vec2ang(v_max - v_stations)  # direction of shower maximum relative to stations
     # distance to shower maximum in [g/cm^2]
     dX = atm.get_atmosphere(
         zen,  # zenith angle of shower maximum seen from each station
@@ -90,30 +93,31 @@ def rand_shower_geometry(N, logE, mass):
     """ Generate random shower geometries: Xmax, axis, core, maximum. """
     N = len(logE)
 
-    # shower axis
+    # 1) shower axis
     phi = 2 * np.pi * (np.random.rand(N) - 0.5)
-    zenith = rand_zenith(N)
-    v_axis = ang2vec(phi, zenith)
-#### TESTIN ######################################################
-    phi = 2 *0* np.pi * (np.random.rand(N) - 0.5)
-    zenith = 60*np.pi/180*rand_zenith(N)
-
-    # shower core on ground (random offset w.r.t. grid origin)
-    v_core = 0*SPACING * (np.random.rand(N, 3) - 0.5)
+    zenith = utils.rand_zenith(N)
+    # phi = 0 * np.ones(N)
+    # zenith = np.deg2rad(60) * np.ones(N)
+    v_axis = utils.ang2vec(phi, zenith)
+
+    # 2) shower core on ground (random offset w.r.t. grid origin)
+    v_core = SPACING * (np.random.rand(N, 3) - 0.5)
     v_core[:, 2] = HEIGHT  # core is always at ground level
 
-    # shower maximum, require h > hmin
+    # 3) shower maximum, require h > hmin
     Xmax = np.empty(N)
     hmax = np.empty(N)
     i = 0
     while i < N:
         Xmax[i] = gumbel.rand_gumbel(logE[i], mass[i])
-        hmax[i] = atm.get_vertical_height(zenith[i], Xmax[i])
+        try:
+            hmax[i] = atm.get_vertical_height(zenith[i], Xmax[i])
+        except:
+            continue
         if hmax[i] < HEIGHT + 200:
             continue  # shower maximum too low, repeat
         i += 1
-
-    # position of shower maximum
+        print(i)
     dmax = (hmax - HEIGHT) / np.cos(zenith)  # distance to Xmax
     v_max = v_core + v_axis * dmax[:, np.newaxis]
 
@@ -122,26 +126,29 @@ def rand_shower_geometry(N, logE, mass):
 
 if __name__ == '__main__':
     # detector array
-    nb_stations = 11**2  # number of stations
-    v_stations = triangular_array(nb_stations**.5)  # x,y,z coordinates of SD stations
+    n = 11
+    v_stations = utils.station_coordinates(n, layout='cartesian')  # x,y,z coordinates of SD stations
+    nb_stations = n**2  # number of stations
 
     # showers
-    nb_events = 2
+    print('simulating showers')
+    nb_events = 100
     logE = 18.5 + 1.5 * np.random.rand(nb_events)
-    logE = 20 *np.ones(nb_events)
+    # logE = 20 * np.ones(nb_events)
     mass = 1 * np.ones(nb_events)
     Xmax, v_axis, v_core, v_max = rand_shower_geometry(nb_events, logE, mass)
 
     # detector response for each event
+    print('simulating detector response')
     T = np.zeros((nb_events, nb_stations))
     S1 = np.zeros((nb_events, nb_stations, 80))
     S2 = np.zeros((nb_events, nb_stations, 80))
     for i in range(nb_events):
         print(i)
-        s1, s2 = detector_response(logE[i], mass[i], v_axis[i], v_max[i], v_stations)
+        s1, s2 = detector_response(logE[i], mass[i], v_axis[i], v_core[i], v_max[i], v_stations)
         S1[i] = s1
         S2[i] = s2
-        T[i] = arrival_time_planar(v_stations, v_core[i], v_axis[i])
+        T[i] = utils.arrival_time_planar(v_stations, v_core[i], v_axis[i])
 
     # add per ton noise on arrival time
     T += 20E-9 * np.random.randn(*T.shape)
@@ -156,7 +163,6 @@ if __name__ == '__main__':
 #    S1 += 0.5 + 0.2 * np.random.randn(*S1.shape)
 #    S2 += 0.5 + 0.2 * np.random.randn(*S2.shape)
 
-    plotting.plot_traces_of_array_for_one_event(Smu=S1[0], Sem=S2[0], v_axis=v_axis[0], v_stations=v_stations, arraysize = 5)
     # trigger threshold: use only stations with sufficient signal-to-noise
 ##    c = S.sum(axis=-1) < 80 * 0.55
 ##    T[c] = np.NaN
@@ -176,3 +182,12 @@ if __name__ == '__main__':
 #        showeraxis=v_axis,
 #        showermax=v_max,
 #        detector=v_stations)
+
+    plotting.plot_traces_of_array_for_one_event(
+        Smu=S1[0], Sem=S2[0], v_axis=v_axis[0], v_stations=v_stations, n=5)
+
+    # r = utils.distance2showeraxis(v_stations, v_core[0], v_axis[0])
+    # plotting.plot_array(v_stations, r)
+
+    import matplotlib.pyplot as plt
+    plt.show()

stuff/test_lognormal_param.py

 from pylab import *
 
-
 close('all')
-# r = linspace(100, 3000)
-# figure()
-# plot(r, 220 * (r / 750)**1.2)
-# plot(r, 420 * (r / 750)**1.2)
-# plot([750, 850, 1250], [220, 260, 500], 'C0+')
-# plot([750, 850, 1250], [420, 500, 1000], 'C1+')
-# xlabel('$r$ [m]')
-# ylabel('$\Delta t$ [ns]')
-# grid()
-
-
-# for r in [750, 850, 1250]:
-#     mean1 = np.log(200 * (r / 750)**1.2)  # * (1 - 0.2 * dX / 1000.)
-#     mean2 = np.log(320 * (r / 750)**1.2)  # * (1 - 0.1 * dX / 1000.)
+
+
+def mean1(r, dx):
+    return 30 + 180 * (r / 750)**1.4 * (1 - 0.2 * dX / 1000.)
+
+def mean2(r, dx):
+    return 50 + 300 * (r / 750)**1.4 * (1 - 0.1 * dX / 1000.)
+
+
+figure()
+r = linspace(100, 3000)
+dX = 600
+plot(r, mean1(r, dX))
+plot(r, mean2(r, dX))
+plot([750, 850, 1250], [220, 260, 500], 'C0+')
+plot([750, 850, 1250], [420, 500, 1000], 'C1+')
+xlabel('$r$ [m]')
+ylabel('$\Delta t$ [ns]')
+grid()
+
+
+for r in [750, 850, 1250]:
+    m1 = np.log(mean1(r, dX))
+    m2 = np.log(mean2(r, dX))
+
+    figure()
+    bins = np.arange(0, 2001, 25)  # time bins [ns]
+    axvline(exp(m1), color='C0')
+    axvline(exp(m2), color='C1')
+    hist(np.random.lognormal(m1, 0.7, size=10000), bins=bins, alpha=0.7)
+    hist(np.random.lognormal(m2, 0.7, size=10000) + exp(m2) - exp(m1), bins=bins, alpha=0.7)
+    xlim(0, 1500)
+    show()
+
+
+# r = 750
+# for dX in [300, 600, 900]:
+#     mean1 = np.log(200 * (r / 750)**1.2 * (1 - 0.2 * dX / 1000.))
+#     mean2 = np.log(320 * (r / 750)**1.2 * (1 - 0.1 * dX / 1000.))
 #     sigma1 = 0.7  # * (r / 1000.)**0.15 * (1 - 0.2 * dX / 1000.)
 #     sigma2 = 0.7  # * (r / 1000.)**0.15
 
@@ -23,25 +47,8 @@ close('all')
 #     bins = np.arange(0, 2001, 25)  # time bins [ns]
 #     axvline(exp(mean1), color='C0')
 #     axvline(exp(mean2), color='C1')
+#     shift = exp(mean2)-exp(mean1)
 #     hist(np.random.lognormal(mean1, sigma1, size=10000), bins=bins, alpha=0.7)
-#     hist(np.random.lognormal(mean2, sigma2, size=10000) + exp(mean2)-exp(mean1), bins=bins, alpha=0.7)
+#     hist(np.random.lognormal(mean2, sigma2, size=10000) + shift, bins=bins, alpha=0.7)
 #     xlim(0, 1500)
 #     show()
-
-
-r = 750
-for dX in [300, 600, 900]:
-    mean1 = np.log(200 * (r / 750)**1.2 * (1 - 0.2 * dX / 1000.))
-    mean2 = np.log(320 * (r / 750)**1.2 * (1 - 0.1 * dX / 1000.))
-    sigma1 = 0.7  # * (r / 1000.)**0.15 * (1 - 0.2 * dX / 1000.)
-    sigma2 = 0.7  # * (r / 1000.)**0.15
-
-    figure()
-    bins = np.arange(0, 2001, 25)  # time bins [ns]
-    axvline(exp(mean1), color='C0')
-    axvline(exp(mean2), color='C1')
-    shift = exp(mean2)-exp(mean1)
-    hist(np.random.lognormal(mean1, sigma1, size=10000), bins=bins, alpha=0.7)
-    hist(np.random.lognormal(mean2, sigma2, size=10000) + shift, bins=bins, alpha=0.7)
-    xlim(0, 1500)
-    show()

utils.py

@@ -3,34 +3,38 @@ from __future__ import division
 import numpy as np
 
 
-HEIGHT = 1400  # detector height in [m]
-SPACING = 1500  # detector spacing in [m]
-
-
 def rectangular_array(n=11):
-    """ Coordinates for rectangular array with n^2 stations and given spacing.
-        Returns: n^2 x 3 array of x,y,z coordinates for each station.
-    """
+    """ Return x,y coordinates for rectangular array with n^2 stations. """
     n0 = (n - 1) / 2
-    x, y = (np.mgrid[0:n, 0:n].astype(float) - n0) * SPACING
-    z = np.ones_like(x) * HEIGHT  # z-position
-    return np.dstack([x, y, z]).reshape(-1, 3)
+    return (np.mgrid[0:n, 0:n].astype(float) - n0)
 
 
 def triangular_array(n=11, offset=True):
-    """ Coordinates for triangular array with n^2 stations and given spacing.
-        Returns: n^2 x 3 array of x,y,z coordinates for each station.
-    """
+    """ Return x,y coordinates for triangular array with n^2 stations. """
     n0 = (n - 1) / 2
     x, y = np.mgrid[0:n, 0:n].astype(float) - n0
     if offset:  # offset coordinates
         x += 0.5 * (y % 2)
     else:  # axial coordinates
         x += 0.5 * y
-    x *= SPACING
-    y *= np.sin(np.pi / 3) * SPACING
-    z = np.ones_like(x) * HEIGHT  # z-position
-    return np.dstack([x, y, z]).reshape(-1, 3)
+    y *= np.sin(np.pi / 3)
+    return x, y
+
+
+def station_coordinates(n=11, layout='axial', spacing=1500, height=1400):
+    """ Return array of n^2*(x,y,z) coordinates of SD stations for given layout. """
+    if layout == 'axial':
+        x, y = triangular_array(n, offset=False)
+    elif layout == 'offset':
+        x, y = triangular_array(n, offset=True)
+    elif layout == 'cartesian':
+        x, y = rectangular_array(n)
+    else:
+        raise ValueError('layout must be one of axial, offset, cartesian')
+    x = x.reshape(n**2) * spacing
+    y = y.reshape(n**2) * spacing
+    z = np.ones_like(x) * height
+    return np.c_[x, y, z]
 
 
 def rand_zenith(N=1, zmax=np.pi / 3):
@@ -71,23 +75,25 @@ def vec2ang(v):
     return phi, zenith
 
 
-def distance2showerplane(v, va):
+def distance2showerplane(v, vc, va):
     """ Get shortest distance to shower plane
     Args:
         v (Nx3 array): array of positions
         va (3 array): shower axis = normal vector of the shower plane
+        vc (3 array): shower core
     """
-    return np.dot(v, va)
+    return np.dot(v - vc, va)
 
 
-def distance2showeraxis(v, va):
+def distance2showeraxis(v, vc, va):
     """ Shortest distance to shower axis.
     Args:
         v (Nx3 array): array of positions
         va (3 array): shower axis
+        vc (3 array): shower core
     """
-    d = distance2showerplane(v, va)
-    vp = v - np.outer(d, va)
+    d = distance2showerplane(v, vc, va)
+    vp = v - vc - np.outer(d, va)
     return np.linalg.norm(vp, axis=-1)
 
 
@@ -110,5 +116,5 @@ def arrival_time_planar(v, vc, va):
     Return:
         array: arrival times [s]
     """
-    d = distance2showerplane(v - vc, -va)
+    d = distance2showerplane(v, vc, -va)
     return d / 3E8