import numpy as np import matplotlib.pyplot as plt """ Script originally written by Ben Hokanson-Fasig """ ########################## Generate the data to plot ########################## ant1 = (0, 0, -2) ant2 = (0, 0, -60) nu_vertex = (-175, 0, -100) nu_direction = np.array([2.0, 0.0, -1.0]) nu_direction /= np.linalg.norm(nu_direction) theta_c = 0.926603721063296 x1 = np.array( [-1.75000000e+02, -1.73703777e+02, -1.72404337e+02, -1.71101614e+02, -1.69795544e+02, -1.68486059e+02, -1.67173089e+02, -1.65856562e+02, -1.64536406e+02, -1.63212546e+02, -1.61884903e+02, -1.60553399e+02, -1.59217952e+02, -1.57878477e+02, -1.56534889e+02, -1.55187099e+02, -1.53835014e+02, -1.52478542e+02, -1.51117584e+02, -1.49752041e+02, -1.48381811e+02, -1.47006786e+02, -1.45626858e+02, -1.44241913e+02, -1.42851836e+02, -1.41456507e+02, -1.40055800e+02, -1.38649589e+02, -1.37237741e+02, -1.35820118e+02, -1.34396581e+02, -1.32966981e+02, -1.31531169e+02, -1.30088987e+02, -1.28640273e+02, -1.27184858e+02, -1.25722569e+02, -1.24253224e+02, -1.22776636e+02, -1.21292610e+02, -1.19800943e+02, -1.18301425e+02, -1.16793837e+02, -1.15277951e+02, -1.13753531e+02, -1.12220329e+02, -1.10678087e+02, -1.09126538e+02, -1.07565401e+02, -1.05994382e+02, -1.04413177e+02, -1.02821463e+02, -1.01218906e+02, -9.96051552e+01, -9.79798407e+01, -9.63425762e+01, -9.46929550e+01, -9.30305497e+01, -9.13549103e+01, -8.96655625e+01, -8.79620062e+01, -8.62437129e+01, -8.45101241e+01, -8.27606484e+01, -8.09946590e+01, -7.92114905e+01, -7.74104358e+01, -7.55907420e+01, -7.37516065e+01, -7.18921721e+01, -7.00115218e+01, -6.81086731e+01, -6.61825706e+01, -6.42320789e+01, -6.22559738e+01, -6.02529319e+01, -5.82215196e+01, -5.61601795e+01, -5.40672152e+01, -5.19407732e+01, -4.97788218e+01, -4.75791266e+01, -4.53392209e+01, -4.30563706e+01, -4.07275318e+01, -3.83492996e+01, -3.59178458e+01, -3.34288407e+01, -3.08773572e+01, -2.82577485e+01, -2.55634925e+01, -2.27869891e+01, -1.99192948e+01, -1.69497642e+01, -1.38655612e+01, -1.06509714e+01, -7.28641062e+00, -3.74694605e+00, 7.95807864e-13] ) z1 = np.array( [-100., -99., -98., -97., -96., -95., -94., -93., -92., -91., -90., -89., -88., -87., -86., -85., -84., -83., -82., -81., -80., -79., -78., -77., -76., -75., -74., -73., -72., -71., -70., -69., -68., -67., -66., -65., -64., -63., -62., -61., -60., -59., -58., -57., -56., -55., -54., -53., -52., -51., -50., -49., -48., -47., -46., -45., -44., -43., -42., -41., -40., -39., -38., -37., -36., -35., -34., -33., -32., -31., -30., -29., -28., -27., -26., -25., -24., -23., -22., -21., -20., -19., -18., -17., -16., -15., -14., -13., -12., -11., -10., -9., -8., -7., -6., -5., -4., -3., -2.] ) x2 = np.array( [-1.75000000e+02, -1.72030338e+02, -1.69033398e+02, -1.66008064e+02, -1.62953146e+02, -1.59867380e+02, -1.56749413e+02, -1.53597801e+02, -1.50410997e+02, -1.47187340e+02, -1.43925046e+02, -1.40622192e+02, -1.37276701e+02, -1.33886328e+02, -1.30448633e+02, -1.26960964e+02, -1.23420427e+02, -1.19823853e+02, -1.16167766e+02, -1.12448330e+02, -1.08661305e+02, -1.04801979e+02, -1.00865096e+02, -9.68447564e+01, -9.27343117e+01, -8.85262176e+01, -8.42118592e+01, -7.97813255e+01, -7.52231208e+01, -7.05237863e+01, -6.56673985e+01, -6.06348894e+01, -5.54031047e+01, -4.99434665e+01, -4.42200142e+01, -3.81864309e+01, -3.17813280e+01, -2.49203503e+01, -1.74820122e+01, -9.27983644e+00, -1.16955334e-10] ) z2 = np.array( [-100., -99., -98., -97., -96., -95., -94., -93., -92., -91., -90., -89., -88., -87., -86., -85., -84., -83., -82., -81., -80., -79., -78., -77., -76., -75., -74., -73., -72., -71., -70., -69., -68., -67., -66., -65., -64., -63., -62., -61., -60.] ) ################ Define a helper function for rotations in x/z ################ def rot(theta): """Calculates rotation matrix for a given angle in the x-z plane.""" return np.array(((np.cos(theta), 0, -np.sin(theta)), (0, 1, 0 ), (np.sin(theta), 0, np.cos(theta)))) ############################ Draw the event diagram ########################### # Use the x-axis and y-axis ranges to determine the aspect ratio of the image xlims = (-210, 20) ylims = (-130, 10) ydim = 5 xdim = round(ydim*(np.diff(xlims)/np.diff(ylims))[0], 1) print("figsize:", (xdim, ydim)) # Set the figure size to use our calculated dimensions plt.figure(figsize=(xdim, ydim)) # Set a horizontal line at zero to represent the ice surface plt.axhline(0, c='k') # Add a filled circle at the neutrino interaction vertex position plt.scatter(nu_vertex[0], nu_vertex[2], marker='o', s=75, color='C1', zorder=10, label="Shower Vertex") # Draw an arrow pointing in the direction of the particle shower arrow_scale = 10 plt.arrow(nu_vertex[0], nu_vertex[2], nu_direction[0]*arrow_scale, nu_direction[2]*arrow_scale, width=1.2, head_width=4, color='C1', label="Shower Direction") # Add a text label near the neutrino event plt.text(-190, -105, "Neutrino\nEvent", fontsize=13, horizontalalignment='center') # To represent the approximate radio signal power coming from the shower # we'll draw many overlapping, nearly-transparent lines. The density of lines # should increase as we get closer to the Cherenkov angle so that the image # gives the impression of higher signal power nearer the Cherenkov angle kwargs = {'color': 'C0', 'alpha': 0.015} for dtheta in np.sqrt(np.logspace(-5, -1, 400)): for sign in (-1, 1): # Calculate the Cherenkov angle directions +/- some offset dir1 = np.dot(rot(theta_c+sign*dtheta), nu_direction) dir2 = np.dot(nu_direction, rot(theta_c+sign*dtheta)) # Plot short lines in the directions defined above, and pass the # keyword arguments defined above for color and opacity plt.plot([nu_vertex[0], nu_vertex[0]+dir1[0]*2*arrow_scale], [nu_vertex[2], nu_vertex[2]+dir1[2]*2*arrow_scale], **kwargs) plt.plot([nu_vertex[0], nu_vertex[0]+dir2[0]*2*arrow_scale], [nu_vertex[2], nu_vertex[2]+dir2[2]*2*arrow_scale], **kwargs) # Add diamonds at the two antenna locations plt.scatter(ant1[0], ant1[2], marker='D', s=75, color='C2', zorder=10, label="Antenna") plt.scatter(ant2[0], ant2[2], marker='D', s=75, color='C2', zorder=10) # Add text labels near the antennas plt.text(0, -22, "Antenna\n(z = -2 m)", fontsize=13, horizontalalignment='center') plt.text(0, -80, "Antenna\n(z = -60 m)", fontsize=13, horizontalalignment='center') # Plot the paths the radio signal takes from the shower vertex to the antennas plt.plot(x1, z1, color='C0', label="Signal Path") plt.plot(x2, z2, color='C0') # Add a text label between the signal paths plt.text(-105, -70, "Signal Paths", fontsize=13, horizontalalignment='center') # Set the axis limits to the previously set values plt.xlim(*xlims) plt.ylim(*ylims) # Turn off axis ticks and labels on both the x and y axes plt.tick_params(bottom=False, labelbottom=False, left=False, labelleft=False) # Turn off the axis lines around the plot plt.gca().axis('off') # Save the image as a PDF so it scales nicely in presentations plt.savefig('event_diagram.pdf')