Python CameraGeometry.guess方法代码示例

# 需要导入模块: from ctapipe.instrument import CameraGeometry [as 别名]
# 或者: from ctapipe.instrument.CameraGeometry import guess [as 别名]
def display_event(event):
    """an extremely inefficient display. It creates new instances of
    CameraDisplay for every event and every camera, and also new axes
    for each event. It's hacked, but it works
    """
    print("Displaying... please wait (this is an inefficient implementation)")
    global fig
    ntels = len(event.r0.tels_with_data)
    fig.clear()

    plt.suptitle("EVENT {}".format(event.r0.event_id))

    disps = []

    for ii, tel_id in enumerate(event.r0.tels_with_data):
        print("\t draw cam {}...".format(tel_id))
        nn = int(ceil(sqrt(ntels)))
        ax = plt.subplot(nn, nn, ii + 1)

        x, y = event.inst.pixel_pos[tel_id]
        geom = CameraGeometry.guess(x, y, event.inst.optical_foclen[tel_id])
        disp = CameraDisplay(geom, ax=ax, title="CT{0}".format(tel_id))
        disp.pixels.set_antialiaseds(False)
        disp.autoupdate = False
        disp.cmap = random.choice(cmaps)
        chan = 0
        signals = event.r0.tel[tel_id].adc_sums[chan].astype(float)
        signals -= signals.mean()
        disp.image = signals
        disp.set_limits_percent(95)
        disp.add_colorbar()
        disps.append(disp)

    return disps

# 需要导入模块: from ctapipe.instrument import CameraGeometry [as 别名]
# 或者: from ctapipe.instrument.CameraGeometry import guess [as 别名]
def test_FitGammaHillas():
    '''
    a test of the complete fit procedure on one event including:
    • tailcut cleaning
    • hillas parametrisation
    • GreatCircle creation
    • direction fit
    • position fit

    in the end, proper units in the output are asserted '''

    filename = get_dataset("gamma_test.simtel.gz")

    fit = HillasReconstructor()

    cam_geom = {}
    tel_phi = {}
    tel_theta = {}

    source = hessio_event_source(filename)

    for event in source:

        hillas_dict = {}
        for tel_id in event.dl0.tels_with_data:

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

                tel_phi[tel_id] = event.mc.tel[tel_id].azimuth_raw * u.rad
                tel_theta[tel_id] = (np.pi/2-event.mc.tel[tel_id].altitude_raw)*u.rad

            pmt_signal = event.r0.tel[tel_id].adc_sums[0]

            mask = tailcuts_clean(cam_geom[tel_id], pmt_signal,
                                  picture_thresh=10., boundary_thresh=5.)
            pmt_signal[mask == 0] = 0

            try:
                moments = hillas_parameters(event.inst.pixel_pos[tel_id][0],
                                            event.inst.pixel_pos[tel_id][1],
                                            pmt_signal)
                hillas_dict[tel_id] = moments
            except HillasParameterizationError as e:
                print(e)
                continue

        if len(hillas_dict) < 2: continue

        fit_result = fit.predict(hillas_dict, event.inst, tel_phi, tel_theta)

        print(fit_result)
        fit_result.alt.to(u.deg)
        fit_result.az.to(u.deg)
        fit_result.core_x.to(u.m)
        assert fit_result.is_valid
        return

