Monitoring the optical quality of the FACT Cherenkov Telescope

Master’s Thesis

Laurits Tani

Supervisor:

Prof. Dr. Adrian Biland

2019

abs.sectionabs.sectionAbstractAbstractabs.section

## Abstract

In gamma ray astronomy muon events have a distinct feature of casting ring-
like images on the sensor plane, thus forming a well known signal class for
Cherenkov telescopes. These ring-like images can then be used to deduce the
optical point spread function (PSF) which is an important measure of the
optical quality of the imaging-reflector. In this thesis the observed
’fuzziness’ of muon rings is used as a measure to infer the PSF. However to
have a good estimate for this ’fuzziness’ parameter, the reconstruction of the
ring center and ring radius itself needs to be accurate, so different methods
of ring feature extraction are studied. To check for the accuracy of the
methods a simulation and analysis is performed. Measuring the evolution of the
PSF over time allows to identify its effects and take them into account for
the reconstruction of gamma-rays postliminary. As a further benefit of the
methods presented here no additional observations are needed to measure the
PSF nor any human activity on site is required. The accuracy of the method,
and the PSF of FACT vs. time are presented.

###### Contents

  1. 1 Introduciton
    1. 1.1 Astroparticle Physics
      1. 1.1.1 Astrophysical messengers
    2. 1.2 Gamma-ray astronomy
    3. 1.3 FACT Cherenkov Telescope
      1. 1.3.1 Photon-stream representation
    4. 1.4 Creation of muon rings
      1. 1.4.1 Cherenkov effect
      2. 1.4.2 Muons and their properties
      3. 1.4.3 Muon rings
  2. 2 Simulation
    1. 2.1 Individual muon event
    2. 2.2 Multiple muon events
      1. 2.2.1 Choice of parameters
      2. 2.2.2 Simulation process
  3. 3 Methods of investigation
    1. 3.1 Cherenkov photon detection
      1. 3.1.1 Spatial clustering
      2. 3.1.2 DBSCAN (DB)
      3. 3.1.3 DBSCAN’s effectiveness in detecting Cherenkov photons
      4. 3.1.4 DBSCAN’s effectiveness with different night-sky-background rate
    2. 3.2 Extraction of muon ring properties
    3. 3.3 Circle Model (RM)
    4. 3.4 Hough transform(HT)
      1. 3.4.1 Hough transform for circle
    5. 3.5 Used combinations of detection and extraction and their properties
      1. 3.5.1 DB-RM
      2. 3.5.2 DB-RM-HT
    6. 3.6 Conclusions of the ring reconstruction methods
    7. 3.7 Classification of the event
      1. 3.7.1 Acceptance vs opening angle
      2. 3.7.2 Acceptance vs point spread function
    8. 3.8 PSF vs fuzziness
      1. 3.8.1 Standard deviation as fuzziness parameter
      2. 3.8.2 Hough response as fuzziness parameter
    9. 3.9 Comparison of fuzz parameters
  4. 4 Real observations
    1. 4.1 Fuzziness parameters and PSF vs time
      1. 4.1.1 Standard deviation as the fuzziness parameter
      2. 4.1.2 Hough response as the fuzziness parameter
  5. 5 Conclusion and Outlook
    1. 5.1 Conclusion
    2. 5.2 Outlook
  6. A Other muon ring feature extraction methods
    1. A.1 DB-RM-SD
    2. A.2 DB-RM-KC
  7. B Found distributions opening angle

##  Chapter 1 Introduciton

###  1.1 Astroparticle Physics

Astroparticle physics investigates the origins of particles of astronomical
origin. Insight into the structure and history of our universe is gained by
researching the mechanisms responsible for e.g. accelerating particles to high
energies, origins of charged cosmic rays and the nature of dark matter.

####  1.1.1 Astrophysical messengers

Particles coming from all around the Universe that are constantly bombarding
Earth are used as messengers to explore the cosmos. These particles can be
divided into three classes: (electrically charged) cosmic rays, neutrinos and
photons.

Figure 1.1.1: Hadronic shower [1]

##### (Charged) Cosmic Rays

Cosmic rays consist mainly ($\sim$ 90$\%$) out of protons. Upon impact with
the Earth’s atmosphere, they can produce showers of secondary particles that
in case of high enough energy will also reach the ground. This class of
airshowers is called hadronic airshower (see figure 1.1.1). Despite being the
easiest to detect with Cherenkov teleskopes, they do not carry information
about their origins. The energy spectrum of cosmic rays is very featureless,
though there are 3 features: two knees and ankle (respectively at energies
$\sim 10^{15}$ eV, $\sim 10^{17}$ eV and $\sim 10^{19}$ eV). [2]

However cosmic rays do have a theoretical upper limit for energy. It is called
Greizen-Zatsepin-Kuzmin limit (or GZK-cutoff for short). The limit comes from
the fact that universe is not transparent anymore for most energetic cosmic
rays ($\sim 10^{21}$ eV) because they interact with the very low-energetic
cosmic microwave background (CMB) photons.

$p+\gamma_{CMB}\rightarrow\Delta^{+}\rightarrow n+\pi^{+}$

or

$p+\gamma_{CMB}\rightarrow\Delta^{+}\rightarrow p+\pi^{0}$

Figure 1.1.2: Charged cosmic-ray energyspectrum [2]

##### Neutrinos

Neutrinos are uncharged leptons with a very small cross section, thus being
highly uninteractive. In addition to gravity which we can neglect due to
almost zero mass of the neutrino, neutrinos interact with matter only via weak
interactions - charged current and neutral current. Neutrinos are often times
called the ultimate messengers of the universe, because their trajectories are
not bent by cosmic magnetic fields, their signal will not get attenuated due
to their very uninteractive character, thereby being highly penetrating thus
allowing us to see into the core of their origin. Unfortunately being that
uninteractive also poses a problem if one were to detect them. For neutrinos
to interact with detector material, large detectors are needed. Otherwise they
would only fly through the detector without interacting. In interactions they
produce leptons ($e^{\pm},\mu^{\pm},\tau^{\pm})$, which then produce Cherenkov
photons that are easy to detect. As examples for neutrino observatories, there
are IceCube and Super Kamiokande. In those experiments the main background is
the atmospheric neutrinos.

##### Photons

