Thesis: SingleSine: rewrite algorithm part

documents/thesis/chapters/single_sine_interferometry.tex

@@ -122,7 +122,6 @@ This falls into the dynamic setup previously mentioned.
 The best period defects must thus be recovered from a single event.
 \\
 When doing the interferometric analysis for a sine beacon synchronised array, waveforms can only be delayed by an integer amount of periods, thereby giving discrete solutions to maximising the interferometric signal.
-\Todo{add size of shower at plane vs period defects in meters}
 
 \clearpage
 \section{Air Shower simulation}
@@ -212,6 +211,61 @@ The small clock defect $\tClockPhase$ is then finally calculated from the beacon
 From the above, we now have a set of reconstructed air shower waveforms with corresponding clock defects smaller than one beacon period $T$.
 Shifting the waveforms to remove these small clocks defects, we are left with resolving the correct number of periods $k$ per waveform (see Figure~\ref{fig:grid_power:repair_phases}).
 \\
+% >>>>
+\section{\textit{k}-finding} % <<<
+
+% unknown origin of air shower signal
+Up until now, the shower axis and thus the origin of the air shower signal have not been resolved.
+This means that the unknown propagation time delays for the air shower ($\tProp$) affect the alignment of the signals in Figure~\ref{fig:beacon_sync:period_alignment} in addition to the unknown clock period defects ($k T$).
+As such, both this origin and the clock defects have to be determined simultaneously.
+\\
+% radio interferometry
+If the antennas had been fully synchronised, radio interferometry as introduced in Chapter~\ref{sec:interferometry} can be applied to find the origin of the air shower signal, thus resolving the shower axis.
+Still, a (rough) first estimate of the shower axis might be made using this technique or by employing other detection techniques such as those using surface or fluorescence detectors.
+% (see Figure~\ref{fig:dynamic-resolve}).
+\\
+
+On the true shower axis, the waveforms would sum most coherently when the correct $k$'s are used.
+Therefore, around the estimated shower axis, we define a grid search to both optimise this sum and the location of the maximum power.
+In this process each waveform of the array is allowed to shift by a restricted amount of periods with respect to a reference waveform (taken to be the waveform with the highest maximum).
+\\
+
+Note that these grids are defined here in shower plane coordinates with $\vec{v}$ the true shower axis and $\vec{B}$ the local magnetic field.
+Searching a grid that is slightly misaligned with the true shower axis is expected to give comparable results.
+\\
+
+The steps of this $k$-finding algorithm are listed shortly below.
+Discussion is found in the next Chapter.
+\begin{enumerate}[start=1, label={Step \arabic*.}, ref=\arabic*, topsep=6pt, widest={Step 1.}]
+	\label{algo:kfinding}
+	\item \label{algo:kfinding:grid}
+		Define a grid around the estimated shower axis, zooming in on each iteration.
+
+	\item  \label{algo:kfinding:optimisation}
+		$k$-optimisation: per grid point, optimise the $k$'s to maximise the sum of the waveforms (see Figure~\ref{fig:single:k-correlation}).
+
+	\item  \label{algo:kfinding:kfinding}
+		$k$-finding: find the grid point with the maximum overall sum (see Figure~\ref{fig:findks:maxima}) and select its set of $k$'s.
+	
+	\item \label{algo:kfinding:break}
+		Stop when the set of $k$'s is equal to the set of the previous iteration.
+
+	\item \label{algo:kfinding:powermapping}
+		Finally, make a power mapping with the obtained $k$'s to determine re-estimate the shower axis (location with maximum power) (see Figure~\ref{fig:findks:reconstruction}), and return to Step~\ref{algo:kfinding:grid}.
+\end{enumerate}
+\vspace*{2pt}
+
+Here, Step~\ref{algo:kfinding:optimisation} has been implemented by summing each waveform to the reference waveform (see above) with different time delays $kT$ and selecting the $k$ that maximises the amplitude of a waveform combination.\footnote{%<<<
+	Note that one could use a correlation method instead of a maximum to select the best time delay.
+	However, for simplicity and ease of computation, this has not been implemented.
+	Other improvements are discussed in the next Section.
+} %>>>
+As shown in Figure~\ref{fig:single:k-correlation}, the maximum possible period shift has been limited to $\pm 3\,\mathrm{periods}$.
+This corresponds to the maximum expected time delay between two antennas with a clock randomisation up to $30\ns$ for the considered beacon frequency.%
+\footnote{
+	Figure~\ref{fig:simu:error:periods} shows this is not completely true.
+	However, overall, it still applies.
+}
 
 \begin{figure}%<<<
 	\centering
