Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
input_file = resource_filename(
'eventio',
'resources/gamma_20deg_0deg_run103___cta-prod4-sst-astri_desert-2150m-Paranal-sst-astri.simtel.gz'
)
with EventIOFile(input_file) as f:
cameras = {}
for o in f:
if isinstance(o, CameraSettings):
cam_data = o.parse()
pix_type = 'square'
pix_rotation = 0 * u.deg
cameras[o.telescope_id] = CameraGeometry(
cam_id='CAM-{}'.format(o.telescope_id),
pix_id=np.arange(cam_data['n_pixels']),
pix_x=cam_data['pixel_x'] * u.m,
pix_y=cam_data['pixel_y'] * u.m,
pix_area=cam_data['pixel_area'] * u.m**2,
pix_type=pix_type,
cam_rotation=cam_data['cam_rot'] * u.rad,
pix_rotation=pix_rotation,
)
if isinstance(o, TelescopeData):
for subo in o:
if isinstance(subo, PhotoElectrons):
pe = subo.parse()
plt.figure()
pix_type = 'hexagonal'
if cam_data['n_pixels'] == 1855:
pix_rotation = 0 * u.deg
else:
pix_rotation = 30 * u.deg
elif cam_data['pixel_shape'][0] == -1:
if cam_data['n_pixels'] > 2000:
pix_type = 'square'
pix_rotation = 0 * u.deg
else:
pix_type = 'hexagonal'
pix_rotation = 0 * u.deg
cameras[o.telescope_id] = CameraGeometry(
cam_id='CAM-{}'.format(o.telescope_id),
pix_id=np.arange(cam_data['n_pixels']),
pix_x=cam_data['pixel_x'] * u.m,
pix_y=cam_data['pixel_y'] * u.m,
pix_area=cam_data['pixel_area'] * u.m**2,
pix_type=pix_type,
cam_rotation=cam_data['cam_rot'] * u.rad,
pix_rotation=pix_rotation,
)
if isinstance(o, Event):
for subo in o:
if isinstance(subo, TelescopeEvent):
for subsubo in subo:
if isinstance(subsubo, ADCSamples):
data = subsubo.parse()
def get_pixel_width(self, tel_id):
"""Guesstimate fov radius for telescope with id `tel_id`"""
# memoize fov calculation
if tel_id not in self.pixel_widths:
x, y = self.get_pixel_coords(tel_id)
self.pixel_widths[tel_id] = CameraGeometry.guess_pixel_width(x, y)
return self.pixel_widths[tel_id]
# Get calibrated image (low gain channel only)
pmt_signal = event.dl1.tel[tel_id].image[0]
if len(event.dl1.tel[tel_id].image) >1:
print(event.dl1.tel[tel_id].image[1][pmt_signal>100])
pmt_signal[pmt_signal > 100] = \
event.dl1.tel[tel_id].image[1][pmt_signal > 100]
# Create nominal system for the telescope (this should later used telescope
# pointing)
nom_system = NominalFrame(array_direction=array_pointing,
pointing_direction=array_pointing)
# Create camera system of all pixels
pix_x, pix_y = event.inst.pixel_pos[tel_id]
fl = event.inst.optical_foclen[tel_id]
if tel_id not in self.geoms:
self.geoms[tel_id] = CameraGeometry.guess(pix_x, pix_y,
event.inst.optical_foclen[
tel_id],
apply_derotation=False)
# Transform the pixels positions into nominal coordinates
camera_coord = CameraFrame(x=pix_x, y=pix_y, z=np.zeros(pix_x.shape) * u.m,
focal_length=fl,
rotation= -1* self.geoms[tel_id].cam_rotation)
nom_coord = camera_coord.transform_to(nom_system)
tx, ty, tz = event.inst.tel_pos[tel_id]
# ImPACT reconstruction is performed in the tilted system,
# so we need to transform tel positions
grd_tel = GroundFrame(x=tx, y=ty, z=tz)
tilt_tel = grd_tel.transform_to(tilted_system)
for tel_id in event.dl0.tels_with_data:
# Get calibrated image (low gain channel only)
pmt_signal = event.dl1.tel[tel_id].image[0]
# Create nominal system for the telescope (this should later used telescope
# pointing)
nom_system = NominalFrame(array_direction=array_pointing,
pointing_direction=array_pointing)
# Create camera system of all pixels
pix_x, pix_y = event.inst.pixel_pos[tel_id]
fl = event.inst.optical_foclen[tel_id]
if tel_id not in self.geoms:
self.geoms[tel_id] = CameraGeometry.guess(pix_x, pix_y,
event.inst.optical_foclen[
tel_id])
# Transform the pixels positions into nominal coordinates
camera_coord = CameraFrame(x=pix_x, y=pix_y, z=np.zeros(pix_x.shape) * u.m,
focal_length=fl,
rotation=-1 * self.geoms[tel_id].cam_rotation)
nom_coord = camera_coord.transform_to(nom_system)
tx, ty, tz = event.inst.tel_pos[tel_id]
# ImPACT reconstruction is performed in the tilted system,
# so we need to transform tel positions
grd_tel = GroundFrame(x=tx, y=ty, z=tz)
tilt_tel = grd_tel.transform_to(tilted_system)
# Clean image using split level cleaning
from ctapipe.instrument import CameraGeometry
from matplotlib import pyplot as plt
geom = CameraGeometry.from_name("LSTCam")
plt.figure(figsize=(8, 3))
plt.subplot(1, 2, 1)
plt.imshow(geom.neighbor_matrix, origin="lower")
plt.title("Pixel Neighbor Matrix")
plt.subplot(1, 2, 2)
plt.scatter(geom.pix_x, geom.pix_y)
plt.title("Pixel Positions")
plt.show()
def build_camera(cam_settings, pixel_settings, telescope):
pixel_shape = cam_settings["pixel_shape"][0]
try:
pix_type, pix_rotation = CameraGeometry.simtel_shape_to_type(pixel_shape)
except ValueError:
warnings.warn(
f"Unkown pixel_shape {pixel_shape} for camera_type {telescope.camera_name}",
UnknownPixelShapeWarning,
)
pix_type = "hexagon"
pix_rotation = "0d"
geometry = CameraGeometry(
telescope.camera_name,
pix_id=np.arange(cam_settings["n_pixels"]),
pix_x=u.Quantity(cam_settings["pixel_x"], u.m),
pix_y=u.Quantity(cam_settings["pixel_y"], u.m),
pix_area=u.Quantity(cam_settings["pixel_area"], u.m ** 2),
pix_type=pix_type,
pix_rotation=pix_rotation,
# instead of blindly enumerating all pixels, let's instead
# store a list of all valid -- i.e. picked by the mask -- 2D
# indices
for i, row in enumerate(square_mask):
for j, val in enumerate(row):
if val is True:
ids.append((i, j))
# the area of the pixels (note that this is still a deformed
# image)
pix_area = (
np.ones_like(grid_x) * (x_edges[1] - x_edges[0]) * (y_edges[1] - y_edges[0])
)
# creating a new geometry object with the attributes we just determined
new_geom = CameraGeometry(
camera_name=geom.camera_name + "_rect",
pix_id=ids, # this is a list of all the valid coordinate pairs now
pix_x=u.Quantity(grid_x.ravel(), u.meter),
pix_y=u.Quantity(grid_y.ravel(), u.meter),
pix_area=pix_area * u.meter ** 2,
neighbors=geom.neighbors,
pix_type="rectangular",
apply_derotation=False,
)
# storing the pixel mask for later use
new_geom.mask = square_mask
# create a transfer map by enumerating all pixel positions in a 2D histogram
hex_to_rect_map = np.histogramdd(
[rot_y.to_value(u.m), rot_x.to_value(u.m)],
'eventio',
'resources/gamma_20deg_0deg_run103___cta-prod4-sst-astri_desert-2150m-Paranal-sst-astri.simtel.gz'
)
with EventIOFile(input_file) as f:
cameras = {}
for o in f:
if isinstance(o, CameraSettings):
cam_data = o.parse()
if cam_data['pixel_shape'][0] == -1:
pixel_shape = 'hexagonal' if cam_data['n_pixels'] < 2000 else 'square'
else:
pixel_shape = 'square' if cam_data['pixel_shape'][0] else 'hexagonal'
cameras[o.telescope_id] = CameraGeometry(
cam_id='CAM-{}'.format(o.telescope_id),
pix_id=np.arange(cam_data['n_pixels']),
pix_x=cam_data['pixel_x'] * u.m,
pix_y=cam_data['pixel_y'] * u.m,
pix_area=cam_data['pixel_area'] * u.m**2,
pix_type=pixel_shape,
cam_rotation=cam_data['cam_rot'] * u.rad,
)
if isinstance(o, Event):
assert len(cameras) > 0
for subo in o:
if isinstance(subo, TelescopeEvent):
for subsubo in subo:
if isinstance(subsubo, ADCSums):
data = subsubo.parse()