Photons are uncharged elementary particles that have no rest mass. They play a
major role in astroparticle physics due to the fact that their trajectories
are not bent by the intergalactic/cosmic magnetic fields, thus they carry the
information of the place of their origin. Photons that are of interest to us
are highly energetic and are called $\gamma$-rays. The detected photon
energies cover energy range of more than twenty orders of magnitude from $\sim
1$ $\mu$eV (radio waves) up to several hundred TeV ($\gamma$-rays).

###  1.2 Gamma-ray astronomy

Gamma-ray astronomy is the observation of the previously mentioned gamma rays,
the most energetic electromagnetic radiation that has energies above 100 keV.
$\gamma$-rays originating from the solar flares have energies mostly in the
MeV range, although they can be created also in the GeV range. [3] There are
several mechanisms to produce these highest energetic photons - inverse
Compton effect, annihilation of electron and positron, decay of radioactive
material (also known as gamma decay) or the decay of high energetic neutral
pion. The highest photon energies recorded so far are in the TeV range.

The detection of $\gamma$-rays is diffucult, since they are much more rare
than the low-energetic X-rays. For gamma-ray energies higher than 100 GeV the
flux is so rare, that the effective area of space-bound experiments is too
small to have any statistical significance. To overcome this, ground-based
telescopes are used, that use the atmosphere as calorimeter. In contrast to
charged cosmic-rays, gamma-rays produce only electromagnetic airshowers, that
are then detected by ground-based telescopes. The cascade arises from pair-
production and Bremsstrahlung of the electrons and positrons. Electrons and
positrons that move faster than light in the atmosphere, will also emit
Cherenkov radiation.

Although the Imaging Atmospheric Cherenkov Teleskope techniques have currently
the highest sensitivity, $\gamma$-ray astronomy is still limited by the fact,
that at low energies non-gamma-ray background dominates and at higher energies
the $\gamma$-ray flux is too low.

###  1.3 FACT Cherenkov Telescope

The first G-APD Cherenkov Telescope (FACT) is a ground based Imaging
Atmospheric Cherenkov Telescope (IACT) that started operating in October 2011.
It is the first IACT to test novel Geiger-mode avalance photo diodes (G-APDs)
as photosensors. The focal length of the telescope is 4.889 m and field of
view is 4.5°, thus for each individual pixel the field of view is $\sim
0.11\degree$ [4]

The telescope is located next to two 17m MAGIC telescopes at the Observatorio
del Roque de los muchachos on the Canary Island La Palma in Spain at the
altitude $\sim 2200$ m. It constantly monitors several bright sources, among
them for example Mrk421, Mrk501 and the Crab nebula.

Imaging atmospheric-Cherenkov method is an indirect measurement technique,
where Cherenkov photons emitted by charged particles are used for imaging.
There are 1440 photosensors, each equipped with 3600 cells, with individual
readouts in the camera to detect these photons. Using silicon photo
multipliers (SiPM’s which is a different name for G-APD) allows observations
also under bright light conditions.

The reflector is comprised of 30 hexogonal mirrors, with the total reflective
surface of 9.51 m2. These 30 mirrors are arranged in the Davies-Cotton
arrangement and placed at their focal distance to their focal point and
oriented in the direction of a point at twice the focal length. [4] [5]

Figure 1.3.1: Photo of the FACT telescope at the Observatorio del Roque de los
muchachos on La Palma (Canary Islands, Spain) [4]

####  1.3.1 Photon-stream representation

We describe incoming photons as a list of arrival times for each pixel. This
list is called ”photon-stream”. Due to low electronics-noise and fast read-out
of the SiPM sensors, FACT is able to provide the accuracy to take full
advantage of the ”photon-stream”. [6] [1]

On figure 1.3.2 one can see the depiction of an event in that representation.

Figure 1.3.2: Recorded event in photon-stream representation. Each blue dot
corresponds to a single photon.

###  1.4 Creation of muon rings

####  1.4.1 Cherenkov effect

Charged particles travelling in a medium with velocities larger than the local
speed of light emit Cherenkov photons. These photons are emitted in spherical
waves, thus forming a cone of photons with an opening angle that is dependent
on the refractive index of the medium and the velocity of the charged
particle.

As can be seen on right side figure of figure 1.4.1, the opening angle of the
cone is:

$\theta_{C}=arccos\left(\frac{1}{n\cdot\beta}\right)$ (1.1)

The refractive index of air at altitude 2200 m is
$n=\frac{c_{0}}{c}=1+\epsilon=1.00022$ where $c_{0}$ is the speed of light in
vacuum, c is the speed of light in medium, and $\epsilon$ is the fractional
refractive index.

In order to produce Cherenkov emission, the following relation must hold:

$\displaystyle\begin{aligned} v&>c\\\ \implies\beta&>\frac{1}{n}\\\
\implies\beta_{min}&=0.99977\end{aligned}$ (1.2)

Figure 1.4.1: Cherenkov effect. On the left: charged particle travelling in
medium with speed less than lightspeed in that medium ($v<c$). On the right:
charged particle travelling in medium with speed larger than lightspeed in
that medium ($v>c$). Figure taken from [7]

####  1.4.2 Muons and their properties

Muons have a mass of $m=105.6583745\pm 0.0000024$ MeV and a mean lifetime in
their reference frame of $\tau=2.1969811\pm 0.0000022\cdot 10^{-6}$ seconds.
Muons are produced in the atmosphere in hadronic airshowers as the decay
product of charged pions:

$\pi^{-}\to\mu^{-}+\overline{\nu}_{\mu}$

and

$\pi^{+}\to\mu^{+}+\nu_{\mu}$

Muons produced in hadronic showers are usually relativistic particles, thus
many reach ground before decaying due to time dilation. The decay products of
muons are:

$\mu^{-}\rightarrow e^{-}+\overline{\nu_{e}}+\nu_{\mu}$ (1.3)

or

$\mu^{+}+\rightarrow e^{+}+\overline{\nu_{\mu}}+\nu_{e}$ (1.4)

Since the Lorentz factor is defined as

$\gamma=\frac{1}{\sqrt{1-\beta^{2}}}=\frac{E}{E_{0}}>1$ (1.5)

then we find that the minimal energy for muons to produce Cherenkov radiation
is:

$E_{min}=\frac{E_{0}}{\sqrt{1-\beta^{2}_{min}}}\approx 4.93GeV$ (1.6)

