Python CameraCalibrator.calibrate方法代码示例

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

# 需要导入模块: from ctapipe.calib import CameraCalibrator [as 别名]
# 或者: from ctapipe.calib.CameraCalibrator import calibrate [as 别名]
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
#.........这里部分代码省略.........

# 需要导入模块: from ctapipe.calib import CameraCalibrator [as 别名]
# 或者: from ctapipe.calib.CameraCalibrator import calibrate [as 别名]
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()

# 需要导入模块: from ctapipe.calib import CameraCalibrator [as 别名]
# 或者: from ctapipe.calib.CameraCalibrator import calibrate [as 别名]
class ImageSumDisplayerTool(Tool):
    description = Unicode(__doc__)
    name = "ctapipe-display-imagesum"

    infile = Unicode(
        help='input simtelarray file',
        default="/Users/kosack/Data/CTA/Prod3/gamma.simtel.gz"
    ).tag(config=True)

    telgroup = Integer(
        help='telescope group number', default=1
    ).tag(config=True)

    max_events = Integer(
        help='stop after this many events if non-zero', default_value=0, min=0
    ).tag(config=True)

    output_suffix = Unicode(
        help='suffix (file extension) of output '
             'filenames to write images '
             'to (no writing is done if blank). '
             'Images will be named [EVENTID][suffix]',
        default_value=""
    ).tag(config=True)

    aliases = Dict({
        'infile': 'ImageSumDisplayerTool.infile',
        'telgroup': 'ImageSumDisplayerTool.telgroup',
        'max-events': 'ImageSumDisplayerTool.max_events',
        'output-suffix': 'ImageSumDisplayerTool.output_suffix'
    })

    classes = List([CameraCalibrator, SimTelEventSource])

    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 start(self):
        geom = None
        imsum = None
        disp = None

        for event in self.reader:

            self.calibrator.calibrate(event)

            if geom is None:
                geom = event.inst.subarray.tel[self._base_tel].camera
                imsum = np.zeros(shape=geom.pix_x.shape, dtype=np.float)
                disp = CameraDisplay(geom, title=geom.cam_id)
                disp.add_colorbar()
                disp.cmap = 'viridis'

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

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

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

            if self.output_suffix is not "":
                filename = "{:020d}{}".format(
                    event.r0.event_id, self.output_suffix
                )
                self.log.info(f"saving: '{filename}'")
#.........这里部分代码省略.........

# 需要导入模块: from ctapipe.calib import CameraCalibrator [as 别名]
# 或者: from ctapipe.calib.CameraCalibrator import calibrate [as 别名]
class MuonDisplayerTool(Tool):
    name = 'ctapipe-display-muons'
    description = t.Unicode(__doc__)

    events = t.Unicode("",
                       help="input event data file").tag(config=True)

    outfile = t.Unicode("muons.hdf5", help='HDF5 output file name').tag(
        config=True)

    display = t.Bool(
        help='display the camera events', default=False
    ).tag(config=True)

    classes = t.List([
        CameraCalibrator, EventSourceFactory
    ])

    aliases = t.Dict({
        'input': 'MuonDisplayerTool.events',
        'outfile': 'MuonDisplayerTool.outfile',
        'display': 'MuonDisplayerTool.display',
        'max_events': 'EventSourceFactory.max_events',
        'allowed_tels': 'EventSourceFactory.allowed_tels',
    })

    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 start(self):

        numev = 0
        self.num_muons_found = defaultdict(int)

        for event in tqdm(self.source, desc='detecting muons'):

            self.calib.calibrate(event)
            muon_evt = analyze_muon_event(event)

            if numev == 0:
                _exclude_some_columns(event.inst.subarray, self.writer)

            numev += 1

            if not muon_evt['MuonIntensityParams']:
                # No telescopes  contained a good muon
                continue
            else:
                if self.display:
                    plot_muon_event(event, muon_evt)

                for tel_id in muon_evt['TelIds']:
                    idx = muon_evt['TelIds'].index(tel_id)
                    intens_params = muon_evt['MuonIntensityParams'][idx]

                    if intens_params is not None:
                        ring_params = muon_evt['MuonRingParams'][idx]
                        cam_id = str(event.inst.subarray.tel[tel_id].camera)
                        self.num_muons_found[cam_id] += 1
                        self.log.debug("INTENSITY: %s", intens_params)
                        self.log.debug("RING: %s", ring_params)
                        self.writer.write(table_name=cam_id,
                                          containers=[intens_params,
                                                      ring_params])

                self.log.info(
                    "Event Number: %d, found %s muons",
                    numev, dict(self.num_muons_found)
                )

    def finish(self):
        Provenance().add_output_file(self.outfile,
                                     role='dl1.tel.evt.muon')
        self.writer.close()

