Clean some functions

vuillaut · Jul 17, 2017 · e093c95 · e093c95
1 parent 5651712
commit e093c95
Show file tree

Hide file tree

Showing 4 changed files with 118 additions and 146 deletions.
diff --git a/notebooks/example.ipynb b/notebooks/example.ipynb
diff --git a/pschitt/camera_image.py b/pschitt/camera_image.py
@@ -39,10 +39,7 @@ def find_closest_pixel(pos, pixel_tab):
     -------
     Pixel index in the given pixel_tab, corresponding distance to the pixel
     """
-    #D = np.linalg.norm(pixel_tab-pos, axis=1)
-    #linalg.norm is surprisingly slow
     x = pixel_tab - pos
-    # D = np.sqrt(x[:, 0]**2 + x[:, 1]**2)
     D2 = x[:, 0] ** 2 + x[:, 1] ** 2
     return D2.argmin()
 
@@ -65,24 +62,20 @@ def photons_to_signal(photon_pos_tab, pixel_tab):
     d_max2 = (pixel_tab[:, 0]**2 + pixel_tab[:, 1]**2).max()
     for photon in photon_pos_tab:
         if photon[0]**2 + photon[1]**2 < d_max2:
-            # pxi = find_closest_pixel(photon, pixel_tab)
             x = pixel_tab - photon
             D2 = x[:, 0] ** 2 + x[:, 1] ** 2
-            # pxi = D2.argmin()
             count[D2.argmin()] += 1
 
-    # for photon in photon_pos_tab:
-    #     pxi = find_closest_pixel(photon, pixel_tab)
-    #     count[pxi] += photon[0]**2 + photon[1]**2 < d_max2
-
     return count
 
 
 def write_camera_image(pix_hist, filename="data/camera_image.txt"):
     """
     Save camera image in a file
-    :param pix_hist: array - pixel histogram
-    :param filename: string
+    Parameters
+    ----------
+    pix_hist : array - pixels signal
+    filename : string
     """
     np.savetxt(filename,pix_hist,fmt='%.5f')
 
@@ -103,9 +96,15 @@ def shower_image_in_camera_old(telescope, photon_pos_tab, pixel_pos_filename):
 def add_noise_poisson(signal, lam=100):
     """
     Add Poisson noise to the image
