Skip to content

Commit

Permalink
added details of the model
Browse files Browse the repository at this point in the history
  • Loading branch information
clara-escanuela committed Mar 16, 2023
1 parent c2a40d0 commit e994d6b
Show file tree
Hide file tree
Showing 2 changed files with 16 additions and 10 deletions.
17 changes: 10 additions & 7 deletions ctapipe/image/ellipsoid.py
Original file line number Diff line number Diff line change
Expand Up @@ -52,7 +52,7 @@ def create_initial_guess(geometry, image, pdf):
scale_skew = hillas.length.to_value(u.m) / np.sqrt(1 - 2 * delta ** 2 / np.pi)

if pdf == "Gaussian":
initial_guess["amplitude"] = hillas.intensity/(np.sqrt(2*np.pi)*hillas.length.value*hillas.width.value)
initial_guess["amplitude"] = hillas.intensity/(2*np.pi*hillas.length.value*hillas.width.value)
if pdf == "Cauchy":
initial_guess["amplitude"] = hillas.intensity/(np.pi*scale_skew*hillas.width.value)
if pdf == "Skewed":
Expand Down Expand Up @@ -115,7 +115,7 @@ def boundaries(geometry, image, cleaning_mask, clean_row_mask, x0, pdf):
pix_area = geometry.pix_area.value[0]
area = pix_area * np.count_nonzero(row_image)

leakage = leakage_parameters(geometry, image, clean_row_mask)
leakage = leakage_parameters(geometry, image, clean_row_mask) #TODO: length boundaries dependent on leakage or distance of centroid to centre of camera
fract_pix_border = leakage.pixels_width_2
fract_int_border = leakage.intensity_width_2

Expand All @@ -124,12 +124,12 @@ def boundaries(geometry, image, cleaning_mask, clean_row_mask, x0, pdf):
longitudinal = delta_x * np.cos(x0["psi"].value) + delta_y * np.sin(x0["psi"].value)
transverse = delta_x * -np.sin(x0["psi"].value) + delta_y * np.cos(x0["psi"].value)

cogx_min, cogx_max = np.min(geometry.pix_x.value[cleaned_image>0]), np.max(geometry.pix_x.value[cleaned_image>0])
cogy_min, cogy_max = np.min(geometry.pix_y.value[cleaned_image>0]), np.max(geometry.pix_y.value[cleaned_image>0])
cogx_min, cogx_max = np.min(geometry.pix_x.value[cleaned_image>0.3*np.max(cleaned_image)]), np.max(geometry.pix_x.value[cleaned_image>0.3*np.max(cleaned_image)])
cogy_min, cogy_max = np.min(geometry.pix_y.value[cleaned_image>0.3*np.max(cleaned_image)]), np.max(geometry.pix_y.value[cleaned_image>0.3*np.max(cleaned_image)])

x_dis = np.max(longitudinal[row_image>0]) - np.min(longitudinal[row_image>0])
y_dis = np.max(transverse[row_image>0]) - np.min(transverse[row_image>0])
length_min, length_max = np.sqrt(pix_area), x_dis/(1 - fract_pix_border)
length_min, length_max = np.sqrt(pix_area), x_dis

width_min, width_max = np.sqrt(pix_area), y_dis
scale = length_min/ np.sqrt(1 - 2 / np.pi)
Expand All @@ -142,7 +142,7 @@ def boundaries(geometry, image, cleaning_mask, clean_row_mask, x0, pdf):
if pdf == "Cauchy":
return [(cogx_min, cogx_max), (cogy_min, cogy_max), (-np.pi/2, np.pi/2), (length_min, length_max), (width_min, width_max), (skew_min, skew_max), (0, np.sum(row_image) * 1/scale * 1/(np.pi*width_min))]

def image_fit_parameters(geom, image, n, cleaned_mask, spe_width, pedestal, pdf, bounds=None):
def image_fit_parameters(geom, image, n, cleaned_mask, spe_width, pedestal, pdf="Cauchy", bounds=None):
"""
Computes image parameters for a given shower image.
Expand All @@ -165,7 +165,7 @@ def image_fit_parameters(geom, image, n, cleaned_mask, spe_width, pedestal, pdf,
pedestal: ndarray
Width of pedestal (:math:`σ_p`).
pdf : str
name of the prob distrib to use for the fit
name of the prob distrib to use for the fit, options = "Gaussian", "Cauchy", "Skewed"
Returns
-------
Expand Down Expand Up @@ -216,6 +216,9 @@ def fit_gauss(cog_x, cog_y, psi, length, width, amplitude):

bounds = boundaries(geom, image, cleaned_mask, mask, x0, pdf)

if bounds == None:
bounds = boundaries(geom, image, cleaned_mask, mask, x0, pdf)
m.limits = bounds
if bounds != None:
m.limits = bounds

Expand Down
9 changes: 6 additions & 3 deletions ctapipe/image/toymodel.py
Original file line number Diff line number Diff line change
Expand Up @@ -247,7 +247,7 @@ def expected_signal(self, camera, intensity):

class Gaussian(ImageModel):
@u.quantity_input(x=u.m, y=u.m, length=u.m, width=u.m)
def __init__(self, x, y, length, width, psi, ampl):
def __init__(self, x, y, length, width, psi, amplitude):
"""Create 2D Gaussian model for a shower image in a camera.
Parameters
Expand All @@ -260,6 +260,7 @@ def __init__(self, x, y, length, width, psi, ampl):
length of shower (major axis)
psi : convertable to `astropy.coordinates.Angle`
rotation angle about the centroid (0=x-axis)
amplitude : normalization amplitude
Returns
-------
Expand All @@ -271,15 +272,15 @@ def __init__(self, x, y, length, width, psi, ampl):
self.width = width
self.length = length
self.psi = psi
self.ampl = ampl
self.amplitude = amplitude

@u.quantity_input(x=u.m, y=u.m)
def pdf(self, x, y):
"""2d probability for photon electrons in the camera plane"""
long = (x.to_value(u.m) - self.x.to_value(u.m)) * np.cos(Angle(self.psi)) + (y.to_value(u.m) - self.y.to_value(u.m)) * np.sin(Angle(self.psi))
trans = (x.to_value(u.m) - self.x.to_value(u.m)) * -np.sin(Angle(self.psi)) + (y.to_value(u.m) - self.y.to_value(u.m)) * np.cos(Angle(self.psi))

gaussian_pdf = self.ampl * np.exp(-1/2*(long)**2/self.length.to_value(u.m)**2 - 1/2*(trans)**2/self.width.to_value(u.m)**2)
gaussian_pdf = self.amplitude * np.exp(-1/2*(long)**2/self.length.to_value(u.m)**2 - 1/2*(trans)**2/self.width.to_value(u.m)**2)

return gaussian_pdf

Expand All @@ -306,6 +307,7 @@ def __init__(self, x, y, length, width, psi, skewness, amplitude):
length of shower (major axis)
psi : convertable to `astropy.coordinates.Angle`
rotation angle about the centroid (0=x-axis)
amplitude : normalization amplitude
Returns
-------
Expand Down Expand Up @@ -393,6 +395,7 @@ def __init__(self, x, y, length, width, psi, skewness, amplitude):
length of shower (major axis)
psi : convertable to `astropy.coordinates.Angle`
rotation angle about the centroid (0=x-axis)
amplitude : normalization amplitude
Returns
-------
Expand Down

0 comments on commit e994d6b

Please sign in to comment.