# 需要导入模块: from ctapipe.calib import CameraCalibrator [as 别名]
# 或者: from ctapipe.calib.CameraCalibrator import calibrate [as 别名]
class DisplayDL1Calib(Tool):
    name = "ctapipe-display-dl1"
    description = __doc__

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

    extractor_product = tool_utils.enum_trait(
        ImageExtractor,
        default='NeighborPeakWindowSum'
    )

    aliases = Dict(
        dict(
            max_events='EventSource.max_events',
            extractor='DisplayDL1Calib.extractor_product',
            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(
        [
            EventSource,
            CameraDL1Calibrator,
            ImagePlotter
        ] + tool_utils.classes_with_traits(ImageExtractor)
    )

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

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

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

        self.plotter = ImagePlotter(parent=self)

    def start(self):
        for event in self.eventsource:
            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()

# 需要导入模块: from ctapipe.calib import CameraCalibrator [as 别名]
# 或者: from ctapipe.calib.CameraCalibrator import calibrate [as 别名]
"""
very simple example that loads a single event into memory, for exploration
purposes
"""
import sys

from ctapipe.calib import CameraCalibrator
from ctapipe.io import event_source
from ctapipe.utils import get_dataset_path

if __name__ == '__main__':

    calib = CameraCalibrator(r1_product="HESSIOR1Calibrator")

    if len(sys.argv) >= 2:
        filename = sys.argv[1]
    else:
        filename = get_dataset_path("gamma_test_large.simtel.gz")

    with event_source(filename, max_events=1) as source:
        for event in source:
            calib.calibrate(event)

    print(event)

# 需要导入模块: from ctapipe.calib import CameraCalibrator [as 别名]
# 或者: from ctapipe.calib.CameraCalibrator import calibrate [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.calib import CameraCalibrator [as 别名]
# 或者: from ctapipe.calib.CameraCalibrator import calibrate [as 别名]
    'ASTRICam': (5, 7, 2),
}


input_url = get_dataset_path('gamma_test_large.simtel.gz')
event_source = EventSourceFactory.produce(input_url=input_url)

calibrator = CameraCalibrator(
    eventsource=event_source,
)

reco = HillasReconstructor()

for event in event_source:
    print('Event', event.count)
    calibrator.calibrate(event)

    # mapping of telescope_id to parameters for stereo reconstruction
    hillas_containers = {}
    pointing_azimuth = {}
    pointing_altitude = {}
    time_gradients = {}

    for telescope_id, dl1 in event.dl1.tel.items():
        camera = event.inst.subarray.tels[telescope_id].camera

        image = dl1.image[0]
        peakpos = dl1.peakpos[0]

        # cleaning
        boundary, picture, min_neighbors = cleaning_level[camera.cam_id]

# 需要导入模块: from ctapipe.calib import CameraCalibrator [as 别名]
# 或者: from ctapipe.calib.CameraCalibrator import calibrate [as 别名]
class SingleTelEventDisplay(Tool):
    name = "ctapipe-display-televents"
    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': 'SingleTelEventDisplay.infile',
        'tel': 'SingleTelEventDisplay.tel',
        'max-events': 'EventSource.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([EventSource, CameraCalibrator])

    def __init__(self, **kwargs):
        super().__init__(**kwargs)

    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, eventsource=self.event_source
        )

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

    def start(self):

        disp = None

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

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

            self.calibrator.calibrate(event)

            if disp is None:
                geom = event.inst.subarray.tel[self.tel].camera
                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].waveform[self.channel]
                for ii in range(data.shape[1]):
                    disp.image = data[:, ii]
                    disp.set_limits_percent(70)
                    plt.suptitle(f"Sample {ii:03d}")
                    if self.display:
                        plt.pause(self.delay)
                    if self.write:
                        plt.savefig(
#.........这里部分代码省略.........

# 需要导入模块: from ctapipe.calib import CameraCalibrator [as 别名]
# 或者: from ctapipe.calib.CameraCalibrator import calibrate [as 别名]
class ImageSumDisplayerTool(Tool):
    description = Unicode(__doc__)
    name = "ctapipe-image-sum-display"

    infile = Unicode(help='input simtelarray file',
                     default="/Users/kosack/Data/CTA/Prod3/gamma.simtel.gz"
                     ).tag(config=True)
    telgroup = Integer(help='telescope group number', default=1).tag(
        config=True)

    max_events = Integer(help='stop after this many events if non-zero',
                         default_value=0, min=0).tag(config=True)

    output_suffix=Unicode(help='suffix (file extension) of output '
                               'filenames to write images '
                               'to (no writing is done if blank). '
                               'Images will be named [EVENTID][suffix]' ,
                          default_value="").tag(config=True)

    aliases = Dict({'infile': 'ImageSumDisplayerTool.infile',
                    'telgroup': 'ImageSumDisplayerTool.telgroup',
                    'max-events': 'ImageSumDisplayerTool.max_events',
                    'output-suffix': 'ImageSumDisplayerTool.output_suffix'})
    classes = List([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
        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 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.calib import CameraCalibrator [as 别名]
# 或者: from ctapipe.calib.CameraCalibrator import calibrate [as 别名]
def extract_images(simtel_file_path,
                   tel_id_filter_list=None,
                   event_id_filter_list=None,
                   output_directory=None,
                   crop=True):

    # 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 (ASTRI ONLY) ###################

                    # TODO

                    #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 != "rectangular") or (geom.cam_id not in ("ASTRICam", "ASTRI")):
                        raise ValueError("Telescope {}: error (the input image is not a valide ASTRI 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,:] <= ASTRI_CAM_CHANNEL_THRESHOLD] = 0
                    calibrated_image[0, calibrated_image[0,:] >  ASTRI_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) ##########################

                    pe_image_2d = geometry_converter.astri_to_2d_array(pe_image, crop=crop)
                    calibrated_image_2d = geometry_converter.astri_to_2d_array(calibrated_image, crop=crop)
                    uncalibrated_image_2d = geometry_converter.astri_to_3d_array(uncalibrated_image, crop=crop)
                    pedestal_2d = geometry_converter.astri_to_3d_array(pedestal, crop=crop)
                    gains_2d = geometry_converter.astri_to_3d_array(gain, crop=crop)
                    pixel_pos_2d = geometry_converter.astri_to_3d_array(pixel_pos, crop=crop)

                    #print(pe_image_2d.shape)
                    #print(calibrated_image_2d.shape)
                    #print(uncalibrated_image_2d.shape)
                    #print(pedestal_2d.shape)
                    #print(gains_2d.shape)
                    #print(pixel_pos_2d.shape)
                    #sys.exit(0)

                    # GET PIXEL MASK ##########################################

                    pixel_mask = geometry_converter.astri_pixel_mask(crop)  # 1 for pixels with actual data, 0 for virtual (blank) pixels

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

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

# 需要导入模块: from ctapipe.calib import CameraCalibrator [as 别名]
# 或者: from ctapipe.calib.CameraCalibrator import calibrate [as 别名]
class SimpleEventWriter(Tool):
    name = 'ctapipe-simple-event-writer'
    description = Unicode(__doc__)

    infile = Unicode(help='input file to read', default='').tag(config=True)
    outfile = Unicode(help='output file name', default_value='output.h5').tag(config=True)
    progress = Bool(help='display progress bar', default_value=True).tag(config=True)

    aliases = Dict({
        'infile': 'EventSourceFactory.input_url',
        'outfile': 'SimpleEventWriter.outfile',
        'max-events': 'EventSourceFactory.max_events',
        'progress': 'SimpleEventWriter.progress'
    })
    classes = List([EventSourceFactory, CameraCalibrator, CutFlow])

    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 start(self):
        self.log.info('Loop on events...')

        for event in tqdm(
                self.event_source,
                desc='EventWriter',
                total=self.event_source.max_events,
                disable=~self.progress):

            self.event_cutflow.count('no_sel')
            self.calibrator.calibrate(event)

            for tel_id in event.dl0.tels_with_data:
                self.image_cutflow.count('no_sel')

                camera = event.inst.subarray.tel[tel_id].camera
                dl1_tel = event.dl1.tel[tel_id]

                # Image cleaning
                image = dl1_tel.image[0]  # Waiting for automatic gain selection
                mask = tailcuts_clean(camera, image, picture_thresh=10, boundary_thresh=5)
                cleaned = image.copy()
                cleaned[~mask] = 0

                # Preselection cuts
                if self.image_cutflow.cut('n_pixel', cleaned):
                    continue
                if self.image_cutflow.cut('image_amplitude', np.sum(cleaned)):
                    continue

                # Image parametrisation
                params = hillas_parameters(camera, cleaned)

                # Save Ids, MC infos and Hillas informations
                self.writer.write(camera.cam_id, [event.r0, event.mc, params])

    def finish(self):
        self.log.info('End of job.')

        self.image_cutflow()
        self.event_cutflow()
        self.writer.close()