# 需要导入模块: from ctapipe.instrument import CameraGeometry [as 别名]
# 或者: from ctapipe.instrument.CameraGeometry import guess [as 别名]
def display_telescope(calibrated_data,cam,tel_id,pdffilename):
    logger.info("Tel_ID %d\n"%tel_id)
    pp = PdfPages("%s"%pdffilename)
    global fig
    for event in calibrated_data:
        tels = list(event.r0.tels_with_data)
        logger.debug(tels)
        # Select telescope
        if tel_id not in tels:
            continue

        fig.clear()

        # geom = cam['CameraTable_VersionFeb2016_TelID%s'%tel_id]
        geom = CameraGeometry.guess(*event.inst.pixel_pos[tel_id],
                                       event.inst.optical_foclen[tel_id])
        # Select number of pads to display. It depends on the numberr of gains:
        nchan = event.inst.num_channels[tel_id]

        plt.suptitle("EVENT {} {:.1e} TeV @({:.1f},{:.1f})deg @{:.1f} m".format(
            event.r0.event_id, event.mc.energy,
            event.mc.alt, event.mc.az,
            np.sqrt(pow(event.mc.core_x, 2) +
                    pow(event.mc.core_y, 2))))
        print("\t draw cam {} (gains={})...".format(tel_id,nchan))
        ax=[]
        disp=[]
        signals=[]
        npads=0
        for i in range(nchan):
            npads += 1
            # Display the camera charge (HG/LG)
            ax.append(plt.subplot(nchan, 2, npads))
            disp.append(visualization.CameraDisplay(geom, ax=ax[-1],
                                                    title="CT{} [{} ADC cts]".format(tel_id,gain_label[i])))
            disp[-1].pixels.set_antialiaseds(False)
            signals.append(event.r0.tel[tel_id].adc_sums[i])
            disp[-1].image = signals[-1]
            disp[-1].pixels.set_cmap(get_cmap(disp[-1].image))
            disp[-1].add_colorbar(label=" [{} ADC cts]".format(gain_label[i]))

            # Display the camera charge for significant pixels (HG/LG)
            npads += 1
            ax.append(plt.subplot(nchan, 2, npads))
            disp.append(visualization.CameraDisplay(geom, ax=ax[-1],
                                                    title="CT{} [{} ADC cts]".format(tel_id,gain_label[i])))
            disp[-1].pixels.set_antialiaseds(False)
            signals.append(event.r0.tel[tel_id].adc_sums[i])
            m = (get_zero_sup_mode(tel_id) & 0x001) or (get_significant(tel_id) & 0x020)
            disp[-1].image = signals[-1]*(m/m.max())
            disp[-1].pixels.set_cmap(get_cmap(disp[-1].image))
            disp[-1].add_colorbar(label=" [{} ADC cts]".format(gain_label[i]))

            
        plt.pause(0.1)
        pp.savefig(fig)

    pp.close()

# 需要导入模块: from ctapipe.instrument import CameraGeometry [as 别名]
# 或者: from ctapipe.instrument.CameraGeometry import guess [as 别名]
    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]

# 需要导入模块: from ctapipe.instrument import CameraGeometry [as 别名]
# 或者: from ctapipe.instrument.CameraGeometry import guess [as 别名]
 def write_camera_geometries(self):
     cam_types = get_camera_types(self.inst)
     print_camera_types(self.inst, 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)
         table = geom.to_table()
         table.meta['SOURCE'] = self.infile
         table.write("{}.camgeom.{}".format(cam_name, ext), **args)

# 需要导入模块: from ctapipe.instrument import CameraGeometry [as 别名]
# 或者: from ctapipe.instrument.CameraGeometry import guess [as 别名]
    def get_geometry(event, telid):
        """
        Obtain the geometry for this telescope.

        Parameters
        ----------
        event : container
            A `ctapipe` event container
        telid : int
            The telescope id.
            The neighbours are calculated once per telescope.

        Returns
        -------
        `CameraGeometry`
        """
        return CameraGeometry.guess(*event.inst.pixel_pos[telid],
                                    event.inst.optical_foclen[telid])

# 需要导入模块: from ctapipe.instrument import CameraGeometry [as 别名]
# 或者: from ctapipe.instrument.CameraGeometry import guess [as 别名]
def test_nb_peak_integration():
    telid = 11
    event = get_test_event()
    data = event.r0.tel[telid].adc_samples
    nsamples = data.shape[2]
    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])
    nei = geom.neighbors

    integrator = NeighbourPeakIntegrator(None, None)
    integrator.neighbours = nei
    integration, peakpos, window = integrator.extract_charge(data_ped)

    assert_almost_equal(integration[0][0], -64, 0)
    assert_almost_equal(integration[1][0], -64, 0)
    assert peakpos[0][0] == 20
    assert peakpos[1][0] == 20

