Merge branch 'develop' into parameter_storage

.github/workflows/run_tests.yaml

@@ -190,6 +190,12 @@ jobs:
         export GSLDIR=$(gsl-config --prefix)
         export PYTHONPATH=$(pwd):$PYTHONPATH
         NuRadioReco/test/trigger_tests/run_trigger_test.sh
+    - name: "Test channelGalacticNoiseAdder"
+      if: always()
+      run: |
+        export PYTHONPATH=$(pwd):$PYTHONPATH
+        python install_dev.py --dev galacticnoise --no-interactive
+        python NuRadioReco/examples/StandAloneScripts/B01galacticnoiseadder.py
     - name: "Test all examples"
       if: always()
       run: |

NuRadioMC/SignalProp/propagation_base_class.py

@@ -25,7 +25,10 @@ class initilization
         medium: medium class
             class describing the index-of-refraction profile
         attenuation_model: string
-            signal attenuation model (so far only "SP1" is implemented)
+            if this parameter is also defined in the config file, the value from the config file
+            will be used. If not, the value from this parameter will be used.
+
+            signal attenuation model
         log_level: logging object
             specify the log level of the ray tracing class
 
@@ -36,10 +39,16 @@ class describing the index-of-refraction profile
 
             default is NOTSET (ie global control)
         n_frequencies_integration: int
-            the number of frequencies for which the frequency dependent attenuation
+            if this parameter is also defined in the config file, the value from the config file
+            will be used. If not, the value from this parameter will be used.
+
+            This parameter specifies the number of frequencies for which the frequency dependent attenuation
             length is being calculated. The attenuation length for all other frequencies
             is obtained via linear interpolation.
         n_reflections: int (default 0)
+            if this parameter is also defined in the config file, the value from the config file
+            will be used. If not, the value from this parameter will be used.
+
             in case of a medium with a reflective layer at the bottom, how many reflections should be considered
         config: nested dictionary
             loaded yaml config file
@@ -60,6 +69,19 @@ class describing the index-of-refraction profile
         self._n_reflections = n_reflections
 
         self._config = config
+        if self._config is not None:
+            if 'n_freq' in self._config['propagation']:
+                if n_frequencies_integration is not None:
+                    self.__logger.warning(f"overriding n_frequencies_integration from config file from {n_frequencies_integration} to {self._config['propagation']['n_freq']}")
+                self._n_frequencies_integration = self._config['propagation']['n_freq']
+            if 'n_reflections' in self._config['propagation']:
+                if n_reflections is not None:
+                    self.__logger.warning(f"overriding n_reflections from config file from {n_reflections} to {self._config['propagation']['n_reflections']}")
+                self._n_reflections = self._config['propagation']['n_reflections']
+            if 'attenuation_model' in self._config['propagation']:
+                if attenuation_model is not None:
+                    self.__logger.warning(f"overriding attenuation_model from config file from {attenuation_model} to {self._config['propagation']['attenuation_model']}")
+                self._attenuation_model = self._config['propagation']['attenuation_model']
         self._detector = detector
         self._max_detector_frequency = None
         if self._detector is not None:

NuRadioMC/examples/Sensitivities/E2_fluxes3.py

@@ -528,40 +528,11 @@ def get_ice_cube_hese_range():
 ara_4year_limit *= units.eV * units.cm ** -2 * units.second ** -1 * units.sr ** -1
 ara_4year_limit *= energyBinsPerDecade
 
