#  On the Boundedness solutions of the difference equation
$x_{n+1}=ax^{\alpha}_{n}+bx^{\alpha}_{n-1},0<\alpha\leq 2$ and its application
in medicine

Zeraoulia Rafik
University of batna2.Algeria
Departement of mathematics
yabous ,khenchela
<EMAIL_ADDRESS>
&Alvaro humberto Salas
Universidad national de Colombia
departement of physics
Bogota Colombia
<EMAIL_ADDRESS>
&Lorenzo Martinez
Universidad national de Colombia
departement of mathematics
Manizales,Caldas
<EMAIL_ADDRESS>

###### Abstract

Recently ,mathematicians have been interested in studying the theory of
discrete dynamical system, specifically difference equation, such that
considerable works about discussing the behavior properties of its solutions
(boundedness and unboundedness) are discussed and published in many areas of
mathematics which involves several interesting results and applications in
applied mathematics and physics ,One of the most important discrete dynamics
which is become of interest for researchers in the field is the rational
dynamical system .In this paper we may discuss qualitative behavior and
properties of the difference equation $x_{n+1}=ax^{2}_{n}+bx^{2}_{n-1}$ with
$a$ and $b$ are two parameters and we shall show its application to medicine

_K_ eywords difference equation $\cdot$ boundedness $\cdot$ number theory
$\cdot$ dynamical system

##  1 Introduction

The theory of difference equations finds many applications in almost all areas
of natural science [Daniel J. Duffy(2006)]. increasingly clearly emerges the
fundamental role that difference equations with discrete and continuous
argument is played for understanding nonlinear dynamics and phenomena also it
is used for combinatorics and in the approximation of solutions of partial
differential equations [Y. Ordokhan1 ,S. Davaei far(2017)]. The increased
interest in difference equations is partly due to their ease of handling. A
minimum is enough computing and graphical tools to see how the solution of
difference equations trace their bifurcations with changing parameters [Josef
Diblik ,Miroslava Ruzickova , Barbora Vaclavikova (2008)]. Thus opens a
complex understanding as well invariant manifolds for linear and nonlinear
dynamical systems.nonlinear difference equations and systems are of wide
interest due to their applications in real life. Such equations appear
naturally as the mathematical models which describe biological, physical and
economical phenomena.

Although difference equations have very simple forms, however, it is extremely
difficult to understand completely the global behavior of their solutions.y
the global behaviors of their solutions. One can refer to [Charlie y
Routh(1992)],[Camouzis and G.Ladas(2008)],[E.A.Grove and G. Ladas(2005)] and
the references therein. Difference equations have always played an important
role in the construction and analysis of mathematical models of biology,
ecology, physics, economic processes, etc. The study of nonlinear rational
difference equations of higher order is of paramount importance, since we
still know so little about such equations. In this paper [A. Q. Khan,S. M.
Qureshi (2020)] A. Q. Khan and S. M. Qureshi discussed Dynamical properties of
some rational systems of difference equations such as they explored the
equilibrium points, local and global dynamics, rate of convergence,
instability and boundedness of positive solution of some rational systems of
difference equations .As application to modern science ,namely, in
mathematical biology they also explored the local dynamics about equilibrium
points of the discrete-time Levin’s model .In the meanwhile A. Q. Khan studied
and discussed Global dynamics of a $3\times 6$ exponential system of
difference equations defined as:

$\begin{cases}x_{n+1}=\frac{\alpha_{1}+\beta_{1}\exp(-x_{n})}{\gamma_{1}+y_{n-1}}\\\
y_{n+1}=\frac{\alpha_{2}+\beta_{2}\exp(-y_{n})}{\gamma_{2}+z_{n-1}}\\\
z_{n+1}=\frac{\alpha_{3}+\beta_{3}\exp(-z_{n})}{\gamma_{3}+y_{n-1}}\end{cases}$

$n=0,1,....$

where parameters $\alpha_{i},\beta_{i},\gamma_{i}(i=1,2,3)$ and initial
conditions $x_{i},y_{i},z_{i}(i=0,-1)$ are nonnegative real numbers.

In [R. Abo-Zeid (2017)] .R. Abo Zeid has discussed the global behavior of all
solutions of the difference equation:

$x_{n+1}=\frac{x_{n}x_{n-1}}{ax_{n}+bx_{n-1}},n\in\mathbb{N}_{\textpdfrender{TextRenderingMode=Stroke,LineWidth=.1pt,}{0}}$
(1)

where $a,b$ are real numbers and the initial conditions $x_{-1},x_{0}$ are
real numbers .In this paper, we discuss the global behavior of the difference
equation :

$x_{n+1}=Ax_{n}^{\alpha}+Bx_{n-1}^{\alpha}$ (2)

with :$A,B$ are two parameters real numbers and $\alpha$ is a real number such
that :

$0<\alpha\leq 2$, For $\alpha=1$ the dynamics defined in (2) is deduced from
(1) by the substitution $y=\frac{1}{x_{n}}$ , the global behavior of all
solutions of (1) is discussed as well in [R. Abo-Zeid (2017)] by R. Abo-Zeid
,The difference equation (2) for $\alpha=1$ ,namely ,

$y_{n+1}=by_{n}+ay_{n-1},n\in\mathbb{N}$ (3)

The characteristic equation of equation (3) is :

$\lambda^{2}-b\lambda-a=0$ (4)

Equation (4) has two roots :
$\lambda_{1}=\frac{b-\sqrt{b^{2}-4ac}}{2a},\lambda_{2}=\frac{b+\sqrt{b^{2}-4ac}}{2a}$
The form of the solution should be according to the value of the quantity
$b^{2}+a^{2}$,The following theorem [S.Elaydi(2005)] is useful in studying the
solutions of the difference equation (3)

###### Theorem 1.1

The following statements holds:

  * •

1) All solutions of (3) oscillates (about zero) if and only if the
characteristic equation has no positive roots.

  * •

2) All solutions of (3) converge to zero if and only if
$\max\\{\lambda_{1},\lambda_{2},\\}<1$

For boundedness solutions of (3) one can refer to [R. Abo-Zeid (2017)] .Now
for $\alpha=2$ which it is the aim of this paper we are ready to do some
analysis and discussion about the global behavior of solutions of this
difference equation:

$y_{n+1}=by_{n}^{2}+ay_{n-1}^{2},n\in\mathbb{N}$ (5)

##  2 Analysis and discussion:

Case:1 $|a|,|b|<1$ ,for this case we may use a nice trick just we assume $a,b$
are two bounded functions such that :
$a=\sin(\theta),b=\cos(\theta),\theta\in\mathbb{R}$ then the dynamics defined
in (5) become :

$y_{n+1}=\cos(t)y_{n}^{2}+\sin(t)y_{n-1}^{2},n\in\mathbb{N},t\in\mathbb{R}$
(6)

Now ,we may ask for which values of $\theta$ does this equation
$x_{n+1}=\cos(\theta)x^{2}_{n}+\sin(\theta)x^{2}_{n-1}$ have bounded
solutions?.

We ran a small computation: The below plot,see figure (1) is created as
follows: For each point $r(\cos\theta,\sin\theta)$, We use $x_{-1}=0$,
$x_{0}=r$ as initial values, and $\theta$ as the parameter.

The white area is where the iterations is less than $2$ (bounded), the first
30 iterations. As we see, a small circle around the origin is white, meaning
that there are small initial values that is bounded for every $\theta$. Now,
this is not a proof, but the picture suggests this.

Changing the cut-off to $40$ iterations does not change the picture much.

Figure 1: Bounded solution for $x_{-1}=0,r=x_{0}$

For analytical proof we have :Let $r=1/\sqrt{2}$. If $x_{-1},x_{0}\leq r$,
then the sequence is bounded for all $\theta$.

We note that: $|\cos(\theta)x_{n-1}^{2}+\sin(\theta)x_{n}^{2}|\leq
r^{2}(|\cos\theta|+|\sin\theta|)\leq r^{2}\sqrt{2}\leq 1/\sqrt{2}$, and the
statement follows from induction.

Case:2

$|a|,|b|=1$, for the case $|a|=|b|=1$ we note after runing with the same
initial conditions for each point $r(\cos\theta,\sin\theta)$ (iteration of the
plot (20 iteration) ) that the obtained picture changed ( Rotation of picture)
and the white region become larger than it for the precedent case ,namely,
Case:1 this mean that we have another new initial values that is bounded for
every $\theta$.(The number of initial values is increasing).see figure 2

Figure 2: new intial values that bounded solutions for every
$\theta,x_{-1}=0,r=x_{0}$

Case3 In this case we may assume
:$a^{\prime}=|a|k,b^{\prime}=|b|k,k\in\mathbb{R}$ and
$|a^{\prime}|,|b^{\prime}|>1$, then our dynamics become :

$y_{n+1}=|\cos(t)|y_{n}^{2}+|\sin(t)|y_{n-1}^{2},n\in\mathbb{N},t\in\mathbb{R}$
(7)

With the same data and conditions that we have used in both of case1 and case
2, we noticed after runing using mathematica the white region become smaller
whenever $k$ being larger this indicate that the dyanmics defined in (7) would
be unbounded,for numerical evidence we just take , $k=5,20,-15$ as shown in
below plots (figure (3+4+5)).

Figure 3: boundedness solutions for the dynamic (7),$k=5$ Figure 4:
boundedness solutions for the dynamic (7),$k=15$ Figure 5: boundedness
solutions for the dynamic (7),$k=-150$

We noticed For the remainder case when $0<\alpha\leq 1$ the iteration of the
dynamics (2) show us much intinial values that is bounded for every $\theta$ ,
we took $\alpha=\frac{1}{2},\frac{1}{3},\alpha=1$ as comparative examples as
shown in below figures:

Figure 6: bounded solutions for $\alpha=1$ (white circle arround origin)

Figure 7: bounded solution for $\alpha=\frac{1}{2},\frac{1}{3},\cdots$ (The
circle become white at all this indicate all solution of our dynamics are
bounded for any arbitary initial values )

##  3 Stability and analysis

Dynamical modeling is the study of change and changes take place everywhere in
life. As a result dynamical systems have a wide range of application areas in
applied science and engineering. With these systems, real life situations can
be turned into the language of mathematics. If one can create a good model for
a real life situation, we will be able to predict the future states of the
system by simply iterations according to this model.

