From 4fd919bd5aa36f36bdef36409b50e0f2de65364f Mon Sep 17 00:00:00 2001 From: Sjoerd Bouma Date: Fri, 8 Mar 2024 20:57:49 +0100 Subject: [PATCH 01/12] remove deprecated interface for opening GenericDetector --- NuRadioReco/detector/detector_browser/detector_provider.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/NuRadioReco/detector/detector_browser/detector_provider.py b/NuRadioReco/detector/detector_browser/detector_provider.py index f842a09b4..b0238df9c 100644 --- a/NuRadioReco/detector/detector_browser/detector_provider.py +++ b/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): """ From 7800c56da2dfbb55c99ce6f4cfd76c7a3b0dc93f Mon Sep 17 00:00:00 2001 From: Sjoerd Bouma Date: Tue, 12 Mar 2024 15:15:48 +0100 Subject: [PATCH 02/12] correct description of fft power normalization --- NuRadioReco/utilities/fft.py | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/NuRadioReco/utilities/fft.py b/NuRadioReco/utilities/fft.py index 164fea1bb..5aa0f3b5f 100644 --- a/NuRadioReco/utilities/fft.py +++ b/NuRadioReco/utilities/fft.py @@ -1,17 +1,24 @@ -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 + +Note that this will have 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``. """ +import numpy as np def time2freq(trace, sampling_rate): """ From 6a26eab651e5804ca50decd07ff2d2fb3354f240 Mon Sep 17 00:00:00 2001 From: Sjoerd Bouma Date: Tue, 12 Mar 2024 15:17:26 +0100 Subject: [PATCH 03/12] other minor fixes to docstrings --- NuRadioReco/framework/trigger.py | 2 +- NuRadioReco/utilities/logging.py | 4 ++-- documentation/make_docs.py | 2 +- documentation/source/Introduction/pages/installation.rst | 4 ++-- documentation/source/conf.py | 2 +- 5 files changed, 7 insertions(+), 7 deletions(-) diff --git a/NuRadioReco/framework/trigger.py b/NuRadioReco/framework/trigger.py index 2ab08c57f..203ac1f8d 100644 --- a/NuRadioReco/framework/trigger.py +++ b/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 ---------- diff --git a/NuRadioReco/utilities/logging.py b/NuRadioReco/utilities/logging.py index 0b5aae15b..6801f795b 100644 --- a/NuRadioReco/utilities/logging.py +++ b/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') diff --git a/documentation/make_docs.py b/documentation/make_docs.py index b7909c628..a4fc0fe18 100644 --- a/documentation/make_docs.py +++ b/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 diff --git a/documentation/source/Introduction/pages/installation.rst b/documentation/source/Introduction/pages/installation.rst index 62fd29dda..8ad0c8600 100644 --- a/documentation/source/Introduction/pages/installation.rst +++ b/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 `_. -It can be downloaded manually from the `repository website `_, +The most recent version of ``NuRadioMC`` is available on `github `__. +It can be downloaded manually from the `repository website `__, or cloned using ``git`` .. code-block:: Bash diff --git a/documentation/source/conf.py b/documentation/source/conf.py index b5242e833..9f2d666dc 100644 --- a/documentation/source/conf.py +++ b/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 From 095efc5138b64eee75d3073310905978eb12a0fa Mon Sep 17 00:00:00 2001 From: Sjoerd Bouma Date: Tue, 12 Mar 2024 15:18:12 +0100 Subject: [PATCH 04/12] compile C++ raytracer before building docs (hopefully cleaner output) --- .github/workflows/deploy_documentation.yaml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/deploy_documentation.yaml b/.github/workflows/deploy_documentation.yaml index e7696517e..fce5310ba 100644 --- a/.github/workflows/deploy_documentation.yaml +++ b/.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: | From 7b050f23640124592b2ed973a0dc4650c7901e71 Mon Sep 17 00:00:00 2001 From: Sjoerd Bouma Date: Tue, 12 Mar 2024 15:27:50 +0100 Subject: [PATCH 05/12] hopefully make default debug output a bit quieter --- .github/workflows/deploy_documentation.yaml | 2 +- documentation/make_docs.py | 11 +++++++---- 2 files changed, 8 insertions(+), 5 deletions(-) diff --git a/.github/workflows/deploy_documentation.yaml b/.github/workflows/deploy_documentation.yaml index fce5310ba..332f74ac9 100644 --- a/.github/workflows/deploy_documentation.yaml +++ b/.github/workflows/deploy_documentation.yaml @@ -60,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'}} diff --git a/documentation/make_docs.py b/documentation/make_docs.py index a4fc0fe18..208833dda 100644 --- a/documentation/make_docs.py +++ b/documentation/make_docs.py @@ -55,10 +55,12 @@ " 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: @@ -66,9 +68,8 @@ 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: From 02dbeeeb186d47bf5b3a013ba5154fea1ae91b33 Mon Sep 17 00:00:00 2001 From: Karen Terveer <94610272+karenterveer@users.noreply.github.com> Date: Thu, 4 Apr 2024 11:15:52 +0200 Subject: [PATCH 06/12] added LOFAR site longitude and latitude to get_site_coordinates() --- NuRadioReco/detector/detector_base.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/NuRadioReco/detector/detector_base.py b/NuRadioReco/detector/detector_base.py index 6c4e3fd72..417cb2f2b 100644 --- a/NuRadioReco/detector/detector_base.py +++ b/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(): From a4c76bf6db8e8cfcec217f2a25faf2f9e3969605 Mon Sep 17 00:00:00 2001 From: Mitja Desmet Date: Thu, 4 Apr 2024 14:01:04 +0200 Subject: [PATCH 07/12] Update changelog.txt --- changelog.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/changelog.txt b/changelog.txt index c60346b12..1db8d719c 100644 --- a/changelog.txt +++ b/changelog.txt @@ -23,6 +23,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 From e1034eb7937171247d20a6cad7f2539e6e420c51 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Felix=20Schl=C3=BCter?= Date: Tue, 16 Apr 2024 18:19:06 +0200 Subject: [PATCH 08/12] Extend functionality of __call__ method of response class --- NuRadioReco/detector/response.py | 31 +++++++++++++++++++++++++------ 1 file changed, 25 insertions(+), 6 deletions(-) diff --git a/NuRadioReco/detector/response.py b/NuRadioReco/detector/response.py index 646c4fe0f..e27f87d3e 100644 --- a/NuRadioReco/detector/response.py +++ b/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 \ No newline at end of file + return resp From 2ba8828f6a3d7355a922c316ef7436d330c4a5a6 Mon Sep 17 00:00:00 2001 From: Sjoerd Bouma Date: Mon, 29 Apr 2024 12:55:20 +0200 Subject: [PATCH 09/12] import hann window from scipy.signal.windows as import from scipy.signal.hann was deprecated --- NuRadioReco/detector/antennapattern.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NuRadioReco/detector/antennapattern.py b/NuRadioReco/detector/antennapattern.py index 1402a9726..2b0a6dd56 100644 --- a/NuRadioReco/detector/antennapattern.py +++ b/NuRadioReco/detector/antennapattern.py @@ -1467,7 +1467,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] From 8a125f360a47df899c473bb0af6219c0ef4f64b0 Mon Sep 17 00:00:00 2001 From: Sjoerd Bouma Date: Mon, 29 Apr 2024 13:04:41 +0200 Subject: [PATCH 10/12] add link to latest release on github --- documentation/source/Introduction/pages/installation.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/documentation/source/Introduction/pages/installation.rst b/documentation/source/Introduction/pages/installation.rst index 8ad0c8600..79896b5cb 100644 --- a/documentation/source/Introduction/pages/installation.rst +++ b/documentation/source/Introduction/pages/installation.rst @@ -28,7 +28,7 @@ The pip installation will also install all core dependencies. Development version --------------------------- The most recent version of ``NuRadioMC`` is available on `github `__. -It can be downloaded manually from the `repository website `__, +It can be downloaded manually from the `repository website `__, or cloned using ``git`` .. code-block:: Bash From e9ed08d3ef54381fcd229dd32627f3ee32dfc46b Mon Sep 17 00:00:00 2001 From: Sjoerd Bouma Date: Mon, 29 Apr 2024 13:18:27 +0200 Subject: [PATCH 11/12] add some more detail to FFT normalization --- NuRadioReco/utilities/fft.py | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/NuRadioReco/utilities/fft.py b/NuRadioReco/utilities/fft.py index 5aa0f3b5f..825d9232f 100644 --- a/NuRadioReco/utilities/fft.py +++ b/NuRadioReco/utilities/fft.py @@ -14,8 +14,18 @@ np.sum(trace**2) * dt = np.sum(spectrum**2) * df -Note that this will have 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``. +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 From 800dbe92c7b60d440cb0592868e580129f92f8f8 Mon Sep 17 00:00:00 2001 From: Christian Glaser Date: Thu, 16 May 2024 00:51:11 +0200 Subject: [PATCH 12/12] refactor ARZ module to accept custom charge_excess profile (#645) --- NuRadioMC/SignalGen/ARZ/ARZ.py | 74 ++++++++++++++++++++-------------- 1 file changed, 43 insertions(+), 31 deletions(-) diff --git a/NuRadioMC/SignalGen/ARZ/ARZ.py b/NuRadioMC/SignalGen/ARZ/ARZ.py index 92827124c..d4db40849 100644 --- a/NuRadioMC/SignalGen/ARZ/ARZ.py +++ b/NuRadioMC/SignalGen/ARZ/ARZ.py @@ -500,7 +500,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 @@ -543,6 +544,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 @@ -558,34 +566,6 @@ 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: @@ -593,8 +573,40 @@ def get_time_trace(self, shower_energy, theta, N, dt, shower_type, n_index, R, s 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)]