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

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

{\it Missing coding: 
\begin{itemize}
\item Real fit to the expected pulse shape \ldots Hendrik, Wolfgang ???
\end{itemize}
}

\subsection{Implementation of Signal Extractors in MARS}

All signal extractor classes are stored in the MARS-directory {\textit{\bf msignal/}}. 
There, the base classes {\textit{\bf MExtractor}}, {\textit{\bf MExtractTime}}, {\textit{\bf MExtractTimeAndCharge}} and 
all individual extractors can be found. Figure~\ref{fig:extractorclasses} gives a sketch of the 
inheritances of each class and what each class calculates.

\begin{figure}[htp]
\includegraphics[width=0.99\linewidth]{ExtractorClasses.eps}
\caption{Sketch of the inheritances of three examplary MARS signal extractor classes: MExtractFixedWindow, MExtractTimeFastSpline and MExtractTimeAndChargeDigitalFilter}
\label{fig:extractorclasses}
\end{figure}

The following base classes for the extractor tasks are used:
\begin{description}
\item[MExtractor:\xspace] This class provides the basic data members equal for all extractors which are:
  \begin{enumerate}
  \item Global extraction ranges, parameterized by the variables 
    {\textit{\bf fHiGainFirst, fHiGainLast, fLoGainFirst, fLoGainLast}} and the function {\textit{\bf SetRange()}}. 
    The ranges always {\textit{\bf include}} the edge slices. 
  \item An internal variable {\textit{\bf fHiLoLast}} regulating the overlap of the desired high-gain
    extraction range into the low-gain array.
  \item The maximum possible FADC value, before the slice is declared as saturated, parameterized 
    by the variable {\textit{\bf fSaturationLimit}} (default:\,254). 
  \item The typical delay between high-gain and low-gain slices, expressed in FADC slices and parameterized
    by the variable {\textit{\bf fOffsetLoGain}} (default:\,1.51)
  \item Pointers to the used storage containers {\textit{\bf MRawEvtData, MRawRunHeader, MPedestalCam}} 
    and~{\textit{\bf MExtractedSignalCam}}, parameterized by the variables 
    {\textit{\bf fRawEvt, fRunHeader, fPedestals}} and~{\textit{\bf fSignals}}.
  \item Names of the used storage containers to be searched for in the parameter list, parameterized 
    by the variables {\textit{\bf fNamePedestalCam}} and~{\textit{\bf fNameSignalCam}} (default: ``MPedestalCam'' 
    and~''MExtractedSignalCam'').
  \item The equivalent number of FADC samples, used for the calculation of the pedestal RMS and then the 
    number of photo-electrons with the F-Factor method (see eq.~\ref{eq:rmssubtraction} and 
    section~\ref{sec:photo-electrons}). This number is parameterized by the variables 
    {\textit{\bf fNumHiGainSamples}} and~{\textit{\bf fNumLoGainSamples}}.
  \end{enumerate}

  {\textit {\bf MExtractor}} is able to loop over all events, if the {\textit{\bf Process()}}-function is not overwritten. 
  It uses the following (virtual) functions, to be overwritten by the derived extractor class:

  \begin{enumerate}
  \item void {\textit {\bf FindSignalHiGain}}(Byte\_t* firstused, Byte\_t* logain, Float\_t\& sum, Byte\_t\& sat) const
  \item void {\textit {\bf FindSignalLoGain}}(Byte\_t* firstused, Float\_t\& sum, Byte\_t\& sat) const
  \end{enumerate}

  where the pointers ``firstused'' point to the first used FADC slice declared by the extraction ranges,
  the pointer ``logain'' points to the beginning of the ``low-gain'' FADC slices array (to be used for 
  pulses reaching into the low-gain array) and the variables ``sum'' and ``sat'' get filled with the 
  extracted signal and the number of saturating FADC slices, respectively.
  \par
  The pedestals get subtracted automatically {\textit {\bf after}} execution of these two functions.