# 需要导入模块: from ctapipe.instrument import CameraGeometry [as 别名]
# 或者: from ctapipe.instrument.CameraGeometry import guess [as 别名]
    def start(self):
        geom = None
        imsum = None
        disp = None

        for data in hessio_event_source(self.infile,
                                        allowed_tels=self._selected_tels,
                                        max_events=self.max_events):

            self.calibrator.calibrate(data)

            if geom is None:
                x, y = data.inst.pixel_pos[self._base_tel]
                flen = data.inst.optical_foclen[self._base_tel]
                geom = CameraGeometry.guess(x, y, flen)
                imsum = np.zeros(shape=x.shape, dtype=np.float)
                disp = CameraDisplay(geom, title=geom.cam_id)
                disp.add_colorbar()
                disp.cmap = 'viridis'

            if len(data.dl0.tels_with_data) <= 2:
                continue

            imsum[:] = 0
            for telid in data.dl0.tels_with_data:
                imsum += data.dl1.tel[telid].image[0]

            self.log.info("event={} ntels={} energy={}" \
                          .format(data.r0.event_id,
                                  len(data.dl0.tels_with_data),
                                  data.mc.energy))
            disp.image = imsum
            plt.pause(0.1)

            if self.output_suffix is not "":
                filename = "{:020d}{}".format(data.r0.event_id,
                                              self.output_suffix)
                self.log.info("saving: '{}'".format(filename))
                plt.savefig(filename)

# 需要导入模块: from ctapipe.instrument import CameraGeometry [as 别名]
# 或者: from ctapipe.instrument.CameraGeometry import guess [as 别名]
def plot_muon_event(event, muonparams, geom_dict=None, 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

        for tel_id in event.dl0.tels_with_data:
            npads = 2
            # 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
            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

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

            muon_phi = np.arctan(muonparams[0].ring_center_y /
                                 muonparams[0].ring_center_x)

            rotr_angle = geom.pix_rotation
            # if event.inst.optical_foclen[tel_id] > 10.*u.m and
            # event.dl0.tel[tel_id].num_pixels != 1764:
            if geom.cam_id == 'LSTCam' or geom.cam_id == 'NectarCam':
                #print("Resetting the rotation angle")
                rotr_angle = 0. * u.deg

            # Convert to camera frame (centre & radius)
            altaz = HorizonFrame(alt=event.mc.alt, az=event.mc.az)

            ring_nominal = NominalFrame(x=muonparams[0].ring_center_x,
                                        y=muonparams[0].ring_center_y,
                                        array_direction=altaz,
                                        pointing_direction=altaz)

            # embed()
            ring_camcoord = ring_nominal.transform_to(CameraFrame(
                pointing_direction=altaz,
                focal_length=event.inst.optical_foclen[tel_id],
                rotation=rotr_angle))

            centroid_rad = np.sqrt(ring_camcoord.y**2 + ring_camcoord.x**2)
            centroid = (ring_camcoord.x.value, ring_camcoord.y.value)

            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,
                                       z=np.zeros(px.shape) * u.m,
                                       focal_length=flen,
                                       rotation=geom.pix_rotation)

            nom_coord = camera_coord.transform_to(
                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)

#.........这里部分代码省略.........