Stability, like one of the most important concepts in Discrete Dynamical
Theory, can tell much about the behavior of the dynamical system. In this
section a symbolic Mathematica package for analysis and control of chaos in
discrete two dimensional dynamical nonlinear systems, is presented. There are
constructed some computer codes to find stability types of the fixed points,
covering the stability of the one-dimensional nonlinear dynamical systems.
Applications are taken from chemical kinetics and population dynamics
(logistic model).Since our dynamics is two dimonsional discret dynamical
system it is enough to run the below code to get satibility and oscillation
points arround origing…

example one Let us apply that code for our dynamics: let $a=2,b=9$

        twoDimStab[2 x^2 + y + 1, 9x^2 - 1]

           1     112
    {Null, {--, -(---)} -> oscillatory source}
            11    121

One can do more examples just applying the above code with change of variable
$a$ and $b$.

For the periodicity of our dynamic for $\alpha=2$ it is already discussed in
[R. Abo-Zeid (2017)],The dynamics is two period solution

##  4 PHASE PLANE DIAGRAMS

Continuous systems are often approximated as discrete processes, meaning that
we look only at the solutions for positive integer inputs. Difference
equations are recurrence relations, and first order difference equations only
depend on the previous value. Using difference equations, we can model
discrete dynamical systems. The observations we can determine, by analyzing
phase plane diagrams of difference equations, are if we are modeling decay or
growth, convergence, stability, and equilibrium points.In this section we may
give a diagram plot of our dynamics for some values $a$ and $b$ and $\alpha=2$
, The dynamics defined in (2) can be written as sytem of two difference
equations :

$\begin{cases}x(t+1)=ax_{t}^{2}+y(t)+1\\\ y(t+1)=bx_{t}^{2}-1\end{cases}$ (8)

The phase portrait plot of the dynamic (8) is shown in below figures, Here is
the mathematica code one can change values of a,b to analyse the slop field

Figure 8: diagram plot of nonlinare dynamics (8), for $a=0.5,b=2$

##  5 Time series analysis and modeling (application to medicin)

Nowadays mathematics is being successfully applied to a number of important
fields in medicine including biofluids, cardiovascular diseases, clinical
schedules and tests, data analysis, drug design and discovery, epidemiology,
genetics, image processing, immunology, instrumentation, microbiology [José M.
Amigó , Michael Small (2017)], neuroscience, oncology, virology and more. The
list of tools includes virtually the whole of applied mathematics. To cite the
most familiar ones: difference equations and discrete-time dynamical systems,
information and coding theory, graph and network theory, integral transforms,
numerical and computational mathematics, ordinary differential equations and
continuous-time dynamical systems, partial differential equations, stochastic
and time-delay differential equations, statistics, probability and time-series
analysis. All this research has contributed to and continues to increasingly
contribute both to better understand medical phenomena and to finding
practical ways of action.

Time series modeling for predictive purpose has been an active research area
of machine learning for many years [Fatoumata Dama and Christine Sinoquet
(2019)]. However, no sufficiently comprehensive and meanwhile substantive
survey was offered so far. This survey strives to meet this need. A unified
presentation has been adopted for entire parts of this compilation. Time
series data are amongst the most ubiquitous data types that capture
information and record activity in most aeras. In any domain involving
temporal measurements via sensors, censuses, transaction records, the capture
of a sequence of observations indexed by time stamps first allows to provide
insights on the past evolution of some measurable quantity. Beyond this goal,
the pervasiveness of time series has generated an increasing demand for
performing various tasks on time series data (visualization, discovery of
recurrent patterns, correlation discovery, classification, clustering, outlier
detection, segmentation, forecasting, data simulation).

Time series classification (TSC) is one of data minings persistent challenges.
Applications of TSC abound in fields including agriculture, medicine, and
engine prognostics [Abdoli, A., Murillo, A. C., Yeh, C. C. M., Gerry, A. C.,
and Keogh, E. J. (2018, December)],and [Yu, W., Kim, I. Y., Mechefske, C.
(2021)] a common goal being to detect instances of sub optimal behavior or
decreased health (biologically or mechanically) as just one real-world
example. Dozens of new TSC algorithms were introduced in the last four years
alone [Bagnall, A., Lines, J., Bostrom, A., Large, J., Keogh, E.(2017)]. This
trend has been intensified by the increasing availability of real-world
datasets. In fact, the classification of any inherently ordered data
(temporally or otherwise) can be treated as a TSC problem [Gamboa, J. C. B.
(2017)] , making for a vast breadth of real-world applications. Deep learning
methods have shown suitability for time series classification in the health
and medical domain, with promising results for electrocardiogram data
classification. Successful identification of myocardial infarction holds life
saving potential and any meaningful improvement upon deep learning models in
this area is of great interest.In this section we show that the discussed
dynamics ,namely , (8) present new classes of heartbeat ,namely, A new
dynamical model for Generating Synthetic Electrocardiogram Signals,The
electrocardiogram (ECG) is a time-varying signal reflecting the ionic current
flow which causes the cardiac fibers to contract and subsequently relax. The
surface ECG is obtained by recording the potential difference between two
electrodes placed on the surface of the skin. A single normal cycle of the ECG
represents the successive atrial depolarization/repolarization and ventricular
depolarization/repolarization which occurs with every heartbeat. These can be
approximately associated with the peaks and troughs of the ECG waveform
labeled $P,Q,R,S$, and $T$ as shown in figure 9

Figure 9: Morphology of a mean PQRST-complex of an ECG recorded from a normal
human

Extracting useful clinical information from the real (noisy) ECG requires
reliable signal processing techniques [AL Goldberger,E Goldberger (1977)].
These include R-peak detection [Jiapu Pan; Willis J. Tompkins (1985)], [C. Li,
C. Zheng, and C. Tai(1995)], QT-interval detection and the derivation of heart
rate and respiration rate from the ECG [P.E. McSharry , G.D. Clifford, L.
Tarassenko, L.A. Smith (2003)]. The RR-interval is the time between successive
R-peaks, the inverse of this time interval gives the instantaneous heart rate.
A series of RR-intervals is known as a RR tachogram and variability of these
RR-intervals reveals important information about the physiological state of
the ECG signal. The ECG may be divided into the following sections.

  * •

Q-wave: A small low-voltage deflection away from the baseline caused by the
depolarization of the atria prior to atrial contraction as the activation
(depolarization) wavefront propagates from the SA node through the atria.

  * •

PQ-interval: The time between the beginning of atrial depolarization and the
beginning of ventricular depolarization.

  * •

QRS-complex: The largest-amplitude portion of the ECG, caused by currents
generated when the ventricles depolarize prior to their contraction. Although
atrial repolarization occurs before ventricular depolarization, the latter
waveform (i.e. the QRS-complex) is of much greater amplitude and atrial
repolarization is therefore not seen on the ECG.

  * •

QT-interval: The time between the onset of ventricular depolarization and the
end of ventricular repolarization. Clinical studies have demonstrated that the
QT-interval increases linearly as the RR-interval increases . Prolonged QT-
interval may be associated with delayed ventricular repolarization which may
cause ventricular tachyarrhythmias leading to sudden cardiac death

  * •

ST-interval: The time between the end of S-wave and the beginning of T-wave.
Significantly elevated or depressed amplitudes away from the baseline are
often associated with cardiac illness.

  * •

T-wave: Ventricular repolarization, whereby the cardiac muscle is prepared for
the next cycle of the ECG.

Analysis of variations in the instantaneous heart rate time series using the
beat-to-beat RR-intervals (the RR tachogram) is known as HRV analysis [M Malik
, AJ camm(1995)], [Task force of the european society of cardiology(1996)].
HRV analysis has been shown to provide an assessment of cardiovascular disease
[M H Crawford,S Bernstein and P Deedwania(1999)].A dynamical model based on
three coupled ordinary differential equations is introduced which is capable
of generating realistic synthetic electrocardiogram (ECG) signals.The
dynamical equations of motion are given by a set of three ordinary
differential equations defined as the following :

$\begin{cases}\dot{x}=\alpha x-\omega y\\\ \dot{y}=\alpha y+\omega x\\\
\dot{z}=-\sum_{i\in\\{P,Q,R,S,T\\}}\exp\bigg{(}-\frac{-\Delta\theta_{i}^{2}}{2b_{i}^{2}}\bigg{)}-(z-z_{0})\end{cases}$
(9)

where $\alpha=1-\sqrt{x^{2}+y^{2}},\Delta\theta_{i}=(\theta-\theta_{i})\bmod
2\pi,\theta=\arctan 2(y,x)$, (the four quadrant arctangent of the real parts
of elements of $x$ and $y$ , with $-\pi\leq\arctan 2(y,x)\leq\pi$) and
$\omega$ is the angular velocity of the trajectory as it moves around the
limit cycle. Baseline wander was introduced by coupling the baseline value
$z_{0}$ in (9) to the respiratory frequency $f_{2}$ using
:$z_{0}(t)=A\sin(2\pi f_{2}t)$ where $A=0.15mV$

Figure 10:  ECG generated by dynamical model: (a) 10 s and (b) 50 s Figure 11:
Comparison between (a) synthetic ECG with additive normally distributed
measurement errors and (b) real ECG signal from a normal human.

##  6 Introduction to Electrocardiography

An electrocardiogram is a recording of the electrical activity of the
heart.The heart can be viewed as a three-dimensional vector. Therefore, its
electrical activity can, in theory, be recorded by three orthogonal leads. In
practice, however, a standard clinical EKG is recorded with 12 leads: 6 limb
leads and 6 precordical leads.A normal EKG reflects the electrical activity of
each of the four heart chambers: left and right ventricles and left and right
atria (Fig. 1).

Figure 12: A normal EKG reflects the electrical activity of each of the four
heart chambers: left and right ventricles and left and right atria

The P wave marks the beginning of atrial depolarization. The onset of the Q
wave is produced by the stimulus from the Purkinje system. The depolarization
of the ventricle produces the R wave. The S wave is initiated when the
excitation reaches the base of the ventricle. The T wave is produced by
ventricular repolarization (Tompkins 1993; Goldman 1992). The synthetic ECG
(Figure 10) illustrates the modulation of the QRS-complex due to RSA and Mayer
waves. Observational uncertainty is incorporated by adding normally
distributed measurement errors with mean zero and standard deviation $0.025mV$
[Figure 11], yielding a similar signal to a segment of real ECG from a normal
human [Figure 11]. Electrocardiogram (ECG) detection is currently the most
effective and direct way to detect ECG signals [Z. Tang, G. Zhao, and T.
Ouyang(2021)]. At present, the diagnosis of cardiac diseases is mainly
determined by medical doctors and clinicians through manual detection and ECG
analysis.ECG is a diagnostic technology that records the electrocardiography
activities of the heart in a certain time unit through the chest of biological
objects. It collects and records the electrodes connected to the skin of
specific parts of biological objects and preserves the relevant contents in a
certain form [D. A. Winter, P. M. Rautaharju, and H. K. Wolf(1997)].

