Python calib.CameraCalibrator类代码示例

    def setup(self):
        # load up the telescope types table (need to first open a file, a bit of
        # a hack until a proper insturment module exists) and select only the
        # telescopes with the same camera type

        self.reader = SimTelEventSource(
            input_url=self.infile, max_events=self.max_events
        )

        for event in self.reader:
            camtypes = event.inst.subarray.to_table().group_by('camera_type')
            event.inst.subarray.info(printer=self.log.info)
            break

        group = camtypes.groups[self.telgroup]
        self._selected_tels = list(group['id'].data)
        self._base_tel = self._selected_tels[0]
        self.log.info(
            "Telescope group %d: %s", self.telgroup,
            str(event.inst.subarray.tel[self._selected_tels[0]])
        )
        self.log.info(f"SELECTED TELESCOPES:{self._selected_tels}")

        self.calibrator = CameraCalibrator(
            parent=self, eventsource=self.reader
        )

        self.reader.allowed_tels = self._selected_tels

    def setup(self):
        self.log.info('Configure EventSourceFactory...')

        self.event_source = EventSourceFactory.produce(
            config=self.config, tool=self, product='SimTelEventSource'
        )
        self.event_source.allowed_tels = self.config['Analysis']['allowed_tels']

        self.calibrator = CameraCalibrator(
            config=self.config, tool=self, eventsource=self.event_source
        )

        self.writer = HDF5TableWriter(
            filename=self.outfile, group_name='image_infos', overwrite=True
        )

        # Define Pre-selection for images
        preselcuts = self.config['Preselect']
        self.image_cutflow = CutFlow('Image preselection')
        self.image_cutflow.set_cuts(dict(
            no_sel=None,
            n_pixel=lambda s: np.count_nonzero(s) < preselcuts['n_pixel']['min'],
            image_amplitude=lambda q: q < preselcuts['image_amplitude']['min']
        ))

        # Define Pre-selection for events
        self.event_cutflow = CutFlow('Event preselection')
        self.event_cutflow.set_cuts(dict(
            no_sel=None
        ))

def test_basic_muon_reco(example_event):
    """
    Really simplistic test: just run the analyze_muon_event code, to make
    sure it doesn't crash. The input event is so far not a muon, so no output
    is generated.

    Parameters
    ----------
    test_event - a sample event (fixture)

    """

    calib = CameraCalibrator()
    calib.calibrate(example_event)

    muon_params = muon.analyze_muon_event(example_event)
    assert muon_params is not None

    def setup(self):
        print('TOLLES INFILE', self.infile)
        self.event_source = EventSource.from_url(self.infile, parent=self)
        self.event_source.allowed_tels = {self.tel, }

        self.calibrator = CameraCalibrator(parent=self)

        self.log.info(f'SELECTING EVENTS FROM TELESCOPE {self.tel}')

 def setup(self):
     if self.events == '':
         raise ToolConfigurationError("please specify --input <events file>")
     self.log.debug("input: %s", self.events)
     self.source = EventSourceFactory.produce(input_url=self.events)
     self.calib = CameraCalibrator(
         config=self.config, tool=self, eventsource=self.source
     )
     self.writer = HDF5TableWriter(self.outfile, "muons")

    def setup(self):
        self.eventsource = EventSource.from_url(
            get_dataset_path("gamma_test_large.simtel.gz"),
            parent=self,
        )

        self.calibrator = CameraCalibrator(parent=self)

        self.plotter = ImagePlotter(parent=self)

    def setup(self):
        kwargs = dict(config=self.config, tool=self)

        reader_factory = EventFileReaderFactory(**kwargs)
        reader_class = reader_factory.get_class()
        self.reader = reader_class(**kwargs)

        self.calibrator = CameraCalibrator(origin=self.reader.origin, **kwargs)

        self.plotter = ImagePlotter(**kwargs)

    def setup(self):

        reader_factory = EventFileReaderFactory(None, self)
        reader_class = reader_factory.get_class()
        self.reader = reader_class(None, self)

        self.calibrator = CameraCalibrator(config=None, tool=self,
                                           origin=self.reader.origin)
        self.source = self.reader.read(allowed_tels=[self.tel, ])

        self.log.info('SELECTING EVENTS FROM TELESCOPE {}'.format(self.tel))

    def setup(self):
        kwargs = dict(config=self.config, tool=self)

        self.eventsource = EventSourceFactory.produce(
            input_url=get_dataset_path("gamma_test.simtel.gz"), **kwargs
        )

        self.calibrator = CameraCalibrator(
            eventsource=self.eventsource, **kwargs
        )

        self.plotter = ImagePlotter(**kwargs)

    def setup(self):

        self.event_source = EventSourceFactory.produce(
            config=self.config, tool=self
        )
        self.event_source.allowed_tels = [
            self.tel,
        ]

        self.calibrator = CameraCalibrator(
            config=self.config, tool=self, eventsource=self.event_source
        )

        self.log.info(f'SELECTING EVENTS FROM TELESCOPE {self.tel}')

 def setup(self):
     # load up the telescope types table (need to first open a file, a bit of
     # a hack until a proper insturment module exists) and select only the
     # telescopes with the same camera type
     data = next(hessio_event_source(self.infile, max_events=1))
     camtypes = get_camera_types(data.inst)
     print_camera_types(data.inst, printer=self.log.info)
     group = camtypes.groups[self.telgroup]
     self._selected_tels = group['tel_id'].data
     self._base_tel = self._selected_tels[0]
     self.log.info("Telescope group %d: %s with %s camera", self.telgroup,
                   group[0]['tel_type'], group[0]['cam_type'])
     self.log.info("SELECTED TELESCOPES:{}".format(self._selected_tels))
     self.calibrator = CameraCalibrator(self.config, self)

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

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

