Merge remote-tracking branch 'origin/develop' into add_fallback_anten…

…naserver

.github/workflows/deploy_documentation.yaml

@@ -37,6 +37,7 @@ jobs:
           python install_dev.py --install --dev documentation proposal --no-interactive
           export PYTHONPATH=$PWD:$PYTHONPATH
           echo $PYTHONPATH
+          NuRadioMC/SignalProp/install.sh
 
       - name: Debugging
         run: |
@@ -59,7 +60,7 @@ jobs:
         run: |
           export PYTHONPATH=$(pwd):$PYTHONPATH
           cd documentation/
-          python make_docs.py
+          python make_docs.py --debug
 
       - name: No Jekyll
         if: ${{ github.ref == 'refs/heads/develop'}}

NuRadioMC/SignalGen/ARZ/ARZ.py

@@ -494,7 +494,8 @@ def get_shower_profile(self, shower_energy, shower_type, iN):
 
 
     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):
+                       same_shower=False, iN=None, output_mode='trace', maximum_angle=20 * units.deg,
+                       profile_depth=None, profile_ce=None):
         """
         calculates the electric-field Askaryan pulse from a charge-excess profile
 
@@ -537,6 +538,13 @@ def get_time_trace(self, shower_energy, theta, N, dt, shower_type, n_index, R, s
             Maximum angular difference allowed between the observer angle and the Cherenkov angle.
             If the difference is greater, the function returns an empty trace.
 
+        profile_depth: (optional) array of floats
+            shower depth values of the charge excess profile
+            if provided, profile_ce must also be provided
+            if provided, the function will not use the library to get the profile but use the provided profile
+        profile_ce: (optional) array of floats
+            charge-excess values of the charge excess profile
+
         Returns
         -------
         efield_trace: array of floats
@@ -552,43 +560,47 @@ def get_time_trace(self, shower_energy, theta, N, dt, shower_type, n_index, R, s
         # should not trigger, we return an empty trace for angular differences > 20 degrees.
         cherenkov_angle = np.arccos(1 / n_index)
 
-        # determine closes available energy in shower library
-        energies = np.array([*self._library[shower_type]])
-        iE = np.argmin(np.abs(energies - shower_energy))
-        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)):
-            if(same_shower):
-                if(shower_type in self._random_numbers):
-                    iN = self._random_numbers[shower_type]
-                    logger.info("using previously used shower {}/{}".format(iN, N_profiles))
-                else:
-                    logger.warning("no previous random number for shower type {} exists. Generating a new random number.".format(shower_type))
-                    iN = self._random_generator.randint(N_profiles)
-                    self._random_numbers[shower_type] = iN
-                    logger.info("picking profile {}/{} randomly".format(iN, N_profiles))
-            else:
-                iN = self._random_generator.randint(N_profiles)
-                self._random_numbers[shower_type] = iN
-                logger.info("picking profile {}/{} randomly".format(iN, N_profiles))
-        else:
-            iN = int(iN)  # saveguard against iN being a float
-            logger.info("using shower {}/{} as specified by user".format(iN, N_profiles))
-            self._random_numbers[shower_type] = iN
-
         # we always need to generate a random shower realization. The second ray tracing solution might be closer
         # to the cherenkov angle, but NuRadioMC will reuse the shower realization of the first ray tracing solution.
         if np.abs(theta - cherenkov_angle) > maximum_angle:
             logger.info(f"viewing angle {theta/units.deg:.1f}deg is more than {maximum_angle/units.deg:.1f}deg away from the cherenkov cone. Returning zero trace.")
             empty_trace = np.zeros((3, N))
             return empty_trace
 
-        profile_depth = profiles['depth']
-        profile_ce = profiles['charge_excess'][iN] * rescaling_factor
+        # determine closes available energy in shower library
+        if profile_depth is None:
+            energies = np.array([*self._library[shower_type]])
+            iE = np.argmin(np.abs(energies - shower_energy))
+            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)):
+                if(same_shower):
+                    if(shower_type in self._random_numbers):
+                        iN = self._random_numbers[shower_type]
+                        logger.info("using previously used shower {}/{}".format(iN, N_profiles))
+                    else:
+                        logger.warning("no previous random number for shower type {} exists. Generating a new random number.".format(shower_type))
+                        iN = self._random_generator.randint(N_profiles)
+                        self._random_numbers[shower_type] = iN
+                        logger.info("picking profile {}/{} randomly".format(iN, N_profiles))
+                else:
+                    iN = self._random_generator.randint(N_profiles)
+                    self._random_numbers[shower_type] = iN
+                    logger.info("picking profile {}/{} randomly".format(iN, N_profiles))
+            else:
+                iN = int(iN)  # saveguard against iN being a float
+                logger.info("using shower {}/{} as specified by user".format(iN, N_profiles))
+                self._random_numbers[shower_type] = iN
+            profile_depth = profiles['depth']
+            profile_ce = profiles['charge_excess'][iN] * rescaling_factor
+        else: # if profile_depth is provided, we don't need to use the library
+            if profile_ce is None:
+                raise ValueError("if profile_depth is provided, profile_ce must also be provided")
+            logger.info("using provided charge-excess profile, shower_energy and iN will be ignored.")
 
         xmax = profile_depth[np.argmax(profile_ce)]
 

NuRadioReco/detector/antennapattern.py

@@ -1456,7 +1456,7 @@ def _get_antenna_response_vectorized_raw(self, freq, theta, phi, group_delay='fr
 
             index = np.argmax(freq > self._cutoff_freq)
             Gain = np.ones_like(freq)
-            from scipy.signal import hann
+            from scipy.signal.windows import hann
             gain_filter = hann(2 * index)
             Gain[:index] = gain_filter[:index]
 

NuRadioReco/detector/detector_base.py

@@ -606,7 +606,8 @@ def get_site_coordinates(self, station_id):
             'auger': (-35.10, -69.55),
             'mooresbay': (-78.74, 165.09),
             'southpole': (-90., 0.),
-            'summit': (72.57, -38.46)
+            'summit': (72.57, -38.46),
+            'lofar': (52.92, 6.87)
         }
         site = self.get_site(station_id)
         if site in sites.keys():

NuRadioReco/detector/detector_browser/detector_provider.py

@@ -47,9 +47,9 @@ def set_generic_detector(self, filename, default_station, default_channel, assum
             ID of the channel to be set as default channel
         """
         import NuRadioReco.detector.generic_detector