# 需要导入模块: from ctapipe.instrument import CameraGeometry [as 别名]
# 或者: from ctapipe.instrument.CameraGeometry import guess [as 别名]
def analyze_muon_event(event, params=None, geom_dict=None):
    """
    Generic muon event analyzer. 

    Parameters
    ----------
    event : ctapipe dl1 event container


    Returns
    -------
    muonringparam, muonintensityparam : MuonRingParameter and MuonIntensityParameter container event

    """
    # Declare a dict to define the muon cuts (ASTRI and SCT missing)
    muon_cuts = {}

    names = ['LST:LSTCam','MST:NectarCam','MST:FlashCam','MST-SCT:SCTCam','SST-1M:DigiCam','SST-GCT:CHEC','SST-ASTRI:ASTRICam']
    TailCuts = [(5,7),(5,7),(10,12),(5,7),(5,7),(5,7),(5,7)] #10,12?
    impact = [(0.2,0.9),(0.1,0.95),(0.2,0.9),(0.2,0.9),(0.1,0.95),(0.1,0.95),(0.1,0.95)]
    ringwidth = [(0.04,0.08),(0.02,0.1),(0.01,0.1),(0.02,0.1),(0.01,0.5),(0.02,0.2),(0.02,0.2)]
    TotalPix = [1855.,1855.,1764.,11328.,1296.,2048.,2368.]#8% (or 6%) as limit
    MinPix = [148.,148.,141.,680.,104.,164.,142.]
    #Need to either convert from the pixel area in m^2 or check the camera specs
    AngPixelWidth = [0.1,0.2,0.18,0.067,0.24,0.2,0.17] #Found from TDRs (or the pixel area)
    hole_rad = []#Need to check and implement
    cam_rad = [2.26,3.96,3.87,4.,4.45,2.86,5.25]#Found from the field of view calculation
    sec_rad = [0.*u.m,0.*u.m,0.*u.m,2.7*u.m,0.*u.m,1.*u.m,1.8*u.m]
    sct = [False,False,False,True,False,True,True]


    muon_cuts = {'Name':names,'TailCuts':TailCuts,'Impact':impact,'RingWidth':ringwidth,'TotalPix':TotalPix,'MinPix':MinPix,'CamRad':cam_rad,'SecRad':sec_rad,'SCT':sct,'AngPixW':AngPixelWidth}
    #print(muon_cuts)

    muonringlist = []#[None] * len(event.dl0.tels_with_data)
    muonintensitylist = []#[None] * len(event.dl0.tels_with_data)
    tellist = []
    #for tid in event.dl0.tels_with_data:
    #    tellist.append(tid)
    muon_event_param = {'TelIds':tellist,'MuonRingParams':muonringlist,'MuonIntensityParams':muonintensitylist}
    #muonringparam = None
    #muonintensityparam = None

    for telid in event.dl0.tels_with_data:

        #print("Analysing muon event for tel",telid)
        muonringparam = None
        muonintensityparam = None
        #idx = muon_event_param['TelIds'].index(telid)

        x, y = event.inst.pixel_pos[telid]

        #image = event.dl1.tel[telid].calibrated_image
        image = event.dl1.tel[telid].image[0]

        # Get geometry
        geom = None
        if geom_dict is not None and telid in geom_dict:
            geom = geom_dict[telid]
        else:
            log.debug("[calib] Guessing camera geometry")
            geom = CameraGeometry.guess(*event.inst.pixel_pos[telid],
                                        event.inst.optical_foclen[telid])
            log.debug("[calib] Camera geometry found")
            if geom_dict is not None:
                geom_dict[telid] = geom


        teldes = event.inst.subarray.tel[telid]
        dict_index = muon_cuts['Name'].index(str(teldes))
        #print('found an index of',dict_index,'for camera',geom.cam_id)

        #tailcuts = (5.,7.)
        tailcuts = muon_cuts['TailCuts'][dict_index]

        #print("Tailcuts are",tailcuts[0],tailcuts[1])

        #rot_angle = 0.*u.deg
        #if event.inst.optical_foclen[telid] > 10.*u.m and event.dl0.tel[telid].num_pixels != 1764:
            #rot_angle = -100.14*u.deg

        clean_mask = tailcuts_clean(geom,image,picture_thresh=tailcuts[0],boundary_thresh=tailcuts[1])
        camera_coord = CameraFrame(x=x,y=y,z=np.zeros(x.shape)*u.m, focal_length = event.inst.optical_foclen[telid], rotation=geom.pix_rotation)

        #print("Camera",geom.cam_id,"focal length",event.inst.optical_foclen[telid],"rotation",geom.pix_rotation)
        #TODO: correct this hack for values over 90
        altval = event.mcheader.run_array_direction[1]
        if (altval > np.pi/2.):
            altval = np.pi/2.

        altaz = HorizonFrame(alt=altval*u.rad,az=event.mcheader.run_array_direction[0]*u.rad)
        nom_coord = camera_coord.transform_to(NominalFrame(array_direction=altaz,pointing_direction=altaz))

        
        x = nom_coord.x.to(u.deg)
        y = nom_coord.y.to(u.deg)

        img = image*clean_mask
        noise = 5.
        weight = img / (img+noise)
#.........这里部分代码省略.........

