// @(#)root/hist:$Name:  $:$Id: TMultiDimFit.cxx,v 1.10 2003/01/16 18:07:52 brun Exp $
// Author: Christian Holm Christensen 07/11/2000

//____________________________________________________________________
//
//
// 
/*



Multidimensional Fits in ROOT







Overview





A common problem encountered in different fields of applied science is
to find an expression for one physical quantity in terms of several
others, which are directly measurable.



An example in high energy physics is the evaluation of the momentum of
a charged particle from the observation of its trajectory in a magnetic
field.  The problem is to relate the momentum of the particle to the
observations, which may consists of of positional measurements at
intervals along the particle trajectory.



The exact functional relationship between the measured quantities
(e.g., the space-points) and the dependent quantity (e.g., the
momentum) is in general not known, but one possible way of solving the
problem, is to find an expression which reliably approximates the
dependence of the momentum on the observations.



This explicit function of the observations can be obtained by a
least squares fitting procedure applied to a representive
sample of the data, for which the dependent quantity (e.g., momentum)
and the independent observations are known. The function can then be
used to compute the quantity of interest for new observations of the
independent variables.



This class TMultiDimFit implements such a procedure in
ROOT. It is largely based on the CERNLIB MUDIFI package
[2]. Though the basic concepts are still sound, and
therefore kept, a few implementation details have changed, and this
class can take advantage of MINUIT [4] to improve the errors
of the fitting, thanks to the class TMinuit.



In [5] and [6] H. Wind demonstrates the utility
of this procedure in the context of tracking, magnetic field
parameterisation, and so on. The outline of the method used in this
class is based on Winds discussion, and I refer these two excellents
text for more information.



And example of usage is given in
$ROOTSYS/tutorials/multidimfit.C.








The Method





Let $ D$ by the dependent quantity of interest, which depends smoothly
on the observable quantities 
$ x_1, \ldots, x_N$, which we'll denote by

$ \mathbf{x}$. Given a training sample of $ M$ tuples of the form,
(TMultiDimFit::AddRow)



$\displaystyle \left(\mathbf{x}_j, D_j, E_j\right)\quad,
$


where 
$ \mathbf{x}_j = (x_{1,j},\ldots,x_{N,j})$ are $ N$ independent
variables, $ D_j$ is the known, quantity dependent at 
$ \mathbf{x}_j$,
and $ E_j$ is the square error in $ D_j$, the class
TMultiDimFit
will
try to find the parameterization








$\displaystyle D_p(\mathbf{x}) = \sum_{l=1}^{L} c_l \prod_{i=1}^{N} p_{li}\left(x_i\right) = \sum_{l=1}^{L} c_l F_l(\mathbf{x})$

(1)





such that








$\displaystyle S \equiv \sum_{j=1}^{M} \left(D_j - D_p\left(\mathbf{x}_j\right)\right)^2$

(2)





is minimal. Here 
$ p_{li}(x_i)$ are monomials, or Chebyshev or Legendre
polynomials, labelled 
$ l = 1, \ldots, L$, in each variable
$ x_i$, 
$ i=1, \ldots, N$.



So what TMultiDimFit does, is to determine the number of
terms $ L$, and then $ L$ terms (or functions) $ F_l$, and the $ L$
coefficients $ c_l$, so that $ S$ is minimal
(TMultiDimFit::FindParameterization).



Of course it's more than a little unlikely that $ S$ will ever become
exact zero as a result of the procedure outlined below. Therefore, the
user is asked to provide a minimum relative error $ \epsilon$
(TMultiDimFit::SetMinRelativeError), and $ S$
will be considered minimized when



$\displaystyle R = \frac{S}{\sum_{j=1}^M D_j^2} < \epsilon
$





Optionally, the user may impose a functional expression by specifying
the powers of each variable in $ L$ specified functions 
$ F_1, \ldots,
F_L$ (TMultiDimFit::SetPowers). In that case, only the
coefficients $ c_l$ is calculated by the class.








Limiting the Number of Terms





As always when dealing with fits, there's a real chance of
over fitting. As is well-known, it's always possible to fit an
$ N-1$ polynomial in $ x$ to $ N$ points $ (x,y)$ with 
$ \chi^2 = 0$, but
the polynomial is not likely to fit new data at all
[1]. Therefore, the user is asked to provide an upper
limit, $ L_{max}$ to the number of terms in $ D_p$
(TMultiDimFit::SetMaxTerms).



However, since there's an infinite number of $ F_l$ to choose from, the
user is asked to give the maximum power. $ P_{max,i}$, of each variable
$ x_i$ to be considered in the minimization of $ S$
(TMultiDimFit::SetMaxPowers).



One way of obtaining values for the maximum power in variable $ i$, is
to perform a regular fit to the dependent quantity $ D$, using a
polynomial only in $ x_i$. The maximum power is $ P_{max,i}$ is then the
power that does not significantly improve the one-dimensional
least-square fit over $ x_i$ to $ D$ [5].



There are still a huge amount of possible choices for $ F_l$; in fact
there are 
$ \prod_{i=1}^{N} (P_{max,i} + 1)$ possible
choices. Obviously we need to limit this. To this end, the user is
asked to set a power control limit, $ Q$
(TMultiDimFit::SetPowerLimit), and a function
$ F_l$ is only accepted if



$\displaystyle Q_l = \sum_{i=1}^{N} \frac{P_{li}}{P_{max,i}} < Q
$


where $ P_{li}$ is the leading power of variable $ x_i$ in function
$ F_l$. (TMultiDimFit::MakeCandidates). So the number of
functions increase with $ Q$ (1, 2 is fine, 5 is way out).






Gram-Schmidt Orthogonalisation





To further reduce the number of functions in the final expression,
only those functions that significantly reduce $ S$ is chosen. What
`significant' means, is chosen by the user, and will be
discussed below (see 2.3).



The functions $ F_l$ are generally not orthogonal, which means one will
have to evaluate all possible $ F_l$'s over all data-points before
finding the most significant [1]. We can, however, do
better then that. By applying the modified Gram-Schmidt
  orthogonalisation algorithm [5] [3] to the
functions $ F_l$, we can evaluate the contribution to the reduction of
$ S$ from each function in turn, and we may delay the actual inversion
of the curvature-matrix
(TMultiDimFit::MakeGramSchmidt).



So we are let to consider an $ M\times L$ matrix 
$ \mathsf{F}$, an
element of which is given by








$\displaystyle f_{jl} = F_j\left(x_{1j} , x_{2j}, \ldots, x_{Nj}\right) = F_l(\mathbf{x}_j) $   with$\displaystyle  j=1,2,\ldots,M,$

(3)





where $ j$ labels the $ M$ rows in the training sample and $ l$ labels
$ L$ functions of $ N$ variables, and $ L \leq M$. That is, $ f_{jl}$ is
the term (or function) numbered $ l$ evaluated at the data point
$ j$. We have to normalise 
$ \mathbf{x}_j$ to $ [-1,1]$ for this to
succeed [5]
(TMultiDimFit::MakeNormalized). We then define a
matrix 
$ \mathsf{W}$ of which the columns 
$ \mathbf{w}_j$ are given by







$\displaystyle \mathbf{w}_1$
$\displaystyle =$
$\displaystyle \mathbf{f}_1 = F_1\left(\mathbf x_1\right)$

(4)


$\displaystyle \mathbf{w}_l$
$\displaystyle =$
$\displaystyle \mathbf{f}_l - \sum^{l-1}_{k=1} \frac{\mathbf{f}_l \bullet
\mathbf{w}_k}{\mathbf{w}_k^2}\mathbf{w}_k .$

(5)





and 
$ \mathbf{w}_{l}$ is the component of 
$ \mathbf{f}_{l}$ orthogonal
to 
$ \mathbf{w}_{1}, \ldots, \mathbf{w}_{l-1}$. Hence we obtain
[3],








$\displaystyle \mathbf{w}_k\bullet\mathbf{w}_l = 0$   if$\displaystyle  k \neq l\quad.$

(6)








We now take as a new model 
$ \mathsf{W}\mathbf{a}$. We thus want to
minimize








$\displaystyle S\equiv \left(\mathbf{D} - \mathsf{W}\mathbf{a}\right)^2\quad,$

(7)





where 
$ \mathbf{D} = \left(D_1,\ldots,D_M\right)$ is a vector of the
dependent quantity in the sample. Differentiation with respect to
$ a_j$ gives, using (6),








$\displaystyle \mathbf{D}\bullet\mathbf{w}_l - a_l\mathbf{w}_l^2 = 0$

(8)





or








$\displaystyle a_l = \frac{\mathbf{D}_l\bullet\mathbf{w}_l}{\mathbf{w}_l^2}$

(9)





Let $ S_j$ be the sum of squares of residuals when taking $ j$ functions
into account. Then








$\displaystyle S_l = \left[\mathbf{D} - \sum^l_{k=1} a_k\mathbf{w}_k\right]^2 = ...
...2 - 2\mathbf{D} \sum^l_{k=1} a_k\mathbf{w}_k + \sum^l_{k=1} a_k^2\mathbf{w}_k^2$

(10)





Using (9), we see that







$\displaystyle S_l$
$\displaystyle =$
$\displaystyle \mathbf{D}^2 - 2 \sum^l_{k=1} a_k^2\mathbf{w}_k^2 +
\sum^j_{k=1} a_k^2\mathbf{w}_k^2$

 


 
$\displaystyle =$
$\displaystyle \mathbf{D}^2 - \sum^l_{k=1} a_k^2\mathbf{w}_k^2$

 


 
$\displaystyle =$
$\displaystyle \mathbf{D}^2 - \sum^l_{k=1} \frac{\left(\mathbf D\bullet \mathbf
w_k\right)}{\mathbf w_k^2}$

(11)








So for each new function $ F_l$ included in the model, we get a
reduction of the sum of squares of residuals of 
$ a_l^2\mathbf{w}_l^2$,
where 
$ \mathbf{w}_l$ is given by (4) and $ a_l$ by
(9). Thus, using the Gram-Schmidt orthogonalisation, we
can decide if we want to include this function in the final model,
before the matrix inversion.








Function Selection Based on Residual





Supposing that $ L-1$ steps of the procedure have been performed, the
problem now is to consider the 
$ L^{\mbox{th}}$ function.



The sum of squares of residuals can be written as








$\displaystyle S_L = \textbf{D}^T\bullet\textbf{D} - \sum^L_{l=1}a^2_l\left(\textbf{w}_l^T\bullet\textbf{w}_l\right)$

(12)





where the relation (9) have been taken into account. The
contribution of the 
$ L^{\mbox{th}}$ function to the reduction of S, is
given by








$\displaystyle \Delta S_L = a^2_L\left(\textbf{w}_L^T\bullet\textbf{w}_L\right)$

(13)








Two test are now applied to decide whether this 
$ L^{\mbox{th}}$
function is to be included in the final expression, or not.








Test 1





Denoting by $ H_{L-1}$ the subspace spanned by

$ \textbf{w}_1,\ldots,\textbf{w}_{L-1}$ the function 
$ \textbf {w}_L$ is
by construction (see (4)) the projection of the function
$ F_L$ onto the direction perpendicular to $ H_{L-1}$. Now, if the
length of 
$ \textbf {w}_L$ (given by 
$ \textbf{w}_L\bullet\textbf{w}_L$)
is very small compared to the length of 
$ \textbf {f}_L$ this new
function can not contribute much to the reduction of the sum of
squares of residuals. The test consists then in calculating the angle
$ \theta $ between the two vectors 
$ \textbf {w}_L$ and 
$ \textbf {f}_L$
(see also figure 1) and requiring that it's
greater then a threshold value which the user must set
(TMultiDimFit::SetMinAngle).










Figure 1:
(a) Angle $ \theta $ between 
$ \textbf {w}_l$ and
      
$ \textbf {f}_L$, (b) angle $ \phi $ between 
$ \textbf {w}_L$ and
      
$ \textbf {D}$

\begin{figure}\begin{center}
\begin{tabular}{p{.4\textwidth}p{.4\textwidth}}
\...
... \put(80,100){$\mathbf{D}$}
\end{picture} \end{tabular} \end{center}\end{figure}










 

Test 2





Let 
$ \textbf {D}$ be the data vector to be fitted. As illustrated in
figure 1, the 
$ L^{\mbox{th}}$ function 
$ \textbf {w}_L$
will contribute significantly to the reduction of $ S$, if the angle

$ \phi^\prime$ between 
$ \textbf {w}_L$ and 
$ \textbf {D}$ is smaller than
an upper limit $ \phi $, defined by the user
(TMultiDimFit::SetMaxAngle)



However, the method automatically readjusts the value of this angle
while fitting is in progress, in order to make the selection criteria
less and less difficult to be fulfilled. The result is that the
functions contributing most to the reduction of $ S$ are chosen first
(TMultiDimFit::TestFunction).



In case $ \phi $ isn't defined, an alternative method of
performing this second test is used: The 
$ L^{\mbox{th}}$ function

$ \textbf {f}_L$ is accepted if (refer also to equation (13))








$\displaystyle \Delta S_L > \frac{S_{L-1}}{L_{max}-L}$

(14)





where  $ S_{L-1}$ is the sum of the $ L-1$ first residuals from the
$ L-1$ functions previously accepted; and $ L_{max}$ is the total number
of functions allowed in the final expression of the fit (defined by
user).



From this we see, that by restricting $ L_{max}$ -- the number of
terms in the final model -- the fit is more difficult to perform,
since the above selection criteria is more limiting.



The more coefficients we evaluate, the more the sum of squares of
residuals $ S$ will be reduced. We can evaluate $ S$ before inverting

$ \mathsf{B}$ as shown below.






Coefficients and Coefficient Errors





Having found a parameterization, that is the $ F_l$'s and $ L$, that
minimizes $ S$, we still need to determine the coefficients
$ c_l$. However, it's a feature of how we choose the significant
functions, that the evaluation of the $ c_l$'s becomes trivial
[5]. To derive 
$ \mathbf{c}$, we first note that
equation (4) can be written as








$\displaystyle \mathsf{F} = \mathsf{W}\mathsf{B}$

(15)





where








$\displaystyle b_{ij} = \left\{\begin{array}{rcl} \frac{\mathbf{f}_j \bullet \ma...
...f} & i < j\  1 & \mbox{if} & i = j\  0 & \mbox{if} & i > j \end{array}\right.$

(16)





Consequently, 
$ \mathsf{B}$ is an upper triangle matrix, which can be
readily inverted. So we now evaluate








$\displaystyle \mathsf{F}\mathsf{B}^{-1} = \mathsf{W}$

(17)





The model 
$ \mathsf{W}\mathbf{a}$ can therefore be written as



$\displaystyle (\mathsf{F}\mathsf{B}^{-1})\mathbf{a} =
\mathsf{F}(\mathsf{B}^{-1}\mathbf{a}) .
$


The original model 
$ \mathsf{F}\mathbf{c}$ is therefore identical with
this if








$\displaystyle \mathbf{c} = \left(\mathsf{B}^{-1}\mathbf{a}\right) = \left[\mathbf{a}^T\left(\mathsf{B}^{-1}\right)^T\right]^T .$

(18)





The reason we use 
$ \left(\mathsf{B}^{-1}\right)^T$ rather then

$ \mathsf{B}^{-1}$ is to save storage, since

$ \left(\mathsf{B}^{-1}\right)^T$ can be stored in the same matrix as

$ \mathsf{B}$
(TMultiDimFit::MakeCoefficients). The errors in
the coefficients is calculated by inverting the curvature matrix
of the non-orthogonal functions $ f_{lj}$ [1]
(TMultiDimFit::MakeCoefficientErrors).








Considerations





It's important to realize that the training sample should be
representive of the problem at hand, in particular along the borders
of the region of interest. This is because the algorithm presented
here, is a interpolation, rahter then a extrapolation
[5].



Also, the independent variables $ x_{i}$ need to be linear
independent, since the procedure will perform poorly if they are not
[5]. One can find an linear transformation from ones
original variables $ \xi_{i}$ to a set of linear independent variables
$ x_{i}$, using a Principal Components Analysis
(see TPrincipal), and
then use the transformed variable as input to this class [5]
[6].



H. Wind also outlines a method for parameterising a multidimensional
dependence over a multidimensional set of variables. An example
of the method from [5], is a follows (please refer to
[5] for a full discussion):





    

  1. Define 
$ \mathbf{P} = (P_1, \ldots, P_5)$ are the 5 dependent
  quantities that define a track.

    2. 

  2. Compute, for $ M$ different values of 
$ \mathbf{P}$, the tracks
  through the magnetic field, and determine the corresponding
  
$ \mathbf{x} = (x_1, \ldots, x_N)$.

    3. 

  3. Use the simulated observations to determine, with a simple
  approximation, the values of 
$ \mathbf{P}_j$. We call these values
  
$ \mathbf{P}^\prime_j, j = 1, \ldots, M$.

    4. 

  4. Determine from 
$ \mathbf{x}$ a set of at least five relevant
  coordinates 
$ \mathbf{x}^\prime$, using contrains, or
    alternative:

    5. 

  5. Perform a Principal Component Analysis (using
  TPrincipal), and use

to get a linear transformation
  
$ \mathbf{x} \rightarrow \mathbf{x}^\prime$, so that
  
$ \mathbf{x}^\prime$ are constrained and linear independent.

    6. 

  6. Perform a Principal Component Analysis on
  
$ Q_i = P_i / P^prime_i  i = 1, \ldots, 5$, to get linear
  indenpendent (among themselves, but not independent of
  
$ \mathbf{x}$) quantities 
$ \mathbf{Q}^\prime$

    7. 

  7. For each component 
$ Q^\prime_i$ make a mutlidimensional fit,
  using 
$ \mathbf{x}^\prime$ as the variables, thus determing a set of
  coefficents 
$ \mathbf{c}_i$.

    8. 





To process data, using this parameterisation, do


    

  1. Test wether the observation 
$ \mathbf{x}$ within the domain of
  the parameterization, using the result from the Principal Component
  Analysis.

    2. 

  2. Determine 
$ \mathbf{P}^\prime$ as before.

    3. 

  3. Detetmine 
$ \mathbf{x}^\prime$ as before.

    4. 

  4. Use the result of the fit to determind 
$ \mathbf{Q}^\prime$.

    5. 

  5. Transform back to 
$ \mathbf{P}$ from 
$ \mathbf{Q}^\prime$, using
  the result from the Principal Component Analysis.

    6. 










Testing the parameterization





The class also provides functionality for testing the, over the
training sample, found parameterization
(TMultiDimFit::Fit). This is done by passing
the class a test sample of $ M_t$ tuples of the form 
$ (\mathbf{x}_{t,j},
D_{t,j}, E_{t,j})$, where 
$ \mathbf{x}_{t,j}$ are the independent
variables, $ D_{t,j}$ the known, dependent quantity, and $ E_{t,j}$ is
the square error in $ D_{t,j}$
(TMultiDimFit::AddTestRow).



The parameterization is then evaluated at every 
$ \mathbf{x}_t$ in the
test sample, and



$\displaystyle S_t \equiv \sum_{j=1}^{M_t} \left(D_{t,j} -
D_p\left(\mathbf{x}_{t,j}\right)\right)^2
$


is evaluated. The relative error over the test sample



$\displaystyle R_t = \frac{S_t}{\sum_{j=1}^{M_t} D_{t,j}^2}
$


should not be to low or high compared to $ R$ from the training
sample. Also, multiple correlation coefficient from both samples should
be fairly close, otherwise one of the samples is not representive of
the problem. A large difference in the reduced $ \chi^2$ over the two
samples indicate an over fit, and the maximum number of terms in the
parameterisation should be reduced.



It's possible to use Minuit
[4] to further improve the fit, using the test sample.





Christian Holm

  November 2000, NBI









Bibliography

1


Philip R. Bevington and D. Keith Robinson.

Data Reduction and Error Analysis for the Physical Sciences.

McGraw-Hill, 2 edition, 1992.


2


René Brun et al.

Mudifi.

Long writeup DD/75-23, CERN, 1980.


3


Gene H. Golub and Charles F. van Loan.

Matrix Computations.

John Hopkins Univeristy Press, Baltimore, 3 edition, 1996.


4


F. James.

Minuit.

Long writeup D506, CERN, 1998.


5


H. Wind.

Function parameterization.

In Proceedings of the 1972 CERN Computing and Data Processing
  School, volume 72-21 of Yellow report. CERN, 1972.


6


H. Wind.

1. principal component analysis, 2. pattern recognition for track
  finding, 3. interpolation and functional representation.

Yellow report EP/81-12, CERN, 1981.



 */