\item[MExtractTime:\xspace] This class provides - additionally to those already declared in {\textit{\bf MExtractor}} - 
  the basic data members equal for all time extractors which are:
  \begin{enumerate}
  \item Pointer to the used storage container {\textit{\bf MArrivalTimeCam}}
    parameterized by the variables 
    {\textit{\bf fArrTime}}.
  \item The name of the used ``MArrivalTimeCam''-container to be searched for in the parameter list, 
    parameterized by the variables {\textit{\bf fNameTimeCam}} (default: ``MArrivalTimeCam'' ).
  \end{enumerate}

  {\textit {\bf MExtractTime}} is able to loop over all events, if the {\textit{\bf Process()}}-function is not 
  overwritten. 
  It uses the following (virtual) functions, to be overwritten by the derived extractor class:

  \begin{enumerate}
  \item void {\textit {\bf FindTimeHiGain}}(Byte\_t* firstused, Float\_t\& time, Float\_t\& dtime, Byte\_t\& sat, const MPedestlPix \&ped) const
  \item void {\textit {\bf FindTimeLoGain}}(Byte\_t* firstused, Float\_t\& time, Float\_t\& dtime, Byte\_t\& sat, const MPedestalPix \&ped) const
  \end{enumerate}

  where the pointers ``firstused'' point to the first used FADC slice declared by the extraction ranges,
  and the variables ``time'', ``dtime'' and ``sat'' get filled with the 
  extracted arrival time, its error and the number of saturating FADC slices, respectively.
  \par
  The pedestals can be used for the arrival time extraction via the reference ``ped''.

\item[MExtractTimeAndCharge:\xspace] This class provides - additionally to those already declared in 
  {\textit{\bf MExtractor}} and {\textit{\bf MExtractTime}} - 
  the basic data members equal for all time and charge extractors which are:
  \begin{enumerate}
  \item The actual extraction window sizes, parameterized by the variables 
    {\textit{\bf fWindowSizeHiGain}} and {\textit{\bf fWindowSizeLoGain}}.
  \item The shift of the low-gain extraction range start w.r.t. to the found high-gain arrival 
    time, parameterized by the variable {\textit{\bf fLoGainStartShift}} (default: -2.8)
  \end{enumerate}

  {\textit {\bf MExtractTimeAndCharge}} is able to loop over all events, if the {\textit{\bf Process()}}-function is not 
  overwritten. 
  It uses the following (virtual) functions, to be overwritten by the derived extractor class:

  \begin{enumerate}
  \item void {\textit {\bf FindTimeAndChargeHiGain}}(Byte\_t* firstused, Byte\_t* logain, Float\_t\& sum, Float\_t\& dsum, 
    Float\_t\& time, Float\_t\& dtime, Byte\_t\& sat, 
    const MPedestlPix \&ped, const Bool\_t abflag) const
  \item void {\textit {\bf FindTimeAndChargeLoGain}}(Byte\_t* firstused, Float\_t\& sum, Float\_t\& dsum, 
    Float\_t\& time, Float\_t\& dtime, Byte\_t\& sat, 
    const MPedestalPix \&ped, const Bool\_t abflag) const
  \end{enumerate}

  where the pointers ``firstused'' point to the first used FADC slice declared by the extraction ranges,
  the pointer ``logain'' point to the beginning of the low-gain FADC slices array (to be used for 
  pulses reaching into the ``low-gain'' array),
  the variables ``sum'', ``dsum'' get filled with the 
  extracted signal and its error. The variables ``time'', ``dtime'' and ``sat'' get filled with the 
  extracted arrival time, its error and the number of saturating FADC slices, respectively. 
  \par
  The pedestals can be used for the extraction via the reference ``ped'', also the AB-flag is given 
  for AB-clock noise correction.
\end{description}
  
  
\subsection{Pure Signal Extractors}

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

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

The following  adjustable parameters have to be set from outside: 
\begin{description}
\item[Global extraction limits:\xspace] Limits in between which the extractor is allowed 
to extract the signal, for high gain and low gain, respectively. 
\end{description}

As the pulses jitter by about one FADC slice, 
not every pulse lies exactly within the optimal limits, especially if one takes small 
extraction windows. 
Moreover, the readout position with respect to the trigger position has changed a couple 
of times during last year, therefore a very careful adjustment of the extraction limits
is mandatory before using these extractors. 