class SingleTelEventDisplay(Tool):
    name = "ctapipe-display-single-tel"
    description = Unicode(__doc__)

    infile = Unicode(help="input file to read", default='').tag(config=True)
    tel = Int(help='Telescope ID to display', default=0).tag(config=True)
    channel = Integer(help="channel number to display", min=0, max=1).tag(
        config=True)
    write = Bool(help="Write out images to PNG files", default=False).tag(
        config=True)
    clean = Bool(help="Apply image cleaning", default=False).tag(config=True)
    hillas = Bool(help="Apply and display Hillas parametrization",
                  default=False).tag(config=True)
    samples = Bool(help="Show each sample", default=False).tag(config=True)
    display = Bool(help="Display results in interactive window",
                   default_value=True).tag(config=True)
    delay = Float(help='delay between events in s', default_value=0.01,
                  min=0.001).tag(config=True)
    progress = Bool(help='display progress bar', default_value=True).tag(
        config=True)

    aliases = Dict({'infile': 'EventFileReaderFactory.input_path',
                    'tel': 'SingleTelEventDisplay.tel',
                    'max-events': 'EventFileReaderFactory.max_events',
                    'channel': 'SingleTelEventDisplay.channel',
                    'write': 'SingleTelEventDisplay.write',
                    'clean': 'SingleTelEventDisplay.clean',
                    'hillas': 'SingleTelEventDisplay.hillas',
                    'samples': 'SingleTelEventDisplay.samples',
                    'display': 'SingleTelEventDisplay.display',
                    'delay': 'SingleTelEventDisplay.delay',
                    'progress': 'SingleTelEventDisplay.progress'
                    })

    classes = List([EventFileReaderFactory, CameraCalibrator])

    def setup(self):

        reader_factory = EventFileReaderFactory(None, self)
        reader_class = reader_factory.get_class()
        self.reader = reader_class(None, self)

        self.calibrator = CameraCalibrator(config=None, tool=self,
                                           origin=self.reader.origin)
        self.source = self.reader.read(allowed_tels=[self.tel, ])

        self.log.info('SELECTING EVENTS FROM TELESCOPE {}'.format(self.tel))

    def start(self):

        disp = None

        for event in tqdm(self.source,
                          desc='Tel{}'.format(self.tel),
                          total=self.reader.max_events,
                          disable=~self.progress):

            self.log.debug(event.trig)
            self.log.debug("Energy: {}".format(event.mc.energy))

            self.calibrator.calibrate(event)

            if disp is None:
                x, y = event.inst.pixel_pos[self.tel]
                focal_len = event.inst.optical_foclen[self.tel]
                geom = CameraGeometry.guess(x, y, focal_len)
                self.log.info(geom)
                disp = CameraDisplay(geom)
                # disp.enable_pixel_picker()
                disp.add_colorbar()
                if self.display:
                    plt.show(block=False)

            # display the event
            disp.axes.set_title('CT{:03d} ({}), event {:06d}'.format(
                self.tel, geom.cam_id, event.r0.event_id)
            )

            if self.samples:
                # display time-varying event
                data = event.dl0.tel[self.tel].pe_samples[self.channel]
                for ii in range(data.shape[1]):
                    disp.image = data[:, ii]
                    disp.set_limits_percent(70)
                    plt.suptitle("Sample {:03d}".format(ii))
                    if self.display:
                        plt.pause(self.delay)
                    if self.write:
                        plt.savefig('CT{:03d}_EV{:10d}_S{:02d}.png'
                                    .format(self.tel, event.r0.event_id, ii))
            else:
                # display integrated event:
                im = event.dl1.tel[self.tel].image[self.channel]

                if self.clean:
                    mask = tailcuts_clean(geom, im, picture_thresh=10,
                                          boundary_thresh=7)
                    im[~mask] = 0.0

                disp.image = im