The pre-ejection-period (PEP) is the time span between the depolarization of
the left ventricle (R onset) and opening of the aortic valve (B point). The R
onset is the beginning of the Q-wave signal; it indicates the beginning of the
depolarization and can be picked up from the ECG signal. As signature for the
beginning of the Q wave, we take the minimum of the ECG’s second derivative.
Its peak indicates the maximum curvature at the transition of the ECG signal
into the Q wave. However, as the Q wave is relatively small, other signatures
in the ECG signal can be misinterpreted as the R onset in an automated
evaluation of the data. In general, first and higher-order derivatives of
noisy signals suffer from containing spurious peaks. We therefore restrict the
possible occurrences of R onset to a time window after the, easily
identifiable, R peak (peak of the QRS complex). Within that window, the R
onset is typically seen as a clear negative peak of the ECG signal’s second
derivative, which can be located reliably and with high precision and thus
allows a reliable identification of the Q wave onset. Heart rate (HR) is then
calculated from the time difference of subsequent R points.[M. Thomas, M. K.
Das, and S. Ari(2015)]

The time point for the opening of the aortic valve (B point) is derived from
the impedance cardiogram (ICG). The impedance Z, and thus the ICG, is
sensitive to a variation of blood volume in the thorax. The first derivative,
$\frac{dZ}{dt}$, corresponds to blood flow. The second derivative
$\frac{d^{2}Z}{dt^{2}}$, in turn, corresponds to a change of the blood flow
and is thus indicative for the opening of the heart valves. The B point is the
onset of the aortic valve’s opening, indicated by a negative peak in the third
derivative, $\frac{d^{3}Z}{dt^{3}}$. While, compared to ECG, the ICG signal is
smooth and devoid of characteristic spikes, its first, second, and third
derivative show distinct features. As selection criterion for picking the
correct peak of the third derivative, we use the, easily identifiable, peak of
the first derivative, $\frac{dZ}{dt}$. The B point is obtained as the minimum
of $\frac{d^{3}Z}{dt^{3}}$ that occurs just before the maximum in
$\frac{dZ}{dt}$. This strategy allows for an automated evaluation of the PEP
interval for the large data sets, with few outliers and the required
precision. To calculate the derivatives of the measured signals, we use the
Savitzky-Golay filter [Jianwen LuoKui ,YingLijing BaiLijing Bai(2005)]. This
method allows data smoothing, while keeping intact signatures like peaks, and
the simultaneous determination of derivatives. Similar to a moving average, a
moving section of the data is selected. However, instead of a simple
averaging, the algorithm fits a polynomial of given degree to the selected
sequence. Then, one point of the fitted polynomial (usually the central point)
is taken as value for the smoothed curve. Higher derivatives are taken from
the corresponding derivatives of the fitted polynomial at the respective
point. The Savitzky-Golay filter is implemented numerically by a list
convolution with a kernel. That kernel is calculated in advance for the number
of points for the moving fitting, the order of the polynomial, and the order
of the derivative. We use a kernel length of 100 points, corresponding to a
time interval of 50 ms, and a 3rd-order polynomial for all kernels and for the
ICG and ECG signals. The third derivative of the ICG signal is calculated from
the first derivative of the ICG signal, which, together with Z, is provided by
the Biopac MP36 system (i.e., the system we used to measure ICG/ECG). We
ensure that no time lag gets introduced between the ICG and ECG signals and
their derivatives by the Savitzky-Golay filter. Thus, PEP and LVET data get
extracted from the ICG and ECG measurements in a semi-automated way and with a
by-heartbeat resolution. The Mathematica Notebook output is stored in a text
file, with every row containing a timestamp, the corresponding length of
cardiac PEP, LVET, and HR for each heartbeat. Here is the Graphical display of
$Z,\frac{dZ}{dt}$,EKG

Figure 13: Graphical display of $Z,\frac{dZ}{dt}$ ,EKG , Istart=1000,Iend
=16000

##  7 Medical interpretation of our discret dynamics

To analyze ECG signal, the most necessary step is to extract its QRS wave
group. The QRS complex reflects changes in depolarization potential and time
of the left and right ventricles. Considering the robustness and stability,
For analysis of ECG signal using delay differential equation (DDE’s) is
computationally fast one can refer to the system of ODE in (9), and could be
the basis for a real time diagnostic system. DDEs reveal non-linear properties
as well as spectral properties of the data as shown in this paper [Claudia
Lainscsek and Terrence J. Sejnowski(2013)] by Claudia Lainscsek and Terrence
J. Sejnowski they analyzed ECG signal using delay differential equation which
leads to the good classification of Electrocardiogram ,may that classification
wouldn’t work as well using discret map ,some authors discussed
Electrocardiogram Signal Classification in the Diagnosis of Heart Disease
Based on RBF Neural Network [Yan Fang , Jianshe Shi, Yifeng Huang, Taisheng
Zeng,Yuguang Ye , Lianta Su, Daxin Zhu,and Jianlong Huang(2022)] such that
they extracted QRS wave using discret map (they used difference equations to
analyze ECG ),Recall that in our paper we are interested to the behavior of
the following dynamics :

$x_{n+1}=ax^{2}_{n}-bx^{2}_{n-1},n=0,1,....$

,$a,b$ two real parameters ,It is known that ECG signal has a very obvious
periodicity, Taking $T$ as sampling period we may rewrite our dynamics using
$T$ as :

