ZH: update rit.py from Harm

This commit is contained in:
Eric Teunis de Boone 2022-11-22 13:41:08 +01:00
parent f88799dd95
commit 8b210514d3

View file

@ -60,7 +60,7 @@ def pow_and_time(test_loc,ev,dt=1.0):
P = 0
return P,t_,a_,a_sum,t_sum
def shower_axis_slice(e,Xb=200,Xe=1200,dX=2,zgr=1400):
def shower_axis_slice(e,Xb=200,Xe=1200,dX=2,zgr=0):
zgr = zgr + e.core[2]
N = int((Xe-Xb)/dX)
Xs = np.array(np.linspace(Xb,Xe,N+1))
@ -79,7 +79,7 @@ def shower_axis_slice(e,Xb=200,Xe=1200,dX=2,zgr=1400):
p = np.asanyarray(p)
return ds,Xs,locs,p
def shower_plane_slice(e,X=750.,Nx=10,Ny=10,wx=1e3,wy=1e3,xoff=0,yoff=0,zgr=1400):
def shower_plane_slice(e,X=750.,Nx=10,Ny=10,wx=1e3,wy=1e3,xoff=0,yoff=0,zgr=0):
zgr = zgr + e.core[2]
dX = atm.distance_to_slant_depth(np.deg2rad(e.zenith),X,zgr)
x = np.linspace(-wx,wx,Nx)
@ -106,7 +106,7 @@ def slice_figure(e,X,xx,yy,p,mode='horizontal'):
fig.suptitle(r'E = %.1f EeV, $\theta$ = %.1f$^\circ$, $\phi$ = %.1f$^\circ$ X = %.f'%(e.energy,e.zenith,e.azimuth,X))
sc = axs.scatter(xx/1e3,yy/1e3,c=p,cmap='Spectral_r',alpha=0.6)
fig.colorbar(sc,ax=axs)
zgr = 1400 + e.core[2]
zgr = 0 + e.core[2]
dX = atm.distance_to_slant_depth(np.deg2rad(e.zenith),X,zgr)
xc = np.sin(np.deg2rad(e.zenith))*np.cos(np.deg2rad(e.azimuth))* dX
yc = np.sin(np.deg2rad(e.zenith))*np.sin(np.deg2rad(e.azimuth))* dX
@ -146,7 +146,7 @@ def dist_to_line_sum(param,data,weights):
# print('%.2e %.2e %.2e %.2e %.2e'%(x0,y0,theta,phi,dsum))
return dsum/len(data)
def get_axis_points(e,savefig=True,path="",zgr=1400):
def get_axis_points(e,savefig=True,path="",zgr=0):
axis_points = []
max_vals = []
Xsteps = np.linspace(300,1000,15)
@ -344,3 +344,15 @@ def reconstruction(e,outfile=''):
res.zenith_rec.append(e.zenith)
res.dcore_rec.append(core)
return res
if __name__ == "__main__":
file = '../ZH_airshower/mysim.sry'
ev = REvent(file)
set_pol_and_bp(ev)
dXref = atm.distance_to_slant_depth(np.deg2rad(ev.zenith),750,0)
scale2d = dXref*np.tan(np.deg2rad(2.))
xx,yy,p,km= shower_plane_slice(ev,700,21,21,scale2d,scale2d)
plt.scatter(xx,yy,c=p)
plt.colorbar()
plt.show()