source: trunk/MagicSoft/TDAS-Extractor/Algorithms.tex@ 5578

\section{Signal Reconstruction Algorithms \label{sec:algorithms}}

\ldots {\it In this section, the extractors are described, especially w.r.t. which free parameters are left to play, 
how they subtract the pedestal, how they compare between calibration and cosmics pulses and how an 
extraction in case of a pure pedestal event takes place. }
\newline
\newline
{\it Missing coding: 
\begin{itemize}
\item Implementing a low-gain extraction based on the high-gain information \ldots Arnau
\item Real fit to the expected pulse shape \ldots Hendrik, Wolfgang ???
\end{itemize}
}

\subsection{Pure signal extractors}

The pure signal extractors have in common that they compute only the 
signal, but no arrival time. All treated extractors here derive from the MARS-base
class {\textit{MExtractor}} which provides the following facilities:

\begin{itemize}
\item The global extraction limits can be set from outside
\item FADC saturation is kept track off
\end{itemize}

The following free adjustable parameters have to be set from outside: 
\begin{description}
\item[Global extraction limits:\xspace] Limits in between the extractor is allowed 
to search. 
\end{description}

\subsubsection{Fixed Window}

This extractor is implemented in the MARS-class {\textit{MExtractFixedWindow}}.
It simply adds the FADC contents in the allowed ranges. 
As it does not correct for the clock-noise, only an even number of samples is allowed.

\subsubsection{Fixed Window with global Peak Search}

This extractor is implemented in the MARS-class {\textit{MExtractFixedWindowPeakSearch}}.
It first fixes a reference point defined as the highest sum of 
consecutive non-saturating FADC slices in a (smaller) peak-search window. This reference 
point removes the coherent movement of the arrival times over whole camera due to the trigger jitter.

Then, simply adds the FADC contents around the reference point in a fixed window manner.
It loops twice over the all pixels every event, because it first has to find the reference point. 
As it does not correct for the clock-noise, only an even number of samples is allowed.

The following free adjustable parameters have to be set from outside: 
\begin{description}
\item[Peak Search Window:\xspace] Defines the ``sliding window'' in which the peaking sum is 
searched for (default: 4 slices)
\item[Offset from Window:\xspace] Defines the offset from the found reference point to start 
extracting the fixed window (default: 1 slice)
\item[Low-Gain Peak shift:\xspace] Defines the shift in the low-gain with respect to the peak found 
in the high-gain (default: 1 slice)
\end{description}

\subsubsection{Fixed Window with integrated cubic spline}

This extractor is implemented in the MARS-class {\textit{MExtractFixedWindowSpline}}.
It uses a cubic spline algorithm, adapted from \cite{NUMREC}. 
It integrated the 
spline interpolated FADC slice values, counting the edge slices as half.
As it does not correct for the clock-noise, only an odd number of samples is allowed.

\subsection{Combined extractors}

The combined extractors have in common that they compute the arrival time and 
the signal in one step. All treated combined extractors here derive from the MARS-base
class {\textit{MExtractTimeAndCharge}} which provides the following facilities:

\begin{itemize}
\item Only one loop over all pixels is performed
\item The individual FADC slice values get the clock-noise-corrected pedestals immediately subtracted.
\item The low-gain extraction range is adapted dynamically, based on the computed arrival time 
    from the high-gain samples
\item Extracted times from the low-gain samples get corrected for the intrinsic time delay of the low-gain
    pulse
\item The global extraction limits can be set from outside
\item FADC saturation is kept track off
\end{itemize}

