Merge pull request #463 from kosack/master

more cleanups for instrument class usage

ctapipe/image/muon/muon_diagnostic_plots.py

@@ -43,7 +43,6 @@ def plot_muon_efficiency(outputpath):
 
     print('Gaussian fit with mu=', mu, 'sigma=', sigma)
 
-    #ax = figeff.add_subplot(1,3,1)
     conteff = ax.hist(t['MuonEff'], nbins)
     ax.set_xlim(0.2 * min(t['MuonEff']), 1.2 * max(t['MuonEff']))
 
@@ -54,35 +53,31 @@ def plot_muon_efficiency(outputpath):
 
     ax.set_ylim(0., 1.2 * max(conteff[0]))
     ax.set_xlabel('Muon Efficiency')
-    # plt.figure(fig.number)
     plt.draw()
 
-    #axip = figeff.add_subplot(1,3,1)
     yimp = np.linspace(min(t['ImpactP']), max(t['ImpactP']), nbins)
     contimp = axip.hist(t['ImpactP'], nbins)
     axip.set_xlim(0.2 * min(t['ImpactP']), 1.2 * max(t['ImpactP']))
     axip.set_ylim(0., 1.2 * max(contimp[0]))
     axip.set_xlabel('Impact Parameter (m)')
-    # plt.figure(figip.number)
+
     plt.draw()
 
     heffimp = Histogram(nbins=[16, 16],
                         ranges=[(min(t['MuonEff']), max(t['MuonEff'])),
                                 (min(t['ImpactP']), max(t['ImpactP']))])  # ,axisNames=["MuonEfficiency","ImpactParameter"])
-    # embed()
-    # heffimp.bin_lower_edges([xtest,yimp])
+
     heffimp.fill([t['MuonEff'], t['ImpactP']])
     heffimp.draw_2d()
 
-    #axrw = figeff.add_subplot(1,3,1)
     yrw = np.linspace(min(t['RingWidth']), max(t['RingWidth']), nbins)
     contrw = axrw.hist(t['RingWidth'], nbins)
     axrw.set_xlim(0.2 * min(t['RingWidth']), 1.2 * max(t['RingWidth']))
     axrw.set_ylim(0., 1.2 * max(contrw[0]))
     axrw.set_xlabel('Ring Width ($^\circ$)')
-    # plt.figure(figrw.number)
+
     plt.draw()
-    # plt.show()
+
     if outputpath is not None:
         print("saving figure at", outputpath)
         fig.savefig(str(outputpath) + '_MuonEff.png')
@@ -93,16 +88,12 @@ def plot_muon_efficiency(outputpath):
         plt.show()
 
 
-def plot_muon_event(event, muonparams, geom_dict=None, args=None):
+def plot_muon_event(event, muonparams, args=None):
 
     if muonparams[0] is not None:
 
         # Plot the muon event and overlay muon parameters
         fig = plt.figure(figsize=(16, 7))
-        # if args.display:
-        #    plt.show(block=False)
-        #pp = PdfPages(args.output_path) if args.output_path is not None else None
-        # pp = None #For now, need to correct this
 
         colorbar = None
         colorbar2 = None
@@ -112,34 +103,21 @@ def plot_muon_event(event, muonparams, geom_dict=None, args=None):
             # Only create two pads if there is timing information extracted
             # from the calibration
             ax1 = fig.add_subplot(1, npads, 1)
-            plotter = CameraPlotter(event, geom_dict)
-            #image = event.dl1.tel[tel_id].calibrated_image
+            plotter = CameraPlotter(event)
             image = event.dl1.tel[tel_id].image[0]
