Commit 89690194 authored by JGlombitza's avatar JGlombitza
Browse files

v1.1

New Path
Correct description
parent 584ff2f2
...@@ -19,8 +19,8 @@ nbatch = 132 ...@@ -19,8 +19,8 @@ nbatch = 132
# load and prepare data # load and prepare data
# ---------------------------------------------------------------------- # ----------------------------------------------------------------------
Folder = '/home/JGlombitza/DataLink/Public/Paper/PreProData' Folder = '../Data_preprocessed'
filenames=glob.glob(Folder + "/Data_PlanarWaveShower_PrePro_NewMC_Mixed_Composition*") filenames=glob.glob(Folder + "/PlanarWaveShower_PrePro*")
N = 10000*len(filenames) N = 10000*len(filenames)
a = 10000 # packagesize a = 10000 # packagesize
...@@ -30,10 +30,10 @@ showeraxis = np.zeros(N*3).reshape(N,3) ...@@ -30,10 +30,10 @@ showeraxis = np.zeros(N*3).reshape(N,3)
Energy = np.zeros(N) Energy = np.zeros(N)
for i in range(0,len(filenames)): for i in range(0,len(filenames)):
data = np.load(filenames[i]) data = np.load(filenames[i])
Input1[a*i:a*(i+1)] = data['Input1'].reshape(a,9,9,80,1) Input1[a*i:a*(i+1)] = data['Input1']
showeraxis[a*i:a*(i+1)] = data['showeraxis'].reshape(a,3) showeraxis[a*i:a*(i+1)] = data['showeraxis']
Input2[a*i:a*(i+1)] = data['Input2'].reshape(a,9,9,2) Input2[a*i:a*(i+1)] = data['Input2']
Energy[a*i:a*(i+1)] = data['Energy'].reshape(a) Energy[a*i:a*(i+1)] = data['Energy']
print('Loaded %i showers' %N ) print('Loaded %i showers' %N )
......
...@@ -18,8 +18,8 @@ nbatch = 132 ...@@ -18,8 +18,8 @@ nbatch = 132
# load and prepare data # load and prepare data
# ---------------------------------------------------------------------- # ----------------------------------------------------------------------
Folder = '/home/JGlombitza/DataLink/Public/Paper/PreProData' Folder = '../Data_preprocessed'
filenames=glob.glob(Folder + "/Data_PlanarWaveShower_PrePro_NewMC_Mixed_Composition*") filenames=glob.glob(Folder + "/PlanarWaveShower_PrePro*")
N = 10000*len(filenames) N = 10000*len(filenames)
a = 10000 # packagesize a = 10000 # packagesize
...@@ -28,9 +28,9 @@ showeraxis = np.zeros(N*3).reshape(N,3) ...@@ -28,9 +28,9 @@ showeraxis = np.zeros(N*3).reshape(N,3)
Energy = np.zeros(N) Energy = np.zeros(N)
for i in range(0,len(filenames)): for i in range(0,len(filenames)):
data = np.load(filenames[i]) data = np.load(filenames[i])
showeraxis[a*i:a*(i+1)] = data['showeraxis'].reshape(a,3) showeraxis[a*i:a*(i+1)] = data['showeraxis']
Input2[a*i:a*(i+1)] = data['Input2'][:,:,:,0].reshape(a,9,9,1) Input2[a*i:a*(i+1)] = data['Input2'][:,:,:,0].reshape(a,9,9,1)
Energy[a*i:a*(i+1)] = data['Energy'].reshape(a) Energy[a*i:a*(i+1)] = data['Energy']
## Reshape, Split to test and train sets ## Reshape, Split to test and train sets
ntest = 40000 ntest = 40000
......
...@@ -13,16 +13,12 @@ lr = 0.001 ...@@ -13,16 +13,12 @@ lr = 0.001
log_dir="." log_dir="."
nbatch = 132 nbatch = 132
def DenselyConnectedSepConv(z, nfilter, **kwargs):
c = SeparableConv2D(nfilter, (3, 3), padding = 'same', depth_multiplier=1, **kwargs)(z)
return concatenate([z, c], axis=-1)
# ---------------------------------------------------------------------- # ----------------------------------------------------------------------
# load and prepare data # load and prepare data
# ---------------------------------------------------------------------- # ----------------------------------------------------------------------
Folder = '/home/JGlombitza/DataLink/Public/Paper/PreProData' Folder = '../Data_preprocessed'
filenames=glob.glob(Folder + "/Data_PlanarWaveShower_PrePro_NewMC_Mixed_Composition*") filenames=glob.glob(Folder + "/PlanarWaveShower_PrePro*")
N = 10000*len(filenames) N = 10000*len(filenames)
a = 10000 # packagesize a = 10000 # packagesize
...@@ -31,9 +27,9 @@ Input2 = np.zeros(N*9*9*2).reshape(N,9,9,2) ...@@ -31,9 +27,9 @@ Input2 = np.zeros(N*9*9*2).reshape(N,9,9,2)
Energy = np.zeros(N) Energy = np.zeros(N)
for i in range(0,len(filenames)): for i in range(0,len(filenames)):
data = np.load(filenames[i]) data = np.load(filenames[i])
Input1[a*i:a*(i+1)] = data['Input1'].reshape(a,9,9,80,1) Input1[a*i:a*(i+1)] = data['Input1']
Energy[a*i:a*(i+1)] = data['Energy'] Energy[a*i:a*(i+1)] = data['Energy']
Input2[a*i:a*(i+1)] = data['Input2'].reshape(a,9,9,2) Input2[a*i:a*(i+1)] = data['Input2']
## Reshape, Split to test and train sets ## Reshape, Split to test and train sets
ntest = 40000 ntest = 40000
...@@ -67,7 +63,7 @@ x = concatenate([TimeProcess, inputT0]) ...@@ -67,7 +63,7 @@ x = concatenate([TimeProcess, inputT0])
# Densely connected convolutions # Densely connected convolutions
nfilter =12 nfilter =12
for i in range(4): for i in range(4):
x = DenselyConnectedSepConv(x, nfilter, **kwargs) x = tools.DenselyConnectedSepConv(x, nfilter, **kwargs)
nfilter = 2*nfilter nfilter = 2*nfilter
x = SeparableConv2D(nfilter, (3, 3), padding = 'same', depth_multiplier=1, **kwargs)(x) x = SeparableConv2D(nfilter, (3, 3), padding = 'same', depth_multiplier=1, **kwargs)(x)
x = Flatten()(x) x = Flatten()(x)
......
...@@ -18,8 +18,8 @@ nbatch = 132 ...@@ -18,8 +18,8 @@ nbatch = 132
# load and prepare data # load and prepare data
# ---------------------------------------------------------------------- # ----------------------------------------------------------------------
Folder = '/home/JGlombitza/DataLink/Public/Paper/PreProData' Folder = '../Data_preprocessed'
filenames=glob.glob(Folder + "/Data_PlanarWaveShower_PrePro_NewMC_Mixed_Composition*") filenames=glob.glob(Folder + "/PlanarWaveShower_PrePro*")
N = 10000*len(filenames) N = 10000*len(filenames)
a = 10000 # packagesize a = 10000 # packagesize
......
...@@ -17,8 +17,8 @@ nbatch = 132 ...@@ -17,8 +17,8 @@ nbatch = 132
# load and prepare data # load and prepare data
# ---------------------------------------------------------------------- # ----------------------------------------------------------------------
Folder = '/home/JGlombitza/DataLink/Public/Paper/PreProData' Folder = '../Data_preprocessed'
filenames=glob.glob(Folder + "/Data_PlanarWaveShower_PrePro_NewMC_Mixed_Composition*") filenames=glob.glob(Folder + "/PlanarWaveShower_PrePro*")
N = 10000*len(filenames) N = 10000*len(filenames)
a = 10000 # packagesize a = 10000 # packagesize
...@@ -28,10 +28,10 @@ Xmax = np.zeros(N) ...@@ -28,10 +28,10 @@ Xmax = np.zeros(N)
Energy = np.zeros(N) Energy = np.zeros(N)
for i in range(0,len(filenames)): for i in range(0,len(filenames)):
data = np.load(filenames[i]) data = np.load(filenames[i])
Input1[a*i:a*(i+1)] = data['Input1'].reshape(a,9,9,80,1) Input1[a*i:a*(i+1)] = data['Input1']
Xmax[a*i:a*(i+1)] = data['Xmax'] Xmax[a*i:a*(i+1)] = data['Xmax']
Energy[a*i:a*(i+1)] = data['Energy'] Energy[a*i:a*(i+1)] = data['Energy']
Input2[a*i:a*(i+1)] = data['Input2'].reshape(a,9,9,2) Input2[a*i:a*(i+1)] = data['Input2']
## Reshape, Split to test and train sets ## Reshape, Split to test and train sets
ntest = 40000 ntest = 40000
......
...@@ -13,13 +13,17 @@ def Distance(y_true, y_pred): ...@@ -13,13 +13,17 @@ def Distance(y_true, y_pred):
mean, var = tf.nn.moments((y_true-y_pred), axes=[0]) mean, var = tf.nn.moments((y_true-y_pred), axes=[0])
return tf.sqrt(var) return tf.sqrt(var)
def MeanElong(y_true, y_pred): ##inGrad
def MeanElong(y_true, y_pred):
''' Metric to control the mean of the angulardistance distribution '''
return 57.29577951308232*K.mean(tf.minimum(tf.acos((y_pred[:,0]*y_true[:,0] + y_pred[:,1]*y_true[:,1] + y_pred[:,2]*y_true[:,2])/(tf.sqrt(y_pred[:,0]*y_pred[:,0]+y_pred[:,1]*y_pred[:,1]+y_pred[:,2]*y_pred[:,2])*tf.sqrt(y_true[:,0]*y_true[:,0]+y_true[:,1]*y_true[:,1]+y_true[:,2]*y_true[:,2]))),1)) return 57.29577951308232*K.mean(tf.minimum(tf.acos((y_pred[:,0]*y_true[:,0] + y_pred[:,1]*y_true[:,1] + y_pred[:,2]*y_true[:,2])/(tf.sqrt(y_pred[:,0]*y_pred[:,0]+y_pred[:,1]*y_pred[:,1]+y_pred[:,2]*y_pred[:,2])*tf.sqrt(y_true[:,0]*y_true[:,0]+y_true[:,1]*y_true[:,1]+y_true[:,2]*y_true[:,2]))),1))
def Percentile68(x): def Percentile68(x):
''' Calculates 68% quantile of distribution''' ''' Calculates 68% quantile of distribution'''
return np.percentile(x, 68) return np.percentile(x, 68)
def vec2ang(v): def vec2ang(v):
'''Calculates phi and theta for given shower axis''' '''Calculates phi and theta for given shower axis'''
x, y, z = np.asarray(v).T x, y, z = np.asarray(v).T
...@@ -27,11 +31,13 @@ def vec2ang(v): ...@@ -27,11 +31,13 @@ def vec2ang(v):
theta = np.pi/2 - np.arctan2(z, (x**2 + y**2)**.5) theta = np.pi/2 - np.arctan2(z, (x**2 + y**2)**.5)
return np.rad2deg(phi), np.rad2deg(theta) return np.rad2deg(phi), np.rad2deg(theta)
def DenselyConnectedSepConv(z, nfilter, **kwargs): def DenselyConnectedSepConv(z, nfilter, **kwargs):
''' Densely Connected SeparableConvolution2D Layer''' ''' Densely Connected SeparableConvolution2D Layer'''
c = SeparableConv2D(nfilter, (3, 3), padding = 'same', depth_multiplier=1, **kwargs)(z) c = SeparableConv2D(nfilter, (3, 3), padding = 'same', depth_multiplier=1, **kwargs)(z)
return concatenate([z, c], axis=-1) return concatenate([z, c], axis=-1)
def PlotHistosXmaxReconstruction(y_true, y_pred, log_dir, name, Energy): def PlotHistosXmaxReconstruction(y_true, y_pred, log_dir, name, Energy):
''' Plottingfunction to evaluate the reconstruction of the shower depths Xmax ''' Plottingfunction to evaluate the reconstruction of the shower depths Xmax
plots: predicted Xmax against MC Xmax and the energy dependence of the reconstruction (embedded)''' plots: predicted Xmax against MC Xmax and the energy dependence of the reconstruction (embedded)'''
...@@ -41,7 +47,7 @@ def PlotHistosXmaxReconstruction(y_true, y_pred, log_dir, name, Energy): ...@@ -41,7 +47,7 @@ def PlotHistosXmaxReconstruction(y_true, y_pred, log_dir, name, Energy):
ax1.scatter(y_true, y_pred, s = 2, color='royalblue', alpha=.90) ax1.scatter(y_true, y_pred, s = 2, color='royalblue', alpha=.90)
ax1.set_xlabel(r'$X_{max,\;true}$ $[g/cm^{2}]$') ax1.set_xlabel(r'$X_{max,\;true}$ $[g/cm^{2}]$')
ax1.set_ylabel(r'$X_{max,\;DNN}$ $[g/cm^{2}]$') ax1.set_ylabel(r'$X_{max,\;DNN}$ $[g/cm^{2}]$')
left, bottom, width, height = [0.23, 0.6, 0.3, 0.26] # Links Oben left, bottom, width, height = [0.23, 0.6, 0.3, 0.26]
ax1.text(0.95,0.2, r"$\mu = $ %.3f " % np.mean(reco) + r"$[g/cm^{2}]$" + "\n" + "$\sigma = $ %.3f " % np.std(reco) + r"$[g/cm^{2}]$", verticalalignment = 'top', horizontalalignment = 'right', transform=ax1.transAxes, backgroundcolor='w') ax1.text(0.95,0.2, r"$\mu = $ %.3f " % np.mean(reco) + r"$[g/cm^{2}]$" + "\n" + "$\sigma = $ %.3f " % np.std(reco) + r"$[g/cm^{2}]$", verticalalignment = 'top', horizontalalignment = 'right', transform=ax1.transAxes, backgroundcolor='w')
ax1.plot([np.min(y_true), np.max(y_true)], [np.min(y_true), np.max(y_true)], color='red') ax1.plot([np.min(y_true), np.max(y_true)], [np.min(y_true), np.max(y_true)], color='red')
# Energy dependence plot # Energy dependence plot
...@@ -59,14 +65,13 @@ def PlotHistosXmaxReconstruction(y_true, y_pred, log_dir, name, Energy): ...@@ -59,14 +65,13 @@ def PlotHistosXmaxReconstruction(y_true, y_pred, log_dir, name, Energy):
fig.savefig(log_dir + '/Xmax_embedding_%s.png' %name) fig.savefig(log_dir + '/Xmax_embedding_%s.png' %name)
return reco return reco
def PlotHistosEnergyReconstruction(y_true, y_pred, log_dir, name): def PlotHistosEnergyReconstruction(y_true, y_pred, log_dir, name):
''' Plottingfunction to evaluate the reconstruction of the cosmic ray energy ''' Plottingfunction to evaluate the reconstruction of the cosmic ray energy
plots: predicted Energy against MC Energy and the energy dependence of the relative energy resolution (embedded)''' plots: predicted Energy against MC Energy and the energy dependence of the relative energy resolution (embedded)'''
y_true = y_true.reshape(y_true.shape[0],1) y_true = y_true.reshape(y_true.shape[0],1)
reco = y_true-y_pred reco = y_true-y_pred
# Embedding plot
# Embedding plot
fig, ax1 = plt.subplots(1) fig, ax1 = plt.subplots(1)
ax1.scatter(y_true, y_pred, s = 2, color='royalblue', alpha=.90) ax1.scatter(y_true, y_pred, s = 2, color='royalblue', alpha=.90)
ax1.set_xlabel(r'$Energy_{true}\;[EeV]$') ax1.set_xlabel(r'$Energy_{true}\;[EeV]$')
...@@ -87,17 +92,18 @@ def PlotHistosEnergyReconstruction(y_true, y_pred, log_dir, name): ...@@ -87,17 +92,18 @@ def PlotHistosEnergyReconstruction(y_true, y_pred, log_dir, name):
fig.savefig(log_dir + '/Energy_embedding_%s.png' %name) fig.savefig(log_dir + '/Energy_embedding_%s.png' %name)
return reco return reco
def PlotHistosAngularReconstruction(y_true, y_pred, log_dir, name, Energy, Signal=0): def PlotHistosAngularReconstruction(y_true, y_pred, log_dir, name, Energy, Signal=0):
''' Plotting function to evaluate the reconstruction of the shower axis ''' Plotting function to evaluate the reconstruction of the shower axis
plots: angular distance distribution and the energy dependence of the angular reconstruction (embedded)''' plots: angular distance distribution and the energy dependence of the angular reconstruction (embedded)'''
phi, theta = vec2ang(y_true) phi, theta = vec2ang(y_true)
angulardistance = 180 / np.pi * np.arccos(np.clip(np.sum(y_true*y_pred, axis=-1)/np.linalg.norm(y_pred, axis=-1)/np.linalg.norm(y_true, axis=-1),-1,1)) angulardistance = 180 / np.pi * np.arccos(np.clip(np.sum(y_true*y_pred, axis=-1)/np.linalg.norm(y_pred, axis=-1)/np.linalg.norm(y_true, axis=-1),-1,1))
resolution = Percentile68(angulardistance) resolution = Percentile68(angulardistance)
# Embedding plot # Embedding plot
fig, ax1 = plt.subplots(1) fig, ax1 = plt.subplots(1)
ax1.hist(angulardistance, bins=np.linspace(0, 4, 51)) ax1.hist(angulardistance, bins=np.linspace(0, 4, 51))
ax1.set(xlabel=r'$Angulardistance\;\delta\;[^\circ]$', ylabel=r'$N$') ax1.set(xlabel=r'$Angulardistance\;\delta\;[^\circ]$', ylabel=r'$N$')
left, bottom, width, height = [0.56, 0.6, 0.31, 0.27] # Rechts Oben left, bottom, width, height = [0.56, 0.6, 0.31, 0.27]
ax2 = fig.add_axes([left, bottom, width, height]) ax2 = fig.add_axes([left, bottom, width, height])
x_bins = 12 x_bins = 12
logE = 18 + np.log10(Energy) logE = 18 + np.log10(Energy)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment