shower.py 7.59 KB
Newer Older
DavidWalz's avatar
DavidWalz committed
1
2
3
from __future__ import division

import numpy as np
DavidWalz's avatar
DavidWalz committed
4
import utils
DavidWalz's avatar
DavidWalz committed
5
6
import atmosphere
import gumbel
JGlombitza's avatar
JGlombitza committed
7
import plotting
DavidWalz's avatar
DavidWalz committed
8
9
10
11
12
13


HEIGHT = 1400
SPACING = 1500


DavidWalz's avatar
tune    
DavidWalz committed
14
def station_response(r, dX, logE=19, zenith=0, mass=1):
DavidWalz's avatar
DavidWalz committed
15
16
17
18
19
    """ Simulate time trace f(t) for a given SD station and shower.
    Args:
        r (float): radial distance to shower axis in [m]
        dX (float): distance to Xmax in [g/cm^2]
        logE (float): log10(E/eV)
DavidWalz's avatar
tune    
DavidWalz committed
20
21
        zenith (float): zenith angle
        mass (float): mass number
DavidWalz's avatar
DavidWalz committed
22
23
24
    Returns:
        time traces f_muon(t) and f_em(t) in [VEM] for t = 0 - 2000 ns in bins of 25 ns
    """
DavidWalz's avatar
DavidWalz committed
25
    # signal strength, scaling with energy
JGlombitza's avatar
JGlombitza committed
26
    S0 = 900 * 10**(logE - 19)
DavidWalz's avatar
tune    
DavidWalz committed
27
28
29
    # # scaling with zenith angle (similar to S1000 --> S38 relation, CIC)
    # x = np.cos(zenith)**2 - np.cos(np.deg2rad(38))**2
    # S0 *= 1 + 0.95 * x - 1.4 * x**2 - 1.1 * x**3
30

DavidWalz's avatar
tune    
DavidWalz committed
31
    # relative scaling mu/em and scaling with mass number
32
33
34
35
36
    r1 = 0.3 * mass**0.15
    r2 = 0.7
    S1 = S0 * r1 / (r1 + r2)
    S2 = S0 * r2 / (r1 + r2)

DavidWalz's avatar
DavidWalz committed
37
    # scaling with distance of station to shower axis
DavidWalz's avatar
tune    
DavidWalz committed
38
39
    S1 *= (np.maximum(r, 150) / 1000)**-4.7
    S2 *= (np.maximum(r, 150) / 1000)**-6.1
40

DavidWalz's avatar
DavidWalz committed
41
    # scaling with traversed atmosphere to station
DavidWalz's avatar
DavidWalz committed
42
43
    S1 *= np.minimum((dX / 100)**-0.1, 10)
    S2 *= np.minimum((dX / 100)**-0.4, 10)
DavidWalz's avatar
DavidWalz committed
44

45
    # limit total signal (saturation / memory constraints)
DavidWalz's avatar
DavidWalz committed
46
    Stot = S1 + S2
DavidWalz's avatar
DavidWalz committed
47
    Smax = 1E7
DavidWalz's avatar
DavidWalz committed
48
49
50
51
52
53
54
55
56
    if Stot > Smax:
        S1 *= Smax / Stot
        S2 *= Smax / Stot

    # draw number of detected muons and em particles
    N1 = np.maximum(int(S1 + S1**.5 * np.random.randn()), 0)  # number of muons
    N2 = np.maximum(int(S2 + S2**.5 * np.random.randn()), 0)  # number of em particles

    # parameters of log-normal distributions
DavidWalz's avatar
DavidWalz committed
57
58
    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.))
DavidWalz's avatar
DavidWalz committed
59
60
61
62
    sigma1 = 0.7
    sigma2 = 0.7

    # draw samples from distributions and create histograms
DavidWalz's avatar
DavidWalz committed
63
    shift = (np.exp(mean2) - np.exp(mean1)) / 1.5
DavidWalz's avatar
DavidWalz committed
64
65
    bins = np.arange(0, 2001, 25)  # time bins [ns]
    h1 = np.histogram(np.random.lognormal(mean1, sigma1, size=N1), bins=bins)[0]
DavidWalz's avatar
DavidWalz committed
66
    h2 = np.histogram(np.random.lognormal(mean2, sigma2, size=N2) + shift, bins=bins)[0]
DavidWalz's avatar
DavidWalz committed
67
68
69
70
71

    # total signal (simplify: 1 particle = 1 VEM)
    return h1, h2