#.........这里部分代码省略.........

The most basic pipeline, using no special features of the framework other 
than a for-loop. This is useful for debugging and profiling of speed.
"""

import sys

import numpy as np

from ctapipe.calib import CameraCalibrator
from ctapipe.io.hessio import hessio_event_source

if __name__ == '__main__':

    filename = sys.argv[1]

    source = hessio_event_source(filename, max_events=None)

    cal = CameraCalibrator(None, None)

    for data in source:

        print("EVENT: {}, ENERGY: {:.2f}, TELS:{}"
              .format(data.r0.event_id,
                      data.mc.energy,
                      len(data.dl0.tels_with_data))
              )

        cal.calibrate(data)

        # now the calibrated images are in data.dl1.tel[x].image

class DisplayDL1Calib(Tool):
    name = "DisplayDL1Calib"
    description = "Calibrate dl0 data to dl1, and plot the photoelectron " \
                  "images."

    telescope = Int(None, allow_none=True,
                    help='Telescope to view. Set to None to display all '
                         'telescopes.').tag(config=True)

    aliases = Dict(dict(f='EventFileReaderFactory.input_path',
                        r='EventFileReaderFactory.reader',
                        max_events='EventFileReaderFactory.max_events',
                        extractor='ChargeExtractorFactory.extractor',
                        window_width='ChargeExtractorFactory.window_width',
                        t0='ChargeExtractorFactory.t0',
                        window_shift='ChargeExtractorFactory.window_shift',
                        sig_amp_cut_HG='ChargeExtractorFactory.sig_amp_cut_HG',
                        sig_amp_cut_LG='ChargeExtractorFactory.sig_amp_cut_LG',
                        lwt='ChargeExtractorFactory.lwt',
                        clip_amplitude='CameraDL1Calibrator.clip_amplitude',
                        T='DisplayDL1Calib.telescope',
                        O='ImagePlotter.output_path'
                        ))
    flags = Dict(dict(D=({'ImagePlotter': {'display': True}},
                         "Display the photoelectron images on-screen as they "
                         "are produced.")
                      ))
    classes = List([EventFileReaderFactory,
                    ChargeExtractorFactory,
                    CameraDL1Calibrator,
                    ImagePlotter
                    ])

    def __init__(self, **kwargs):
        super().__init__(**kwargs)
        self.reader = None
        self.calibrator = None
        self.plotter = None

    def setup(self):
        kwargs = dict(config=self.config, tool=self)

        reader_factory = EventFileReaderFactory(**kwargs)
        reader_class = reader_factory.get_class()
        self.reader = reader_class(**kwargs)

        self.calibrator = CameraCalibrator(origin=self.reader.origin, **kwargs)

        self.plotter = ImagePlotter(**kwargs)

    def start(self):
        source = self.reader.read()
        for event in source:
            self.calibrator.calibrate(event)

            tel_list = event.r0.tels_with_data

            if self.telescope:
                if self.telescope not in tel_list:
                    continue
                tel_list = [self.telescope]
            for telid in tel_list:
                self.plotter.plot(event, telid)

    def finish(self):
        self.plotter.finish()