-        self.__detector = NuRadioReco.detector.generic_detector.GenericDetector.__new__(
-            NuRadioReco.detector.generic_detector.GenericDetector, filename, default_station, 
-            default_channel, assume_inf=assume_inf, antenna_by_depth=antenna_by_depth)
+        self.__detector = NuRadioReco.detector.generic_detector.GenericDetector(
+            filename, default_station, default_channel, assume_inf=assume_inf,
+            antenna_by_depth=antenna_by_depth, create_new=True)
 
     def set_event_file(self, filename):
         """

NuRadioReco/detector/response.py

@@ -174,7 +174,7 @@ def __init__(self, frequency, y, y_unit, time_delay=0, weight=1,
                         f'{datetime.datetime.utcnow().strftime("%Y%m%dT%H%M%S")}_debug.png', transparent=False)
             plt.close()
 
-    def __call__(self, freq, component_name=None):
+    def __call__(self, freq, component_names=None, blacklist=True):
         """
         Returns the complex response for a given frequency.
 
@@ -184,17 +184,36 @@ def __call__(self, freq, component_name=None):
         freq: list(float)
             The frequencies for which to get the response.
 
+        component_names: list(str) or None (Default: None)
+            Only return the (combined) response of components selected by their names.
+            List of names to consider or not to consider (depends on `blacklist`).
+            `None` mean no selection.
+
+        blacklist: bool (Default: True)
+            If True (and `component_names is not None`), ignore components selected with `component_names`.
+            If False, only consider components selected with `component_names`.
+
         Returns