Since with higher energy the opening angle increases then we can find out the
maximal opening angle corresponding to the maximum speed of the particle
($\beta=1$). Using equation 1.1 we find:

$\theta_{max}=\frac{1}{n}\approx 1.22\degree$ (1.7)

####  1.4.3 Muon rings

Photons that are hitting the aperture of the telescope are reflected to the
sensor. Since the aperture is an imaging mirror, then all photons with the
same direction-vector are reflected onto the same place in the sensor. Thus
all Cherenkov photons that are emitted in a cone with a certain fixed angle
will form a ring-like shape in the sensor plane. Unfortunately muon events are
very dim, so they are difficult to detect. For FACT they are among the dimmest
events that can be triggered.

##### Zero inclination angle and distance within aperture radius

In this case muon rings that are forming on the sensor plane will be evenly
illuminated and always form in the center of the camera. On figure 1.4.2 one
can see how muon rings formation in case of zero inclination angle and with an
impact radius of zero.

Figure 1.4.2: Formation of the muon ring in case of $\alpha=0$ and $r=0$.
Cherenkov cone opening angle $\theta_{C}=\theta_{max}$ Figure taken from [8]

##### Non-zero inclination angle

On figure 1.4.3 one can see the case when the inclination angle of the muon is
not zero. Depending on the inclination angle of the muon, the center position
of the ring moves on the sensorplane.

Figure 1.4.3: $\alpha=2.5^{o}$, $r=-0.25$m, $\theta_{C}=\theta_{max}$. Figure
taken from [8]

##### Greater distance than the aperture radius

If muons are hitting the aperture plane further away than the radius of the
aperture, then the Cherenkov photons do not produce full rings anymore - some
amount of photons, depending on the distance, won’t hit the mirror anymore.
Consequently the muon ring is not illuminated evenly and the ring might be
incomplete. On figure 1.4.4 one can see the case when muon is not hitting the
aperture of the telescope anymore.

Figure 1.4.4: $\alpha=0$, $r=1.2$ m, $\theta_{C}=\theta_{max}$ Figure taken
from [8]

##  Chapter 2 Simulation

In order to benchmark the effectiveness of muon detection and ring feature
extraction methods, we simulate FACT muon event observations. For this purpose
a custom muon simulation was created. Having a simulation allows us to relate
observed ’fuzziness’ (the spread of photons around ring line) of the muon ring
with the optical point spread function. (For scripts see [9])

###  2.1 Individual muon event

A single muon will be simulated, given the following parameters:

  * •

direction vector of the muon

  * •

support position of the muon

  * •

night-sky-background rate

  * •

opening angle of the Cherenkov cone

  * •

Cherenkov photon emission rate

  * •

point spread function of the reflector

  * •

arrival time standard deviation of Cherenkov photons

A muon will emit Cherenkov photons when traversing to the ground. Distance
travelled before emitting another Cherenkov photon is Poisson distributed. In
case the emitted Cherenkov photons hit the aperture plane within the aperture
radius, they are reflected onto the sensor plane. Also additional night-sky-
background photons are simulated. On figure 2.1.1 is a simulated muon event
with non-zero inclination angle and impact parameter that is larger than
aperture radius.

Figure 2.1.1: Example of a muon event with $\alpha\neq 0$ and
$r_{impact}>r_{aperture}$ and night-sky-background rate of $35\times
10^{6}s^{-1}$ which corresponds to the dark night

###  2.2 Multiple muon events

####  2.2.1 Choice of parameters

In order to have all possible muon events within a given range muon directions
on the skydome and its support position are drawn randomly from uniform
distribution. In addition also the opening angle of the Cherenkov cone is
drawn from an uniform distribution. These parameters are then given to the
function that simulates all these individual muon events.

####  2.2.2 Simulation process

There are 13 possible input parameters:

  * •

Number of muons to be simulated

  * •

Maximum inclination angle

  * •

Maximum aperture radius of muon support vector

  * •

Minimum opening angle of the Cherenkov cone [default: 0.4 degrees]

  * •

Maximum opening angle of the Cherenkov cone [default: 1.6 degrees]

  * •

Night-sky-background-photon rate per pixel [default $35\times
10^{6}s^{-1}pixel^{-1}$ which corresponds to dark nights] [5]

  * •

Standard deviation of the arrival times of the photons [default: $0.5\times
10^{-9}$ seconds]

  * •

Rate of Cherenkov photons to be generated per meter [default: 3.0] (This value
was chosen based on observed muon events. Otherwise Frank-Tamm formula,
refractive index, mirror reflectivity, photo-sensor-efficiency and possibly
even more should be taken into account to get this value.) [10]

  * •

FACT aperture radius [default: 1.965 m]

  * •

Standard deviation of the Point Spread Function (PSF) [default: 0] (for
perfect imaging)

##  Chapter 3 Methods of investigation

###  3.1 Cherenkov photon detection

Due to the fact that muon events are very dim, the biggest obstacle in seeing
these is the triggering the telescope. In case of high night-sky-background
rate the trigger threshold is higher and the dim muon event will not trigger
the telescope.

Thus the next step is the identification of Cherenkov photons (detection for
short) from night-sky-background photons. Because Cherenkov photons reside
much more densly than night-sky-background photons, we need to use a density
based algorithm like Two-Level-Time-Neighbor-Cleaning or density clustering.

####  3.1.1 Spatial clustering

Spatial clustering is a technique used in data mining, which groups objects
into classes or clusters based on their similarities within the cluster. One
example of these spatial clustering techniques is DBSCAN ([11]), that will be
used in this thesis.

####  3.1.2 DBSCAN (DB)

DBSCAN (Density-Based Spatial Clustering of Applications with Noise)
recognizes Cherenkov photon clusters by taking into account the fact that
within this cluster the typical density of photons is much higher than outside
of the cluster. Since it can discover arbitrary shaped Cherenkov clusters, it
is very powerful. Basic ideas:

  * •

The neighbourhood (on the $c_{x}$ and $c_{y}$ (the directional cosines) plane)
that is within a range $\epsilon$ of a given photon is called the
$\epsilon$-neighborhood

  * •

For a photon to be considered to be the core point of a cluster, there has to
be a minimum number of photons, MinPts, inside its $\epsilon$-neighborhood

  * •

