<reponame>fermi-lat/CalRecon
\documentclass[a4paper,11pt]{article}

\title{Calorimeter moments analysis}
\author{<NAME> (<EMAIL>),
  <NAME> (<EMAIL>)}


\usepackage{ifthen}
\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{hyperref}
\usepackage{bbm}
\usepackage[margin=2cm,bottom=3cm,top=4cm,marginparsep=0pt,marginparwidth=0pt]%
	   {geometry}
\usepackage[margin=1cm, font=small, labelfont=bf, labelsep=endash]{caption}

\newcommand{\pder}[2]{\frac{\partial#1}{\partial#2}}
\newcommand{\pdersec}[3]%
           {\ifthenelse{\equal{#2}{#3}}%
             {\frac{\partial^2#1}{\partial#2^2}}%
             {\frac{\partial^2#1}{\partial#2\partial#3}}%
           }
\newcommand{\itm}{\mathbbm I}
\newcommand{\itc}[1]{\itm_{#1}}
\newcommand{\firstder}[2]{\frac{{\rm d}#1}{{\rm d}#2}}
\newcommand{\secder}[2]{\frac{{\rm d}^2#1}{{\rm d}#2^2}}
\newcommand{\tmax}{t_{\rm max}}
\newcommand{\diff}{{\rm d}}
\newcommand{\xdir}{\ensuremath{x_{\rm dir}}}
\newcommand{\ydir}{\ensuremath{y_{\rm dir}}}
\newcommand{\zdir}{\ensuremath{z_{\rm dir}}}

\begin{document}

\maketitle

\abstract{These are some (more or less) random notes about the calorimeter
  moments analysis. The first part is a brief description of the basic
  code that I put together while trying to understand it.
  There's also some stuff dealing with possible improvements of the
  moments analysis, namely the measurement of the shower skeweness and the
  the error analysis aimed at the projection of the calorimeter clusters
  into the ACD.
}


\section{Basics of the CAL moments analysis}

The code for the moments analysis basically calculates the moment of inertia
tensor (using energy as the weight instead of mass) and then diagonalizes this
to get the three principal axes. The basic definitions can be found in
\cite{goldstein,landau} and in our case they read:
\begin{align}
  \itc{xx}  = \sum_{i=1}^n w_i(r_i^2 - x_i^2),\quad
  \itc{yy} &= \sum_{i=1}^n w_i(r_i^2 - y_i^2),\quad
  \itc{zz}  = \sum_{i=1}^n w_i(r_i^2 - z_i^2)\\
  \itc{xy}  = -\sum_{i=1}^n w_ix_iy_i,\quad
  \itc{xz} &= -\sum_{i=1}^n w_ix_iz_i,\quad
  \itc{yz}  = -\sum_{i=1}^n w_iy_iz_i
\end{align}
where the index $i$ runs over the $n$ hits in the calorimeter and the $w_i$ are
the weights associated with the hits (essentially the energy release).
In addition to the moment of inertia, the sum of weights and the coordinates
of the energy centroids are also used:
\begin{align}
  W &= \sum_{i=1}^n w_i\\
  \mathbf{r}_c &= \frac{\sum_{i=1}^n w_i\mathbf{r}_i}{W}
\end{align}

In order to reduce the tensor of inertia to the principal axes we do have to
solve the secular equation:
\begin{equation}
  \det(\itm - \lambda{\mathbbm 1}) =
  \det\begin{pmatrix}
  \itc{xx} - \lambda & \itc{xy} & \itc{xz}\\
  \itc{xy} & \itc{yy} - \lambda & \itc{yz}\\
  \itc{xz} & \itc{yz} & \itc{zz} - \lambda
  \end{pmatrix} = 0
\end{equation}
which is a cubic equation  in $\lambda$ yielding the three eigenvalues.
By working out the tedious algebra we can write the equation as:
$$
\lambda^3 + c_2\lambda^2 + c_1\lambda + c_0 = 0
$$
where:
\begin{align}
  c_2 &= -(\itc{xx} + \itc{yy} + \itc{zz})\\
  c_1 &= \itc{xx}\itc{yy} + \itc{yy}\itc{zz} + \itc{xx}\itc{zz} -
  (\itc{xy}^2 + \itc{yz}^2 + \itc{xz}^2)\\
  c_0 &= -\itc{xx}\itc{yy}\itc{zz} - 2\itc{xy}\itc{yz}\itc{xz} +
  \itc{xx}\itc{yz}^2 + \itc{yy}\itc{xz}^2 + \itc{zz}\itc{xy}^2
\end{align}
If we now define a new variable $\lambda' = \lambda + c_2/3$, the previous
equation becomes:
$$
\lambda'^3 + a\lambda' + b = 0
$$
where:
\begin{align}
  a &= \left(\frac{3c_1 - c_2^2}{3}\right)\\
  b &= \left(\frac{27c_0 + 2c_2^2 - 9c_1c_2}{27}\right)
\end{align}
(the algebra is fully worked out in \cite{wolfram}).
We now set:
\begin{align}
  m &= 2\sqrt{\frac{-a}{3}}\\
  \psi &= \frac{1}{3}\arccos\left(\frac{3b}{am}\right)
\end{align}
and, finally we get the three real solutions (guranteed by the fact that the
tensor of inertia is symmetruc):
\begin{align}
  \lambda_0 &= m\cos(\psi) - c_2/3\\
  \lambda_1 &= m\cos(\psi + 2\pi/3) - c_2/3\\
  \lambda_2 &= m\cos(\psi + 4\pi/3) - c_2/3
\end{align}
Once we have the three eigenvalues we can work out the calcultaion of the
eigenvectors $\mathbf{e}^i$ ($i = 1\ldots3$) defined by:
$$
\itm\mathbf{e}^i = \lambda_i\mathbf{e}^i
$$
Following these conventions, $\lambda_1$ is the largest eigenvalue and,
as a consequence, $\mathbf{e}^1$ is the principal axis of the cluster.

Once the three principal axis of the cluster have been found, the
cluster $\chi^2$ (normalized to the number of \emph{degree of freedom}) is
calculated as:
\begin{equation}
  \chi^2 = \frac{\sum_{i=1}^n w_i d_i^2}{nW}
\end{equation}
where $d_i$ are the distances from each of the calorimeter hits to the
axis parallel to $\mathbf{e}^1$ and passing throught the cluster
centroid. Finally the some well know Merit quantities are calculated:
\begin{align}
  \texttt{CalTransRms} &= \sqrt{\frac{|\lambda_1|}{W}}\\
  \texttt{CalLongRms}  &= \sqrt{\frac{|\lambda_0| + |\lambda_2|}{2W\log L}}\\
  \texttt{CalLRmsAsym} &= \sqrt{\frac{|\lambda_0| - |\lambda_2|}
    {|\lambda_0| + |\lambda_2|}}
\end{align}
where $L$ is the number of radiation lengths transversed.

\subsection{Outline of the iterative moments analysis}

Put here some details about the iteration scheme, as they are
relevant for the calculation of the skewness (i.e. you get significantly
different answers at different steps).



\section{The (toy) problem in 2D}

I thought the problem of the diagonalization of the inertia tensor in 2D
could be useful for working out the error analysis, so here is a quick look at
it. The basic definitions read as:
\begin{align}
  \itc{xx}  = \sum_{i=1}^n w_i y_i^2,\quad
  \itc{yy}  = \sum_{i=1}^n w_i x_i^2,\quad
  \itc{xy}  = -\sum_{i=1}^n w_ix_iy_i
\end{align}
and the secular equation is:
\begin{equation}
  \det(\itm - \lambda{\mathbbm 1}) =
  \det\begin{pmatrix}
  \itc{xx} - \lambda & \itc{xy}\\
  \itc{xy} & \itc{yy} - \lambda
  \end{pmatrix} = 0
\end{equation}
The eigenvalues are readily found:
\begin{align}
  \lambda_0 &= \frac{(\itc{xx} + \itc{yy}) -
    \sqrt{(\itc{xx} - \itc{yy})^2 + 4\itc{xy}^2}}{2}\\
  \lambda_1 &= \frac{(\itc{xx} + \itc{yy}) +
    \sqrt{(\itc{xx} - \itc{yy})^2 + 4\itc{xy}^2}}{2}\\
\end{align}

At this point, assuming that $\lambda_0$ is the smallest eigenvalue, we're
left with the problem of calculating the corresponding eigenvector
$\mathbf{e}^0$, obeying the equation:
$$
\itm \mathbf{e}^0 = \lambda_0\mathbf{e}^0
$$
Being the problem bi-dimensional, the two eigenvectors can be parametrized in
terms of a single rotation angle $\phi$ (to be determined), representing the
angle of the principal axis with respect to the original reference frame i.e.:
\begin{equation}
  \begin{pmatrix}
    e^0_x\\
    e^0_y
  \end{pmatrix} = 
  \begin{pmatrix}
    \cos\phi\\
    \sin\phi
  \end{pmatrix},\quad
  \begin{pmatrix}
    e^1_x\\
    e^1_y
  \end{pmatrix} = 
  \begin{pmatrix}
    -\sin\phi\\
    \cos\phi
  \end{pmatrix}
\end{equation}
By putting everything together, from the definition of the principal axis
we get:
$$
\begin{pmatrix}
  \itc{xx} & \itc{xy}\\
  \itc{xy} & \itc{yy}
\end{pmatrix} 
\begin{pmatrix}
  \cos\phi\\
  \sin\phi
\end{pmatrix} = \lambda_1
\begin{pmatrix}
  \cos\phi\\
  \sin\phi
\end{pmatrix}
$$
The first of the two equations in the system can be squared, yielding:
$$
\frac{(\itc{xx} - \itc{yy})}{\itc{xy}}\tan\phi = \tan^2\phi - 1 
$$
and eventually, through the trigonometric equation
$$
\tan(2\phi) = \frac{2\tan\phi}{1 - \tan^2\phi}
$$
we get:
\begin{equation}
  \phi = -\frac{1}{2}
  \arctan \left( \frac{2\itc{xy}}{\itc{yy} - \itc{xx}} \right)
\end{equation}
The rotation matrix between the original system and the one defined by the
principal axis has the (transpose of) the two eigenvectors as its rows:
\begin{equation}
  S =
  \begin{pmatrix}
    e^0_x & e^0_y\\
    e^1_x & e^1_y
  \end{pmatrix} = 
  \begin{pmatrix}
    \cos\phi & \sin\phi\\
    -\sin\phi & \cos\phi
  \end{pmatrix}
\end{equation}
and obviously we have:
\begin{equation}\label{eq:rotation}
  \lambda = 
  \begin{pmatrix}
    \lambda_0 & 0\\
    0 & \lambda_1
  \end{pmatrix} = S\;\itm\;S^{-1} = S\;\itm\;S^{\rm T}
\end{equation}


\subsection{Error analysis in 2 dimensions}

As a first ingredient we'll need the derivatives of the components of the
tensor of inertia with respect to the coordinates and the weights:
\begin{align}
\pder{\itc{xx}}{x_i} &= 0,\quad
\pder{\itc{xx}}{y_i}  = 2w_iy_i,\quad
\pder{\itc{xx}}{w_i}  = y_i^2\\
\pder{\itc{yy}}{x_i} &= 2w_ix_i,\quad
\pder{\itc{yy}}{y_i}  = 0,\quad
\pder{\itc{yy}}{w_i}  = x_i^2\\
\pder{\itc{xy}}{x_i} &= -w_iy_i,\quad
\pder{\itc{xy}}{y_i}  = -w_ix_i,\quad
\pder{\itc{xy}}{w_i}  = -x_iy_i
\end{align}
We can then calculate the full covariance matrix of the errors using the
usual formula (see \cite{pdg}, section 32.1.4 for instance). Assuming that
the errors on the two positions and on the weights are mutually not correlated
(i.e. their covariance matrix is diagonal), we have:
\begin{equation}
  \Sigma_{k-l} = \sum_{i=1}^n 
  \pder{\itc{k}}{x_i}\pder{\itc{l}}{x_i}(\Delta x_i)^2 +
  \pder{\itc{k}}{y_i}\pder{\itc{l}}{y_i}(\Delta y_i)^2 +
  \pder{\itc{k}}{w_i}\pder{\itc{l}}{w_i}(\Delta w_i)^2
\end{equation}
where $k$ and $l$ run over the three (double) indexes $xx$, $yy$ and $xy$.
We're ready to work out the details:
\begin{align}
  \Sigma_{xx-xx} &= (\Delta\itc{xx})^2 =
  \sum_{i=1}^n \left[4w^2_iy^2_i(\Delta y_i)^2 + y_i^4(\Delta w_i)^2\right]\\
  \Sigma_{yy-yy} &= (\Delta\itc{yy})^2 =
  \sum_{i=1}^n \left[4w^2_ix^2_i(\Delta x_i)^2 + x_i^4(\Delta w_i)^2\right]\\
  \Sigma_{xy-xy} &= (\Delta\itc{xy})^2 =
  \sum_{i=1}^n \left[w_i^2y_i^2(\Delta x_i)^2 + w_i^2x_i^2(\Delta y_i)^2 +
    x_i^2y_i^2(\Delta w_i)^2 \right]\\
  \Sigma_{xx-yy} &= 
  \sum_{i=1}^n \left[ x_i^2y_i^2 (\Delta w_i)^2\right]\\
  \Sigma_{xx-xy} &= 
  -\sum_{i=1}^n \left[ 2w_i^2x_iy_i(\Delta y_i)^2 +
    x_iy_i^3(\Delta w_i)^2 \right]\\
  \Sigma_{xy-yy} &= 
  -\sum_{i=1}^n \left[ 2w_i^2x_iy_i(\Delta x_i)^2 +
    x_i^3y_i(\Delta w_i)^2 \right]
\end{align}

The rest of this section follows closely the prescription described in
\cite{errors} for the error propagation. We can slice the $2 \times 2$ tensor
of inertia and define a 4-component vector with the two columns one on top 
of the other (this is what we call the $vec$ operator):
\begin{equation}
  vec(\itm) = 
  \begin{pmatrix}
    \itc{xx}\\
    \itc{xy}\\
    \itc{xy}\\
    \itc{yy}
  \end{pmatrix}
\end{equation}
That said, we can rewite the equation (\ref{eq:rotation}) using the
Kronecker product of the rotation matrix
$$
T = S \otimes S
$$
as:
\begin{equation}\label{eq:tensor_transform}
  vec(\lambda) = 
  \begin{pmatrix}
    \lambda_0\\
    0\\
    0\\
    \lambda_1
  \end{pmatrix} = T vec(\itm)
\end{equation}
It is useful to rearrange the the elements of the $vec$ operator in such a way
that the diagonal elements of the tensor come first, followed by the
non-diagonal ones, getting rid of the duplicated terms.
This is accomplished by introducing the matrix:
\begin{equation}
  D = 
  \begin{pmatrix}
    1 & 0 & 0 & 0\\
    0 & 0 & 0 & 1\\
    0 & \frac{1}{2} & \frac{1}{2} & 0
  \end{pmatrix}
\end{equation}
which allows to define a new operator $v_d$:
\begin{equation}
  v_d(\itm) = D vec(\itm) =
  \begin{pmatrix}
    1 & 0 & 0 & 0\\
    0 & 0 & 0 & 1\\
    0 & \frac{1}{2} & \frac{1}{2} & 0
  \end{pmatrix}
  \begin{pmatrix}
    \itc{xx}\\
    \itc{xy}\\
    \itc{xy}\\
    \itc{yy}
  \end{pmatrix} = 
  \begin{pmatrix}
    \itc{xx}\\
    \itc{yy}\\
    \itc{xy}
  \end{pmatrix}
\end{equation}
Talking about which, the covariance matrix of $v_d(\itm)$ reads:
\begin{equation}
  \Sigma_{v_d(\itm)} = 
  \begin{pmatrix}
    \Sigma_{xx-xx} & \Sigma_{xx-yy} & \Sigma_{xx-xy}\\
    \Sigma_{xx-yy} & \Sigma_{yy-yy} & \Sigma_{xy-yy}\\
    \Sigma_{xx-xy} & \Sigma_{xy-yy} & \Sigma_{xy-xy}
  \end{pmatrix}
\end{equation}
in terms of the quantities we have calculated a few lines above.
In order to go back from the $v_d$ to the $vec$ representation we need
the so called pseudo-inverse of $D$:
\begin{equation}
  D^+ = 
  \begin{pmatrix}
    1 & 0 & 0\\
    0 & 0 & 1\\
    0 & 0 & 1\\
    0 & 1 & 0
  \end{pmatrix}
\end{equation}
The paper \cite{errors} is wrong to this respect (see \cite{errors_corr}).
Using these definitions, equation (\ref{eq:tensor_transform}) can be rewritten
as:
\begin{equation}
  v_d(\lambda) =
  \begin{pmatrix}
    \lambda_0\\
    \lambda_1\\
    0
  \end{pmatrix} =
  D vec(\lambda) = DT vec(\itm) = DTD^+ v_d(\itm)
\end{equation}
If we define:
\begin{equation}
  V = DTD^+
\end{equation}
we can rewrite the previous equation as:
\begin{equation}
  v_d(\lambda) = V v_d(\itm)
\end{equation}

The infinitesimal change $\diff S$ in the rotation matrix $S$ when we change
the rotation angle $\phi$ by an infinitesimal amount $\diff\phi$ can be easily
calculated by differentiating $S$ itself:
$$
\diff S =
\begin{pmatrix}
  -\sin\phi\;\diff\phi & \cos\phi\;\diff\phi\\
  -\cos\phi\;\diff\phi & -\sin\phi\;\diff\phi\\
\end{pmatrix}
$$
If we introduce the antisimmetric tensor:
\begin{equation}
  \Omega = 
  \begin{pmatrix}
    0 & -\diff\phi\\
    \diff\phi & 0
  \end{pmatrix}
\end{equation}
we can rewrite the previous equation as:
\begin{equation}
  \diff S = -S\Omega
\end{equation}
and, along the same lines, we have:
\begin{equation}
  \diff S^{\rm T} = \Omega S^{\rm T}
\end{equation}
as can be verified by direct matrix multiplication. It seems a bit odd, here,
to call $\Omega$ an infinitesimal quantity (I would have rather named it
$\diff\Omega$), but we'll bravely follow the conventions used in the paper to
avoid confusion.
We now define the quantity:
\begin{equation}
  \Omega^p = S\Omega S^{\rm T} =
\begin{pmatrix}
  0 & -\diff\phi\\
  \diff\phi & 0\\
\end{pmatrix} = \Omega
\end{equation}
By putting all together we have:
\begin{equation}
  S \diff \itm S^{\rm T} = \diff\mu =
  \Omega^p\lambda - \lambda\Omega^p + \diff\lambda = 
  \begin{pmatrix}
    \diff\lambda_0 & (\lambda_0 - \lambda_1)\diff\phi\\
    (\lambda_0 - \lambda_1)\diff\phi & \diff\lambda_1
  \end{pmatrix}
\end{equation}
and we can write:
\begin{equation}
  vec(\diff\mu) = G\beta
\end{equation}
where:
\begin{equation}
  G =
  \begin{pmatrix}
    1 & 0 & 0\\
    0 & 0 & (\lambda_1 - \lambda_0)\\
    0 & 0 & (\lambda_1 - \lambda_0)\\
    0 & 1 & 0
  \end{pmatrix}
\end{equation}
and
\begin{equation}
  \beta =
  \begin{pmatrix}
    \diff\lambda_0\\
    \diff\lambda_1\\
    -\diff\phi
  \end{pmatrix}
\end{equation}
(again, $\beta$ is a differential quantities, so it is a bit odd to name
it without a $\diff$ in front). Further on through the paper:
\begin{equation}
  v_d(\diff\itm) =
  \begin{pmatrix}
    d\itc{xx}\\
    d\itc{yy}\\
    d\itc{xy}
  \end{pmatrix} =
  D (S^{\rm T} \otimes S^{\rm T}) vec(d\mu) =
  D (S^{\rm T} \otimes S^{\rm T}) G\beta = F\beta
\end{equation}
where we have defined:
\begin{equation}
  F = D (S^{\rm T} \otimes S^{\rm T}) G
\end{equation}
All we have to do is invert this equation, namely find $F^{-1}$. If we were
dealing with square matrices, all we had to do would be:
$$
F^{-1} = \left( D (S^{\rm T} \otimes S^{\rm T}) G \right)^{-1} =
G^{-1} \left( S^{\rm T} \otimes S^{\rm T} \right)^{-1} D^{-1} = 
G^{-1} \left( S \otimes S \right) D^{-1} = G^{-1} T D^{-1}
$$
But in fact $G$ and $D$ are rectangular matrices, so we need again the pseudo
inverses. We have already the one for $D$, while for $G$, since
$G^TG$ is not singular, the solution is even easier:
\begin{equation}
  G^+ = \left( G^{\rm T}G \right)^{-1} G^{\rm T} =
  \begin{pmatrix}
    1 & 0 & 0 & 0\\
    0 & 0 & 0 & 1\\
    0 & \frac{1}{2(\lambda_1-\lambda_0)} & \frac{1}{2(\lambda_1-\lambda_0)} & 0
  \end{pmatrix}
\end{equation}
(this satisfies $G^+G = I_{3\times3}$, as can be easily verified by direct
multiplication). Summarizing:
\begin{equation}
  \beta = F^{-1}v_d(\diff\itm)
\end{equation}
and, more important, we get the full covariance matrix of the eigenvalues
and rotation angle:
\begin{equation}
  \Sigma_\beta = F^{-1} \Sigma_{v_d(\itm)} (F^{-1})^T
\end{equation}
Once we have the error on the rotation angle:
\begin{equation}
  (\Delta\phi)^2 = (\Sigma_\beta)_{33}
\end{equation}
the covariance matrix of the two components of the principal axis is:
\begin{equation}
  \Sigma_{\mathbf{e}^0} =
  \begin{pmatrix}
    \pder{e^0_x}{\phi}\pder{e^0_x}{\phi}(\Delta\phi)^2 &
    \pder{e^0_x}{\phi}\pder{e^0_y}{\phi}(\Delta\phi)^2\\
    \pder{e^0_x}{\phi}\pder{e^0_y}{\phi}(\Delta\phi)^2 &
    \pder{e^0_y}{\phi}\pder{e^0_y}{\phi}(\Delta\phi)^2
  \end{pmatrix} =
  \begin{pmatrix}
    \sin^2\phi & -\sin\phi\cos\phi\\
    -\sin\phi\cos\phi & \cos^2\phi
  \end{pmatrix}(\Delta\phi)^2
\end{equation}


\section{Error analysis in 3 dimensions}

Along the same lines of the previous subsection the derivatives of the
components of the inertia tensor are:
\begin{align}
\pder{\itc{xx}}{x_i} &= 0,\quad
\pder{\itc{xx}}{y_i}  = 2w_iy_i,\quad
\pder{\itc{xx}}{z_i}  = 2w_iz_i,\quad
\pder{\itc{xx}}{w_i}  = (y_i^2 + z_i^2)\\
\pder{\itc{yy}}{x_i} &= 2w_ix_i,\quad
\pder{\itc{yy}}{y_i}  = 0,\quad
\pder{\itc{yy}}{z_i}  = 2w_iz_i,\quad
\pder{\itc{yy}}{w_i}  = (x_i^2 + z_i^2)\\
\pder{\itc{zz}}{x_i} &= 2w_ix_i,\quad
\pder{\itc{zz}}{y_i}  = 2w_iy_i,\quad
\pder{\itc{zz}}{z_i}  = 0,\quad
\pder{\itc{zz}}{w_i}  = (x_i^2 + y_i^2)\\
\pder{\itc{xy}}{x_i} &= -w_iy_i,\quad
\pder{\itc{xy}}{y_i}  = -w_ix_i,\quad
\pder{\itc{xy}}{z_i}  = 0,\quad
\pder{\itc{xy}}{w_i}  = -x_iy_i\\
\pder{\itc{xz}}{x_i} &= -w_iz_i,\quad
\pder{\itc{xz}}{y_i}  = 0,\quad
\pder{\itc{xz}}{z_i}  = -w_ix_i,\quad
\pder{\itc{xz}}{w_i}  = -x_iz_i\\
\pder{\itc{yz}}{x_i} &= 0,\quad
\pder{\itc{yz}}{y_i}  = -w_iz_i,\quad
\pder{\itc{yz}}{z_i}  = -w_iy_i,\quad
\pder{\itc{yz}}{w_i}  = -y_iz_i
\end{align}
The elements of the covariance matrix are:
\begin{equation}
  \Sigma_{k-l} = \sum_{i=1}^n 
  \pder{\itc{k}}{x_i}\pder{\itc{l}}{x_i}(\Delta x_i)^2 +
  \pder{\itc{k}}{y_i}\pder{\itc{l}}{y_i}(\Delta y_i)^2 +
  \pder{\itc{k}}{z_i}\pder{\itc{l}}{z_i}(\Delta z_i)^2 +
  \pder{\itc{k}}{w_i}\pder{\itc{l}}{w_i}(\Delta w_i)^2
\end{equation}
where now $k$ and $l$ run over the 6 independent (double) indexes of the
indertia tensor $xx$, $yy$, $zz$, $xy$, $xz$ and $yz$. Therefore the
$6\times6$ covariancewe will have $6(6+1)/2 = 21$ components:
\begin{align}
  \Sigma_{xx-xx} &= \sum_{i=1}^n \left\{
  4w_i^2 \left[ y_i^2(\Delta y_i)^2 + z_i^2(\Delta z_i)^2 \right] +
  \left(y_i^2 + z_i^2\right)^2(\Delta w_i)^2
  \right\}\\
  \Sigma_{xx-yy} &= \sum_{i=1}^n \left\{
  4w_i^2z_i^2 (\Delta z_i)^2 + (x_i^2 + z_i^2)(y_i^2 + z_i^2)(\Delta w_i)^2
  \right\}\\
  \Sigma_{xx-zz} &= \sum_{i=1}^n \left\{
  4w_i^2y_i^2 (\Delta y_i)^2 + (x_i^2 + y_i^2)(y_i^2 + z_i^2)(\Delta w_i)^2
  \right\}\\
  \Sigma_{xx-xy} &= -\sum_{i=1}^n x_iy_i  \left\{
  2w_i^2 (\Delta y_i)^2 + (y_i^2 + z_i^2)(\Delta w_i)^2
  \right\}\\
  \Sigma_{xx-xz} &= -\sum_{i=1}^n x_iz_i \left\{
  2w_i^2 (\Delta z_i)^2 + (y_i^2 + z_i^2)(\Delta w_i)^2
  \right\}\\
  \Sigma_{xx-yz} &= -\sum_{i=1}^n y_iz_i \left\{
  2w_i^2 \left[ (\Delta y_i)^2 + (\Delta z_i)^2 \right] +
  (y_i^2 + z_i^2)(\Delta w_i)^2
  \right\}\\
  \Sigma_{yy-yy} &= \sum_{i=1}^n \left\{
  4w_i^2 \left[ x_i^2(\Delta x_i)^2 + z_i^2(\Delta z_i)^2 \right] +
  \left(x_i^2 + z_i^2\right)^2(\Delta w_i)^2
  \right\}\\
  \Sigma_{yy-zz} &= \sum_{i=1}^n \left\{
  4w_i^2x_i^2 (\Delta x_i)^2 + (x_i^2 + y_i^2)(x_i^2 + z_i^2)(\Delta w_i)^2
  \right\}\\
  \Sigma_{yy-xy} &= -\sum_{i=1}^n x_iy_i \left\{
  2w_i^2 (\Delta x_i)^2 + (x_i^2 + z_i^2)(\Delta w_i)^2
  \right\}\\
  \Sigma_{yy-xz} &= -\sum_{i=1}^n x_iz_i \left\{
  2w_i^2 \left[ (\Delta x_i)^2 + (\Delta z_i)^2 \right] +
  (x_i^2 + z_i^2)(\Delta w_i)^2
  \right\}
\end{align}
\begin{align}
  \Sigma_{yy-yz} &= -\sum_{i=1}^n y_iz_i \left\{
  2w_i^2 (\Delta z_i)^2 + (x_i^2 + z_i^2)(\Delta w_i)^2
  \right\}\\
  \Sigma_{zz-zz} &= \sum_{i=1}^n \left\{
  4w_i^2 \left[ x_i^2(\Delta x_i)^2 + y_i^2(\Delta y_i)^2 \right] +
  \left(x_i^2 + y_i^2\right)^2(\Delta w_i)^2
  \right\}\\
  \Sigma_{zz-xy} &= -\sum_{i=1}^n x_iy_i \left\{
  2w_i^2 \left[ (\Delta x_i)^2 + (\Delta y_i)^2 \right] +
  (x_i^2 + y_i^2)(\Delta w_i)^2
  \right\}\\
  \Sigma_{zz-xz} &= -\sum_{i=1}^n x_iz_i \left\{
  2w_i^2 (\Delta x_i)^2 + (x_i^2 + y_i^2)(\Delta w_i)^2
  \right\}\\
  \Sigma_{zz-yz} &= -\sum_{i=1}^n y_iz_i \left\{
  2w_i^2 (\Delta y_i)^2 + (x_i^2 + y_i^2)(\Delta w_i)^2
  \right\}\\
  \Sigma_{xy-xy} &= \sum_{i=1}^n \left\{
  w_i^2 \left[ y_i^2 (\Delta x_i)^2 + x_i^2 (\Delta y_i)^2 \right] +
  x_i^2y_i^2 (\Delta w_i)^2
  \right\}\\
  \Sigma_{xy-xz} &= \sum_{i=1}^n y_iz_i \left\{
  w_i^2 (\Delta x_i)^2 + x_i^2 (\Delta w_i)^2
  \right\}\\
  \Sigma_{xy-yz} &= \sum_{i=1}^n x_iz_i \left\{
  w_i^2 (\Delta y_i)^2 + y_i^2 (\Delta w_i)^2
  \right\}\\
  \Sigma_{yz-yz} &= \sum_{i=1}^n \left\{
  w_i^2 \left[ z_i^2 (\Delta y_i)^2 + y_i^2 (\Delta z_i)^2 \right] +
  y_i^2z_i^2 (\Delta w_i)^2
  \right\}\\
  \Sigma_{yz-xz} &= \sum_{i=1}^n x_iy_i \left\{
  w_i^2 (\Delta z_i)^2 + z_i^2 (\Delta w_i)^2
  \right\}\\
  \Sigma_{xz-xz} &= \sum_{i=1}^n \left\{
  w_i^2 \left[ z_i^2 (\Delta x_i)^2 + x_i^2 (\Delta z_i)^2 \right] +
  x_i^2z_i^2 (\Delta w_i)^2
  \right\}
\end{align}
In analogy with the two-dimensional case we define:
\begin{equation}
  D = 
  \begin{pmatrix}
    1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\
    0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0\\
    0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1\\
    0 & \frac{1}{2} & 0 & \frac{1}{2} & 0 & 0 & 0 & 0 & 0\\
    0 & 0 & \frac{1}{2} & 0 & 0 & 0 & \frac{1}{2} & 0 & 0\\
    0 & 0 & 0 & 0 & 0 & \frac{1}{2} & 0 & \frac{1}{2} & 0
  \end{pmatrix}
\end{equation}
so that 
\begin{equation}
  v_d(\itm) = D vec(\itm) =
  \begin{pmatrix}
    \itc{xx}\\
    \itc{yy}\\
    \itc{zz}\\
    \itc{xy}\\
    \itc{xz}\\
    \itc{yz}
  \end{pmatrix}
\end{equation}
and the $6\times6$ covariance matrix of the components of the inertia tensor
is:
\begin{equation}
  \Sigma_{v_d(\itm)} = 
  \begin{pmatrix}
    \Sigma_{xx-xx} & \Sigma_{xx-yy} & \Sigma_{xx-zz} &
    \Sigma_{xx-xy} & \Sigma_{xx-xz} & \Sigma_{xx-yz}\\
    \Sigma_{xx-yy} & \Sigma_{yy-yy} & \Sigma_{yy-zz} &
    \Sigma_{yy-xy} & \Sigma_{yy-xz} & \Sigma_{yy-yz}\\
    \Sigma_{xx-zz} & \Sigma_{yy-zz} & \Sigma_{zz-zz} &
    \Sigma_{zz-xy} & \Sigma_{zz-xz} & \Sigma_{zz-yz}\\
    \Sigma_{xx-xy} & \Sigma_{yy-xy} & \Sigma_{zz-xy} &
    \Sigma_{xy-xy} & \Sigma_{xy-xz} & \Sigma_{xy-yz}\\
    \Sigma_{xx-xz} & \Sigma_{yy-xz} & \Sigma_{zz-xz} &
    \Sigma_{xy-xz} & \Sigma_{xz-xz} & \Sigma_{xz-yz}\\
    \Sigma_{xx-yz} & \Sigma_{yy-yz} & \Sigma_{zz-yz} &
    \Sigma_{xy-yz} & \Sigma_{xz-yz} & \Sigma_{yz-yz}\\
  \end{pmatrix}
\end{equation}
In the three-dimensional case $D^+$ reads:
\begin{equation}
  D^+ =
  \begin{pmatrix}
    1 & 0 & 0 & 0 & 0 & 0\\
    0 & 0 & 0 & 1 & 0 & 0\\
    0 & 0 & 0 & 0 & 1 & 0\\
    0 & 0 & 0 & 1 & 0 & 0\\
    0 & 1 & 0 & 0 & 0 & 0\\
    0 & 0 & 0 & 0 & 0 & 1\\
    0 & 0 & 0 & 0 & 1 & 0\\
    0 & 0 & 0 & 0 & 0 & 1\\
    0 & 0 & 1 & 0 & 0 & 0
  \end{pmatrix}
\end{equation}
and the rotation matrix is obviously:
\begin{equation}
  S =
  \begin{pmatrix}
    e^0_x & e^0_y & e^0_z\\
    e^1_x & e^1_y & e^1_z\\
    e^2_x & e^2_y & e^2_z
  \end{pmatrix}
\end{equation}
The last pieces of information that we need are:
\begin{equation}
  S_p =
  \begin{pmatrix}
    0 & S_{20} & -S_{10}\\
    -S_{20} & 0 & S_{00}\\
    S_{10} & -S_{00} & 0\\
    0 & S_{21} & -S_{11}\\
    -S_{21} & 0 & S_{01}\\
    S_{11} & -S_{01} & 0\\
    0 & S_{22} & -S_{12}\\
    -S_{22} & 0 & S_{02}\\
    S_{12} & -S_{02} & 0
  \end{pmatrix}
\end{equation}
and:
\begin{equation}
  G^+ =
  \begin{pmatrix}
    1 & 0  & 0  & 0  & 0 & 0  & 0  & 0  & 0\\
    0 & 0  & 0  & 0  & 1 & 0  & 0  & 0  & 0\\
    0 & 0  & 0  & 0  & 0 & 0  & 0  & 0  & 1\\
    0 & 0  & 0  & 0  & 0 & \frac{1}{2(\lambda_2 - \lambda_1)} &
    0  & \frac{1}{2(\lambda_2 - \lambda_1)} & 0\\
    0 & 0  & \frac{1}{2(\lambda_0 - \lambda_2)} &
    0  & 0 & 0  & \frac{1}{2(\lambda_0 - \lambda_2)} & 0  & 0\\
    0 & \frac{1}{2(\lambda_1 - \lambda_0)} & 0 &
    \frac{1}{2(\lambda_1 - \lambda_0)} & 0 & 0  & 0  & 0  & 0
  \end{pmatrix}
\end{equation}
And at this point we just follow the paper (the new one, not the old one):
\begin{align}
  F^{-1} &= G^+ (S \otimes S) D^+\\
  K &=
  \begin{pmatrix}
    0_{9\times3} & S_p\\
    I_{3\times3} & 0_{3\times3}
  \end{pmatrix}
  F^{-1}
\end{align}
and eventually:
\begin{equation}
  \Sigma_{vec(S), \lambda} = K \Sigma_{v_d(\itm)} K^T
\end{equation}


\subsection{Simplified error analysis}

The main purpose in this section is to recover the expression for the
covariance matrix in the $(\xdir, \ydir, \zdir)$ representation when error
estimates on $\theta$ and $\phi$ are available---either on an event-by-event
basis or through some phenomenological parametrization. We will assume for
simplicity that the errors on $\theta$ and $\phi$ are uncorrelated.

The basic transformation equations are:
\begin{align}
  \xdir & = \sin\theta \cos\phi\\
  \ydir & = \sin\theta \sin\phi\\
  \zdir & = \cos\theta
\end{align}
and the elements of the covariance matrix, as a function of $\theta$
and $\phi$, read:
\begin{equation}
  \Sigma_{ij} = E[C_i C_j]
\end{equation}
where, denoting the central values with the subscript zero:
\begin{align}
  C_x & = (\sin\theta \cos\phi - \sin\theta_0 \cos\phi_0)\\
  C_y & = (\sin\theta \sin\phi - \sin\theta_0 \sin\phi_0)\\
  C_z & = (\cos\theta - \cos\theta_0)
\end{align}
We will take advantage of the fact that in the series expansion of the 
expectation values:
\begin{align}
  E[f(\theta, \phi)] & \approx f(\theta_0, \phi_0) + 
  \left[\pder{f(\theta, \phi)}{\theta}\right]_{\theta_0, \phi_0}
  \!\!\!\!\!\!\!\!\!\!E[(\theta - \theta_0)] +
  \left[\pder{f(\theta, \phi)}{\phi}\right]_{\theta_0, \phi_0}
  \!\!\!\!\!\!\!\!\!\!E[(\phi - \phi_0)] +\\\nonumber
  & \frac{1}{2}\left[\pdersec{f(\theta, \phi)}{\theta}{\theta}
  \right]_{\theta_0, \phi_0}
  \!\!\!\!\!\!\!\!\!\!E[(\theta - \theta_0)^2] +
  \frac{1}{2}\left[\pdersec{f(\theta, \phi)}{\phi}{\phi}
  \right]_{\theta_0, \phi_0}
  \!\!\!\!\!\!\!\!\!\!E[(\phi - \phi_0)^2] +\\\nonumber
  & \left[\pdersec{f(\theta, \phi)}{\theta}{\phi}\right]_{\theta_0, \phi_0}
  \!\!\!\!\!\!\!\!\!\!E[(\theta - \theta_0)(\phi - \phi_0)]
\end{align}
the first three terms vanish for obvious reasons, and so does the last,
if the errors on $\theta$ and $\phi$ are not correlated. Therefore:
\begin{align}
  E[f(\theta, \phi)] \approx
  \frac{1}{2}\left[\pdersec{f(\theta, \phi)}{\theta}{\theta}
  \right]_{\theta_0, \phi_0}
  \!\!\!\!\!\!\!\!\!\!\sigma^2_\theta + 
  \frac{1}{2}\left[\pdersec{f(\theta, \phi)}{\phi}{\phi}
    \right]_{\theta_0, \phi_0}
  \!\!\!\!\!\!\!\!\!\!\sigma^2_\phi
\end{align}
In the following we will need:
\begin{align}
  \nonumber
  \left[\pdersec{C_i(\theta, \phi)C_j(\theta, \phi)}{\theta}{\theta}
    \right]_{\theta_0, \phi_0}
  \quad {\rm and}\quad
  \left[\pdersec{C_i(\theta, \phi)C_j(\theta, \phi)}{\phi}{\phi}
    \right]_{\theta_0, \phi_0}
\end{align}
The basic blocks are:
\begin{align}
  \nonumber
  \pder{C_x}{\theta} &= \cos\theta\cos\phi ,\quad 
  \pder{C_x}{\phi} = -\sin\theta\sin\phi\\\nonumber
  \pder{C_y}{\theta} &= \cos\theta\sin\phi ,\quad
  \pder{C_y}{\phi} = \sin\theta\cos\phi\\\nonumber
  \pder{C_z}{\theta} &= -\sin\theta ,\quad
  \pder{C_z}{\phi} = 0\\\nonumber
  \pdersec{C_x}{\theta}{\theta} &= -\sin\theta\cos\phi ,\quad 
  \pdersec{C_x}{\phi}{\phi} = -\sin\theta\cos\phi\\\nonumber
  \pdersec{C_y}{\theta}{\theta} &= -\sin\theta\sin\phi ,\quad
  \pdersec{C_y}{\phi}{\phi} = -\sin\theta\sin\phi\\\nonumber
  \pdersec{C_z}{\theta}{\theta} &= -\cos\theta ,\quad
  \pdersec{C_z}{\phi}{\phi} = 0\\\nonumber
\end{align}
And the actual ingredients for the covariance matrix can be easily
calculated remembering that all the terms containing a $C_i$ multiplicative
factor are zero when calculated in $(\theta_0, \phi_0)$:
\begin{align}
  \left[\pdersec{C^2_x}{\theta}{\theta}\right]_{\theta_0, \phi_0}
  \!\!\!\!\!\!\!\! =
  \left[\pder{}{\theta}\left(2C_x\pder{C_x}{\theta}\right)
    \right]_{\theta_0, \phi_0}
  \!\!\!\!\!\!\!\! =
  \left[2\left(\pder{C_x}{\theta}\right)^2 +
    2C_x\pdersec{C_x}{\theta}{\theta}\right]_{\theta_0, \phi_0}
  \!\!\!\!\!\!\!\! = 2\cos^2\theta_0\cos^2\phi_0
\end{align}
So, along the same lines (and defining for compactness
$s_{\theta_0} = \sin\theta_0$, $c_{\theta_0} = \cos\theta_0$,
$s_{\phi_0} = \sin\phi_0$, $c_{\phi_0} = \cos\phi_0$):
\begin{equation}
  %m[0, 0] = ct2*cp2*dt2 + st2*sp2*dp2
  %m[1, 1] = ct2*sp2*dt2 + st2*cp2*dp2
  %m[2, 2] = st2*dt2
  %m[0, 1] = ct2*cp*sp*dt2 - st2*cp*sp*dp2
  %m[0, 2] = -st*ct*cp*dt2
  %m[1, 2] = -st*ct*sp*dt2
  %m[1, 0] = m[0, 1]
  %m[2, 0] = m[0, 2]
  %m[2, 1] = m[1, 2]
  \Sigma =
  \begin{pmatrix}
    c_{\theta_0}^2c_{\phi_0}^2\sigma^2_\theta +
    s_{\theta_0}^2s_{\phi_0}^2\sigma^2_\phi
    & c_{\theta}^2 c_{\phi}s_{\phi}\sigma^2_\theta -
    s_{\theta}^2c_{\phi}s_{\phi}\sigma^2_\phi
    & -s_{\theta}c_{\theta}c_{\phi} \sigma^2_\theta\\
    c_{\theta}^2 c_{\phi}s_{\phi}\sigma^2_\theta -
    s_{\theta}^2c_{\phi}s_{\phi}\sigma^2_\phi
    & c_{\theta_0}^2s_{\phi_0}^2\sigma^2_\theta +
    s_{\theta_0}^2c_{\phi_0}^2\sigma^2_\phi
    & -s_{\theta}c_{\theta}s_{\phi} \sigma^2_\theta\\
    -s_{\theta}c_{\theta}c_{\phi} \sigma^2_\theta
    & -s_{\theta}c_{\theta}s_{\phi} \sigma^2_\theta
    & s_{\theta_0}^2\sigma^2_\theta\\
  \end{pmatrix}
\end{equation}


\section{Shower development: basic formul\ae}

The longitudinal profile of an electromagnetic shower is described by:
\begin{equation}
\firstder{E}{t} = E_0 p(t) = E_0 k t^\alpha e^{-bt}
\end{equation}
where
$$
k = \frac{b^{\alpha + 1}}{\Gamma(\alpha + 1)}
$$
(with this definition $p(t)$ is normalized to 1 and is therefore a probability
density) and the Euler $\Gamma$ function, defined by:
$$
\Gamma(\alpha) = \int_{0}^\infty t^{\alpha - 1} e^{-t} \diff t
$$
satisfies the well know relation:
$$
\Gamma(\alpha + 1) = \alpha \Gamma(\alpha)
$$
The position of the shower maximum is given by the condition:
$$
\left.\firstder{p}{t}\right|_{\tmax} =
k \tmax^{\alpha - 1} e^{-b\tmax} (\alpha  - b\tmax) = 0
$$
and therefore:
\begin{equation}
\tmax = \frac{\alpha}{b}
\end{equation}
The other two pieces of necessary information are the dependences of $\alpha$
and $b$ on the energy. These are given by the relations:
\begin{equation}
b \approx 0.5
\end{equation}
and:
\begin{equation}
\tmax = \frac{\alpha}{b} = \ln\left(\frac{E_0}{E_c}\right) + C
\end{equation}
where $C=0.5$ for photons and $C=-0.5$ for electrons and $E_c$ is the critical
energy for the material.


\section{Longitudinal moments}

Let's start from the calculation of the lowest order moments of the shower
longitudinal profile around $t=0$. The first one is the mean:
\begin{align}
\left< t \right> &= \mu = \int_{0}^\infty t p(t) ] \diff t =
k \int_{0}^\infty t^{\alpha + 1} e^{-bt} \diff t =\nonumber\\
&= \frac{b^{\alpha + 1}}{\Gamma(\alpha + 1)}
\frac{\Gamma(\alpha + 2)}{b^{\alpha + 2}} = \frac{(\alpha + 1)}{b}
\end{align}
(i.e. the mean of the profile is exactly $1/b$ radiation lengths to the
right of the shower maximum). Along the same lines:
\begin{align}
\left< t^2 \right> = \frac{b^{\alpha + 1}}{\Gamma(\alpha + 1)}
\frac{\Gamma(\alpha + 3)}{b^{\alpha + 3}} =
\frac{(\alpha + 2)(\alpha + 1)}{b^2}
\end{align}
and:
\begin{align}
\left< t^3 \right> = \frac{b^{\alpha + 1}}{\Gamma(\alpha + 1)}
\frac{\Gamma(\alpha + 4)}{b^{\alpha + 4}} =
\frac{(\alpha + 3)(\alpha + 2)(\alpha + 1)}{b^3}
\end{align}

We can apply the usual formul\ae\ for the moments $M_n$ centered around the
mean (as opposed to the ones centered around 0):
\begin{equation}
M_2 = \sigma^2 = \left< t^2 \right> - \mu^2 =
\frac{(\alpha + 1)}{b^2}
\end{equation}
and
\begin{equation}
M_3 = \left< t^3 \right> - 3\mu\sigma^2 - \mu^3 = \frac{2(\alpha + 1)}{b^3}
\end{equation}
The skewness $\gamma$ is given by:
\begin{equation}
\gamma = \frac{M_3}{\sigma^3} = \frac{2}{\sqrt{\alpha + 1}}
\end{equation}

Let's look at the problem from a different perspective, which will hopefully
turn out to be handy in the following. Integrating by parts, we get:
\begin{align*}
\left< t^n \right> & = k \int_{0}^\infty t^n \cdot t^\alpha e^{-bt} \diff t =
k \int_{0}^\infty t^\alpha e^{-bt} \diff\left(\frac{t^{n+1}}{n+1}\right) =\\
&= k\left.\frac{t^{n+1}}{n+1}  t^\alpha e^{-bt}\right|_0^\infty -
k \int_0^\infty \frac{t^{n+1}}{n+1} \left(
\alpha t^{\alpha - 1}e^{-bt} - bt^\alpha e^{-bt}
\right) \diff t =\\
&= \frac{kb}{n+1} \int_{0}^\infty  t^{\alpha + n + 1} e^{-bt} \diff t -
\frac{k\alpha}{n+1} \int_{0}^\infty  t^{\alpha + n} e^{-bt} \diff t =
\frac{b \left< t^{n+1} \right> - \alpha\left< t^n \right>}{n+1}
\end{align*}
from which it follows that:
\begin{equation}
\left< t^{n+1} \right> = \frac{(\alpha + n + 1)}{b}\left< t^n \right>
\end{equation}
For $n = 1$ we get:
$$
\left< t^2 \right> = \frac{(\alpha + 2)}{b}\left< t \right>
$$
or:
\begin{equation}
\sigma^2 = \frac{(\alpha + 2)}{b}\mu - \mu^2
\end{equation}
Whereas for $n = 2$:
$$
\left< t^3 \right> = \frac{(\alpha + 3)}{b}\left< t^2 \right>
$$
which translates into:
\begin{equation}
\gamma = \frac{\mu}{\sigma^3}\left[
  \frac{(\alpha + 3)(\alpha + 2)}{b^2} - 3\sigma^2 - \mu^2
  \right]
\end{equation}

All this equations can be directly verified by plugging in the expressions
for $\mu$, $\sigma$ and $\gamma$ explicitly obtained before, but the hope is
to generalize them to the case in which we don't sample the entire shower
(see the following section).


\section{Longitudinal moments over a finite interval}

We can generalize the previous relations to the case in which we only
sample a finite fraction of the longitudinal shower development,
say between $t_1$ and $t_2$. The formalism is essentially identical, except
for the fact that now we're dealing with a probability density function
over a finite interval:
$$
p_{\rm f}(t) = k_{\rm f} t^\alpha e^{-bt}
$$
with $k_{\rm f}$ being:
$$
k_{\rm f} = \frac{1}{\int_{t_1}^{t_2} t^\alpha e^{-bt} \diff t}
$$
(physically $k_{\rm f}$ is the ratio between the raw energy deposited in the
calorimeter and the true energy of the particle). So now we have:
\begin{equation}
\left< t^{n+1} \right> = \frac{(\alpha + n + 1)}{b}\left< t^n \right> -
\left.\frac{k_{\rm f}}{b} t^{(\alpha + n + 1)} e^{-bt}\right|_{t_1}^{t_2}
\end{equation}
and therefore:
\begin{equation}
\left< t^2 \right> = \frac{(\alpha + 2)}{b}\left< t \right> -
\frac{k_{\rm f}}{b}
\left[t_2^{(\alpha + 2)} e^{-bt_2} - t_1^{(\alpha + 2)} e^{-bt_1}\right]
\end{equation}
and:
\begin{equation}
\left< t^3 \right> = \frac{(\alpha + 3)}{b}\left< t^2 \right> -
\frac{k_{\rm f}}{b}
\left[t_2^{(\alpha + 3)} e^{-bt_2} - t_1^{(\alpha + 3)} e^{-bt_1}\right]
\end{equation}

Some more formula that might turn out to be useful for the normalization
of the skewness to the expected value for electronmagnetic showers.
The moments of the longitudinal distribution can be written as a function
of the incomplete gamma function, defined as:
\begin{equation}
\gamma(\alpha, t) = \frac{1}{\Gamma(\alpha)}
\int_{0}^{t} t^{\alpha - 1} e^{-t} \diff t
\end{equation}
from which it follows that:
\begin{equation}
\int_{t_1}^{t_2} t^{\alpha} e^{-bt} \diff t =
\frac{\Gamma(\alpha+1)}{b^{\alpha+1}}
\left( \gamma(\alpha+1, bt_2) - \gamma(\alpha+1, bt_1) \right)
\end{equation}
If we define:
$$
{\mathcal G}(\alpha, b, t_1, t_2) = \frac{\Gamma(\alpha)}{b^{\alpha}}
\left( \gamma(\alpha, bt_2) - \gamma(\alpha, bt_1) \right)
$$
we have:
\begin{align}
\left< t^n \right> =
\frac{{\mathcal G}(\alpha + n + 1, b, t_1, t_2)}
     {{\mathcal G}(\alpha + 1, b, t_1, t_2)}
\end{align}



\clearpage
\appendix

\emph{Caution: the stuff in the appendix is mostly crap, at this time.
I'll move it into appropriate sections as soon as it's in a reasonable shape
(and, of course, this does not mean that people should not take a look).}



Let's go back to the basic equation for the principal eigenvector:
$$
\itm\mathbf{e}^1 = \lambda_1\mathbf{e}^1
$$
Doing a full error propagation is not easy, since in this equation we do have
error on the six independent components of the inertia tensor, as well as on
the eigenvalue $\lambda_1$ we've just calculated. The errors on the
$\itc{ij}$ are reasonably easy to calculate, starting from the errors
associated with the finite dimensions of the crystals. On the other side
the propagation of the errors to $\lambda_1$ is not trivial, as the
expression is complicated. On top of that, these different error are not
indipendent from each other, as $\lambda_1$ is calculated starting from the
component of the inertia tensor.

The solution to this equation is:
\begin{align}
e^1_x &= \frac{1}{\sqrt{1 + \frac{A^2}{B^2} + \frac{A^2}{C^2}}}\\
e^1_y &= \frac{1}{\sqrt{1 + \frac{B^2}{A^2} + \frac{B^2}{C^2}}}\\
e^1_z &= \frac{1}{\sqrt{1 + \frac{C^2}{A^2} + \frac{C^2}{B^2}}}
\end{align}
where:
\begin{align}
A &= \itc{yz}(\itc{xx} - \lambda_1) - \itc{xy}\itc{xz}\\
B &= \itc{xz}(\itc{yy} - \lambda_1) - \itc{xy}\itc{yz}\\
C &= \itc{xy}(\itc{zz} - \lambda_1) - \itc{xz}\itc{yz}
\end{align}



\begin{thebibliography}{100}
\bibitem{goldstein}H.~Goldstein, \emph{Classical mechanics}.
\bibitem{landau}L.~D.~Landau, E.~M.~Lif\^sic, \emph{Mechanics}.
\bibitem{wolfram}\url{http://mathworld.wolfram.com/CubicFormula.html}
\bibitem{pdg}PDG Review of Particle Physics,
  \emph{Physics Letters B.} (2004) {\bf 592}
\bibitem{errors}T.~Soler and B.~H.~W.~van~Gelder,
  \emph{Geophys. J. Int.} (1991) {\bf 105}, 537--546.
\bibitem{errors_corr}T.~Soler and B.~H.~W.~van~Gelder,
  \emph{Geophys. J. Int.} (2006) {\bf 165}, 382.
\end{thebibliography}



\end{document}
