Merge branch 'develop' of github.com:nu-radio/NuRadioMC into feature/…

…block_offsets

NuRadioMC/EvtGen/generate_unforced.py

@@ -328,7 +328,7 @@ def points_in_cylinder(pt1, pt2, r, q):
                 'flavors': [],
                 'energies': []}
     # calculate rotation matrix to transform position on area to 3D
-    mask_int = np.zeros_like(mask, dtype=np.bool)
+    mask_int = np.zeros_like(mask, dtype=bool)
     t0 = time.perf_counter()
     n_cylinder = 0
     for j, i in enumerate(np.arange(n_events, dtype=int)[mask]):
@@ -586,7 +586,7 @@ def points_in_cylinder(pt1, pt2, r, q):
     data_sets['flavors'] = np.ones(np.sum(mask_int))
     data_sets["event_ids"] = np.arange(np.sum(mask_int)) + start_event_id
     data_sets["n_interaction"] = np.ones(np.sum(mask_int), dtype=int)
-    data_sets["vertex_times"] = np.zeros(np.sum(mask_int), dtype=np.float)
+    data_sets["vertex_times"] = np.zeros(np.sum(mask_int), dtype=float)
 
     data_sets["interaction_type"] = inelasticities.get_ccnc(np.sum(mask_int))
     data_sets["inelasticity"] = inelasticities.get_neutrino_inelasticity(np.sum(mask_int))

NuRadioMC/EvtGen/generator.py