\subsubsection{Fixed Window}

This extractor is implemented in the MARS-class {\textit{\bf MExtractFixedWindow}}.
It simply adds the FADC slice contents in the assigned ranges. 
As it does not correct for the clock-noise, only an even number of samples is allowed.
Figure~\ref{fig:fixedwindowsketch} gives a sketch of the used extraction ranges for this 
paper and two typical calibration pulses. 

\begin{figure}[htp]
  \includegraphics[width=0.49\linewidth]{MExtractFixedWindow_5Led_UV.eps}
  \includegraphics[width=0.49\linewidth]{MExtractFixedWindow_23Led_Blue.eps}
\caption[Sketch extraction ranges MExtractFixedWindow]{%
Sketch of the extraction ranges for the extractor {\textit{\bf MExtractFixedWindow}} 
for two typical calibration pulses (pedestals have been subtracted) and a typical inner pixel.
The pulse would be shifted half a slice to the right for an outer pixel. }
\label{fig:fixedwindowsketch}
\end{figure}


\subsubsection{Fixed Window with Integrated Cubic Spline}

This extractor is implemented in the MARS-class {\textit{\bf MExtractFixedWindowSpline}}. It
uses a cubic spline algorithm, adapted from \cite{NUMREC} and integrates the 
spline interpolated FADC slice values from a fixed extraction range. The edge slices are counted as half.
As it does not correct for the clock-noise, only an odd number of samples is allowed.
Figure~\ref{fig:fixedwindowsplinesketch} gives a sketch of the used extraction ranges for this 
paper and two typical calibration pulses. 

\begin{figure}[htp]
  \includegraphics[width=0.49\linewidth]{MExtractFixedWindowSpline_5Led_UV.eps}
  \includegraphics[width=0.49\linewidth]{MExtractFixedWindowSpline_23Led_Blue.eps}
\caption[Sketch extraction ranges MExtractFixedWindowSpline]{%
Sketch of the extraction ranges for the extractor {\textit{\bf MExtractFixedWindowSpline}} 
for two typical calibration pulses (pedestals have been subtracted) and a typical inner pixel. 
The pulse would be shifted half a slice to the right for an outer pixel. }
\label{fig:fixedwindowsplinesketch}
\end{figure}

\subsubsection{Fixed Window with Global Peak Search}

This extractor is implemented in the MARS-class {\textit{\bf MExtractFixedWindowPeakSearch}}.
The basic idea of this extractor is to correct for coherent movements in arrival time for all pixels, 
as e.g. caused by the trigger jitter.  
In a first loop, it fixes a reference point defined as the highest sum of 
consecutive non-saturating FADC slices in a (smaller) peak-search window.
\par
In a second loop over the pixels, 
it adds the FADC contents starting from a pre-defined offset from the obtained peak-search window 
over an extraction window of a pre-defined window size.
It loops twice over all pixels in every event, because it has to find the reference point, first. 
As it does not correct for the clock-noise, only an even number of samples is allowed.
For a high intensity calibration run causing high-gain saturation in the whole camera, this 
extractor apparently fails since only dead pixels which cannot produced a saturated signal, are taken into account 
in the peak search. For this special case, we modified {\textit{\bf MExtractFixedWindowPeakSearch}} 
such to define the peak search window as the one starting from the mean position of the first saturating slice.
\par
The following  adjustable parameters have to be set from outside: 
\begin{description}
\item[Peak Search Window:\xspace] Defines the ``sliding window'' size within which the peaking sum is 
searched for (default: 4 slices)
\item[Offset from Window:\xspace] Defines the offset of the start of the extraction window w.r.t. the 
starting point of the obtained peak search 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}

Figure~\ref{fig:fixedwindowpeaksearchsketch} gives a sketch of the possible peak-search and extraction 
window positions in two typical calibration pulses. 