-            # Get geometry
-            geom = None
-            if geom_dict is not None and tel_id in geom_dict:
-                geom = geom_dict[tel_id]
-            else:
-                #log.debug("[calib] Guessing camera geometry")
-                geom = CameraGeometry.guess(*event.inst.pixel_pos[tel_id],
-                                            event.inst.optical_foclen[tel_id])
-                #log.debug("[calib] Camera geometry found")
-                if geom_dict is not None:
-                    geom_dict[tel_id] = geom
+            geom = event.inst.subarray.tel[tel_id].camera
 
             tailcuts = (5., 7.)
             # Try a higher threshold for
             if geom.cam_id == 'FlashCam':
                 tailcuts = (10., 12.)
 
-            #print("Using Tail Cuts:",tailcuts)
             clean_mask = tailcuts_clean(geom, image,
                                         picture_thresh=tailcuts[0],
                                         boundary_thresh=tailcuts[1])
 
             signals = image * clean_mask
 
-            #print("Ring Centre in Nominal Coords:",muonparams[0].ring_center_x,muonparams[0].ring_center_y)
             muon_incl = np.sqrt(muonparams[0].ring_center_x**2. +
                                 muonparams[0].ring_center_y**2.)
 
@@ -173,10 +151,6 @@ def plot_muon_event(event, muonparams, geom_dict=None, args=None):
             ringrad_camcoord = muonparams[0].ring_radius.to(u.rad) \
                                * event.inst.optical_foclen[tel_id] * 2.  # But not FC?
 
-            #rot_angle = 0.*u.deg
-            # if event.inst.optical_foclen[tel_id] > 10.*u.m and event.dl0.tel[tel_id].num_pixels != 1764:
-            #rot_angle = -100.14*u.deg
-
             px, py = event.inst.pixel_pos[tel_id]
             flen = event.inst.optical_foclen[tel_id]
             camera_coord = CameraFrame(x=px, y=py,
@@ -188,7 +162,6 @@ def plot_muon_event(event, muonparams, geom_dict=None, args=None):
                 NominalFrame(array_direction=altaz,
                              pointing_direction=altaz)
             )
-            #,focal_length = event.inst.optical_foclen[tel_id])) # tel['TelescopeTable_VersionFeb2016'][tel['TelescopeTable_VersionFeb2016']['TelID']==telid]['FL'][0]*u.m))
 
             px = nom_coord.x.to(u.deg)
             py = nom_coord.y.to(u.deg)
@@ -224,8 +197,6 @@ def plot_muon_event(event, muonparams, geom_dict=None, args=None):
             camera1.add_ellipse(centroid, ringrad_camcoord.value,
                                 ringrad_camcoord.value, 0., 0., color="red")
 
-#            ax1.set_title("CT {} ({}) - Mean pixel charge"
-#                          .format(tel_id, geom_dict[tel_id].cam_id))
 
             if muonparams[1] is not None:
                 # continue #Comment this...(should ringwidthfrac also be *0.5?)
@@ -271,9 +242,7 @@ def plot_muon_event(event, muonparams, geom_dict=None, args=None):
                 camera2.update(True)
                 plt.pause(1.)  # make shorter
 
-            # plt.pause(0.1)
-            # if pp is not None:
-            #    pp.savefig(fig)
+
             fig.savefig(str(args.output_path) + "_" +
                         str(event.dl0.event_id) + '.png')
 

ctapipe/image/tests/test_charge_extraction.py

@@ -1,12 +1,11 @@
-from ctapipe.io.hessio import hessio_event_source
-from ctapipe.utils import get_dataset
-from ctapipe.instrument import CameraGeometry
 import numpy as np
 from numpy.testing import assert_almost_equal
 
 from ctapipe.image.charge_extractors import FullIntegrator, \
     SimpleIntegrator, GlobalPeakIntegrator, LocalPeakIntegrator, \
     NeighbourPeakIntegrator, ChargeExtractorFactory, AverageWfPeakIntegrator
+from ctapipe.io.hessio import hessio_event_source
+from ctapipe.utils import get_dataset
 
 
 def get_test_event():
