Skip to content

pyhiperta.hillas

Hillas (and some additional) parameters computation from integrated waveforms.

Functions:

Name Description
hillas

Compute hillas parameters (ellipse position, orientation, skewness, kurtosis) from shower images.

hillas 

hillas(waveform, position, nan_threshold=0.0001)

Compute hillas parameters (ellipse position, orientation, skewness, kurtosis) from shower images.

The Hillas parametrization consists in estimating the ellipse formed by a gamma-ray shower over the background noise in the telescope integrated and cleaned images (R1 data level).

The ellipse parameters are found by performing a principal component analysis (PCA) of the images pixels weighted by their charge. The ellipse center is defined as the image center of gravity ; the ellipse axis and angle with respect to the x-axis are given by the PCA eigenvectors ; the ellipse length and width are defined as twice the standard deviation along the ellipse axis. In addition, the skewness and kurtosis of the image are computed along the major axis of the ellipse.

This method works with a single image, of with a batch of images. The implementation first computes the 1st and 2nd order image moments to build the covariance matrix (PCA is computed via the covariance matrix). The ellipse parameters are deduced from the covariance matrix and finally the skewness and kurtosis are computed as 3rd and 4rd along the ellipse major axis.

This function also returns the longitudinal coordinate of the pixels along the ellipsis greater axis since it is usefull for timing parameters computation.

Notes

A gamma-ray shower is radially symmetric in the direction of arrival, therefore the distribution along the small ellipse axis should be gaussian and the skewness and kurtosis are not computed along this axis.

