From 9ad92dfec96809e2e4af0d801a3316e61b77e4a1 Mon Sep 17 00:00:00 2001 From: Francesca Capel Date: Mon, 12 Apr 2021 20:37:13 +0200 Subject: [PATCH 1/7] Add total flux density calc to PowerLawFlux --- icecube_tools/source/flux_model.py | 28 +++++++++++++++++++++++----- 1 file changed, 23 insertions(+), 5 deletions(-) diff --git a/icecube_tools/source/flux_model.py b/icecube_tools/source/flux_model.py index aedb54d..c5168cf 100644 --- a/icecube_tools/source/flux_model.py +++ b/icecube_tools/source/flux_model.py @@ -40,7 +40,7 @@ def __init__( upper_energy=np.inf, ): """ - Power law flux models. + Power law flux models. :param normalisation: Flux normalisation [GeV^-1 cm^-2 s^-1 sr^-1] or [GeV^-1 cm^-2 s^-1] for point sources. :param normalisation energy: Energy at which flux is normalised [GeV]. @@ -82,7 +82,7 @@ def spectrum(self, energy): def integrated_spectrum(self, lower_energy_bound, upper_energy_bound): """ \int spectrum dE over finite energy bounds. - + :param lower_energy_bound: [GeV] :param upper_energy_bound: [GeV] """ @@ -108,11 +108,29 @@ def integrated_spectrum(self, lower_energy_bound, upper_energy_bound): - np.power(lower_energy_bound, 1 - self._index) ) + def total_flux_density(self): + """ + Total flux density in units of + [GeV cm^-2 s^-1] + """ + + norm = self._normalisation + index = self._index + lower, upper = self._lower_energy, self._upper_energy + + if index == 2: + # special case + int_norm = norm / (np.power(self._normalisation_energy, -index)) + return int_norm * (np.log(upper / lower)) + + int_norm = norm / (np.power(self._normalisation_energy, -index) * (2 - index)) + return int_norm * (np.power(upper, 2 - index) - np.power(lower, 2 - index)) + def sample(self, N): """ Sample energies from the power law. Uses inverse transform sampling. - + :param min_energy: Minimum energy to sample from [GeV]. :param N: Number of samples. """ @@ -216,7 +234,7 @@ def spectrum(self, energy): def integrated_spectrum(self, lower_energy_bound, upper_energy_bound): """ \int spectrum dE over finite energy bounds. - + :param lower_energy_bound: [GeV] :param upper_energy_bound: [GeV] """ @@ -272,7 +290,7 @@ def sample(self, N): """ Sample energies from the power law. Uses inverse transform sampling. - + :param N: Number of samples. """ From 8f9c4d5bc6fdc1f20be90f98e663de40f24680e8 Mon Sep 17 00:00:00 2001 From: Francesca Capel Date: Wed, 21 Apr 2021 17:43:26 +0200 Subject: [PATCH 2/7] Add a reader for the 2013 HESE effective areas --- icecube_tools/detector/effective_area.py | 82 +++++++++++++++++++++++- 1 file changed, 81 insertions(+), 1 deletion(-) diff --git a/icecube_tools/detector/effective_area.py b/icecube_tools/detector/effective_area.py index 1cbcabb..99867e5 100644 --- a/icecube_tools/detector/effective_area.py +++ b/icecube_tools/detector/effective_area.py @@ -1,4 +1,5 @@ import numpy as np +import os from abc import ABC, abstractmethod from icecube_tools.utils.data import ( @@ -12,11 +13,12 @@ effective area information. """ +R2013_AEFF_FILENAME = "effective_areas" R2015_AEFF_FILENAME = "effective_area.h5" R2018_AEFF_FILENAME = "TabulatedAeff.txt" BRAUN2008_AEFF_FILENAME = "AeffBraun2008.csv" -_supported_dataset_ids = ["20150820", "20181018"] +_supported_dataset_ids = ["20131121", "20150820", "20181018"] class IceCubeAeffReader(ABC): @@ -60,6 +62,84 @@ def read(self): pass +class R2013AeffReader(IceCubeAeffReader): + """ + Reader for the 2013 Nov 21 release. + Link: https://icecube.wisc.edu/data-releases/2013/11/ + search-for-contained-neutrino-events-at-energies-above + -30-tev-in-2-years-of-data. + """ + + def __init__(self, filename, **kwargs): + """ + Here, filename is the path to the folder, as + the effective areas are given in a bunch of + separate files. + """ + + if "nu_flavors" in kwargs: + + self.nu_flavors = kwargs["nu_flavors"] + + else: + + self.nu_flavors = ["numu", "nue", "nutau"] + + self._cosz_range = np.linsapce(-1, 1, 21) + + self._fname_str = "_cosZenRange_from_%.1f_to_%.1f.txt" + + super().__init__(filename) + + def read(self): + + self.cos_zenith_bins = self._cosz_range + + self.true_energy_bins = self._get_true_energy_bins() + + self.effective_area_values = np.zeros( + (self.true_energy_bins.size - 1, self.cos_zenith_bins.size - 1) + ) + + for i, bounds in enumerate( + zip(self.cos_zenith_bins[:-1], self.cos_zenith_bins[1:]) + ): + + l, u = bounds + + for nu_flavor in self.nu_flavors: + + tmp_file_name = nu_flavor + self._fname_str % (l, u) + + file_name = os.path.join(self.filename, tmp_file_name) + + tmp_read = np.loadtxt(file_name, skiprows=2).T + + self.effective_area_values.T[i] += tmp_read[2] + + # Assume equal flavour ratio in flux + self.effective_area_values.T[i] /= 3 + + def _get_true_energy_bins(self): + """ + These are the same in all files, so can + just read out once. + """ + + tmp_file_name = self.nu_flavors[0] + self._fname_str % ( + self._cosz_range[0], + self._cosz_range[1], + ) + + file_name = os.path.join(self.filename, tmp_file_name) + + tmp_read = np.loadtxt(file_name, skiprows=2).T + + energy_bins = np.append(tmp_read[0], tmp_read[1][-1]) + + return energy_bins + + class R2015AeffReader(IceCubeAeffReader): """ Reader for the 2015 Aug 20 release. From e9b27647831bbf0d302aa81346108668d13e68a0 Mon Sep 17 00:00:00 2001 From: Francesca Capel Date: Wed, 21 Apr 2021 17:53:26 +0200 Subject: [PATCH 3/7] Test new R2013AeffReader and fix bugs --- icecube_tools/detector/effective_area.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/icecube_tools/detector/effective_area.py b/icecube_tools/detector/effective_area.py index 99867e5..0d4d8bd 100644 --- a/icecube_tools/detector/effective_area.py +++ b/icecube_tools/detector/effective_area.py @@ -85,9 +85,9 @@ def __init__(self, filename, **kwargs): self.nu_flavors = ["numu", "nue", "nutau"] - self._cosz_range = np.linsapce(-1, 1, 21) + self._cosz_range = np.linspace(-1, 1, 21) - self._fname_str = "_cosZenRange_from_%.1f_to_%.1f.txt" + self._fname_str = "_cosZenRange_from_%+.1f_to_%+.1f.txt" super().__init__(filename) @@ -111,7 +111,7 @@ def read(self): tmp_file_name = nu_flavor + self._fname_str % (l, u) - file_name = os.path.join(self.filename, tmp_file_name) + file_name = os.path.join(self._filename, tmp_file_name) tmp_read = np.loadtxt(file_name, skiprows=2).T @@ -131,7 +131,7 @@ def _get_true_energy_bins(self): self._cosz_range[1], ) - file_name = os.path.join(self.filename, tmp_file_name) + file_name = os.path.join(self._filename, tmp_file_name) tmp_read = np.loadtxt(file_name, skiprows=2).T From 36ac698d4ecd78840cf27bcce0c089039f471475 Mon Sep 17 00:00:00 2001 From: Francesca Capel Date: Wed, 21 Apr 2021 18:16:06 +0200 Subject: [PATCH 4/7] Add extraction of inner zip files for downloaded datasets and find_folders function --- icecube_tools/utils/data.py | 31 ++++++++++++++++++++++++++++++- 1 file changed, 30 insertions(+), 1 deletion(-) diff --git a/icecube_tools/utils/data.py b/icecube_tools/utils/data.py index bb03da2..56a5961 100644 --- a/icecube_tools/utils/data.py +++ b/icecube_tools/utils/data.py @@ -161,14 +161,22 @@ def fetch(self, datasets, overwrite=False, write_to=None): # Delete zipfile os.remove(local_path) - # Check for further tarfiles in the extraction + # Check for further compressed files in the extraction tar_files = find_files(dataset_dir, ".tar") + zip_files = find_files(dataset_dir, ".zip") + for tf in tar_files: tar = tarfile.open(tf) tar.extractall(os.path.splitext(tf)[0]) + for zf in zip_files: + + with ZipFile(zf, "r") as zip_ref: + + zip_ref.extractall(zf[:-4]) + crawl_delay() if write_to: @@ -225,3 +233,24 @@ def find_files(directory, keyword): found_files.append(os.path.join(root, f)) return found_files + + +def find_folders(directory, keyword): + """ + Find subfolders in a directory that + contain a keyword. + """ + + found_folders = [] + + for root, dirs, files in os.walk(directory): + + if dirs: + + for d in dirs: + + if keyword in d: + + found_folders.append(os.path.join(root, d)) + + return found_folders From a72e9533174846f97931a2de23ce40876bbedbd3 Mon Sep 17 00:00:00 2001 From: Francesca Capel Date: Wed, 21 Apr 2021 18:20:32 +0200 Subject: [PATCH 5/7] Add high-level interface to 2013 HESE Aeffs --- icecube_tools/detector/effective_area.py | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/icecube_tools/detector/effective_area.py b/icecube_tools/detector/effective_area.py index 0d4d8bd..0eb46a7 100644 --- a/icecube_tools/detector/effective_area.py +++ b/icecube_tools/detector/effective_area.py @@ -5,6 +5,7 @@ from icecube_tools.utils.data import ( IceCubeData, find_files, + find_folders, data_directory, ) @@ -277,6 +278,10 @@ def get_reader(self, **kwargs): Define an IceCubeAeffReader based on the filename. """ + if R2013_AEFF_FILENAME in self._filename: + + return R2013AeffReader(self._filename, **kwargs) + if R2015_AEFF_FILENAME in self._filename: return R2015AeffReader(self._filename, **kwargs) @@ -372,7 +377,6 @@ def from_dataset(cls, dataset_id, fetch=True, **kwargs): dataset_dir = data_directory - # Find filename if dataset_id == "20181018": files = find_files(dataset_dir, R2018_AEFF_FILENAME) @@ -387,4 +391,11 @@ def from_dataset(cls, dataset_id, fetch=True, **kwargs): # Latest dataset aeff_file_name = files[0] + elif dataset_id == "20131121": + + folders = find_folders(dataset_dir, R2013_AEFF_FILENAME) + + # Folder containing all Aeff info + aeff_file_name = folders[0] + return cls(aeff_file_name, **kwargs) From d0aa2eba18e86fbb94f9e0de0d5aa03a2daa6610 Mon Sep 17 00:00:00 2001 From: Francesca Capel Date: Wed, 21 Apr 2021 18:42:02 +0200 Subject: [PATCH 6/7] Move data location to somewhere sensible and add test for Aeff load --- icecube_tools/utils/data.py | 2 +- tests/test_aeff_load.py | 11 +++++++++++ tests/test_data_interface.py | 2 +- 3 files changed, 13 insertions(+), 2 deletions(-) create mode 100644 tests/test_aeff_load.py diff --git a/icecube_tools/utils/data.py b/icecube_tools/utils/data.py index 56a5961..52ed348 100644 --- a/icecube_tools/utils/data.py +++ b/icecube_tools/utils/data.py @@ -9,7 +9,7 @@ from tqdm import tqdm icecube_data_base_url = "https://icecube.wisc.edu/data-releases" -data_directory = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "data")) +data_directory = os.path.abspath(os.path.join(os.path.expanduser("~"), ".icecube_data")) class IceCubeData: diff --git a/tests/test_aeff_load.py b/tests/test_aeff_load.py new file mode 100644 index 0000000..c024168 --- /dev/null +++ b/tests/test_aeff_load.py @@ -0,0 +1,11 @@ +from icecube_tools.detector.effective_area import ( + EffectiveArea, + _supported_dataset_ids, +) + + +def test_aeff_load(): + + for dataset_id in _supported_dataset_ids: + + _ = EffectiveArea.from_dataset(dataset_id) diff --git a/tests/test_data_interface.py b/tests/test_data_interface.py index ac9b656..8f7ee0a 100644 --- a/tests/test_data_interface.py +++ b/tests/test_data_interface.py @@ -5,7 +5,7 @@ def test_data_scan(): - assert my_data.datasets[0] == "20080911_AMANDA_7_Year_Data.zip" + assert my_data.datasets[1] == "20080911_AMANDA_7_Year_Data.zip" def test_file_download(output_directory): From c6c8f3672619dfbfa5526df166d450d6c1d1b812 Mon Sep 17 00:00:00 2001 From: Francesca Capel Date: Wed, 21 Apr 2021 18:47:51 +0200 Subject: [PATCH 7/7] Add h5py to fix tests --- setup.cfg | 1 + 1 file changed, 1 insertion(+) diff --git a/setup.cfg b/setup.cfg index a61dbec..aa9da19 100644 --- a/setup.cfg +++ b/setup.cfg @@ -28,6 +28,7 @@ install_requires = versioneer vMF pandas + h5py [versioneer] VCS=git