-# ARA 2023 projection
-'''
-This estimate is built by using the actual recorded livetime for the ARA stations 
-through June 2021. Specifically:
-1747 days of A1
-5627 days of A2 + A3 + A4
-826 days of A5
-
-And then adding projected livetime
-A1: 7/12 of a year for 2021, then 1 year of 2022, and 1 year of 2023
-A2: no more data for 2021, no data for 2022, 1 year of 2023
-A3: 7/12 of a year for 2021, then 1 year of 2022, and 1 year of 2023
-A4: no more data for 2021, no data for 2022, 1 year of 2023
-A5: 7/12 of a year for 2021, then 1 year of 2022, and 1 year of 2023
-
-We do include different effective areas for A1, A2/3/4, and A5,
-since A1 is smaller (only being at 100m), while A5 is larger (having the phased array).
-
-We also included the trigger level and analysis level estimate.
-The analysis level option assumes the A2 analysis efficiency for stations A1-4,
-and a (preliminary) analysis efficiency estimate from the phased-array analysis.
-
-'''
-ara_2023_E_TL, ara_2023_limit_TL, t1, t2 = np.loadtxt(os.path.join(os.path.dirname(__file__), 'data', "limit_ara_2023_projected_trigger.txt"), unpack=True)
-ara_2023_E_TL *= units.GeV
-ara_2023_limit_TL *= units.GeV * units.cm ** -2 * units.second ** -1 * units.sr ** -1
-ara_2023_limit_TL *= energyBinsPerDecade
-ara_2023_limit_TL *= 2.44  # convert to 90%CL limit to be comparable with information from other experiments
-
-ara_2023_E, ara_2023_limit, t1, t2 = np.loadtxt(os.path.join(os.path.dirname(__file__), 'data', "limit_ara_2023_projected_analysis.txt"), unpack=True)
-ara_2023_E *= units.GeV
-ara_2023_limit *= units.GeV * units.cm ** -2 * units.second ** -1 * units.sr ** -1
-ara_2023_limit *= energyBinsPerDecade
-ara_2023_limit *= 2.44  # convert to 90%CL limit to be comparable with information from other experiments
+# ARA phased array 6-month limit
+ara_PA_6month_E, ara_PA_6month_limit, t1, t2 = np.loadtxt(os.path.join(os.path.dirname(__file__), 'data', "limit_ARA_PA_6months.txt"), unpack=True)
+ara_PA_6month_E *= units.eV
+ara_PA_6month_limit *= units.eV * units.cm ** -2 * units.second ** -1 * units.sr ** -1
+ara_PA_6month_limit *= energyBinsPerDecade
 
 ARIANNA_HRA = np.array([[1.00000003e+07, 3.16228005e+07, 9.99999984e+07, 3.16227997e+08,
                          9.99999984e+08, 3.16228010e+09, 9.99999998e+09, 3.16228010e+10,
@@ -631,8 +602,7 @@ def get_E2_limit_figure(diffuse=True,
                         show_anita_I_IV_limit=True,
                         show_auger_limit=True,
                         show_ara=True,
-                        show_ara_2023=False,
-                        show_ara_2023_TL=False,
+                        show_ara_PA=False,
                         show_arianna=True,
                         show_neutrino_best_fit=True,
                         show_neutrino_best_case=True,
@@ -773,8 +743,8 @@ def get_E2_limit_figure(diffuse=True,
                             horizontalalignment='left', color='saddlebrown', rotation=50, fontsize=legendfontsize)
         else:
             ax.annotate('GRAND 10k',
-                xy=(1.3e19 * units.eV / plotUnitsEnergy, 2.6e-8), xycoords='data',
-                horizontalalignment='left', va="bottom", color='saddlebrown', rotation=38, fontsize=legendfontsize)
+                xy=(1.5e19 * units.eV / plotUnitsEnergy, 3.4e-8), xycoords='data',
+                horizontalalignment='left', va="bottom", color='saddlebrown', rotation=40, fontsize=legendfontsize)
 
     if show_grand_200k:
         ax.plot(GRAND_energy / plotUnitsEnergy, GRAND_200k / plotUnitsFlux, linestyle=":", color='saddlebrown',
@@ -830,10 +800,10 @@ def get_E2_limit_figure(diffuse=True,
                     horizontalalignment='center', color='dodgerblue', rotation=0, fontsize=legendfontsize)
 
         # Extrapolation
-        energy_placeholder = np.array(([1e14, 1e19])) * units.eV
-        plt.plot(energy_placeholder / plotUnitsEnergy,
-                 ice_cube_nu_fit(energy_placeholder) * energy_placeholder ** 2 / plotUnitsFlux,
-                 color='dodgerblue', linestyle=':')
+        # energy_placeholder = np.array(([1e14, 1e19])) * units.eV
+        # plt.plot(energy_placeholder / plotUnitsEnergy,
+        #          ice_cube_nu_fit(energy_placeholder) * energy_placeholder ** 2 / plotUnitsFlux,
+        #          color='dodgerblue', linestyle=':')
 
         uplimit = np.copy(nu_mu_data[:, 3])
         uplimit[np.where(nu_mu_data[:, 3] == 0)] = 1
@@ -913,27 +883,16 @@ def get_E2_limit_figure(diffuse=True,
             ax.annotate('ARA',
                     xy=(1e18* units.eV / plotUnitsEnergy, 0.7e-6), xycoords='data',
                     horizontalalignment='left', color='indigo', rotation=0, fontsize=legendfontsize)
-    if show_ara_2023:
-        ax.plot(ara_2023_E / plotUnitsEnergy, ara_2023_limit / plotUnitsFlux, color='grey', linestyle='--')
-        if energyBinsPerDecade == 2:
-            ax.annotate('ARA 2023',
-                        xy=(2E16 * units.eV / plotUnitsEnergy, 6e-7), xycoords='data',
-                        horizontalalignment='left', color='grey', rotation=0, fontsize=legendfontsize)
-        else:
-            ax.annotate('ARA 2023',
-                    xy=(4E17 * units.eV / plotUnitsEnergy, 6e-8), xycoords='data',
-                    horizontalalignment='left', color='grey', rotation=0, fontsize=legendfontsize)
-
-    if show_ara_2023_TL:
-        ax.plot(ara_2023_E_TL / plotUnitsEnergy, ara_2023_limit_TL / plotUnitsFlux, color='grey', linestyle='--')
+    if show_ara_PA:
+        ax.plot(ara_PA_6month_E / plotUnitsEnergy, ara_PA_6month_limit / plotUnitsFlux, color='blue', linestyle='-')
         if energyBinsPerDecade == 2:
-            ax.annotate('ARA 2023 \n(TL)',
-                        xy=(1E16 * units.eV / plotUnitsEnergy, 6e-7), xycoords='data',
+            ax.annotate('ARA PA',
+                        xy=(2.95E1 * units.eV / plotUnitsEnergy, 6e-7), xycoords='data',
                         horizontalalignment='left', color='grey', rotation=0, fontsize=legendfontsize)
         else:
-            ax.annotate('ARA 2023 \n(TL)',
-                    xy=(1E16 * units.eV / plotUnitsEnergy, 6e-8), xycoords='data',
-                    horizontalalignment='left', color='grey', rotation=0, fontsize=legendfontsize)
+            ax.annotate('ARA PA',
+                    xy=(2.9E17 * units.eV / plotUnitsEnergy, 5.2e-6), xycoords='data',
+                    horizontalalignment='left', color='blue', rotation=0, fontsize=legendfontsize)
 
     if show_arianna:
         ax.plot(ARIANNA_HRA[:, 0] / plotUnitsEnergy, ARIANNA_HRA[:, 1] / plotUnitsFlux, color='red')

NuRadioMC/examples/Sensitivities/data/limit_ARA_PA_6months.txt

@@ -0,0 +1,20 @@
+# Source: Phased Array Six-Month Real Limit
+# Data type: Upper limit
+# Energy bins per decade: 1
+# Column values:
+# Energy [eV]
+# Flux [eV cm^-2 s^-1 sr^-1]
+# Flux band minimum
+# Flux band maximum
+
+1.00E+16  89771.61706761  0 0
+3.16E+16  27085.56086071  0 0
+1.00E+17  8786.44941194   0 0
+3.16E+17  4822.58991603   0 0
+1.00E+18  3210.43840531   0 0
+3.16E+18  3120.34963572     0 0
+1.00E+19  3255.37433081 0 0
+3.16E+19  4246.17439525 0 0
+1.00E+20  6213.49250206 0 0
+3.16E+20  8835.75610193 0 0
+1.00E+21  17992.49247862  0 0

NuRadioMC/simulation/simulation.py

@@ -1245,10 +1245,8 @@ def __init__(
 
         prop = propagation.get_propagation_module(self._config['propagation']['module'])
         self._propagator = prop(
-            self._ice, self._config['propagation']['attenuation_model'],
+            self._ice,
             log_level=self._log_level_ray_propagation,
-            n_frequencies_integration=int(self._config['propagation']['n_freq']),
-            n_reflections=int(self._config['propagation']['n_reflections']),
             config=self._config,
             detector=self._det
         )