@@ -223,51 +277,12 @@ Shifting the waveforms to remove these small clocks defects, we are left with re
 	\label{fig:single:k-correlation}
 \end{figure}%>>>
 
-% >>>>
-\section{\textit{k}-finding} % <<<
-
-% unknown origin of air shower signal
-Up until now, the shower axis and thus the origin of the air shower signal have not been resolved.
-This means that the unknown propagation time delays for the air shower ($\tProp$) affect the alignment of the signals in Figure~\ref{fig:beacon_sync:period_alignment} in addition to the unknown clock period defects ($k_j T$).
-As such, both this origin and the clock defects have to be determined simultaneously.
-\\
-% radio interferometry
-If the antennas had been fully synchronised, radio interferometry as introduced in Section~\ref{sec:interferometry} can be applied to find the origin of the air shower signal, thus resolving the shower axis.
-Still, a rough first estimate of the shower axis might be made using this or other techniques.% (see Figure~\ref{fig:dynamic-resolve}).
-\\
-
-Starting with an initial grid around this estimated axis, a two-step process zooms in on the shower axis while optimising the interferometric signal.
-In this process each waveform of the array is allowed to shift by a restricted amount of periods with respect to a reference waveform (taken to be the waveform with the highest maximum).
-\\
-Note that these grids are defined in shower plane coordinates with $\vec{v}$ the true shower axis and $\vec{B}$ the local magnetic field.
-Searching a grid that is slightly misaligned with the true shower axis is expected to give comparable results.
-\\
-
-The first step consists of finding the set of period shifts $k_j$ that maximises the overall maximum amplitude on the grid.
-At each location, after removing propagation delays, each waveform and the reference waveform are summed with different time delays $kT$ to find the maximum attainable amplitude of the combined trace.%
-\footnote{%<<<
-	Note that one could use a correlation method instead of a maximum to select the best time delay.
-	However, for simplicity and ease of computation, this has not been implemented.
-%As shown in Figure~\ref{fig:single:annotated_full_waveform}, the air shower signal has a length in the order of a few nanoseconds.
-%Since it is this peak that is of interest, it would have been possible to cut the waveforms such to only correlate the peaks.
-} %>>>
-\\
-As shown in Figure~\ref{fig:single:k-correlation}, the maximum possible period shift has been limited to $\pm 3\,\mathrm{periods}$.
-This corresponds to the maximum expected time delay between two antennas with a clock randomisation up to $30\ns$ for the considered beacon frequency.
-\\
-%
-This amplitude optimisation is iterated over the grid (see Figure~\ref{fig:findks:maxima}) resulting in a grid measurement of the maximum attainable amplitude and the corresponding set of period defects $k$.
-\\
-
-The second step applies the obtained period defects and measures the interferometric power on the same grid (see Figure~\ref{fig:findks:reconstruction}).
-Afterwards, a new grid zooms in on the power maximum and the process is repeated (Figures~\ref{fig:findks:maxima:zoomed},~\ref{fig:findks:reconstruction:zoomed}) until the set of period defects is constant between zooming grids.
-\\
 
 \begin{figure}%<<< findks
 	\begin{subfigure}[t]{0.45\textwidth}
 		\includegraphics[width=\textwidth]{ZH_simulation/findks/ca_period_from_shower.py.maxima.run0.pdf}
 		\caption{
-			Maximum amplitudes obtainable by shifting the waveforms.
+			$k$-finding: optimise the $k$'s by shifting waveforms to find the maximum amplitude obtainable at each point.
 		}
 		\label{fig:findks:maxima}
 	\end{subfigure}