\begin{figure}[htp]
  \includegraphics[width=0.49\linewidth]{MExtractFixedWindowPeakSearch_5Led_UV.eps}
  \includegraphics[width=0.49\linewidth]{MExtractFixedWindowPeakSearch_23Led_Blue.eps}
\caption[Sketch extraction ranges MExtractFixedWindowPeakSearch]{%
Sketch of the extraction ranges for the extractor {\textit{\bf MExtractFixedWindowPeakSearch}} 
for two typical calibration pulses (pedestals have been subtracted) and a typical inner pixel.
The pulse would be shifted half a slice to the right for an outer pixel. }
\label{fig:fixedwindowpeaksearchsketch}
\end{figure}

\subsection{Combined Extractors}

The combined extractors have in common that they reconstruct the arrival time and 
the charge at the same time and for the same pulse. 
All treated combined extractors here derive from the MARS-base
class {\textit{\bf MExtractTimeAndCharge}} which itself derives from MExtractor and MExtractTime. 
It 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 of.
\end{itemize}

The following adjustable parameters have to be set from outside, additionally to those declared in the 
base classes MExtractor and MExtractTime:
 
\begin{description}
\item[Global extraction limits:\xspace] Limits in between which 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 (default: -2.8).
\end{description}

\subsubsection{Sliding Window with Amplitude-Weighted Time}

This extractor is implemented in the MARS-class {\textit{\bf 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) consecutive FADC slice contents.
\par
The amplitude-weighted arrival time is calculated from the window with 
the highest FADC slice contents integral using the following formula:

\begin{equation}
  t = \frac{\sum_{i=i_0}^{i_0+ws} s_i \cdot i}{\sum_{i=i_0}^{i_0t+ws} i} 
\end{equation}
where $i$ denotes the FADC slice index, starting from $i_0$
window and running over a window of size $ws$. $s_i$ the clock-noise and 
pedestal-corrected FADC slice contents at slice position $i$. 
\par
The following adjustable parameters have to be set from outside: 
\begin{description}
\item[Window sizes:\xspace] Independently for high-gain and low-gain (default: 6,6)
\end{description}

\begin{figure}[htp]
  \includegraphics[width=0.49\linewidth]{MExtractTimeAndChargeSlidingWindow_5Led_UV.eps}
  \includegraphics[width=0.49\linewidth]{MExtractTimeAndChargeSlidingWindow_23Led_Blue.eps}
\caption[Sketch calculated arrival times MExtractTimeAndChargeSlidingWindow]{%
Sketch of the calculated arrival times for the extractor {\textit{\bf MExtractTimeAndChargeSlidingWindow}} 
for two typical calibration pulses (pedestals have been subtracted) and a typical inner pixel. 
The extraction window sizes modify the position of the (amplitude-weighted) mean FADC-slices slightly.
The pulse would be shifted half a slice to the right for an outer pixel. }
\label{fig:slidingwindowsketch}
\end{figure}

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

This extractor is implemented in the MARS-class {\textit{\bf MExtractTimeAndChargeSpline}}.
It interpolates the FADC contents using a cubic spline algorithm, adapted from \cite{NUMREC}. 
In a second step, it searches for the position of the spline maximum. From then on, two 
possibilities are offered:

\begin{enumerate}
\item[Extraction Type Amplitude:\xspace] The amplitude of the spline maximum is taken as charge signal
and the (precisee) position of the maximum is returned as arrival time. This type is faster, since it 
performs not spline intergraion . 
\item[Extraction Type Integral:\xspace] The integrated spline between maximum position minus 
rise time (default: 1.5 slices) and maximum position plus fall time (default: 4.5 slices) 
is taken as charge signal and the position of the half maximum left from the position of the maximum 
is returned as arrival time (default). 
The low-gain signal stretches the rise and fall time by a stretch factor (default: 1.5). This type 
is slower, but yields more precise results (see section~\ref{sec:performance}) . 
The charge integration resolution is set to 0.1 FADC slices.
\end{enumerate}

The following adjustable parameters have to be set from outside: 

\begin{description}
\item[Charge Extraction Type:\xspace] The amplitude of the spline maximum can be chosen while the position 
of the maximum is returned as arrival time. This type is fast. \\
Otherwise, the integrated spline between maximum position minus rise time (default: 1.5 slices) 
and maximum position plus fall time (default: 4.5 slices) is taken as signal and the position of the 
half maximum is returned as arrival time (default). 
The low-gain signal stretches the rise and fall time by a stretch factor (default: 1.5). This type 
is slower, but more precise. The charge integration resolution is 0.1 FADC slices.
\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.
\item[Low Gain Stretch:\xspace] Can be adjusted to account for the larger rise and fall times in the 
low-gain as compared to the high gain pulses (default: 1.5)
\end{description}

\begin{figure}[htp]
  \includegraphics[width=0.49\linewidth]{MExtractTimeAndChargeSpline_5Led_UV.eps}
  \includegraphics[width=0.49\linewidth]{MExtractTimeAndChargeSpline_23Led_Blue.eps}
\caption[Sketch calculated arrival times MExtractTimeAndChargeSpline]{%
Sketch of the calculated arrival times for the extractor {\textit{\bf MExtractTimeAndChargeSpline}} 
for two typical calibration pulses (pedestals have been subtracted) and a typical inner pixel. 
The extraction window sizes modify the position of the (amplitude-weighted) mean FADC-slices slightly.
The pulse would be shifted half a slice to the right for an outer pixel. }
\label{fig:splinesketch}
\end{figure}

\subsubsection{Digital Filter}

This extractor is implemented in the MARS-class {\textit{\bf 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, three 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.}
\item{The noise auto-correlation matrix does not change its form significantly with time.}
\end{itemize}

\par
\ldots {\textit{\bf IS THIS TRUE FOR MAGIC???? }} \ldots
\par

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 contribution. For small time shifts $\tau$ (usually smaller than 
one FADC slice width), 
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$ can be expressed in the 
noise autocorrelation matrix $\boldsymbol{B}$:

\begin{equation} 
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 the excess noise contribution with respect to the known noise 
auto-correlation:

\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}_j) \\
&=& (\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. $\chi^2$ is a continuous function of $\tau$ and will have to be discretized itself for a 
desired resolution.  $\chi^2$ is also proportional to the auto noise-correlation matrix where increases in the noise level lead to 
a multiplicative factor for all matrix elements and thus do not affect the position of the maximum of $\chi^2$.
The minimum of $\chi^2$ 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}


