merged from develop

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++'
                 ),