removing the setting of AOD track references for TPC only tracks

[u/mrichter/AliRoot.git] / PWG2 / FORWARD / doc / doc.tex

\documentclass[11pt]{article}
\usepackage[margin=2cm,twoside,a4paper]{geometry}
\usepackage{amstext}
\usepackage{amsmath}
\usepackage[ruled,vlined,linesnumbered]{algorithm2e}
\usepackage{graphicx}
\usepackage{color}
\usepackage{units}
\usepackage{listings}
\def\AlwaysText#1{\ifmmode\relax\text{#1}\else #1\fi}
\newcommand{\AbbrName}[1]{\AlwaysText{{\scshape #1}}}
\newcommand{\SPD}{\AbbrName{spd}}
\newcommand{\ESD}{\AbbrName{esd}}
\newcommand{\AOD}{\AbbrName{aod}}
\newcommand{\INEL}{\AbbrName{inel}}
\newcommand{\INELONE}{$\AbbrName{inel}>0$}
\newcommand{\NSD}{\AbbrName{nsd}}
\newcommand{\FMD}[1][]{\AbbrName{fmd\ifx|#1|\else#1\fi}}
\newcommand{\dndetadphi}[1][]{{\ensuremath% 
    \ifx|#1|\else\left.\fi%
    \frac{d^2N_{ch}}{d\eta\,d\varphi}%
    \ifx|#1|\else\right|_{#1}\fi%
}}
\newcommand{\landau}[1]{{\ensuremath% 
    \text{landau}\left(#1\right)}}
\newcommand{\dndeta}[1][]{{\ensuremath% 
    \ifx|#1|\else\left.\fi%
    \frac{1}{N}\frac{dN_{ch}}{d\eta}%
    \ifx|#1|\else\right|_{#1}\fi%
}}
\newcommand{\MC}{\AlwaysText{MC}}
\newcommand{\Nsel}{{\ensuremath N_{v,\text{trigger,vertex}}}}
\newcommand{\GeV}[1]{\unit[#1]{\AlwaysText{GeV}}}
\newcommand{\cm}[1]{\unit[#1]{\AlwaysText{cm}}}

\setlength{\parskip}{1ex}
\setlength{\parindent}{0em}

\title{Analysing the FMD data for $\dndeta$}
\author{Christian Holm
  Christensen\thanks{\texttt{$\langle$cholm@nbi.dk$\rangle$}}\quad\&\quad
  Hans Hjersing Dalsgaard\thanks{\texttt{$\langle$canute@nbi.dk$\rangle$}}\\ 
  Niels Bohr Institute\\
  University of Copenhagen}
\date{\today}
\begin{document}
\maketitle 

\tableofcontents 
\section{Introduction}

This document describes the steps performed in the analysis of the
charged particle multiplicity in the forward pseudo--rapidity
regions. 

The analysis is performed as a two--step process.  
\begin{enumerate}
\item The Event--Summary--Data (\ESD{}) is processed event--by--event
  and passed through a number of algorithms, and
  $\dndetadphi$ for each event is output to an Analysis--Object--Data
  (\AOD{}) tree.
\item The \AOD{} data is read in and the sub--sample of the data under
  investigation is selected (e.g., \INEL{}, \INELONE{}, \NSD{}, or
  some centrality class) and the $\dndetadphi$ histogram read in for
  those events to build up $\dndeta$
\end{enumerate}
The details of each step above will be expanded upon in the
following. 

\section{Generating $\dndetadphi[i]$ event--by--event}

When reading in the \ESD{}s and generating the $\dndetadphi$
event--by--event the following steps are taken (in order) for each
event $i$
\begin{description}
\item[Event inspection] The global properties of the event is
  determined, including the trigger type, vertex $z$ coordinate, and
  whether this is a low--flux event or not. 
\item[Sharing filter] The \ESD{} object is read in and corrected for
  sharing.  The result is a new \ESD{} object.
\item[Density calculator] The (possibly un--corrected) \ESD{} object
  is then inspected and an inclusive, per--ring charged particle
  density $\dndetadphi[incl,r,v,i]$  is made.  This
  calculation depends in general upon the interaction
  vertex\footnote{In the following simply labelled 'primary vertex' or
    'vertex'.} position along the $z$ axis ($v_z$).  
\item[Corrections] The 5 $\dndetadphi[incl,r,v,i]$ are
  corrected for secondary production, event selection efficiency, and
  possibly the sharing efficiency.  These corrections are highly
  dependent on the vertex $z$ coordinate.  The result is an per--ring,
  charged primary particle density $\dndetadphi[r,v,i]$
\item[Histogram collector] Finally, the 5 $\dndetadphi[r,v,i]$ are
  summed into a single $\dndetadphi[v,i]$ histogram, taking care of
  the overlaps between the detector rings.  In principle, this
  histogram is independent of the vertex, except that the
  pseudo--rapidity range, and possible holes in that range, depends on
  $v_z$ --- or rather the bin in which the $v_z$ falls. 
\end{description}

Each of these steps will be detailed in the following. 

\subsection{Event inspection}

The first thing to do, is to inspect the event for triggers.  A number
of trigger bits, like \INEL{}, \INELONE{}, \NSD{}, and so on is then
propagated to the \AOD{} output.  

Next, the number of \emph{tracklets} reconstructed in the Silicon
Pixel Detector (\SPD{}) compared to a threshold.  If the number of
track--lets falls belows this threshold, the event is consider a
low--flux event.   

Just after the sharing filter (described below) but before any more
processing, the vertex information is queried.  If there is no vertex
information, or if the vertex $z$ coordinate is outside the
pre--defined range, then no further processing takes place. 

\subsection{Sharing filter}

The \FMD{} \ESD{} object contains the scaled energy deposited $\Delta
E/\Delta E_{mip}$ for each of the 51,200 strips.  The \FMD{} is
organised in 3 \emph{sub--detectors} \FMD{1}, \FMD{2}, and \FMD{3}, each
consisting of 1 (\FMD{1}) or 2 (\FMD{2} and \FMD{3}) \emph{rings}.
The rings fall into two types: \emph{Inner} or \emph{outer} rings.
Each ring is in turn  azimuthal divided into \emph{sectors}, and each
sector is radially divided into \emph{strips}.  How many sectors,
strips, as well as the $\eta$ coverage is given in
\tablename~\ref{tab:fmd:overview}. 

\begin{table}[htbp]
  \begin{center}
    \caption{Physical dimensions of Si segments and strips.}
    \label{tab:fmd:overview}
    \vglue0.2cm
    \begin{tabular}{|c|cc|cr@{\space--\space}l|r@{\space--\space}l|}
      \hline
      \textbf{Sub--detector/} &
      \textbf{Azimuthal}&
      \textbf{Radial} &
      $z$ &
      \multicolumn{2}{c|}{\textbf{$r$}} &
      \multicolumn{2}{c|}{\textbf{$\eta$}} \\
      \textbf{Ring}& 
      \textbf{sectors} &
      \textbf{strips} & 
      \textbf{[cm]} &
      \multicolumn{2}{c|}{\textbf{range [cm]}} &
      \multicolumn{2}{c|}{\textbf{coverage}} \\
      \hline
      FMD1i & 20& 512& 320  &  4.2& 17.2& 3.68&  5.03\\
      FMD2i & 20& 512&  83.4&  4.2& 17.2& 2.28&  3.68\\
      FMD2o & 40& 256&  75.2& 15.4& 28.4& 1.70&  2.29\\
      FMD3i & 20& 512& -75.2&  4.2& 17.2&-2.29& -1.70\\
      FMD3o & 40& 256& -83.4& 15.4& 28.4&-3.40& -2.01\\
      \hline
    \end{tabular}
  \end{center}
\end{table}

A particle originating from the vertex can, because of it's incident
angle on the \FMD{} sensors traverse more than one strip.  That means
that the energy loss of the particle is distributed over 1 or more
strips.  The signal in each strip should therefore possibly merged
with it's neighbor strip signals to properly reconstruct the energy
loss of a single particle.  

The effect is most pronounced in low--flux events, like proton--proton
collisions or peripheral Pb--Pb collisions, while in high--flux events
the hit density is so high that most likely each and every strip will
be hit and the effect cancel out on average. 

Since the particles travel more or less in straight lines toward the
\FMD{} sensors, the sharing effect predominantly in the $r$ or
\emph{strip} direction.  Only neighboring strips in a given sector is
therefor investigated for this effect.  

Algorithm~\ref{algo:sharing} is applied to the signals in a given
sector.

\begin{algorithm}[htpb]
  \SetKwData{usedThis}{current strip used}
  \SetKwData{usedPrev}{previous strip used}
  \SetKwData{Output}{output}
  \SetKwData{Input}{input}
  \SetKwData{Nstr}{\# strips}
  \SetKwData{Signal}{current}
  \SetKwData{Eta}{$\eta$}
  \SetKwData{prevE}{previous strip signal} 
  \SetKwData{nextE}{next strip signal} 
  \SetKwData{lowFlux}{low flux flag} 
  \SetKwFunction{SignalInStrip}{SignalInStrip}
  \SetKwFunction{MultiplicityOfStrip}{MultiplicityOfStrip}
  \usedThis $\leftarrow$ false\;
  \usedPrev $\leftarrow$ false\;
  \For{$t\leftarrow1$ \KwTo \Nstr}{ 
    \Output${}_t\leftarrow 0$\;
    \Signal $\leftarrow$ \SignalInStrip($t$)\;

    \uIf{\Signal is not valid}{ 
      \Output${}_t \leftarrow$ invalid\;
    }
    \uElseIf{\Signal is 0}{ 
      \Output${}_t \leftarrow$ 0\;
    }
    \Else{
      \Eta$\leftarrow$ $\eta$ of \Input${}_t$\;
      \prevE$\leftarrow$ 0\;
      \nextE$\leftarrow$ 0\;
      \lIf{$t \ne 1$}{ 
        \prevE$\leftarrow$ \SignalInStrip($t-1$)\;
      }
      \lIf{$t \ne $\Nstr}{ 
        \nextE$\leftarrow$ \SignalInStrip($t+1$)\;
      }
      \Output${}_t\leftarrow$
      \MultiplicityOfStrip(\Signal,\Eta,\prevE,\nextE,\\
      \hfill\lowFlux,$t$,\usedPrev,\usedThis)\;
    }   
  }
  \caption{Sharing correction}
  \label{algo:sharing}
\end{algorithm}

Here the function \FuncSty{SignalInStrip}($t$) returns the properly
path--length corrected signal in strip $t$.  The function
\FuncSty{MultiplicityInStrip} is where the real processing takes
place (see page \pageref{func:MultiplicityInStrip}). 

\begin{function}[htbp]
  \caption{MultiplicityInStrip(\DataSty{current},$\eta$,\DataSty{previous},\DataSty{next},\DataSty{low
      flux flag},\DataSty{previous signal used},\DataSty{this signal
      used})} 
  \label{func:MultiplicityInStrip}
  \SetKwData{Current}{current} 
  \SetKwData{Next}{next} 
  \SetKwData{Previous}{previous} 
  \SetKwData{lowFlux}{low flux flag}
  \SetKwData{usedPrev}{previous signal used}
  \SetKwData{usedThis}{this signal used}
  \SetKwData{lowCut}{low cut}
  \SetKwData{total}{Total}
  \SetKwData{highCut}{high cut}
  \SetKwData{Eta}{$\eta$}  
  \SetKwFunction{GetHighCut}{GetHighCut}
  \If{\Current is very large or \Current $<$ \lowCut} {
    \usedThis $\leftarrow$ false\;
    \usedPrev $\leftarrow$ false\;
    \Return{0}
  }
  \If{\usedThis}{ 
    \usedThis $\leftarrow$ false\;
    \usedPrev $\leftarrow$ true\;
    \Return{0}
  }
  \highCut $\leftarrow$ \GetHighCut($t$,\Eta)\;
  \If{\Current $<$ \Next and \Next $>$ \highCut and \lowFlux set}{ 
    \usedThis $\leftarrow$ false\;
    \usedPrev $\leftarrow$ false\;
    \Return{0}
  }
  \total $\leftarrow$ \Current\;
  \lIf{\lowCut $<$ \Previous $<$ \highCut and not \usedPrev}{ 
    \total $\leftarrow$ \total + \Previous\;
  }
  \If{\lowCut $<$ \Next $<$ \highCut}{ 
    \total $\leftarrow$ \total + \Next\;  
    \usedThis $\leftarrow$ true\;
  }
  \eIf{\total $>$ 0}{ 
    \usedPrev $\leftarrow$ true\;
    \Return{\total}
  }{
    \usedPrev $\leftarrow$ false\;
    \usedThis $\leftarrow$ false\;
    \Return{0}
  }
\end{function}
Here, the function \FuncSty{GetHighCut} evaluates a Landau
distribution fitted to the energy spectrum in the $\eta$ bin
specified.  It returns 
$$
\Delta_{mp} - 2 w
$$
where $\Delta_{mp}$ is the most probable energy loss, and $w$ is the
width of the Landau distribution.  

The \KwSty{if} in line 5, says that if the previous strip was merged
with current one, and the signal of the current strip was added to
that, then we the current signal is set to 0, and we mark it as used
for the next iteration (\DataSty{previous signal
  used}$\leftarrow$true). 

The \KwSty{if} in line 10 checks if the current signal is smaller than
the next signal, the next signal is larger than the upper cut defined
above, and if we have a low--flux event.  If that condition is met,
then the current signal is the smaller of two possible candidates for
merging, and it should be merged into the next signal.  Note, that
this \emph{only} applies in low--flux events.  

On line 15, we test if the previous signal lies between our low and
high cuts, and if it has not been marked as being used.  If so, we add
it to our current signal.  

The next \KwSty{if} on line 16 checks if the next signal is within our
cut bounds.  If so, we add that signal to the current signal and mark
it as used for the next iteration (\DataSty{this signal
  used}$\leftarrow$true).  It will then be zero'ed on the next
iteration by the condition on line 6.

Finally, if our signal is still larger than 0, we return the signal
and mark this signal as used (\DataSty{previous signal
  used}$\leftarrow$true) so that it will not be used in the next
iteration. Otherwise, we mark the current signal and the next signal
as unused and return a 0. 


\subsection{Density calculator}

The density calculator loops over all the strip signals and calculates
the inclusive (primary + secondary) charged particle density in
pre-defined $(\eta,\varphi)$ bins.  

\subsubsection{Inclusive number of charged particles} 

If the event is classified as a low--flux event, then the number of
charged particles in a given by a simple threshold: 
\begin{align}
  N_{ch,t} &= \left\{
    \begin{array}{cl}
      0 & \Delta_t < \text{low cut}\\ 
      1 & \Delta_t \ge \text{low cut}\\ 
    \end{array}\right.\quad,
\end{align}
where $t$ is the strip identifier, $\Delta_t$ is the scaled energy
deposition in that strip, and 'low cut' is a predefined
cut\footnote{This low--flux mode is perhaps deprecated.}.  

For high flux events, the number charged particles in a strip is
calculated using multiple Landau distributions fitted to the energy
loss spectrum at a given $\eta$ value.
\begin{align}
  \Delta_{i,mp} &= i (\Delta_{1,mp}+ \xi_1 \log(i))\nonumber\\
  \xi_i         &= i\xi_1\nonumber\\
  \sigma_i      &= \sqrt{i}\sigma_1\nonumber\\
  N_{ch,t} &= \frac{\sum_i^{N_{max}}
    i\,a_i\,F(\Delta_t;\Delta_{i,mp},\xi_i,\sigma_i)}{
    \sum_i^{N_{max}}\,a_i\,F(\Delta_t;\Delta_{i,mp},\xi_i,\sigma_i)}\quad,
\end{align}
where $F(x;\Delta_{mp},\xi,\sigma)$ is the evaluation of the Landau
distribution $f_L$ with most probable value $\Delta_{mp}$ and width
$\xi$, folded with a Gaussian distribution with spread $\sigma$
\cite{nim:b1:16,phyrev:a28:615}.
$$
F(x;\Delta_{mp},\xi,\sigma) = \frac{1}{\sigma \sqrt{2 \pi}}
\int_{-\infty}^{+\infty} d\Delta' f_{L}(x;\Delta',\xi)
\exp{-\frac{(\Delta_{mp}-\Delta')^2}{2\sigma^2}}
$$
$\Delta_{1,mp}$, $\xi_1$, and $\sigma_1$ are the parameters for the
first MIP peak, $a_1=1$, and $a_i$ is the relative weight of the
$i^{\text{th}}$ MIP peak.   The parameters $\Delta_{1,mp}, \xi_1,
\sigma_1, a_2, \ldots a_{N_{max}}$ are obtained by fitting 
$$
\sum_{i=1}^{j} F(x;\Delta_{i,mp},\xi_{i},\sigma_i) 
$$
for increasing $j$ to the energy loss spectra in separate $\eta$
bins. 

\subsubsection{Acceptance and double-hit corrections}

Before the signal $N_{ch,t}$ can be added to the $(\eta,\varphi)$
bin in one of the 5 per--ring histograms, it needs to be corrected for
the $\varphi$ acceptance of the strip, as well as a correction for
double hits in low--flux events.   

The acceptance correction is only applicable where the strip length
does not cover the full sector.  This is the case for the outer strips
in both the inner and outer type rings.  The acceptance correction is
then simply 
\begin{align}
  \label{eq:acc_corr}
  a_t &= \frac{l_t}{\Delta\varphi}\quad
\end{align}
where $l_t$ is the strip length in radians at constant $r$, and
$\Delta\varphi$ is $2\pi$ divided by the number of sectors in the
ring (20 for inner type rings, and 40 for outer type rings). 

Even in low--flux events it is possible that more than one particle
hits a strip.  However, for low--flux events, it is not possible to
reconstruct the 3\textsuperscript{rd} nor even the
2\textsuperscript{nd} MIP peak in the energy loss spectrum.
Therefore, the strip signal needs to be corrected to the average
number of particle impinging on a strip at a given $\eta$\footnote{As
  before, this low--flux mode is deprecated and this correction is not
  applied}.  
\begin{align}
  d_{t,r}(\eta) &= \left\{
    \begin{array}{cl} \langle
      n_{t,r}(\eta)\rangle & \text{low flux}\\
      1 & \text{high flux}
    \end{array}\right.\quad.
\end{align}
This correction is calculated separately for each ring in $\eta$ bins
from simulated events like
\begin{align}
  \label{eq:double_hit_corr} 
  \langle n_{t,r}(\eta)\rangle &= \frac{
    \sum_i N_{\text{strips},r,i}(\eta)}{
    \sum_i N_{ch,r,i}(\eta)}\quad,
\end{align}
where $N_{\text{strips},r,i}(\eta)$ is the number of strips in ring
$r$ in a given $\eta$ bin that have been hit in the $i^{\text{th}}$
event, and $N_{ch,r,i}(\eta)$ is the number of charged particles that
fell within the same $\eta$ bin in the $i^{\text{th}}$ event.

The final $(\eta,\varphi)$ content of the 5 output vertex dependent,
per--ring histograms of the inclusive charged particle density is then
given by
\begin{align}
  \label{eq:density}
  \dndetadphi[incl,r,v,i(\eta,\varphi)] &= \sum_t^{t\in(\eta,\varphi)}
  N_{ch,t}\,a_t\,d_{t,r}(\eta)
\end{align}
where $t$ runs over the strips in the $(\eta,\varphi)$ bin. 

\subsection{Corrections}

The corrections code receives the five vertex dependent,
per--ring histograms of the inclusive charged particle density
$\dndetadphi[incl,r,v,i]$ from the density calculator and applies
three corrections 

\subsubsection{Secondary correction}
%%
%%                hHits_FMD<d><r>_vtx<v> 
%% hCorrection = -----------------------
%%                hPrimary_FMD_<r>_vtx<v>
%%
%% where 
%% - hPrimary_FMD_<r>_vtx<vtx> is 2D of eta,phi for all primary ch
%%   particles
%% - hHits_FMD<d><r>_vtx<v>  is 2D of eta,phi for all track-refs that
%%   hit the FMD - The 2D version of hMCHits_nocuts_FMD<d><r>_vtx<v>
%%   used below. 
  This is a 2 dimensional histogram generated from simulations of the
  ratio of primary particles to the total number of particles that
  fall within an $(\eta,\varphi)$ bin for a given vertex bin

  \begin{align}
    \label{eq:secondary}
    s_v(\eta,\varphi) &=
    \frac{\sum_i^{\Nsel}N_{ch,\text{primary},i}(\eta,\varphi)}{
      \sum_i^{\Nsel}N_{ch,\text{\FMD{}},i}(\eta,\varphi)}\quad,
  \end{align}
  where $\Nsel$ is the number of events with a valid trigger and a
  vertex in bin $v$, and $N_{ch,\FMD{},i}$ is the total number of
  charged particles that hit the \FMD{} in event $i$ in the specified
  $(\eta,\varphi)$ bin and $N_{ch,\text{primary},i}$ is number of
  primary charged particles in event $i$ within the specified
  $(\eta,\varphi)$ bin.

  $N_{ch,primary}(\eta,\varphi)$ is given by summing over the charged
  particles labelled as primaries \emph{at the time of the collision}
  as defined in the simulation code.  That is, it is the number of
  primaries within the bin at the collision point --- not at the
  \FMD{}.  

\subsubsection{Event selection efficiency}
%% AliFMDAnalysisTaskGenerateCorrection
%% 
%% Analysed_FMD<r>_vtx<v> 2D histogram (eta,phi)
%% EventsSelected[<v>]    bin content 
%% Inel_FMD<r>_vtx<v>     2D histogram (eta,phi) 
%% EventsAll[<v>]         bin content 
%%
%% correction is 
%%
%%      Analysed_FMD<d><r>_vtx<v> / EventsSelected[<v>]
%%    ---------------------------------------------------
%%      Inel_FMD<r>_vtx<v> / EventsAll[<v>] 
%% 
%% Inel_FMD<r>_vtx<v>        is d^2N_incl,ch/d\eta d\varphi for all
%%                           events  
%% Analysed_FMD<d><r>_vtx<v> is d^2N_{incl,ch}/d\eta d\varphi for
%%                           events with trigger _and_ vertex
%% EventsSelected            is dN/dv_z for events with trigger _and_
%%                           vertex
%% EventsAll                 is dN/dv_z for all events 
  This correction is made from simulations.  It is a 2--dimensional
  histogram of the ratio of average number of charged particles in
  events that are triggered as interactions \emph{and} have vertex, to
  the average number of charged particles in all events in
  $(\eta,\varphi)$ bins
  \begin{align}
    \label{eq:event_sel_eff}
    e_{v}(\eta,\varphi) &= 
    \frac{\frac{1}{N_{v}} \sum_{i}^{N_{v}} N_{ch,\MC{},i}(\eta,\varphi)}{
      \frac{1}{\Nsel{}}\sum_{i}^{\Nsel{}}
      N_{ch,\MC{},i}(\eta,\varphi)}\nonumber\\
    &= \frac{\Nsel{}}{N_{v}} \frac{\sum_{i}^{N_{v}} 
      N_{ch,\MC{},i}(\eta,\varphi)}{\sum_{i}^{\Nsel{}}
      N_{ch,\MC{},i}(\eta,\varphi)}\quad,
  \end{align}
  where $N_{v}$ is the total number of events in vertex bin $v$, and
  $\Nsel{}$ is the number of events with a collision trigger and
  vertex within the predefined range, for vertex bin $v$.
  $N_{ch,\MC{},i}$ is the number of charged 'MC truth' particles
  within the specified $(\eta,\varphi)$ range for event $i$.

\subsubsection{Sharing correction efficiency}
%% hits_NoCuts_FMD<d><r>_vtxbin<v>_proj 
%% hMCHits_nocuts_FMD<d><r>_vtxbin<v>
%% 
%% Correction is 
%% 
%%     hits_NoCuts_FMD<d><r>_vtxbin<v>_proj / hEvents[<v>]
%%    -----------------------------------------------------
%%     hMChits_nocuts_FMD<d><r>_vtxbin<v> / nMCEventsNoCuts[<v>] 
%%
%% Or 
%%
%%     sum_eta hits_FMD<d><r>_vtxbin<v> / EventSel1D / hEvents[<v>]
%%   ---------------------------------------------------------------
%%     hMCHits_nocuts_FMD<d><r>_vtxbin<v> / nMCEventsNoCuts[<v>] 
%%   
%% - hits_NoCuts_FMD<d><r>_vtxbin<v>_proj 
%%   A projection of hits_FMD<d><r>_vtxbin<v> onto the eta axis and
%%   possibly scaled to EventSel1D. 
%% - EventSel1D - the 1D event selection efficiency (per vertex) 
%%   This is the ratio of hEventAll divided hEventsSelected, with
%%   errors calculated as 
%% 
%%      da = sqrt(e_sel^2 + e_all^2)
%%      a  = c_all - c_sel 
%%      e  = sqrt(da^2 + a^2 / (c_ell*c_sel)^2)
%%         = sqrt(e_sel^2 + e_all^2 + (c_all-c_sel)^2 /
%%         (c_ell*c_sel)^2)
%%
%%   where hEventsAll[<v>] is the total number of events in the v bin,
%%   and hEventsSelected[<v>] is the notoal number of events with
%%   trigger and vertex in the v bin. Both for MC
%% - hits_FMD<d><r>_vtxbin<v>
%%   from AliFMDAnalysisTaskBackgroundCorrection - the uncorrected
%%   multiplicity in each ring - i.e.,
%%   \dndetadphi[incl,r,v,i(\eta,\varphi)] 
%% - hMCHits_nocuts_FMD<d><r>_vtxbin<v> - from
%%   AliFMDAnalysisTaskSharing - a 1D histogram in eta - sum of eta of
%%   track-refs that hit FMD with some requirements on the neighbors  
%% - hEvents[<v>] - bin content of dN/dv_{z} 
%% - nMCEventsNoCuts -  hPrimEvents[<v>] - bin content of dN/dv_{z}
%%   for MC events - from AliFMDAnalysisTaskSharing - number of not
%%   empty mc events per vertex bin 
%% 
This is a one dimension histogram in $\eta$, generated from running
the full analysis on simulated data\footnote{With improved energy loss
  fits, this correction is redundant and not used} .
  \begin{align}
    \label{eq:sharing_corr}
    m_v(\eta) &=
    \frac{\frac1{\Nsel}\sum_i^{\Nsel{}}N_{ch,\text{\FMD{}},i}(\eta)}{
      \frac{1}{e_{v}}\frac1{\Nsel{}}\left.\frac{dN_{ch}}{d\eta}
      \right|_{incl,r,v(\eta)}^{\MC{}}}
    \nonumber\\
    &=
    \frac{\frac1{\Nsel}\sum_i^{\Nsel{}}N_{ch,\text{\FMD{}},i}(\eta)}{
      \frac1{N_v}\left.\frac{dN_{ch}}{d\eta}\right|_{incl,r,v(\eta)}^{\MC}}
    \nonumber\\
    &= 
    \frac{N_v}{\Nsel{}}
    \frac{\sum_i^{\Nsel{}}N_{ch,\text{\FMD{}},i}(\eta)}{
      \left.\frac{dN_{ch}}{d\eta}\right|_{incl,r,v(\eta)}^{\MC}}\\
    \intertext{since}
    e_{v} &= \frac{N_{v}}{\Nsel{}}\quad,
%% \pm 
%%    \sqrt{\sqrt{\Nsel{}}^2 + \sqrt{N_v^2} +
%%      \frac{(N_v-\Nsel{})^2}{(\Nsel{}N_v)^2}}\nonumber \\
%%    &= \frac{N_{v}}{\Nsel{}} \pm 
%%    \sqrt{\Nsel{} + N_v + \frac1{\Nsel{}^2} + \frac1{N_v^2} -
%%      \frac1{\Nsel{}N_v}}\quad, 
  \end{align}
  and $N_v$ is the total number of simulated events with vertex in bin
  $v$, $\Nsel{}$ is the number of simulated events with a collision
  trigger and a reconstructed vertex, $N_{ch,\text{\FMD{}},i}(\eta)$ is
  the total number of 'MC truth' charged particles that impinge on the
  \FMD{} in the given eta  bin in event $i$, and 
  \begin{align}
    \left.\frac{dN_{ch}}{d\eta}\right|_{incl,r,v(\eta)}^{\MC{}} &= 
    \int_0^{2\pi}d\varphi\dndetadphi[incl,r,v(\eta,\varphi)]
    \nonumber\\
    &=  \int_0^{2\pi}d\varphi\sum_i^{\Nsel{}}\dndetadphi[incl,i,r,v(\eta,\varphi)]\quad,
  \end{align}
  where the integrand is \eqref{eq:density} for simulated data i.e.,
  this is the inclusive total charge particle density in a given $v$
  bin for ring $r$ in event $i$. 

Putting equations \eqref{eq:secondary}, \eqref{eq:event_sel_eff}, and
\eqref{eq:sharing_corr} together we get the total correction
$c_v(\eta,\varphi)$ to be 
\begin{align}
  \label{eq:total_corr}
  c_v(\eta,\varphi) =\ &
  s_v(\eta,\varphi)\,e_v(\eta,\varphi)\,m_v(\eta)\nonumber\\ 
  =\ & \frac{\sum_i^{\Nsel}N_{ch,\text{primary},i}(\eta,\varphi)}{
      \sum_i^{\Nsel}N_{ch,\text{\FMD{}},i}(\eta,\varphi)}
    \nonumber\\
    & \times\ \frac{\Nsel{}}{N_{v}} \frac{\sum_{i}^{N_{v}} 
      N_{ch,\MC{},i}(\eta,\varphi)}{\sum_{i}^{\Nsel{}}
      N_{ch,\MC{},i}(\eta,\varphi)}\nonumber\\
    & \times\ \frac{N_v}{\Nsel{}}
    \frac{\sum_i^{\Nsel{}}N_{ch,\text{\FMD{}},i}(\eta)}{
      \left.\frac{dN_{ch}}{d\eta}\right|_{incl,r,v(\eta)}^{\MC}}\nonumber\\
  =\ & \frac{\sum_i^{\Nsel}N_{ch,\text{primary},i}(\eta,\varphi)}{
      \sum_i^{\Nsel}N_{ch,\text{\FMD{}},i}(\eta,\varphi)}
    \frac{\sum_{i}^{N_{v}} N_{ch,\MC{},i}(\eta,\varphi)}{
      \sum_{i}^{\Nsel{}} N_{ch,\MC{},i}(\eta,\varphi)}
    \frac{\sum_i^{\Nsel{}}N_{ch,\text{\FMD{}},i}(\eta)}{
      \left.\frac{dN_{ch}}{d\eta}\right|_{incl,r,v(\eta)}^{\MC}}\quad,
\end{align}
which\footnote{Note that $\sum_{i}^{\Nsel{}} N_{ch,\FMD{},i}(\eta) =
  \int_0^{2\pi}d\varphi\,\sum_{i}^{\Nsel{}}
  N_{ch,\FMD{},i}(\eta,\varphi)$} is independent of the overall
number of events $N_v$ and $\Nsel{}$.

The 5 output vertex dependent, per--ring histograms of the primary
charged particle density is then given by
\begin{align}
  \dndetadphi[r,v,i(\eta,\varphi)] &=
  c_v(\eta,\varphi)\dndetadphi[incl,r,v,i(\eta,\varphi)]
\end{align}

\subsection{Histogram collector}

The histogram collector collects the information from the 5 vertex
dependent, per--ring histograms of the primary charged particle
density $\dndetadphi[r,v,i]$ into a single vertex dependent histogram
of the charged particle density $\dndetadphi[v,i]$.  

To do this, it first calculates, for each vertex bin, the $\eta$ bin
range to use for each ring.  It investigates the secondary correction
maps $s_v(\eta,\varphi)$ to find the edges of the map.  The edges are
given by the $\eta$ range where $s_v(\eta,\varphi)$ is larger than
some threshold\footnote{Typically $t_s\approx 0.1$.}  $t_s$. The code
applies safety margin of a $N_{cut}$ bins\footnote{Typically
  $N_{cut}=1$.}, to ensure that the data selected does not have too
large corrections associated with it.

It then loops over the bins in the defined $\eta$ range and sums the
contributions from each of the 5 histograms.  In the $\eta$ ranges
where two rings overlap, the collector calculates the average and adds
the errors in quadrature.

The output vertex dependent histogram of the primary
charged particle density is then given by
\begin{align}
  \label{eq:superhist}
  \dndetadphi[v,i(\eta,\varphi)] &=
  \frac{1}{N_{r\in(\eta,\varphi)}}\sum_{r}^{r\in(\eta,\varphi)}  
  \dndetadphi[r,v,i(\eta,\varphi)]\\
  \delta\left[\dndetadphi[v,i(\eta,\varphi)]\right] &=
  \frac{1}{N_{r\in(\eta,\varphi)}}\sqrt{\sum_{r}^{r\in(\eta,\varphi)}   
    \delta\left[\dndetadphi[r,v,i(\eta,\varphi)]\right]^2}
  \quad,
\end{align}
where $N_{r\in(\eta,\varphi)}$ is the number of overlapping histograms
in the given $(\eta,\varphi)$ bin. 

The histogram collector stores the found $\eta$ ranges in the
underflow bin of the histogram produced.  The content of the overflow
bins are 
\begin{align}
  \label{eq:overflow}
  I_{v,i}(\eta) &= 
  \frac{1}{N_{r\in(\eta)}}
  \sum_{r}^{r\in(\eta)} \left\{\begin{array}{cl} 
      0 & \eta \text{\ bin not selected}\\ 
      1 & \eta \text{\ bin selected}
      \end{array}\right.\quad,
\end{align}
where $N_{r\in(\eta)}$ is the number of overlapping histograms in the
given $\eta$ bin.  The subscript $v$ indicates that the content
depends on the current vertex bin of event $i$.

\section{Building the final $\dndeta$}

To build the final $\dndeta$ distribution it is enough to sum
\eqref{eq:superhist} and \eqref{eq:overflow} over all interesting
events and then scale it to the number of events that had a collision
trigger
\begin{align}
  \dndetadphi[(\eta,\varphi)] &= \sum_i^{N_{\text{selected}}}
  \dndetadphi[i,v(\eta,\varphi)]\\ 
  I(\eta) &= \sum_i^{N_{\text{selected}}}I_{i,v}(\eta)\\
  \intertext{With $N_{\text{trigger}}$ equal to the number of events with a
    collision trigger, we get}
  \dndeta[(\eta)] &=
  \frac{\Nsel{}}{N_{\text{trigger}}}\frac{1}{\Delta\eta}
  \frac{1}{I(\eta)}\sum_{\varphi}\dndetadphi[(\eta,\varphi)]\quad,
\end{align}
where $\Delta\eta$ is the $\eta$ bin width.  Note that in general
$\Nsel{}\ne N_{\text{selected}}$. 

\section{Using the per--event $\dndetadphi[i,v]$ histogram for other
  analysis} 

\subsection{Multiplicity distribution} 

To build the multiplicity distribution for a given $\eta$ range
$[\eta_1,\eta_2]$, one needs to find the total multiplicity in that
$\eta$ range for each event. To do so, one should sum the
$\dndetadphi[i,v]$ histogram over all $\varphi$ and in the selected
$\eta$ range.
\begin{align}
  n'_{i[\eta_1,\eta_2]}, &= \int_{\eta_1}^{\eta_2}d\eta\int_0^{2\pi}d\varphi
  \dndetadphi[i,v]\quad.\nonumber
\end{align}
However, $n'_i$ is not corrected for the coverage in $\eta$ for the
particular vertex range $v$.  One therefor needs to correct for the
number of missing bins in the range $[\eta_1,\eta_2]$.  Suppose
$[\eta_1,\eta_2]$ covers $N_{[\eta_1,\eta_2]}$ $\eta$ bins, then the acceptance
correction is given by 
\begin{align}
  A_{i,[\eta_1,\eta_2]} = \frac{N_{[\eta_1,\eta_2]}}{\int_{\eta_1}^{\eta_2}d\eta\,
    I_{i,v}(\eta)}\quad.\nonumber
\end{align}
The per--event multiplicity is then given by 
\begin{align}
  n_{i,[\eta_1,\eta_2]} &= A_{i,[\eta_1,\eta_2]}\,n'_{i,[\eta_1,\eta_2]}\nonumber\\
  &= \frac{N_{[\eta_1,\eta_2]}}{\int_{\eta_1}^{\eta_2}\eta
    I_{i,v}(\eta)} \int_{\eta_1}^{\eta_2}d\eta\int_0^{2\pi}d\varphi
  \dndetadphi[i,v]
  \label{eq:event_n}
\end{align}

\subsection{Forward--Backward correlations} 

To do forward--backward correlations, one need to calculate
$n_{i,[\eta_1,\eta_2]}$ as shown in \eqref{eq:event_n} in two bins
$n_{i,[\eta_1,\eta_2]}$ and $n_{i,[-\eta_2,-\eta_1]}$ \textit{e.g.},
$n_{i,f}=n_{i,[-3,-1]}$ and $n_{i,b}=n_{i,[1,3]}$. 

\section{Some results}

\figurename{}s \ref{fig:1} to \ref{fig:3} shows some results.

\begin{figure}[tbp]
  \centering
  \includegraphics[keepaspectratio,width=\textwidth]{%
    dndeta_0900GeV_m10-p10cm_rb05_inel}
  \caption{$\dndeta$ for pp for \INEL{} events at $\sqrt{s}=\GeV{900}$,
    $\cm{-10}\le v_z\le\cm{10}$, rebinned by a factor 5.  Comparisons
    to other measurements shown where applicable}
  \label{fig:1}
\end{figure} 
\begin{figure}[tbp]
  \centering
  \includegraphics[keepaspectratio,width=\textwidth]{%
    dndeta_0900GeV_m10-p10cm_rb05_inelgt0}
  \caption{$\dndeta$ for pp for \INELONE{} events at
    $\sqrt{s}=\GeV{900}$, $\cm{-10}\le v_z\le\cm{10}$, rebinned by a
    factor 5.  Comparisons to other measurements shown where
    applicable}
  \label{fig:2}
\end{figure} 
\begin{figure}[tbp]
  \centering
  \includegraphics[keepaspectratio,width=\textwidth]{%
    dndeta_0900GeV_m10-p10cm_rb05_nsd}
  \caption{$\dndeta$ for pp for \NSD{} events at $\sqrt{s}=\GeV{900}$,
    $\cm{-10}\le v_z\le\cm{10}$, rebinned by a factor 5.  Comparisons
    to other measurements shown where applicable}
  \label{fig:3}
\end{figure} 

\clearpage
\appendix 
\section{Second pass example code}
\lstset{basicstyle=\small\ttfamily,% 
  keywordstyle=\color[rgb]{0.627,0.125,0.941}\bfseries,% 
  identifierstyle=\color[rgb]{0.133,0.545,0.133}\itshape,%
  commentstyle=\color[rgb]{0.698,0.133,0.133},%
  stringstyle=\color[rgb]{0.737,0.561,0.561},
  emph={TH2D,TFile,TTree,AliAODForwardMult},emphstyle=\color{blue},%
  emph={[2]dndeta,sum,norm},emphstyle={[2]\bfseries\underbar},%
  language=c++,%
}
\begin{lstlisting}[caption={Example 2\textsuperscript{nd} pass code to
    do $\dndeta$},label={lst:example},frame=single,captionpos=b]
void Analyse()
{ 
  gSystem->Load("libANALYSIS.so");      // Load analysis libraries
  gSystem->Load("libANALYSISalice.so"); // General ALICE stuff
  gSystem->Load("libPWG2forward2.so");  // New code 

  TH2D*              sum        = 0;                  // Summed hist
  TFile*             file       = TFile::Open("AliAODs.root","READ");
  TTree*             tree       = static_cast<TTree*>(file->Get("aodTree"));
  AliAODForwardMult* mult       = 0;                  // AOD object
  int                nTriggered = 0;                  // # of triggered ev.
  int                nWithVertex= 0;                  // # of ev. w/vertex
  int                nAvailable = tree->GetEntries(); // How many entries
  float              vzLow      = -10;                // Lower ip cut
  float              vzHigh     = +10;                // Upper ip cut
  int                mask       = AliAODForwardMult::kInel;// Trigger mask
  tree->SetBranchAddress("Forward", &forward);        // Set the address

  for (int i = 0; i < nAvailable; i++) { 
    // Read the i'th event 
    tree->GetEntry(i);

    // Create sum histogram on first event - to match binning to input
    if (!sum) 
      sum = static_cast<TH2D*>(mult->GetHistogram()->Clone("d2ndetadphi"));
    
    // Other trigger/event requirements could be defined 
    if (!mult->IsTriggerBits(mask)) continue; 
    nTriggered++;

    // Check if we have vertex 
    if (!mult->HasIpZ()) continue;
    nWithVertex++;

    // Select vertex range (in centimeters) 
    if (!mult->InRange(vzLow, vzHigh) continue; 

    // Add contribution from this event
    sum->Add(&(mult->GetHistogram()));
  }

  // Get acceptance normalisation from underflow bins 
  TH1D* norm   = sum->Projection("norm", 0, 1, "");
  // Project onto eta axis - _ignoring_underflow_bins_!
  TH1D* dndeta = sum->Projection("dndeta", 1, -1, "e");
  // Normalize to the acceptance, and scale by the vertex efficiency 
  dndeta->Divide(norm);
  dndeta->Scale(double(nWithVertex)/nTriggered, "width");
  // And draw the result
  dndeta->Draw();
}
\end{lstlisting}

\begin{thebibliography}{99}
\bibitem{nim:b1:16} Nucl.Instrum.Meth.B1:16
\bibitem{phyrev:a28:615} Phys.Rev.A28:615
\end{thebibliography}
\end{document}