Taking into account that $\boldsymbol{B}$ is a symmetric matrix, 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 following solutions:

\begin{equation}
\overline{E}(\tau) = \boldsymbol{w}_{\text{amp}}^T (\tau)\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}(\tau)= \boldsymbol{w}_{\text{time}}^T(\tau) \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}}(\tau)$, and time shift, $w_{\text{time}}(\tau)$
where the time dependency gets discretized once again leading to a set of weights samples depending on the 
discretized time $\tau$.

\par
\ldots {\textit{\bf IS THIS CORRECT LIKE THIS???}} \ldots
\par

Note the remaining time dependence of the two weights which follow from the dependency of $\boldsymbol{g}$ and 
$\dot{\boldsymbol{g}}$ on the position of the pulse with respect to the FADC bin positions.
\par
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}
\left(\boldsymbol{V}^{-1}\right)_{i,j}=\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=\boldsymbol{V}_{E,E}=\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=\boldsymbol{V}_{E\tau,E\tau}=\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 interval of the MAGIC FADCs.


For an IACT there are two types of background noise. On the one hand there is the constantly present electronics noise, on the other hand the light of the night sky introduces a sizeable background noise to the measurement of Cherenkov photons from air showers.

The electronics noise is largely white, uncorrelated in time. The noise from the night sky background photons is the superposition of the detector response to single photo electrons arriving randomly distributed in time. Figure \ref{fig:noise_autocorr_AB_36038_TDAS} shows the noise autocorrelation matrix for an open camera. The large noise autocorrelation in time of the current FADC system is due to the pulse shaping with a shaping constant of 6 ns. 