+        -------
 
         response: np.array(np.complex128)
             The complex response at the desired frequencies
         """
         response = np.ones_like(freq, dtype=np.complex128)
 
+        if component_names is not None:
+            if isinstance(component_names, str):
+                component_names = [component_names]
+
         for gain, phase, weight, name in zip(self.__gains, self.__phases, self.__weights, self.__names):
 
-            if component_name is not None and component_name != name:
-                continue
+            if component_names is not None:
+                if blacklist:
+                    if name in component_names:  # if name in blacklist skip
+                        continue
+                else:
+                    if name not in component_names:  # if name *not* in whitelist skip
+                        continue
 
             _gain = gain(freq / units.GHz)
 
@@ -205,9 +224,9 @@ def __call__(self, freq, component_name=None):
             response *= (_gain * np.exp(1j * phase(freq / units.GHz))) ** weight
 
         if np.allclose(response, np.ones_like(freq, dtype=np.complex128)):
-            if component_name is not None:
+            if component_names is not None:
                 raise ValueError("Returned response is equal to 1. "
-                                f"Did you requested a non-existing component ({component_name})? "
+                                f"Did you requested a non-existing component ({component_names})? "
                                 f"Options are: {self.__names}")
             else:
                 self.logger.warning("Returned response is equal to 1.")
@@ -405,4 +424,4 @@ def subtract_time_delay_from_response(frequencies, resp, phase=None, time_delay=
 
     resp = gain * np.exp(1j * (phase + 2 * np.pi * time_delay * frequencies))
 
-    return resp
+    return resp

NuRadioReco/framework/trigger.py

@@ -125,7 +125,7 @@ def set_pre_trigger_times(self, pre_trigger_times):
         Set the pre-trigger times
 
         This parameter should only be set if this trigger 
-        determines the readout windows (e.g. by :py:`NuRadioReco.modules.triggerTimeAdjuster`)
+        determines the readout windows (e.g. by :class:`NuRadioReco.modules.triggerTimeAdjuster`)
 
         Parameters
         ----------

NuRadioReco/utilities/fft.py

@@ -1,17 +1,34 @@
-import numpy as np
-
 """
 A wrapper around the numpy fft routines to achive a coherent normalization of the fft
+
 As we have real valued data in the time domain, we use the 'real ffts' that omit the negative frequencies in Fourier
 space. To account for the missing power in the frequency domain, we multiply the frequency spectrum by sqrt(2),
 and divide the iFFT with 1/sqrt(2) accordingly. The frequency spectrum is divided by the sampling rate so that the
 units if the spectrum are volts/GHz instead of volts/bin and the voltage in the frequency domain is independent of the
 sampling rate.
-Then, a calculation of the power leads the same result in
+
+Then, a calculation of the energy leads the same result in
 the time and frequency domain, i.e.
-np.sum(trace**2) * dt = np.sum(spectrum**2/dt**2) * df
+
+.. code-block::
+
+    np.sum(trace**2) * dt = np.sum(spectrum**2) * df
+
+with units of :math:`V^2/GHz`. In order to obtain the correct units
+(i.e. energy in eV) one additionally has to divide by the impedance ``Z``.
+
+Note that in our convention, the above equality holds only **approximately**, as also the zero-frequency
+and the Nyquist frequency are multiplied by sqrt(2). This effect is however very small and negligible
+in practice.
+
+Finally, we remark that our normalization ensures implies the FFT describes the **energy** spectral density;
+this is in contrast to another common convention where the Fourier transform describes the
+**power** spectral density. One should keep this in mind especially when working with e.g. noise temperatures
+which are defined using the latter convention.
+
 """
 
+import numpy as np
 
 def time2freq(trace, sampling_rate):
     """

NuRadioReco/utilities/logging.py

