shower.py 7.59 KB
Newer Older
 DavidWalz committed May 04, 2017 1 2 3 ``````from __future__ import division import numpy as np `````` DavidWalz committed May 04, 2017 4 ``````import utils `````` DavidWalz committed May 04, 2017 5 6 ``````import atmosphere import gumbel `````` JGlombitza committed May 04, 2017 7 ``````import plotting `````` DavidWalz committed May 04, 2017 8 9 10 11 12 13 `````` HEIGHT = 1400 SPACING = 1500 `````` DavidWalz committed May 05, 2017 14 ``````def station_response(r, dX, logE=19, zenith=0, mass=1): `````` DavidWalz committed May 04, 2017 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 committed May 05, 2017 20 21 `````` zenith (float): zenith angle mass (float): mass number `````` DavidWalz committed May 04, 2017 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 committed May 16, 2017 25 `````` # signal strength, scaling with energy `````` JGlombitza committed May 09, 2017 26 `````` S0 = 900 * 10**(logE - 19) `````` DavidWalz committed May 05, 2017 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 committed May 05, 2017 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 committed May 04, 2017 37 `````` # scaling with distance of station to shower axis `````` DavidWalz committed May 05, 2017 38 39 `````` S1 *= (np.maximum(r, 150) / 1000)**-4.7 S2 *= (np.maximum(r, 150) / 1000)**-6.1 `````` 40 `````` `````` DavidWalz committed May 04, 2017 41 `````` # scaling with traversed atmosphere to station `````` DavidWalz committed May 16, 2017 42 43 `````` S1 *= np.minimum((dX / 100)**-0.1, 10) S2 *= np.minimum((dX / 100)**-0.4, 10) `````` DavidWalz committed May 04, 2017 44 `````` `````` 45 `````` # limit total signal (saturation / memory constraints) `````` DavidWalz committed May 04, 2017 46 `````` Stot = S1 + S2 `````` DavidWalz committed May 04, 2017 47 `````` Smax = 1E7 `````` DavidWalz committed May 04, 2017 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 committed May 04, 2017 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 committed May 04, 2017 59 60 61 62 `````` sigma1 = 0.7 sigma2 = 0.7 # draw samples from distributions and create histograms `````` DavidWalz committed May 04, 2017 63 `````` shift = (np.exp(mean2) - np.exp(mean1)) / 1.5 `````` DavidWalz committed May 04, 2017 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 committed May 04, 2017 66 `````` h2 = np.histogram(np.random.lognormal(mean2, sigma2, size=N2) + shift, bins=bins)[0] `````` DavidWalz committed May 04, 2017 67 68 69 70 71 `````` # total signal (simplify: 1 particle = 1 VEM) return h1, h2 `````` DavidWalz committed May 04, 2017 72 ``````def detector_response(logE, mass, v_axis, v_core, v_max, v_stations): `````` DavidWalz committed May 16, 2017 73 `````` """ Simulate the detector response for all SD stations and one event. """ `````` DavidWalz committed May 05, 2017 74 75 `````` _, zenith = utils.vec2ang(v_axis) `````` DavidWalz committed May 04, 2017 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 committed May 05, 2017 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 committed May 04, 2017 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 committed May 05, 2017 88 `````` h1, h2 = station_response(r=r[j], dX=dX[j], logE=logE, zenith=zenith, mass=mass) `````` DavidWalz committed May 04, 2017 89 90 91 92 93 94 `````` S1[j] = h1 S2[j] = h2 return S1, S2 `````` DavidWalz committed May 16, 2017 95 96 97 ``````def rand_shower_geometry(logE, mass): """ Generate random shower geometries: Xmax, axis, core, maximum, virtual origin """ nb_showers = len(logE) `````` DavidWalz committed May 04, 2017 98 `````` `````` DavidWalz committed May 05, 2017 99 `````` # 1) random shower axis `````` DavidWalz committed May 16, 2017 100 101 `````` phi = 2 * np.pi * (np.random.rand(nb_showers) - 0.5) zenith = utils.rand_zenith(nb_showers) `````` DavidWalz committed May 04, 2017 102 103 `````` v_axis = utils.ang2vec(phi, zenith) `````` DavidWalz committed May 05, 2017 104 `````` # 2) random shower core on ground (offset w.r.t. grid origin) `````` DavidWalz committed May 19, 2017 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 committed May 04, 2017 111 `````` `````` DavidWalz committed May 16, 2017 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 committed May 04, 2017 115 `````` `````` DavidWalz committed May 16, 2017 116 `````` while not np.all(Xmax < Xmax_max): `````` DavidWalz committed May 19, 2017 117 `````` idx = Xmax > Xmax_max # resample shower maxima less than 200 above ground `````` DavidWalz committed May 16, 2017 118 `````` Xmax[idx] = gumbel.rand_gumbel(logE[idx], mass[idx]) `````` DavidWalz committed May 04, 2017 119 `````` `````` DavidWalz committed May 16, 2017 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 committed May 04, 2017 124 `````` `````` 125 126 `````` # 5) virtual origin at 50% of Xmax h = atmosphere.height_at_slant_depth(0.5 * Xmax, zenith) `````` DavidWalz committed May 16, 2017 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 committed May 16, 2017 134 `````` """Simulate events for given energy, mass, and detector positions.""" `````` DavidWalz committed May 04, 2017 135 `````` `````` DavidWalz committed May 16, 2017 136 137 138 139 `````` nb_events = len(logE) nb_stations = len(v_stations) # simulation showers `````` DavidWalz committed May 04, 2017 140 `````` print('simulating showers') `````` DavidWalz committed May 16, 2017 141 `````` Xmax, v_axis, v_core, v_max, v_origin = rand_shower_geometry(logE, mass) `````` DavidWalz committed May 04, 2017 142 143 `````` # detector response for each event `````` DavidWalz committed May 04, 2017 144 `````` print('simulating detector response') `````` DavidWalz committed May 04, 2017 145 `````` T = np.zeros((nb_events, nb_stations)) `````` JGlombitza committed May 04, 2017 146 147 `````` S1 = np.zeros((nb_events, nb_stations, 80)) S2 = np.zeros((nb_events, nb_stations, 80)) `````` DavidWalz committed May 04, 2017 148 `````` for i in range(nb_events): `````` 149 `````` # print('%4i, logE = %.2f' % (i, logE[i])) `````` DavidWalz committed May 16, 2017 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 committed May 04, 2017 157 `````` `````` DavidWalz committed May 05, 2017 158 159 `````` # total signal S = S1 + S2 `````` DavidWalz committed May 04, 2017 160 `````` # add per ton noise on arrival time `````` DavidWalz committed May 18, 2017 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 committed May 04, 2017 165 `````` T += 100E-9 * (np.random.rand(nb_events, 1) - 0.5) `````` DavidWalz committed May 05, 2017 166 `````` # add relative noise to signal pulse `````` 167 `````` S += 0.05 * S * np.random.randn(*S.shape) `````` DavidWalz committed May 05, 2017 168 `````` # add absolute noise to signal pulse `````` DavidWalz committed May 05, 2017 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 committed May 04, 2017 172 `````` # trigger threshold: use only stations with sufficient signal-to-noise `````` DavidWalz committed May 05, 2017 173 `````` c = S.sum(axis=-1) < 80 * noise_level * 1.2 `````` DavidWalz committed May 05, 2017 174 175 `````` T[c] = np.NaN S[c] = np.NaN `````` DavidWalz committed May 04, 2017 176 `````` `````` DavidWalz committed May 16, 2017 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')``````