DavidWalz's avatar
DavidWalz committed
72
def detector_response(logE, mass, v_axis, v_core, v_max, v_stations):
DavidWalz's avatar
DavidWalz committed
73
    """ Simulate the detector response for all SD stations and one event. """
DavidWalz's avatar
tune    
DavidWalz committed
74
75
    _, zenith = utils.vec2ang(v_axis)

DavidWalz's avatar
DavidWalz committed
76
77
    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
DavidWalz's avatar
DavidWalz committed
78
79
80

    # distance from each station to shower maximum in [g/cm^2]
    dX = atmosphere.slant_depth(v_stations[:, 2], zen) - atmosphere.slant_depth(v_max[2], zen)
DavidWalz's avatar
DavidWalz committed
81
82
83
84
85
86
87

    # time traces for each station
    n = len(v_stations)
    S1 = np.zeros((n, 80))  # muon traces
    S2 = np.zeros((n, 80))  # em traces

    for j in range(n):
DavidWalz's avatar
tune    
DavidWalz committed
88
        h1, h2 = station_response(r=r[j], dX=dX[j], logE=logE, zenith=zenith, mass=mass)
DavidWalz's avatar
DavidWalz committed
89
90
91
92
93
94
        S1[j] = h1
        S2[j] = h2

    return S1, S2


DavidWalz's avatar
DavidWalz committed
95
96
97
def rand_shower_geometry(logE, mass):
    """ Generate random shower geometries: Xmax, axis, core, maximum, virtual origin """
    nb_showers = len(logE)
DavidWalz's avatar
DavidWalz committed
98

99
    # 1) random shower axis
DavidWalz's avatar
DavidWalz committed
100
101
    phi = 2 * np.pi * (np.random.rand(nb_showers) - 0.5)
    zenith = utils.rand_zenith(nb_showers)
DavidWalz's avatar
DavidWalz committed
102
103
    v_axis = utils.ang2vec(phi, zenith)

104
    # 2) random shower core on ground (offset w.r.t. grid origin)
105
106
107
108
109
110
    r = SPACING / 2 * np.random.rand(nb_showers)**.5
    p = np.random.rand(nb_showers) * 2 * np.pi
    x = r * np.cos(p)
    y = r * np.sin(p)
    z = HEIGHT * np.ones_like(p)
    v_core = np.c_[x, y, z]
DavidWalz's avatar
DavidWalz committed
111

DavidWalz's avatar
DavidWalz committed
112
113
114
    # 3) random shower maximum, require Xmax 200m above ground
    Xmax = gumbel.rand_gumbel(logE, mass)
    Xmax_max = atmosphere.slant_depth(HEIGHT + 200, zenith)
DavidWalz's avatar
DavidWalz committed
115

DavidWalz's avatar
DavidWalz committed
116
    while not np.all(Xmax < Xmax_max):
117
        idx = Xmax > Xmax_max  # resample shower maxima less than 200 above ground
DavidWalz's avatar
DavidWalz committed
118
        Xmax[idx] = gumbel.rand_gumbel(logE[idx], mass[idx])
DavidWalz's avatar
DavidWalz committed
119

DavidWalz's avatar
DavidWalz committed
120
121
122
123
    # 4) point of shower maximum
    h = atmosphere.height_at_slant_depth(Xmax, zenith)
    d = (h - HEIGHT) / np.cos(zenith)
    v_max = v_core + v_axis * d[:, np.newaxis]
DavidWalz's avatar
DavidWalz committed
124

125
126
    # 5) virtual origin at 50% of Xmax
    h = atmosphere.height_at_slant_depth(0.5 * Xmax, zenith)
DavidWalz's avatar
DavidWalz committed
127
128
129
130
131
132
    d = (h - HEIGHT) / np.cos(zenith)
    v_origin = v_core + v_axis * d[:, np.newaxis]

    return Xmax, v_axis, v_core, v_max, v_origin


133
def rand_events(logE, mass, v_stations, fname=None, wavefront='planar'):
DavidWalz's avatar
DavidWalz committed
134
    """Simulate events for given energy, mass, and detector positions."""
DavidWalz's avatar
DavidWalz committed
135

DavidWalz's avatar
DavidWalz committed
136
137
138
139
    nb_events = len(logE)
    nb_stations = len(v_stations)

    # simulation showers
DavidWalz's avatar
DavidWalz committed
140
    print('simulating showers')