A photon is directly density reachable if it is within $\epsilon$-neighborhood
of a core photon

  * •

For a photon to be density reachable from photon P to photon Q, there has to
be a chain of photons beween that are directly density reachable

  * •

For a photon P to be density-connected to photon Q both have to to be density
reachable with respect to a photon O

  * •

All photons that are not inside a density based cluster are considered to be
noise.

A new cluster is created if a photon P has within its $\epsilon$-neighborhood
more than MinPts photons. During iteration some clusters may be connected and
lasts until no more photon can be added to any clusters. [12]

####  3.1.3 DBSCAN’s effectiveness in detecting Cherenkov photons

A benchmark for the sensitivity and precision of this method was set by
simulating muon events with a set night-sky-background (corresponding to dark
night in this case).

Sensitivity ([13]) is calculated as:

$\frac{true\\_positives}{true\\_positives+false\\_negatives}$ (3.1)

and precision ([13]) is calculated as:

$\frac{true\\_positives}{true\\_positives+false\\_positives}$ (3.2)

where:

  * •

true_positives is the number of photons that were classified correctly to be
Cherenkov photons,

  * •

false_positives is the number of photons that were night-sky-background
photons that were classified to be Cherenkov photons,

  * •

true_negatives is the number of photons that were classified correctly to be
night-sky-ackground photons,

  * •

false_negatives is the number of photons that were Cherenkov photons that were
classified to be night-sky-background photons.

The used parameters for DBSCAN were:

Parameter | Value
---|---
min-samples | 20
$\epsilon$-radius | 0.45 degrees

With these parameters the found sensitivity and precision corresponding to
dark night night-sky-background were $98.9\pm 0.3$ and $82.7\pm 0.4$
respectively.

Naturally it is expected that the precision of detecting muons can not be 100
percent due to the fact that night-sky-background photons also will reside in
close proximity to the cluster center and thus be considered to part of the
cluster (for 2D projection of an event see figure 3.1.1). This happens even
with the night-sky-background rate corresponding to dark night. Furthermore,
detecting Cherenkov photons takes only into account the spatial positioning
and time of an individual photon, thus in order to increase precision one
needs to apply some other cuts after the density clustering.

Figure 3.1.1: Inside Cherenkov cluster (black dots) reside also some night-
sky-background photons (red dots). Before clustering

####  3.1.4 DBSCAN’s effectiveness with different night-sky-background rate

Since observations are done under different night-sky-bacground-light-
conditions 10 different night-sky-background rates were simulated. Although
adjusting DBSCANs parameters might improve the precision, it was not
implemented. Regardless one can see from figures 3.1.2 and 3.1.3 that only
precision drops with bigger night-sky-background rate while sensitivity stays
rouhly the same. This was of course expected, since increasing the night-sky-
background rate only increases the chances of a photon residing in the
vicinity of the actual Cherenkov photon cluster and thus being included as one
of the elements. At rougly 2 times the dark night rate (dark night night-sky-
background rate taken to be $35\times 10^{6}$ photons per second per pixel)
the precision already drops to roughly 60%.

Figure 3.1.2: Precision of the DBSCAN algorithm with different night-sky-
background rates Figure 3.1.3: Sensitivity of the DBSCAN algorithm with
different night-sky-background rates

Due to inhomogeneities in the night-sky-background rate and cosmic rays,
branches as depicted on figure 3.1.4 can emerge due to the fact that only
densities are taken into account when searching for Cherenkov photons.

Figure 3.1.4: Having much night-sky-background might cause branches consiting
out of night-sky-background photons (red cluster on figure)

Taking into account muon ring geometry, we can implement a cut, discarding
photons that reside from the reconstructed ringline more than 2 standard
deviations away (See figure 3.1.5), thus ending up with a muon ring similar to
the one on the figure 3.1.6.

Figure 3.1.5: Implementing cut on the $c_{x}$-$c_{y}$ plane (DBSCAN+) when the
distance of a photon from reconstructed ring radius is larger than 2 standard
deviation Figure 3.1.6: Implementing cut on the $c_{x}$-$c_{y}$ plane
(DBSCAN+) when the distance of a photon from reconstructed ring radius is
larger than 2 standard deviation

This cut (DBSCAN+) notably increased the precision of the detection while
almost not decreasing sensitivity at all as one can see on figures 3.1.7 and
3.1.8. With DBSCAN+ at 2 times the dark night night-sky-background rate the
precision has risen to $\sim 70\%$.

Figure 3.1.7: Precisions of the DBSCAN algorithm with and DBSCAN+ different
night-sky-background rates Figure 3.1.8: Sensitivities of the DBSCAN algorithm
and DBSCAN+ with different night-sky-background rates

Although for night-sky-background rate of $70\times 10^{6}$ photons
$s^{-1}pixel^{-1}$ the precision is $\sim 70\%$, we can see from 3.1.9, then
most of the observations are done with less night-sky-background, thus having
higher precision.

Figure 3.1.9: Observation time for different night-sky-background rates (dark
night = $35\times 10^{6}$ photons $s^{-1}pixel^{-1}$) [1]

###  3.2 Extraction of muon ring properties

The goal of extracting the muon ring features (extraction for short) is to get
the radius and the center position of the muon ring. Cherenkov photons are
expected to form a sharp ring on the sensor plane. Due to camera’s pixellation
and misaligned reflector facets the ring image is blurred (see figure 3.2.1).
To describe the spread, a fuzziness parameter is used.

Figure 3.2.1: Sharp ring (dashed grey line) is blurred (solid black line) due
to pixellation of the camera and misaligned reflector facets

Accurate reconstruction of the ring parameters is highly preferable to avoid
large uncertainties in the fuzziness parameter. In contrast to detection, the
choice for extraction methods is much wider Methods of extraction include for
example a simple Circle Model and Hough Transformation.

###  3.3 Circle Model (RM)

The Circle model used in this thesis is taken directly from skimage.measure
package ([14]). It estimates 2D circles using total least squares method. The
functional model for the circle is

$r^{2}=(x-x_{c})^{2}+(y-y_{c})^{2}$ (3.3)

This model minimizes the squared distances from all the points to the circle:

$min\\{\sum(r-\sqrt[]{(x-x_{c})^{2}+(y-y_{c})^{2}})\\}$ (3.4)

