diff --git a/ctapipe/image/muon/muon_diagnostic_plots.py b/ctapipe/image/muon/muon_diagnostic_plots.py index 28f0080c2e5..cfaa928c0d9 100644 --- a/ctapipe/image/muon/muon_diagnostic_plots.py +++ b/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') diff --git a/ctapipe/image/tests/test_charge_extraction.py b/ctapipe/image/tests/test_charge_extraction.py index 5255efe5b96..4e4a703b45c 100644 --- a/ctapipe/image/tests/test_charge_extraction.py +++ b/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) diff --git a/ctapipe/instrument/__init__.py b/ctapipe/instrument/__init__.py index f9984235a59..ef1c33ef71b 100644 --- a/ctapipe/instrument/__init__.py +++ b/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'] \ No newline at end of file diff --git a/ctapipe/instrument/camera.py b/ctapipe/instrument/camera.py index 174f4e7fc7e..42481354040 100644 --- a/ctapipe/instrument/camera.py +++ b/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,6 +61,27 @@ 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 @@ -70,30 +89,7 @@ class CameraGeometry: 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))) diff --git a/ctapipe/instrument/subarray.py b/ctapipe/instrument/subarray.py index 137b1cddd8d..2b9e27fc549 100644 --- a/ctapipe/instrument/subarray.py +++ b/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() diff --git a/ctapipe/instrument/telescope.py b/ctapipe/instrument/telescope.py index 13d969d4d3f..f67e1c5977e 100644 --- a/ctapipe/instrument/telescope.py +++ b/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): diff --git a/ctapipe/plotting/camera.py b/ctapipe/plotting/camera.py index b0730248e85..833bacc208c 100644 --- a/ctapipe/plotting/camera.py +++ b/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): """ diff --git a/ctapipe/plotting/event_viewer.py b/ctapipe/plotting/event_viewer.py index 8c4351dc079..c8044400178 100644 --- a/ctapipe/plotting/event_viewer.py +++ b/ctapipe/plotting/event_viewer.py @@ -2,7 +2,6 @@ import matplotlib.pyplot as plt from ctapipe import visualization -from ctapipe.instrument import CameraGeometry from ctapipe.plotting.array import ArrayPlotter, NominalPlotter from numpy import ceil, sqrt from ctapipe.core import Component @@ -11,9 +10,10 @@ class EventViewer(Component): """ - Event viewer class built on top of the other plotters to allow a single view of both the camera images - and the projected Hillas parameters for a single event. Can be further modified to show the reconstructed - shower direction and core position if needed. Plus further info + Event viewer class built on top of the other plotters to allow a single + view of both the camera images and the projected Hillas parameters for a + single event. Can be further modified to show the reconstructed shower + direction and core position if needed. Plus further info """ name = 'EventViewer' test = Bool(True, help='').tag(config=True) @@ -29,7 +29,6 @@ def __init__(self, draw_hillas_planes=False): self.array_view = None self.nominal_view = None - self.geom = dict() self.cam_display = dict() self.draw_hillas_planes = draw_hillas_planes @@ -97,15 +96,12 @@ def draw_event(self, event, hillas_parameters=None): # Loop over camera images of all telescopes and create plots for ii, tel_id in zip(range(ntels), tel_list): - # Cache of camera geometries, this may go away soon - if tel_id not in self.geom: - self.geom[tel_id] = CameraGeometry.guess( - event.inst.pixel_pos[tel_id][0], - event.inst.pixel_pos[tel_id][1], - event.inst.optical_foclen[tel_id]) - ax = plt.subplot(camera_grid[ii]) - self.get_camera_view(tel_id, images.tel[tel_id].image[0], ax) + geom = event.inst.subarray.tel[tel_id].camera + self.get_camera_view(tel_id, + image=images.tel[tel_id].image[0], + ax=ax, + geom=geom) # If we want to draw the Hillas parameters in different planes we need to make a couple more viewers if self.draw_hillas_planes: @@ -129,7 +125,7 @@ def draw_event(self, event, hillas_parameters=None): return - def get_camera_view(self, tel_id, image, axis): + def get_camera_view(self, tel_id, image, axis, geom): """ Create camera viewer for a given camera image @@ -141,6 +137,8 @@ def get_camera_view(self, tel_id, image, axis): Array of calibrated pixel intensities axis: matplotlib axis Axis on which to draw plot + geom: ctapipe.instrument.CameraGeometry + Camera geometry for this tel_id Returns ------- @@ -148,7 +146,8 @@ def get_camera_view(self, tel_id, image, axis): """ #if tel_id not in self.cam_display: # Argh this is annoying, for some reason we cannot cahe the displays - self.cam_display[tel_id] = visualization.CameraDisplay(self.geom[tel_id], title="CT{0}".format(tel_id)) + self.cam_display[tel_id] = visualization.CameraDisplay(geom, title="CT{" + "0}".format(tel_id)) self.cam_display[tel_id].add_colorbar() self.cam_display[tel_id].pixels.set_antialiaseds(False) self.cam_display[tel_id].autoupdate = True diff --git a/ctapipe/plotting/tests/test_array.py b/ctapipe/plotting/tests/test_array.py index 4322daa00bb..9c880d691af 100644 --- a/ctapipe/plotting/tests/test_array.py +++ b/ctapipe/plotting/tests/test_array.py @@ -9,14 +9,13 @@ from ctapipe.coordinates import NominalFrame, CameraFrame from ctapipe.image.cleaning import tailcuts_clean from ctapipe.image.hillas import hillas_parameters, HillasParameterizationError -from ctapipe.instrument import CameraGeometry from ctapipe.io.hessio import hessio_event_source from ctapipe.plotting.array import NominalPlotter from ctapipe.utils import get_dataset + @pytest.mark.skip def test_array_draw(): - filename = get_dataset("gamma_test.simtel.gz") cam_geom = {} @@ -27,60 +26,58 @@ def test_array_draw(): calibrator = CameraDL1Calibrator(None, None) for event in source: - array_pointing=SkyCoord(event.mcheader.run_array_direction[1]*u.rad, - event.mcheader.run_array_direction[0]*u.rad, - frame=AltAz ) - #array_view = ArrayPlotter(instrument=event.inst, - # system=TiltedGroundFrame(pointing_direction=array_pointing)) + array_pointing = SkyCoord(event.mcheader.run_array_direction[1] * u.rad, + event.mcheader.run_array_direction[0] * u.rad, + frame=AltAz) + # array_view = ArrayPlotter(instrument=event.inst, + # system=TiltedGroundFrame( + # pointing_direction=array_pointing)) hillas_dict = {} r1.calibrate(event) dl0.reduce(event) - calibrator.calibrate(event) # calibrate the events + calibrator.calibrate(event) # calibrate the events # store MC pointing direction for the array for tel_id in event.dl0.tels_with_data: pmt_signal = event.dl1.tel[tel_id].image[0] - - x, y = event.inst.pixel_pos[tel_id] - fl = event.inst.optical_foclen[tel_id] - - if tel_id not in cam_geom: - cam_geom[tel_id] = CameraGeometry.guess( - event.inst.pixel_pos[tel_id][0], - event.inst.pixel_pos[tel_id][1], - event.inst.optical_foclen[tel_id]) - + geom = event.inst.subarray.tel[tel_id].camera + fl = event.inst.subarray.tel[tel_id].optics.effective_focal_length # Transform the pixels positions into nominal coordinates - camera_coord = CameraFrame(x=x, y=y, z=np.zeros(x.shape) * u.m, + camera_coord = CameraFrame(x=geom.pix_x, y=geom.pix_y, + z=np.zeros(geom.pix_x.shape) * u.m, focal_length=fl, - rotation= 90*u.deg - cam_geom[tel_id].cam_rotation ) + rotation=90 * u.deg - geom.cam_rotation) - nom_coord = camera_coord.transform_to(NominalFrame(array_direction=array_pointing, - pointing_direction=array_pointing)) + nom_coord = camera_coord.transform_to( + NominalFrame(array_direction=array_pointing, + pointing_direction=array_pointing)) - mask = tailcuts_clean(cam_geom[tel_id], pmt_signal, + mask = tailcuts_clean(geom, pmt_signal, picture_thresh=10., boundary_thresh=5.) try: moments = hillas_parameters(nom_coord.x, nom_coord.y, - pmt_signal*mask) + pmt_signal * mask) hillas_dict[tel_id] = moments - nom_coord = NominalPlotter(hillas_parameters=hillas_dict, draw_axes=True) + nom_coord = NominalPlotter(hillas_parameters=hillas_dict, + draw_axes=True) nom_coord.draw_array() except HillasParameterizationError as e: print(e) continue - #array_view.background_image(np.ones((4,4)), ((-1500,1500),(-1500,1500))) - #array_view.overlay_hillas(hillas_dict, draw_axes=True) - #array_view.draw_array(range=((-1000,1000),(-1000,1000))) - #return + # array_view.background_image(np.ones((4,4)), ((-1500,1500), + # (-1500,1500))) + # array_view.overlay_hillas(hillas_dict, draw_axes=True) + # array_view.draw_array(range=((-1000,1000),(-1000,1000))) + # return + if __name__ == "__main__": test_array_draw() diff --git a/ctapipe/plotting/tests/test_camera.py b/ctapipe/plotting/tests/test_camera.py index 95ca47e46b4..24e39ac3d96 100644 --- a/ctapipe/plotting/tests/test_camera.py +++ b/ctapipe/plotting/tests/test_camera.py @@ -8,14 +8,18 @@ def test_eventplotter(): dataset = get_dataset("gamma_test.simtel.gz") source = hessio_event_source(dataset, max_events=1) event = next(source) - data = event.r0.tel[38].adc_samples[0] + del source + + telid = list(event.r0.tels_with_data)[0] + + data = event.r0.tel[telid].adc_samples[0] plotter = CameraPlotter(event) - camera = plotter.draw_camera(38, data[:, 0]) + camera = plotter.draw_camera(telid, data[:, 0]) assert camera is not None np.testing.assert_array_equal(camera.image, data[:, 0]) - plotter.draw_camera_pixel_ids(38, [0, 1, 2]) + plotter.draw_camera_pixel_ids(telid, [0, 1, 2]) waveform = plotter.draw_waveform(data[0, :]) assert waveform is not None @@ -24,3 +28,7 @@ def test_eventplotter(): line = plotter.draw_waveform_positionline(0) assert line is not None np.testing.assert_array_equal(line.get_xdata(), [0, 0]) + + +if __name__ == '__main__': + test_eventplotter() \ No newline at end of file diff --git a/ctapipe/reco/HillasReconstructor.py b/ctapipe/reco/HillasReconstructor.py index b881eed481f..84552149ed0 100644 --- a/ctapipe/reco/HillasReconstructor.py +++ b/ctapipe/reco/HillasReconstructor.py @@ -224,7 +224,7 @@ class are set to np.nan "need at least two telescopes, have {}" .format(len(hillas_dict))) - self.get_great_circles(hillas_dict, inst, tel_phi, tel_theta) + self.get_great_circles(hillas_dict, inst.subarray, tel_phi, tel_theta) # algebraic direction estimate dir = self.fit_origin_crosses()[0] @@ -261,17 +261,18 @@ class are set to np.nan return result - def get_great_circles(self, hillas_dict, inst, tel_phi, tel_theta): + def get_great_circles(self, hillas_dict, subarray, tel_phi, tel_theta): """ - creates a dictionary of :class:`.GreatCircle` from a dictionary of hillas + creates a dictionary of :class:`.GreatCircle` from a dictionary of + hillas parameters Parameters ---------- hillas_dict : dictionary dictionary of hillas moments - inst : ctapipe instrument container - the container you get from event.inst + subarray : ctapipe.instrument.SubarrayDescription + subarray information tel_phi, tel_theta : dictionaries dictionaries of the orientation angles of the telescopes needs to contain at least the same keys as in `hillas_dict` @@ -282,6 +283,7 @@ def get_great_circles(self, hillas_dict, inst, tel_phi, tel_theta): p2_x = moments.cen_x + moments.length * np.cos(moments.psi) p2_y = moments.cen_y + moments.length * np.sin(moments.psi) + foclen = subarray.tel[tel_id].optics.effective_focal_length circle = GreatCircle( guess_pix_direction(np.array([moments.cen_x / u.m, @@ -289,10 +291,10 @@ def get_great_circles(self, hillas_dict, inst, tel_phi, tel_theta): np.array([moments.cen_y / u.m, p2_y / u.m]) * u.m, tel_phi[tel_id], tel_theta[tel_id], - inst.optical_foclen[tel_id]), + foclen), moments.size * (moments.length / moments.width) ) - circle.pos = inst.tel_pos[tel_id] + circle.pos = subarray.positions[tel_id] self.circles[tel_id] = circle def fit_origin_crosses(self): diff --git a/ctapipe/tools/dump_instrument.py b/ctapipe/tools/dump_instrument.py index d2ee9069a2e..98d70cf9273 100644 --- a/ctapipe/tools/dump_instrument.py +++ b/ctapipe/tools/dump_instrument.py @@ -5,13 +5,35 @@ automatically generated. """ -from ctapipe.core.traits import (Unicode, Dict, Bool, Enum) +from collections import defaultdict + +import numpy as np +from astropy import units as u +from astropy.table import Table + from ctapipe.core import Tool +from ctapipe.core.traits import (Unicode, Dict, Enum) from ctapipe.io.hessio import hessio_event_source -from ctapipe.instrument import CameraGeometry, get_camera_types, print_camera_types -from astropy.table import Table -from astropy import units as u -import numpy as np + + +def get_camera_types(subarray): + """ return dict of camera names mapped to a list of tel_ids + that use that camera + + Parameters + ---------- + subarray: ctapipe.instrument.SubarrayDescription + + """ + + cam_types = defaultdict(list) + + for telid in subarray.tel: + geom = subarray.tel[telid].camera + cam_types[geom.cam_id].append(telid) + + return cam_types + class DumpInstrumentTool(Tool): description = Unicode(__doc__) @@ -55,24 +77,22 @@ def _get_file_format_info(self, format_name, table_type, cam_name): raise NameError("format not supported") def write_camera_geometries(self): - cam_types = get_camera_types(self.inst) - print_camera_types(self.inst, printer=self.log.info) + cam_types = get_camera_types(self.inst.subarray) + self.inst.subarray.info(printer=self.log.info) for cam_name in cam_types: ext, args = self._get_file_format_info(self.format, 'CAMGEOM', cam_name) self.log.debug("writing {}".format(cam_name)) tel_id = cam_types[cam_name].pop() - pix = self.inst.pixel_pos[tel_id] - flen = self.inst.optical_foclen[tel_id] - geom = CameraGeometry.guess(*pix, flen) + geom = self.inst.subarray.tel[tel_id].camera table = geom.to_table() table.meta['SOURCE'] = self.infile table.write("{}.camgeom.{}".format(cam_name, ext), **args) def write_optics_descriptions(self): # we'll use the camera names for the optics names... - cam_types = get_camera_types(self.inst) + cam_types = get_camera_types(self.inst.subarray) cam_names = list(cam_types.keys()) optical_foclens = [] mirror_dish_areas = [] diff --git a/examples/display_summed_images.py b/examples/display_summed_images.py index 4c2b200b16f..1022f1de0b7 100644 --- a/examples/display_summed_images.py +++ b/examples/display_summed_images.py @@ -8,19 +8,12 @@ from astropy.table import Table from ctapipe.core import Tool from ctapipe.core.traits import * -from ctapipe.instrument import CameraGeometry, SubarrayDescription from ctapipe.io.hessio import hessio_event_source from ctapipe.visualization import CameraDisplay from matplotlib import pyplot as plt from ctapipe.calib import CameraCalibrator -def get_camera_types(subarray: SubarrayDescription): - """ return dict of group_id to list of telescopes in group, - where a group is defined as similar telescopes""" - return - - class ImageSumDisplayerTool(Tool): description = Unicode(__doc__) name = "ctapipe-image-sum-display"