@@ -97,8 +96,7 @@ def test_nb_peak_integration():
     ped = event.mc.tel[telid].pedestal
     data_ped = data - np.atleast_3d(ped/nsamples)
     data_ped = np.array([data_ped[0], data_ped[0]])  # Test LG functionality
-    geom = CameraGeometry.guess(*event.inst.pixel_pos[telid],
-                                event.inst.optical_foclen[telid])
+    geom = event.inst.subarray.tel[telid].camera
     nei = geom.neighbors
 
     integrator = NeighbourPeakIntegrator(None, None)

ctapipe/instrument/__init__.py

@@ -1,10 +1,10 @@
-from .camera import CameraGeometry, get_camera_types, print_camera_types
+from .camera import CameraGeometry
 from .atmosphere import get_atmosphere_profile_table, get_atmosphere_profile_functions
 from .telescope import TelescopeDescription
 from .optics import OpticsDescription
 from .subarray import SubarrayDescription
 
 
-__all__ = ['CameraGeometry', 'get_camera_types','print_camera_types',
+__all__ = ['CameraGeometry',
            'get_atmosphere_profile_functions', 'TelescopeDescription',
            'OpticsDescription','SubarrayDescription']

ctapipe/instrument/camera.py

@@ -15,9 +15,7 @@
 from ctapipe.utils import get_dataset, find_all_matching_datasets
 from ctapipe.utils.linalg import rotation_matrix_2d
 
-__all__ = ['CameraGeometry',
-           'get_camera_types',
-           'print_camera_types']
+__all__ = ['CameraGeometry',]
 
 logger = logging.getLogger(__name__)
 
@@ -63,37 +61,35 @@ class CameraGeometry:
     of the data. Note that this function is memoized, so calling it multiple 
     times with the same inputs will give back the same object (for speed).
 
+    Parameters
+    ----------
+    self: type
+        description
+    cam_id: camera id name or number
+        camera identification string
+    pix_id: array(int)
+        pixels id numbers
+    pix_x: array with units
+        position of each pixel (x-coordinate)
+    pix_y: array with units
+        position of each pixel (y-coordinate)
+    pix_area: array(float)
+        surface area of each pixel, if None will be calculated
+    neighbors: list(arrays)
+        adjacency list for each pixel
+    pix_type: string
+        either 'rectangular' or 'hexagonal'
+    pix_rotation: value convertable to an `astropy.coordinates.Angle`
+        rotation angle with unit (e.g. 12 * u.deg), or "12d"
+    cam_rotation: overall camera rotation with units
     """
 
     _geometry_cache = {}  # dictionary CameraGeometry instances for speed
 
     def __init__(self, cam_id, pix_id, pix_x, pix_y, pix_area, pix_type,
                  pix_rotation="0d", cam_rotation="0d",
                  neighbors=None, apply_derotation=True):
-        """
-        Parameters
-        ----------
-        self: type
-            description
-        cam_id: camera id name or number
-            camera identification string
-        pix_id: array(int)
-            pixels id numbers
-        pix_x: array with units
-            position of each pixel (x-coordinate)
-        pix_y: array with units
-            position of each pixel (y-coordinate)
-        pix_area: array(float)
-            surface area of each pixel, if None will be calculated
-        neighbors: list(arrays)
-            adjacency list for each pixel
-        pix_type: string
-            either 'rectangular' or 'hexagonal'
-        pix_rotation: value convertable to an `astropy.coordinates.Angle`
-            rotation angle with unit (e.g. 12 * u.deg), or "12d"
-        cam_rotation: overall camera rotation with units
 
-        """
         self.cam_id = cam_id
         self.pix_id = pix_id
         self.pix_x = pix_x
@@ -480,47 +476,4 @@ def _neighbor_list_to_matrix(neighbors):
     return neigh2d
 
 
-def get_camera_types(inst):
-    """ return dict of camera names mapped to a list of tel_ids
-     that use that camera
-
-     Parameters
-     ----------
-     inst: instument Container
-
-     """
-
-    cam_types = defaultdict(list)
-
-    for telid in inst.pixel_pos:
-        x, y = inst.pixel_pos[telid]
-        f = inst.optical_foclen[telid]
-        geom = CameraGeometry.guess(x, y, f)
-
-        cam_types[geom.cam_id].append(telid)
-
-    return cam_types
-
-
-def print_camera_types(inst, printer=print):
-    """
-    Print out a friendly table of which camera types are registered in the 
-    inst dictionary (from a hessio file), along with their starting and 
-    stopping tel_ids.
-
-    Parameters
-    ----------
-    inst: ctapipe.io.containers.InstrumentContainer
-        input container
-    printer: func
-        function to call to output the text (default is the standard python 
-        print command, but you can give for example logger.info to have it 
-        write to a logger) 
-    """
-    camtypes = get_camera_types(inst)
 
-    printer("              CAMERA  Num IDmin  IDmax")
-    printer("=====================================")
-    for cam, tels in camtypes.items():
-        printer("{:>20s} {:4d} {:4d} ..{:4d}".format(cam, len(tels), min(tels),
-                                                     max(tels)))

ctapipe/instrument/subarray.py

@@ -22,11 +22,20 @@ class SubarrayDescription:
         dict of telescope positions by tel_id
     tel_descriptions: dict(TelescopeDescription)
         array of TelescopeDescriptions by tel_id
+
+    Attributes
+    ----------
+    name
+       name of subarray
+    positions
+       x,y position of each telescope as length-2 arrays of unit quantities
+    tels
+       dict of TelescopeDescription for each telescope in the subarray
     """
 
     def __init__(self, name, tel_positions=None, tel_descriptions=None):
 
-        self.name = name
+        self.name = name #: name of telescope
         self.positions = tel_positions or dict()
         self.tels = tel_descriptions or dict()
 

ctapipe/instrument/telescope.py

@@ -42,9 +42,18 @@ def __init__(self,
                  optics: OpticsDescription,
                  camera: CameraGeometry):
 
-        self.optics = optics
-        self.camera = camera
-
+        self._optics = optics
+        self._camera = camera
+
+    @property
+    def optics(self):
+        """ OpticsDescription for this telescope """
+        return self._optics
+
+    @property
+    def camera(self):
+        """ CameraGeometry for this telescope"""
+        return self._camera
 
     @classmethod
     def guess(cls, pix_x, pix_y, effective_focal_length):

ctapipe/plotting/camera.py

@@ -4,7 +4,6 @@
 
 from matplotlib import pyplot as plt
 
-from ctapipe.instrument import CameraGeometry
 from ctapipe.visualization import CameraDisplay
 
 from astropy import units as u
@@ -21,35 +20,20 @@ class CameraPlotter:
     ----------
     event : container
         A `ctapipe` event container
-    geom_dict : dict
-        Dictionary to store the geometries of cameras
     """
 
-    def __init__(self, event, geom_dict=None):
+    def __init__(self, event):
         """
         Parameters
         ----------
         event : container
             A `ctapipe` event container
-        geom_dict : dict
-            A pre-build geom_dict, or an empty dict to store any geoms
-            calculated
-            dict[(num_pixels, focal_length)] = 
-            `ctapipe.instrument.CameraGeometry`
         """
         self.event = event
-        self.geom_dict = {} if geom_dict is None else geom_dict
         self.cameradisplay_dict = {}
 
     def get_geometry(self, tel):
-        npix = len(self.event.r0.tel[tel].adc_sums[0])
-        cam_dimensions = (npix, self.event.inst.optical_foclen[tel])
-
-        if tel not in self.geom_dict:
-            self.geom_dict[tel] = \
-                CameraGeometry.guess(*self.event.inst.pixel_pos[tel],
-                                     self.event.inst.optical_foclen[tel])
-        return self.geom_dict[tel]
+        return self.event.inst.subarray.tel[tel].camera
 
     def draw_camera(self, tel, data, axes=None):
         """