//
//

#include "Riostream.h"
#include "TMultiDimFit.h"
#include "TMath.h"
#include "TH1.h"
#include "TH2.h"
#include "TROOT.h"
#include "TBrowser.h"

#define RADDEG (180. / TMath::Pi())
#define DEGRAD (TMath::Pi() / 180.)
#define HIST_XORIG     0
#define HIST_DORIG     1
#define HIST_XNORM     2
#define HIST_DSHIF     3
#define HIST_RX        4
#define HIST_RD        5
#define HIST_RTRAI     6
#define HIST_RTEST     7
#define PARAM_MAXSTUDY 1
#define PARAM_SEVERAL  2
#define PARAM_RELERR   3
#define PARAM_MAXTERMS 4


//____________________________________________________________________
static void mdfHelper(int&, double*, double&, double*, int);

//____________________________________________________________________
ClassImp(TMultiDimFit);

//____________________________________________________________________
// Static instance. Used with mdfHelper and TMinuit
TMultiDimFit* TMultiDimFit::fgInstance = 0;


//____________________________________________________________________
 TMultiDimFit::TMultiDimFit()
{
  // Empty CTOR. Do not use
  fMeanQuantity            = 0;
  fMaxQuantity             = 0;
  fMinQuantity             = 0;
  fSumSqQuantity           = 0;
  fSumSqAvgQuantity        = 0;
  fPowerLimit              = 1;

  fMaxAngle                = 0;
  fMinAngle                = 1;

  fNVariables              = 0;
  fMaxVariables            = 0;
  fMinVariables            = 0;
  fSampleSize              = 0;

  fMaxAngle                = 0;
  fMinAngle                = 0;

  fPolyType                = kMonomials;
  fShowCorrelation         = kFALSE;

  fIsUserFunction          = kFALSE;

  fPowers                  = 0;
  fMaxPowers               = 0;
  fMaxPowersFinal          = 0;

  fHistograms              = 0;
  fHistogramMask           = 0;

  fQuantity(0);
  fVariables(0);
  fMaxVariables(0);
  fMinVariables(0);
  fMeanVariables(0);

  fFitter                  = 0;
}