For this model a minimum of 3 photons is required. (See examples in [15])

###  3.4 Hough transform(HT)

Hough transform ([16]) is a technique used in computer vision to identify
objects from a background. This technique finds imperfect instances of objects
that lie within a certain class. The best candidate is chosen via a
integrating procedure that is carried out in a parameter space (called Hough
accumulator or Hough space).

Dimension of parameter space depends on the amount of variables needed to
describe a certain object. For example a line can be described by 2
parameters: slope and distance from the origin, thus the parameter space is
two dimensional. Visually it would be represented like a 2D-plane (see figure
3.4.1)

Figure 3.4.1: Example of 2D parameter space (Hough accumulator) for a line.
Darkest bin corresponds to best parameters for a line

Collecting all votes from the parameter space we obtain the local maxima,
whose coordinates in this accumulator correspond to the best fit of the
sought-after object (be it a line, circle or something else).

####  3.4.1 Hough transform for circle

Since in this thesis the goal is to infer true optical pointspread function
using muon rings, then in order to reconstruct the rings we need Hough
transform for recognizing rings. As a first attempt the accumulator was filled
with votes that come from astropy Ring2D fit ([17]. Unfortunately to achieve
good enough accuracy the grid of the accumulator needed to be very fine.
Consequently this increased also the computational time needed.

To overcome this problem a Hough transform script was written, that had some
extra features: instead of the Heaviside weighing function ([18]) (see figure
3.4.2) triangular weighing function was used (see figure 3.4.3), thus avoiding
same value bin-plateaus in the Hough accumulator, thereby making it possible
to choose one best parametrization for the ring. This increases the accuracy
when deciding which of the accumulator bins had the most votes.

Furthermore our custom Hough transform is iterative, so in every interation-
step the area to be scanned gets smaller but with the same bin-count. Thus
computing time is reduced and higher accuracy is reached.

Finally, an initial guess can be given to the transformation together with the
total uncertainty, thus decreasing the amount of parameter space to be scanned
by the Hough transformation. In case when no initial guess is available, one
should adjust the uncertainty according to the amount of area to be scanned.
(See more details in [19])

Figure 3.4.2: Heaviside weighing function Figure 3.4.3: Triangular weighing
function

###  3.5 Used combinations of detection and extraction and their properties

In the following one can see the results of 300 000 muons that were simulated
and different detection and extraction method combinations were tested.

####  3.5.1 DB-RM

The most straight-forward method to test was density clustering with a simple
ring model. As the name suggests then in order to detect Cherenkov photons
from the night-sky-background photons, this method uses density clustering
(DBSCAN algorithm to be more exact). Found Cherenkov photons are then used to
extract muon rings and their features. The first estimate is done by using a
regular CircleModel from skimage which is wrapped with RanSaC (Random Sampling
Consensus) ([20]) to get rid of the outliers.

As one can see from figure 3.5.2 and figure 3.5.3 at first glance it would
seem that the center positions of the muon ring are reconstructed fairly well.
Unfortunately this cannot be said about the ring radius. When one looks closer
at the scale, then it becomes clear that ring center positions are
reconstructed with an error of rougly three quarters of a degree. This error
causes also the ring radius to be reconstructed very inaccurately and one can
even see a major bias - this method underestimates the ring radius very often.

When subtracting true opening angle from reconstructed one, one sees how
inaccurately the radius is reconstructed (see figure 3.5.1). The standard
deviation of the distribution for the difference of reconstructed opening
angle and true opening angle is 0.183 degrees.

On figure 3.5.5 we can see that indeed the simulated muons support position on
apperture plane is distributed uniformly over the disk with the given aperture
radius (4 m). As one would expect then the muons which hit the aperture plane
closer to the center were detected more efficiently. On figure 3.5.6 we see,
that the simulated muons had directions uniformly distributed and smaller
inclination angles were more likely to be detected.

Figure 3.5.1: Reconstructed opening angle minus true opening angle Figure
3.5.2: Confusion matrix for cx with DB-RM method Figure 3.5.3: Confusion
matrix for cy with DB-RM method Figure 3.5.4: Confusion matrix for opening
angle with DB-RM method Figure 3.5.5: Impact radius squared with DB-RM method
Figure 3.5.6: Ring center coordinates squared with DB-RM method

####  3.5.2 DB-RM-HT

The custom Hough transform script uses initial guess for improving the
accuracy of the reconstruction. As a byproduct it reduces computing time. The
initial guess for ring parameters is taken from the DB-RM results. To increase
accuracy even more, one performs multiple Hough transforms (iterative Hough,
decreasing bin sizes with every step), taking the results from the previous
iteration to be the initial guess and thus reducing the amount of Hough space
to be scanned each iteration.

As one can see from figures 3.5.8 and 3.5.9 it is clear that the extracted
ring center position is very accurate. Furthermore, when one compares
reconstructed ring radius with Hough transform and the one from Circle Model
(DB-RM) on figure 3.5.7 it is obvious that Hough transform version is
superior, having much narrower and symmetric distribution.

In addition in order to reach this accuracy without iterative Hough, one needs
$\sim 100$ bins for all of the 3 variable, thus $10^{6}$ bins needed to be
checked. With iterative Hough only $\sim 8000$ bins needs to be checked, thus
the computing time reduced by a factor of 125 at least.

Figure 3.5.7: Opening angle of the Cherenkov cone Figure 3.5.8: Confusion
matrix for cx with DB-RM-HT method Figure 3.5.9: Confusion matrix for cy with
DB-RM-HT method Figure 3.5.10: Confusion matrix for opening angle with DB-RM-
HT method Figure 3.5.11: Ring center coordinates squared with DB-RM-HT method

Comparison of DB-RM and DM-RM-HT ring center reconstructions are found on
figures 3.5.12 and 3.5.12. The standard deviations of the distributions are in
table 3.1.

Figure 3.5.12: Reconstructed cx minus true cx Figure 3.5.13: Reconstructed cy minus true cy Method | $\sigma(c_{x})(degrees)$ | $\sigma(c_{y})(degrees)$
---|---|---
DB-RM-HT | 0.0712 | 0.0703
DB-RM | 0.183 | 0.180
Table 3.1: Standard deviation of the difference of true and reconstructed ring
center

###  3.6 Conclusions of the ring reconstruction methods

In table 3.2 are listed the standard deviations of the difference of true and
reconstructed opening angles. It is evident, that DB-RM-HT method is superior
to other tried methods, thereby reconstructing rings $\sim 2.5$ times more
accurately.

Method | $\sigma(degrees)$
---|---
DB-KC-SD 111see appendix A.2 | 0.0202
DB-RM-HT | 0.0730
DB-RM | 0.183
DB-RM-SD 222see appendix A.1 | 0.187
Table 3.2: Standard deviation of the difference of true and reconstructed
opening angle

###  3.7 Classification of the event

For classification of the events the following cuts are applied after
Cherenkov photon detection and ring feature extraction:

Parameter | cut
---|---
number-of-photons-min | 3
muon-ring-radius-min | 0.45 deg
muon-ring-radius-max | 1.6 deg
muon-ring-overlapp-with-field-of-view-min | 20%
arrival-time-stddev-max | 5e-9 s
initial-circle-model-photon-ratio-min | 0.6
visible-ring-circumfance-min | 1.5 deg
off-density | 0
on-off-ratio-min | 3.5
density-circle-model-inner-ratio-max | 0.25
evenly-population | True
Table 3.3: Current cut parameters to find muons

where:

  * •

number-of-photons-min: Minimum number of photons required to perform Circle
Model.

  * •

muon-ring-radius-min: Minimum muon ring radius for it to be classified as muon
ring.

  * •

muon-ring-radius-max: Maximum muon ring radius for it to be classified as muon
ring.

  * •

muon-ring-overlapp-with-field-of-view-min: Minimum muon ring overlap with the
field of view (fraction).

  * •

arrival-time-stddev-max: Maximum arrival time standard deviation.

  * •

initial-circle-model-photon-ratio-min: Minimum fraction of inlier photons
divided by all photons.

  * •

visible-ring-circumfance-min: Minimum visible ring circumference in degrees.

  * •

off-density: average density of photons offset from the ring circumference.

  * •

on-off-ratio-min: Minimal fraction of on-density (points on ring circle)
divided by off-density

  * •

density-circle-model-inner-ratio-max: Maximum fraction of number of photons
off ring inside circle divided by number of all photons.

  * •

evenly-population: Check whether ring is populated evenly. Threshold for
evenly-population to be true is if the standard deviation is less than 0.8

With these cuts, the performance of classifying events to be muons had the
following performance:

Precision | 99.3%
---|---
Sensitivity | 67.2%

####  3.7.1 Acceptance vs opening angle

For the best found method, DB-RM-HT, also muon acceptance dependent on the
opening angle of the Cherenkov cone was found by simulating muons with a set
opening angle. As one can see from figure 3.7.1 then the effective area peaks
near $\sim 1.2$ degrees, which is also the maximum angle under which muons in
at the height of the Telescope can emit Cherenkov emission.

Figure 3.7.1: Acceptance vs opening angle

####  3.7.2 Acceptance vs point spread function

Also muon acceptance dependence was studied for different point spread
functions. As one can see from figure 3.7.2 then the acceptance reduces when
muon ring image gets fuzzier.

Figure 3.7.2: Acceptance vs point spread function

###  3.8 PSF vs fuzziness

Since the goal of the thesis is to measure PSF using muon rings, then we need
to relate muon ring fuzziness parameter to the actual PSF. Since perfect muon
rings should be sharp, then the task is easily achievable by simulating muons
with different point spread functions and relating them to the calculated
fuzziness. In total 300 000 muons were simulated for every different point
spread function between 0 and 0.20 degrees PSF with 20 steps.

We expect the fuzz vs. PSF plot to have a 3 regions that can be approximated
with lines with different slopes. This is expected because in the region with
very small point spread function it will be saturated due to the pixellation
([21]). The theoretical limit for the point spread function caused by
pixellation would be thus
$\sigma_{pixel}=\frac{1}{\sqrt{12}}\times\theta=0.0289$. The second region is
expected to be linear because of the relation between fuzziness PSF should be
linear. Third region starts roughly when the PSF is the size of the pixel and
thus unexpected behaviour is expected. Thus it makes sense to use a third
degree polynomial fit, because this is the lowest order polynomial that has
two curvatures. Because we want to have still some info about saturated parts,
we will keep some of it and fit a polynomial of third degree.

####  3.8.1 Standard deviation as fuzziness parameter

This fuzziness parameter is the standard deviation of the distances of all the
photons from the ringline.

On figure 3.8.1 one can see a clear fuzziness dependence on point spread
function, although there is a slight saturation with very low point spread
function caused by the pixellation of the camera.

Figure 3.8.1: Fuzziness vs point spread function

Fit for standard deviation as fuzziness parameter:

$-53.5x^{3}+9.14x^{2}-1.32\times 10^{-1}x+1.35\times 10^{-1}$ (3.5)

Figure 3.8.2: Standard deviation as fuzziness parameter

####  3.8.2 Hough response as fuzziness parameter

The fuzziness parameter called ‘Hough response’ or just ‘response’ is the
normed value of the highest value bin in the Hough accumulator (Hough space).
It is calculated as:

$\frac{highest\\_value\\_bin}{all\\_photons\\_count}*100$ (3.6)

Thus for the sharpest ring the response would be 100

Fit for response as the fuzziness parameter:

$2.77\times 10^{4}x^{3}-2.76\times 10^{3}x^{2}-2.86\times 10^{2}x+51.3$ (3.7)

Figure 3.8.3: Response as the fuzziness parameter Figure 3.8.4: Response as
the fuzziness parameter

Comparing DB-RM (ringModel) and DB-RM-HT (Hough) on figure 3.8.4 we have an
indication that Hough transform would outperform ringModel by having steeper
slope in the correlation of fuzziness vs. point spread function.

###  3.9 Comparison of fuzz parameters

As one can see from figure 3.9.1 an important advantage of the fuzziness
parameter ’response’ over ’stdev’ (see subsection 3.8.1) is that the relative
standard deviation is smaller manyfold, thus one would expect much more
accurate results when reconstructiong PSF, which we will see in the next
chapter.

Figure 3.9.1: Relative standard error for fuzziness parameter value vs psf

##  Chapter 4 Real observations

###  4.1 Fuzziness parameters and PSF vs time

As we saw in chapter 3 that Hough transform (DB-RM-HT) outperforms ringM (DB-
RM) in reconstructing ring parameters. In the following only results using
Hough transform as the extraction method are shown. To arrive at our goal to
see how point spread function varies in time we have to go through four steps:

  * •

Simulate event with different point spread functions to get the relation
between our chosen fuzziness parameter and true PSF (See section 3.8)

  * •

Fit a curve to infer PSF (see also section 3.8)

  * •

Find fuzziness parameter over time

  * •

Convert fuzziness to PSF

####  4.1.1 Standard deviation as the fuzziness parameter

As the first step we calculate the fuzziness parameter (currently standard
deviation of the photon from the reconstructed ring) for all the observations.
This gives us figure 4.1.1. The dimmer a datapoint is, the less significant
the date is, because the less muons were detected. The size of the errorbar
around the datapoint is indicates how big was the standard error of the
fuzziness parameter at that date. Already from this plot we can see rougly how
true point spread function would behave. The green lines on the plot show
timestamps where mirror alignments were done. Furthermore there is a hatched
area on the figure, which shows when the data taken was not reliable due to
bad DRS (Domino Ring Sampler) calibration.

Figure 4.1.1: Standard deviation as the fuzziness parameter over time

####  4.1.2 Hough response as the fuzziness parameter

As seen on the figure 4.1.2, it is clear that this fuzziness parameter
distinguishes different point spread functions very well, thus changes in PSF
should be visible. One should not forget, that the higher the response, the
smaller the point spread function. Furthermore, as was shown on figure 3.9.1,
then we would expect not much noise in the PSF vs. time figure. Those
assumptions turn out to be true, as one can see from figure 4.1.3.

Green dashed lines on figures 4.1.2, 4.1.3 indicate times when a mirror
alignment was done. We can see that after the most thorough mirror realignment
at 14th of May 2015 (fourth green dashed line from the left) the point spread
function dropped significantly as was also expected from the recorded pictures
(see figures 4.1.4 and 4.1.5). In the beginning of the year 2016 the PSF
degraded very suddenly and it was most probably caused by the high wind speeds
and snow. The dimmer the datapoint, the less important the night was due to
less muons detected.

Figure 4.1.2: Response as the fuzziness parameter over time Figure 4.1.3:
Reconstructed PSF using response as the fuzziness parameter vs time

As one might notice when comparing figures 4.1.3 and 4.1.5, the measured point
spread function is much larger than the reconstructed one after the mirror
alignment on 14th of May 2014. The difference probably comes from the fact,
that if even one mirror is misaligned a lot then the image cleaning will not
find the photons reflected by this particular facet because it will not be in
a dense enough cluster. In contrast, the equipment measuring point spread
function is able to spot also the photons reflected by the misaligned facet.
This is however not a big shortcoming of the reconstruction of the image
cleaning because it takes into account only the differences, not the overall
value, thus it is possible to scale it accordingly if PSF at one point is
known.

Figure 4.1.4: PSF before 10.05.2015 alignment [22] Figure 4.1.5: PSF after
10.05.2015 alignment [22]

##  Chapter 5 Conclusion and Outlook

###  5.1 Conclusion

Finding muon rings and using them to infer the true optical point spread
function turned out to be a good method that can be used efficiently also on
FACT.

When using photon stream, using density clustering (in our case DBSCAN) turned
out to be very efficient and precise method for identifying Cherenkov photons.
We managed to even increase the precision by making cuts to get rid of the
branch-like structures which emerged due to inhomogeneities in the night-sky-
background.

After finding the Cherenkov photons, an accurate way to reconstruct the muon
ring was needed. The most accurate way was found out to be DB-RM-HT, which
made use of scikit-learn’s CircleModel by taking it’s results as a rough
estimate for the ring parameters. These ring parameters were used as an
initial guess for Hough transform. In contrast with AstroPy’s Circle2D model,
our Hough transform used triangular weighing function instead of the Heaviside
function.

Third obstacle to overcome was finding the best parameter for the fuzziness of
the muon ring, which then could be used as a tool to infer the true optical
point spread function. The fuzziness parameter ’response’ turned out to be the
best due to it’s smaller uncertainty.

As a final step the fuzziness parameter was to be calculated using real
observations. After getting the fuzziness over time it was possible to infer
the point spread function fit that we got by fitting the datapoints that we
got from simulating different night-sky-backgrounds.

It was clearly visible from the final PSF-vs-time plot how the reconstructed
PSF correlated very well with the mirror alignments. Due to stormy weather in
the beginning of 2016 some of the mirror facets were misaligned thus worsening
the point spread function, which is clearly visible from the presented
figures. Using the results it is clear that mirrors should be again aligned.

All these found methods turned out to be very good when one wants to infer the
PSF from the past, thus allowing one to reconstruct gamma-ray flux with better
instrument response functions.

###  5.2 Outlook

Although found methods were good in achieving the goals, there is still room
to improve the given methods, that were not implemented in this thesis.

Finding Cherenkov photons from the night-sky-background photons might benefit
when one would adjust DBSCAN’s $\epsilon$-neighborhood depending on the night-
sky-background rate. Also the cut made after DBSCAN can be adjusted depending
on the night-sky-background rate.

After finding the Cherenkov photons and reconstructing the rings one might
want to find better cut values for classifying events to be muon-like or not,
thus increasing the sensitivity and precision of the detection method. For the
evaluation of the detection method a bigger manually classified set of events
would be beneficial. This can be achieved by using Zooniverse as did VERITAS
(see [23]).

## Acknowledgements

I wish to thank my supervisor prof. Dr. Adrian Biland, who read through my
drafts multiple times and gave very good advice on how to continue my thesis.
Furthermore Sebastian A. Müller and Dominik Neise were always there to give
advice and help out in a difficult situation. Last but not least Axel Arbet-
Engels was very helpful giving comments and opinions when I needed them.

## References

  * [1] S. A. Mueller, “Cherenkov-plenoscope,” arXiv preprint arXiv:1904.13368, 2019\.
  * [2] K. A. Olive, P. D. Group, et al., “Review of particle physics,” Chinese Physics C, vol. 38, no. 9, p. 090001, 2014.
  * [3] “Solar Flares: Gamma Rays.” http://cds.cern.ch/record/367281/files/9810089.pdf.
  * [4] H. Anderhub, M. Backes, A. Biland, V. Boccone, I. Braun, T. Bretz, J. Buß, F. Cadoux, V. Commichau, L. Djambazov, et al., “Design and operation of fact–the first g-apd cherenkov telescope,” Journal of Instrumentation, vol. 8, no. 06, p. P06008, 2013.
  * [5] A. Biland, T. Bretz, J. Buß, V. Commichau, L. Djambazov, D. Dorner, S. Einecke, D. Eisenacher, J. Freiwald, O. Grimm, et al., “Calibration and performance of the photon sensor response of fact—the first g-apd cherenkov telescope,” Journal of Instrumentation, vol. 9, no. 10, p. P10012, 2014.
  * [6] “Photon Stream.” https://github.com/fact-project/photon_stream.
  * [7] “Cherenkov Radiation.” https://icecube.wisc.edu/~kjero/Bootcamp/2015/Notebooks/Reconstruction_Introduction.html.
  * [8] M. Nöthe, “Extraktion und analyse von myon-ereignissen am fact-teleskop,” Bachelorarbeit. TU Dortmund, 2014.
  * [9] “Muon Ring Simulation.” https://github.com/fact-project/muons/tree/master/muons/muon_ring_simulation.
  * [10] I. Frank and I. Tamm, “Coherent visible radiation of fast electrons passing through matter,” in Selected Papers, pp. 29–35, Springer, 1991.
  * [11] M. Ester, H.-P. Kriegel, J. Sander, X. Xu, et al., “A density-based algorithm for discovering clusters in large spatial databases with noise.,” in Kdd, vol. 96, pp. 226–231, 1996.
  * [12] B. Borah and D. Bhattacharyya, “An improved sampling-based dbscan for large spatial databases,” in Intelligent Sensing and Information Processing, 2004\. Proceedings of International Conference on, pp. 92–96, IEEE, 2004.
  * [13] D. M. Powers, “Evaluation: from precision, recall and f-measure to roc, informedness, markedness and correlation,” 2011.
  * [14] “skimage.measure.” https://scikit-image.org/docs/0.11.x/api/skimage.measure.html.
  * [15] “CircleModel.” http://scikit-image.org/docs/dev/api/skimage.measure.html#skimage.measure.CircleModel.
  * [16] R. O. Duda and P. E. Hart, “Use of the hough transformation to detect lines and curves in pictures,” Communications of the ACM, vol. 15, no. 1, pp. 11–15, 1972.
  * [17] “Ring2D.” http://docs.astropy.org/en/stable/api/astropy.modeling.functional_models.Ring2D.html.
  * [18] “Heaviside weighing function.” http://mathworld.wolfram.com/HeavisideStepFunction.html.
  * [19] “Circlehough.” https://github.com/Laurits7/circlehough.
  * [20] M. A. Fischler and R. C. Bolles, “Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography,” Communications of the ACM, vol. 24, no. 6, pp. 381–395, 1981\.
  * [21] B. Kamgar-Parsi and W. Sander, “Quantization error in spatial sampling: comparison between square and hexagonal pixels,” in Computer Vision and Pattern Recognition, 1989. Proceedings CVPR’89., IEEE Computer Society Conference on, pp. 604–611, IEEE, 1989.
  * [22] S. Müller, “Clear sight in cherenkov-astronomy. investigations of reflector geometry and mirror alignment for the fact-telescope,” Magisterarb. TU-Dortmund, 2014.
  * [23] “Muon Hunter.” https://www.zooniverse.org/projects/zooniverse/muon-hunter/.

##  Appendix A Other muon ring feature extraction methods

###  A.1 DB-RM-SD

In hopes that ring center position is reconstructed accurately enough, method
DB-RM-SD was tested.

Once again in order to identify Cherenkov photons, density clustering is used.
Now instead of directly using the results from RingModel, we use the found
ring-center positions and estimate ring radius once again by calculating the
distance of every single Cherenkov photon from the previously found center and
taking the median of the distances to be the new radius.

As one can see figure A.1.1, then the reconstructed radius still had very
similar characteristic - it once again underestimated the ring radius very
often.

When subtracting true opening angle from reconstructed one one sees how bad
the radius is reconstructed (see figure A.1.2). The standard deviation of the
difference of reconstructed opening angle and true opening angle for this
method is even worse than before 0.187 degrees. So, this method doesn’t imrove
the results from DB-RM.

Figure A.1.1: Confusion matrix for opening angle with DB-RM-SD method Figure
A.1.2: Reconstructed opening angle minus true opening angle

###  A.2 DB-RM-KC

Inaccuracy in finding correct ring center motivated to test the case when ring
center was as accurate as possible - it was taken from the simulation truth.
Although one cannot know the true ring center in reality, this check showed
that guessing accurately the ring center has a huge impact in reconstructing
the ring radius. As one can see in figure A.2.1 then in that case the muon
ring radius was also reconstructed very accurately as one would have expected.
Of course this was only a test how much the choice of the center position
influences the reconstrucion of the ring radius.

As one can see from figure A.2.2 the ring radius is reconstructed very
accurately. The reason for this not being exactly one is probably the
pixellation of the camera. The standard deviation of the difference of
extracted opening angle and true opening angle is 0.0202 degrees. This method
showed that the ring center and radius from only using circle model from
skimage is not good enough even though ring center seemed to be reconstructed
well.

Figure A.2.1: Confusion matrix for opening angle with known center Figure
A.2.2: Reconstructed opening angle of the Cherenkov cone

##  Appendix B Found distributions opening angle

As seen before from figure 3.5.10 and 3.5.7 we see that Hough Transform
reconstructs ring parameters much more accurately than using the ringM method.
Comparing the distributions on real data (figure B.0.1), we can once again
see, that the distributions are different - ringM method reconstructs rings to
be with smaller radii while with the added Hough transform the reconstructed
radius is shifted towards the bigger radii.

Furthermore, from figure B.0.1 we can see that the reconstructed ring radius
for Hough transform peaks at opening angle $\sim 1.22$ degrees, which was also
expected (see 1.7). For ringM it peaks at $\sim 1$ degree.

The reason for Hough transform (DB-RM-HT) not being as smoothly distributed as
for DB-RM is most probably because of binning effects when performing the
Hough transform.

Figure B.0.1: Comparison of two models opening angle distribution. Both
distributions are normed