An ellipse with given orientation and skewness is equivalent to another ellipse with orientation + pi and opposite skewness (psi + pi, -skewness). This method consistently returns the orientation in the ]-pi/2, pi/2[ range and corresponding skewness.

The skewness and kurtosis are computed using the 3rd and 4th moments, which are biased estimators.

Hillas, A. (1985). Cerenkov Light Images of EAS Produced by Primary Gamma Rays and by Nuclei. In 19th International Cosmic Ray Conference (ICRC19), Volume 3 (pp. 445).

Parameters:

Name Type Description Default
waveform ndarray

1D shower image(s). If a single image is provided the shape must be (N_pixels,). If a batch of images is provided, the shape should be (N_batch, N_pixels).

 required
position ndarray

Image pixel position in (x, y) coordinates. The shape must be (2, N_pixels): position[0,:] is the x coordinate, and position[1, :] the y coordinate, of the image pixels.

 required
nan_threshold float

Threshold to set some computed values to NaN instead of leaving them with problematic values. Images with intensity smaller than the threshold will have all output set to NaN. Images with lengths smaller than the threshold will have length, skewness and kurtosis set to NaN.

 0.0001

Returns:

Name Type Description
hillas_parameters ndarray

Image(s) hillas parameters. This is a 1D array with shape (10,) if a single image was provided, or a 2D array with shape (N_batch, 10) if a batch of images was provided. The parameters are ordered like so: hillas_parameters[..., 0]: intensity hillas_parameters[..., 1]: center of gravity x coordinate hillas_parameters[..., 2]: center of gravity y coordinate hillas_parameters[..., 3]: center of gravity radial coordinate radius (r) hillas_parameters[..., 4]: center of gravity radial coordinate angle (phi) hillas_parameters[..., 5]: ellipse length hillas_parameters[..., 6]: ellipse width hillas_parameters[..., 7]: ellipse angle with respect to x-axis (psi) in ]-pi/2, pi/2[ hillas_parameters[..., 8]: ellipse skewness hillas_parameters[..., 9]: ellipse kurtosis
longitudinal_coordinates ndarray

Coordinates of the pixels along the long axis of the ellipsis.
Source code in src/pyhiperta/hillas.py 
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
def hillas(
    waveform: np.ndarray, position: np.ndarray, nan_threshold: float = 1.0e-4
) -> tuple[np.ndarray, np.ndarray]:
    """Compute hillas parameters (ellipse position, orientation, skewness, kurtosis) from shower images.

    The Hillas parametrization consists in estimating the ellipse formed by a gamma-ray shower over the
    background noise in the telescope integrated and cleaned images (R1 data level).

    The ellipse parameters are found by performing a principal component analysis (PCA) of the images pixels
    weighted by their charge. The ellipse center is defined as the image center of gravity ; the ellipse
    axis and angle with respect to the x-axis are given by the PCA eigenvectors ; the ellipse length and
    width are defined as twice the standard deviation along the ellipse axis.
    In addition, the skewness and kurtosis of the image are computed along the major axis of the ellipse.

    This method works with a single image, of with a batch of images. The implementation first computes
    the 1st and 2nd order image moments to build the covariance matrix (PCA is computed via the
    covariance matrix). The ellipse parameters are deduced from the covariance matrix and finally the
    skewness and kurtosis are computed as 3rd and 4rd along the ellipse major axis.

    This function also returns the longitudinal coordinate of the pixels along the ellipsis greater axis
    since it is usefull for timing parameters computation.

    Notes
    -----
    A gamma-ray shower is radially symmetric in the direction of arrival, therefore the distribution
    along the small ellipse axis should be gaussian and the skewness and kurtosis are not computed
    along this axis.

    An ellipse with given orientation and skewness is equivalent to another ellipse with
    orientation + pi and opposite skewness (psi + pi, -skewness). This method consistently returns
    the orientation in the ]-pi/2, pi/2[ range and corresponding skewness.

    The skewness and kurtosis are computed using the 3rd and 4th moments, which are biased estimators.

    Hillas, A. (1985). Cerenkov Light Images of EAS Produced by Primary Gamma Rays and by Nuclei.
    In 19th International Cosmic Ray Conference (ICRC19), Volume 3 (pp. 445).

    Parameters
    ----------
    waveform : np.ndarray
        1D shower image(s). If a single image is provided the shape must be (N_pixels,). If a
        batch of images is provided, the shape should be (N_batch, N_pixels).
    position : np.ndarray
        Image pixel position in (x, y) coordinates. The shape must be (2, N_pixels): position[0,:]
        is the x coordinate, and position[1, :] the y coordinate, of the image pixels.
    nan_threshold : float
        Threshold to set some computed values to NaN instead of leaving them with problematic values.
        Images with intensity smaller than the threshold will have all output set to NaN.
        Images with lengths smaller than the threshold will have length, skewness and kurtosis set
        to NaN.

    Returns
    -------
    hillas_parameters : np.ndarray
        Image(s) hillas parameters. This is a 1D array with shape (10,) if a single image was
        provided, or a 2D array with shape (N_batch, 10) if a batch of images was provided. The
        parameters are ordered like so:
        hillas_parameters[..., 0]: intensity
        hillas_parameters[..., 1]: center of gravity x coordinate
        hillas_parameters[..., 2]: center of gravity y coordinate
        hillas_parameters[..., 3]: center of gravity radial coordinate radius (r)
        hillas_parameters[..., 4]: center of gravity radial coordinate angle (phi)
        hillas_parameters[..., 5]: ellipse length
        hillas_parameters[..., 6]: ellipse width
        hillas_parameters[..., 7]: ellipse angle with respect to x-axis (psi) in ]-pi/2, pi/2[
        hillas_parameters[..., 8]: ellipse skewness
        hillas_parameters[..., 9]: ellipse kurtosis
    longitudinal_coordinates : np.ndarray
        Coordinates of the pixels along the long axis of the ellipsis.
    """
    # some computation could have division by 0 or sqrt(negative_value)
    # we don't want warnings and actually want the NaNs in the result, so
    # simply disable these warnings for this block
    with np.errstate(divide="ignore", invalid="ignore"):
        # intensity is simply the sum of all pixels (for each image in batch)
        intensity = waveform.sum(axis=-1)
        # If intensity is too small, set it to NaN (all results will be meaningless)
        intensity = nanify(intensity, intensity < nan_threshold)

        # The PCA can be computed with the singular value decomposition (SVD) or covariance matrix. Here we
        # use the covariance matrix.
        # We want to compute the variance of x and y, as well as the covariance(x,y) along the distribution
        # given by the image charge (pixel values). The distribution will be normalized by the total charge.
        # To do so we will compute the 1st and 2nd order moments of x and y.
        # (Notation: E[x] is expected value of x, w is pixel value)

        # Compute all moments: E[w*x], E[w*y], E[w*x**2], E[w*y**2], E[w*x*y] and normalize
        # (stack all positions powers: x, y, x**2, y**2, x*y) and broadcast multiply with waveform
        moments = (
            waveform[..., np.newaxis, :]  # shape (N_batch, 1, N_pixels) to broadcast for all moments
            * np.stack([position[0], position[1], position[0] ** 2, position[1] ** 2, position.prod(axis=0)])
        ).sum(axis=-1) / intensity[..., np.newaxis]

        # E[x]**2 and E[y]**2 (will be re-used for center of gravity radial coordinates computation)
        E_x_y_square = moments[..., 0:2] ** 2
        # E[w*x**2] - E[w*x]**2 and E[w*y**2] - E[w*y]**2 (shape: (N_batch, 2))
        variance = moments[..., 2:4] - E_x_y_square
        # E[w*x*y] - E[w*x]*E[w*y]  (shape: (N_batch,))
        covariance = moments[..., 4] - moments[..., 0:2].prod(axis=-1)

        # The covariance matrix would be: (Var(x)    Cov(x,y))
        #                                 (Cov(x,y)  Var(y)  )
        # Now we compute the eigenvalues lambda as the roots of the characteristic polynomial
        #     (lambda - Var(x)) * (lambda - Var(y)) - Cov(x,y)**2 = 0
        varx_plus_vary = variance.sum(axis=-1)
        char_pol_delta = np.sqrt(varx_plus_vary**2 + 4.0 * (covariance**2 - variance.prod(axis=-1)))
        eigenValueHigh = (varx_plus_vary + char_pol_delta) / 2.0
        eigenValueLow = (varx_plus_vary - char_pol_delta) / 2.0

        # The length and width are twice the square roots of the eigenvalues (*2 is performed later)
        # (<=> std along the major and minor axis)
        # Note: an ellipse parametrized by A and B is defined by
        # x^T (A   B/2) x = 1
        #     (B/2 A  )
        # the semi-axis are given by the square root of the ellipse matrix eigenvalues.
        length = np.sqrt(eigenValueHigh)
        width = np.sqrt(eigenValueLow)
        # If the length is too small, set it to NaN to avoid crazy values later on.
        length = nanify(length, length < nan_threshold)

        # The ellipse axis are along the eigenvectors of the covariance matrix, and the ellipse
        # orientation angle (psi) is defined as the angle between the major axis and the x-axis.
        # However, we can define the axis with 2 possible directions: for a given eigenvector e,
        # -e is also an eigenvector with the same direction, but reversed sense.
        # The convention here is to chose the eigenvector so that psi is in ]-pi/2, pi/2[
        # Note: eigenvectors e are defined by lamba e = (Var(x)    Cov(x, y)) e
        #                                               (Cov(x, y) Var(y)   )
        # To solve this system e[0] and e[1] must respect a fixed ratio, that is most simply
        # expressed as e = (Cov(x, y), (lambda - Var(x))) or e = (lambda - Var(y), Cov(x,y))
        # We are free to chose the direction of e, so we chose the one giving psi in ]-pi/2,pi/2[
        psi = np.arctan2(covariance, eigenValueHigh - variance[..., 1])

        # Finally, we compute the 3rd and 4th order moments along the ellipse major axis
        # First we compute the pixels coordinates in the PCA coordinates: (cos(psi),  sin(psi)) (p - cog)
        #                                                                 (-sin(psi), cos(psi))
        # but we only care about the major axis, so only the 1st coordinate is computed
        # longitudinal shape (N_batch, N_pixels) (the sum is along coordinate axis to do p * cos + p * sin)
        longitudinal = (
            (position - moments[..., 0:2, np.newaxis])  # shape (N_batch, 2, N_pixels)
            * np.stack([np.cos(psi), np.sin(psi)], axis=-1)[..., np.newaxis]  # shape (N_batch, 2, 1)
        ).sum(axis=-2)

        # we want to compute longitudinal **3 and **4
        # to do this with broadcasting, we build an array containing [3, 4] that broadcasts
        # to (2, *longitudinal.shape). See np.power()
        exponents = np.array([3, 4], dtype=np.float32)[:, *([np.newaxis] * longitudinal.ndim)]
        # Compute 3rd and 4th moments E[w * l**3] and E[w * l**4], shape: (2, N_batch)
        longitudinal_moments = ((longitudinal**exponents) * waveform).sum(axis=-1)
        # Normalize by ellipse length
        skewness, kurtosis = longitudinal_moments / ((length ** exponents.squeeze(-1)) * intensity)

        # Radial coordinates of the center of gravity
        r = np.sqrt(E_x_y_square.sum(axis=-1))
        phi = np.arctan2(moments[..., 1], moments[..., 0])

    # Stack everything in columns so it is ready for scikit-learn
    return np.stack(
        [
            intensity,
            moments[..., 0],
            moments[..., 1],
            r,
            phi,
            length * 2.0,
            width * 2.0,
            psi,
            skewness,
            kurtosis,
        ],
        axis=-1,
    ), longitudinal