$x((n+1)T)=ax^{2}(nT)-bx^{2}(T(n-1),n=0,1,...$ (10)

The high-frequency characteristics can be enhanced by a nonlinear square
function [Yan Fang , Jianshe Shi, Yifeng Huang, Taisheng Zeng,Yuguang Ye ,
Lianta Su, Daxin Zhu,and Jianlong Huang(2022)],see page 3 equation 4, whose
equation can be expressed as :

$y((nT)=x^{2}(nT)$ (11)

This means strongly that our discret dynamics which is defined in (10) can be
interpreted in the medical point of view as :Enhanced high frequencies
regarding the behavior of the heartbeat in short time, we may consider the
discret time defined in (10) as a discretized boundary value problem (delay
differential equation) using some standrad numerical methods like finite
difference method in 1D (dimension) , we may give a short analysis of the ECG
signal which is produced by the dynamics defined in (10), Assume the
correspending non linear boundary value problem which we wanted to solve
satisfies $x_{1}=0$ and $x_{n+1}=1$ using the Finit difference technique [R.
P. AGARWAL and Y. M. CHOW(1985)]. The first step is to partition the domain
$[0,1]$ into a number of sub-domains or intervals of length $h$. So, if the
number of intervals is equal to $n$, then $nh=1$. We denote by $x_{i}$ the
interval end points or nodes, . In general, we have $x_{i}=(i-1)h$,
$i=1,2,\cdots n+1$. Let us denote the concentration at the ith node by
$C_{i}$, for short enough time (say between $t_{i},t_{i+1}$) the dynamics (10)
become as a simple linear differential equation such that : for $h\to 0$ we
have $x((n+1)T)\to\dot{x}\leq 0$ implies that we have a decreasing frequencies
coming up to a constant signal $x(t)$ thus ,the dynamic (10) can be
interpreted as heart attack or heart failure ,In general that case indicate
Cardiac diseases ,(see figure in example 1), Enhanced high frequency appear
whenever $\dot{x}>0$ which means derivative of the analyzed signal is always
positive (one can refer to figures in examples (2+3+4)).It is hard to analyze
ECG signal using coupled signal which is defined in (12)

we have plotted many figures with some values of $a$ and $b$ such that we call
for medical interpretation $a$ is a factor of life (Electrocardiogram (EKG))
,identification of that factor(EKG) present attempt to improve patient
survival and $b$ may present blood loss which causes the cardiac arrest.we may
introduce a new parameter here $\sigma$ as a toxine factor this for good
modelization of the discussed phenomena. Now,let us try showing heartbeat
classes ,The dynamics (8) for $\alpha=2$ can be written as :

$\begin{cases}x(t+1)=ax_{t}^{2}+\sigma y(t)+1\\\
y(t+1)=bx_{t}^{2}-1\end{cases}$ (12)

We have noticed that For heartbeat classes plots we should have : $a\leq
0.15$,we may try playing with $b$ and $\sigma<0$ (should be negative) and for
$x_{0}=0.07,y_{0}=0.08$ initial values ,increasing values of $\sigma$ means
attempt to reduce toxine factor and increasing values of $b$ means raise
factor of EKG thus a good attempt to improve patient survival ,for the reverse
assumption we would have cardiac arrest record.

Example1

we start with the first heartplot , $a=0.15,\sigma=-0.45,b=0.45$

Figure 14: heartbeat classe for
$a=0.15,\sigma=-0.45,b=0.45,x_{0}=0.07,y_{0}=0.08$

Example2 let :$a=0.15,\sigma=-0.45,b=0.45$

Figure 15: heartbeat case for $a=0.15,\sigma=-0.6,b=0.5,x_{0}=0.07,y_{0}=0.08$

Example3 let :$a=0.15,\sigma=-0.65,b=0.58$

Figure 16: heartbeat case for
$a=0.15,\sigma=-0.65,b=0.58,x_{0}=0.07,y_{0}=0.08$

Example4 let :$a=0.15,\sigma=-0.65,b=0.57$

Figure 17: heartbeat case for
$a=0.15,\sigma=-0.75,b=0.6,x_{0}=0.07,y_{0}=0.08$

##  8 Conclusion:

This study of difference equations with public health applications to science
develops the methodology for the solution of the general kth order linear
difference equation using the generating function approach and computers tools
like mathematica. It includes an examination of the dynamics of disease spread
and containment in populations using illness-death models ,Cardiovascular
disease is one of the major hazard to human health today. ECG stands for
electrocardiogram and it is an important way for clinical diagnosis
cardiovascular disease. The ECG refers to the small voltages ( 1mv) found on
the skin as a result of electrical activity of the heart. These electrical
actions trigger various electrical and muscular activity in the heart. The
health and function of the heart can be measured by the shape of the ECG
waveform. Typical heart problems are leaking valves and blocked coronary
arteries.in our paper we have discussed a new discret dynamics which lead to
discover a new illness-death model(Enhanced high frequencies model ) using
time series diagram and ECG analyses with some data which are defined by
iterative discret dynamics such as heartbeat classes ,in particulary attempts
to improve patient for being survival,we may call it model of life for more
time.

##  9 Data Availability

The data that support the findings of this study can be obtained from the
corresponding author upon reasonable request.

##  10 Conflicts of Interest

The author declare that they have no conflicts of interest.

##  11 Authors’ Contributions

All authors contributed to the study conception and design of that research .
Material preparation, data collection, analysis and medical interpretation
were performed by Zeraoulia Rafik . The first draft of the manuscript was
written by Zeraoulia Rafik and all authors commented on previous versions of
the manuscript. All authors read and approved the revised manuscript.

## References

  * [Daniel J. Duffy(2006)] Daniel J. Duffy, 2006.Finite Difference Methods in Financial Engineering: A Partial Differential Equation Approach, ISBN:=9780470858820, 9780470858820
  * [Y. Ordokhan1 ,S. Davaei far(2017) ] Able, B., 1956.Approximate Solutions of Differential Equations by Using the Bernstein Polynomials
DOI: 10.5402/2011/787694

  * [Josef Diblik ,Miroslava Ruzickova , Barbora Vaclavikova (2008)] Josef Diblik ,Miroslava Ruzickova , Barbora Vaclavikova, 2008. Bounded Solutions of Dynamic Equations on Time Scales,International Journal of Difference Equations (IJDE). ,ISSN 0973-6069 Volume 3 Number 1 (2008), pp. 61–69
  * [Charlie y Routh(1992)] Agarwal,Difference equations and inequalities,1992.first edition, Marcel Dekker
  * [ Camouzis and G.Ladas(2008)] Dynamics of Third-Order Rational Difference Equations;With Open Problems and Conjectures, 2008. Chapman and Hall/HRC Boca Raton,
  * [E.A.Grove and G. Ladas(2005)] E.A. Grove and G. Ladas,Periodicities in Nonlinear Difference Equations,2005. Chapmanand Hall/CRC, 2005
  * [R. Abo-Zeid (2017)] R. Abo-Zeid On the solutions of a second order difference equation,2017.Mathematica Moravica Vol. 21, No. 2 (61-73
  * [S.Elaydi(2005)] S.Elaydi, An Introduction to Difference Equations,2005. Third Edition, Springer, New York.
  * [Fatoumata Dama and Christine Sinoquet (2019)] Fatoumata Dama and Christine Sinoquet,2019. Time Series Analysis and Modeling to Forecast: a Survey .LS2N / UMR CNRS 6004, Nantes University, France
  * [Abdoli, A., Murillo, A. C., Yeh, C. C. M., Gerry, A. C., and Keogh, E. J. (2018, December)] Abdoli, A., Murillo, A. C., Yeh, C. C. M., Gerry, A. C., Keogh, E. J. 2018. Time series classification to improve poultry welfare. In 2018 17TH IEEE International conference on machine learning and applications (ICMLA) (pp. 635, 642). IEEE
  * [Yu, W., Kim, I. Y., Mechefske, C. (2021)] Yu, W., Kim, I. Y., Mechefske, C. ,2021. Analysis of different RNN autoencoder variants for time series classification and machine prognostics. Mechanical Systems and Signal Processing, 149, 107322
  * [Bagnall, A., Lines, J., Bostrom, A., Large, J., Keogh, E.(2017)] Bagnall, A., Lines, J., Bostrom, A., Large, J., Keogh, E. ,2017. The great time series classification bake off: a review and experimental evaluation of recent algorithmic advances. Data mining and knowledge discovery, 31(3), 606-660
  * [Gamboa, J. C. B. (2017)] Gamboa, J. C. B. 2017. Deep learning for time series analysis.
  * [A. Q. Khan,S. M. Qureshi (2020)] A. Q. Khan,S. M. Qureshi (2020) . Dynamical properties of some rational systems of difference equations.Mathematical method in the applied science
  * [José M. Amigó , Michael Small (2017)] A. Q. Khan,S. M. Qureshi (2017) .Mathematical methods in medicine: neuroscience, cardiology and pathology National library of medcin
  * [AL Goldberger,E Goldberger (1977)] AL Goldberger,E Goldberger (1977) .Clinical electrocardiac ,Louis MO
  * [Jiapu Pan; Willis J. Tompkins (1985)] Jiapu Pan; Willis J. Tompkins (1985).A Real-Time QRS Detection Algorithm
  * [P.E. McSharry , G.D. Clifford, L. Tarassenko, L.A. Smith (2003)] P.E. McSharry , G.D. Clifford, L. Tarassenko, L.A. Smith (2003).A dynamical model for generating synthetic electrocardiogram signals .IEEE Transactions on Biomedical Engineering ( Volume: 50, Issue: 3)
  * [M Malik , AJ camm(1995)] M Malik , AJ camm(1995). Heart rate variability ,Armonk,NY :futura
  * [Task force of the european society of cardiology(1996)] Task force of the european society of cardiology(1996). Heart rate variability :standard of measurement physiological interpretation and clinical use the north americain society of pacing and electrophysiology,sophia antipolis,France
  * [M H Crawford,S Bernstein and P Deedwania(1999)] M H Crawford,S Bernstein and P Deedwania(1999). ACC AHA , guidlines for ambulatory ,electrocardiography.circulation vol 100, pp 886-893
  * [Jianwen LuoKui ,YingLijing BaiLijing Bai(2005)] Jianwen LuoKui ,YingLijing BaiLijing Bai(2005). Savitzky–Golay smoothing and differentiation filter for even number data
  * [C. Li, C. Zheng, and C. Tai(1995)] C. Li, C. Zheng, and C. Tai(1995).Detection of ECG characteristic points using wavelet transforms. IEEE Transactions on Biomedical Engineering, vol. 42, no. 1, pp. 21–28
  * [Z. Tang, G. Zhao, and T. Ouyang(2021)] Z. Tang, G. Zhao, and T. Ouyang(2021).Two-phase deep learning model for short-term wind direction forecasting,Renewable Energy, vol. 173, pp. 1005–1016
  * [D. A. Winter, P. M. Rautaharju, and H. K. Wolf(1997)] , D. A. Winter, P. M. Rautaharju, and H. K. Wolf(1997).Measurement and characteristics of over-all noise content in exercise electrocardiograms,American Heart Journal, vol. 74, no. 3, pp. 324–331, 1967.
  * [M. Thomas, M. K. Das, and S. Ari(2015)] M. Thomas, M. K. Das, and S. Ari(2015). Automatic ECG arrhythmia classification using dual tree complex wavelet based features, AEU - International Journal of Electronics and Communications, vol. 69, no. 4, pp. 715–721, 2015
  * [Yan Fang , Jianshe Shi, Yifeng Huang, Taisheng Zeng,Yuguang Ye , Lianta Su, Daxin Zhu,and Jianlong Huang(2022)] Yan Fang , Jianshe Shi, Yifeng Huang, Taisheng Zeng,Yuguang Ye , Lianta Su, Daxin Zhu,and Jianlong Huang(2022).Electrocardiogram Signal Classification in the Diagnosis of Heart Disease Based on RBF Neural Network .Hindawi Computational and Mathematical Methods in Medicine Volume 2022, Article ID 9251225, 9 pages
  * [Claudia Lainscsek and Terrence J. Sejnowski(2013)] Claudia Lainscsek and Terrence J. Sejnowski(2013).Electrocardiogram classification using delay differential equations doi: 10.1063/1.4811544
  * [R. P. AGARWAL and Y. M. CHOW(1985)] R. P. AGARWAL and Y. M. CHOW(1985). FINITE-DIFFERENCE METHODS FOR BOUNDARY-VALUE PROBLEMS OF DIFFERENTIAL EQUATIONS WITH DEVIATING ARGUMENTS


# Augmenting momentum resolution with well tuned probability distributions

Gregorio Landia Corresponding author. Giovanni E. Landib

a Dipartimento di Fisica e Astronomia  Universita’ di Firenze
Largo E. Fermi 2 50125 Firenze Italy
and INFN  Sezione di Firenze  Firenze  Italy
E-mail:

b E.T.T. S.a.g.l<EMAIL_ADDRESS>
Via Balestra 33
6900 Lugano  Switzerland.

###### Abstract

The realistic probability distributions of a previous article are applied to
the reconstruction of tracks in constant magnetic field. The complete forms
and their schematic approximations produce excellent momentum estimations,
drastically better than standard fits. A simplified derivation of one of our
probability distributions is illustrated. The momentum reconstructions are
compared with standard fits (least squares) with two different position
algorithms: the $\eta_{2}$-algorithm and the two-strip center of gravity. The
quality of our results are expressed as the increase of the magnetic field and
signal-to-noise ratio that overlap the standard fit reconstructions with ours
best distributions. The data and the simulations are tuned on the tracker of a
running experiment and its double sided microstrip detectors, here each
detector side is simulated to measure the magnetic bending. To overlap with
our best distributions, the magnetic field must be increased by a factor 1.5
for the least squares based on the $\eta_{2}$-algorithm and 1.8 for the two-
strip center of gravity for the low noise side, and 1.8 and 2.0 for the high
noise side. The signal-to-noise ratio must be increased by 1.6 for the low
noise side and 2.2 for the high noise side ($\eta_{2}$-algorithms). The fits,
built on the positioning with the center of gravity, are not modified by a
reduction of the signal-to-noise ratio.

###### keywords:

Particle tracking detectors, Performance of High Energy Physics Detectors, Si
microstrip and pad detectors, Analysis and statistical methods

##  1 Introduction

Some properties of our well tuned probability density functions (PDFs) were
described in ref. [1] and were tested on fits of simulated straight tracks.
The limitation to straight tracks was mainly due to the complexity of the
method we applied for the first time. Hence, working with two parameters, the
debugging and testing of the application can be followed on a surface.
Furthermore, the data, elaborated in our approach, were collected in a CERN
test beam [2] in the absence of a magnetic field. The detectors we used, were
few samples of double-sided silicon microstrip detector [3, 4] as those
composing the PAMELA tracker [5]. The results of ref. [1] showed a drastic
improvement of the fitted track parameters respect to the results of the least
squares methods. We observed excellent reconstructions even in presence of
very noisy hits, often called outliers, that generally produce worse fits.
This achievement is almost natural for our non gaussian PDFs. In fact, an
outlier hit is an event that is incompatible with the gaussian distribution
that is always assumed as the error distribution of the hit positions. Thus,
being the least squares a method strictly optimum for a gaussian PDF, it must
be modified in a somewhat arbitrary way to handle these pathological hits. On
the contrary, our PDFs are essentially different from a gaussian, and those
non standard cases are allowed (with low probability) by the non-gaussian
tails of the PDFs.

We have to recall that the perception of a rough handling of the hit
properties is well present in literature, and Kalman filter modifications are
often studied to accept extended deviations from a pure gaussian model. These
extensions imply the use of linear combinations of different gaussians [6],
but in a small number to avoid the intractability of the equations. The
unknown parameters are optimized from an improvements of the fits. In late
sense, our schematic approximations, used to initialize the maximum likelihood
search, could be reconnected to those extensions. In fact, to speed the
convergence, we calculate a different gaussian PDF (or more precisely an
effective variance) for each hit. Our comparisons take as references the least
squares methods. The Gauss-Markov theorem states their optimality among the
linear methods, thus no other linear method can be better.

The confidence gained in ref. [1] with straight tracks allow us to deal with
more complex tasks, i.e. to reconstruct the tracks in a magnetic field and
determine their momenta. To face this extension with the minimum modifications
of our previous work, we will utilize the simulated data and parameters used
in ref. [1], adapting them to the geometry and average magnetic field of the
PAMELA tracker. In fact, the relatively low magnetic field introduces small
modifications of the parameters obtained in its absence. The average signal
distribution of a minimum ionizing particle (MIP) is slightly distorted by a
Lorentz angle, that for the average magnetic field of the PAMELA tracker (
$0.43\,T$) is about $0.7^{\circ}$. Around this incidence angle, the average
signal distribution of a MIP is practically identical to that at orthogonal
incidence in the absence of a magnetic field, thus without further notice we
will assume a rotation of the detectors of their Lorentz angle. With these
assumptions, we simulate high momentum tracks for each side of the double
sided detectors. In the PAMELA tracker, the low noise side of the detector
(the junction side) is used for momentum measurements. This side has an
excellent resolution, a strip each two is connected to the read out system and
the unconnected strip distributes the incident signal on the nearby ones
optimizing the detector resolution. The other side (the ohmic side) has the
strips oriented perpendicularly to the junction side. Each strip is connected
to the readout system and has, by the complexity of construction, an higher
noise and a small, if any, signal spread to nearby strips. For those
characteristics, this side responds in a way similar to the types of
microstrip arrays used in the large trackers of the CERN LHC [7, 8, 9]. Thus,
the simulations on this side, as bending side, can give a glimpse of our
approach for other type of trackers even for the (small angle stereo) double-
sided detectors of the ALICE [7] experiment.

##  2 A compact derivation of the PDF

As always, the MIP incidence of our geometry imposes the use of the minimum
number of signal strips to reduce the noise, hence only two signal strips will
be used. In ref. [1] we indicated the principal steps required to obtain the
PDF for the center of gravity (COG) with two-strips, those steps followed the
standard method described in the books about probability. That procedure uses
the essential tools of the calculus: integrals and derivatives. In our cases,
many pieces of integrals over geometrical domains must be calculated to build
the cumulative probability, and its derivative gives our final PDF. It is
evident the length and the complexity of this development. Here we follow a
very different approach (Quantum Mechanics style) that reaches identical
results in few steps. The two strip COG (COG2 in the following) is calculated
with the origin in the center of the maximum-signal strip. The strip signals
are indicated with: $x_{1}$, $x_{2}$, and $x_{3}$, respectively the signal of
the right strip, central strip (with the maximum signal) and left strip. If
$x_{1}>x_{3}$ the COG2 is $x=x_{1}/(x_{1}+x_{2})$, if $x_{3}>x_{1}$ it is
$x=-x_{3}/(x_{3}+x_{2})$, other minor details will not be discussed here.
Thus:

$\displaystyle P_{x_{g2}}(x)=$
$\displaystyle\int_{-\infty}^{+\infty}\mathrm{d}\,x_{1}\int_{-\infty}^{+\infty}\mathrm{d}\,x_{2}\int_{-\infty}^{+\infty}\mathrm{d}\,x_{3}\,P(x_{1},x_{2},x_{3})$
(1)
$\displaystyle\,\Big{[}\theta(x_{1}-x_{3})\delta(x-\frac{x_{1}}{x_{1}+x_{2}})+\theta(x_{3}-x_{1})\delta(x+\frac{x_{3}}{x_{3}+x_{2}})\Big{]}$

where $P(x_{1},x_{2},x_{3})$ is the probability to have the signals
$x_{1},x_{2},x_{3}$ from the strips $1,2,3$. The signals $x_{i}$ are at their
final elaboration stage and ready to be used for position reconstruction of
the hit. The function $\theta(x_{j})$ is the Heaviside $\theta$-function:
$\theta(x_{j})=1$ for $x_{j}>0$, $\theta(x_{j})=0$ for $x_{j}\leq 0$, and
$\delta(x)$ is the Dirac $\delta$-function. It is immediate to verify the
normalization of $P_{x_{g2}}(x)$ by direct integration. Splitting the sum of
eq. 1 in two independent integrals and transforming the variables $x_{1}=\xi$,
$x_{1}+x_{2}=z_{1}$ and $x_{3}=\beta$, $x_{3}+x_{2}=z_{2}$, the jacobian of
the transformation is one and the integrals in $z_{1}$ and $z_{2}$ can be
performed with the rule:

$\int_{-\infty}^{+\infty}\mathrm{d}\,z\,F(z-\nu)\,\delta(x\mp\frac{\nu}{z})=F(\frac{\pm\nu}{x}-\nu)\,\frac{|\nu|}{x^{2}}\,.$
(2)

Applying eq. 2 to eq. 1 and using the limitations of the two
$\theta$-functions, eq. 1 becomes:

$\displaystyle P_{x_{g2}}(x)=\frac{1}{x^{2}}\Big{[}$
$\displaystyle\int_{-\infty}^{+\infty}\mathrm{d}\,\xi\int_{-\infty}^{\xi}\,\mathrm{d}\beta\,P\big{(}\xi,\,\xi\frac{1-x}{x},\beta\big{)}\,|\xi|\,+$
(3)
$\displaystyle\int_{-\infty}^{+\infty}\mathrm{d}\,\beta\,\int_{-\infty}^{\beta}\,\mathrm{d}\xi\,P\big{(}\xi,\,\beta\frac{-1-x}{x},\beta\big{)}\,|\beta|\Big{]}\,.$

This form underlines very well the similarity with the Cauchy PDF; in the
limit of $x\rightarrow\infty$ the $x$-part of the $P$ arguments are $-1$ and
$P_{x_{g2}}(x)\propto 1/x^{2}$ for large $x$.

###  2.1 The probability $P_{x_{g2}}(x)$ for small $x$

The probability $P(x_{1},x_{2},x_{3})$ can handle a strict correlation among
the arguments, we will release this strict correlation in a weakest one: the
mean values of the strip signals are correlated, but the fluctuations around
the mean values are independent. Thus the probability $P(x_{1},x_{2},x_{3})$
becomes the product of three functions $\\{P_{i}(x_{i}),i=1,2,3\\}$. Each
$P_{i}(x_{i})$ is assumed to be a gaussian PDF with mean values $a_{i}$ and
standard deviation $\sigma_{i}$. To simplify, the constants $a_{i},i=1,2,3$
are the noiseless signals released by a MIP with impact point $\varepsilon$:

$P_{i}(x_{i})=\exp\big{[}-\frac{(x_{i}-a_{i})^{2}}{2\sigma_{i}^{2}}\big{]}\frac{1}{\sqrt{2\pi}\sigma_{i}}.$
(4)

Even with the gaussian functions, the integrals of eq. 3 have no analytical
expressions and effective approximations must be constructed. We will not
report our final forms that are very long, instead we will illustrate a
limiting case which gives a simple approximation and eliminates a disturbing
singularity for the numerical integrations. It easy to show that for
$|x|\rightarrow 0$, $x^{-1}P_{2}\big{(}\xi(1-x)/x\big{)}$ and
$x^{-1}P_{2}\big{(}\beta(-1-x)/x\big{)}$ converge to two Dirac
$\delta$-functions. Hence, for small $|x|$, the integrals of eq. 3 can be
expressed as:

$\displaystyle P_{x_{g2}}(x)=$
$\displaystyle\frac{|a_{2}|}{\sqrt{2\pi}}\Big{\\{}\frac{\exp\big{[}-(x-\frac{a_{1}}{a_{1}+a_{2}})^{2}\frac{(a_{1}+a_{2})^{2}}{2\sigma_{1}^{2}(1-x)^{2}}\big{]}\big{[}1-\mathrm{erf}\big{(}(\frac{a_{3}}{a_{2}+a_{3}}-x)\frac{a_{2}+a_{3}}{\sqrt{2}\sigma_{3}(1-x)}\big{)}\big{]}}{2\sigma_{1}(1-x)^{2}}+$
(5)
$\displaystyle\frac{\exp\big{[}-(x+\frac{a_{3}}{a_{3}+a_{2}})^{2}\frac{(a_{3}+a_{2})^{2}}{2\sigma_{3}^{2}(1+x)^{2}}\big{]}\big{[}1-\mathrm{erf}\big{(}(\frac{a_{1}}{a_{2}+a_{1}}+x)\frac{a_{2}+a_{1}}{\sqrt{2}\sigma_{1}(1+x)}\big{)}\big{]}}{2\sigma_{3}(1+x)^{2}}\Big{\\}}\,.$

Equation 5 is correct only for $|x|\rightarrow 0$, but it is useful beyond
this limit and contains many ingredients of more complex expressions. An
example can be seen in fig.2 of ref. [1], the essential elements are the two
bumps (gaussian-like) centered in the possible noiseless two-strip COG:
$a_{1}/(a_{1}+a_{2})$ and $-a_{3}/(a_{3}+a_{2})$. Their effective standard
deviations are modulated by the signal-to-noise ratio
$\sigma_{1}/(a_{2}+a_{1})$ and $\sigma_{3}/(a_{2}+a_{3})$ and by the $1\pm x$
factor. More complete expressions contain terms very similar to Cauchy PDFs
that are $\propto 1/x^{2}$ for large $x$. The dimensions of the constants
$a_{j}$ must be those of the $\sigma_{j}$, for both of them we take directly
the ADC counts. The $x$-variable (the COG2) is a pure number expressed as a
fraction of the strip size, or more precisely, the strip size is the scale of
lengths.

The form of the COG2 inserted in eq. 1 synthesizes well our positioning
algorithm: each strip signal around the strip with the maximum signal is used
for position reconstruction. The strategy of cluster detection is supposed to
be optimized to suppress the false hits and its results are used only for the
maximum-signal-strip selection, and the signals of the two lateral strips are
used in any case even for small positive or not too negative values. It must
be underlined that the tiny amount of information is relevant for positioning
or to maintain the noise distribution. In fact, the track reconstruction
algorithms reach their optimum with realistic probability distributions, and
any unaccounted distortion induces a loss of resolution.

For construction, $P_{x_{g2}}(x)$ gives to the probability of an $x$ value
with constant $\\{a_{j}\\}$. For the track fitting this simple probability is
useless, we need a PDF that, for any $x$-value, gives the probabilities of the
impact points $\\{\varepsilon\\}$. To insert this functional dependence on
$P_{x_{g2}}(x)$, we must know the average variations of the energies
$\\{a_{1},a_{2},a_{3}\\}$ with $\varepsilon$. For this task, it is essential a
theorem illustrated in ref. [1] and the details of refs. [10, 11, 12], this
theorem allows the extraction of these functional dependencies from the data.
Once extended in $\varepsilon$, the PDF can be rewritten as:

$P_{x_{g2}}(x,E_{t},\varepsilon)=\frac{F(a_{1}(\varepsilon),a_{2}(\varepsilon),a_{3}(\varepsilon),E_{t},\sigma_{1},\sigma_{2},\sigma_{3},x)}{x^{2}}\,.$
(6)

Where
$F(a_{1}(\varepsilon),a_{2}(\varepsilon),a_{3}(\varepsilon),E_{t},\sigma_{1},\sigma_{2},\sigma_{3},x)$
is the term in square brackets of eq. 3. The functions
$a_{1}(\varepsilon),a_{2}(\varepsilon),a_{3}(\varepsilon)$ are the average
noiseless fractions of signal collected by each strip and $E_{t}$ is the total
signal of the three strips: the central one with the maximum signal and the
two lateral. We would need the noiseless $E_{t}$, but the measured one is the
only possibility. The extraction of the functions $\\{a_{i}(\varepsilon)\\}$
from real data of ref. [1] necessitates further refinements and numerical
recipes. But, once defined and checked their results, the procedure can be
standardized. In any case, slight variations of the $\\{a_{i}(\varepsilon)\\}$
around the best one give almost identical track parameter distributions, thus
their selection is very important but non critical. The functions
$\\{a_{i}(\varepsilon)\\}$ have a definition similar to the "templates" of
ref. [13]. The parameters $\sigma_{1},\sigma_{2},\sigma_{3}$ are the noise
(eq. 4) of the three strips considered. Our PDF can easily accomodate strips
with different $\sigma$, but, in the simulations, we will use a single
$\sigma$ for all the strips of same detector side and these parameters will
not be reported in the future expressions. Equation 6 is normalized as
function of $x$ but is not normalized in $\varepsilon$, due to its irrelevance
for the fit, this normalization will be neglected here. A set of normalized
PDFs from eq. 6 are illustrated in ref. [1] as $\varepsilon$ functions.

###  2.2 Track definition

Given our exploration of a new approach, our track definition must be the
easiest allowed. The tracks are circles with a large radius to simulate the
high momentum MIPs where the multiple scattering is negligible. The relation
of the track parameters to the momenta is $p=0.3\,\mathrm{B}\,R$, the form,
adapted to our needs, from ref. [14] where $p$ is in $\mathrm{GeV/c}$,
$\mathrm{B}$ in Tesla and $R$, the track radius, in meters.
The tracker model is formed by six parallel equidistant ($89\,mm$) detector
layers as for the PAMELA tracker, with constant magnetic field of $0.43\,T$
perpendicular to the track plane ($\xi\,,z$). The center of the tracker has
coordinates $\xi=0\,,z=0$, the $\xi$ axis is parallel to the layer planes and
the $z$ axis is perpendicular. The tracks are circles with center in $\xi=-R$
and $z=0$, and the magnetic field is parallel to the analyzing strips. To
simplify the geometry, the overall small rotation of the Lorentz angle
($0.7^{\circ}$) is neglected now, but it will be introduced in the following.
The simulated hits are generated (ref. [1]) with a uniform random distribution
on a strip, they are collected in groups of six to produce the tracks. For the
least squares, the exact impact position $\varepsilon$ of each hit is
subtracted from its reconstructed $\eta_{2}(x_{g2})$ position (as defined in
refs. [1, 11, 15]) and it is added the value of the fiducial track for the
corresponding detector layer. In this way each group of six hits defines a
track with our geometry and the error distribution of $\eta_{2}$-algorithm.
Identically for the COG2 positioning algorithm. This hit collection simulates
a set of tracks populating a large portion of the tracker system with slightly
non parallel strips on different layers (as it is always in real detectors).
All these tracks are possible realizations of our model track. At our high
momenta the track bending is very small and the track sagitta is smaller than
the strip width. Thus, the bunch of tracks has a transversal section around a
strip size, on this width we have to consider the Lorentz angle but its effect
is clearly negligible. Our preferred position reconstruction is the $\eta_{2}$
algorithm as in ref. [1] because it gives parameter distributions better than
those obtained with the simplest COG2 positions, but even the results for the
COG2 will be reported in the following.
In the $\\{{\xi},{z}\\}$-plane, the circular tracks are approximated, as
usual, with parabolas linear in the track parameters:

$\displaystyle\xi=\beta+\gamma{z}-\alpha{z}^{2}=\varphi(z)\,$ (7)
$\displaystyle\xi_{n}=\beta_{n}+\gamma_{n}{z}-\alpha_{n}{z}^{2}=\varphi_{n}(z)\,.$

The first line of eq. 7 is the model track, the second line is the fitted one.
The circular track is the osculating circle of the parabola $\varphi(z)$, at
our high momenta and tracker size the differences are negligible. Our model
track has $\gamma=0,\,\beta=0$, and $1/\alpha$ is proportional to the track
momentum. Due to the noise, the reconstructed track has equation
$\varphi_{n}(z)$ and the fitted parameters
$\\{\alpha_{n},\beta_{n},\gamma_{n}\\}$ are distributed around the model
values $\alpha,\beta,\gamma$. For their non gaussian forms, our PDFs must be
used with the non-linear search of the likelihood maxima, but, as always, the
search will be transformed in a minimization. The momentum and the other
parameters of the track $n$ are obtained minimizing, respect to $\alpha_{n}$,
$\beta_{n}$ and $\gamma_{n}$ the function $L(\alpha_{n},\beta_{n},\gamma_{n})$
defined as the negative logarithm of the likelihood with the PDFs of eq. 6:

$\displaystyle
L(\alpha_{n},\beta_{n},\gamma_{n})=-\sum_{j=6n+1}^{6n+6}\,\ln[P_{x_{g2}}(x(j),E_{t}(j),\psi_{j}(\alpha_{n},\beta_{n},\gamma_{n})]$
(8)
$\displaystyle\psi_{j}(\alpha_{n},\beta_{n},\gamma_{n})=\varepsilon(j)-\varphi(z_{j})+\varphi_{n}(z_{j})\,.$

The parameters $x(j)$, $E_{t}(j)$ (introduced in eq. 6) are respectively: the
COG${}_{2}(j)$ position, the sum of signal in the three strips for the hit $j$
and the track $n$, ${z}_{j}$ is the position of the detector plane $j$ of the
$n$th track. The $\varepsilon$-dependence in the $\\{a_{i}(\varepsilon)\\}$ is
modified to $\psi_{j}(\alpha_{n},\beta_{n},\gamma_{n})$, to place the impact
points on the track. In real data,
$\varepsilon(j)-\varphi(\alpha,\beta,\gamma,z_{j})$ is absent (and unknown)
but the data are supposed to be on a track. We can easily use a non linear
form for the function $\psi_{j}(\alpha_{n},\beta_{n},\gamma_{n})$, but in this
case is of scarce meaning. In more complex cases, non linearities of various
origin can be easily implemented.

We will reserve the definition of maximum likelihood evaluation (MLE) to the
results of eq. 8 even if the least squares method can be derived by a maximum
likelihood. The minimum search routine for the MLE is initialized as in ref.
[1]. The initial track parameters are given by a weighted least squares with
weights given by an effective variance ($\sigma_{eff}(i)^{2}$) for each hit
$(i)$ and $\eta_{2}(i)$ as hit position. The $\sigma_{eff}(i)^{2}$ is obtained
from our PDF, but, for the form of eq. 6, the variance is an hill defined
parameter even in $\varepsilon$, and this hill definition must be eliminated
with cuts on the integration ranges that suppress the PDF tails. We use two
sizes of cuts, one for the low noise side and one for the high noise side of
our double side detector. The cuts are optimized to obtain gaussian
distributions, with standard deviation $\sigma_{eff}(i)$, reproducing our PDF
for a good hit. For a set of hits (good hits), eq. 6 has the form of a narrow
high peak and a gaussian (centered in $\eta_{2}(i)$) with standard deviation
$\sigma_{eff}(i)$ is built to reproduce well the
$P_{x_{g2}}(x(i),E_{t}(i),\varepsilon)$ on few of them. These gaussian
approximations look good in linear plots, the logarithmic plots show marked
differences even in these happy cases, the tails are non-gaussian. The cuts,
so defined, are used for the $\sigma_{eff}(i)$ extraction in all the other
hits, even where $P_{x_{g2}}(x(i),E_{t}(i),\varepsilon)$ is poorly reproduced
by a gaussian.

The $\\{\alpha_{n},\beta_{n},\gamma_{n}\\}$ given by the weighted least
squares are almost always near those given by the minimization of eq. 8, thus
rendering less heavy the minimum search to the MATLAB [17] $fminsearch$
routine. The closeness of these approximations to the MLE supports the non
criticality of the extraction of the functions $\\{a_{i}(\varepsilon)\\}$. In
fact, the approximate gaussian distributions are often very different from the
$P_{x_{g2}}(x(i),E_{t}(i),\varepsilon)$, but the realistic pieces of
information about the hits are sufficient to produce near optimal results.
When the tails of the PDFs are important the MLE results are better. The
strong variations of the set $\\{\sigma_{eff}\\}$ on the strip are illustrated
in figs. 6 and 11 of ref. [1], and the patterns of those variations support a
complex interplay among the physical properties of the strips and the
reconstruction algorithm. The time-consuming extraction of $\sigma_{eff}(i)$
(as indicated in ref. [1]) can be reduced building the surface
$\sigma_{eff}(x,E_{t})$ and calculating $\sigma_{eff}(i)$ with an
interpolation.

##  3 Low noise, high resolution, floating strip side

The floating strip side is the best of the two sides of this type of strip
detector. It is just this side that measures the track bending in the PAMELA
magnetic spectrometer. In the test beam [2], the noise of the "average" strips
is well reproduced by a gaussian with a standard deviation of 4 ADC counts,
and the PDF for the sum of three strip signals has its maximum at 142 ADC
counts with the most probable signal-to-noise ratio of $35.5$ for a three
strip cluster ($SNR(n)=\sum_{i=1}^{3}x_{i}(n)/\sigma_{i}$). The functions
$a_{j}(\varepsilon)$ are those of ref. [1] for this strip type. In the
simulations we will use a high momentum of $350\,GeV/c$. For this momentum and
similar geometry, we have a report with some histograms of a CERN test beam of
the PAMELA tracker before its installation in the satellite. The tracks were
reconstructed with the $\eta_{2}$ algorithm, and this allows a positive check
of our simulations.

At $350\,GeV/c$, the hits of the tracks are largely contained

Figure 1: Left plot. Blue line: distribution of the differences of hit
$\\#\,3$ minus hit $\\#\,1$ of a track for the $\eta_{2}$-algorithm; magenta
line: the same for the COG2. Right plot. True residuals of the reconstructed
tracks: from MLE (red), from $\sigma_{eff}(i)$ (black), from $\eta_{2}$
(blue), from COG2 (magenta), and the position errors of the $\eta_{2}$ (cyan)
and COG2 (green).

in a strip width ($51\,\mu m$), the left side of fig. 1 shows the
distributions of the differences between two hits on a track (the $\\#\,3$
minus the $\\#\,1$) for two different hit reconstruction algorithm. They peak
around 8 $\mu m$, but, as any COG reconstruction algorithm, the COG2 algorithm
produces a wider distribution due to a systematic error of refs. [10]. The
$\eta_{2}$-algorithm is built to be free of this systematic error, here we
corrected even the asymmetry errors discussed in ref. [11, 12]. Even if they
are now very small, the clear symmetry of the plots is the product of this
accuracy. For the absence of the COG systematic error, the PDFs of the track
parameters given by the $\eta_{2}$ algorithm are better than those given by
the COG2. The results of our MLE are compared with those obtained by three
different least squares. Having to plot the results of four type of
reconstructions we will use the following color convention:

  * •

red lines refer to our MLE (eq. 8),

  * •

black lines are the weighted least squares with weight $1/\sigma_{eff}(i)^{2}$
and $\eta_{2}$ position,

  * •

blue lines for the least squares with the $\eta_{2}$ position algorithm

  * •

magenta lines for the least squares with the COG2 position algorithm.

In the right plot of fig. 1 the distributions (or more precisely the histogram
values divided by the number of entries and the step size) of the differences
of the fitted positions respect to the exact ones are reported with the above
color convention. In the following we will call these differences as true
residuals (allowed only in simulations) to distinguish from the residuals
generally defined as the differences of the fitted positions from the
reconstructed positions. For further tests, two other lines are added, the
cyan and green lines are the error distributions of the $\eta_{2}$ and COG2
hit reconstruction algorithms. The PDFs of the residuals are not reported,
but, those for the COG2 and $\eta_{2}$ are almost identical to the reported
PDFs of the true residual. The residuals for other two approaches are very
different from the true residuals of fig. 1, they have very high peaks around
zero. Our PDF of eq. 6 and $\sigma_{eff}(i)$ allow the recognition of the good
hits and the fit optimization selects to pass near to them, giving an high
frequency of small residuals.

The general expectancy of an improvement of the position reconstruction by
redundant data is evidently disappointed for the $\eta_{2}$ least squares: the
cyan distribution is higher than the blue one. This mismatch is very similar
to the case of a least squares on data from a Cauchy distribution [16].
Instead, the true residuals of COG2 least squares look much better than their
error distribution, that is very large. The least squares looks able to round
the error distribution. But this effort consumes all its power, in fact if the
statistical noise is suppressed, the COG systematic error does not allow any
modification of the true residuals. On the contrary, the true residual PDF of
the $\eta_{2}$ least squares grows toward a Dirac $\delta$ function in the
absence of the statistical noise. We have to signal a light inconsistency
among the histograms reported here and those of ref. [1]. There we worked
always with the strip pitch as scale for the length (as defined for the
probability of eq. 5), in the plots the horizontal scale was turned to the
appropriate scales ($\mu m$). Here, the conversion to the correct dimensions
was performed before building the histograms, thus the vertical scales turn
out different respect to those of ref. [1].

Figure 2: Left plot. Distributions of the curvature in $(GeV/c)^{-1}$ (the
$\alpha$ parameters) for the four types of fit. Right plot: distributions of
the reconstructed momenta.

As illustrated in fig. 2, the MLEs give the best results for the momentum
reconstruction, and the weighted least squares, are very near to them. The
fits of the standard least squares with $\eta_{2}$ or COG2 positioning
algorithms show a drastic decrease in resolution. The use of the simple COG2
algorithm is the worst one. Often the distributions of the left side of fig. 2
are reported as resolution of the momentum reconstruction, the k-value of ref.
[14]. For the $\eta_{2}$ and COG2 least squares, the plots of the momentum
distributions have appreciable shifts of the maxima (most probable value)
respect to the fiducial value of 350 $GeV/c$, the shifts are negligible for
the other two fits. These shifts are mathematical consequences of the change
of variable from curvature to momentum.

###  3.1 Other track parameters

The complete track reconstruction must consider even the other two parameters
of a track, the $\beta_{n}$ and $\gamma_{n}$. Their fits give very similar
distributions to those plotted in ref. [1]. The maxima of the distributions
are now a little lower, in particular for the $\beta$ parameter. This is not
unexpected, in ref. [1] we had 3 degree of freedom for two parameters, here we
have 3 degree of freedom for three parameters, an effective reduction of the
redundancy that has a slight effect on the results.

##  4 High noise, low resolution, normal strip side

The other side of the double sided silicon microstrip detector has very
different properties respect to the floating strip side (the junction side).
We will not recall the special treatments required to transform a ohmic side
in strip detector and all the other particular setups necessary to its
functioning. From the point of view of their data, this side produces data
very similar to a normal strip with a gaussian noise of 8 ADC counts, twice of
the other side (most probable signal-to-noise ratio $SNR(n)=18.2$). The
absence of the floating strips gives to the histograms of the COG2 the normal
aspect with a high central density and a drop around $x_{g2}=0$. No additional
rises around $x_{g2}=\pm 1/2$ are present, they are typical of the charge
spreads given by the floating strips. This absence reduces the efficiency of
the positioning algorithms that gain substantially from the charge spread.
Now, the functions $\\{a_{i}(\varepsilon)\\}$ of ref. [1] are similar to those
of an interval function with a weak rounding to the borders, this rounding is
mainly due to the convolution of the strip response with the charge spread
produced by the drift in the collecting field. If a residual capacitive
coupling is present, it is very small. In any case, it is just this rounding
that renders very good the resolution of the hits on the strip borders with a
small $\sigma_{eff}(i)$. Due to the differences respect to the other side, we
have to implement the simulations with a lower momentum of $150\,GeV/c$.

Figure 3: Left plot. Blue line: distribution of the differences of hit
$\\#\,3$ minus hit $\\#\,1$ of a track for the $\eta_{2}$-algorithm, magenta
line: the same for the COG2. Right plot. True residuals of the reconstructed
tracks. Color cades defined above for the fits, the cyan and green PDFs are
the positioning errors of the $\eta_{2}$ and COG2 algorithms.

The left side of fig. 3 shows the distributions of the hit differences
($\\#\,3$ minus $\\#\,1$) to test the containment of the track on a strip
width, the largest noise has evident effects on their forms. As usual the COG2
distribution is larger than the $\eta_{2}$ distribution. The true residuals of
the fits are reported in the right side of fig. 3 with the cyan and green
lines for the error distributions of the $\eta_{2}$ and COG2 hit
reconstruction algorithms. The other color codes are the standard ones. Even
here the data redundancy for the $\eta_{2}$-fit does not improve the position
reconstruction. The $\eta_{2}$ hit error distribution is drastically better
than that of the true residuals for the fit. The COG2 position error
distribution has two maxima due to the systematic error that is positive for
the first half of the strip and negative in the second half, this difference
of sign combines with the noise random error to round the two maxima. The fit,
based on the COG2 positions, has a good probability to partially average these
sign differences and produce a single wide maximum. The momentum distributions
given by the four fits are reported in fig. 4.

Figure 4: Left plot. Distributions of the curvature in $(GeV/c)^{-1}$ (the
$\alpha$ parameters) for the four different fits. Right plot: distributions of
the reconstructed momenta.

The results of our MLE and the weighted least squares are drastically better
than the standard least squares. Now the shift of the maxima for the momentum
distributions of two least squares are more evident than in fig. 2.

##  5 Discussion

A meaningful comparison among the different fits is somewhat difficult. Often
the standard deviation of the PDF is extracted interpolating a gaussian
function on non gaussian distributions, or two (or more) gaussian PDF are
fitted hoping to describe two (or more) independent gaussian processes. In
ref. [1] we used the full width at half maximum, but, apart from its complex
extraction, it does not characterize very well the distributions. In any case
those differences are not of easy interpretation. Here we will use two very
precious resources contributing to the momentum measurements: the magnetic
field and the signal-to-noise ratio. The simulated magnetic field intensity
and the signal-to-noise ratio are increased in the $\eta_{2}$ and COG2 least
squares reconstructions to reach our best distributions (the red lines).

###  5.1 Increasing the magnetic field

Figure 5: Top plots. Low noise case, Curvature in $(GeV/c)^{-1}$ and momentum
distributions for our MLE (red lines) and $\eta_{2}$ least squares with a
magnetic field increased by a factor 1.5. Bottom plots. Higher noise case,
here the $\eta_{2}$ least square has a magnetic field 1.8 times greater than
the MLE.

With fixed momentum, the magnetic field is increased in the fit for the
$\eta_{2}$ least squares. The upper sector of fig. 5 illustrates these results
for the low noise side and the overlaps with our MLE. The red lines are
identical to those of fig. 2, the blue line are the $\eta_{2}$ least squares
for a magnetic field 1.5 times higher. The plots for the COG2 least squares
are not reported to render easily legible the figures, in this case the
magnetic field must be increased of factor 1.8 to overlap our red lines. The
lowest sector of fig. 5 reports the higher noise side, the red lines are
identical to those of fig. 4, the blue lines ($\eta_{2}$ least squares)
overlap the red lines with a magnetic field increased of 1.8 times. Even here,
the overlap with the COG2 least squares is not reported, now, to obtain a
reasonable overlap, the magnetic field must be doubled. Clearly, the increment
of the magnetic field produces other slight differences in the tails of the
distributions, but the main results are well evidenced.

###  5.2 Increase of the signal-to-noise ratio

Let us see now the effects of the increase of the signal-to-noise ratio. These
modifications reproduce in part the results of the magnetic increase, but it
is different in other parts. In fact, higher values of the magnetic field do
not modify the form of the distributions of the curvature $\alpha$ (with
dimension ${legth}^{-1}$), the curvature distributions translate toward higher
values for a reduction of the track radius. The dimensional transformations to
$(\mathrm{GeV/c})^{-1}$ and $\mathrm{GeV/c}$ have the magnetic field intensity
as scaling factor, that can rise the distributions to the forms of fig. 5 even
for the COG2 case. But, as for the $\alpha$ parameters (with dimension
${length}^{-1}$), the distributions of the true residuals of fig. 1 and fig. 3
are not influenced by the rising of the magnetic, a part small effects on the
Lorentz angle.

Keeping the magnetic field to $0.43\,T$, the raise of the signal-to-noise
ratio modifies the distributions of the $\alpha$ parameters and those of the
true residuals and can bring them to overlap our red distributions. In our
simulations, this is accomplished scaling the amplitude of the random noise,
added to the strip signals, and rerunning the least squares fits. The plots of
these results for the momentum (in $\mathrm{GeV/c}$) and curvature (in
$(\mathrm{GeV/c})^{-1}$) are not reported because they are practically
identical to those of fig. 5. For the low noise side the gaussian strip noise
of $\sigma=4$ ADC counts must be reduced to 2.5 ADC counts, raising the
signal-to-noise ratio by a factor 1.6 to a huge most probable value of
$SNR(n)=56.8$. Now, even the blue line in fig. 1 of the true residuals reaches
the red line, but the redundancy continues to be unable to move the true
residual distribution beyond that of the new $\eta_{2}$ errors. A special
mention must be devoted to the COG2 distributions, the magenta lines. They
remain essentially unchanged in any plot for a (any) reduction of the strip
noise. The COG2 error distribution of fig. 1 (the green line) shows a slight
modification becoming less rounded for the hard presence of the systematic
error, unaffected by the random noise.

For the higher noise side of the detector, a reduction to less of one half of
the noise (from $8$ to $3.6\ \ \mathrm{ADC}$-counts) is required to reproduce
the low parts of fig. 5. With this doubling (2.2) of the signal-to-noise ratio
the blue line (the $\eta_{2}$ best fit) overlaps the red line in fig. 3 rising
the most probable signal-to-noise ratio to 40.5. As for the other side, no
noise reduction moves the distributions given by the COG2 least-squares. The
magenta lines remain practically unchanged a part a reduction of the
fluctuations and a sharpening of the green line of fig. 3.

As just said above, the insensitiveness to the noise reduction is due to the
COG systematic error of ref. [10]. This fact renders of weak relevance the
efforts to improve the signal-to-noise ratio in the detectors if the
positioning algorithm remains the COG. On the other side, this method is
stable respect to a decrease of the signal-to-noise ratio due to ageing or
radiation damages, as far as the COG systematic error dominates.

A last point remains to be analyzed, i.e. the huge difference between the
curvature PDF for the $\eta_{2}$ least squares in the low part fig. 5 and the
blue line of fig. 2. Now the two detector types have very similar signal-to-
noise ratio (3.6 ADC counts in the first and 4 ADC counts in the second), but,
to reach the overlap of the two PDFs, another factor of around 2.4 is needed.
This factor is too large to be due to the $20\,\%$ of differences between the
sizes of their strips. The main difference must be due to the beneficial
charge sharing of the floating strip and the nearing of its detector
architecture to the ideal detector defined in ref. [10].

##  6 Conclusions

We extended the simulated track reconstructions of ref. [1] with the study of
curved tracks in a constant magnetic field of $0.43\,T$. All the other
effects, $\delta$-rays, multiple scattering (negligible at our high momenta),
energy loss, etc., are explicitly excluded by the simulation, focusing on the
differences among various fit methods. Each type of detectors of our double
sided detector are examined as momentum analyzer. The higher noise side has
signal properties very similar to single sided detectors of large use in
running trackers. We use our well tuned PDFs in two form: complete in the MLE
or schematic with weight parameters $1/\sigma_{eff}(i)^{2}$. Each one of these
two form is very effective for this task respect to the results of two other
type of track fitting explored i.e. the least squares with $\eta_{2}$ as
position algorithm and the least squares with the COG2 positioning. To
establish a comparison among these fitting methods, we skip the usual ways of
comparing standard deviations or full width at half maximum. Instead, we
modify two different tracker properties to reach an overlap among the fit
outputs: the magnetic field and the signal-to-noise ratio. For the overlap of
the two standard fits with our best distributions, the magnetic field must be
increased by a factor 1.5 for the $\eta_{2}$ fit and 1.8 for the COG2 in the
low noise side and 1.8 and 2 for the higher noise side. The increase of
signal-to-noise ratio is effective only for the $\eta_{2}$ based least
squares, the overlaps are obtained with factors 1.6 and 2.2 for the two
detector sides. Any increase of the signal-to-noise ratio is unable to improve
the curvature and momentum distributions in the case of COG2 positioning, the
COG intrinsic systematic error survives untouched to any reduction of the
detector random noise. In any case, the drastic improvement of our well tuned
PDFs in the momentum reconstructions is evident. We have to remind that these
are simulations and come with all their uncertainties. Assuming an optimistic
view, this increase in resolution can be spent in different ways either for
better results on running experiments or in reducing the complexity of future
experiments if the baseline fits (almost always based on the COG positioning)
are estimated sufficient. With eq. 5, we start to give a glimpse of the
analytical forms of our PDFs. Even if this expression is of limited validity,
its combination with the appropriate $\\{a_{j}(\varepsilon)\\}$ allows a
faithful reproduction of the green and cyan distributions of figs. 1 and 3
from the data. Those distributions are exclusive products of simulations, but,
with slight redefinitions, the simplified PDF is able to reproduce them. The
complex task of MLE requires more advanced expressions, future papers will be
devoted to their discussion.

## References

  * [1] G. Landi, G.E. Landi, Improvement of track reconstruction with well tuned probability distributions 9 2014 P10006.
  * [2] O. Adriani et al., "In-flight performance of the PAMELA magnetic spectrometer" $16^{th}$ _International Workshop on Vertex Detectors_ September 2007 NY. USA. PoS(Vertex 2007)048
  * [3] G. Batignani et al., Double sided read out silicon strip detectors for the ALEPH minivertex http://dx.doi.org/10.1016/0168-9002(89)90546-9 _Nucl. Instrum. and Meth. A_ 277 (1989) 147.
  * [4] O. Adriani et al., [L3 SMD Collaboration] The New double sided silicon microvertex detector for the L3 experiment. http://dx.doi.org/10.1016/0168-9002(94)90774-9 _Nucl. Instrum. Meth. A_ 348 (1994) 431.
  * [5] P. Picozza et al., PAMELA-A Payload for Matter Antmatter Exploration and Light-nuclei Astrophysics http://dx.doi.org/10.1016/j.astropartphys.2006.12.002 _Astropart. Phys._ 27 (2007) 296 [0608697], O. Adriani et. al., An anomalous positron abundance in cosmic rays with energies 1.5-100 GeV http://dx.doi.org/10.1038/nature07942 _Nature_ 458 (2009) 607
  * [6] R. Fr$\mathrm{\ddot{u}}$hwirth Track fitting with non-Gaussian noise http://dx.doi.org/10.1016/S0010-4655(96)00155-5 _Computer Physics Communications_ 100 (1997) 1. R. Fr$\mathrm{\ddot{u}}$hwirth, T. Speer, A Gaussian-sum filter for vertex reconstruction http://dx.doi.org/10.1016/j.nima.2004.07.090 _Nucl. Instr. and Meth. A_ 534 (2004) 217
  * [7] ALICE Collaboration, The ALICE experiment at the CERN LHC. 3 2008 S08002
  * [8] ATLAS Collaboration, The ATLAS experiment at the CERN Large Hadron Collider. 3 2008 S08003
  * [9] CMS Collaboration, The CMS experiment at the CERN Large Hadron Collider. 3 2008 S08004.
  * [10] G. Landi, Properties of the center of gravity as an algorithm for position measurements http://dx.doi.org/10.1016/S0168-9002(01)02071-X _Nucl. Instrum. and Meth. A_ 485 (2002) 698. G. Landi, _Properties of the center of gravity as an algorithm for position measurements: two-dimensional geometry_ http://dx.doi.org/10.1016/S0168-9002(02)01822-3 _Nucl. Instrum. and Meth. A_ 497 (2003) 511.
  * [11] G. Landi, _Problems of position reconstruction in silicon microstrip detectors_ , http://dx.doi.org/10.1016/j.nima.2005.08.094 Nucl. Instr. and Meth. A 554 (2005) 226.
  * [12] G. Landi and G.E. Landi, "Asymmetries in Silicon Microstrip Response Function and Lorentz Angle"  [http://arxiv.org/abs/1403.4273physics/arXiv:1403.4273]
  * [13] The CMS Collaboration _Description and performance of track and primary-vertex reconstruction with the CMS tracker_ 9 2014 P10009
  * [14] K.A. Olive et al. _Particle Data Book_ Chin. Phys. C 38 090001 (2014)
  * [15] E. Belau et al., _Charge collection in silicon strip detector_ http://dx.doi.org/10.1016/0167-5087(83)90591-4 _Nucl. Instrum. and Meth. A_ 214 (1983) 253.
  * [16] G. Landi Tracks of detector theory Slides for a seminar to the Firenze Theory Group, October 2014
  * [17] MATLAB 8. The MathWorks Inc.