DavidWalz's avatar
DavidWalz committed
141
    Xmax, v_axis, v_core, v_max, v_origin = rand_shower_geometry(logE, mass)
DavidWalz's avatar
DavidWalz committed
142
143

    # detector response for each event
DavidWalz's avatar
DavidWalz committed
144
    print('simulating detector response')
DavidWalz's avatar
DavidWalz committed
145
    T = np.zeros((nb_events, nb_stations))
JGlombitza's avatar
JGlombitza committed
146
147
    S1 = np.zeros((nb_events, nb_stations, 80))
    S2 = np.zeros((nb_events, nb_stations, 80))
DavidWalz's avatar
DavidWalz committed
148
    for i in range(nb_events):
149
        # print('%4i, logE = %.2f' % (i, logE[i]))
DavidWalz's avatar
DavidWalz committed
150
        S1[i], S2[i] = detector_response(logE[i], mass[i], v_axis[i], v_core[i], v_max[i], v_stations)
151
152
153
154
155
156
        if wavefront == 'planar':
            T[i] = utils.arrival_time_planar(v_stations, v_core[i], v_axis[i])
        elif wavefront == 'spherical':
            T[i] = utils.arrival_time_spherical(v_stations, v_origin[i], v_core[i])
        else:
            raise ValueError('wavefront must be "planar" or "spherical"')
DavidWalz's avatar
DavidWalz committed
157

DavidWalz's avatar
DavidWalz committed
158
159
    # total signal
    S = S1 + S2
DavidWalz's avatar
DavidWalz committed
160
    # add per ton noise on arrival time
161
    Stot = S.sum(axis=-1)
162
163
164
    sigma = 60E-9 * 8. / (1 + np.log10(Stot + 1))  # varies from ~(1-8)*60 ns
    T += sigma * np.random.randn(*T.shape)
    # add time offset per event (to account for random core position)
DavidWalz's avatar
DavidWalz committed
165
    T += 100E-9 * (np.random.rand(nb_events, 1) - 0.5)
DavidWalz's avatar
DavidWalz committed
166
    # add relative noise to signal pulse
167
    S += 0.05 * S * np.random.randn(*S.shape)
DavidWalz's avatar
DavidWalz committed
168
    # add absolute noise to signal pulse
DavidWalz's avatar
tune    
DavidWalz committed
169
170
    noise_level = 1.2
    S += noise_level * (1 + 0.5 * np.random.randn(*S.shape))
171
    S = np.clip(S, 0, None)
DavidWalz's avatar
DavidWalz committed
172
    # trigger threshold: use only stations with sufficient signal-to-noise
DavidWalz's avatar
tune    
DavidWalz committed
173
    c = S.sum(axis=-1) < 80 * noise_level * 1.2
DavidWalz's avatar
DavidWalz committed
174
175
    T[c] = np.NaN
    S[c] = np.NaN
DavidWalz's avatar
DavidWalz committed
176

DavidWalz's avatar
DavidWalz committed
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
    return {
        'logE': logE,
        'mass': mass,
        'Xmax': Xmax,
        'time': T,
        'signal': S,
        'signal1': S1,
        'signal2': S2,
        'showercore': v_core,
        'showeraxis': v_axis,
        'showermax': v_max,
        'detector': v_stations}


if __name__ == '__main__':
    # detector array, vector of (x,y,z) positions
    v_stations = utils.station_coordinates(11, layout='offset')
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210

    # simulate events
    n = 1000
    logE = 18.5 + 1.5 * np.random.rand(n)
    mass = 1
    data = rand_events(logE, mass, v_stations)

    phi, zenith = utils.vec2ang(data['showeraxis'])
    plotting.plot_time_distribution(data['time'], fname='time_distribution.png')
    plotting.plot_signal_distribution(data['signal'], fname='signal_distribution.png')
    plotting.plot_energy_distribution(data['logE'], fname='energy_distribution.png')
    plotting.plot_xmax_distribution(data['Xmax'], fname='xmax_distribution.png')
    plotting.plot_zenith_distribution(zenith, fname='zenith_distribution.png')
    plotting.plot_phi_distribution(phi, fname='phi_distribution.png')
    plotting.plot_stations_vs_energy(data['logE'], data['signal'], fname='stations_vs_energy.png')
    plotting.plot_stations_vs_zenith(zenith, data['signal'], fname='stations_vs_zenith.png')
    plotting.plot_array_traces(Smu=data['signal1'][1], Sem=data['signal2'][1], v_stations=data['detector'], n=5, fname='example-trace.png')