//____________________________________________________________________
 TMultiDimFit::TMultiDimFit(Int_t dimension,
			   EMDFPolyType type,
			   Option_t *option)
  : TNamed("multidimfit","Multi-dimensional fit object"),
    fQuantity(dimension),
    fSqError(dimension),
    fVariables(dimension*100),
    fMeanVariables(dimension),
    fMaxVariables(dimension),
    fMinVariables(dimension)
{
  // Constructor
  // Second argument is the type of polynomials to use in
  // parameterisation, one of:
  //      TMultiDimFit::kMonomials
  //      TMultiDimFit::kChebyshev
  //      TMultiDimFit::kLegendre
  //
  // Options:
  //   K      Compute (k)correlation matrix
  //   V      Be verbose
  //
  // Default is no options.
  //

  fgInstance = this;

  fMeanQuantity           = 0;
  fMaxQuantity            = 0;
  fMinQuantity            = 0;
  fSumSqQuantity          = 0;
  fSumSqAvgQuantity       = 0;
  fPowerLimit             = 1;

  fMaxAngle               = 0;
  fMinAngle               = 1;

  fNVariables             = dimension;
  fMaxVariables           = 0;
  fMinVariables           = 0;
  fSampleSize             = 0;
  fMinRelativeError       = 0.01;
  fError                  = 0;
  fTestError              = 0;
  fPrecision              = 0;
  fTestPrecision          = 0;
  fParameterisationCode   = 0;

  fPolyType               = type;
  fShowCorrelation        = kFALSE;
  fIsVerbose              = kFALSE;

  TString opt             = option;
  opt.ToLower();

  if (opt.Contains("k")) fShowCorrelation = kTRUE;
  if (opt.Contains("v")) fIsVerbose       = kTRUE;

  fIsUserFunction         = kFALSE;

  fHistograms             = 0;
  fHistogramMask          = 0;

  fPowers                 = 0;
  fMaxPowers              = new Int_t[dimension];
  fMaxPowersFinal         = new Int_t[dimension];
  fFitter                 = 0;
}


//____________________________________________________________________
 TMultiDimFit::~TMultiDimFit()
{
  // DTOR
  delete [] fPowers;
  delete [] fMaxPowers;
  delete [] fMaxPowersFinal;
  delete [] fPowerIndex;
  delete [] fFunctionCodes;
  delete fHistograms;
}


//____________________________________________________________________
 void TMultiDimFit::AddRow(const Double_t *x, Double_t D, Double_t E)
{
  // Add a row consisting of fNVariables independent variables, the
  // known, dependent quantity, and optionally, the square error in
  // the dependent quantity, to the training sample to be used for the
  // parameterization.
  // The mean of the variables and quantity is calculated on the fly,
  // as outlined in TPrincipal::AddRow.
  // This sample should be representive of the problem at hand.
  // Please note, that if no error is given Poisson statistics is
  // assumed and the square error is set to the value of dependent
  // quantity.  See also the
  // class description
  if (!x)
    return;

  if (++fSampleSize == 1) {
    fMeanQuantity  = D;
    fMaxQuantity   = D;
    fMinQuantity   = D;
  }
  else {
    fMeanQuantity  *= 1 - 1./Double_t(fSampleSize);
    fMeanQuantity  += D / Double_t(fSampleSize);
    fSumSqQuantity += D * D;

    if (D >= fMaxQuantity) fMaxQuantity = D;
    if (D <= fMinQuantity) fMinQuantity = D;
  }


  // If the vector isn't big enough to hold the new data, then
  // expand the vector by half it's size.
  Int_t size = fQuantity.GetNrows();
  if (fSampleSize > size) {
    fQuantity.ResizeTo(size + size/2);
    fSqError.ResizeTo(size + size/2);
  }

  // Store the value
  fQuantity(fSampleSize-1) = D;
  fSqError(fSampleSize-1) = (E == 0 ? D : E);

  // Store data point in internal vector
  // If the vector isn't big enough to hold the new data, then
  // expand the vector by half it's size
  size = fVariables.GetNrows();
  if (fSampleSize * fNVariables > size)
    fVariables.ResizeTo(size + size/2);


   // Increment the data point counter
  Int_t i,j;
  for (i = 0; i < fNVariables; i++) {
    if (fSampleSize == 1) {
      fMeanVariables(i) = x[i];
      fMaxVariables(i)  = x[i];
      fMinVariables(i)  = x[i];
    }
    else {
      fMeanVariables(i) *= 1 - 1./Double_t(fSampleSize);
      fMeanVariables(i) += x[i] / Double_t(fSampleSize);

      // Update the maximum value for this component
      if (x[i] >= fMaxVariables(i)) fMaxVariables(i)  = x[i];

      // Update the minimum value for this component
      if (x[i] <= fMinVariables(i)) fMinVariables(i)  = x[i];

    }

    // Store the data.
    j = (fSampleSize-1) * fNVariables + i;
    fVariables(j) = x[i];
  }
}


//____________________________________________________________________
 void TMultiDimFit::AddTestRow(const Double_t *x, Double_t D, Double_t E)
{
  // Add a row consisting of fNVariables independent variables, the
  // known, dependent quantity, and optionally, the square error in
  // the dependent quantity, to the test sample to be used for the
  // test of the parameterization.
  // This sample needn't be representive of the problem at hand.
  // Please note, that if no error is given Poisson statistics is
  // assumed and the square error is set to the value of dependent
  // quantity.  See also the
  // class description
  if (fTestSampleSize++ == 0) {
    fTestQuantity.ResizeTo(fNVariables);
    fTestSqError.ResizeTo(fNVariables);
    fTestVariables.ResizeTo(fNVariables * 100);
  }

  // If the vector isn't big enough to hold the new data, then
  // expand the vector by half it's size.
  Int_t size = fTestQuantity.GetNrows();
  if (fTestSampleSize > size) {
    fTestQuantity.ResizeTo(size + size/2);
    fTestSqError.ResizeTo(size + size/2);
  }

  // Store the value
  fTestQuantity(fTestSampleSize-1) = D;
  fTestSqError(fTestSampleSize-1) = (E == 0 ? D : E);

  // Store data point in internal vector
  // If the vector isn't big enough to hold the new data, then
  // expand the vector by half it's size
  size = fTestVariables.GetNrows();
  if (fTestSampleSize * fNVariables > size)
    fTestVariables.ResizeTo(size + size/2);


   // Increment the data point counter
  Int_t i,j;
  for (i = 0; i < fNVariables; i++) {
    j = fNVariables * (fTestSampleSize - 1) + i;
    fTestVariables(j) = x[i];

    if (x[i] > fMaxVariables(i))
      Warning("AddTestRow", "variable %d (row: %d) too large: %f > %f",
	       i, fTestSampleSize, x[i], fMaxVariables(i));
    if (x[i] < fMinVariables(i))
      Warning("AddTestRow", "variable %d (row: %d) too small: %f < %f",
	       i, fTestSampleSize, x[i], fMinVariables(i));
  }
}


//____________________________________________________________________
 void TMultiDimFit::Browse(TBrowser* b)
{
  // Browse the TMultiDimFit object in the TBrowser.
  if (fHistograms) {
    TIter next(fHistograms);
    TH1* h = 0;
    while ((h = (TH1*)next()))
      b->Add(h,h->GetName());
  }
  if (fVariables.IsValid())
    b->Add(&fVariables, "Variables (Training)");
  if (fQuantity.IsValid())
    b->Add(&fQuantity, "Quantity (Training)");
  if (fSqError.IsValid())
    b->Add(&fSqError, "Error (Training)");
  if (fMeanVariables.IsValid())
    b->Add(&fMeanVariables, "Mean of Variables (Training)");
  if (fMaxVariables.IsValid())
    b->Add(&fMaxVariables, "Mean of Variables (Training)");
  if (fMinVariables.IsValid())
    b->Add(&fMinVariables, "Min of Variables (Training)");
  if (fTestVariables.IsValid())
    b->Add(&fTestVariables, "Variables (Test)");
  if (fTestQuantity.IsValid())
    b->Add(&fTestQuantity, "Quantity (Test)");
  if (fTestSqError.IsValid())
    b->Add(&fTestSqError, "Error (Test)");
  if (fFunctions.IsValid())
    b->Add(&fFunctions, "Functions");
  if(fCoefficients.IsValid())
    b->Add(&fCoefficients,"Coefficients");
  if(fCoefficientsRMS.IsValid())
    b->Add(&fCoefficientsRMS,"Coefficients Errors");
  if (fOrthFunctions.IsValid())
    b->Add(&fOrthFunctions, "Orthogonal Functions");
  if (fOrthFunctionNorms.IsValid())
    b->Add(&fOrthFunctionNorms, "Orthogonal Functions Norms");
  if (fResiduals.IsValid())
    b->Add(&fResiduals, "Residuals");
  if(fOrthCoefficients.IsValid())
    b->Add(&fOrthCoefficients,"Orthogonal Coefficients");
  if (fOrthCurvatureMatrix.IsValid())
    b->Add(&fOrthCurvatureMatrix,"Orthogonal curvature matrix");
  if(fCorrelationMatrix.IsValid())
    b->Add(&fCorrelationMatrix,"Correlation Matrix");
  if (fFitter)
    b->Add(fFitter, fFitter->GetName());
}


//____________________________________________________________________
 void TMultiDimFit::Clear(Option_t *option)
{
  // Clear internal structures and variables
  Int_t i, j, n = fNVariables, m = fMaxFunctions;

  // Training sample, dependent quantity
  fQuantity.Zero();
  fSqError.Zero();
  fMeanQuantity                 = 0;
  fMaxQuantity                  = 0;
  fMinQuantity                  = 0;
  fSumSqQuantity                = 0;
  fSumSqAvgQuantity             = 0;

  // Training sample, independent variables
  fVariables.Zero();
  fNVariables                   = 0;
  fSampleSize                   = 0;
  fMeanVariables.Zero();
  fMaxVariables.Zero();
  fMinVariables.Zero();

  // Test sample
  fTestQuantity.Zero();
  fTestSqError.Zero();
  fTestVariables.Zero();
  fTestSampleSize               = 0;

  // Functions
  fFunctions.Zero();
  for (i = 0; i < fMaxTerms; i++) {
    fFunctionCodes[i]           = 0;
    fPowerIndex[i]              = 0;
  }
  fMaxFunctions                 = 0;
  fMaxStudy                     = 0;
  fOrthFunctions.Zero();
  fOrthFunctionNorms.Zero();

  // Control parameters
  fMinRelativeError             = 0;
  fMinAngle                     = 0;
  fMaxAngle                     = 0;
  fMaxTerms                     = 0;

  // Powers
  for (i = 0; i < n; i++) {
    fMaxPowers[i]               = 0;
    fMaxPowersFinal[i]          = 0;
    for (j = 0; j < m; j++)
      fPowers[i * n + j]        = 0;
  }
  fPowerLimit                   = 0;

  // Residuals
  fMaxResidual                  = 0;
  fMinResidual                  = 0;
  fMaxResidualRow               = 0;
  fMinResidualRow               = 0;
  fSumSqResidual                = 0;

  // Fit
  fNCoefficients                = 0;
  fOrthCoefficients             = 0;
  fOrthCurvatureMatrix          = 0;
  fRMS                          = 0;
  fCorrelationMatrix.Zero();
  fError                        = 0;
  fTestError                    = 0;
  fPrecision                    = 0;
  fTestPrecision                = 0;

  // Coefficients
  fCoefficients.Zero();
  fCoefficientsRMS.Zero();
  fResiduals.Zero();
  fHistograms->Clear(option);

  // Options
  fPolyType                     = kMonomials;
  fShowCorrelation              = kFALSE;
  fIsUserFunction               = kFALSE;
}