# 需要导入模块: from ctapipe.instrument import CameraGeometry [as 别名]
# 或者: from ctapipe.instrument.CameraGeometry import guess [as 别名]
    def draw_event(self, event, hillas_parameters=None):
        """
        Draw display for a given event

        Parameters
        ----------
        event: ctapipe event object
        hillas_parameters: dict
            Dictionary of Hillas parameters (in nominal system)

        Returns
        -------
            None
        """
        tel_list = event.r0.tels_with_data
        images = event.dl1

        # First close any plots that already exist
        plt.close()
        ntels = len(tel_list)

        fig = plt.figure(figsize=(20, 20 * 0.66))

        # If we want to draw the Hillas parameters in different planes we need to split our figure
        if self.draw_hillas_planes:
            y_axis_split = 2
        else:
            y_axis_split = 1

        outer_grid = gridspec.GridSpec(1, y_axis_split, width_ratios=[y_axis_split, 1])
        # Create a square grid for camera images
        nn = int(ceil(sqrt(ntels)))
        nx = nn
        ny = nn

        while nx * ny >= ntels:
            ny-=1
        ny+=1
        while nx * ny >= ntels:
            nx -= 1
        nx += 1

        camera_grid = gridspec.GridSpecFromSubplotSpec(ny, nx, subplot_spec=outer_grid[0])

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

        # If we want to draw the Hillas parameters in different planes we need to make a couple more viewers
        if self.draw_hillas_planes:
            # Split the second sub figure into two further figures
            reco_grid = gridspec.GridSpecFromSubplotSpec(2, 1, subplot_spec=outer_grid[1])
            # Create plot of telescope positions at ground level
            array = ArrayPlotter(telescopes=tel_list, instrument=event.inst,# system=tilted_system,
                                ax=plt.subplot(reco_grid[0]))
            # Draw MC position (this should change later)
            array.draw_position(event.mc.core_x, event.mc.core_y, use_centre=True)
            array.draw_array(((-300,300),(-300,300)))

            # If we have valid Hillas parameters we should draw them in the Nominal system
            if hillas_parameters is not None:
                array.overlay_hillas(hillas_parameters, draw_axes=True)

                nominal =  NominalPlotter(hillas_parameters=hillas_parameters, draw_axes=True, ax=plt.subplot(reco_grid[1]))
                nominal.draw_array()

        plt.show()

        return

# 需要导入模块: from ctapipe.instrument import CameraGeometry [as 别名]
# 或者: from ctapipe.instrument.CameraGeometry import guess [as 别名]
def extract_images(simtel_file_path,
                   tel_id_filter_list=None,
                   event_id_filter_list=None,
                   output_directory=None):

    # EXTRACT IMAGES ##########################################################

    # hessio_event_source returns a Python generator that streams data from an
    # EventIO/HESSIO MC data file (e.g. a standard CTA data file).
    # This generator contains ctapipe.core.Container instances ("event").
    # 
    # Parameters:
    # - max_events: maximum number of events to read
    # - allowed_tels: select only a subset of telescope, if None, all are read.
    source = hessio_event_source(simtel_file_path, allowed_tels=tel_id_filter_list)

    # ITERATE OVER EVENTS #####################################################

    calib = CameraCalibrator(None, None)

    for event in source:

        calib.calibrate(event)  # calibrate the event

        event_id = int(event.dl0.event_id)

        if (event_id_filter_list is None) or (event_id in event_id_filter_list):

            #print("event", event_id)

            # ITERATE OVER IMAGES #############################################

            for tel_id in event.trig.tels_with_trigger:

                tel_id = int(tel_id)

                if tel_id in tel_id_filter_list:

                    #print("telescope", tel_id)

                    # CHECK THE IMAGE GEOMETRY ################################

                    #print("checking geometry")

                    x, y = event.inst.pixel_pos[tel_id]
                    foclen = event.inst.optical_foclen[tel_id]
                    geom = CameraGeometry.guess(x, y, foclen)

                    if (geom.pix_type != "hexagonal") or (geom.cam_id != "LSTCam"):
                        raise ValueError("Telescope {}: error (the input image is not a valide LSTCam telescope image) -> {} ({})".format(tel_id, geom.pix_type, geom.cam_id))

                    # GET IMAGES ##############################################

                    pe_image = event.mc.tel[tel_id].photo_electron_image   # 1D np array

                    #uncalibrated_image = event.dl0.tel[tel_id].adc_sums  # ctapipe 0.3.0
                    uncalibrated_image = event.r0.tel[tel_id].adc_sums    # ctapipe 0.4.0
                    pedestal = event.mc.tel[tel_id].pedestal
                    gain = event.mc.tel[tel_id].dc_to_pe
                    pixel_pos = event.inst.pixel_pos[tel_id]

                    calibrated_image = event.dl1.tel[tel_id].image

                    calibrated_image[1, calibrated_image[0,:] <= LST_CAM_CHANNEL_THRESHOLD] = 0
                    calibrated_image[0, calibrated_image[0,:] >  LST_CAM_CHANNEL_THRESHOLD] = 0
                    calibrated_image = calibrated_image.sum(axis=0)

                    #print(pe_image.shape)
                    #print(calibrated_image.shape)
                    #print(uncalibrated_image.shape)
                    #print(pedestal.shape)
                    #print(gain.shape)
                    #print(pixel_pos.shape)
                    #print(pixel_pos[0])
                    #print(pixel_pos[1])

                    # CONVERTING GEOMETRY (1D TO 2D) ##########################

                    buffer_id_str = geom.cam_id + "0"

                    geom2d, pe_image_2d =           ctapipe_geom_converter.convert_geometry_1d_to_2d(geom, pe_image,           buffer_id_str, add_rot=0)
                    geom2d, calibrated_image_2d =   ctapipe_geom_converter.convert_geometry_1d_to_2d(geom, calibrated_image,   buffer_id_str, add_rot=0)

                    geom2d, uncalibrated_image_2d_ch0 = ctapipe_geom_converter.convert_geometry_1d_to_2d(geom, uncalibrated_image[0], buffer_id_str, add_rot=0)
                    geom2d, uncalibrated_image_2d_ch1 = ctapipe_geom_converter.convert_geometry_1d_to_2d(geom, uncalibrated_image[1], buffer_id_str, add_rot=0)
                    geom2d, pedestal_2d_ch0 =           ctapipe_geom_converter.convert_geometry_1d_to_2d(geom, pedestal[0],           buffer_id_str, add_rot=0)
                    geom2d, pedestal_2d_ch1 =           ctapipe_geom_converter.convert_geometry_1d_to_2d(geom, pedestal[1],           buffer_id_str, add_rot=0)
                    geom2d, gains_2d_ch0 =              ctapipe_geom_converter.convert_geometry_1d_to_2d(geom, gain[0],               buffer_id_str, add_rot=0)
                    geom2d, gains_2d_ch1 =              ctapipe_geom_converter.convert_geometry_1d_to_2d(geom, gain[1],               buffer_id_str, add_rot=0)

                    # Make a mock pixel position array...
                    pixel_pos_2d = np.array(np.meshgrid(np.linspace(pixel_pos[0].min(), pixel_pos[0].max(), pe_image_2d.shape[0]),
                                                        np.linspace(pixel_pos[1].min(), pixel_pos[1].max(), pe_image_2d.shape[1])))

                    ###########################################################

                    # The ctapipe geometry converter operate on one channel
                    # only and then takes and return a 2D array but datapipe
                    # fits files keep all channels and thus takes 3D arrays...

