From fdd2517111fb68443a57d394dd17ad04cf686e79 Mon Sep 17 00:00:00 2001 From: jorge Date: Mon, 17 Mar 2025 11:17:53 -0600 Subject: [PATCH 1/6] Added chasm run scripts. Pending further integration --- src/nuspacesim/chasm_run.py | 363 +++++++++++++++++ src/nuspacesim/chasm_to_offline_original.py | 420 ++++++++++++++++++++ 2 files changed, 783 insertions(+) create mode 100644 src/nuspacesim/chasm_run.py create mode 100644 src/nuspacesim/chasm_to_offline_original.py diff --git a/src/nuspacesim/chasm_run.py b/src/nuspacesim/chasm_run.py new file mode 100644 index 0000000..4183035 --- /dev/null +++ b/src/nuspacesim/chasm_run.py @@ -0,0 +1,363 @@ +#!@Python_EXECUTABLE@ + +import sys +import numpy as np +from scipy.spatial.transform import Rotation as R +from scipy import interpolate +from astropy.table import QTable, Table, Column +from astropy import units as u + +try: + import CHASM as ch +except ImportError as e: + print("Module CHASM could not be imported. It is most likely not installed. \n \ + It can be installed through pip: python -m pip install CHASM-NuSpacesim. \n The detailed error message is: \n",e) + sys.exit() + + +def get_vectors(zenith, azimuth, altitude, grid_angle, npoints): + + # get distance along axis corresponding to detector altitude + r = sim.ingredients['axis'].h_to_axis_R_LOC(altitude*1.e3, zenith) + grid_width = r*np.tan(np.radians(grid_angle)) + orig = [0, 0, r] + rcoord_count = np.linspace(0, grid_width, npoints) + ang_step=np.arange(0,355,10) + nperiods=(npoints // ang_step.size) + 1 + # Radial random angle + ang = np.tile(np.radians(ang_step),nperiods)[:npoints] + xall = rcoord_count*np.cos(ang) + yall = rcoord_count*np.sin(ang) + zall = np.full_like(xall, r) + + # Square grid + # x = np.linspace(-grid_width, grid_width, n_side) + # y = np.linspace(-grid_width, grid_width, n_side) + # xx, yy = np.meshgrid(x,y) + # zz = np.full_like(xx, r) #convert altitude to m + + # vecs = np.vstack((xx.flatten(),yy.flatten(),zz.flatten())).T + vecs = np.vstack((xall.flatten(), yall.flatten(), zall.flatten())).T + + theta_rot_axis = np.array([0, 1, 0]) + theta_rotation = R.from_rotvec(theta_rot_axis * zenith) + + z_rot_axis = np.array([0, 0, 1]) + z_rotation = R.from_rotvec(z_rot_axis * np.pi/2) + + phi_rot_axis = np.array([0, 0, 1]) + phi_rot = R.from_rotvec(phi_rot_axis*azimuth) + + vecs = z_rotation.apply(vecs) + vecs = theta_rotation.apply(vecs) + + tel_vecs = phi_rot.apply(vecs) + orig = z_rotation.apply(orig) + orig = theta_rotation.apply(orig) + orig = phi_rot.apply(orig) + return tel_vecs, orig + + +def run_chasm(sim, orig): + sig = sim.run(mesh=False, att=True) + del sim + phot = np.array(sig.photons) + phot_lambda = phot.sum(axis=2).sum(axis=0) + phot_plane = phot.sum(axis=1) + source_vec = sig.source_points #Vectors to axis points + costheta=sig.cos_theta # This returns the cosine of the angle between the z-axis and the vector from the axis to the counter + times = np.array(sig.times) + times = times-np.min(times) + phot_plane[phot_plane < 0] = 0 + + total_propagation_times = times + photons_on_plane = phot_plane + wlbins = sig.wavelengths + wlhist = phot_lambda + return source_vec, total_propagation_times, costheta, photons_on_plane, wlbins, wlhist + +def signal_to_astropy(sig: ShowerSignal) -> QTable: + '''This function outputs the data in a shower signal object to an astropy table. + ''' + column_list = [] + column_list.append(Column(sig.source_points, name='source points',unit=u.m)) + column_list.append(Column(sig.charged_particles, name='charged particles',unit=u.ct)) + column_list.append(Column(sig.depths, name='depths',unit=u.g/u.cm**2)) + for i in range(sig.photons.shape[0]): + column_list.append(Column(sig.photons[i].T, name=f'counter {i} photons',unit=u.ct)) + column_list.append(Column(sig.times[i].T, name=f'counter {i} arrival times',unit=u.nanosecond)) + column_list.append(Column(sig.cos_theta[i].T, name=f'counter {i} cos zenith',unit=u.rad)) + return QTable(column_list) + +def weighted_percentile(data, percents, weights=None): + ''' percents in units of 1% + weights specifies the frequency (count) of data. + ''' + if weights is None: + return np.percentile(data, percents) + ind = np.argsort(data) + d = data[ind] + w = weights[ind] + p = 1.*w.cumsum()/w.sum()*100 + y = np.interp(percents, p, d) + return y + + +def theta_D_from_theta_EE(det_alt, theta_EE): + """ + det_alt : float + detector altitude in km + theta_EE : float + earth emergence angle in rad + + Returns + ------- + float + detector viewing angle in degrees + """ + theta_EE = theta_EE * (np.pi / 180) + Re = 6371 + theta_D = (180 / np.pi) * np.arcsin((Re / (Re + det_alt)) + * np.sin(np.pi / 2 + theta_EE)) + + return theta_D + + +def h_to_axis_R_LOC(h, theta) -> np.ndarray: + earth_radius = 6371e3 + ground_level = 0 + """Return the length along the shower axis from the point of Earth + emergence to the height above the surface specified + + Parameters: + h: array of heights (m above ground level) + theta: polar angle of shower axis (radians) + + returns: r (m) (same size as h), an array of distances along the shower + axis_sp. + """ + cos_EM = np.cos(np.pi-theta) + R = earth_radius + ground_level + r_CoE = h + R # distance from the center of the earth to the specified height + rs = R*cos_EM + np.sqrt(R**2*cos_EM**2-R**2+r_CoE**2) + # rs -= rs[0] + # rs[0] = 1. + return rs + +zen90 = tshower["zenith"].array() +phi = np.array(np.radians(tshower["azimuth"].array())) +Hfirst = tshower["Hfirst"].array()/1000 +Xfirst = tshower["Xfirst"].array() +theta = np.array(np.radians(180-zen90)) +nEvents = np.size(zen90) +lgenergy = tshower["lgE"].array() +X = tshower["X"].array() +N = tshower["N"].array() + + + +# add shower axis +zenith = theta +azimuth = phi + +# add grid of detectors +n_side = 30 # For square grid +npoints = args.n # For radial +detector_grid_alt = args.alt # km +grid_angle = 5 # degrees +filename = args.outfile +eng = np.zeros(1, dtype=float) +zen = np.zeros(1, dtype=float) +azi = np.zeros(1, dtype=float) +startdist = np.zeros(1, dtype=float) +X0 = np.zeros(1, dtype=float) + + + +for i in range(nEvents): + sim = ch.ShowerSimulation() + sim.add(ch.UpwardAxis(zenith[i], azimuth[i], curved=True)) + tel_vecs, orig = get_vectors( + zenith[i], azimuth[i], detector_grid_alt, grid_angle, npoints) + + sim.add(ch.UserShower(np.array(X[i]), np.array(N[i]))) + sim.add(ch.SphericalCounters(tel_vecs, np.sqrt(1/np.pi))) + sim.add(ch.Yield(270, 1000, N_bins=100)) + print('Running CHASM simulation shower {}'.format(i)) + + plane_vec, total_propagation_times, photon_arrival_phZenith, photons_on_plane, wlbins, wlhist, dist_counter = run_chasm( + sim, orig) + r_coord = np.sqrt((plane_vec**2).sum(axis=-1)) + r_coordinates = (r_coord/1000).flatten() + + #rotate plane_vec to have it in z=0 + + theta_rot_axis = np.array([0,1,0]) + theta_rotation = R.from_rotvec(theta_rot_axis * -zenith[i]) + + z_rot_axis = np.array([0,0,1]) + z_rotation = R.from_rotvec(z_rot_axis * -np.pi/2) + + phi_rot_axis=np.array([0,0,1]) + phi_rot=R.from_rotvec(phi_rot_axis*-azimuth[i]) + + plane_vec_rot=phi_rot.apply(plane_vec) + plane_vec_rot=theta_rotation.apply(plane_vec_rot) + plane_vec_rot=z_rotation.apply(plane_vec_rot) + + x_plane=plane_vec_rot[:,0] + y_plane=plane_vec_rot[:,1] + + azimuth_angles = 0*(180/np.pi)*np.arctan2(y_plane, x_plane) + azi_below = azimuth_angles<-180 + azi_above = azimuth_angles>180 + azimuth_angles[azi_below] = 360+azimuth_angles[azi_below] + azimuth_angles[azi_above] = -360+azimuth_angles[azi_above] + + energy = 10**lgenergy[i] + azimdeg = np.degrees(azimuth[i]) + theta_D = theta_D_from_theta_EE(detector_grid_alt, 90-np.degrees(theta[i])) + L_particle = Hfirst[i] # meters + X0_First=Xfirst[i] + + photons_on_plane[photons_on_plane < 0] = 0 + + bins, start, stop = 100, 0, 4 # For time + time_bins = np.logspace(float(start), float(stop), num=int(bins)) + bins, start, stop = 100, -3, 1 # For Zenith + phZenith_bins = np.logspace(float(start), float(stop), num=int(bins)) + bins, start, stop = 1440, -180, 180 # For Azimuth + bin_size = (float(stop) - float(start)) / float(bins) + phAzi_bins = np.arange(float(start), float(stop), bin_size) + spatial_bins = 100 + pct_cap = 100 # 30% recommended in easchersim, but 50 used + r_percentiles = weighted_percentile( + r_coordinates, pct_cap, weights=photons_on_plane) + r_bounds = weighted_percentile( + r_coordinates, np.array([1, 99.7]), weights=photons_on_plane) + r_bins = np.linspace(np.min(r_coordinates)-1e-7, + r_percentiles, spatial_bins) + ncounter_rbin = np.zeros(spatial_bins-1) + + time_bounds, zenith_bounds = [], [] + offset_propagation_times = total_propagation_times + offset_arrival_phZenith = photon_arrival_phZenith + for j in range(0, len(r_bins) - 1): + # Find coordinates which are inside each spatial bin + r_lower, r_upper = r_bins[j], r_bins[j + 1] + inbin = (dist_counter <= r_upper) & (dist_counter > r_lower) + ncounter_rbin[j] = int(np.sum(inbin)) + in_range = (r_coordinates >= r_lower) & (r_coordinates <= r_upper) + in_range_time, in_range_zenith, in_range_photons = total_propagation_times[ + in_range], photon_arrival_phZenith[in_range], photons_on_plane[in_range] + # Calculate the 1% and 99% weighted percentiles of the arrival time and arrival angle inside the bin + time_percentiles = weighted_percentile( + in_range_time, np.array([1, 99]), weights=in_range_photons) + zenith_percentiles = weighted_percentile( + in_range_zenith, np.array([1, 99]), weights=in_range_photons) + + time_bounds.append(time_percentiles) + zenith_bounds.append(zenith_percentiles) + + time_offset = interpolate.interp1d( + r_bins[:-1], np.array(time_bounds)[:, 0], kind='cubic', bounds_error=False, fill_value=0) + phZenith_offset = interpolate.interp1d( + r_bins[:-1], np.array(zenith_bounds)[:, 0], kind='cubic', bounds_error=False, fill_value=0) + + offset_propagation_times = total_propagation_times - \ + time_offset(r_coordinates).reshape(total_propagation_times.shape) + onephotonmask = (photons_on_plane >= 1) + offset_propagation_times -= np.min(offset_propagation_times[onephotonmask]) + + offset_arrival_phZenith = offset_arrival_phZenith - \ + phZenith_offset(r_coordinates).reshape(photon_arrival_phZenith.shape) + + # Calculating the 1D histogram for the spatial counts and the 2D histogram for the arrival times and arrival angles + spatial_counts = np.histogram( + r_coordinates, r_bins, weights=photons_on_plane)[0] + time_counts = np.histogram2d(r_coordinates, offset_propagation_times, bins=( + r_bins, time_bins), weights=photons_on_plane)[0] + phZenith_counts = np.histogram2d(r_coordinates, offset_arrival_phZenith, bins=( + r_bins, phZenith_bins), weights=photons_on_plane)[0] + phAzi_counts = np.histogram2d(r_coordinates.flatten(), azimuth_angles.flatten(), bins=( + r_bins, phAzi_bins), weights=photons_on_plane.flatten())[0] + + spatial_counts = spatial_counts / ncounter_rbin + + # create branches for tshower + + eng[0] = energy + zen[0] = theta_D + azi[0] = azimdeg + startdist[0] = L_particle + X0[0]=X0_First + tshower.Fill() + + # vec_time = vector("TH1D")() + # vec_phZenith = vector("TH1D")() + + # define histograms + + num_wl_bin = len(wlbins) - 1 + num_dist_bin = len(r_bins) - 1 + num_time_bin = len(time_bins) - 1 + num_phZenith_bin = len(phZenith_bins) - 1 + num_phAzi_bin = len(phAzi_bins) - 1 + + + # histo_w = TH1D() + # histo_r = TH1D() + # histo_t_off = TH1D() + # histo_phZenith_off = TH1D() + + # histo_t=[TH1D()for i in range(num_dist_bin)] + # histo_phZenith=[TH1D()for i in range(num_dist_bin)] + + # define histograms + histo_w = TH1D("wl", "wavelength", num_wl_bin, wlbins) + histo_r = TH1D("r", "distance", num_dist_bin, r_bins) + histo_t_off = TH1D("t_off", "time offset", + num_dist_bin, r_bins) + histo_phZenith_off = TH1D("phZenith_off", "phZenith offset", + num_dist_bin, r_bins) + histo_t = [TH1D("t_dist_" + str(i), "time_dist_" + str(i), + num_time_bin, time_bins) for i in range(num_dist_bin)] + histo_phZenith = [TH1D("phZenith_dist_" + str(i), "phZenith_dist_" + str(i), + num_phZenith_bin, phZenith_bins) for i in range(num_dist_bin)] + histo_phAzi = [TH1D("phAzi_dist_" + str(i), "phAzi_dist_" + str(i), + num_phAzi_bin, phAzi_bins) for i in range(num_dist_bin)] + # fill histograms + for wl_bin, counts in enumerate(wlhist): + histo_w.SetBinContent(wl_bin + 1, counts) + for r_bin, counts in enumerate(spatial_counts): + histo_r.SetBinContent(r_bin + 1, counts) + histo_t_off.SetBinContent( + r_bin + 1, np.array(time_bounds)[:, 0][r_bin]) + histo_phZenith_off.SetBinContent( + r_bin + 1, np.array(zenith_bounds)[:, 0][r_bin]) + for t_bin, counts_t in enumerate(time_counts[r_bin]): + histo_t[r_bin].SetBinContent(t_bin + 1, counts_t) + for phZenith_bin, counts_phZenith in enumerate(phZenith_counts[r_bin]): + histo_phZenith[r_bin].SetBinContent(phZenith_bin + 1, counts_phZenith) + for phAzi_bin, counts_phAzi in enumerate(phAzi_counts[r_bin]): + histo_phAzi[r_bin].SetBinContent(phAzi_bin + 1, counts_phAzi) + + # set branch for histograms + tcher_ph.SetBranchAddress("wavelength", histo_w) + tcher_ph.SetBranchAddress("distance", histo_r) + tcher_ph.SetBranchAddress("phZenith_offset", histo_phZenith_off) + tcher_ph.SetBranchAddress("time_offset", histo_t_off) + vec_time.assign(histo_t) + vec_phZenith.assign(histo_phZenith) + vec_phAzi.assign(histo_phAzi) + + + tcher_ph.Fill() + # , vec_phZenith, vec_time + del histo_w, histo_r, histo_t, histo_phZenith, histo_phAzi, histo_t_off, histo_phZenith_off +tshower.Write("", TFile.kOverwrite) +tcher_ph.Write("", TFile.kOverwrite) +print('Output written to {}'.format(filename)) + + +file.Close() diff --git a/src/nuspacesim/chasm_to_offline_original.py b/src/nuspacesim/chasm_to_offline_original.py new file mode 100644 index 0000000..7288264 --- /dev/null +++ b/src/nuspacesim/chasm_to_offline_original.py @@ -0,0 +1,420 @@ +#!@Python_EXECUTABLE@ + +import sys +import uproot +import numpy as np +import argparse +from scipy.spatial.transform import Rotation as R +from scipy import interpolate +from ROOT import TFile, TTree, TH1D, vector +try: + import CHASM as ch +except ImportError as e: + print("Module CHASM could not be imported. It is most likely not installed. \n \ + It can be installed through pip: python -m pip install CHASM-NuSpacesim. \n The detailed error message is: \n",e) + sys.exit() + + +def get_vectors(zenith, azimuth, altitude, grid_angle, npoints): + + # get distance along axis corresponding to detector altitude + r = sim.ingredients['axis'].h_to_axis_R_LOC(altitude*1.e3, zenith) + grid_width = r*np.tan(np.radians(grid_angle)) + orig = [0, 0, r] + rcoord_count = np.linspace(0, grid_width, npoints) + + # Radial random angle + ang = 2*np.pi*np.random.rand(npoints) + xall = rcoord_count*np.cos(ang) + yall = rcoord_count*np.sin(ang) + zall = np.full_like(xall, r) + + # Square grid + # x = np.linspace(-grid_width, grid_width, n_side) + # y = np.linspace(-grid_width, grid_width, n_side) + # xx, yy = np.meshgrid(x,y) + # zz = np.full_like(xx, r) #convert altitude to m + + # vecs = np.vstack((xx.flatten(),yy.flatten(),zz.flatten())).T + vecs = np.vstack((xall.flatten(), yall.flatten(), zall.flatten())).T + + theta_rot_axis = np.array([0, 1, 0]) + theta_rotation = R.from_rotvec(theta_rot_axis * zenith) + + z_rot_axis = np.array([0, 0, 1]) + z_rotation = R.from_rotvec(z_rot_axis * np.pi/2) + + phi_rot_axis = np.array([0, 0, 1]) + phi_rot = R.from_rotvec(phi_rot_axis*azimuth) + + vecs = z_rotation.apply(vecs) + vecs = theta_rotation.apply(vecs) + + tel_vecs = phi_rot.apply(vecs) + orig = z_rotation.apply(orig) + orig = theta_rotation.apply(orig) + orig = phi_rot.apply(orig) + return tel_vecs, orig + + +def run_chasm(sim, orig): + sig = sim.run(mesh=False, att=True) + del sim + phot = np.array(sig.photons) + phot_lambda = phot.sum(axis=2).sum(axis=0) + phot_plane = phot.sum(axis=1) + counters_vec = sig.counters.vectors + dist_counter = np.sqrt(((counters_vec-orig)**2).sum(axis=-1))/1000 # km + source_vec = sig.source_points + to_orig = orig-source_vec + travel_vec = sig.counters.travel_vectors(source_vec) + dist_travel = np.sqrt((travel_vec**2).sum(axis=-1)) + plane_vec = travel_vec-to_orig + plane_vec = plane_vec.reshape([-1, 3]) + r_coord = np.sqrt((plane_vec**2).sum(axis=-1)) + + anglesin = 180/np.pi*np.arcsin(r_coord/dist_travel.flatten()) + times = np.array(sig.times) + times = times-np.min(times) + + r_coordinates = (r_coord/1000).flatten() + total_propagation_times = times.flatten() + photon_arrival_phZenith = anglesin.flatten() + photons_on_plane = phot_plane.flatten() + wlbins = sig.wavelengths + wlhist = phot_lambda.flatten() + return r_coordinates, total_propagation_times, photon_arrival_phZenith, photons_on_plane, wlbins, wlhist, dist_counter + + +def weighted_percentile(data, percents, weights=None): + ''' percents in units of 1% + weights specifies the frequency (count) of data. + ''' + if weights is None: + return np.percentile(data, percents) + ind = np.argsort(data) + d = data[ind] + w = weights[ind] + p = 1.*w.cumsum()/w.sum()*100 + y = np.interp(percents, p, d) + return y + + +def theta_D_from_theta_EE(det_alt, theta_EE): + """ + det_alt : float + detector altitude in km + theta_EE : float + earth emergence angle in rad + + Returns + ------- + float + detector viewing angle in degrees + """ + theta_EE = theta_EE * (np.pi / 180) + Re = 6371 + theta_D = (180 / np.pi) * np.arcsin((Re / (Re + det_alt)) + * np.sin(np.pi / 2 + theta_EE)) + + return theta_D + + +def h_to_axis_R_LOC(h, theta) -> np.ndarray: + earth_radius = 6371e3 + ground_level = 0 + """Return the length along the shower axis from the point of Earth + emergence to the height above the surface specified + + Parameters: + h: array of heights (m above ground level) + theta: polar angle of shower axis (radians) + + returns: r (m) (same size as h), an array of distances along the shower + axis_sp. + """ + cos_EM = np.cos(np.pi-theta) + R = earth_radius + ground_level + r_CoE = h + R # distance from the center of the earth to the specified height + rs = R*cos_EM + np.sqrt(R**2*cos_EM**2-R**2+r_CoE**2) + # rs -= rs[0] + # rs[0] = 1. + return rs + + +TH1D.AddDirectory(False) + +parser = argparse.ArgumentParser() +parser.add_argument("-nss", "--nuspacesim", dest="nssfile", + default=False, help="NuSpaceSim Conex-like input") +parser.add_argument("-p", "--profile", dest="profile", default=False, + help="Text file with slant depth (g/cm2) and N charged particles") +parser.add_argument("-gh", "--GaisserHillas", dest="gh", type=float, nargs=4, + default=False, help="Gaisser Hilas parameters as in -gh X N X0 Lambda") +parser.add_argument("-o", "--outputfile", dest="outfile", + default='chasm_to_offline.root', help="Output file name") +parser.add_argument("-e", "--showerenergy", dest="senergy", + type=float, default=False, help="Shower Energy in log10(E/eV)") +parser.add_argument("-z", "--zenith", dest="theta", type=float, default=False, + help="Zenith emergence angle in degrees (89deg is very horizontal, 0 deg is upward vertical)") +parser.add_argument("-a", "--azimuth", dest="phi", type=float, + default=0, help="Azimuth angle in degrees") +parser.add_argument("-hf", "--heightfirst", dest="hfirst", + type=float, default=0, help="Height of first interaction in km") +parser.add_argument("-alt", "--altitude", dest="alt", type=float, + default=33, help="Detector altitude in km (default 33km)") +parser.add_argument("-n", "--npoints", dest="n", type=int, default=2500, + help="Number of sampled points on detection plane (default 2500)") + +args = parser.parse_args() +if [args.nssfile, args.profile, args.gh].count(False) != 2: + print('Use only one input flag. Use either a nuspacesim input, a text profile input, or a Gaisser Hillas input') + exit() +if args.nssfile: + try: + with uproot.open(args.nssfile) as file: + print('Reading NuSpaceSim file') + + tshower = file["Shower"] + theader = file["Header"] + + zen90 = tshower["zenith"].array() + phi = np.array(np.radians(tshower["azimuth"].array())) + Hfirst = tshower["Hfirst"].array()/1000 + theta = np.array(np.radians(180-zen90)) + nEvents = np.size(zen90) + lgenergy = tshower["lgE"].array() + X = tshower["X"].array() + N = tshower["N"].array() + except: + print('Could not read the file {}'.format(args.nssfile)) + exit() + +else: + nEvents = 1 + if args.senergy: + if args.senergy < 30 and args.senergy > 9: + lgenergy = [args.senergy] + else: + print('Shower energy must be in log10(E/eV)') + exit() + else: + print('Introduce shower energy in log10(E/eV) using -e') + exit() + if args.theta: + if args.theta <= 90 and args.theta >= 0: + theta = [np.radians(args.theta)] + else: + print( + 'Shower must be upward. Zenith angle must be between 0(vertical) and 90(horizontal) degrees') + exit() + else: + print('Introduce zenith angle between 0(vertical) and 90 (horizontal) degrees using -z') + exit() + phi = [args.phi] + Hfirst = [args.hfirst] + if args.profile: + try: + X, N = np.loadtxt(args.profile) + X = [X] + N = [N] + print('Reading profile from text file') + except: + print('Could not read the file {}'.format(args.profile)) + exit() + else: + Xgh = args.gh[0] + Ngh = args.gh[1] + X0gh = args.gh[2] + Lgh = args.gh[3] + + +# add shower axis +zenith = theta +azimuth = phi + +# add grid of detectors +n_side = 30 # For square grid +npoints = args.n # For radial +detector_grid_alt = args.alt # km +grid_angle = 5 # degrees +file_name = args.outfile +filename = file_name +eng = np.zeros(1, dtype=float) +zen = np.zeros(1, dtype=float) +azi = np.zeros(1, dtype=float) +startdist = np.zeros(1, dtype=float) +file = TFile(filename, 'recreate') +tcher_ph = TTree("cherPhProp", "Cherenkov Photon Properties") +tshower = TTree("showerProp", "Shower Properties") +tshower.Branch("energy", eng, 'eng/D') +tshower.Branch("zenith", zen, 'zen/D') +tshower.Branch("azimuth", azi, 'azi/D') +tshower.Branch("startdist", startdist, 'startdist/D') + + +vec_time = vector("TH1D")() +vec_phZenith = vector("TH1D")() + +histo_w = TH1D() +histo_r = TH1D() +histo_t_off = TH1D() +histo_phZenith_off = TH1D() + +tcher_ph.Branch("wavelength", 'TH1D', histo_w) +tcher_ph.Branch("distance", 'TH1D', histo_r) +tcher_ph.Branch("time_offset", 'TH1D', histo_t_off) +tcher_ph.Branch("phZenith_offset", 'TH1D', histo_phZenith_off) +tcher_ph.Branch("time_dist", vec_time) +tcher_ph.Branch("phZenith_dist", vec_phZenith) + + +for i in range(nEvents): + sim = ch.ShowerSimulation() + sim.add(ch.UpwardAxis(zenith[i], azimuth[i], curved=True)) + tel_vecs, orig = get_vectors( + zenith[i], azimuth[i], detector_grid_alt, grid_angle, npoints) + if args.gh: + print('Generating Gaisser Hillas shower') + sim.add(ch.GHShower(Xgh, Ngh, X0gh, Lgh)) + else: + sim.add(ch.UserShower(np.array(X[i]), np.array(N[i]))) + sim.add(ch.SphericalCounters(tel_vecs, np.sqrt(1/np.pi))) + sim.add(ch.Yield(270, 1000, N_bins=100)) + print('Running CHASM simulation shower {}'.format(i)) + + r_coordinates, total_propagation_times, photon_arrival_phZenith, photons_on_plane, wlbins, wlhist, dist_counter = run_chasm( + sim, orig) + + energy = 10**lgenergy[i] + azimdeg = np.degrees(azimuth[i]) + theta_D = theta_D_from_theta_EE(detector_grid_alt, 90-np.degrees(theta[i])) + L_particle = Hfirst[i] # meters + + photons_on_plane[photons_on_plane < 0] = 0 + + bins, start, stop = 100, 0, 4 # For time + time_bins = np.logspace(float(start), float(stop), num=int(bins)) + bins, start, stop = 100, -3, 1 # For Zenith + angle_bins = np.logspace(float(start), float(stop), num=int(bins)) + spatial_bins = 100 + pct_cap = 100 # 30% recommended in easchersim, but 50 used + r_percentiles = weighted_percentile( + r_coordinates, pct_cap, weights=photons_on_plane) + r_bounds = weighted_percentile( + r_coordinates, np.array([1, 99.7]), weights=photons_on_plane) + r_bins = np.linspace(np.min(r_coordinates)-1e-7, + r_percentiles, spatial_bins) + ncounter_rbin = np.zeros(spatial_bins-1) + + time_bounds, zenith_bounds = [], [] + offset_propagation_times = total_propagation_times + offset_arrival_phZenith = photon_arrival_phZenith + for j in range(0, len(r_bins) - 1): + # Find coordinates which are inside each spatial bin + r_lower, r_upper = r_bins[j], r_bins[j + 1] + inbin = (dist_counter <= r_upper) & (dist_counter > r_lower) + ncounter_rbin[j] = int(np.sum(inbin)) + in_range = (r_coordinates >= r_lower) & (r_coordinates <= r_upper) + in_range_time, in_range_zenith, in_range_photons = total_propagation_times[ + in_range], photon_arrival_phZenith[in_range], photons_on_plane[in_range] + # Calculate the 1% and 99% weighted percentiles of the arrival time and arrival angle inside the bin + time_percentiles = weighted_percentile( + in_range_time, np.array([1, 99]), weights=in_range_photons) + zenith_percentiles = weighted_percentile( + in_range_zenith, np.array([1, 99]), weights=in_range_photons) + + time_bounds.append(time_percentiles) + zenith_bounds.append(zenith_percentiles) + + time_offset = interpolate.interp1d( + r_bins[:-1], np.array(time_bounds)[:, 0], kind='cubic', bounds_error=False, fill_value=0) + phZenith_offset = interpolate.interp1d( + r_bins[:-1], np.array(zenith_bounds)[:, 0], kind='cubic', bounds_error=False, fill_value=0) + + offset_propagation_times = total_propagation_times - \ + time_offset(r_coordinates).reshape(total_propagation_times.shape) + onephotonmask = (photons_on_plane >= 1) + offset_propagation_times -= np.min(offset_propagation_times[onephotonmask]) + + offset_arrival_phZenith = offset_arrival_phZenith - \ + phZenith_offset(r_coordinates).reshape(photon_arrival_phZenith.shape) + + # Calculating the 1D histogram for the spatial counts and the 2D histogram for the arrival times and arrival angles + spatial_counts = np.histogram( + r_coordinates, r_bins, weights=photons_on_plane)[0] + time_counts = np.histogram2d(r_coordinates, offset_propagation_times, bins=( + r_bins, time_bins), weights=photons_on_plane)[0] + angle_counts = np.histogram2d(r_coordinates, offset_arrival_phZenith, bins=( + r_bins, angle_bins), weights=photons_on_plane)[0] + spatial_counts = spatial_counts / ncounter_rbin + + # create branches for tshower + + eng[0] = energy + zen[0] = theta_D + azi[0] = azimdeg + startdist[0] = L_particle + + tshower.Fill() + + # vec_time = vector("TH1D")() + # vec_phZenith = vector("TH1D")() + + # define histograms + + num_wl_bin = len(wlbins) - 1 + num_dist_bin = len(r_bins) - 1 + num_time_bin = len(time_bins) - 1 + num_ang_bin = len(angle_bins) - 1 + + # histo_w = TH1D() + # histo_r = TH1D() + # histo_t_off = TH1D() + # histo_phZenith_off = TH1D() + + # histo_t=[TH1D()for i in range(num_dist_bin)] + # histo_phZenith=[TH1D()for i in range(num_dist_bin)] + + # define histograms + histo_w = TH1D("wl", "wavelength", num_wl_bin, wlbins) + histo_r = TH1D("r", "distance", num_dist_bin, r_bins) + histo_t_off = TH1D("t_off", "time offset", + num_dist_bin, r_bins) + histo_phZenith_off = TH1D("ang_off", "angle offset", + num_dist_bin, r_bins) + histo_t = [TH1D("t_dist_" + str(i), "time_dist_" + str(i), + num_time_bin, time_bins) for i in range(num_dist_bin)] + histo_phZenith = [TH1D("ang_dist_" + str(i), "angle_dist_" + str(i), + num_ang_bin, angle_bins) for i in range(num_dist_bin)] + # fill histograms + for wl_bin, counts in enumerate(wlhist): + histo_w.SetBinContent(wl_bin + 1, counts) + for r_bin, counts in enumerate(spatial_counts): + histo_r.SetBinContent(r_bin + 1, counts) + histo_t_off.SetBinContent( + r_bin + 1, np.array(time_bounds)[:, 0][r_bin]) + histo_phZenith_off.SetBinContent( + r_bin + 1, np.array(zenith_bounds)[:, 0][r_bin]) + for t_bin, counts_t in enumerate(time_counts[r_bin]): + histo_t[r_bin].SetBinContent(t_bin + 1, counts_t) + for ang_bin, counts_ang in enumerate(angle_counts[r_bin]): + histo_phZenith[r_bin].SetBinContent(ang_bin + 1, counts_ang) + + # set branch for histograms + tcher_ph.SetBranchAddress("wavelength", histo_w) + tcher_ph.SetBranchAddress("distance", histo_r) + tcher_ph.SetBranchAddress("phZenith_offset", histo_phZenith_off) + tcher_ph.SetBranchAddress("time_offset", histo_t_off) + vec_time.assign(histo_t) + vec_phZenith.assign(histo_phZenith) + + tcher_ph.Fill() + # , vec_phZenith, vec_time + del histo_w, histo_r, histo_t, histo_phZenith, histo_t_off, histo_phZenith_off +tshower.Write("", TFile.kOverwrite) +tcher_ph.Write("", TFile.kOverwrite) +print('Writing output to {}'.format(filename)) + + +file.Close() From be5a1e17ad3bc2ffa383bc53b51b9bf51fd0f2a6 Mon Sep 17 00:00:00 2001 From: jorge Date: Sat, 26 Apr 2025 20:53:43 -0600 Subject: [PATCH 2/6] Chasm implementation in cphotang.py --- sample_input_file.toml | 4 +- src/nuspacesim/compute.py | 56 ++++++++++++++++++- .../simulation/eas_optical/cphotang.py | 19 ++++++- 3 files changed, 74 insertions(+), 5 deletions(-) diff --git a/sample_input_file.toml b/sample_input_file.toml index 00f850d..3e0f79b 100644 --- a/sample_input_file.toml +++ b/sample_input_file.toml @@ -4,7 +4,7 @@ title = "NuSpaceSim" name = "Default Name" [detector.initial_position] -altitude = "525.0 km" +altitude = "500.0 km" latitude = "0.0 deg" longitude = "0.0 deg" @@ -30,7 +30,7 @@ gain = "1.8 dB" [simulation] mode = "Diffuse" -thrown_events = 100000 +thrown_events = 4 max_cherenkov_angle = "3.0000000000000004 deg" max_azimuth_angle = "360.0 deg" angle_from_limb = "7.0 deg" diff --git a/src/nuspacesim/compute.py b/src/nuspacesim/compute.py index 7826950..d7ef71f 100644 --- a/src/nuspacesim/compute.py +++ b/src/nuspacesim/compute.py @@ -62,7 +62,8 @@ from .simulation.eas_radio.radio import EASRadio from .simulation.eas_radio.radio_antenna import calculate_snr from .simulation.geometry.region_geometry import RegionGeom, RegionGeomToO - +from .chasm_testing import * +from .simulation.eas_optical.shower_properties import path_length_tau_atm # from .simulation.geometry.too import * from .simulation.spectra.spectra import Spectra from .simulation.taus.taus import Taus @@ -223,6 +224,59 @@ def add_meta(self, name: str, value: Any, comment: str): logv("Computing [green] Decay Altitudes.[/]") altDec, lenDec = eas.altDec(beta_tr, tauBeta, tauLorentz, store=sw) + #CHASM STUFF + #Need to figure out how to send azimuth, detcoords to cphotang function + #Need to check that coord system of detcoords is the same as azimTrSubN (should be correct) + zenith=np.pi/2-beta_tr + + from scipy.spatial.transform import Rotation as R + + #azimuth= + + + thetaTrSubV=geom.thetas() + phiTrSubV=geom.phis() + costhetaNSubV=geom.valid_costhetaNSubV() + det_azi=geom.valid_aziAngVSubN() + det_elev=geom.valid_elevAngVSubN() + thetaNSubV=np.arccos(costhetaNSubV) + Pvecx=np.sin(thetaTrSubV)*np.cos(phiTrSubV) + Pvecy=np.sin(thetaTrSubV)*np.sin(phiTrSubV) + Pvecz=np.cos(thetaTrSubV) + + P=np.vstack((Pvecx,Pvecy,Pvecz)).T + n_up_rot_axis = np.array([0, 1, 0]) + rot_vectors = n_up_rot_axis * thetaNSubV[:, np.newaxis] # Shape (N, 3) + theta_rotation = R.from_rotvec(rot_vectors) + + z_rot_axis = np.array([0, 0, 1]) # z-axis + z_rot_vectors = z_rot_axis * (det_azi)[:, np.newaxis] # Shape (N, 3) + z_rotation = R.from_rotvec(z_rot_vectors) + + Pn = theta_rotation.apply(P) + tr_n=z_rotation.apply(Pn) + azimuth=np.arctan2(tr_n[:,1],tr_n[:,0]) + + + path_len=geom.pathLens() + #path_len2=path_length_tau_atm(altitude,geom.valid_elevAngVSubN()) This is the same as path_len + x = path_len * np.cos(det_elev) * np.cos(det_azi) # East + y = path_len * np.cos(det_elev) * np.sin(det_azi) # North + z = path_len * np.sin(det_elev) # Up + detcoords=np.vstack((x,y,z)).T + print(path_len,np.degrees(beta_tr)) + #print('Zen, zenith, beta,detazi,detzen',np.degrees(zen),np.degrees(zenith),np.degrees(beta_tr),np.degrees(det_azi),np.degrees(det_zen)) + #print(det_coord) + for i in range(0, len(beta_tr)): + sim = ch.ShowerSimulation() + sim.add(ch.UpwardAxis(zenith[i], azimuth[i], curved=True)) + sim.add(ch.UserShower(np.array(X[i]), np.array(N[i]))) + sim.add(ch.SphericalCounters(detcoords, np.sqrt(1/np.pi))) + sim.add(ch.Yield(270, 1000, N_bins=100)) + r_coordinates, total_propagation_times, photon_arrival_phZenith, photons_on_plane, wlbins, wlhist, dist_counter = run_chasm(sim, orig) + + + # if config.detector.method == "Optical" or config.detector.method == "Both": if config.detector.optical.enable: logv("Computing [green] EAS Optical Cherenkov light.[/]") diff --git a/src/nuspacesim/simulation/eas_optical/cphotang.py b/src/nuspacesim/simulation/eas_optical/cphotang.py index 1c31b49..e4a8de5 100644 --- a/src/nuspacesim/simulation/eas_optical/cphotang.py +++ b/src/nuspacesim/simulation/eas_optical/cphotang.py @@ -47,7 +47,8 @@ from numpy.polynomial import Polynomial from .detector_geometry import distance_to_detector - +from .chasm_geom import detector_coordinates_and_tr_azimuth +import CHASM as ch # Wrapped in try-catch block as a hack to enable sphinx documentation to be generated # on ReadTheDocs without pre-compiling. try: @@ -591,7 +592,21 @@ def cherenkov_area(self, AveCangI, DistStep, izRNmax): CherArea *= DistStep[izRNmax] CherArea = self.pi * np.power(CherArea, 2, dtype=self.dtype) return CherArea - + + def chasm(betaE, gramsum, RN): + detcoord, azimuth_tr= detector_coordinates_and_tr_azimuth() + sim = ch.ShowerSimulation() + sim.add(ch.UpwardAxis(np.pi/2-betaE, azimuth_tr, curved=True)) + sim.add(ch.UserShower(np.array(gramsum), np.array(RN))) + sim.add(ch.SphericalCounters(detcoord, np.sqrt(1/np.pi))) + sim.add(ch.Yield(270, 1000, N_bins=100)) + sig = sim.run(mesh=False, att=True) + del sim + ch_photons=np.array(sig.photons) + times = np.array(sig.times) + costheta=np.array(sig.cos_theta) + return ch_photons, times, costheta + def run(self, betaE, alt, Eshow100PeV, lat, long, cloudf=None): """Main body of simulation code.""" From 0305f1a891df24acc89b7e23545ce8b941fb0900 Mon Sep 17 00:00:00 2001 From: jorge Date: Sat, 26 Apr 2025 21:30:42 -0600 Subject: [PATCH 3/6] Chasm Implementation continuation --- src/nuspacesim/compute.py | 11 ++++++-- .../simulation/eas_optical/cphotang.py | 27 ++++++++++--------- src/nuspacesim/simulation/eas_optical/eas.py | 4 +++ 3 files changed, 27 insertions(+), 15 deletions(-) diff --git a/src/nuspacesim/compute.py b/src/nuspacesim/compute.py index d7ef71f..4b73b5d 100644 --- a/src/nuspacesim/compute.py +++ b/src/nuspacesim/compute.py @@ -62,11 +62,12 @@ from .simulation.eas_radio.radio import EASRadio from .simulation.eas_radio.radio_antenna import calculate_snr from .simulation.geometry.region_geometry import RegionGeom, RegionGeomToO -from .chasm_testing import * from .simulation.eas_optical.shower_properties import path_length_tau_atm # from .simulation.geometry.too import * from .simulation.spectra.spectra import Spectra from .simulation.taus.taus import Taus +from .chasm_geom import detector_coordinates_and_tr_azimuth + __all__ = ["compute"] @@ -224,6 +225,7 @@ def add_meta(self, name: str, value: Any, comment: str): logv("Computing [green] Decay Altitudes.[/]") altDec, lenDec = eas.altDec(beta_tr, tauBeta, tauLorentz, store=sw) + """ #CHASM STUFF #Need to figure out how to send azimuth, detcoords to cphotang function #Need to check that coord system of detcoords is the same as azimTrSubN (should be correct) @@ -274,9 +276,12 @@ def add_meta(self, name: str, value: Any, comment: str): sim.add(ch.SphericalCounters(detcoords, np.sqrt(1/np.pi))) sim.add(ch.Yield(270, 1000, N_bins=100)) r_coordinates, total_propagation_times, photon_arrival_phZenith, photons_on_plane, wlbins, wlhist, dist_counter = run_chasm(sim, orig) + """ + detcoords, azimuth = detector_coordinates_and_tr_azimuth( + geom.thetas(), geom.phis(), geom.valid_costhetaNSubV(), geom.valid_aziAngVSubN(), geom.valid_elevAngVSubN(), geom.pathLens() - + ) # if config.detector.method == "Optical" or config.detector.method == "Both": if config.detector.optical.enable: logv("Computing [green] EAS Optical Cherenkov light.[/]") @@ -287,6 +292,8 @@ def add_meta(self, name: str, value: Any, comment: str): showerEnergy, init_lat, init_long, + detcoords, + azimuth, cloudf=cloud, store=sw, plot=to_plot, diff --git a/src/nuspacesim/simulation/eas_optical/cphotang.py b/src/nuspacesim/simulation/eas_optical/cphotang.py index e4a8de5..631554b 100644 --- a/src/nuspacesim/simulation/eas_optical/cphotang.py +++ b/src/nuspacesim/simulation/eas_optical/cphotang.py @@ -47,7 +47,6 @@ from numpy.polynomial import Polynomial from .detector_geometry import distance_to_detector -from .chasm_geom import detector_coordinates_and_tr_azimuth import CHASM as ch # Wrapped in try-catch block as a hack to enable sphinx documentation to be generated # on ReadTheDocs without pre-compiling. @@ -546,8 +545,9 @@ def valid_arrays(self, zsave, delgram, gramsum, gramz, ZonZ, ThetPrpA, Eshow): s = s[mask] RN = RN[mask] e2hill = e2hill[mask] + gramsum = gramsum[mask] - return zs, delgram, ZonZ, ThetPrpA, AirN, s, RN, e2hill + return zs, delgram, ZonZ, ThetPrpA, AirN, s, RN, e2hill, gramsum def e0(self, shape, s): """Hillas Energy Paramaterization. @@ -593,21 +593,20 @@ def cherenkov_area(self, AveCangI, DistStep, izRNmax): CherArea = self.pi * np.power(CherArea, 2, dtype=self.dtype) return CherArea - def chasm(betaE, gramsum, RN): - detcoord, azimuth_tr= detector_coordinates_and_tr_azimuth() + def chasm(self, betaE, gramsum, RN, detcoord, azimuth_tr): sim = ch.ShowerSimulation() sim.add(ch.UpwardAxis(np.pi/2-betaE, azimuth_tr, curved=True)) sim.add(ch.UserShower(np.array(gramsum), np.array(RN))) sim.add(ch.SphericalCounters(detcoord, np.sqrt(1/np.pi))) sim.add(ch.Yield(270, 1000, N_bins=100)) sig = sim.run(mesh=False, att=True) - del sim ch_photons=np.array(sig.photons) - times = np.array(sig.times) - costheta=np.array(sig.cos_theta) - return ch_photons, times, costheta + ch_times = np.array(sig.times) + ch_costheta=np.array(sig.cos_theta) #cosine of the angle between the z-axis (zenith of spot on ground) and the vector from the axis to the detector + + return ch_photons, ch_times, ch_costheta - def run(self, betaE, alt, Eshow100PeV, lat, long, cloudf=None): + def run(self, betaE, alt, Eshow100PeV, lat, long, detcoords, azimuth, cloudf=None): """Main body of simulation code.""" # Should we just skip these with a mask in valid_arrays? @@ -623,10 +622,10 @@ def run(self, betaE, alt, Eshow100PeV, lat, long, cloudf=None): # Shower # - zs, delgram, ZonZ, ThetPrpA, AirN, s, RN, e2hill = self.valid_arrays( + zs, delgram, ZonZ, ThetPrpA, AirN, s, RN, e2hill, gramsum = self.valid_arrays( *self.slant_depth(alt, sinThetView), Eshow ) - + ch_photons, ch_times, ch_costheta = self.chasm(betaE, gramsum, RN, detcoords, azimuth) # Cloud top height cloud_top_height = cloudf(lat, long) if cloudf else -np.inf @@ -684,10 +683,11 @@ def run(self, betaE, alt, Eshow100PeV, lat, long, cloudf=None): photonDen *= altitude_scaling Cang = np.degrees(AveCangI + CangsigI) + print('Cang',Cang.shape, Cang) return photonDen, Cang - def __call__(self, betaE, alt, Eshow100PeV, init_lat, init_long, cloudf=None): + def __call__(self, betaE, alt, Eshow100PeV, init_lat, init_long, detcoords, azimuth, cloudf=None): """ Iterate over the list of events and return the result as pair of numpy arrays. @@ -707,5 +707,6 @@ def __call__(self, betaE, alt, Eshow100PeV, init_lat, init_long, cloudf=None): zip(betaE, alt, Eshow100PeV, init_lat, init_long), partition_size=100 ) with ProgressBar(): - Dphots, Cang = zip(*b.map(lambda x: self.run(*x, cloudf)).compute()) + Dphots, Cang = zip(*b.map(lambda x: self.run(*x, detcoords, azimuth, cloudf)).compute()) + print('Cang',Cang.shape, Cang) return np.asarray(Dphots), np.array(Cang) diff --git a/src/nuspacesim/simulation/eas_optical/eas.py b/src/nuspacesim/simulation/eas_optical/eas.py index 2033fc9..06fdf3a 100644 --- a/src/nuspacesim/simulation/eas_optical/eas.py +++ b/src/nuspacesim/simulation/eas_optical/eas.py @@ -86,6 +86,8 @@ def __call__( showerEnergy, init_lat, init_long, + detcoords, + azimuth, *args, cloudf=None, **kwargs @@ -110,6 +112,8 @@ def __call__( showerEnergy[mask], init_lat[mask], init_long[mask], + detcoords, + azimuth, cloudf, ) From 7abbdbb27450d7b24fd22b3aad89c04c23fb37a7 Mon Sep 17 00:00:00 2001 From: jorge Date: Mon, 12 May 2025 16:31:31 -0600 Subject: [PATCH 4/6] Progress on CHASM implementation checks --- sample_input_file.toml | 6 +- src/nuspacesim/compute.py | 78 +++++++++++++++++-- .../simulation/eas_optical/cphotang.py | 27 ++++--- src/nuspacesim/simulation/eas_optical/eas.py | 10 +-- 4 files changed, 95 insertions(+), 26 deletions(-) diff --git a/sample_input_file.toml b/sample_input_file.toml index 3e0f79b..8c395a9 100644 --- a/sample_input_file.toml +++ b/sample_input_file.toml @@ -4,7 +4,7 @@ title = "NuSpaceSim" name = "Default Name" [detector.initial_position] -altitude = "500.0 km" +altitude = "33.0 km" latitude = "0.0 deg" longitude = "0.0 deg" @@ -30,7 +30,7 @@ gain = "1.8 dB" [simulation] mode = "Diffuse" -thrown_events = 4 +thrown_events = 100000 max_cherenkov_angle = "3.0000000000000004 deg" max_azimuth_angle = "360.0 deg" angle_from_limb = "7.0 deg" @@ -48,7 +48,7 @@ table_version = "3" [simulation.spectrum] id = "monospectrum" -log_nu_energy = 8.0 +log_nu_energy = 12 [simulation.cloud_model] id = "no_cloud" diff --git a/src/nuspacesim/compute.py b/src/nuspacesim/compute.py index 4b73b5d..cc23537 100644 --- a/src/nuspacesim/compute.py +++ b/src/nuspacesim/compute.py @@ -51,6 +51,7 @@ from typing import Any, Iterable import numpy as np +import awkward as ak from astropy.table import Table as AstropyTable from numpy.typing import ArrayLike from rich.console import Console @@ -66,7 +67,7 @@ # from .simulation.geometry.too import * from .simulation.spectra.spectra import Spectra from .simulation.taus.taus import Taus -from .chasm_geom import detector_coordinates_and_tr_azimuth +from .chasm_geom import detector_coordinates_and_tr_azimuth, point_to_line_distances, angle_at_decay __all__ = ["compute"] @@ -279,25 +280,90 @@ def add_meta(self, name: str, value: Any, comment: str): """ detcoords, azimuth = detector_coordinates_and_tr_azimuth( - geom.thetas(), geom.phis(), geom.valid_costhetaNSubV(), geom.valid_aziAngVSubN(), geom.valid_elevAngVSubN(), geom.pathLens() - + geom.thetas(), geom.phis()%(2*np.pi), geom.valid_costhetaNSubV(), geom.valid_aziAngVSubN()%(2*np.pi), geom.valid_elevAngVSubN(), geom.pathLens() ) + altDec1km= np.full_like(altDec, 1.05) # if config.detector.method == "Optical" or config.detector.method == "Both": if config.detector.optical.enable: logv("Computing [green] EAS Optical Cherenkov light.[/]") - numPEs, costhetaChEff = eas( + numPEs, costhetaChEff, ch_photons, ch_times, ch_costheta = eas( beta_tr, - altDec, + altDec1km, showerEnergy, init_lat, init_long, detcoords, azimuth, cloudf=cloud, - store=sw, + #store=sw, plot=to_plot, ) + import matplotlib.pyplot as plt + + mask = (altDec1km < 1) | (altDec1km > 1.1) | (showerEnergy*2 < 100) | (showerEnergy*2 > 200) + mask= ~mask + distance= point_to_line_distances(detcoords[mask], beta_tr[mask], azimuth[mask])/1000 #km + angleatdecay=angle_at_decay(detcoords[mask],altDec1km[mask], beta_tr[mask], azimuth[mask]) + logtauenergy = np.log10(tauEnergy)+9 + print(logtauenergy[mask]) + print('HOLA',ch_photons, ch_times, ak.num(ch_photons,axis=1), ak.num(ch_times)) + ch_phot_integrated=ak.sum(ch_photons, axis=1) + ch_theta=np.arccos(ch_costheta) + #print(len(ch_photons), len(ch_photons[0]), len(ch_photons[0][0])) + #print(len(ch_phot_integrated), len(ch_phot_integrated[0]), ak.num(ch_phot_integrated)) + """for i in range(len(ch_photons)): + maskphotons=(ch_phot_integrated[i]>0.1) + if np.sum(maskphotons) == 0: + continue + + plt.figure() + plt.plot(np.degrees(ch_theta[i][maskphotons]),ch_phot_integrated[i][maskphotons],'.',label=f'Tau Energy={logtauenergy[mask][i]:.2f} \n Beta_tr ={np.degrees(beta_tr)[mask][i]:.2f}') + plt.legend() + plt.grid() + plt.xlabel('Theta (degrees)') + plt.xscale('log') + plt.yscale('log') + plt.ylabel('Integrated Photons at detector') + plt.title('Integrated Photons vs Theta') + plt.savefig(f'IntegratedPhotons_{i}.png') + print('TOTAL CHERENKOV PHOTONS AT DETECTOR: ', ak.sum(ch_phot_integrated[i])) + print(ch_phot_integrated[i])""" + totphot=ak.sum(ch_phot_integrated,axis=1) + data = np.column_stack((totphot, distance, np.degrees(angleatdecay))) + #n = data.shape[0] # Get the number of rows + np.savetxt('data100000.txt', data, delimiter=',', header='totphot,distance,angleatdecay_degrees', comments='', fmt='%f') + + plt.figure() + plt.plot(distance, totphot,'.') + plt.xlabel('Distance to shower axis (km)') + plt.ylabel('Total Cherenkov Photons at detector for that shower') + plt.yscale('log') + plt.grid() + plt.title('Total photons vs distance to shower axis, tau energy 100-200 EeV') + plt.savefig('photonsvsdistance.png') + + plt.figure() + plt.plot(np.degrees(angleatdecay), totphot,'.') + plt.xlabel('Angle between shower axis and detector at decay point (degrees)') + plt.ylabel('Total Cherenkov Photons at detector for that shower') + plt.yscale('log') + plt.grid() + plt.title('Total photons vs angle to shower axis, tau energy 100-200 EeV') + plt.savefig('photonsvsangle.png') + + plt.figure() + plt.plot(np.degrees(angleatdecay), totphot,'.') + plt.xlabel('Angle between shower axis and detector at decay point (degrees)') + plt.ylabel('Total Cherenkov Photons at detector for that shower') + plt.gca().set_ylim(bottom=50) + plt.xlim([0,6]) + plt.yscale('log') + plt.grid() + plt.title('Total photons vs angle to shower axis, tau energy 100-200 EeV') + plt.savefig('photonsvsangle_lims.png') + #NEXT CHECK ANGLE BETWEEN SHOWER AXIS AND DETECTOR AT DECAY POINT OR AT XMAX + # Save to CSV with header logv("Computing [green] Optical Monte Carlo Integral.[/]") mcint, mcintgeo, passEV, mcunc = geom.mcintegral( diff --git a/src/nuspacesim/simulation/eas_optical/cphotang.py b/src/nuspacesim/simulation/eas_optical/cphotang.py index 631554b..ca1fb36 100644 --- a/src/nuspacesim/simulation/eas_optical/cphotang.py +++ b/src/nuspacesim/simulation/eas_optical/cphotang.py @@ -43,10 +43,14 @@ import dask.bag as db import numpy as np +import awkward as ak from dask.diagnostics import ProgressBar from numpy.polynomial import Polynomial from .detector_geometry import distance_to_detector +#from .....chasm.src.CHASM.simulation import ShowerSimulation +#from .....chasm.src.CHASM.user_inputs import UpwardAxis, UserShower, SphericalCounters + import CHASM as ch # Wrapped in try-catch block as a hack to enable sphinx documentation to be generated # on ReadTheDocs without pre-compiling. @@ -597,14 +601,14 @@ def chasm(self, betaE, gramsum, RN, detcoord, azimuth_tr): sim = ch.ShowerSimulation() sim.add(ch.UpwardAxis(np.pi/2-betaE, azimuth_tr, curved=True)) sim.add(ch.UserShower(np.array(gramsum), np.array(RN))) - sim.add(ch.SphericalCounters(detcoord, np.sqrt(1/np.pi))) + sim.add(ch.SphericalCounters(detcoord[np.newaxis, :], np.sqrt(1/np.pi))) sim.add(ch.Yield(270, 1000, N_bins=100)) sig = sim.run(mesh=False, att=True) - ch_photons=np.array(sig.photons) - ch_times = np.array(sig.times) - ch_costheta=np.array(sig.cos_theta) #cosine of the angle between the z-axis (zenith of spot on ground) and the vector from the axis to the detector + ch_photons=np.array(sig.photons.sum(axis=0)) + ch_times = np.array(sig.times.sum(axis=0)) + ch_costheta=np.array(sig.cos_theta.sum(axis=0).sum(axis=0)) #cosine of the angle between the z-axis (zenith of spot on ground) and the vector from the axis to the detector - return ch_photons, ch_times, ch_costheta + return np.array(ch_photons), np.array(ch_times), np.array(ch_costheta) def run(self, betaE, alt, Eshow100PeV, lat, long, detcoords, azimuth, cloudf=None): """Main body of simulation code.""" @@ -626,6 +630,7 @@ def run(self, betaE, alt, Eshow100PeV, lat, long, detcoords, azimuth, cloudf=Non *self.slant_depth(alt, sinThetView), Eshow ) ch_photons, ch_times, ch_costheta = self.chasm(betaE, gramsum, RN, detcoords, azimuth) + # Cloud top height cloud_top_height = cloudf(lat, long) if cloudf else -np.inf @@ -683,9 +688,7 @@ def run(self, betaE, alt, Eshow100PeV, lat, long, detcoords, azimuth, cloudf=Non photonDen *= altitude_scaling Cang = np.degrees(AveCangI + CangsigI) - print('Cang',Cang.shape, Cang) - - return photonDen, Cang + return photonDen, Cang, ch_photons, ch_times, ch_costheta def __call__(self, betaE, alt, Eshow100PeV, init_lat, init_long, detcoords, azimuth, cloudf=None): """ @@ -704,9 +707,9 @@ def __call__(self, betaE, alt, Eshow100PeV, init_lat, init_long, detcoords, azim ####################### b = db.from_sequence( - zip(betaE, alt, Eshow100PeV, init_lat, init_long), partition_size=100 + zip(betaE, alt, Eshow100PeV, init_lat, init_long, detcoords, azimuth), partition_size=100 ) with ProgressBar(): - Dphots, Cang = zip(*b.map(lambda x: self.run(*x, detcoords, azimuth, cloudf)).compute()) - print('Cang',Cang.shape, Cang) - return np.asarray(Dphots), np.array(Cang) + Dphots, Cang, ch_photons, ch_times, ch_costheta = zip(*b.map(lambda x: self.run(*x, cloudf)).compute()) + + return np.asarray(Dphots), np.array(Cang), ak.Array(ch_photons), ak.Array(ch_times), ak.Array(ch_costheta) diff --git a/src/nuspacesim/simulation/eas_optical/eas.py b/src/nuspacesim/simulation/eas_optical/eas.py index 06fdf3a..4d36cdc 100644 --- a/src/nuspacesim/simulation/eas_optical/eas.py +++ b/src/nuspacesim/simulation/eas_optical/eas.py @@ -98,7 +98,7 @@ def __call__( # Mask out-of-bounds events. Do not pass to CphotAng. Instead use # Default values for dphots and thetaCh - mask = (altDec < 0.0) | (altDec > 20.0) + mask = (altDec < 1) | (altDec > 1.1) | (showerEnergy*2 < 100) | (showerEnergy*2 > 200) mask = ~mask # phots and theta arrays with default 0 and 1.5 values. @@ -106,14 +106,14 @@ def __call__( thetaCh100PeV = np.full_like(beta, 1.5) # Run CphotAng on in-bounds events - dphots[mask], thetaCh100PeV[mask] = self.CphotAng( + dphots[mask], thetaCh100PeV[mask], ch_photons, ch_times, ch_costheta = self.CphotAng( beta[mask], altDec[mask], showerEnergy[mask], init_lat[mask], init_long[mask], - detcoords, - azimuth, + detcoords[mask], + azimuth[mask], cloudf, ) @@ -136,7 +136,7 @@ def __call__( costhetaChEff = np.cos(np.radians(thetaChEff)) - return numPEs, costhetaChEff + return numPEs, costhetaChEff, ch_photons, ch_times, ch_costheta def show_plot(sim, simclass, plot): From 4976d0970ac50d42596cdcfc240ade6089907074 Mon Sep 17 00:00:00 2001 From: jorge Date: Tue, 20 May 2025 12:26:16 -0600 Subject: [PATCH 5/6] Detector coordinate and azimuth geometry calculation for Chasm --- src/nuspacesim/chasm_geom.py | 92 ++++++++++++++++++++++++++++++++++++ 1 file changed, 92 insertions(+) create mode 100644 src/nuspacesim/chasm_geom.py diff --git a/src/nuspacesim/chasm_geom.py b/src/nuspacesim/chasm_geom.py new file mode 100644 index 0000000..b3feefe --- /dev/null +++ b/src/nuspacesim/chasm_geom.py @@ -0,0 +1,92 @@ +import numpy as np +from scipy.spatial.transform import Rotation as R + +def detector_coordinates_and_tr_azimuth(thetaTrSubV,phiTrSubV,costhetaNSubV,det_azi,det_elev, path_len): + + + thetaNSubV=np.arccos(costhetaNSubV) + Pvecx=np.sin(thetaTrSubV)*np.cos(phiTrSubV) + Pvecy=np.sin(thetaTrSubV)*np.sin(phiTrSubV) + Pvecz=np.cos(thetaTrSubV) + + P=np.vstack((Pvecx,Pvecy,Pvecz)).T + n_up_rot_axis = np.array([0, 1, 0]) + rot_vectors = n_up_rot_axis * thetaNSubV[:, np.newaxis] # Shape (N, 3) + theta_rotation = R.from_rotvec(rot_vectors) + + z_rot_axis = np.array([0, 0, 1]) # z-axis + z_rot_vectors = z_rot_axis * (det_azi)[:, np.newaxis] # Shape (N, 3) + z_rotation = R.from_rotvec(z_rot_vectors) + + Pn = theta_rotation.apply(P) + tr_n=z_rotation.apply(Pn) + azimuth=np.arctan2(tr_n[:,1],tr_n[:,0])%(2*np.pi) #change to between 0 and 2pi. Does it matter that azim definition doesn't match between NSS and CHASM?? + #path_len2=path_length_tau_atm(altitude,geom.valid_elevAngVSubN()) This is the same as path_len + x = path_len * np.cos(det_elev) * np.cos(det_azi) # East + y = path_len * np.cos(det_elev) * np.sin(det_azi) # North + z = path_len * np.sin(det_elev) # Up + detcoords=1000*np.vstack((x,y,z)).T + return detcoords, azimuth + +def point_to_line_distances(detcoords, betaE, azimuth): + """ + Calculate the shortest distance from each point in detcoords to the corresponding line + defined by betaE (emergence angle) and azimuth. + + Parameters: + - detcoords: numpy array of shape (N, 3), points in 3D space (x, y, z). + - betaE: numpy array of shape (N,), emergence angles in radians. + - azimuth: numpy array of shape (N,), azimuth angles in radians. + + Returns: + - distances: numpy array of shape (N,), distances from each point to its corresponding line. + """ + # Ensure inputs are numpy arrays + detcoords = np.asarray(detcoords) + betaE = np.asarray(betaE) + azimuth = np.asarray(azimuth) + + # Check shapes + N = detcoords.shape[0] + if detcoords.shape != (N, 3) or betaE.shape != (N,) or azimuth.shape != (N,): + raise ValueError("detcoords must be (N, 3), betaE and azimuth must be (N,)") + + # Compute direction vectors of the lines + # Direction vector: (cos(betaE) * cos(azimuth), cos(betaE) * sin(azimuth), sin(betaE)) + dx = np.cos(betaE) * np.cos(azimuth) + dy = np.cos(betaE) * np.sin(azimuth) + dz = np.sin(betaE) + + # Stack into direction vectors of shape (N, 3) + direction = np.stack([dx, dy, dz], axis=-1) # Shape: (N, 3) + + # Normalize the direction vectors (though not strictly necessary since norm will be computed) + # norm_direction = np.linalg.norm(direction, axis=1, keepdims=True) + # direction = direction / norm_direction + + # Compute the cross product: detcoords × direction + # detcoords: (N, 3), direction: (N, 3) + cross_product = np.cross(detcoords, direction) # Shape: (N, 3) + + # Compute the magnitude of the cross product + cross_magnitude = np.linalg.norm(cross_product, axis=1) # Shape: (N,) + + # Compute the magnitude of the direction vector + direction_magnitude = np.linalg.norm(direction, axis=1) # Shape: (N,) + + # Compute the distance + distances = cross_magnitude / direction_magnitude + + #print('DISTANCES',distances, np.shape(distances)) + + return distances + +def angle_at_decay(detcoord,altDec,beta_tr,azimuth): + xdec = 1000*altDec*np.cos(beta_tr) * np.cos(azimuth) + ydec = 1000*altDec*np.cos(beta_tr) * np.sin(azimuth) + zdec = 1000*altDec*np.sin(beta_tr) + decaycoord=np.vstack((xdec,ydec,zdec)).T + decaytodetector=detcoord-decaycoord + cosangle = np.sum(decaytodetector * decaycoord, axis=1) / (np.linalg.norm(decaytodetector, axis=1) * np.linalg.norm(decaycoord, axis=1)) + #print('COSANGLE',cosangle, np.shape(cosangle)) + return np.arccos(cosangle) \ No newline at end of file From eee209052d196431ac7b025ce695800476757df4 Mon Sep 17 00:00:00 2001 From: jorge Date: Tue, 20 May 2025 12:26:59 -0600 Subject: [PATCH 6/6] Chasm implementation in cphotang and compute --- sample_input_file.toml | 4 +- src/nuspacesim/compute.py | 91 ++++--------------- .../simulation/eas_optical/cphotang.py | 4 +- src/nuspacesim/simulation/eas_optical/eas.py | 2 +- 4 files changed, 23 insertions(+), 78 deletions(-) diff --git a/sample_input_file.toml b/sample_input_file.toml index 8c395a9..ba9551b 100644 --- a/sample_input_file.toml +++ b/sample_input_file.toml @@ -30,7 +30,7 @@ gain = "1.8 dB" [simulation] mode = "Diffuse" -thrown_events = 100000 +thrown_events = 10 max_cherenkov_angle = "3.0000000000000004 deg" max_azimuth_angle = "360.0 deg" angle_from_limb = "7.0 deg" @@ -48,7 +48,7 @@ table_version = "3" [simulation.spectrum] id = "monospectrum" -log_nu_energy = 12 +log_nu_energy = 10 [simulation.cloud_model] id = "no_cloud" diff --git a/src/nuspacesim/compute.py b/src/nuspacesim/compute.py index cc23537..81bc31e 100644 --- a/src/nuspacesim/compute.py +++ b/src/nuspacesim/compute.py @@ -226,70 +226,16 @@ def add_meta(self, name: str, value: Any, comment: str): logv("Computing [green] Decay Altitudes.[/]") altDec, lenDec = eas.altDec(beta_tr, tauBeta, tauLorentz, store=sw) - """ - #CHASM STUFF - #Need to figure out how to send azimuth, detcoords to cphotang function - #Need to check that coord system of detcoords is the same as azimTrSubN (should be correct) - zenith=np.pi/2-beta_tr - - from scipy.spatial.transform import Rotation as R - - #azimuth= - - - thetaTrSubV=geom.thetas() - phiTrSubV=geom.phis() - costhetaNSubV=geom.valid_costhetaNSubV() - det_azi=geom.valid_aziAngVSubN() - det_elev=geom.valid_elevAngVSubN() - thetaNSubV=np.arccos(costhetaNSubV) - Pvecx=np.sin(thetaTrSubV)*np.cos(phiTrSubV) - Pvecy=np.sin(thetaTrSubV)*np.sin(phiTrSubV) - Pvecz=np.cos(thetaTrSubV) - - P=np.vstack((Pvecx,Pvecy,Pvecz)).T - n_up_rot_axis = np.array([0, 1, 0]) - rot_vectors = n_up_rot_axis * thetaNSubV[:, np.newaxis] # Shape (N, 3) - theta_rotation = R.from_rotvec(rot_vectors) - - z_rot_axis = np.array([0, 0, 1]) # z-axis - z_rot_vectors = z_rot_axis * (det_azi)[:, np.newaxis] # Shape (N, 3) - z_rotation = R.from_rotvec(z_rot_vectors) - - Pn = theta_rotation.apply(P) - tr_n=z_rotation.apply(Pn) - azimuth=np.arctan2(tr_n[:,1],tr_n[:,0]) - - - path_len=geom.pathLens() - #path_len2=path_length_tau_atm(altitude,geom.valid_elevAngVSubN()) This is the same as path_len - x = path_len * np.cos(det_elev) * np.cos(det_azi) # East - y = path_len * np.cos(det_elev) * np.sin(det_azi) # North - z = path_len * np.sin(det_elev) # Up - detcoords=np.vstack((x,y,z)).T - print(path_len,np.degrees(beta_tr)) - #print('Zen, zenith, beta,detazi,detzen',np.degrees(zen),np.degrees(zenith),np.degrees(beta_tr),np.degrees(det_azi),np.degrees(det_zen)) - #print(det_coord) - for i in range(0, len(beta_tr)): - sim = ch.ShowerSimulation() - sim.add(ch.UpwardAxis(zenith[i], azimuth[i], curved=True)) - sim.add(ch.UserShower(np.array(X[i]), np.array(N[i]))) - sim.add(ch.SphericalCounters(detcoords, np.sqrt(1/np.pi))) - sim.add(ch.Yield(270, 1000, N_bins=100)) - r_coordinates, total_propagation_times, photon_arrival_phZenith, photons_on_plane, wlbins, wlhist, dist_counter = run_chasm(sim, orig) - """ - detcoords, azimuth = detector_coordinates_and_tr_azimuth( geom.thetas(), geom.phis()%(2*np.pi), geom.valid_costhetaNSubV(), geom.valid_aziAngVSubN()%(2*np.pi), geom.valid_elevAngVSubN(), geom.pathLens() ) - altDec1km= np.full_like(altDec, 1.05) # if config.detector.method == "Optical" or config.detector.method == "Both": if config.detector.optical.enable: logv("Computing [green] EAS Optical Cherenkov light.[/]") numPEs, costhetaChEff, ch_photons, ch_times, ch_costheta = eas( beta_tr, - altDec1km, + altDec, showerEnergy, init_lat, init_long, @@ -299,17 +245,20 @@ def add_meta(self, name: str, value: Any, comment: str): #store=sw, plot=to_plot, ) - import matplotlib.pyplot as plt - - mask = (altDec1km < 1) | (altDec1km > 1.1) | (showerEnergy*2 < 100) | (showerEnergy*2 > 200) - mask= ~mask - distance= point_to_line_distances(detcoords[mask], beta_tr[mask], azimuth[mask])/1000 #km - angleatdecay=angle_at_decay(detcoords[mask],altDec1km[mask], beta_tr[mask], azimuth[mask]) - logtauenergy = np.log10(tauEnergy)+9 - print(logtauenergy[mask]) - print('HOLA',ch_photons, ch_times, ak.num(ch_photons,axis=1), ak.num(ch_times)) - ch_phot_integrated=ak.sum(ch_photons, axis=1) - ch_theta=np.arccos(ch_costheta) + + ch_phot_integrated=ak.sum(ch_photons, axis=1) #sum across all wavelengths + ch_theta=np.arccos(ch_costheta) #angle between the z-axis (zenith of spot on ground) and the vector from the axis to the detector + totphot=ak.sum(ch_phot_integrated,axis=1) #total photons arriving at detector for each event (sum across shower axis) + + #Some plots for checks + + #mask = (altDec1km < 1) | (altDec1km > 1.1) | (showerEnergy*2 < 100) | (showerEnergy*2 > 200) + #mask= ~mask + #logtauenergy = np.log10(tauEnergy)+9 + #distance= point_to_line_distances(detcoords[mask], beta_tr[mask], azimuth[mask])/1000 #km + #angleatdecay=angle_at_decay(detcoords[mask],altDec1km[mask], beta_tr[mask], azimuth[mask]) + #data = np.column_stack((totphot, distance, np.degrees(angleatdecay))) + #np.savetxt('cherenkov_data.txt', data, delimiter=',', header='totphot,distance,angleatdecay_degrees', comments='', fmt='%f') #print(len(ch_photons), len(ch_photons[0]), len(ch_photons[0][0])) #print(len(ch_phot_integrated), len(ch_phot_integrated[0]), ak.num(ch_phot_integrated)) """for i in range(len(ch_photons)): @@ -329,11 +278,8 @@ def add_meta(self, name: str, value: Any, comment: str): plt.savefig(f'IntegratedPhotons_{i}.png') print('TOTAL CHERENKOV PHOTONS AT DETECTOR: ', ak.sum(ch_phot_integrated[i])) print(ch_phot_integrated[i])""" - totphot=ak.sum(ch_phot_integrated,axis=1) - data = np.column_stack((totphot, distance, np.degrees(angleatdecay))) - #n = data.shape[0] # Get the number of rows - np.savetxt('data100000.txt', data, delimiter=',', header='totphot,distance,angleatdecay_degrees', comments='', fmt='%f') + ''' plt.figure() plt.plot(distance, totphot,'.') plt.xlabel('Distance to shower axis (km)') @@ -361,9 +307,8 @@ def add_meta(self, name: str, value: Any, comment: str): plt.yscale('log') plt.grid() plt.title('Total photons vs angle to shower axis, tau energy 100-200 EeV') - plt.savefig('photonsvsangle_lims.png') - #NEXT CHECK ANGLE BETWEEN SHOWER AXIS AND DETECTOR AT DECAY POINT OR AT XMAX - # Save to CSV with header + plt.savefig('photonsvsangle_lims.png') ''' + logv("Computing [green] Optical Monte Carlo Integral.[/]") mcint, mcintgeo, passEV, mcunc = geom.mcintegral( diff --git a/src/nuspacesim/simulation/eas_optical/cphotang.py b/src/nuspacesim/simulation/eas_optical/cphotang.py index ca1fb36..d97a3b7 100644 --- a/src/nuspacesim/simulation/eas_optical/cphotang.py +++ b/src/nuspacesim/simulation/eas_optical/cphotang.py @@ -604,10 +604,10 @@ def chasm(self, betaE, gramsum, RN, detcoord, azimuth_tr): sim.add(ch.SphericalCounters(detcoord[np.newaxis, :], np.sqrt(1/np.pi))) sim.add(ch.Yield(270, 1000, N_bins=100)) sig = sim.run(mesh=False, att=True) - ch_photons=np.array(sig.photons.sum(axis=0)) + ch_photons=np.array(sig.photons.sum(axis=0)) # ch_times = np.array(sig.times.sum(axis=0)) ch_costheta=np.array(sig.cos_theta.sum(axis=0).sum(axis=0)) #cosine of the angle between the z-axis (zenith of spot on ground) and the vector from the axis to the detector - + return np.array(ch_photons), np.array(ch_times), np.array(ch_costheta) def run(self, betaE, alt, Eshow100PeV, lat, long, detcoords, azimuth, cloudf=None): diff --git a/src/nuspacesim/simulation/eas_optical/eas.py b/src/nuspacesim/simulation/eas_optical/eas.py index 4d36cdc..02bd106 100644 --- a/src/nuspacesim/simulation/eas_optical/eas.py +++ b/src/nuspacesim/simulation/eas_optical/eas.py @@ -98,7 +98,7 @@ def __call__( # Mask out-of-bounds events. Do not pass to CphotAng. Instead use # Default values for dphots and thetaCh - mask = (altDec < 1) | (altDec > 1.1) | (showerEnergy*2 < 100) | (showerEnergy*2 > 200) + mask = (altDec < 0) | (altDec > 20) # | (showerEnergy*2 < 100) | (showerEnergy*2 > 200) mask = ~mask # phots and theta arrays with default 0 and 1.5 values.