In general the amplitude and timing weights, $\boldsymbol{w}_{\text{amp}}$ and $\boldsymbol{w}_{\text{time}}$, depend on the pulse shape, the derivative of the pulse shape and the noise autocorrelation. In the high gain samples the correlated night sky background noise dominates over the white electronics noise. Thus different noise levels just cause the noise autocorrelation matrix $\boldsymbol{B}$ to change by a factor, which cancels out in the weights calculation. Thus in the high gain the weights are to a very good approximation independent of the night sky background noise level.

Contrary to that in the low gain samples ... .



\begin{figure}[h!]
\begin{center}
\includegraphics[totalheight=7cm]{noise_autocorr_AB_36038_TDAS.eps}
\end{center}
\caption[Noise autocorrelation.]{Noise autocorrelation matrix for open camera including the noise due to night sky background fluctuations.} \label{fig:noise_autocorr_AB_36038_TDAS}
\end{figure}



Using the average reconstructed pulpo pulse shape, as shown in figure \ref{fig:pulpo_shape_low}, and the 
reconstructed noise autocorrelation matrix from a pedestal run with random triggers, the digital filter 
weights are computed. Figures \ref{fig:w_time_MC_input_TDAS} and \ref{fig:w_amp_MC_input_TDAS} show the 
parameterization of the amplitude and timing weights for the MC pulse shape as a function of the ...

\par
\ldots {\textit{\bf MISSING END OF SENTENCE }} \ldots
\par