#.........这里部分代码省略.........

# 需要导入模块: from ctapipe.instrument import CameraGeometry [as 别名]
# 或者: from ctapipe.instrument.CameraGeometry import guess [as 别名]
def test_guess_camera():
    px = np.linspace(-10, 10, 11328) * u.m
    py = np.linspace(-10, 10, 11328) * u.m
    geom = CameraGeometry.guess(px, py,0 * u.m)
    assert geom.pix_type.startswith('rect')

# 需要导入模块: from ctapipe.instrument import CameraGeometry [as 别名]
# 或者: from ctapipe.instrument.CameraGeometry import guess [as 别名]
    def plot(self, input_file, event, telid, chan, extractor_name):
        # Extract required images
        dl0 = event.dl0.tel[telid].pe_samples[chan]

        t_pe = event.mc.tel[telid].photo_electron_image
        dl1 = event.dl1.tel[telid].image[chan]
        max_time = np.unravel_index(np.argmax(dl0), dl0.shape)[1]
        max_charges = np.max(dl0, axis=1)
        max_pix = int(np.argmax(max_charges))
        min_pix = int(np.argmin(max_charges))

        geom = CameraGeometry.guess(*event.inst.pixel_pos[telid],
                                    event.inst.optical_foclen[telid])

        nei = geom.neighbors

        # Get Neighbours
        max_pixel_nei = nei[max_pix]
        min_pixel_nei = nei[min_pix]

        # Get Windows
        windows = event.dl1.tel[telid].extracted_samples[chan]
        length = np.sum(windows, axis=1)
        start = np.argmax(windows, axis=1)
        end = start + length - 1

        # Draw figures
        ax_max_nei = {}
        ax_min_nei = {}
        fig_waveforms = plt.figure(figsize=(18, 9))
        fig_waveforms.subplots_adjust(hspace=.5)
        fig_camera = plt.figure(figsize=(15, 12))

        ax_max_pix = fig_waveforms.add_subplot(4, 2, 1)
        ax_min_pix = fig_waveforms.add_subplot(4, 2, 2)
        ax_max_nei[0] = fig_waveforms.add_subplot(4, 2, 3)
        ax_min_nei[0] = fig_waveforms.add_subplot(4, 2, 4)
        ax_max_nei[1] = fig_waveforms.add_subplot(4, 2, 5)
        ax_min_nei[1] = fig_waveforms.add_subplot(4, 2, 6)
        ax_max_nei[2] = fig_waveforms.add_subplot(4, 2, 7)
        ax_min_nei[2] = fig_waveforms.add_subplot(4, 2, 8)

        ax_img_nei = fig_camera.add_subplot(2, 2, 1)
        ax_img_max = fig_camera.add_subplot(2, 2, 2)
        ax_img_true = fig_camera.add_subplot(2, 2, 3)
        ax_img_cal = fig_camera.add_subplot(2, 2, 4)

        # Draw max pixel traces
        ax_max_pix.plot(dl0[max_pix])
        ax_max_pix.set_xlabel("Time (ns)")
        ax_max_pix.set_ylabel("DL0 Samples (ADC)")
        ax_max_pix.set_title("(Max) Pixel: {}, True: {}, Measured = {:.3f}"
                             .format(max_pix, t_pe[max_pix], dl1[max_pix]))
        max_ylim = ax_max_pix.get_ylim()
        ax_max_pix.plot([start[max_pix], start[max_pix]],
                        ax_max_pix.get_ylim(), color='r', alpha=1)
        ax_max_pix.plot([end[max_pix], end[max_pix]],
                        ax_max_pix.get_ylim(), color='r', alpha=1)
        for i, ax in ax_max_nei.items():
            if len(max_pixel_nei) > i:
                pix = max_pixel_nei[i]
                ax.plot(dl0[pix])
                ax.set_xlabel("Time (ns)")
                ax.set_ylabel("DL0 Samples (ADC)")
                ax.set_title("(Max Nei) Pixel: {}, True: {}, Measured = {:.3f}"
                             .format(pix, t_pe[pix], dl1[pix]))
                ax.set_ylim(max_ylim)
                ax.plot([start[pix], start[pix]],
                        ax.get_ylim(), color='r', alpha=1)
                ax.plot([end[pix], end[pix]],
                        ax.get_ylim(), color='r', alpha=1)

        # Draw min pixel traces
        ax_min_pix.plot(dl0[min_pix])
        ax_min_pix.set_xlabel("Time (ns)")
        ax_min_pix.set_ylabel("DL0 Samples (ADC)")
        ax_min_pix.set_title("(Min) Pixel: {}, True: {}, Measured = {:.3f}"
                             .format(min_pix, t_pe[min_pix], dl1[min_pix]))
        ax_min_pix.set_ylim(max_ylim)
        ax_min_pix.plot([start[min_pix], start[min_pix]],
                        ax_min_pix.get_ylim(), color='r', alpha=1)
        ax_min_pix.plot([end[min_pix], end[min_pix]],
                        ax_min_pix.get_ylim(), color='r', alpha=1)
        for i, ax in ax_min_nei.items():
            if len(min_pixel_nei) > i:
                pix = min_pixel_nei[i]
                ax.plot(dl0[pix])
                ax.set_xlabel("Time (ns)")
                ax.set_ylabel("DL0 Samples (ADC)")
                ax.set_title("(Min Nei) Pixel: {}, True: {}, Measured = {:.3f}"
                             .format(pix, t_pe[pix], dl1[pix]))
                ax.set_ylim(max_ylim)
                ax.plot([start[pix], start[pix]],
                        ax.get_ylim(), color='r', alpha=1)
                ax.plot([end[pix], end[pix]],
                        ax.get_ylim(), color='r', alpha=1)

        # Draw cameras
        nei_camera = np.zeros_like(max_charges, dtype=np.int)
        nei_camera[min_pixel_nei] = 2
#.........这里部分代码省略.........

# 需要导入模块: from ctapipe.instrument import CameraGeometry [as 别名]
# 或者: from ctapipe.instrument.CameraGeometry import guess [as 别名]
 def get_geometry(self, event, telid):
     if telid not in self._geom_dict:
         geom = CameraGeometry.guess(*event.inst.pixel_pos[telid],
                                     event.inst.optical_foclen[telid])
         self._geom_dict[telid] = geom
     return self._geom_dict[telid]