-    :param pix_hist: pixel histogram
-    :param lam: lambda for Poisson law
-    :return: pixel histogram
+
+    Parameters
+    ----------
+    signal : pixels signal - array
+    lam : lambda for Poisson law - float
+
+    Returns
+    -------
+    signal array
     """
     if(lam>=0):
         signal += np.random.poisson(lam, signal.size)
@@ -199,26 +198,3 @@ def array_shower_imaging(shower, tel_array, noise):
         shower_camera_image(shower, tel, noise)
 
 
-def shower_camera_image_argorder(shower, noise, tel):
-    shower_camera_image(shower, tel, noise)
-    return None
-
-
-def array_shower_imaging_multiproc(shower, tel_array, noise):
-    """
-    TEST
-    Given a shower object and an array of telescopes, compute the image of the shower in each camera
-    Background noise can be added
-    The camera signal is registered in all tel.signal_hist
-    Parameters
-    ----------
-    shower: 3D Numpy array with a list of space position points composing the shower
-    tel_array: Numpy array or list of telescope classes
-    noise: float
-    """
-
-    pool = Pool(processes=4)
-    func = partial(shower_camera_image_argorder, shower, noise)
-    pool.map(func, tel_array)
-    pool.close()
-    pool.join()
diff --git a/pschitt/hillas.py b/pschitt/hillas.py
@@ -2,8 +2,6 @@
 
 
 import numpy as np
-from astropy.units import Quantity
-import astropy.units as u
 from . import geometry as geo
 from math import *
 
@@ -14,12 +12,9 @@ def hillas_parameters(pix_x, pix_y, image):
     Hillas parameters calculation.
     Parameters
     ----------
-    pix_x : array_like
-        Pixel x-coordinate
-    pix_y : array_like
-        Pixel y-coordinate
-    image : array_like
-        Pixel values corresponding
+    pix_x : 1D Numpy array - Pixel x-coordinate
+    pix_y : 1D Numpy array - Pixel y-coordinate
+    image : 1D Numpy array - signal in each pixel
     Returns
     -------
     hillas_parameters - 1D Numpy array
@@ -28,14 +23,6 @@ def hillas_parameters(pix_x, pix_y, image):
     amplitude = image.sum()
     assert amplitude > 0
 
-
-    # pix_x = Quantity(np.asanyarray(pix_x, dtype=np.float64)).value
-    # pix_y = Quantity(np.asanyarray(pix_y, dtype=np.float64)).value
-
-    # assert pix_x.shape == image.shape
-    # assert pix_y.shape == image.shape
-
-
     momdata = np.row_stack([pix_x,
                             pix_y,
                             pix_x * pix_x,
@@ -157,6 +144,19 @@ def reco_impact_parameter_test(hillas_parameters_1, tel1, hillas_parameters_2, t
 
 
 def impact_parameter_average(alltel, HillasParameters):
+    """
+    Reconstruction of the impact parameter taking the directions intersection
+    of telescopes two by two and averaging the result
+
+    Parameters
+    ----------
+    alltel: list of telescopes
+    HillasParameters: list of hillas_parameters
+
+    Returns
+    -------
+    Space coordinates of the impact parameter - list of floats
+    """
     P = []
     for i in np.arange(len(alltel)-1):
         for j in np.arange(i+1, len(alltel)):
@@ -171,6 +171,17 @@ def impact_parameter_average(alltel, HillasParameters):
 
 
 def coef_hillas_ponderation(hillas_parameters_1, hillas_parameters_2):
+    """
+    Weight of two telescopes reconstruction
+    Parameters
+    ----------
+    hillas_parameters_1: list of hillas parameters for camera 1
+    hillas_parameters_2: list of hillas parameters for camera 2
+
+    Returns
+    -------
+    Float - weight
+    """
     phi1 = hillas_parameters_1[6]
     intensity1 = hillas_parameters_1[0]
     width1 = hillas_parameters_1[4]
@@ -184,14 +195,26 @@ def coef_hillas_ponderation(hillas_parameters_1, hillas_parameters_2):
     gx2 = hillas_parameters_2[1]
     gy2 = hillas_parameters_2[2]
 
-    r = (gx1*gx1 + gy1*gy1 + gx2*gx2 + gy2*gy2)**6;
-    #print(gx1,gy1,gx2,gy2)
+    r = (gx1*gx1 + gy1*gy1 + gx2*gx2 + gy2*gy2)**6
+
     if r < 1e-6:
-        return fabs(sin(phi1 - phi2)) / ((1.0/intensity1 + 1.0/intensity2) + (width1/length1 + width2/length2));
+        return fabs(sin(phi1 - phi2)) / ((1.0/intensity1 + 1.0/intensity2) + (width1/length1 + width2/length2))
     else:
-        return fabs(sin(phi1 - phi2)) / ( (1.0/intensity1 + 1.0/intensity2 + width1/length1 + width2/length2) * r);
+        return fabs(sin(phi1 - phi2)) / ( (1.0/intensity1 + 1.0/intensity2 + width1/length1 + width2/length2) * r)
+
 
 def impact_parameter_ponderated(alltel, HillasParameters):
+    """
+    Reconstruction of the impact parameter using the weigths from coef_hillas_ponderation
+    Parameters
+    ----------
+    alltel: list of telescopes
+    HillasParameters: list of hillas_parameters for each telescope
+
+    Returns
+    -------
+    Coordinates of the impact parameter - list of floats
+    """
     P = []
     C = []
     for i in np.arange(len(alltel)):
@@ -202,16 +225,10 @@ def impact_parameter_ponderated(alltel, HillasParameters):
                 p = reco_impact_parameter_test(hp1, alltel[i], hp2, alltel[j], 0)
                 c = coef_hillas_ponderation(hp1, hp2)
                 P.append(p)
-                r12 = sqrt(hp1[1]**2 + hp1[2]**2)
-                r22 = sqrt(hp2[1]**2 + hp2[2]**2)
-                #if(r12>0.9*alltel[i].camera_size or r22>0.9*alltel[j].camera_size):
-                #    c = 0
                 C.append(c)
 
     P = np.array(P)
     C = np.array(C)
-    #if(C.min()!=C.max()):
-    #    C = (C - C.min())/(C.max()-C.min())
 
     return [np.average(P[:,0],weights=C), np.average(P[:,1],weights=C), np.average(P[:,2],weights=C)]
 

diff --git a/pschitt/vizualisation.py b/pschitt/vizualisation.py
@@ -6,14 +6,17 @@
 from scipy import stats
 
 
-def plot_shower3d(shower, alltel, density_color=True):
+def plot_shower3d(shower, alltel, **options):
     """
     Display the sky object (shower) and the telescope in a 3D representation
     Parameters
     ----------
     shower: array of points (arrays [x,y,z])
     alltel: array of telescopes (telescope class)
-    density_color: boolean to use density for particles color. Set False to speed-up plotting.
+    options:
+        - density_color = True: use density for particles color. False by default.
+        - display = True: show the plot. False by default
+        - outfile = "file.eps" : save the plot as `file.eps`. False by default.
     """
     fig = plt.figure(figsize=(12,12))
     ax = fig.add_subplot(111, projection='3d')
@@ -28,7 +31,7 @@ def plot_shower3d(shower, alltel, density_color=True):
 
     values = shower.array.T
 
-    if density_color:
+    if opt.get("density_color") == True:
         kde = stats.gaussian_kde(values)
         density = kde(values)
         ax.scatter(values[0] , values[1], values[2], marker='o', c=density)
@@ -39,8 +42,12 @@ def plot_shower3d(shower, alltel, density_color=True):
     ax.set_xlabel("[m]")
     ax.set_zlabel("altitude [m]")
     plt.legend()
-    #plt.axis('equal')
-    #plt.show()
+    if opt.get("display") == True:
+        plt.show()
+    if opt.get("outfile"):
+        outfile = opt.get("outfile")
+        assert isinstance(outfile, str), "The given outfile option should be a string"
+        plt.savefig(outfile + '.eps', format='eps', dpi=200)
 
 
 def display_camera_image(telescope):
@@ -54,6 +61,7 @@ def display_camera_image(telescope):
     plt.scatter(telescope.pixel_tab[:, 0], telescope.pixel_tab[:, 1], c=telescope.signal_hist)
     plt.axis('equal')
     plt.colorbar(label='counts')
+    plt.show()
 
 
 def display_stacked_cameras(telescope_array):
@@ -73,6 +81,7 @@ def display_stacked_cameras(telescope_array):
     plt.scatter(tel0.pixel_tab[:, 0], tel0.pixel_tab[:, 1], c=stacked_hist)
     plt.axis('equal')
     plt.colorbar(label='counts')
+    plt.show()
 
 
 def display_pointing_tel(tel, show=True):