Merge pull request #25 from vuillaut/eff_area

refactoring of the cta irf

Author: vuillaut

ctaplot/__init__.py

@@ -2,3 +2,4 @@
 from .plotcamera import *
 from .ana import *
 from .dataset import *
+from .io import *

ctaplot/ana.py

@@ -20,8 +20,8 @@ def __init__(self):
         self.E = logbin_mean(self.E_bin)
 
         # Area of CTA sites in meters
-        self.ParanalArea = 19.63e6
-        self.LaPalmaArea = 6.61e6
+        self.ParanalArea_prod3 = 19.63e6
+        self.LaPalmaArea_prod3 = 11341149 #6.61e6
 
     def set_E_bin(self, E_bin):
         self.E_bin = E_bin
@@ -355,11 +355,11 @@ def energy_bias(SimuE, RecoE):
 
 
 def get_angles_pipi(angles):
-    return np.mod(angles+pi, 2*pi) - pi
+    return np.mod(angles + np.pi, 2 * np.pi) - np.pi
 
 
 def get_angles_0pi(angles):
-    return np.mod(angles, pi)
+    return np.mod(angles, np.pi)
 
 
 def theta2(RecoAlt, RecoAz, AltSource, AzSource):

ctaplot/io.py

@@ -0,0 +1,69 @@
+import json
+import numpy as np
+
+
+class mc_run_header():
+
+    def __init__(self):
+        pass
+
+    def read(self, filename):
+        """
+        Read a json file containing a Monte-Carlo run header.
+
+        Parameters
+        ----------
+        filename: string
+            path to the json file
+
+        Returns
+        -------
+        dict
+        """
+        with open(filename) as file:
+            self.data = json.load(file)
+
+    def get_e_range(self):
+        """
+        return simulated energy range
+
+        Returns
+        -------
+        [e_min, e_max]
+        """
+        return self.data['E_range']
+
+    def get_core_range(self):
+        """
+        return core range
+
+        Returns
+        -------
+        [r_min, r_max]
+        """
+        return self.data['core_range']
+
+    def get_core_area(self):
+        """
+        return simulated core area
+
+        Returns
+        -------
+        area: float
+        """
+        if self.data['core_pos_mode'] == 1:
+            return np.pi * (self.data['core_range'][1] - self.data['core_range'][0])**2
+        else:
+            raise Exception("Sorry I do not know this simulation mode")
+
+    def get_number_simulated_events(self):
+        """
+        return the number of simulated events (number of simulated showers times number of use)
+
+        Returns
+        -------
+        float
+        """
+        return self.data['n_showers'] * self.data['n_use']
+
+

ctaplot/plots.py

@@ -555,7 +555,7 @@ def plot_effective_area_per_energy(SimuE, RecoE, simuArea, ax=None, **kwargs):
     >>> irf = ctaplot.ana.irf_cta()
     >>> simue = 10**(-2 + 4*np.random.rand(1000))
     >>> recoe = 10**(-2 + 4*np.random.rand(100))
-    >>> ax = ctaplot.plots.plot_effective_area_per_energy(simue, recoe, irf.LaPalmaArea)
+    >>> ax = ctaplot.plots.plot_effective_area_per_energy(simue, recoe, irf.LaPalmaArea_prod3)
     """
 
     ax = plt.gca() if ax is None else ax
@@ -576,7 +576,7 @@ def plot_effective_area_per_energy(SimuE, RecoE, simuArea, ax=None, **kwargs):
     return ax
 
 
-def plot_effective_area_requirement(cta_site, ax=None, **kwargs):
+def plot_effective_area_cta_requirements(cta_site, ax=None, **kwargs):
     """
     Plot the CTA requirement for the effective area
 
@@ -609,7 +609,7 @@ def plot_effective_area_requirement(cta_site, ax=None, **kwargs):
     return ax
 
 
-def plot_effective_area_performances(cta_site, ax=None, **kwargs):
+def plot_effective_area_cta_performances(cta_site, ax=None, **kwargs):
     """
     Plot the CTA performances for the effective area
 
@@ -667,15 +667,15 @@ def saveplot_effective_area_per_energy(SimuE, RecoE, simuArea, ax=None, Outfile=
     ax = plot_effective_area_per_energy(SimuE, RecoE, simuArea, ax=ax, **kwargs)
 
     if cta_site:
-        ax = plot_effective_area_requirement(cta_site, ax=ax, color='black')
+        ax = plot_effective_area_cta_requirements(cta_site, ax=ax, color='black')
 
     plt.savefig(Outfile + ".png", bbox_inches="tight", format='png', dpi=200)
     plt.close()
 
     return ax
 
 
-def plot_sensitivity_requirement(cta_site, ax=None, **kwargs):
+def plot_sensitivity_cta_requirements(cta_site, ax=None, **kwargs):
     """
     Plot the CTA requirement for the sensitivity
     Parameters