//____________________________________________________________________
 Double_t TMultiDimFit::Eval(const Double_t *x, const Double_t* coeff)
{
  // Evaluate parameterization at point x. Optional argument coeff is
  // a vector of coefficients for the parameterisation, fNCoefficients
  // elements long.
  Double_t returnValue = fMeanQuantity;
  Double_t term        = 0;
  Int_t    i, j;

  for (i = 0; i < fNCoefficients; i++) {
    // Evaluate the ith term in the expansion
    term = (coeff ? coeff[i] : fCoefficients(i));
    for (j = 0; j < fNVariables; j++) {
      // Evaluate the factor (polynomial) in the j-th variable.
      Int_t    p  =  fPowers[fPowerIndex[i] * fNVariables + j];
      Double_t X  =  1 + 2. / (fMaxVariables(j) - fMinVariables(j))
	* (x[j] - fMaxVariables(j));
      term        *= EvalFactor(p,X);
    }
    // Add this term to the final result
    returnValue += term;
  }
  return returnValue;
}


//____________________________________________________________________
 Double_t TMultiDimFit::EvalControl(const Int_t *iv)
{
  // PRIVATE METHOD:
  // Calculate the control parameter from the passed powers
  Double_t s = 0;
  Double_t epsilon = 1e-6; // a small number
  for (Int_t i = 0; i < fNVariables; i++) {
    if (fMaxPowers[i] != 1)
      s += (epsilon + iv[i] - 1) / (epsilon + fMaxPowers[i] - 1);
  }
  return s;
}

//____________________________________________________________________
 Double_t TMultiDimFit::EvalFactor(Int_t p, Double_t x)
{
  // PRIVATE METHOD:
  // Evaluate function with power p at variable value x
  Int_t    i   = 0;
  Double_t p1  = 1;
  Double_t p2  = 0;
  Double_t p3  = 0;
  Double_t r   = 0;

  switch(p) {
  case 1:
    r = 1;
    break;
  case 2:
    r =  x;
    break;
  default:
    p2 = x;
    for (i = 3; i <= p; i++) {
      p3 = p2 * x;
      if (fPolyType == kLegendre)
	p3 = ((2 * i - 3) * p2 * x - (i - 2) * p1) / (i - 1);
      else if (fPolyType == kChebyshev)
	p3 = 2 * x * p2 - p1;
      p1 = p2;
      p2 = p3;
    }
    r = p3;
  }

  return r;
}


//____________________________________________________________________
 void TMultiDimFit::FindParameterization(Option_t *)
{
  // Find the parameterization
  //
  // Options:
  //     None so far
  //
  // For detailed description of what this entails, please refer to the
  // class description
  MakeNormalized();
  MakeCandidates();
  MakeParameterization();
  MakeCoefficients();
  MakeCoefficientErrors();
  MakeCorrelation();
}

//____________________________________________________________________
 void TMultiDimFit::Fit(Option_t *option)
{
  // Try to fit the found parameterisation to the test sample.
  //
  // Options
  //     M     use Minuit to improve coefficients
  //
  // Also, refer to
  // class description
  Int_t i, j;
  Double_t*      x    = new Double_t[fNVariables];
  Double_t  sumSqD    = 0;
  Double_t    sumD    = 0;
  Double_t  sumSqR    = 0;
  Double_t    sumR    = 0;

  // Calculate the residuals over the test sample
  for (i = 0; i < fTestSampleSize; i++) {
    for (j = 0; j < fNVariables; j++)
      x[j] = fTestVariables(i * fNVariables + j);
    Double_t res =  fTestQuantity(i) - Eval(x);
    sumD         += fTestQuantity(i);
    sumSqD       += fTestQuantity(i) * fTestQuantity(i);
    sumR         += res;
    sumSqR       += res * res;
    if (TESTBIT(fHistogramMask,HIST_RTEST))
      ((TH1F*)fHistograms->FindObject("res_test"))->Fill(res);
  }
  Double_t dAvg         = sumSqD - (sumD * sumD) / fTestSampleSize;
  Double_t rAvg         = sumSqR - (sumR * sumR) / fTestSampleSize;
  fTestCorrelationCoeff = (dAvg - rAvg) / dAvg;
  fTestError            = sumSqR;
  fTestPrecision        = sumSqR / sumSqD;

  TString opt(option);
  opt.ToLower();

  if (!opt.Contains("m"))
    MakeChi2();

  if (fNCoefficients * 50 > fTestSampleSize)
    Warning("Fit", "test sample is very small");

  if (!opt.Contains("m"))
    return;

  fFitter = TVirtualFitter::Fitter(0,fNCoefficients);
  fFitter->SetFCN(mdfHelper);

  const Int_t  maxArgs = 16;
  Int_t           args = 1;
  Double_t*   arglist  = new Double_t[maxArgs];
  arglist[0]           = -1;
  fFitter->ExecuteCommand("SET PRINT",arglist,args);

  for (i = 0; i < fNCoefficients; i++) {
    Double_t startVal = fCoefficients(i);
    Double_t startErr = fCoefficientsRMS(i);
    fFitter->SetParameter(i, Form("coeff%02d",i),
			  startVal, startErr, 0, 0);
  }

  // arglist[0]           = 0;
  args                 = 1;
  // fFitter->ExecuteCommand("SET PRINT",arglist,args);
  fFitter->ExecuteCommand("MIGRAD",arglist,args);

  for (i = 0; i < fNCoefficients; i++) {
    Double_t val = 0, err = 0, low = 0, high = 0;
    fFitter->GetParameter(i, Form("coeff%02d",i),
			 val, err, low, high);
    fCoefficients(i)    = val;
    fCoefficientsRMS(i) = err;
  }
}


//____________________________________________________________________
 void TMultiDimFit::MakeCandidates()
{
  // PRIVATE METHOD:
  // Create list of candidate functions for the parameterisation. See
  // also
  // class description
  Int_t i = 0;
  Int_t j = 0;
  Int_t k = 0;

  // The temporary array to store the powers in. We don't need to
  // initialize this array however.
  Int_t *powers = new Int_t[fNVariables * fMaxFunctions];

  // store of `control variables'
  Double_t* control  = new Double_t[fMaxFunctions];

  // We've better initialize the variables
  Int_t *iv = new Int_t[fNVariables];
  for (i = 0; i < fNVariables; i++)
    iv[i] = 1;

  if (!fIsUserFunction) {

    // Number of funcs selected
    Int_t     numberFunctions = 0;

    // Absolute max number of functions
    Int_t maxNumberFunctions = 1;
    for (i = 0; i < fNVariables; i++)
      maxNumberFunctions *= fMaxPowers[i];

    while (kTRUE) {
      // Get the control value for this function
      Double_t s = EvalControl(iv);

      if (s <= fPowerLimit) {

	// Call over-loadable method Select, as to allow the user to
	// interfere with the selection of functions.
	if (Select(iv)) {
	  numberFunctions++;

	  // If we've reached the user defined limit of how many
	  // functions we can consider, break out of the loop
	  if (numberFunctions > fMaxFunctions)
	  break;

	  // Store the control value, so we can sort array of powers
	  // later on
	  control[numberFunctions-1] = s;

	  // Store the powers in powers array.
	  for (i = 0; i < fNVariables; i++) {
	    j = (numberFunctions - 1) * fNVariables + i;
	    powers[j] = iv[i];
	  }
	} // if (Select())
      } // if (s <= fPowerLimit)

      for (i = 0; i < fNVariables; i++)
	if (iv[i] < fMaxPowers[i])
	  break;

      // If all variables have reached their maximum power, then we
      // break out of the loop
      if (i == fNVariables) {
	fMaxFunctions = numberFunctions;
	break;
      }

      // Next power in variable i
      iv[i]++;

      for (j = 0; j < i; j++)
	iv[j] = 1;
    } // while (kTRUE)
  }
  else {
    // In case the user gave an explicit function
    for (i = 0; i < fMaxFunctions; i++) {
      // Copy the powers to working arrays
      for (j = 0; j < fNVariables; j++) {
	powers[i * fNVariables + j] = fPowers[i * fNVariables + j];
	iv[j]                 = fPowers[i * fNVariables + j];
      }

      control[i] = EvalControl(iv);
    }
  }

  // Now we need to sort the powers according to least `control
  // variable'
  Int_t *order = new Int_t[fMaxFunctions];
  for (i = 0; i < fMaxFunctions; i++)
    order[i] = i;
  fPowers = new Int_t[fMaxFunctions * fNVariables];

  for (i = 0; i < fMaxFunctions; i++) {
    Double_t x = control[i];
    Int_t    l = order[i];
    k = i;

    for (j = i; j < fMaxFunctions; j++) {
      if (control[j] <= x) {
	x = control[j];
	l = order[j];
	k = j;
      }
    }

    if (k != i) {
      control[k] = control[i];
      control[i] = x;
      order[k]   = order[i];
      order[i]   = l;
    }
  }

  for (i = 0; i < fMaxFunctions; i++)
    for (j = 0; j < fNVariables; j++)
      fPowers[i * fNVariables + j] = powers[order[i] * fNVariables + j];

  delete [] control;
  delete [] powers;
  delete [] order;
  delete [] iv;
}


//____________________________________________________________________
 Double_t TMultiDimFit::MakeChi2(const Double_t* coeff)
{
  // Calculate Chi square over either the test sample. The optional
  // argument coeff is a vector of coefficients to use in the
  // evaluation of the parameterisation. If coeff == 0, then the found
  // coefficients is used.
  // Used my MINUIT for fit (see TMultDimFit::Fit)
  fChi2 = 0;
  Int_t i, j;
  Double_t* x = new Double_t[fNVariables];
  for (i = 0; i < fTestSampleSize; i++) {
    // Get the stored point
    for (j = 0; j < fNVariables; j++)
      x[j] = fTestVariables(i * fNVariables + j);

    // Evaluate function. Scale to shifted values
    Double_t f = Eval(x,coeff);

    // Calculate contribution to Chic square
    fChi2 += 1. / TMath::Max(fTestSqError(i),1e-20)
      * (fTestQuantity(i) - f) * (fTestQuantity(i) - f);
  }

  // Clean up
  delete x;

  return fChi2;
}


//____________________________________________________________________
 void TMultiDimFit::MakeCode(const char* filename, Option_t *option)
{
  // Generate the file <filename> with .C appended if argument doesn't
  // end in .cxx or .C. The contains the implementation of the
  // function:
  //
  //   Double_t <funcname>(Double_t *x)
  //
  // which does the same as TMultiDimFit::Eval. Please refer to this
  // method.
  //
  // Further, the static variables:
  //
  //     Int_t    gNVariables
  //     Int_t    gNCoefficients
  //     Double_t gDMean
  //     Double_t gXMean[]
  //     Double_t gXMin[]
  //     Double_t gXMax[]
  //     Double_t gCoefficient[]
  //     Int_t    gPower[]
  //
  // are initialized. The only ROOT header file needed is Rtypes.h
  //
  // See TMultiDimFit::MakeRealCode for a list of options


  TString outName(filename);
  if (!outName.EndsWith(".C") && !outName.EndsWith(".cxx"))
    outName += ".C";

  MakeRealCode(outName.Data(),"",option);
}