@@ -275,7 +290,7 @@ Afterwards, a new grid zooms in on the power maximum and the process is repeated
 	\begin{subfigure}[t]{0.45\textwidth}
 		\includegraphics[width=\textwidth]{ZH_simulation/findks/ca_period_from_shower.py.reconstruction.run0.power.pdf}
 		\caption{
-			Power measurement with the $k$'s belonging to the overall maximum of the amplitude maxima.
+			Power measurement with the $k$'s belonging to the overall maximum of the tested amplitudes.
 		}
 		\label{fig:findks:reconstruction}
 	\end{subfigure}
@@ -283,7 +298,8 @@ Afterwards, a new grid zooms in on the power maximum and the process is repeated
 	\begin{subfigure}[t]{0.45\textwidth}
 		\includegraphics[width=\textwidth]{ZH_simulation/findks/ca_period_from_shower.py.maxima.run1.pdf}
 		\caption{
-			Maximum amplitudes, zoomed to the location in \ref{fig:findks:reconstruction} with the highest amplitude.
+			$2^\mathrm{nd}$ $k$-finding iteration:
+				Zoom in on the location in \subref{fig:findks:reconstruction} with the highest amplitude and repeat algorithm.
 		}
 		\label{fig:findks:maxima:zoomed}
 	\end{subfigure}
@@ -299,7 +315,8 @@ Afterwards, a new grid zooms in on the power maximum and the process is repeated
 	\begin{subfigure}[t]{0.45\textwidth}
 		\includegraphics[width=\textwidth]{ZH_simulation/findks/ca_period_from_shower.py.maxima.run2.pdf}
 		\caption{
-			Final test grid obtaining the same $k$'s as Figure~\ref{fig:findks:maxima:zoomed}.
+			$3^\mathrm{rd}$ $k$-finding iteration:
+			The same set of $k$'s has been found and we stop the algorithm.
 		}
 		\label{fig:findks:maxima:zoomed2}
 	\end{subfigure}