\begin{figure}[h!]
\begin{center}
\includegraphics[totalheight=7cm]{w_time_MC_input_TDAS.eps}
\end{center}
\caption[Time weights.]{Time weights $w_{\mathrm{time}}(t_0) \ldots w_{\mathrm{time}}(t_5)$ for a window size of 6 FADC slices for the pulse shape used in the MC simulations. The first weight $w_{\mathrm{time}}(t_0)$ is plotted as a function of the relative time $t_{\text{rel}}$ the trigger and the FADC clock in the range $[-0.5;0.5[ \ T_{\text{ADC}}$, the second weight in the range $[0.5;1.5[ \ T_{\text{ADC}}$ and so on. A binning resolution of $0.1  T_{\text{ADC}}$ has been chosen.} \label{fig:w_time_MC_input_TDAS}
\end{figure}

\begin{figure}[h!]
\begin{center}
\includegraphics[totalheight=7cm]{w_amp_MC_input_TDAS.eps}
\end{center}
\caption[Amplitude weights.]{Amplitude weights $w_{\mathrm{amp}}(t_0) \ldots w_{\mathrm{amp}}(t_5)$ for a window size of 6 FADC slices for the pulse shape used in the MC simulations. The first weight $w_{\mathrm{amp}}(t_0)$ is plotted as a function of the relative time $t_{\text{rel}}$ the trigger and the FADC clock in the range $[-0.5;0.5[ \ T_{\text{ADC}}$, the second weight in the range $[0.5;1.5[ \ T_{\text{ADC}}$ and so on. A binning resolution of $0.1  T_{\text{ADC}}$ has been chosen.} \label{fig:w_amp_MC_input_TDAS}
\end{figure}



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 using a window of $n$ slices:

\begin{equation}
e_{i_0}=\sum_{i=i_0}^{i_0+n-1} w_{\mathrm{amp}}(t_i)y(t_{i+i_0}) \qquad (e\tau)_{i_0}=\sum_{i=i_0}^{i_0+n-1} w_{\mathrm{time}}(t_i)y(t_{i+i_0}) 
\end{equation}

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

\begin{equation}
\tau=\frac{(e\tau)_{i_0^*}}{e_{i_0^*}}
\end{equation}

and the weights iterated:

\begin{equation}
E=\sum_{i=i_0^*}^{i_0^*+n-1} w_{\mathrm{amp}}(t_i - \tau)y(t_{i+i_0^*}) \qquad E \theta=\sum_{i=i_0^*}^{i_0^*+n-1} w_{\mathrm{time}}(t_i - \tau)y(t_{i+i_0^*}) \ .
\end{equation}

The reconstructed signal is then taken to be $E$ and the reconstructed arrival time $t_{\text{arrival}}$ is 

\begin{equation}
t_{\text{arrival}} = i_0^* + \tau + \theta \ .
\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}


\begin{figure}[h!]
\begin{center}
\includegraphics[totalheight=7cm]{amp_sliding.eps}
\includegraphics[totalheight=7cm]{time_sliding.eps}
\end{center}
\caption[Digital filter weights applied.]{Digital filter weights applied to the recorded FADC time slices of 
one calibration pulse. The left plot shows the result of the applied amplitude weights 
$e(t_0)=\sum_{i=0}^{i=n-1} w_{\mathrm{amp}}(t_0+i \cdot T_{\text{ADC}})y(t_0+i \cdot T_{\text{ADC}})$ and 
the right plot shows the result of the applied timing weights 
$e\tau(t_0)=\sum_{i=0}^{i=n-1} w_{\mathrm{time}}(t_0+i \cdot T_{\text{ADC}})y(t_0+i \cdot T_{\text{ADC}})$ .} 
\label{fig:amp_sliding}
\end{figure}


Figure \ref{fig:shape_fit_TDAS} shows the FADC slices of a single MC event together with the result of a full 
fit of the input MC pulse shape to the simulated FADC samples together with the result of the numerical fit 
using the digital filter. 


\begin{figure}[h!]
\begin{center}
\includegraphics[totalheight=7cm]{shape_fit_TDAS.eps}
\end{center}
\caption[Shape fit.]{Full fit to the MC pulse shape with the MC input shape and a numerical fit using the 
digital filter.} \label{fig:shape_fit_TDAS}
\end{figure}


\ldots {\it Hendrik ... }

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

\begin{description}
\item[Weights File:\xspace] An ascii-file containing the weights, the binning resolution and 
the window size. Currently, the following weight files have been created:
\begin{itemize}
\item "cosmics\_weights.dat'' with a window size of 6 FADC slices
\item "cosmics\_weights4.dat'' with a window size of 4 FADC slices
\item "calibration\_weights\_blue.dat'' with a window size of 6 FADC slices
\item "calibration\_weights\_UV.dat'' with a window size of 6 FADC slices and in the low-gain the 
calibration weigths obtained from blue pulses\footnote{UV-pulses saturating the high-gain are not yet
available.}.
\item "cosmics\_weights\_logaintest.dat'' with a window size of 6 FADC slices and swapped high-gain and low-gain
weights. This file is only used for stability tests.
\item "cosmics\_weights4\_logaintest.dat'' with a window size of 4 FADC slices and swapped high-gain and low-gain
weights. This file is only used for stability tests.
\end{itemize}
\end{description}

\begin{figure}[htp]
  \includegraphics[width=0.49\linewidth]{MExtractTimeAndChargeDigitalFilter_5Led_UV.eps}
  \includegraphics[width=0.49\linewidth]{MExtractTimeAndChargeDigitalFilter_23Led_Blue.eps}
\caption[Sketch calculated arrival times MExtractTimeAndChargeDigitalFilter]{%
Sketch of the calculated arrival times for the extractor {\textit{MExtractTimeAndChargeDigitalFilter}} 
for two typical calibration pulses (pedestals have been subtracted) and a typical inner pixel. 
The extraction window sizes modify the position of the (amplitude-weighted) mean FADC-slices slightly.
The pulse would be shifted half a slice to the right for an outer pixels. }
\label{fig:dfsketch}
\end{figure}

\subsubsection{Real Fit to the Expected Pulse Shape }

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

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

\begin{figure}[h!]
\begin{center}
\includegraphics[totalheight=7cm]{probability_fit_0ns.eps}
\end{center}
\caption[Fit Probability.]{Probability of the fit with the input signal shape to the simulated FADC samples 
including electronics and NSB noise.} \label{fig:w_amp_MC_input_TDAS.eps}
\end{figure}



\subsection{Used Extractors for this Analysis}

We tested in this TDAS the following parameterized extractors:

\begin{description}
\item[MExtractFixedWindow]: with the following intialization, if {\textit{maxbin}} defines the 
   mean position of the high-gain FADC slice which carries the pulse maximum \footnote{The function 
{\textit{MExtractor::SetRange(higain first, higain last, logain first, logain last)}} sets the extraction 
range with the high gain start bin {\textit{higain first}} to (including) the last bin {\textit{higain last}}. 
Analoguously for the low gain extraction range. Note that in MARS, the low-gain FADC samples start with 
the index 0 again, thus {\textit{maxbin+0.5}} means in reality {\textit{maxbin+15+0.5}}. }
:
\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}}-2,{\textit{maxbin}}+5,{\textit{maxbin}}-0.5,{\textit{maxbin}}+6.5);
\item SetRange({\textit{maxbin}}-3,{\textit{maxbin}}+10,{\textit{maxbin}}-1.5,{\textit{maxbin}}+7.5);
\suspend{enumerate}
\item[MExtractFixedWindowSpline]: with the following initialization, if {\textit{maxbin}} defines the 
   mean position of the high-gain FADC slice carrying the pulse maximum \footnote{The function 
{\textit{MExtractor::SetRange(higain first, higain last, logain first, logain last)}} sets the extraction 
range with the high gain start bin {\textit{higain first}} to (including) the last bin {\textit{higain last}}. 
Analoguously for the low gain extraction range. Note that in MARS, the low-gain FADC samples start with 
the index 0 again, thus {\textit{maxbin+0.5}} means in reality {\textit{maxbin+15+0.5}}.}:
\resume{enumerate}
\item SetRange({\textit{maxbin}}-1,{\textit{maxbin}}+3,{\textit{maxbin}}+0.5,{\textit{maxbin}}+4.5);
\item SetRange({\textit{maxbin}}-1,{\textit{maxbin}}+3,{\textit{maxbin}}-0.5,{\textit{maxbin}}+5.5);
\item SetRange({\textit{maxbin}}-2,{\textit{maxbin}}+4,{\textit{maxbin}}-0.5,{\textit{maxbin}}+5.5);
\item SetRange({\textit{maxbin}}-2,{\textit{maxbin}}+6,{\textit{maxbin}}-0.5,{\textit{maxbin}}+7.5);
\item SetRange({\textit{maxbin}}-3,{\textit{maxbin}}+11,{\textit{maxbin}}-1.5,{\textit{maxbin}}+8.5);
\suspend{enumerate}
\item[MExtractFixedWindowPeakSearch]: with the following initialization: \\
SetRange(0,18,2,14); and:
\resume{enumerate}
\item SetWindows(2,2,2); SetOffsetFromWindow(0);
\item SetWindows(4,4,2); SetOffsetFromWindow(1);
\item SetWindows(4,6,4); SetOffsetFromWindow(0);
\item SetWindows(6,6,4); SetOffsetFromWindow(1);
\item SetWindows(8,8,4); SetOffsetFromWindow(1);
\item SetWindows(14,10,4); SetOffsetFromWindow(2);
\suspend{enumerate}
\item[MExtractTimeAndChargeSlidingWindow]: with the following initialization: \\
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 initialization:
\resume{enumerate}
\item SetChargeType(MExtractTimeAndChargeSpline::kAmplitude); \\
SetRange(0,10,4,11); 
\suspend{enumerate}
SetChargeType(MExtractTimeAndChargeSpline::kIntegral); \\
SetRange(0,18,2,14); \\
and:
\resume{enumerate}
\item SetRiseTime(0.5); SetFallTime(0.5);
\item SetRiseTime(0.5); SetFallTime(1.5);
\item SetRiseTime(1.0); SetFallTime(3.0);
\item SetRiseTime(1.5); SetFallTime(4.5);
\suspend{enumerate}
\item[MExtractTimeAndChargeDigitalFilter]: with the following initialization:
\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'');
\suspend{enumerate}
\item[``Real Fit'']: (not yet implemented, one try)
\resume{enumerate}
\item Real Fit
\end{enumerate}
\end{description}

Note that the extractors \#30, \#31 are used only to test the stability of the extraction against 
changes in the pulse-shape.

References: \cite{OF77,OF94}.


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