//____________________________________________________________________
 void TMultiDimFit::MakeCoefficientErrors()
{
  // PRIVATE METHOD:
  // Compute the errors on the coefficients. For this to be done, the
  // curvature matrix of the non-orthogonal functions, is computed.
  Int_t    i = 0;
  Int_t    j = 0;
  Int_t    k = 0;
  TVectorD iF(fSampleSize);
  TVectorD jF(fSampleSize);
  fCoefficientsRMS.ResizeTo(fNCoefficients);

  TMatrixD curvatureMatrix(fNCoefficients,fNCoefficients);

  // Build the curvature matrix
  for (i = 0; i < fNCoefficients; i++) {
    iF = TMatrixDRow(fFunctions,i);
    for (j = 0; j <= i; j++) {
      jF = TMatrixDRow(fFunctions,j);
      for (k = 0; k < fSampleSize; k++)
	curvatureMatrix(i,j) +=
	  1 / TMath::Max(fSqError(k), 1e-20) * iF(k) * jF(k);
      curvatureMatrix(j,i) = curvatureMatrix(i,j);
    }
  }

  // Calculate Chi Square
  fChi2 = 0;
  for (i = 0; i < fSampleSize; i++) {
    Double_t f = 0;
    for (j = 0; j < fNCoefficients; j++)
      f += fCoefficients(j) * fFunctions(j,i);
    fChi2 += 1. / TMath::Max(fSqError(i),1e-20) * (fQuantity(i) - f)
      * (fQuantity(i) - f);
  }

  // Invert the curvature matrix
  Double_t det = 1;
  curvatureMatrix.Invert(&det);
  if (det == 0)
    Warning("MakeCoefficientErrors", "curvature matrix is singular");

  if (fIsVerbose)
    cout << "Determinant of curvature matrix: " << det << endl;

  for (i = 0; i < fNCoefficients; i++)
    fCoefficientsRMS(i) = TMath::Sqrt(curvatureMatrix(i,i));
}


//____________________________________________________________________
 void TMultiDimFit::MakeCoefficients()
{
  // PRIVATE METHOD:
  // Invert the model matrix B, and compute final coefficients. For a
  // more thorough discussion of what this means, please refer to the
  // class description
  //
  // First we invert the lower triangle matrix fOrthCurvatureMatrix
  // and store the inverted matrix in the upper triangle.

  Int_t i = 0, j = 0;
  Int_t col = 0, row = 0;

  // Invert the B matrix
  for (col = 1; col < fNCoefficients; col++) {
    for (row = col - 1; row > -1; row--) {
      fOrthCurvatureMatrix(row,col) = 0;
      for (i = row; i <= col ; i++)
 	fOrthCurvatureMatrix(row,col) -=
	  fOrthCurvatureMatrix(i,row)
	  * fOrthCurvatureMatrix(i,col);
    }
  }

  // Compute the final coefficients
  fCoefficients.ResizeTo(fNCoefficients);

  for (i = 0; i < fNCoefficients; i++) {
    Double_t sum = 0;
    for (j = i; j < fNCoefficients; j++)
      sum += fOrthCurvatureMatrix(i,j) * fOrthCoefficients(j);
    fCoefficients(i) = sum;
  }

  // Compute the final residuals
  fResiduals.ResizeTo(fSampleSize);
  for (i = 0; i < fSampleSize; i++)
    fResiduals(i) = fQuantity(i);

  for (i = 0; i < fNCoefficients; i++)
    for (j = 0; j < fSampleSize; j++)
      fResiduals(j) -= fCoefficients(i) * fFunctions(i,j);

  // Compute the max and minimum, and squared sum of the evaluated
  // residuals
  fMinResidual = 10e10;
  fMaxResidual = -10e10;
  Double_t sqRes  = 0;
  for (i = 0; i < fSampleSize; i++){
    sqRes += fResiduals(i) * fResiduals(i);
    if (fResiduals(i) <= fMinResidual) {
      fMinResidual     = fResiduals(i);
      fMinResidualRow  = i;
    }
    if (fResiduals(i) >= fMaxResidual) {
      fMaxResidual     = fResiduals(i);
      fMaxResidualRow  = i;
    }
  }

  fCorrelationCoeff = fSumSqResidual / fSumSqAvgQuantity;
  fPrecision        = TMath::Sqrt(sqRes / fSumSqQuantity);

  // If we use histograms, fill some more
  if (TESTBIT(fHistogramMask,HIST_RD) ||
      TESTBIT(fHistogramMask,HIST_RTRAI) ||
      TESTBIT(fHistogramMask,HIST_RX)) {
    for (i = 0; i < fSampleSize; i++) {
      if (TESTBIT(fHistogramMask,HIST_RD))
	((TH2D*)fHistograms->FindObject("res_d"))->Fill(fQuantity(i),
							fResiduals(i));
      if (TESTBIT(fHistogramMask,HIST_RTRAI))
	((TH1D*)fHistograms->FindObject("res_train"))->Fill(fResiduals(i));

      if (TESTBIT(fHistogramMask,HIST_RX))
	for (j = 0; j < fNVariables; j++)
	  ((TH2D*)fHistograms->FindObject(Form("res_x_%d",j)))
	    ->Fill(fVariables(i * fNVariables + j),fResiduals(i));
    }
  } // If histograms

}


//____________________________________________________________________
 void TMultiDimFit::MakeCorrelation()
{
  // PRIVATE METHOD:
  // Compute the correlation matrix
  if (!fShowCorrelation)
    return;

  fCorrelationMatrix.ResizeTo(fNVariables,fNVariables+1);

  Double_t d2      = 0;
  Double_t DdotXi  = 0;
  Double_t XiNorm  = 0;
  Double_t XidotXj = 0;
  Double_t XjNorm  = 0;

  Int_t i, j, k, l;
  for (i = 0; i < fSampleSize; i++)
    d2 += fQuantity(i) * fQuantity(i);

  for (i = 0; i < fNVariables; i++) {
    for (j = 0; j< fSampleSize; j++) {
      // Index of sample j of variable i
      k =  j * fNVariables + i;
      DdotXi += fQuantity(j) * (fVariables(k) - fMeanVariables(i));
      XiNorm += (fVariables(k) - fMeanVariables(i))
	* (fVariables(k) - fMeanVariables(i));
    }
    fCorrelationMatrix(i,0) = DdotXi / TMath::Sqrt(d2 * XiNorm);

    for (j = 0; j < i; j++) {
      for (k = 0; k < fSampleSize; k++) {
	// Index of sample j of variable i
	l =  j * fNVariables + k;
	XidotXj += (fVariables(i) - fMeanVariables(i))
	  * (fVariables(l) - fMeanVariables(j));
	XjNorm  += (fVariables(l) - fMeanVariables(j))
	  * (fVariables(l) - fMeanVariables(j));
      }
      fCorrelationMatrix(i+1,j) = XidotXj / TMath::Sqrt(XiNorm * XjNorm);
    }
  }
}



//____________________________________________________________________
 Double_t TMultiDimFit::MakeGramSchmidt(Int_t function)
{
  // PRIVATE METHOD:
  // Make Gram-Schmidt orthogonalisation. The class description gives
  // a thorough account of this algorithm, as well as
  // references. Please refer to the
  // class description


  // calculate w_i, that is, evaluate the current function at data
  // point i
  Double_t f2                        = 0;
  fOrthCoefficients(fNCoefficients)      = 0;
  fOrthFunctionNorms(fNCoefficients)  = 0;
  Int_t j        = 0;
  Int_t k        = 0;

  for (j = 0; j < fSampleSize; j++) {
    fFunctions(fNCoefficients, j) = 1;
    fOrthFunctions(fNCoefficients, j) = 0;
    // First, however, we need to calculate f_fNCoefficients
    for (k = 0; k < fNVariables; k++) {
      Int_t    p   =  fPowers[function * fNVariables + k];
      Double_t x   =  fVariables(j * fNVariables + k);
      fFunctions(fNCoefficients, j) *= EvalFactor(p,x);
    }

    // Calculate f dot f in f2
    f2 += fFunctions(fNCoefficients,j) *  fFunctions(fNCoefficients,j);
    // Assign to w_fNCoefficients f_fNCoefficients
    fOrthFunctions(fNCoefficients, j) = fFunctions(fNCoefficients, j);
  }

  // the first column of w is equal to f
  for (j = 0; j < fNCoefficients; j++) {
    Double_t fdw = 0;
    // Calculate (f_fNCoefficients dot w_j) / w_j^2
    for (k = 0; k < fSampleSize; k++) {
      fdw += fFunctions(fNCoefficients, k) * fOrthFunctions(j,k)
	/ fOrthFunctionNorms(j);
    }

    fOrthCurvatureMatrix(fNCoefficients,j) = fdw;
    // and subtract it from the current value of w_ij
    for (k = 0; k < fSampleSize; k++)
      fOrthFunctions(fNCoefficients,k) -= fdw * fOrthFunctions(j,k);
  }

  for (j = 0; j < fSampleSize; j++) {
    // calculate squared length of w_fNCoefficients
    fOrthFunctionNorms(fNCoefficients) +=
      fOrthFunctions(fNCoefficients,j)
      * fOrthFunctions(fNCoefficients,j);

    // calculate D dot w_fNCoefficients in A
    fOrthCoefficients(fNCoefficients) += fQuantity(j)
      * fOrthFunctions(fNCoefficients, j);
  }

  // First test, but only if didn't user specify
  if (!fIsUserFunction)
    if (TMath::Sqrt(fOrthFunctionNorms(fNCoefficients) / (f2 + 1e-10))
	< TMath::Sin(fMinAngle*DEGRAD))
      return 0;

  // The result found by this code for the first residual is always
  // much less then the one found be MUDIFI. That's because it's
  // supposed to be. The cause is the improved precision of Double_t
  // over DOUBLE PRECISION!
  fOrthCurvatureMatrix(fNCoefficients,fNCoefficients) = 1;
  Double_t b = fOrthCoefficients(fNCoefficients);
  fOrthCoefficients(fNCoefficients) /= fOrthFunctionNorms(fNCoefficients);

  // Calculate the residual from including this fNCoefficients.
  Double_t dResidur = fOrthCoefficients(fNCoefficients) * b;

  return dResidur;
}