The following free adjustable parameters have to be set from outside: 
\begin{description}
\item[Global extraction limits:\xspace] Limits in between the extractor is allowed 
to search. They are fixed by the extractor for the high-gain, but re-adjusted for 
every event in the low-gain, depending on the arrival time found in the low-gain. 
However, the dynamically adjusted window is not allowed to pass beyond the global 
limits.
\item[Low-gain start shift:\xspace] Global shift between the computed high-gain arrival
time and the start of the low-gain extraction limit (corrected for the intrinsic time offset). 
This variable tells where the extractor is allowed to start searching for the low-gain signal 
if the high-gain arrival time is known. It avoids that the extractor gets confused by possible high-gain 
signals leaking into the ``low-gain'' region.
\end{description}

\ldots {\it Note for the usage of this class together with {\textit{MJCalibration}}: In order to access the 
arrival times computed by these classes, the option: MJCalibration::SetTimeAndCharge() has to 
be chosen} \ldots 

\subsubsection{Sliding Window with amplitude-weighted time}

This extractor is implemented in the MARS-class {\textit{MExtractTimeAndChargeSlidingWindow}}.
It extracts the signal from a sliding window of an adjustable size, for high-gain and low-gain 
individually (default: 6 and 6) The signal is the one which maximizes the summed 
(clock-noise and pedestal-corrected) FADC signal over the window.
\par
The amplitude-weighted arrival time is calculated from the window with 
the highest integral using the following formula:

\begin{equation}
  t = \frac{\sum_{i=0}^{windowsize} s_i \cdot i}{\sum_{i=0}^{windowsize} i} 
\end{equation}
where $i$ denotes the FADC slice index, starting from the beginning of the derived 
window and running over the window and $s_i$ the clock-noise and 
pedestal-corrected FADC value at slice index i. 
\par
The following free adjustable parameters have to be set from outside: 
\begin{description}
\item[Window sizes:\xspace] Independenty for high-gain and low-gain (default: 6,6)
\end{description}

\subsubsection{Cubic Spline with Sliding Window or Amplitude extraction}

This extractor is implemented in the MARS-class {\textit{MExtractTimeAndChargeSpline}}.
It uses a cubic spline algorithm, adapted from \cite{NUMREC}. 
The following free adjustable parameters have to be set from outside: 

\begin{description}
\item[Time Extraction Type:\xspace] The position of the maximum can be chosen (default) or the 
position of the half maximum at the rising edge of the pulse
\item[Charge Extraction Type:\xspace] The amplitude of the maximum can be chosen (default) or the 
integrated spline between maximum position minus rise time (default: 1.5 slices) and maximum position plus
fall time (default: 4.5 slices). The low-gain signal integrates one slice more at the falling part of the 
signal.
\item[Rise Time and Fall Time:\xspace] Can be adjusted for the integration charge extraction type.
\item[Resolution:\xspace] Defined as the maximum allowed difference between the calculated half maximum value and 
the computed spline value at the arrival time position. Can be adjusted for the half-maximum time extraction 
type.
\end{description}

\subsubsection{Digital Filter}

This extractor is implemented in the MARS-class {\textit{MExtractTimeAndChargeDigitalFilter}}.


The goal of the digital filtering method \cite{OF94,OF77} is to optimally reconstruct the amplitude and time origin of a signal with a known signal shape from discrete measurements of the signal. Thereby the noise contribution to the amplitude reconstruction is minimized.

For the digital filtering method two assumptions have to be made:

\begin{itemize}
\item{The normalized signal shape has to be independent of the signal amplitude.}
\item{The noise properties have to be independent of the signal amplitude.}
\end{itemize}

Let $g(t)$ be the normalized signal shape, $E$ the signal amplitude and $\tau$ the time shift of the physical signal from the predicted signal shape. Then the time dependence of the signal, $y(t)$, is given by:

\begin{equation}
y(t)=E \cdot  g(t-\tau) + b(t) \ ,
\end{equation}

where $b(t)$ is the time dependent noise distribution. For small time shifts $\tau$ the time dependence can be linearized by the use of a Taylor expansion:

\begin{equation} \label{shape_taylor_approx}
y(t)=E \cdot  g(t) - E\tau \cdot  \dot{g}(t) + b(t) \ ,
\end{equation}