@@ -19,8 +19,8 @@ def addLoggingLevel(levelName, levelNum, methodName=None):
     raise an `AttributeError` if the level name is already an attribute of the
     `logging` module or if the method name is already present
 
-    Example
-    -------
+    Examples
+    --------
     >>> addLoggingLevel('TRACE', logging.DEBUG - 5)
     >>> logging.getLogger(__name__).setLevel("TRACE")
     >>> logging.getLogger(__name__).trace('that worked')

changelog.txt

@@ -26,6 +26,7 @@ can then iterate over a channel group, which yields all channels for a given gro
 - Added LOFAR antenna pattern and function to parse TXT simulation file
 - The assume_inf and antenna_by_depth Detector keywords are now saved in the nur file, and will be
 loaded again when reading in the Detector from a nur file
+- Added LOFAR coordinates to Detector site coordinates
 
 bugfixes:
 - Fixed bug in get_travel_time in directRayTracing propagation module

documentation/make_docs.py

@@ -6,7 +6,7 @@
 import re
 
 logging.basicConfig()
-logger = logging.getLogger(name="Make_docs")
+logger = logging.getLogger(name="documentation.make_docs")
 logger.setLevel(logging.INFO)
 
 # we do some error classification, because we don't want to fail on all errors
@@ -55,20 +55,21 @@
             " Useful if you are only modifying or adding (not moving/removing) pages.")
     )
     argparser.add_argument(
-        '--debug', default=False, const=True, action='store_const',
+        '--debug', '-v', default=0, action='count',
         help="Store full debugging output in make_docs.log."
         )
     parsed_args = argparser.parse_args()
+
+    pipe_stdout = subprocess.PIPE # by default, hide the stdout, unless we want LOTS of debugging
     if parsed_args.debug: # set up the logger to also write output to make_docs.log
         logfile = 'make_docs.log'
         with open(logfile, 'w') as file:
             pass
         file_logger = logging.FileHandler(logfile)
         file_logger.setLevel(logger.getEffectiveLevel())
         logger.addHandler(file_logger)
-        pipe_stdout = None
-    else:
-        pipe_stdout = subprocess.PIPE # hide the stdout
+        if parsed_args.debug > 1:
+            pipe_stdout = None
 
     doc_path = os.path.dirname(os.path.realpath(__file__))
     os.chdir(doc_path)
@@ -111,7 +112,9 @@
     if not parsed_args.no_clean:
         logger.info('Removing old \'build\' directory...')
         subprocess.check_output(['make', 'clean'])
+    logger.info('Building documentation. This may take a while...')
     sphinx_log = subprocess.run(['make', 'html'], stderr=subprocess.PIPE, stdout=pipe_stdout)
+    logger.info('Finished building documentation.')
 
     # we write out all the sphinx errors to sphinx-debug.log, and parse these
     with open('sphinx-debug.log') as f:

documentation/source/Introduction/pages/installation.rst

@@ -27,8 +27,8 @@ The pip installation will also install all core dependencies.
 
 Development version
 ---------------------------
-The most recent version of ``NuRadioMC`` is available on `github <github.com>`_. 
-It can be downloaded manually from the `repository website <https://github.com/nu-radio/NuRadioMC.git>`_,
+The most recent version of ``NuRadioMC`` is available on `github <https://github.com/nu-radio/NuRadioMC.git>`__.
+It can be downloaded manually from the `repository website <https://github.com/nu-radio/NuRadioMC/releases/latest>`__,
 or cloned using ``git``
 
 .. code-block:: Bash

documentation/source/conf.py

@@ -239,7 +239,7 @@
 autodoc_mock_imports = [
     'ROOT', 'mysql-python', 'pygdsm', 'MySQLdb', 'healpy', 'scripts',
     'uproot', 'radiopropa', 'plotly', 'past',
-    'nifty5', 'mattak'
+    'nifty5', 'mattak', 'MCEq', 'crflux'
     ]
 # Raise warnings if any cross-references are broken
 nitpicky = True