@@ -312,7 +329,7 @@ Afterwards, a new grid zooms in on the power maximum and the process is repeated
 		\label{fig:findks:reconstruction:zoomed2}
 	\end{subfigure}
 	\caption{
-		Iterative $k$-finding algorithm:
+			Iterative $k$-finding algorithm (see page~\pageref{algo:kfinding} for explanation):
 		First \subref{fig:findks:maxima}, find the set of period shifts $k$ per point on a grid that returns the highest maximum amplitude (blue cross).
 		The grid starts as a $8^\circ$ wide shower plane slice at $X=400\,\mathrm{g/cm^2}$, centred at the true shower axis (red cross).
 		Second \subref{fig:findks:reconstruction}, perform the interferometric reconstruction with this set of period shifts.
@@ -332,22 +349,20 @@ Phase synchronising the antennas gives a small increase in observed power, while
 
 The initial grid plays an important role here in finding the correct axis.
 Due to selecting the highest maximum amplitude, and the process above zooming in aggressively, wrong candidate axes are selected when there is no grid-location sufficiently close to the true axis.
-Figure~\ref{fig:findks:reconstruction} shows such a potential point near $(-1, 0.5)$ with a maximum that is comparable to the selected maximum near the true axis.\Todo{no longer shown}
-\\
 % premature optimisation / degeneracy
 Such locations are subject to differences in propagation delays that are in the order of a few beacon periods.
 The restriction of the possible delays is therefore important to limit the number of potential axis locations.
 \\
-
 % fall in local extremum, maximum
-In this analysis, the initial grid is defined as $8^\circ$ wide around the true axis.
+In this analysis, the initial grid is defined as a very wide $8^\circ$ around the true axis.
 As the number of computations scales linearly with the number of grid points ($N = N_x N_y$), it is favourable to minimise the number of grid locations.
-Unfortunately, the above process has been observed to fall into local maxima when a too coarse initial grid ($N_x < 13$ at $X=400\,\mathrm{g/cm^2}$) was used while restricting the time delays to $\left| k \right| \leq 3$.
+Unfortunately, the above process has been observed to fall into local maxima when a too coarse and wide initial grid ($N_x < 13$ at $X=400\,\mathrm{g/cm^2}$) was used while restricting the time delays to $\left| k \right| \leq 3$.
 \\
 
 % Missing power / wrong k
 As visible in the right side of Figure~\ref{fig:grid_power:repair_full}, not all waveforms are in sync after the optimisation.
-In this case, the period defects have been resolved incorrectly for two waveforms, lagging 1 and 3 periods respectively (see Figure~\ref{fig:simu:error:periods}).
+In this case, the period defects have been resolved incorrectly for two waveforms (see Figure~\ref{fig:simu:error:periods}) due to too stringent limits on the allowable $k$'s.
+Looking at Figure~\ref{fig:grid_power:repair_phases}, this was to be predicted since there are two waveforms peaking at $k=4$ from the reference waveform's peak (dashed line).
 As a result, the obtained power for the resolved clock defects is slightly less than the obtained power for the true clocks.
 \\
 
@@ -357,17 +372,6 @@ Figure~\ref{fig:grid_power:axis} shows the power mapping at four different atmos
 Except for the low power case at $X=800\,\mathrm{g/cm^2}$, the shower axis is found to be $<0.1^\circ$ of the true shower axis.
 \\
 
-
-% Future: at multiple depths
-Of course, this algorithm must be evaluated at relevant atmospheric depths where the interferometric technique can resolve the air shower.
-In this case, after manual inspection, the air shower was found to have \Xmax\ at roughly $400\,\mathrm{g/cm^2}$.
-The algorithm is expected to perform as long as a region of strong coherent power is resolved.
-This means that with the power in both Figure~\ref{fig:grid_power:axis:X200} and Figure~\ref{fig:grid_power:axis:X600}, the clock defects and air shower should be identified to the same degree.
-\\
-Additionally, since the true period shifts are static per event, evaluating the $k$-finding algorithm at multiple atmospheric depths allows to compare the obtained sets thereof to further minimise any incorrectly resolved period defect.
-\\
-
-
 \begin{figure}% fig:simu:error
 	\centering
 	%\includegraphics[width=\textwidth]{ZH_simulation/cb_report_measured_antenna_offsets.py.time-amplitudes-missing-k.residuals.pdf}
@@ -438,6 +442,26 @@ Additionally, since the true period shifts are static per event, evaluating the
 	\label{fig:grid_power_time_fixes}
 \end{figure}%>>>
 
+\pagebreak
+% Future: at multiple depths
+Of course, this algorithm must be evaluated at relevant atmospheric depths where the interferometric technique can resolve the air shower.
+In this case, after manual inspection, the air shower was found to have \Xmax\ at roughly $400\,\mathrm{g/cm^2}$.
+The algorithm is expected to perform as long as a region of strong coherent power is resolved.
+This means that with the power in both Figure~\ref{fig:grid_power:axis:X200} and Figure~\ref{fig:grid_power:axis:X600}, the clock defects and air shower should be identified to the same degree.
+\\
+
+Additionally, since the true period shifts are static per event, evaluating the $k$-finding algorithm at multiple atmospheric depths allows to compare the obtained sets thereof to further minimise any incorrectly resolved period defect.
+\\
+
+Further improvements to the algorithm are foreseen in both the definition of the initial grid (Step~\ref{algo:kfinding:grid}) and the optimisation of the $k$'s (Step~\ref{algo:kfinding:optimisation}).
+For example, the $k$-optimisation step currently sums the full waveform for each $k$ to find the maximum amplitude for each sum.
+Instead, the timestamp of the amplitude maxima of each waveform can be compared, directly allowing to compute $k$ from the difference.
+\\
+
+Finally, from the overlapping traces in Figure~\ref{fig:grid_power:repair_full}, it is easily recognisable that some period defects have been determined incorrectly.
+Inspecting Figure~\ref{fig:grid_power:repair_phases}, this was to be expected as there are two waveforms with the peak at $\left|k\right| = 4$ from the reference waveform.
+Therefore, either the $k$-optimisation should have been run with a higher limit on the allowable $k$'s, or, preferably, these waveforms must be optimised after the algorithm is finished with a higher maximum $k$.
+
 \begin{figure}%<<< grid_power:axis:X600
 	\vspace*{-5mm}
 	\centering