//____________________________________________________________________
 void TMultiDimFit::MakeHistograms(Option_t *option)
{
  // Make histograms of the result of the analysis. This message
  // should be sent after having read all data points, but before
  // finding the parameterization
  //
  // Options:
  //     A         All the below
  //     X         Original independent variables
  //     D         Original dependent variables
  //     N         Normalised independent variables
  //     S         Shifted dependent variables
  //     R1        Residuals versus normalised independent variables
  //     R2        Residuals versus dependent variable
  //     R3        Residuals computed on training sample
  //     R4        Residuals computed on test sample
  //
  // For a description of these quantities, refer to
  // class description
  TString opt(option);
  opt.ToLower();

  if (opt.Length() < 1)
    return;

  if (!fHistograms)
    fHistograms = new TList;

  // Counter variable
  Int_t i = 0;

  // Histogram of original variables
  if (opt.Contains("x") || opt.Contains("a")) {
    SETBIT(fHistogramMask,HIST_XORIG);
    for (i = 0; i < fNVariables; i++)
      if (!fHistograms->FindObject(Form("x_%d_orig",i)))
	fHistograms->Add(new TH1D(Form("x_%d_orig",i),
				  Form("Original variable # %d",i),
				  100, fMinVariables(i),
				  fMaxVariables(i)));
  }

  // Histogram of original dependent variable
  if (opt.Contains("d") || opt.Contains("a")) {
    SETBIT(fHistogramMask,HIST_DORIG);
    if (!fHistograms->FindObject("d_orig"))
      fHistograms->Add(new TH1D("d_orig", "Original Quantity",
				100, fMinQuantity, fMaxQuantity));
  }

  // Histograms of normalized variables
  if (opt.Contains("n") || opt.Contains("a")) {
    SETBIT(fHistogramMask,HIST_XNORM);
    for (i = 0; i < fNVariables; i++)
      if (!fHistograms->FindObject(Form("x_%d_norm",i)))
	fHistograms->Add(new TH1D(Form("x_%d_norm",i),
				  Form("Normalized variable # %d",i),
				  100, -1,1));
  }

  // Histogram of shifted dependent variable
  if (opt.Contains("s") || opt.Contains("a")) {
    SETBIT(fHistogramMask,HIST_DSHIF);
    if (!fHistograms->FindObject("d_shifted"))
      fHistograms->Add(new TH1D("d_shifted", "Shifted Quantity",
				100, fMinQuantity - fMeanQuantity,
				fMaxQuantity - fMeanQuantity));
  }

  // Residual from training sample versus independent variables
  if (opt.Contains("r1") || opt.Contains("a")) {
    SETBIT(fHistogramMask,HIST_RX);
    for (i = 0; i < fNVariables; i++)
      if (!fHistograms->FindObject(Form("res_x_%d",i)))
	fHistograms->Add(new TH2D(Form("res_x_%d",i),
				  Form("Computed residual versus x_%d", i),
				  100, -1,    1,
				  35,
				  fMinQuantity - fMeanQuantity,
				  fMaxQuantity - fMeanQuantity));
  }

  // Residual from training sample versus. dependent variable
  if (opt.Contains("r2") || opt.Contains("a")) {
    SETBIT(fHistogramMask,HIST_RD);
    if (!fHistograms->FindObject("res_d"))
      fHistograms->Add(new TH2D("res_d",
				"Computed residuals vs Quantity",
				100,
				fMinQuantity - fMeanQuantity,
				fMaxQuantity - fMeanQuantity,
				35,
				fMinQuantity - fMeanQuantity,
				fMaxQuantity - fMeanQuantity));
  }

  // Residual from training sample
  if (opt.Contains("r3") || opt.Contains("a")) {
    SETBIT(fHistogramMask,HIST_RTRAI);
    if (!fHistograms->FindObject("res_train"))
      fHistograms->Add(new TH1D("res_train",
				"Computed residuals over training sample",
				100, fMinQuantity - fMeanQuantity,
				fMaxQuantity - fMeanQuantity));

  }
  if (opt.Contains("r4") || opt.Contains("a")) {
    SETBIT(fHistogramMask,HIST_RTEST);
    if (!fHistograms->FindObject("res_test"))
      fHistograms->Add(new TH1F("res_test",
				"Distribution of residuals from test",
				100,fMinQuantity - fMeanQuantity,
				fMaxQuantity - fMeanQuantity));
  }
}


//____________________________________________________________________
 void TMultiDimFit::MakeMethod(const Char_t* classname, Option_t* option)
{
  // Generate the file <classname>MDF.cxx which contains the
  // implementation of the method:
  //
  //   Double_t <classname>::MDF(Double_t *x)
  //
  // which does the same as  TMultiDimFit::Eval. Please refer to this
  // method.
  //
  // Further, the public static members:
  //
  //   Int_t    <classname>::fgNVariables
  //   Int_t    <classname>::fgNCoefficients
  //   Double_t <classname>::fgDMean
  //   Double_t <classname>::fgXMean[]       //[fgNVariables]
  //   Double_t <classname>::fgXMin[]        //[fgNVariables]
  //   Double_t <classname>::fgXMax[]        //[fgNVariables]
  //   Double_t <classname>::fgCoefficient[] //[fgNCoeffficents]
  //   Int_t    <classname>::fgPower[]       //[fgNCoeffficents*fgNVariables]
  //
  // are initialized, and assumed to exist. The class declaration is
  // assumed to be in <classname>.h and assumed to be provided by the
  // user.
  //
  // See TMultiDimFit::MakeRealCode for a list of options
  //
  // The minimal class definition is:
  //
  //   class <classname> {
  //   public:
  //     Int_t    <classname>::fgNVariables;     // Number of variables
  //     Int_t    <classname>::fgNCoefficients;  // Number of terms
  //     Double_t <classname>::fgDMean;          // Mean from training sample
  //     Double_t <classname>::fgXMean[];        // Mean from training sample
  //     Double_t <classname>::fgXMin[];         // Min from training sample
  //     Double_t <classname>::fgXMax[];         // Max from training sample
  //     Double_t <classname>::fgCoefficient[];  // Coefficients
  //     Int_t    <classname>::fgPower[];        // Function powers
  //
  //     Double_t Eval(Double_t *x);
  //   };
  //
  // Whether the method <classname>::Eval should be static or not, is
  // up to the user.

  MakeRealCode(Form("%sMDF.cxx", classname), classname, option);
}



//____________________________________________________________________
 void TMultiDimFit::MakeNormalized()
{
  // PRIVATE METHOD:
  // Normalize data to the interval [-1;1]. This is needed for the
  // classes method to work.

  Int_t i = 0;
  Int_t j = 0;
  Int_t k = 0;

  for (i = 0; i < fSampleSize; i++) {
    if (TESTBIT(fHistogramMask,HIST_DORIG))
      ((TH1F*)fHistograms->FindObject("d_orig"))->Fill(fQuantity(i));

    fQuantity(i) -= fMeanQuantity;
    fSumSqAvgQuantity  += fQuantity(i) * fQuantity(i);

    if (TESTBIT(fHistogramMask,HIST_DSHIF))
      ((TH1F*)fHistograms->FindObject("d_shifted"))->Fill(fQuantity(i));

    for (j = 0; j < fNVariables; j++) {
      Double_t range = 1. / (fMaxVariables(j) - fMinVariables(j));
      k              = i * fNVariables + j;

      // Fill histograms of original independent variables
      if (TESTBIT(fHistogramMask,HIST_XORIG))
	((TH1F*)fHistograms->FindObject(Form("x_%d_orig",j)))
	  ->Fill(fVariables(k));

      // Normalise independent variables
      fVariables(k) = 1 + 2 * range * (fVariables(k) - fMaxVariables(j));

      // Fill histograms of normalised independent variables
      if (TESTBIT(fHistogramMask,HIST_XNORM))
	((TH1F*)fHistograms->FindObject(Form("x_%d_norm",j)))
	  ->Fill(fVariables(k));

    }
  }
  // Shift min and max of dependent variable
  fMaxQuantity -= fMeanQuantity;
  fMinQuantity -= fMeanQuantity;

  // Shift mean of independent variables
  for (i = 0; i < fNVariables; i++) {
    Double_t range = 1. / (fMaxVariables(i) - fMinVariables(i));
    fMeanVariables(i) = 1 + 2 * range * (fMeanVariables(i)
					 - fMaxVariables(i));
  }
}


//____________________________________________________________________
 void TMultiDimFit::MakeParameterization()
{
  // PRIVATE METHOD:
  // Find the parameterization over the training sample. A full account
  // of the algorithm is given in the
  // class description

  Int_t     i              = -1;
  Int_t     j              = 0;
  Int_t     k              = 0;
  Int_t     maxPass        = 3;
  Int_t     studied        = 0;
  Double_t  squareResidual = fSumSqAvgQuantity;
  fNCoefficients            = 0;
  fSumSqResidual           = fSumSqAvgQuantity;
  fFunctions.ResizeTo(fMaxTerms,fSampleSize);
  fOrthFunctions.ResizeTo(fMaxTerms,fSampleSize);
  fOrthFunctionNorms.ResizeTo(fMaxTerms);
  fOrthCoefficients.ResizeTo(fMaxTerms);
  fOrthCurvatureMatrix.ResizeTo(fMaxTerms,fMaxTerms);
  fFunctions = 1;

  fFunctionCodes = new Int_t[fMaxFunctions];
  fPowerIndex    = new Int_t[fMaxTerms];

  if (fMaxAngle != 0)  maxPass = 100;
  if (fIsUserFunction) maxPass = 1;

  // Loop over the number of functions we want to study.
  // increment inspection counter
  while(kTRUE) {

    // Reach user defined limit of studies
    if (studied++ >= fMaxStudy) {
      fParameterisationCode = PARAM_MAXSTUDY;
      break;
    }

    // Considered all functions several times
    if (k >= maxPass) {
      fParameterisationCode = PARAM_SEVERAL;
      break;
    }

    // increment function counter
    i++;

    // If we've reached the end of the functions, restart pass
    if (i == fMaxFunctions) {
      if (fMaxAngle != 0)
	fMaxAngle += (90 - fMaxAngle) / 2;
      i = 0;
      studied--;
      k++;
      continue;
    }

    if (studied == 1)
      fFunctionCodes[i] = 0;
    else if (fFunctionCodes[i] >= 2)
      continue;

    // Print a happy message
    if (fIsVerbose && studied == 1)
      cout << "Coeff   SumSqRes    Contrib   Angle      QM   Func"
	   << "     Value        W^2  Powers" << endl;

    // Make the Gram-Schmidt
    Double_t dResidur = MakeGramSchmidt(i);

    if (dResidur == 0) {
      // This function is no good!
      // First test is in MakeGramSchmidt
      fFunctionCodes[i] = 1;
      continue;
    }

    // If user specified function, assume she/he knows what he's doing
    if (!fIsUserFunction) {
      // Flag this function as considered
      fFunctionCodes[i] = 2;

      // Test if this function contributes to the fit
      if (!TestFunction(squareResidual, dResidur)) {
	fFunctionCodes[i] = 1;
	continue;
      }
    }

    // If we get to here, the function currently considered is
    // fNCoefficients, so we increment the counter
    // Flag this function as OK, and store and the number in the
    // index.
    fFunctionCodes[i]          = 3;
    fPowerIndex[fNCoefficients] = i;
    fNCoefficients++;

    // We add the current contribution to the sum of square of
    // residuals;
    squareResidual -= dResidur;


    // Calculate control parameter from this function
    for (j = 0; j < fNVariables; j++) {
      if (fNCoefficients == 1
	  || fMaxPowersFinal[j] <= fPowers[i * fNVariables + j] - 1)
	fMaxPowersFinal[j] = fPowers[i * fNVariables + j] - 1;
    }
    Double_t s = EvalControl(&fPowers[i * fNVariables]);

    // Print the statistics about this function
    if (fIsVerbose) {
#ifndef R__MACOSX
      cout << setw(5)  << fNCoefficients << " "
	   << setw(10) << setprecision(4) << squareResidual << " "
	   << setw(10) << setprecision(4) << dResidur << " "
	   << setw(7)  << setprecision(3) << fMaxAngle << " "
	   << setw(7)  << setprecision(3) << s << " "
	   << setw(5)  << i << " "
	   << setw(10) << setprecision(4)
	   << fOrthCoefficients(fNCoefficients-1) << " "
	   << setw(10) << setprecision(4)
	   << fOrthFunctionNorms(fNCoefficients-1) << " "
	   << flush;
#else
      printf("%d5 %g10.4 %g10.4 %g7.3 %g7.3 %d5 %g10.4 %g10.4 ",
             fNCoefficients, squareResidual, dResidur, fMaxAngle, s, i,
             fOrthCoefficients(fNCoefficients-1),
             fOrthFunctionNorms(fNCoefficients-1));
      fflush(stdout);
#endif
      for (j = 0; j < fNVariables; j++)
	cout << " " << fPowers[i * fNVariables + j] - 1 << flush;
      cout << endl;
    }

    if (fNCoefficients >= fMaxTerms /* && fIsVerbose */) {
      fParameterisationCode = PARAM_MAXTERMS;
      break;
    }

    Double_t err  = TMath::Sqrt(TMath::Max(1e-20,squareResidual) /
				fSumSqAvgQuantity);
    if (err < fMinRelativeError) {
      fParameterisationCode = PARAM_RELERR;
      break;
    }

  }

  fError          = TMath::Max(1e-20,squareResidual);
  fSumSqResidual -= fError;
  fRMS = TMath::Sqrt(fError / fSampleSize);
}