@@ -707,7 +707,7 @@ def plot_sensitivity_requirement(cta_site, ax=None, **kwargs):
     return ax
 
 
-def plot_sensitivity_performances(cta_site, ax=None, **kwargs):
+def plot_sensitivity_cta_performances(cta_site, ax=None, **kwargs):
     """
     Plot the CTA performances for the sensitivity
     Parameters
@@ -900,15 +900,15 @@ def saveplot_angular_res_per_energy(RecoAlt, RecoAz, AltSource, AzSource, SimuE,
     ax = plot_angular_res_per_energy(RecoAlt, RecoAz, AltSource, AzSource, SimuE, ax=ax, **kwargs)
 
     if cta_site:
-        ax = plot_angular_res_requirement(cta_site, ax=ax)
+        ax = plot_angular_res_requirements(cta_site, ax=ax)
 
     plt.savefig(Outfile + ".png", bbox_inches="tight", format='png', dpi=200);
     plt.close()
 
     return ax
 
 
-def plot_angular_res_requirement(cta_site, ax=None, **kwargs):
+def plot_angular_res_cta_requirements(cta_site, ax=None, **kwargs):
     """
     Plot the CTA requirement for the angular resolution
     Parameters
@@ -1244,7 +1244,7 @@ def plot_energy_resolution(SimuE, RecoE, ax=None, bias_correction=False, **kwarg
     return ax
 
 
-def plot_energy_resolution_requirements(cta_site, ax=None, **kwargs):
+def plot_energy_resolution_cta_requirements(cta_site, ax=None, **kwargs):
     """
     Plot the cta requirement for the energy resolution
 
@@ -1321,7 +1321,7 @@ def saveplot_energy_resolution(SimuE, RecoE, Outfile="EnergyResolution.png", cta
     ax = plot_energy_resolution(SimuE, RecoE)
 
     if cta_site != None:
-        ax = plot_energy_resolution_requirements(cta_site, ax=ax)
+        ax = plot_energy_resolution_cta_requirements(cta_site, ax=ax)
 
     plt.savefig(Outfile, bbox_inches="tight", format='png', dpi=200)
     plt.close()
@@ -1600,11 +1600,11 @@ def set_site(self, site):
         if site == 'north' or site == 'lapalma':
             self.site = site
             irf = ana.irf_cta()
-            self.SimuArea = irf.LaPalmaArea
+            self.SimuArea = irf.LaPalmaArea_prod3
         if site == 'south' or site == 'paranal':
             self.site = site
             irf = ana.irf_cta()
-            self.SimuArea = irf.ParanalArea
+            self.SimuArea = irf.ParanalArea_prod3
 
     def plot_all(self, MultiplicityMin=[2, 4]):
         self.multiplicity_hist()
@@ -1700,7 +1700,7 @@ def angular_res_per_energy(self, MultiplicityMin=[2], cta_goal=True, ax=None, **
         ax = plt.gca() if ax is None else ax
 
         if cta_goal:
-            ax = plot_angular_res_requirement(self.site, ax=ax, color='black')
+            ax = plot_angular_res_cta_requirements(self.site, ax=ax, color='black')
 
         if len(MultiplicityMin) > 1 and not 'alpha' in kwargs:
             kwargs['alpha'] = 0.8
@@ -1767,7 +1767,7 @@ def effective_area_per_energy(self, MultiplicityMin=[2], cta_goal=True, ax=None,
         ax = plt.gca() if ax is None else ax
 
         if cta_goal:
-            ax = plot_effective_area_requirement(self.site, ax=ax, color='black')
+            ax = plot_effective_area_cta_requirements(self.site, ax=ax, color='black')
 
         if len(MultiplicityMin) > 1 and not 'alpha' in kwargs:
             kwargs['alpha'] = 0.8
@@ -1868,7 +1868,7 @@ def energy_resolution(self, cta_goal=True, MultiplicityMin=[2], bias_correction=
         ax = plt.gca()
 
         if cta_goal:
-            ax = plot_energy_resolution_requirements(self.site, ax=ax, color='black')
+            ax = plot_energy_resolution_cta_requirements(self.site, ax=ax, color='black')
 
         alpha = 1.0
         if len(MultiplicityMin) > 1:

ctaplot/tests/test_plots.py

@@ -17,9 +17,9 @@ def test_plot_energy_resolution():
     plots.plot_energy_resolution(E, E**2, color='red')
 
 
-def test_plot_energy_resolution_requirements():
+def test_plot_energy_resolution_cta_requirements():
     plt.close('all')
-    plots.plot_energy_resolution_requirements('north', color='green')
+    plots.plot_energy_resolution_cta_requirements('north', color='green')
 
 
 def test_saveplot_energy_resolution():
@@ -39,24 +39,24 @@ def test_plot_angular_resolution_cta_performances():
     plots.plot_angular_res_cta_performance('north', color='green')
 
 
-def test_plot_angular_resolution_requirements():
-    plots.plot_angular_res_requirement('north', color='green')
+def test_plot_angular_resolution_cta_requirements():
+    plots.plot_angular_res_cta_requirements('north', color='green')
 
 
-def test_plot_effective_area_performances():
-    plots.plot_effective_area_performances('north', color='green')
+def test_plot_effective_area_cta_performances():
+    plots.plot_effective_area_cta_performances('north', color='green')
 
 
-def test_plot_effective_area_requirement():
-    plots.plot_effective_area_requirement('north', color='green')
+def test_plot_effective_area_cta_requirements():
+    plots.plot_effective_area_cta_requirements('north', color='green')
 
 
-def test_plot_sensitivity_performances():
-    plots.plot_sensitivity_performances('north', color='green')
+def test_plot_sensitivity_cta_performances():
+    plots.plot_sensitivity_cta_performances('north', color='green')
 
 
-def test_plot_sensitivity_requirement():
-    plots.plot_sensitivity_requirement('north', color='green')
+def test_plot_sensitivity_cta_requirements():
+    plots.plot_sensitivity_cta_requirements('north', color='green')
 
 def test_plot_theta2():
     n = 10

setup.py

@@ -18,7 +18,7 @@ def package_files(directory):
 print("dataset {}".format(dataset))
 
 setup(name='ctaplot',
-      version=0.1,
+      version=0.2,
       description="DESCRIPTION",
       # these should be minimum list of what is needed to run (note
       # don't need to list the sub-dependencies like numpy, since

share/simu_mc/gamma_20deg_0deg_cta-prod3-LaPalma_run_header.json

@@ -0,0 +1,30 @@
+{"shower_prog_id": 1,
+ "shower_prog_vers": 6990,
+ "shower_prog_start": 1482580800,
+ "detector_prog_id": 1,
+ "detector_prog_vers": 1481293431,
+ "detector_prog_start": 1482579523,
+ "obsheight": 2147.0,
+ "n_showers": 25000,
+ "n_use": 20,
+ "core_pos_mode": 1,
+ "core_range": [0.0, 1900.0],
+ "alt_range": [1.2217305, 1.2217305],
+ "az_range": [0.0, 0.0],
+ "diffuse": 1,
+ "viewcone": [ 0.0, 10.0],
+ "E_range": [3.0e-03, 3.3e+02],
+ "spectral_index": -2.0,
+ "B_total": 38.7280158996582,
+ "B_inclination": 0.662263810634613,
+ "B_declination": 0.0,
+ "injection_height": -1.0,
+ "atmosphere": 36,
+ "corsika_iact_options": 187,
+ "corsika_low_E_model": 2,
+ "corsika_high_E_model": 3,
+ "corsika_bunchsize": 5.0,
+ "corsika_wlen_min": 240.0,
+ "corsika_wlen_max": 700.0,
+ "corsika_low_E_detail": 0,
+ "corsika_high_E_detail": 101}

share/simu_mc/gamma_20deg_0deg_cta-prod3_Paranal_run_header.json

@@ -0,0 +1,30 @@
+{   "shower_prog_vers": 6990,
+ "shower_prog_start": 1467288000,
+ "detector_prog_id": 1,
+ "detector_prog_vers": 1462392225,
+ "detector_prog_start": 1467287384,
+ "obsheight": 2150.0,
+ "n_showers": 20000,
+ "n_use": 10,
+ "core_pos_mode": 1,
+ "core_range": [   0.0, 2500.0],
+ "alt_range": [1.2217305, 1.2217305],
+ "az_range": [0.0, 0.0],
+ "diffuse": 0,
+ "viewcone": [0.0, 0.0],
+ "E_range": [3.0e-03, 3.3e+02],
+ "spectral_index": -2.0,
+ "B_total": 23.11772346496582,
+ "B_inclination": -0.39641156792640686,
+ "B_declination": 0.0,
+ "injection_height": -1.0,
+ "atmosphere": 26,
+ "corsika_iact_options": 187,
+ "corsika_low_E_model": 2,
+ "corsika_high_E_model": 3,
+ "corsika_bunchsize": 5.0,
+ "corsika_wlen_min": 240.0,
+ "corsika_wlen_max": 700.0,
+ "corsika_low_E_detail": 0,
+ "corsika_high_E_detail": 303
+}