where $\dot{g}(t)$ is the time derivative of the signal shape. Discrete
measurements $y_i$ of the signal at times $t_i \ (i=1,...,n)$ have the form:

\begin{equation}
y_i=E \cdot g_i- E\tau \cdot \dot{g}_i +b_i \ .
\end{equation}

The correlation of the noise contributions at times $t_i$ and $t_j$ is expressed in the noise autocorrelation matrix $\boldsymbol{B}$:

\begin{equation} 
\boldsymbol{B_{ij}} = \langle b_i b_j \rangle - \langle b_i \rangle \langle b_j
\rangle  \ .
\label{eq:autocorr}
\end{equation}
%\equiv \langle b_i b_j \rangle with $\langle b_i \rangle = 0$.

The signal amplitude, $E$, and the product of amplitude and time shift, $E \tau$, can be estimated from the given set of
measurements $\boldsymbol{y} = (y_1, ... ,y_n)$ by minimizing:

\begin{eqnarray}
\chi^2(E, E\tau) &=& \sum_{i,j}(y_i-E g_i-E\tau \dot{g}_i) \boldsymbol{B}^{-1}_{ij} (y_j - E g_j-E\tau \dot{g}_i) \\
&=& (\boldsymbol{y} - E
\boldsymbol{g} - E\tau \dot{\boldsymbol{g}})^T \boldsymbol{B}^{-1} (\boldsymbol{y} - E \boldsymbol{g}- E\tau \dot{\boldsymbol{g}}) \ ,
\end{eqnarray}

where the last expression is matricial. The minimum is obtained for:

\begin{equation}
\frac{\partial \chi^2(E, E\tau)}{\partial E} = 0 \qquad \text{and} \qquad \frac{\partial \chi^2(E, E\tau)}{\partial(E\tau)} = 0 \ .
\end{equation}

This leads to the following two equations for the estimated amplitude $\overline{E}$ and the estimation for the product of amplitude and time offset $\overline{E\tau}$:

\begin{eqnarray}
0&=&-\boldsymbol{g}^T\boldsymbol{B}^{-1}\boldsymbol{y}+\boldsymbol{g}^T\boldsymbol{B}^{-1}\boldsymbol{g}\overline{E}+\boldsymbol{g}^T\boldsymbol{B}^{-1}\dot{\boldsymbol{g}}\overline{E\tau} 
\\
0&=&-\dot{\boldsymbol{g}}^T\boldsymbol{B}^{-1}\boldsymbol{y}+\dot{\boldsymbol{g}}^T\boldsymbol{B}^{-1}\boldsymbol{g}\overline{E}+\dot{\boldsymbol{g}}^T\boldsymbol{B}^{-1}\dot{\boldsymbol{g}}\overline{E\tau} \ .
\end{eqnarray}

Solving these equations one gets the solutions:

\begin{equation}
\overline{E}= \boldsymbol{w}_{\text{amp}}^T \boldsymbol{y} \quad \mathrm{with} \quad \boldsymbol{w}_{\text{amp}} = \frac{ (\dot{\boldsymbol{g}}^T\boldsymbol{B}^{-1}\dot{\boldsymbol{g}}) \boldsymbol{B}^{-1} \boldsymbol{g} -(\boldsymbol{g}^T\boldsymbol{B}^{-1}\dot{\boldsymbol{g}})  \boldsymbol{B}^{-1} \dot{\boldsymbol{g}}}  {(\boldsymbol{g}^T \boldsymbol{B}^{-1} \boldsymbol{g})(\dot{\boldsymbol{g}}^T\boldsymbol{B}^{-1}\dot{\boldsymbol{g}}) -(\dot{\boldsymbol{g}}^T\boldsymbol{B}^{-1}\boldsymbol{g})^2 } \ ,
\end{equation}