//____________________________________________________________________
 void TMultiDimFit::MakeRealCode(const char *filename,
				const char *classname,
				Option_t *)
{
  // PRIVATE METHOD:
  // This is the method that actually generates the code for the
  // evaluation the parameterization on some point.
  // It's called by TMultiDimFit::MakeCode and TMultiDimFit::MakeMethod.
  //
  // The options are: NONE so far
  Int_t i, j;

  Bool_t  isMethod     = (classname[0] == '0' ? kFALSE : kTRUE);
  const char *prefix   = (isMethod ? Form("%s::", classname) : "");
  const char *cv_qual  = (isMethod ? "" : "static ");

  ofstream outFile(filename,ios::out
  if (!outFile) {
    Error("MakeRealCode","couldn't open output file '%s'",filename);
    return;
  }

  if (fIsVerbose)
    cout << "Writing on file "" << filename << "" ... " << flush;
  //
  // Write header of file
  //
  // Emacs mode line ;-)
  outFile << "// -*- mode: c++ -*-" << endl;
  // Info about creator
  outFile << "// " << endl
          << "// File " << filename
          << " generated by TMultiDimFit::MakeRealCode" << endl;
  // Time stamp
  TDatime date;
  outFile << "// on " << date.AsString() << endl;
  // ROOT version info
  outFile << "// ROOT version " << gROOT->GetVersion()
          << endl << "//" << endl;
  // General information on the code
  outFile << "// This file contains the function " << endl
          << "//" << endl
          << "//    Double_t  " << prefix << "MDF(Double_t *x); " << endl
          << "//" << endl
          << "// For evaluating the parameterization obtained" << endl
	  << "// from TMultiDimFit and the point x" << endl
          << "// " << endl
          << "// See TMultiDimFit class documentation for more "
          << "information " << endl << "// " << endl;
  // Header files
  outFile << "#ifndef __CINT__" << endl;
  if (isMethod)
    // If these are methods, we need the class header
    outFile << "#include "" << classname << ".h"" << endl;
  else
    // otherwise, we need the typedefs of Int_t and Double_t
    outFile << "#include <Rtypes.h> // needed for Double_t etc" << endl;
  // Finish the preprocessor block
  outFile << "#endif" << endl << endl;

  //
  // Now for the data
  //
  outFile << "//" << endl
          << "// Static data variables"  << endl
          << "//" << endl;
  outFile << cv_qual << "Int_t    " << prefix << "gNVariables    = "
          << fNVariables << ";" << endl;
  outFile << cv_qual << "Int_t    " << prefix << "gNCoefficients = "
          << fNCoefficients << ";" << endl;
  outFile << cv_qual << "Double_t " << prefix << "gDMean         = "
          << fMeanQuantity << ";" << endl;

  // Assignment to mean vector.
  outFile << "// Assignment to mean vector." << endl;
  outFile << cv_qual << "Double_t " << prefix
          << "gXMean[] = {" << endl;
  for (i = 0; i < fNVariables; i++)
    outFile << (i != 0 ? ", " : "  ") << fMeanVariables(i) << flush;
  outFile << " };" << endl << endl;

  // Assignment to minimum vector.
  outFile << "// Assignment to minimum vector." << endl;
  outFile << cv_qual << "Double_t " << prefix
          << "gXMin[] = {" << endl;
  for (i = 0; i < fNVariables; i++)
    outFile << (i != 0 ? ", " : "  ") << fMinVariables(i) << flush;
  outFile << " };" << endl << endl;

  // Assignment to maximum vector.
  outFile << "// Assignment to maximum vector." << endl;
  outFile << cv_qual << "Double_t " << prefix
          << "gXMax[] = {" << endl;
  for (i = 0; i < fNVariables; i++)
    outFile << (i != 0 ? ", " : "  ") << fMaxVariables(i) << flush;
  outFile << " };" << endl << endl;

  // Assignment to coefficients vector.
  outFile << "// Assignment to coefficients vector." << endl;
  outFile << cv_qual << "Double_t " << prefix
          << "gCoefficient[] = {" << flush;
  for (i = 0; i < fNCoefficients; i++)
    outFile << (i != 0 ? "," : "") << endl
            << "  " << fCoefficients(i) << flush;
  outFile << endl << " };" << endl << endl;

  // Assignment to powers vector.
  outFile << "// Assignment to powers vector." << endl
	  << "// The powers are stored row-wise, that is" << endl
	  << "//  p_ij = " << prefix
	  << "gPower[i * NVariables + j];" << endl;
  outFile << cv_qual << "Int_t    " << prefix
          << "gPower[] = {" << flush;
  for (i = 0; i < fNCoefficients; i++) {
    for (j = 0; j < fNVariables; j++) {
      if (j != 0) outFile << flush << "  ";
      else        outFile << endl << "  ";
      outFile << fPowers[fPowerIndex[i] * fNVariables + j]
	      << (i == fNCoefficients - 1 && j == fNVariables - 1 ? "" : ",")
	      << flush;
    }
  }
  outFile << endl << "};" << endl << endl;


  //
  // Finally we reach the function itself
  //
  outFile << "// " << endl
          << "// The "
          << (isMethod ? "method " : "function ")
          << "  Double_t " << prefix
          << "MDF(Double_t *x)"
          << endl << "// " << endl;
  outFile << "Double_t " << prefix
          << "MDF(Double_t *x) {" << endl
	  << "  Double_t returnValue = " << prefix << "gDMean;" << endl
	  << "  Int_t    i = 0, j = 0;" << endl
	  << "  for (i = 0; i < " << prefix << "gNCoefficients ; i++) {"
	  << endl
	  << "    // Evaluate the ith term in the expansion" << endl
	  << "    Double_t term = " << prefix << "gCoefficient[i]; <<"
	  << endl
	  << "    for (j = 0; j < " << prefix << "gNVariables; j++) {"
	  << endl
	  << "      // Evaluate the polynomial in the jth variable." << endl
	  << "      Int_t power = "<< prefix << "gPower["
	  << prefix << "gNVariables * i + j]; " << endl
	  << "      Double_t p1 = 1, p2 = 0, p3 = 0, r = 0;" << endl
	  << "      Double_t v =  1 + 2. / ("
	  << prefix << "gXMax[j] - " << prefix
	  << "gXMin[j]) * (x[j] - " << prefix << "gXMax[j]);" << endl
	  << "      // what is the power to use!" << endl
	  << "      switch(power) {" << endl
	  << "      case 1: r = 1; break; " << endl
	  << "      case 2: r = v; break; " << endl
	  << "      default: " << endl
	  << "        p2 = v; " << endl
	  << "        for (i = 3; i <= power; i++) { " << endl
	  << "          p3 = p2 * v;" << endl;
  if (fPolyType == kLegendre)
    outFile << "          p3 = ((2 * i - 3) * p2 * v - (i - 2) * p1)"
	    << " / (i - 1);" << endl;
  if (fPolyType == kChebyshev)
    outFile << "          p3 = 2 * v * p2 - p1; " << endl;
  outFile << "          p1 = p2; p2 = p3; " << endl << "        }" << endl
	  << "        r = p3;" << endl << "      }" << endl
	  << "      // multiply this term by the poly in the jth var" << endl
	  << "      term *= r; " << endl << "    }" << endl
	  << "    // Add this term to the final result" << endl
	  << "    returnValue += term;" << endl << "  }" << endl
	  << "  return returnValue;" << endl << "}" << endl << endl;

  // EOF
  outFile << "// EOF for " << filename << endl;

  // Close the file
  outFile.close();

  if (fIsVerbose)
    cout << "done" << endl;
}


//____________________________________________________________________
 void TMultiDimFit::Print(Option_t *option) const
{
  // Print statistics etc.
  // Options are
  //   P        Parameters
  //   S        Statistics
  //   C        Coefficients
  //   R        Result of parameterisation
  //   F        Result of fit
  //   K        Correlation Matrix
  //   M        Pretty print formula
  //
  Int_t i = 0;
  Int_t j = 0;

  TString opt(option);
  opt.ToLower();

  if (opt.Contains("p")) {
    // Print basic parameters for this object
    cout << "User parameters:" << endl
	 << "----------------" << endl
	 << " Variables:                    " << fNVariables << endl
	 << " Data points:                  " << fSampleSize << endl
	 << " Max Terms:                    " << fMaxTerms << endl
	 << " Power Limit Parameter:        " << fPowerLimit << endl
	 << " Max functions:                " << fMaxFunctions << endl
	 << " Max functions to study:       " << fMaxStudy << endl
	 << " Max angle (optional):         " << fMaxAngle << endl
	 << " Min angle:                    " << fMinAngle << endl
	 << " Relative Error accepted:      " << fMinRelativeError << endl
	 << " Maximum Powers:               " << flush;
    for (i = 0; i < fNVariables; i++)
      cout << " " << fMaxPowers[i] - 1 << flush;
    cout << endl << endl
	 << " Parameterisation will be done using " << flush;
    if (fPolyType == kChebyshev)
      cout << "Chebyshev polynomials" << endl;
    else if (fPolyType == kLegendre)
      cout << "Legendre polynomials" << endl;
    else
      cout << "Monomials" << endl;
    cout << endl;
  }

  if (opt.Contains("s")) {
    // Print statistics for read data
    cout << "Sample statistics:" << endl
	 << "------------------" << endl
	 << "                 D"  << flush;
#ifndef R__MACOSX
    for (i = 0; i < fNVariables; i++)
      cout << " " << setw(10) << i+1 << flush;
    cout << endl << " Max:   " << setw(10) << setprecision(7)
	 << fMaxQuantity << flush;
    for (i = 0; i < fNVariables; i++)
      cout << " " << setw(10) << setprecision(4)
	   << fMaxVariables(i) << flush;
    cout << endl << " Min:   " << setw(10) << setprecision(7)
	 << fMinQuantity << flush;
    for (i = 0; i < fNVariables; i++)
      cout << " " << setw(10) << setprecision(4)
	   << fMinVariables(i) << flush;
    cout << endl << " Mean:  " << setw(10) << setprecision(7)
	 << fMeanQuantity << flush;
    for (i = 0; i < fNVariables; i++)
      cout << " " << setw(10) << setprecision(4)
	   << fMeanVariables(i) << flush;
    cout << endl << " Function Sum Squares:         " << fSumSqQuantity
	 << endl << endl;
#else
    for (i = 0; i < fNVariables; i++) {
      fprintf(stdout," %d10",i+1);
      fflush(stdout);
    }
    fprintf(stdout,"\n Max:   %g10.7",fMaxQuantity);
    fflush(stdout);
    for (i = 0; i < fNVariables; i++) {
      fprintf(stdout," %g10.4",fMaxVariables(i));
      fflush(stdout);
    }
    fprintf(stdout,"\n Min:   %g10.7",fMinQuantity);
    fflush(stdout);
    for (i = 0; i < fNVariables; i++) {
      fprintf(stdout," %g10.4",fMinVariables(i));
      fflush(stdout);
    }
    fprintf(stdout,"\n Mean:  %g10.7",fMeanQuantity);
    fflush(stdout);
    for (i = 0; i < fNVariables; i++) {
      fprintf(stdout," %g10.4",fMeanVariables(i));
      fflush(stdout);
    }
    fprintf(stdout,"\n Function Sum Squares:         %g10.4\n\n",
	    fSumSqQuantity);
#endif
  }

  if (opt.Contains("r")) {
    cout << "Results of Parameterisation:" << endl
	 << "----------------------------" << endl
	 << " Total reduction of square residuals    "
	 << fSumSqResidual << endl
	 << " Relative precision obtained:           "
	 << fPrecision   << endl
	 << " Error obtained:                        "
	 << fError << endl
	 << " Multiple correlation coefficient:      "
	 << fCorrelationCoeff   << endl
	 << " Reduced Chi square over sample:        "
	 << fChi2 / (fSampleSize - fNCoefficients) << endl
	 << " Maximum residual value:                "
	 << fMaxResidual << endl
	 << " Minimum residual value:                "
	 << fMinResidual << endl
	 << " Estimated root mean square:            "
	 << fRMS << endl
	 << " Maximum powers used:                   " << flush;
    for (j = 0; j < fNVariables; j++)
      cout << fMaxPowersFinal[j] << " " << flush;
    cout << endl
	 << " Function codes of candidate functions." << endl
	 << "  1: considered,"
	 << "  2: too little contribution,"
	 << "  3: accepted." << flush;
    for (i = 0; i < fMaxFunctions; i++) {
      if (i % 60 == 0)
	cout << endl << " " << flush;
      else if (i % 10 == 0)
	cout << " " << flush;
      cout << fFunctionCodes[i];
    }
    cout << endl << " Loop over candidates stopped because " << flush;
    switch(fParameterisationCode){
    case PARAM_MAXSTUDY:
      cout << "max allowed studies reached" << endl; break;
    case PARAM_SEVERAL:
      cout << "all candidates considered several times" << endl; break;
    case PARAM_RELERR:
      cout << "wanted relative error obtained" << endl; break;
    case PARAM_MAXTERMS:
      cout << "max number of terms reached" << endl; break;
    default:
      cout << "some unknown reason" << endl;
      break;
    }
    cout << endl;
  }

  if (opt.Contains("f")) {
    cout << "Results of Fit:" << endl
	 << "---------------" << endl
	 << " Test sample size:                      "
	 << fTestSampleSize << endl
	 << " Multiple correlation coefficient:      "
	 << fTestCorrelationCoeff << endl
	 << " Relative precision obtained:           "
	 << fTestPrecision   << endl
	 << " Error obtained:                        "
	 << fTestError << endl
	 << " Reduced Chi square over sample:        "
	 << fChi2 / (fSampleSize - fNCoefficients) << endl
	 << endl;
    if (fFitter) {
      fFitter->PrintResults(1,1);
      cout << endl;
    }
  }

  if (opt.Contains("c")){
    cout << "Coefficients:" << endl
	 << "-------------" << endl
	 << "   #         Value        Error   Powers" << endl
	 << " ---------------------------------------" << endl;
    for (i = 0; i < fNCoefficients; i++) {
#ifndef R__MACOSX
      cout << " " << setw(3) << i << "  "
	   << setw(12) << fCoefficients(i) << "  "
	   << setw(12) << fCoefficientsRMS(i) << "  " << flush;
      for (j = 0; j < fNVariables; j++)
	cout << " " << setw(3)
	     << fPowers[fPowerIndex[i] * fNVariables + j] - 1 << flush;
#else
      fprintf(stdout," %d3 %g12.4 %g12.4 ",
              i,fCoefficients(i),fCoefficientsRMS(i));
      fflush(stdout);
      for (j = 0; j < fNVariables; j++) {
        fprintf(stdout," %d3",fPowers[fPowerIndex[i] * fNVariables + j] - 1);
        fflush(stdout);
      }
#endif
      cout << endl;
    }
    cout << endl;
  }
  if (opt.Contains("k") && fCorrelationMatrix.IsValid()) {
    cout << "Correlation Matrix:" << endl
	 << "-------------------";
    fCorrelationMatrix.Print();
  }
  
  if (opt.Contains("m")) {
    cout << "Parameterization:" << endl
	 << "-----------------" << endl
	 << "  Normalised variables: " << endl;
    for (i = 0; i < fNVariables; i++) 
      cout << "\ty_" << i << "\t= 1 + 2 * (x_" << i << " - " 
	   << fMaxVariables(i) << ") / (" 
	   << fMaxVariables(i) << " - " << fMinVariables(i) << ")" 
	   << endl;
    cout << endl
	 << "  f(";
    for (i = 0; i < fNVariables; i++) {
      cout << "y_" << i;
      if (i != fNVariables-1) cout << ", ";
    }
    cout << ") = ";
    for (Int_t i = 0; i < fNCoefficients; i++) {
      if (i != 0)
	cout << endl << "\t" << (fCoefficients(i) < 0 ? "- " : "+ ") 
	     << TMath::Abs(fCoefficients(i));
      else 
	cout << fCoefficients(i);
      for (Int_t j = 0; j < fNVariables; j++) {
	Int_t p = fPowers[fPowerIndex[i] * fNVariables + j];
	switch (p) { 
	case 1: break;
	case 2: cout << " * y_" << j; break;
	default:
	  switch(fPolyType) {
	  case kLegendre:  cout << " * L_" << p << "(y_" << j << ")"; break;
	  case kChebyshev: cout << " * C_" << p << "(y_" << j << ")"; break;
	  default:         cout << " * y_" << j << "^" << p; break;
	  }
	}
			  
      }
    }
    cout << endl;
  }
}


//____________________________________________________________________
 Bool_t TMultiDimFit::Select(const Int_t *)
{
  // Selection method. User can override this method for specialized
  // selection of acceptable functions in fit. Default is to select
  // all. This message is sent during the build-up of the function
  // candidates table once for each set of powers in
  // variables. Notice, that the argument array contains the powers
  // PLUS ONE. For example, to De select the function
  //     f = x1^2 * x2^4 * x3^5,
  // this method should return kFALSE if given the argument
  //     { 3, 4, 6 }
  return kTRUE;
}

//____________________________________________________________________
 void TMultiDimFit::SetMaxAngle(Double_t ang)
{
  // Set the max angle (in degrees) between the initial data vector to
  // be fitted, and the new candidate function to be included in the
  // fit.  By default it is 0, which automatically chooses another
  // selection criteria. See also
  // class description
  if (ang >= 90 || ang < 0) {
    Warning("SetMaxAngle", "angle must be in [0,90)");
    return;
  }

  fMaxAngle = ang;
}

//____________________________________________________________________
 void TMultiDimFit::SetMinAngle(Double_t ang)
{
  // Set the min angle (in degrees) between a new candidate function
  // and the subspace spanned by the previously accepted
  // functions. See also
  // class description
  if (ang > 90 || ang <= 0) {
    Warning("SetMinAngle", "angle must be in [0,90)");
    return;
  }

  fMinAngle = ang;

}


//____________________________________________________________________
 void TMultiDimFit::SetPowers(const Int_t* powers, Int_t terms)
{
  // Define a user function. The input array must be of the form
  // (p11, ..., p1N, ... ,pL1, ..., pLN)
  // Where N is the dimension of the data sample, L is the number of
  // terms (given in terms) and the first number, labels the term, the
  // second the variable.  More information is given in the
  // class description
  fIsUserFunction = kTRUE;
  fMaxFunctions   = terms;
  fMaxTerms       = terms;
  fMaxStudy       = terms;
  fPowers         = new Int_t[fMaxFunctions * fNVariables];
  Int_t i, j;
  for (i = 0; i < fMaxFunctions; i++)
    for(j = 0; j < fNVariables; j++)
      fPowers[i * fNVariables + j] = powers[i * fNVariables + j]  + 1;
}

//____________________________________________________________________
 void TMultiDimFit::SetPowerLimit(Double_t limit)
{
  // Set the user parameter for the function selection. The bigger the
  // limit, the more functions are used. The meaning of this variable
  // is defined in the
  // class description
  fPowerLimit = limit;
}

//____________________________________________________________________
 void TMultiDimFit::SetMaxPowers(const Int_t* powers)
{
  // Set the maximum power to be considered in the fit for each
  // variable. See also
  // class description
  if (!powers)
    return;

  for (Int_t i = 0; i < fNVariables; i++)
    fMaxPowers[i] = powers[i]+1;
}

//____________________________________________________________________
 void TMultiDimFit::SetMinRelativeError(Double_t error)
{
  // Set the acceptable relative error for when sum of square
  // residuals is considered minimized. For a full account, refer to
  // the
  // class description
  fMinRelativeError = error;
}


//____________________________________________________________________
 Bool_t TMultiDimFit::TestFunction(Double_t squareResidual,
				  Double_t dResidur)
{
  // PRIVATE METHOD:
  // Test whether the currently considered function contributes to the
  // fit. See also
  // class description

  if (fNCoefficients != 0) {
    // Now for the second test:
    if (fMaxAngle == 0) {
      // If the user hasn't supplied a max angle do the test as,
      if (dResidur <
	  squareResidual / (fMaxTerms - fNCoefficients + 1 + 1E-10)) {
	return kFALSE;
      }
    }
    else {
      // If the user has provided a max angle, test if the calculated
      // angle is less then the max angle.
      if (TMath::Sqrt(dResidur/fSumSqAvgQuantity) <
	  TMath::Cos(fMaxAngle*DEGRAD)) {
	return kFALSE;
      }
    }
  }
  // If we get here, the function is OK
  return kTRUE;
}


//____________________________________________________________________
void mdfHelper(int& /*npar*/, double* /*divs*/, double& chi2,
	       double* coeffs, int /*flag*/)
{
  // Helper function for doing the minimisation of Chi2 using Minuit

  // Get pointer  to current TMultiDimFit object.
  TMultiDimFit* mdf = TMultiDimFit::Instance();
  chi2     = mdf->MakeChi2(coeffs);
}

ROOT page - Class index - Top of the page
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.