@@ -891,7 +891,7 @@ def generate_surface_muons(filename, n_events, Emin, Emax,
 
         data_sets["event_group_ids"] = np.arange(i_batch * max_n_events_batch, i_batch * max_n_events_batch + n_events_batch, dtype=int) + start_event_id
         data_sets["n_interaction"] = np.ones(n_events_batch, dtype=int)
-        data_sets["vertex_times"] = np.zeros(n_events_batch, dtype=np.float)
+        data_sets["vertex_times"] = np.zeros(n_events_batch, dtype=float)
 
         # generate neutrino flavors randomly
 
@@ -928,7 +928,7 @@ def generate_surface_muons(filename, n_events, Emin, Emax,
         if('fiducial_rmax' in attributes):
             mask_phi = mask_arrival_azimuth(data_sets, attributes['fiducial_rmax'])  # this currently only works for cylindrical volumes
         else:
-            mask_phi = np.ones(len(data_sets["event_group_ids"]), dtype=np.bool)
+            mask_phi = np.ones(len(data_sets["event_group_ids"]), dtype=bool)
         # TODO: combine with `get_intersection_volume_neutrino` function
         for iE, event_id in enumerate(data_sets["event_group_ids"]):
             if not mask_phi[iE]:
@@ -1230,7 +1230,7 @@ def generate_eventlist_cylinder(filename, n_events, Emin, Emax,
         data_sets["event_group_ids"] = np.arange(i_batch * max_n_events_batch, i_batch * max_n_events_batch + n_events_batch) + start_event_id
         logger.debug("generating number of interactions")
         data_sets["n_interaction"] = np.ones(n_events_batch, dtype=int)
-        data_sets["vertex_times"] = np.zeros(n_events_batch, dtype=np.float)
+        data_sets["vertex_times"] = np.zeros(n_events_batch, dtype=float)
 
         # generate neutrino flavors randomly
         logger.debug("generating flavors")

NuRadioMC/SignalGen/ARZ/ARZ.py

@@ -63,7 +63,7 @@ def get_vector_potential(
     profile_ce: array of floats
         charge-excess values of the charge excess profile
     shower_type: string (default "HAD")
-        type of shower, either "HAD" (hadronic), "EM" (electromagnetic) or "TAU" (tau lepton induced)
+        type of shower, either "HAD" (hadronic) or "EM" (electromagnetic)
     n_index: float (default 1.78)
         index of refraction where the shower development takes place
     distance: float (default 1km)
@@ -239,7 +239,7 @@ def get_dist_shower(X, z):
         v[0] = vperp_x
         v[1] = vperp_y
         v[2] = vperp_z
-#         v = np.array([vperp_x, vperp_y, vperp_z], dtype=np.float64)
+#         v = np.array([vperp_x, vperp_y, vperp_z], dtype=float)
         """
         Function F_p Eq.(15) PRD paper.
         """
@@ -468,6 +468,37 @@ def set_interpolation_factor2(self, interp_factor):
         """
         self._interp_factor2 = interp_factor
 
+    def get_shower_profile(self, shower_energy, shower_type, iN):
+        """
+        Gets a charge-excess profile
+
+        Parameters
+        ----------
+        shower_energy: float
+            the energy of the shower
+        shower_type: string (default "HAD")
+            type of shower, either "HAD" (hadronic) or "EM" (electromagnetic)
+        iN: int
+            specify shower number
+
+        Returns
+        -------
+        depth, excess: two arrays of floats
+            slant depths and charge profile amplitudes
+        """
+
+        energies = np.array([*self._library[shower_type]])
+        iE = np.argmin(np.abs(energies - shower_energy))
+
+        rescaling_factor = shower_energy / energies[iE]
+
+        profiles = self._library[shower_type][energies[iE]]
+        profile_depth = profiles['depth']
+        profile_ce = profiles['charge_excess'][iN] * rescaling_factor
+
+        return profile_depth, profile_ce
+
+
     def get_time_trace(self, shower_energy, theta, N, dt, shower_type, n_index, R, shift_for_xmax=False,
                        same_shower=False, iN=None, output_mode='trace', maximum_angle=20 * units.deg):
         """
@@ -483,12 +514,8 @@ def get_time_trace(self, shower_energy, theta, N, dt, shower_type, n_index, R, s
             number of samples in the time domain
         dt: float
             size of one time bin in units of time
-        profile_depth: array of floats
-            shower depth values of the charge excess profile
-        profile_ce: array of floats
-            charge-excess values of the charge excess profile
         shower_type: string (default "HAD")
-            type of shower, either "HAD" (hadronic), "EM" (electromagnetic) or "TAU" (tau lepton induced)
+            type of shower, either "HAD" (hadronic) or "EM" (electromagnetic)
         n_index: float (default 1.78)
             index of refraction where the shower development takes place
         R: float (default 1km)
@@ -537,6 +564,7 @@ def get_time_trace(self, shower_energy, theta, N, dt, shower_type, n_index, R, s
         rescaling_factor = shower_energy / energies[iE]
         logger.info("shower energy of {:.3g}eV requested, closest available energy is {:.3g}eV. The amplitude of the charge-excess profile will be rescaled accordingly by a factor of {:.2f}".format(shower_energy / units.eV, energies[iE] / units.eV, rescaling_factor))
         profiles = self._library[shower_type][energies[iE]]
+
         N_profiles = len(profiles['charge_excess'])
 
         if(iN is None or np.isnan(iN)):
@@ -597,7 +625,7 @@ def get_time_trace(self, shower_energy, theta, N, dt, shower_type, n_index, R, s
             logger.error("Tau showers are not yet implemented")
             raise NotImplementedError("Tau showers are not yet implemented")
         else:
-            msg = "showers of type {} are not implemented. Use 'HAD', 'EM' or 'TAU'".format(shower_type)
+            msg = "showers of type {} are not implemented. Use 'HAD', 'EM'".format(shower_type)
             logger.error(msg)
             raise NotImplementedError(msg)
         if self._use_numba:
@@ -669,7 +697,7 @@ def get_vector_potential_fast(self, shower_energy, theta, N, dt, profile_depth,
         profile_ce: array of floats
             charge-excess values of the charge excess profile
         shower_type: string (default "HAD")
-            type of shower, either "HAD" (hadronic), "EM" (electromagnetic) or "TAU" (tau lepton induced)
+            type of shower, either "HAD" (hadronic) or "EM" (electromagnetic)
         n_index: float (default 1.78)
             index of refraction where the shower development takes place
         distance: float (default 1km)
@@ -1136,7 +1164,7 @@ def get_time_trace(self, shower_energy, theta, N, dt, shower_type, n_index, R,
         dt: float
             size of one time bin in units of time
         shower_type: string (default "HAD")
-            type of shower, either "HAD" (hadronic), "EM" (electromagnetic) or "TAU" (tau lepton induced)
+            type of shower, either "HAD" (hadronic) or "EM" (electromagnetic)
         n_index: float (default 1.78)
             index of refraction where the shower development takes place
         R: float (default 1km)

NuRadioMC/SignalProp/CPPAnalyticRayTracing/setup.py

@@ -16,7 +16,7 @@
         Extension('wrapper', ['wrapper.pyx'],
                 include_dirs=[numpy.get_include(), '../../utilities/', str(os.environ['GSLDIR']) + '/include/'],
                 library_dirs=[str(os.environ['GSLDIR']) + '/lib/'],
-                extra_compile_args=['-O3',"-mfpmath=sse"],
+                extra_compile_args=['-O3'],
                 libraries=['gsl', 'gslcblas'],
                 language='c++'
                 ),

NuRadioMC/SignalProp/analyticraytracing.py

@@ -639,7 +639,7 @@ def dt(t, C_0, frequency):
                 mask = frequency > 0
                 freqs = self.__get_frequencies_for_attenuation(frequency, max_detector_freq)
                 gamma_turn, z_turn = self.get_turning_point(self.medium.n_ice ** 2 - C_0 ** -2)
-                print("_use_optimized_calculation", self._use_optimized_calculation)
+                self.__logger.info("_use_optimized_calculation", self._use_optimized_calculation)
 
                 if self._use_optimized_calculation:
                     # The integration of the attenuation factor along the path with scipy.quad is inefficient. The 
@@ -1751,8 +1751,8 @@ def set_start_and_end_point(self, x1, x2):
         if(self._X2[2] < self._X1[2]):
             self._swap = True
             self.__logger.debug('swap = True')
-            self._X2 = np.array(x1, dtype =np.float)
-            self._X1 = np.array(x2, dtype =np.float)
+            self._X2 = np.array(x1, dtype =float)
+            self._X1 = np.array(x2, dtype =float)
 
         dX = self._X2 - self._X1
         self._dPhi = -np.arctan2(dX[1], dX[0])

NuRadioMC/SignalProp/propagation_base_class.py

@@ -92,8 +92,8 @@ def set_start_and_end_point(self, x1, x2):
             stop point of the ray
         """
         self.reset_solutions()
-        self._X1 = np.array(x1, dtype =np.float)
-        self._X2 = np.array(x2, dtype = np.float)
+        self._X1 = np.array(x1, dtype =float)
+        self._X2 = np.array(x2, dtype = float)
         if (self._n_reflections):
             if (self._X1[2] < self._medium.reflection or self._X2[2] < self._medium.reflection):
                 self.__logger.error("start or stop point is below the reflective bottom layer at {:.1f}m".format(

NuRadioMC/examples/03_station_coincidences/A05visualization.py

@@ -79,7 +79,7 @@
         x2 = det.get_relative_position(101, iC)
         if(plot):
             ax.plot([x2[0]], [x2[1]], [x2[2]], 'ko')
-        if(j != 0 and (~(np.array(triggered_deep[iE], dtype=np.bool) & mask))[j]):
+        if(j != 0 and (~(np.array(triggered_deep[iE], dtype=bool) & mask))[j]):
             continue
         vertex = np.array([fin['xx'][iE], fin['yy'][iE], fin['zz'][iE]])
 #         print(fin.keys())

NuRadioMC/simulation/scripts/T05visualize_sim_output.py

@@ -42,14 +42,14 @@
         plot_folder = os.path.join(dirname, 'plots', filename, args.trigger_name[0])
         if(not os.path.exists(plot_folder)):
             os.makedirs(plot_folder)
-        triggered = np.zeros(len(fin['multiple_triggers'][:, 0]), dtype=np.bool)
+        triggered = np.zeros(len(fin['multiple_triggers'][:, 0]), dtype=bool)
         for trigger in args.trigger_name[1:]:
             iTrigger = np.squeeze(np.argwhere(fin.attrs['trigger_names'] == trigger))
-            triggered = triggered | np.array(fin['multiple_triggers'][:, iTrigger], dtype=np.bool)
+            triggered = triggered | np.array(fin['multiple_triggers'][:, iTrigger], dtype=bool)
     else:
         trigger_name = args.trigger_name[0]
         iTrigger = np.argwhere(fin.attrs['trigger_names'] == trigger_name)
-        triggered = np.array(fin['multiple_triggers'][:, iTrigger], dtype=np.bool)
+        triggered = np.array(fin['multiple_triggers'][:, iTrigger], dtype=bool)
         print("\tyou selected '{}'".format(trigger_name))
         plot_folder = os.path.join(dirname, 'plots', filename, trigger_name)
         if(not os.path.exists(plot_folder)):
@@ -147,14 +147,14 @@ def a(theta):
             if(len(args.trigger_name) > 1):
                 print("trigger {} selected which is a combination of {}".format(args.trigger_name[0], args.trigger_name[1:]))
                 trigger_name = args.trigger_name[0]
-                triggered = np.zeros(len(station['multiple_triggers'][:, 0]), dtype=np.bool)
+                triggered = np.zeros(len(station['multiple_triggers'][:, 0]), dtype=bool)
                 for trigger in args.trigger_name[1:]:
                     iTrigger = np.squeeze(np.argwhere(fin.attrs['trigger_names'] == trigger))
-                    triggered = triggered | np.array(station['multiple_triggers'][:, iTrigger], dtype=np.bool)
+                    triggered = triggered | np.array(station['multiple_triggers'][:, iTrigger], dtype=bool)
             else:
                 trigger_name = args.trigger_name[0]
                 iTrigger = np.argwhere(fin.attrs['trigger_names'] == trigger_name)
-                triggered = np.array(station['multiple_triggers'][:, iTrigger], dtype=np.bool)
+                triggered = np.array(station['multiple_triggers'][:, iTrigger], dtype=bool)
                 print("\tyou selected '{}'".format(trigger_name))
 
         ###########################