ZH: rit.py: simple edits for saving and c other X

This commit is contained in:
Eric Teunis de Boone 2022-11-22 18:06:12 +01:00
parent e6379e2cd2
commit f905dfb2f3

View file

@ -146,10 +146,10 @@ def dist_to_line_sum(param,data,weights):
# print('%.2e %.2e %.2e %.2e %.2e'%(x0,y0,theta,phi,dsum)) # print('%.2e %.2e %.2e %.2e %.2e'%(x0,y0,theta,phi,dsum))
return dsum/len(data) return dsum/len(data)
def get_axis_points(e,savefig=True,path="",zgr=0): def get_axis_points(e,savefig=True,path="",zgr=0,Xlow=300, Xhigh=1000, N_X=15):
axis_points = [] axis_points = []
max_vals = [] max_vals = []
Xsteps = np.linspace(300,1000,15) Xsteps = np.linspace(Xlow, Xhigh, N_X)
zgr=zgr+e.core[2] #not exact zgr=zgr+e.core[2] #not exact
dXref = atm.distance_to_slant_depth(np.deg2rad(e.zenith),750,zgr) dXref = atm.distance_to_slant_depth(np.deg2rad(e.zenith),750,zgr)
scale2d = dXref*np.tan(np.deg2rad(2.)) scale2d = dXref*np.tan(np.deg2rad(2.))
@ -311,14 +311,14 @@ def fill_stations_propeties(e,res):
z = np.array([a.z for a in e.antennas]) z = np.array([a.z for a in e.antennas])
ids = [a.name for a in e.antennas] ids = [a.name for a in e.antennas]
maxE = np.array([np.max(a.E_AxB) for a in e.antennas]) maxE = np.array([np.max(a.E_AxB) for a in e.antennas])
has_pulse = np.array([np.max(a.has_pulse) for a in e.antennas]) #has_pulse = np.array([np.max(a.has_pulse) for a in e.antennas])
res.station_x.append(x) res.station_x.append(x)
res.station_y.append(y) res.station_y.append(y)
res.station_z.append(z) res.station_z.append(z)
res.station_ids.append(ids) res.station_ids.append(ids)
res.has_pulse.append(has_pulse) #res.has_pulse.append(has_pulse)
def reconstruction(e,outfile=''): def reconstruction(e,outfile='', slice_outdir=None):
res = RITResult() res = RITResult()
res.isMC.append(True) res.isMC.append(True)
res.zenith_ini.append(e.zenith) res.zenith_ini.append(e.zenith)
@ -329,7 +329,7 @@ def reconstruction(e,outfile=''):
#only use signal that have a signal in data #only use signal that have a signal in data
fill_stations_propeties(e,res) fill_stations_propeties(e,res)
Xs,axis_points,max_vals = get_axis_points(e,False) Xs,axis_points,max_vals = get_axis_points(e,slice_outdir is not None, slice_outdir, 200, 1000, 5)
zen,azi,core = fit_track(e,axis_points,max_vals,1e2) zen,azi,core = fit_track(e,axis_points,max_vals,1e2)
fig = figure_3D(axis_points,max_vals,zen,azi,core,res) fig = figure_3D(axis_points,max_vals,zen,azi,core,res)
fig.savefig(outfile) fig.savefig(outfile)
@ -349,10 +349,13 @@ if __name__ == "__main__":
file = '../ZH_airshower/mysim.sry' file = '../ZH_airshower/mysim.sry'
ev = REvent(file) ev = REvent(file)
set_pol_and_bp(ev) set_pol_and_bp(ev)
X = 750
dXref = atm.distance_to_slant_depth(np.deg2rad(ev.zenith),750,0) dXref = atm.distance_to_slant_depth(np.deg2rad(ev.zenith),X,0)
scale2d = dXref*np.tan(np.deg2rad(2.)) scale2d = dXref*np.tan(np.deg2rad(2.))
xx,yy,p,km= shower_plane_slice(ev,700,21,21,scale2d,scale2d) xx,yy,p,km= shower_plane_slice(ev,X,21,21,scale2d,scale2d)
plt.scatter(xx,yy,c=p)
plt.colorbar() slice_figure(ev,X,xx,yy,p,mode='sp')
#plt.scatter(xx,yy,c=p)
#plt.colorbar()
plt.show() plt.show()