\begin{equation}
\overline{E\tau}= \boldsymbol{w}_{\text{time}}^T \boldsymbol{y} \quad \mathrm{with} \quad \boldsymbol{w}_{\text{time}} = \frac{ ({\boldsymbol{g}}^T\boldsymbol{B}^{-1}{\boldsymbol{g}}) \boldsymbol{B}^{-1} \dot{\boldsymbol{g}} -(\boldsymbol{g}^T\boldsymbol{B}^{-1}\dot{\boldsymbol{g}})  \boldsymbol{B}^{-1} {\boldsymbol{g}}}  {(\boldsymbol{g}^T \boldsymbol{B}^{-1} \boldsymbol{g})(\dot{\boldsymbol{g}}^T\boldsymbol{B}^{-1}\dot{\boldsymbol{g}}) -(\dot{\boldsymbol{g}}^T\boldsymbol{B}^{-1}\boldsymbol{g})^2 } \ .
\end{equation}


Thus $\overline{E}$ and $\overline{E\tau}$ are given by a weighted sum of the discrete measurements $y_i$ with the digital filtering weights for the amplitude, $w_{\text{amp},i}$, and time shift, $w_{\text{time},i}$.

Because of the truncation of the Taylor series in equation (\ref{shape_taylor_approx}) the above results are only valid for vanishing time offsets $\tau$. For non-zero time offsets one has to iterate the problem using the time shifted signal shape $g(t-\tau)$.

The covariance matrix $\boldsymbol{V}$ of $\overline{E}$ and $\overline{E\tau}$ is given by:

\begin{equation}
\boldsymbol{V}^{-1}=\frac{1}{2}\left(\frac{\partial^2 \chi^2(E, E\tau)}{\partial \alpha_i \partial \alpha_j} \right) \quad \text{with} \quad \alpha_i,\alpha_j \in [E, E\tau] \ .
\end{equation}

The expected contribution of the noise to the estimated amplitude, $\sigma_E$, is:

\begin{equation}\label{of_noise}
\sigma_E^2=\frac{\dot{\boldsymbol{g}}^T\boldsymbol{B}^{-1}\dot{\boldsymbol{g}}}{(\boldsymbol{g}^T \boldsymbol{B}^{-1} \boldsymbol{g})(\dot{\boldsymbol{g}}^T\boldsymbol{B}^{-1}\dot{\boldsymbol{g}}) -(\dot{\boldsymbol{g}}^T\boldsymbol{B}^{-1}\boldsymbol{g})^2} \ .
\end{equation}

The expected contribution of the noise to the estimated timing ,$\sigma_{\tau}$, is:

\begin{equation}\label{of_noise_time}
E^2 \cdot \sigma_{\tau}^2=\frac{{\boldsymbol{g}}^T\boldsymbol{B}^{-1}{\boldsymbol{g}}}{(\boldsymbol{g}^T \boldsymbol{B}^{-1} \boldsymbol{g})(\dot{\boldsymbol{g}}^T\boldsymbol{B}^{-1}\dot{\boldsymbol{g}}) -(\dot{\boldsymbol{g}}^T\boldsymbol{B}^{-1}\boldsymbol{g})^2} \ .
\end{equation}

For the MAGIC signals, as implemented in the MC simulations, a pedestal RMS of a single FADC slice of 4 FADC counts introduces an error in the reconstructed signal and time of:

\begin{equation}\label{of_noise}
\sigma_E \approx 8.3 \ \mathrm{FADC\ counts} \qquad \sigma_{\tau}  \approx  \frac{6.5\ \Delta T_{\mathrm{FADC}}}{E\ /\ \mathrm{FADC\ counts}} \ ,
\end{equation}

where $\Delta T_{\mathrm{FADC}} = 3.33$ ns is the sampling intervall of the MAGIC FADCs.

In the current implementation a two step procedure is applied to reconstruct the signal. The weight functions $w_{\mathrm{amp}}(t)$ and $w_{\mathrm{time}}(t)$ are computed numerically with a resolution of $1/10$ of an FADC slice. In the first step the quantities $e_{i_0}$ and $e\tau_{i_0}$ are computed:

\begin{equation}
e_{i_0}=\sum_{i=i_0}^{i=i_0+5} w_{\mathrm{amp}(t_i)}y(t_{i+i_0}) \qquad e\tau_{i_0}=\sum_{i=i_0}^{i=i_0+5} w(t_i)_{\mathrm{time}(t_i)}y(t_{i+i_0}) 
\end{equation}

for all possible signal start samples $i_0$. Let $i_0^*$ be the signal start sample with the largest $e_{i_0}$. Then in a seconst step the timing offset $\tau$ is calculated:

\begin{equation}
\tau=\frac{e\tau_{i_0^*}}{e_{i_0^*}}
\end{equation}

and the weigths iterated:

\begin{equation}
E_{i_0^*}=\sum_{i=i_0^*}^{i=i_0^*+5} w_{\mathrm{amp}(t_i - \tau)}y(t_{i+i_0^*}) \qquad E \tau_{i_0^*}=\sum_{i=i_0^*}^{i=i_0^*+5} w(t_i - \tau)_{\mathrm{time}(t_i)}y(t_{i+i_0^*}) \ .
\end{equation}

The reconstructed signal is then taken to be $E_{i_0^*}$ and the reconstructed arrival time $t_{\text{arrival}}$ is 

\begin{equation}
t_{\text{arrival}} = i_0^* + \tau_{i_0^*} + \tau \ .
\end{equation}



% This does not apply for MAGIC as the LONs are giving always a correlated noise (in addition to the artificial shaping)

%In the case of an uncorrelated noise with zero mean the noise autocorrelation matrix is:

%\begin{equation}
%\boldsymbol{B}_{ij}= \langle b_i b_j \rangle \delta_{ij} = \sigma^2(b_i) \ ,
%\end{equation}

%where $\sigma(b_i)$ is the standard deviation of the noise of the discrete measurements. Equation (\ref{of_noise}) than becomes:


%\begin{equation}
%\frac{\sigma^2(b_i)}{\sigma_E^2} = \sum_{i=1}^{n}{g_i^2} - \frac{\sum_{i=1}^{n}{g_i \dot{g}_i}}{\sum_{i=1}^{n}{\dot{g}_i^2}} \ .
%\end{equation}



\ldots {\it Hendrik ... }

The following free adjustable parameters have to be set from outside: 

\begin{description}
\item[Weights File:\xspace] 
\item[Window Sizes:\xspace]
\item[Binning Resolution:\xspace]
\end{description}

\subsubsection{Real fit to the expected pulse shape }

This extractor is not yet implemented as MARS-class...
\par
It fit the pulse shape to a Landau convoluted with a Gaussian using the following 
parameters:...

\ldots {\it Hendrik, Wolfgang ... }


\subsection{Used Extractors for this Analysis}

We propose to test the following parameterized extractors:

\begin{description}
\item[MExtractFixedWindow]: with the following parameters, if {\textit{maxbin}} defines the 
   mean position of the High-Gain FADC slice carrying the pulse maximum:
\begin{enumerate}
\item SetRange({\textit{maxbin}}-1,{\textit{maxbin}}+2,{\textit{maxbin}}+0.5,{\textit{maxbin}}+3.5);
\item SetRange({\textit{maxbin}}-1,{\textit{maxbin}}+2,{\textit{maxbin}}-0.5,{\textit{maxbin}}+4.5);
\item SetRange({\textit{maxbin}}-2,{\textit{maxbin}}+3,{\textit{maxbin}}-0.5,{\textit{maxbin}}+4.5);
\item SetRange({\textit{maxbin}}-3,{\textit{maxbin}}+4,{\textit{maxbin}}-1.5,{\textit{maxbin}}+5.5);
\item SetRange({\textit{maxbin}}-5,{\textit{maxbin}}+8,{\textit{maxbin}}-1.5,{\textit{maxbin}}+7.5);
\suspend{enumerate}
\item[MExtractFixedWindowSpline]: with the following parameters, if {\textit{maxbin}} defines the 
   mean position of the High-Gain FADC slice carrying the pulse maximum:
\resume{enumerate}
\item SetRange({\textit{maxbin}}-1,{\textit{maxbin}}+2,{\textit{maxbin}}+0.5,{\textit{maxbin}}+3.5);
\item SetRange({\textit{maxbin}}-1,{\textit{maxbin}}+2,{\textit{maxbin}}-0.5,{\textit{maxbin}}+4.5);
\item SetRange({\textit{maxbin}}-2,{\textit{maxbin}}+3,{\textit{maxbin}}-0.5,{\textit{maxbin}}+4.5);
\item SetRange({\textit{maxbin}}-3,{\textit{maxbin}}+4,{\textit{maxbin}}-1.5,{\textit{maxbin}}+5.5);
\item SetRange({\textit{maxbin}}-5,{\textit{maxbin}}+8,{\textit{maxbin}}-1.5,{\textit{maxbin}}+7.5);
\suspend{enumerate}
\item[MExtractFixedWindowPeakSearch]: with the following parameters: \\
SetRange(0,18,2,14); and:
\resume{enumerate}
\item SetWindows(2,2,2);
\item SetWindows(4,4,2);
\item SetWindows(6,6,4);
\item SetWindows(4,6,4);
\item SetWindows(8,8,4);
\item SetWindows(14,10,4);
\suspend{enumerate}
\item[MExtractTimeAndChargeSlidingWindow]: with the following parameters: \\
SetRange(0,18,2,14); and:
\resume{enumerate}
\item SetWindowSize(2,2);
\item SetWindowSize(4,4);
\item SetWindowSize(4,6);
\item SetWindowSize(6,6);
\item SetWindowSize(8,8);
\item SetWindowSize(14,10);
\suspend{enumerate}
\item[MExtractTimeAndChargeSpline]: with the following parameters:\\
SetChargeType(MExtractTimeAndChargeSpline::kIntegral); \\
SetRange(0,18,2,14); \\
and:

\resume{enumerate}
\item SetRiseTime(1.5); SetFallTime(4.5);
\item SetRiseTime(0.5); SetFallTime(2.5);
\item SetRiseTime(0.5); SetFallTime(1.5);
\item SetRiseTime(0.5); SetFallTime(0.5);
\suspend{enumerate}
and:
\resume{enumerate}
\item SetChargeType(MExtractTimeAndChargeSpline::kAmplitude); \\
SetRange(0,10,4,11); and:
\suspend{enumerate}
\item[MExtractTimeAndChargeDigitalFilter]: with the following parameters:
\resume{enumerate}
\item SetWeightsFile(``cosmic\_weights6.dat'');
\item SetWeightsFile(``cosmic\_weights4.dat'');
\item SetWeightsFile(``cosmic\_weights\_logain6.dat'');
\item SetWeightsFile(``cosmic\_weights\_logain4.dat'');
\item SetWeightsFile(``calibration\_weights\_UV6.dat'');
\item SetWeightsFile(``calibration\_weights\_UV\_logain6.dat'');
\suspend{enumerate}
\item[``Real Fit'']: (not yet implemented, one try)
\resume{enumerate}
\item Real Fit
\end{enumerate}
\end{description}


References: \cite{OF77,OF94}.


%%% Local Variables: 
%%% mode: latex
%%% TeX-master: "MAGIC_signal_reco"
%%% TeX-master: "Algorithms"
%%% TeX-master: "MAGIC_signal_reco"
%%% TeX-master: "MAGIC_signal_reco"
%%% TeX-master: "MAGIC_signal_reco"
%%% TeX-master: "MAGIC_signal_reco"
%%% End: