Add features

icecube_tools/detector/effective_area.py

@@ -1,9 +1,11 @@
 import numpy as np
+import os
 from abc import ABC, abstractmethod
 
 from icecube_tools.utils.data import (
     IceCubeData,
     find_files,
+    find_folders,
     data_directory,
 )
 
@@ -12,11 +14,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 +63,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.linspace(-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.
@@ -197,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)
@@ -292,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)
@@ -307,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)

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.
         """
 

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:
@@ -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

setup.cfg

@@ -28,6 +28,7 @@ install_requires =
     versioneer
     vMF
     pandas
+    h5py
 
 [versioneer]
 VCS=git

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)

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):