ROOT  6.06/09
Reference Guide
TMultiDimFit.cxx
Go to the documentation of this file.
1 // @(#)root/hist:$Id$
2 // Author: Christian Holm Christensen 07/11/2000
3 
4 /** \class TMultiDimFit
5  \ingroup Hist
6 
7  Multidimensional Fits in ROOT.
8  ## Overview
9  A common problem encountered in different fields of applied science is
10  to find an expression for one physical quantity in terms of several
11  others, which are directly measurable.
12 
13  An example in high energy physics is the evaluation of the momentum of
14  a charged particle from the observation of its trajectory in a magnetic
15  field. The problem is to relate the momentum of the particle to the
16  observations, which may consists of positional measurements at
17  intervals along the particle trajectory.
18 
19  The exact functional relationship between the measured quantities
20  (e.g., the space-points) and the dependent quantity (e.g., the
21  momentum) is in general not known, but one possible way of solving the
22  problem, is to find an expression which reliably approximates the
23  dependence of the momentum on the observations.
24 
25  This explicit function of the observations can be obtained by a
26  <I>least squares</I> fitting procedure applied to a representive
27  sample of the data, for which the dependent quantity (e.g., momentum)
28  and the independent observations are known. The function can then be
29  used to compute the quantity of interest for new observations of the
30  independent variables.
31 
32  This class <TT>TMultiDimFit</TT> implements such a procedure in
33  ROOT. It is largely based on the CERNLIB MUDIFI package [2].
34  Though the basic concepts are still sound, and
35  therefore kept, a few implementation details have changed, and this
36  class can take advantage of MINUIT [4] to improve the errors
37  of the fitting, thanks to the class TMinuit.
38 
39  In [5] and [6] H. Wind demonstrates the utility
40  of this procedure in the context of tracking, magnetic field
41  parameterisation, and so on. The outline of the method used in this
42  class is based on Winds discussion, and I refer these two excellents
43  text for more information.
44 
45  And example of usage is given in $ROOTSYS/tutorials/fit/multidimfit.C.
46 
47  ## The Method
48  Let \f$ D \f$ by the dependent quantity of interest, which depends smoothly
49  on the observable quantities \f$ x_1, \ldots, x_N \f$ which we'll denote by
50  \f$\mathbf{x}\f$. Given a training sample of \f$ M\f$ tuples of the form, (TMultiDimFit::AddRow)
51 
52  \f[
53  \left(\mathbf{x}_j, D_j, E_j\right)\quad,
54  \f]
55  where \f$\mathbf{x}_j = (x_{1,j},\ldots,x_{N,j})\f$ are \f$ N\f$ independent
56  variables, \f$ D_j\f$ is the known, quantity dependent at \f$\mathbf{x}_j\f$ and \f$ E_j\f$ is
57  the square error in \f$ D_j\f$, the class will try to find the parameterization
58  \f[
59  D_p(\mathbf{x}) = \sum_{l=1}^{L} c_l \prod_{i=1}^{N} p_{li}\left(x_i\right)
60  = \sum_{l=1}^{L} c_l F_l(\mathbf{x})
61  \f]
62  such that
63 
64  \f[
65  S \equiv \sum_{j=1}^{M} \left(D_j - D_p\left(\mathbf{x}_j\right)\right)^2
66  \f]
67  is minimal. Here \f$p_{li}(x_i)\f$ are monomials, or Chebyshev or Legendre
68  polynomials, labelled \f$l = 1, \ldots, L\f$, in each variable \f$ x_i\f$,\f$ i=1, \ldots, N\f$.
69 
70  So what TMultiDimFit does, is to determine the number of terms \f$ L\f$, and then \f$ L\f$ terms
71  (or functions) \f$ F_l\f$, and the \f$ L\f$ coefficients \f$ c_l\f$, so that \f$ S\f$ is minimal
72  (TMultiDimFit::FindParameterization).
73 
74  Of course it's more than a little unlikely that \f$ S\f$ will ever become
75  exact zero as a result of the procedure outlined below. Therefore, the
76  user is asked to provide a minimum relative error \f$ \epsilon\f$ (TMultiDimFit::SetMinRelativeError),
77  and \f$ S\f$ will be considered minimized when
78 
79  \f[
80  R = \frac{S}{\sum_{j=1}^M D_j^2} < \epsilon
81  \f]
82  Optionally, the user may impose a functional expression by specifying
83  the powers of each variable in \f$ L\f$ specified functions \f$ F_1, \ldots,F_L\f$ (TMultiDimFit::SetPowers).
84  In that case, only the coefficients \f$ c_l\f$ is calculated by the class.
85 
86  ## Limiting the Number of Terms
87  As always when dealing with fits, there's a real chance of *over fitting*. As is well-known, it's
88  always possible to fit an \f$ N-1\f$ polynomial in \f$ x\f$ to \f$ N\f$ points \f$ (x,y)\f$ with
89  \f$\chi^2 = 0\f$, but the polynomial is not likely to fit new data at all [1].
90  Therefore, the user is asked to provide an upper limit, \f$ L_{max}\f$ to the number of terms in
91  \f$ D_p\f$ (TMultiDimFit::SetMaxTerms).
92 
93  However, since there's an infinite number of \f$ F_l\f$ to choose from, the
94  user is asked to give the maximum power. \f$ P_{max,i}\f$, of each variable
95  \f$ x_i\f$ to be considered in the minimization of \f$ S\f$ (TMultiDimFit::SetMaxPowers).
96 
97  One way of obtaining values for the maximum power in variable \f$ i\f$, is
98  to perform a regular fit to the dependent quantity \f$ D\f$, using a
99  polynomial only in \f$ x_i\f$. The maximum power is \f$ P_{max,i}\f$ is then the
100  power that does not significantly improve the one-dimensional
101  least-square fit over \f$ x_i\f$ to \f$ D\f$ [5].
102 
103  There are still a huge amount of possible choices for \f$ F_l\f$; in fact
104  there are \f$\prod_{i=1}^{N} (P_{max,i} + 1)\f$ possible
105  choices. Obviously we need to limit this. To this end, the user is
106  asked to set a *power control limit*, \f$ Q\f$ (TMultiDimFit::SetPowerLimit), and a function
107  \f$ F_l\f$ is only accepted if
108  \f[
109  Q_l = \sum_{i=1}^{N} \frac{P_{li}}{P_{max,i}} < Q
110  \f]
111  where \f$ P_{li}\f$ is the leading power of variable \f$ x_i\f$ in function \f$ F_l\f$ (TMultiDimFit::MakeCandidates).
112  So the number of functions increase with \f$ Q\f$ (1, 2 is fine, 5 is way out).
113 
114  ## Gram-Schmidt Orthogonalisation</A>
115  To further reduce the number of functions in the final expression,
116  only those functions that significantly reduce \f$ S\f$ is chosen. What
117  `significant' means, is chosen by the user, and will be
118  discussed below (see [2.3](TMultiFimFit.html#sec:selectiondetail)).
119 
120  The functions \f$ F_l\f$ are generally not orthogonal, which means one will
121  have to evaluate all possible \f$ F_l\f$'s over all data-points before
122  finding the most significant [1]. We can, however, do
123  better then that. By applying the *modified Gram-Schmidt
124  orthogonalisation* algorithm [5] [3] to the
125  functions \f$ F_l\f$, we can evaluate the contribution to the reduction of
126  \f$ S\f$ from each function in turn, and we may delay the actual inversion
127  of the curvature-matrix (TMultiDimFit::MakeGramSchmidt).
128 
129  So we are let to consider an \f$ M\times L\f$ matrix \f$\mathsf{F}\f$, an
130  element of which is given by
131  \f[
132  f_{jl} = F_j\left(x_{1j} , x_{2j}, \ldots, x_{Nj}\right)
133  = F_l(\mathbf{x}_j)\, \quad\mbox{with}~j=1,2,\ldots,M,
134  \f]
135  where \f$ j\f$ labels the \f$ M\f$ rows in the training sample and \f$ l\f$ labels
136  \f$ L\f$ functions of \f$ N\f$ variables, and \f$ L \leq M\f$. That is, \f$ f_{jl}\f$ is
137  the term (or function) numbered \f$ l\f$ evaluated at the data point
138  \f$ j\f$. We have to normalise \f$\mathbf{x}_j\f$ to \f$ [-1,1]\f$ for this to
139  succeed [5] (TMultiDimFit::MakeNormalized). We then define a
140  matrix \f$\mathsf{W}\f$ of which the columns \f$\mathbf{w}_j\f$ are given by
141  \f{eqnarray*}{
142  \mathbf{w}_1 &=& \mathbf{f}_1 = F_1\left(\mathbf x_1\right)\\
143  \mathbf{w}_l &=& \mathbf{f}_l - \sum^{l-1}_{k=1} \frac{\mathbf{f}_l \bullet
144  \mathbf{w}_k}{\mathbf{w}_k^2}\mathbf{w}_k\,.
145  \f}
146  and \f$\mathbf{w}_{l}\f$ is the component of \f$\mathbf{f}_{l} \f$ orthogonal
147  to \f$\mathbf{w}_{1}, \ldots, \mathbf{w}_{l-1}\f$. Hence we obtain [3],
148  \f[
149  \mathbf{w}_k\bullet\mathbf{w}_l = 0\quad\mbox{if}~k \neq l\quad.
150  \f]
151  We now take as a new model \f$\mathsf{W}\mathbf{a}\f$. We thus want to
152  minimize
153  \f[
154  S\equiv \left(\mathbf{D} - \mathsf{W}\mathbf{a}\right)^2\quad,
155  \f]
156  where \f$\mathbf{D} = \left(D_1,\ldots,D_M\right)\f$ is a vector of the
157  dependent quantity in the sample. Differentiation with respect to
158  \f$ a_j\f$ gives, using [6], <a name="eq:dS2"></a>
159  \f[
160  \mathbf{D}\bullet\mathbf{w}_l - a_l\mathbf{w}_l^2 = 0
161  \f]
162  or
163  \f[
164  a_l = \frac{\mathbf{D}_l\bullet\mathbf{w}_l}{\mathbf{w}_l^2}
165  \f]
166  Let \f$ S_j\f$ be the sum of squares of residuals when taking \f$ j\f$ functions
167  into account. Then
168  \f[
169  S_l = \left[\mathbf{D} - \sum^l_{k=1} a_k\mathbf{w}_k\right]^2
170  = \mathbf{D}^2 - 2\mathbf{D} \sum^l_{k=1} a_k\mathbf{w}_k
171  + \sum^l_{k=1} a_k^2\mathbf{w}_k^2
172  \f]
173  Using [9], we see that
174  \f{eqnarray*}{
175  S_l &=& \mathbf{D}^2 - 2 \sum^l_{k=1} a_k^2\mathbf{w}_k^2 +
176  \sum^j_{k=1} a_k^2\mathbf{w}_k^2\nonumber\\
177  &=& \mathbf{D}^2 - \sum^l_{k=1} a_k^2\mathbf{w}_k^2\nonumber\\
178  &=& \mathbf{D}^2 - \sum^l_{k=1} \frac{\left(\mathbf D\bullet \mathbf
179  w_k\right)}{\mathbf w_k^2}
180  \f}
181  So for each new function \f$ F_l\f$ included in the model, we get a
182  reduction of the sum of squares of residuals of \f$a_l^2\mathbf{w}_l^2\f$,
183  where \f$\mathbf{w}_l\f$ is given by [4] and \f$ a_l\f$ by [9]. Thus, using
184  the Gram-Schmidt orthogonalisation, we
185  can decide if we want to include this function in the final model,
186  *before* the matrix inversion.
187 
188  ## Function Selection Based on Residual
189  Supposing that \f$ L-1\f$ steps of the procedure have been performed, the
190  problem now is to consider the \f$L^{\mbox{th}}\f$ function.
191 
192  The sum of squares of residuals can be written as
193  \f[
194  S_L = \textbf{D}^T\bullet\textbf{D} -
195  \sum^L_{l=1}a^2_l\left(\textbf{w}_l^T\bullet\textbf{w}_l\right)
196  \f]
197  where the relation [9] have been taken into account. The
198  contribution of the \f$L^{\mbox{th}}\f$ function to the reduction of S, is
199  given by
200  \f[
201  \Delta S_L = a^2_L\left(\textbf{w}_L^T\bullet\textbf{w}_L\right)
202  \f]
203  Two test are now applied to decide whether this \f$L^{\mbox{th}}\f$
204  function is to be included in the final expression, or not.
205 
206  ## Test 1
207  Denoting by \f$ H_{L-1}\f$ the subspace spanned by \f$\textbf{w}_1,\ldots,\textbf{w}_{L-1}\f$
208  the function \d$\textbf{w}_L\d$ is by construction (see 4) the projection of the function
209  \f$ F_L\f$ onto the direction perpendicular to \f$ H_{L-1}\f$. Now, if the
210  length of \f$\textbf{w}_L\f$ (given by \f$\textbf{w}_L\bullet\textbf{w}_L\f$)
211  is very small compared to the length of \f$\textbf{f}_L\f$ this new
212  function can not contribute much to the reduction of the sum of
213  squares of residuals. The test consists then in calculating the angle
214  \f$ \theta \f$ between the two vectors \f$\textbf{w}_L\f$ \f$ \textbf {f}_L\f$
215  (see also figure 1) and requiring that it's
216  *greater* then a threshold value which the user must set (TMultiDimFit::SetMinAngle).
217 
218  \image html multidimfit_img86.gif "Figure 1: (a) angle \f$\theta\f$ between \f$\textbf{w}_l\f$ and \f$\textbf{f}_L\f$, (b) angle \f$ \phi \f$ between \f$\textbf{w}_L\f$ and \f$\textbf{D}\f$"
219 
220  ## Test 2
221  Let \f$\textbf{D}\f$ be the data vector to be fitted. As illustrated in
222  figure 1, the \f$L^{\mbox{th}}\f$ function \f$\textbf{w}_L\f$
223  will contribute significantly to the reduction of \f$ S\f$, if the angle
224  \f$\phi^\prime\f$ between \f$\textbf{w}_L\f$ and \f$\textbf{D}\f$ is smaller than
225  an upper limit \f$ \phi \f$, defined by the user (MultiDimFit::SetMaxAngle)
226 
227  However, the method automatically readjusts the value of this angle
228  while fitting is in progress, in order to make the selection criteria
229  less and less difficult to be fulfilled. The result is that the
230  functions contributing most to the reduction of \f$ S\f$ are chosen first
231  (TMultiDimFit::TestFunction).
232 
233  In case \f$ \phi \f$ isn't defined, an alternative method of
234  performing this second test is used: The \f$L^{\mbox{th}}\f$
235  function \f$\textbf{f}_L\f$ is accepted if (refer also to equation (13))
236  \f[
237  \Delta S_L > \frac{S_{L-1}}{L_{max}-L}
238  \f]
239  where \f$ S_{L-1}\f$ is the sum of the \f$ L-1\f$ first residuals from the
240  \f$ L-1\f$ functions previously accepted; and \f$ L_{max}\f$ is the total number
241  of functions allowed in the final expression of the fit (defined by
242  user).
243 
244  From this we see, that by restricting \f$ L_{max}\f$ -- the number of
245  terms in the final model -- the fit is more difficult to perform,
246  since the above selection criteria is more limiting.
247 
248  The more coefficients we evaluate, the more the sum of squares of
249  residuals \f$ S\f$ will be reduced. We can evaluate \f$ S\f$ before inverting
250  \f$\mathsf{B}\f$ as shown below.
251 
252  ## Coefficients and Coefficient Errors
253  Having found a parameterization, that is the \f$ F_l\f$'s and \f$ L\f$, that
254  minimizes \f$ S\f$, we still need to determine the coefficients
255  \f$ c_l\f$. However, it's a feature of how we choose the significant
256  functions, that the evaluation of the \f$ c_l\f$'s becomes trivial [5]. To derive
257  \f$\mathbf{c}\f$, we first note that
258  equation (4) can be written as
259  \f[
260  \mathsf{F} = \mathsf{W}\mathsf{B}
261  \f]
262  where
263  \f{eqnarray*}{
264  b_{ij} = \frac{\mathbf{f}_j \bullet \mathbf{w}_i}{\mathbf{w}_i^2}
265  & \mbox{if} & i < j\\
266  1 & \mbox{if} & i = j\\
267  0 & \mbox{if} & i > j
268  \f}
269  Consequently, \f$\mathsf{B}\f$ is an upper triangle matrix, which can be
270  readily inverted. So we now evaluate
271  \f[
272  \mathsf{F}\mathsf{B}^{-1} = \mathsf{W}
273  \f]
274  The model \f$\mathsf{W}\mathbf{a}\f$ can therefore be written as
275  \f$(\mathsf{F}\mathsf{B}^{-1})\mathbf{a} = \mathsf{F}(\mathsf{B}^{-1}\mathbf{a})\,.\f$
276 
277  The original model \f$\mathsf{F}\mathbf{c}\f$ is therefore identical with
278  this if
279  \f[
280  \mathbf{c} = \left(\mathsf{B}^{-1}\mathbf{a}\right) =
281  \left[\mathbf{a}^T\left(\mathsf{B}^{-1}\right)^T\right]^T\,.
282  \f]
283  The reason we use \f$\left(\mathsf{B}^{-1}\right)^T\f$ rather then
284  \f$\mathsf{B}^{-1}\f$ is to save storage, since \f$\left(\mathsf{B}^{-1}\right)^T\f$
285  can be stored in the same matrix as \f$\mathsf{B}\f$ (TMultiDimFit::MakeCoefficients).
286  The errors in the coefficients is calculated by inverting the curvature matrix
287  of the non-orthogonal functions \f$ f_{lj}\f$ [1] (TMultiDimFit::MakeCoefficientErrors).
288 
289  ## Considerations
290  It's important to realize that the training sample should be
291  representive of the problem at hand, in particular along the borders
292  of the region of interest. This is because the algorithm presented
293  here, is a *interpolation*, rather then a *extrapolation* [5].
294 
295  Also, the independent variables \f$ x_{i}\f$ need to be linear
296  independent, since the procedure will perform poorly if they are not
297  [5]. One can find an linear transformation from ones
298  original variables \f$ \xi_{i}\f$ to a set of linear independent variables
299  \f$ x_{i}\f$, using a *Principal Components Analysis* (see TPrincipal), and
300  then use the transformed variable as input to this class [5] [6].
301 
302  H. Wind also outlines a method for parameterising a multidimensional
303  dependence over a multidimensional set of variables. An example
304  of the method from [5], is a follows (please refer to
305  [5] for a full discussion):
306 
307  1. Define \f$\mathbf{P} = (P_1, \ldots, P_5)\f$ are the 5 dependent
308  quantities that define a track.
309  2. Compute, for \f$ M\f$ different values of \f$\mathbf{P}\f$, the tracks
310  through the magnetic field, and determine the corresponding
311  \f$\mathbf{x} = (x_1, \ldots, x_N)\f$.
312  3. Use the simulated observations to determine, with a simple
313  approximation, the values of \f$\mathbf{P}_j\f$. We call these values
314  \f$\mathbf{P}^\prime_j, j = 1, \ldots, M\f$.
315  4. Determine from \f$\mathbf{x}\f$ a set of at least five relevant
316  coordinates \f$\mathbf{x}^\prime\f$, using contrains, *or
317  alternative:*
318  5. Perform a Principal Component Analysis (using TPrincipal), and use
319  to get a linear transformation \f$\mathbf{x} \rightarrow \mathbf{x}^\prime\f$, so that
320  \f$\mathbf{x}^\prime\f$ are constrained and linear independent.
321  6. Perform a Principal Component Analysis on
322  \f$Q_i = P_i / P^\prime_i\, i = 1, \ldots, 5\f$, to get linear
323  indenpendent (among themselves, but not independent of \f$\mathbf{x}\f$) quantities
324  \f$\mathbf{Q}^\prime\f$
325  7. For each component \f$Q^\prime_i\f$ make a mutlidimensional fit,
326  using \f$\mathbf{x}^\prime\f$ as the variables, thus determing a set of
327  coefficents \f$\mathbf{c}_i\f$.
328 
329  To process data, using this parameterisation, do
330  1. Test wether the observation \f$\mathbf{x}\f$ within the domain of
331  the parameterization, using the result from the Principal Component
332  Analysis.
333  2. Determine \f$\mathbf{P}^\prime\f$ as before.
334  3. Detetmine \f$\mathbf{x}^\prime\f$ as before.
335  4. Use the result of the fit to determind \f$\mathbf{Q}^\prime\f$.
336  5. Transform back to \f$\mathbf{P}\f$ from \f$\mathbf{Q}^\prime\f$, using
337  the result from the Principal Component Analysis.
338 
339  ## Testing the parameterization
340  The class also provides functionality for testing the, over the
341  training sample, found parameterization (TMultiDimFit::Fit). This is done by passing
342  the class a test sample of \f$ M_t\f$ tuples of the form
343  \f$(\mathbf{x}_{t,j},D_{t,j}, E_{t,j})\f$, where \f$\mathbf{x}_{t,j}\f$ are the independent
344  variables, \f$ D_{t,j}\f$ the known, dependent quantity, and \f$ E_{t,j}\f$ is
345  the square error in \f$ D_{t,j}\f$ (TMultiDimFit::AddTestRow).
346 
347  The parameterization is then evaluated at every \f$\mathbf{x}_t\f$ in the
348  test sample, and
349  \f[
350  S_t \equiv \sum_{j=1}^{M_t} \left(D_{t,j} -
351  D_p\left(\mathbf{x}_{t,j}\right)\right)^2
352  \f]
353  is evaluated. The relative error over the test sample
354  \f[
355  R_t = \frac{S_t}{\sum_{j=1}^{M_t} D_{t,j}^2}
356  \f]
357  should not be to low or high compared to \f$ R\f$ from the training
358  sample. Also, multiple correlation coefficient from both samples should
359  be fairly close, otherwise one of the samples is not representive of
360  the problem. A large difference in the reduced \f$ \chi^2\f$ over the two
361  samples indicate an over fit, and the maximum number of terms in the
362  parameterisation should be reduced.
363 
364  It's possible to use [4] to further improve the fit, using the test sample.
365 
366  Christian Holm
367 
368  ## Bibliography
369  - <a name="bevington"></a> Philip R. Bevington and D. Keith Robinson. *Data Reduction and Error Analysis for
370  the Physical Sciences*. McGraw-Hill, 2 edition, 1992.
371  - <a name="mudifi"></a> R. Brun et al. *Long writeup DD/75-23*, CERN, 1980.
372  - Gene H. Golub and Charles F. van Loan. *Matrix Computations*.
373  John Hopkins Univeristy Press, Baltimore, 3 edition, 1996.
374  - <a name="minuit"></a>F. James. *Minuit*. Long writeup D506, CERN, 1998.
375  - <a name="wind72"></a>H. Wind. *Function parameterization*. Proceedings of the 1972 CERN Computing and Data Processing
376  School, volume 72-21 of Yellow report. CERN, 1972.
377  - <a name="wind81"></a>H. Wind. 1. principal component analysis, 2. pattern recognition for track
378  finding, 3. interpolation and functional representation. Yellow report EP/81-12, CERN, 1981.
379 
380 [1]: classTMultiDimFit.html#bevington
381 [2]: classTMultiDimFit.html#mudifi
382 [4]: classTMultiDimFit.html#minuit
383 [5]: classTMultiDimFit.html#wind72
384 [6]: classTMultiDimFit.html#wind81
385 [9]: classTMultiDimFit.html#eq:dS2
386 */
387 
388 
389 #include "Riostream.h"
390 #include "TMultiDimFit.h"
391 #include "TMath.h"
392 #include "TH1.h"
393 #include "TH2.h"
394 #include "TROOT.h"
395 #include "TBrowser.h"
396 #include "TDecompChol.h"
397 
398 #define RADDEG (180. / TMath::Pi())
399 #define DEGRAD (TMath::Pi() / 180.)
400 #define HIST_XORIG 0
401 #define HIST_DORIG 1
402 #define HIST_XNORM 2
403 #define HIST_DSHIF 3
404 #define HIST_RX 4
405 #define HIST_RD 5
406 #define HIST_RTRAI 6
407 #define HIST_RTEST 7
408 #define PARAM_MAXSTUDY 1
409 #define PARAM_SEVERAL 2
410 #define PARAM_RELERR 3
411 #define PARAM_MAXTERMS 4
412 
413 
414 ////////////////////////////////////////////////////////////////////////////////
415 
416 static void mdfHelper(int&, double*, double&, double*, int);
417 
418 ////////////////////////////////////////////////////////////////////////////////
419 
421 
422 //____________________________________________________________________
423 // Static instance. Used with mdfHelper and TMinuit
424 TMultiDimFit* TMultiDimFit::fgInstance = 0;
425 
426 
427 ////////////////////////////////////////////////////////////////////////////////
428 /// Empty CTOR. Do not use
429 
431 {
432  fMeanQuantity = 0;
433  fMaxQuantity = 0;
434  fMinQuantity = 0;
435  fSumSqQuantity = 0;
436  fSumSqAvgQuantity = 0;
437 
438  fNVariables = 0;
439  fSampleSize = 0;
440  fTestSampleSize = 0;
441 
442  fMinAngle = 1;
443  fMaxAngle = 0;
444  fMaxTerms = 0;
445  fMinRelativeError = 0;
446  fMaxPowers = 0;
447  fPowerLimit = 0;
448 
449  fMaxFunctions = 0;
450  fFunctionCodes = 0;
451  fMaxStudy = 0;
452  fMaxFuncNV = 0;
453 
454  fMaxPowersFinal = 0;
455  fPowers = 0;
456  fPowerIndex = 0;
457 
458  fMaxResidual = 0;
459  fMinResidual = 0;
460  fMaxResidualRow = 0;
461  fMinResidualRow = 0;
462  fSumSqResidual = 0;
463 
464  fNCoefficients = 0;
465  fRMS = 0;
466  fChi2 = 0;
468 
469  fError = 0;
470  fTestError = 0;
471  fPrecision = 0;
472  fTestPrecision = 0;
473  fCorrelationCoeff = 0;
475 
476  fHistograms = 0;
477  fHistogramMask = 0;
478  fBinVarX = 100;
479  fBinVarY = 100;
480 
481  fFitter = 0;
485  fIsVerbose = kFALSE;
486 
487 }
488 
489 
490 ////////////////////////////////////////////////////////////////////////////////
491 /// Constructor
492 /// Second argument is the type of polynomials to use in
493 /// parameterisation, one of:
494 /// TMultiDimFit::kMonomials
495 /// TMultiDimFit::kChebyshev
496 /// TMultiDimFit::kLegendre
497 ///
498 /// Options:
499 /// K Compute (k)correlation matrix
500 /// V Be verbose
501 ///
502 /// Default is no options.
503 ///
504 
507  Option_t *option)
508 : TNamed("multidimfit","Multi-dimensional fit object"),
509 fQuantity(dimension),
510 fSqError(dimension),
511 fVariables(dimension*100),
512 fMeanVariables(dimension),
513 fMaxVariables(dimension),
514 fMinVariables(dimension)
515 {
516  fgInstance = this;
517 
518  fMeanQuantity = 0;
519  fMaxQuantity = 0;
520  fMinQuantity = 0;
521  fSumSqQuantity = 0;
522  fSumSqAvgQuantity = 0;
523 
524  fNVariables = dimension;
525  fSampleSize = 0;
526  fTestSampleSize = 0;
527 
528  fMinAngle = 1;
529  fMaxAngle = 0;
530  fMaxTerms = 0;
531  fMinRelativeError = 0.01;
532  fMaxPowers = new Int_t[dimension];
533  fPowerLimit = 1;
534 
535  fMaxFunctions = 0;
536  fFunctionCodes = 0;
537  fMaxStudy = 0;
538  fMaxFuncNV = 0;
539 
540  fMaxPowersFinal = new Int_t[dimension];
541  fPowers = 0;
542  fPowerIndex = 0;
543 
544  fMaxResidual = 0;
545  fMinResidual = 0;
546  fMaxResidualRow = 0;
547  fMinResidualRow = 0;
548  fSumSqResidual = 0;
549 
550  fNCoefficients = 0;
551  fRMS = 0;
552  fChi2 = 0;
554 
555  fError = 0;
556  fTestError = 0;
557  fPrecision = 0;
558  fTestPrecision = 0;
559  fCorrelationCoeff = 0;
561 
562  fHistograms = 0;
563  fHistogramMask = 0;
564  fBinVarX = 100;
565  fBinVarY = 100;
566 
567  fFitter = 0;
568  fPolyType = type;
571  fIsVerbose = kFALSE;
572  TString opt = option;
573  opt.ToLower();
574 
575  if (opt.Contains("k")) fShowCorrelation = kTRUE;
576  if (opt.Contains("v")) fIsVerbose = kTRUE;
577 }
578 
579 
580 ////////////////////////////////////////////////////////////////////////////////
581 /// Destructor
582 
584 {
585  delete [] fPowers;
586  delete [] fMaxPowers;
587  delete [] fMaxPowersFinal;
588  delete [] fPowerIndex;
589  delete [] fFunctionCodes;
590  if (fHistograms) fHistograms->Clear("nodelete");
591  delete fHistograms;
592 }
593 
594 
595 ////////////////////////////////////////////////////////////////////////////////
596 /// Add a row consisting of fNVariables independent variables, the
597 /// known, dependent quantity, and optionally, the square error in
598 /// the dependent quantity, to the training sample to be used for the
599 /// parameterization.
600 /// The mean of the variables and quantity is calculated on the fly,
601 /// as outlined in TPrincipal::AddRow.
602 /// This sample should be representive of the problem at hand.
603 /// Please note, that if no error is given Poisson statistics is
604 /// assumed and the square error is set to the value of dependent
605 /// quantity. See also the
606 /// Begin_Html<a href="#TMultiDimFit:description">class description</a>End_Html
607 
609 {
610  if (!x)
611  return;
612 
613  if (++fSampleSize == 1) {
614  fMeanQuantity = D;
615  fMaxQuantity = D;
616  fMinQuantity = D;
617  fSumSqQuantity = D * D;// G.Q. erratum on August 15th, 2008
618  }
619  else {
622  fSumSqQuantity += D * D;
623 
624  if (D >= fMaxQuantity) fMaxQuantity = D;
625  if (D <= fMinQuantity) fMinQuantity = D;
626  }
627 
628 
629  // If the vector isn't big enough to hold the new data, then
630  // expand the vector by half it's size.
631  Int_t size = fQuantity.GetNrows();
632  if (fSampleSize > size) {
633  fQuantity.ResizeTo(size + size/2);
634  fSqError.ResizeTo(size + size/2);
635  }
636 
637  // Store the value
638  fQuantity(fSampleSize-1) = D;
639  fSqError(fSampleSize-1) = (E == 0 ? D : E);
640 
641  // Store data point in internal vector
642  // If the vector isn't big enough to hold the new data, then
643  // expand the vector by half it's size
644  size = fVariables.GetNrows();
645  if (fSampleSize * fNVariables > size)
646  fVariables.ResizeTo(size + size/2);
647 
648 
649  // Increment the data point counter
650  Int_t i,j;
651  for (i = 0; i < fNVariables; i++) {
652  if (fSampleSize == 1) {
653  fMeanVariables(i) = x[i];
654  fMaxVariables(i) = x[i];
655  fMinVariables(i) = x[i];
656  }
657  else {
658  fMeanVariables(i) *= 1 - 1./Double_t(fSampleSize);
659  fMeanVariables(i) += x[i] / Double_t(fSampleSize);
660 
661  // Update the maximum value for this component
662  if (x[i] >= fMaxVariables(i)) fMaxVariables(i) = x[i];
663 
664  // Update the minimum value for this component
665  if (x[i] <= fMinVariables(i)) fMinVariables(i) = x[i];
666 
667  }
668 
669  // Store the data.
670  j = (fSampleSize-1) * fNVariables + i;
671  fVariables(j) = x[i];
672  }
673 }
674 
675 
676 ////////////////////////////////////////////////////////////////////////////////
677 /// Add a row consisting of fNVariables independent variables, the
678 /// known, dependent quantity, and optionally, the square error in
679 /// the dependent quantity, to the test sample to be used for the
680 /// test of the parameterization.
681 /// This sample needn't be representive of the problem at hand.
682 /// Please note, that if no error is given Poisson statistics is
683 /// assumed and the square error is set to the value of dependent
684 /// quantity. See also the
685 /// Begin_Html<a href="#TMultiDimFit:description">class description</a>End_Html
686 
688 {
689  if (fTestSampleSize++ == 0) {
693  }
694 
695  // If the vector isn't big enough to hold the new data, then
696  // expand the vector by half it's size.
697  Int_t size = fTestQuantity.GetNrows();
698  if (fTestSampleSize > size) {
699  fTestQuantity.ResizeTo(size + size/2);
700  fTestSqError.ResizeTo(size + size/2);
701  }
702 
703  // Store the value
705  fTestSqError(fTestSampleSize-1) = (E == 0 ? D : E);
706 
707  // Store data point in internal vector
708  // If the vector isn't big enough to hold the new data, then
709  // expand the vector by half it's size
710  size = fTestVariables.GetNrows();
711  if (fTestSampleSize * fNVariables > size)
712  fTestVariables.ResizeTo(size + size/2);
713 
714 
715  // Increment the data point counter
716  Int_t i,j;
717  for (i = 0; i < fNVariables; i++) {
718  j = fNVariables * (fTestSampleSize - 1) + i;
719  fTestVariables(j) = x[i];
720 
721  if (x[i] > fMaxVariables(i))
722  Warning("AddTestRow", "variable %d (row: %d) too large: %f > %f",
723  i, fTestSampleSize, x[i], fMaxVariables(i));
724  if (x[i] < fMinVariables(i))
725  Warning("AddTestRow", "variable %d (row: %d) too small: %f < %f",
726  i, fTestSampleSize, x[i], fMinVariables(i));
727  }
728 }
729 
730 
731 ////////////////////////////////////////////////////////////////////////////////
732 /// Browse the TMultiDimFit object in the TBrowser.
733 
735 {
736  if (fHistograms) {
738  TH1* h = 0;
739  while ((h = (TH1*)next()))
740  b->Add(h,h->GetName());
741  }
742  if (fVariables.IsValid())
743  b->Add(&fVariables, "Variables (Training)");
744  if (fQuantity.IsValid())
745  b->Add(&fQuantity, "Quantity (Training)");
746  if (fSqError.IsValid())
747  b->Add(&fSqError, "Error (Training)");
748  if (fMeanVariables.IsValid())
749  b->Add(&fMeanVariables, "Mean of Variables (Training)");
750  if (fMaxVariables.IsValid())
751  b->Add(&fMaxVariables, "Mean of Variables (Training)");
752  if (fMinVariables.IsValid())
753  b->Add(&fMinVariables, "Min of Variables (Training)");
754  if (fTestVariables.IsValid())
755  b->Add(&fTestVariables, "Variables (Test)");
756  if (fTestQuantity.IsValid())
757  b->Add(&fTestQuantity, "Quantity (Test)");
758  if (fTestSqError.IsValid())
759  b->Add(&fTestSqError, "Error (Test)");
760  if (fFunctions.IsValid())
761  b->Add(&fFunctions, "Functions");
762  if(fCoefficients.IsValid())
763  b->Add(&fCoefficients,"Coefficients");
765  b->Add(&fCoefficientsRMS,"Coefficients Errors");
766  if (fOrthFunctions.IsValid())
767  b->Add(&fOrthFunctions, "Orthogonal Functions");
769  b->Add(&fOrthFunctionNorms, "Orthogonal Functions Norms");
770  if (fResiduals.IsValid())
771  b->Add(&fResiduals, "Residuals");
773  b->Add(&fOrthCoefficients,"Orthogonal Coefficients");
775  b->Add(&fOrthCurvatureMatrix,"Orthogonal curvature matrix");
777  b->Add(&fCorrelationMatrix,"Correlation Matrix");
778  if (fFitter)
779  b->Add(fFitter, fFitter->GetName());
780 }
781 
782 
783 ////////////////////////////////////////////////////////////////////////////////
784 /// Clear internal structures and variables
785 
787 {
788  Int_t i, j, n = fNVariables, m = fMaxFunctions;
789 
790  // Training sample, dependent quantity
791  fQuantity.Zero();
792  fSqError.Zero();
793  fMeanQuantity = 0;
794  fMaxQuantity = 0;
795  fMinQuantity = 0;
796  fSumSqQuantity = 0;
797  fSumSqAvgQuantity = 0;
798 
799  // Training sample, independent variables
800  fVariables.Zero();
801  fNVariables = 0;
802  fSampleSize = 0;
806 
807  // Test sample
809  fTestSqError.Zero();
811  fTestSampleSize = 0;
812 
813  // Functions
814  fFunctions.Zero();
815  //for (i = 0; i < fMaxTerms; i++) fPowerIndex[i] = 0;
816  //for (i = 0; i < fMaxTerms; i++) fFunctionCodes[i] = 0;
817  fMaxFunctions = 0;
818  fMaxStudy = 0;
821 
822  // Control parameters
823  fMinRelativeError = 0;
824  fMinAngle = 0;
825  fMaxAngle = 0;
826  fMaxTerms = 0;
827 
828  // Powers
829  for (i = 0; i < n; i++) {
830  fMaxPowers[i] = 0;
831  fMaxPowersFinal[i] = 0;
832  for (j = 0; j < m; j++)
833  fPowers[i * n + j] = 0;
834  }
835  fPowerLimit = 0;
836 
837  // Residuals
838  fMaxResidual = 0;
839  fMinResidual = 0;
840  fMaxResidualRow = 0;
841  fMinResidualRow = 0;
842  fSumSqResidual = 0;
843 
844  // Fit
845  fNCoefficients = 0;
846  fOrthCoefficients = 0;
848  fRMS = 0;
850  fError = 0;
851  fTestError = 0;
852  fPrecision = 0;
853  fTestPrecision = 0;
854 
855  // Coefficients
858  fResiduals.Zero();
859  fHistograms->Clear(option);
860 
861  // Options
865 }
866 
867 
868 ////////////////////////////////////////////////////////////////////////////////
869 /// Evaluate parameterization at point x. Optional argument coeff is
870 /// a vector of coefficients for the parameterisation, fNCoefficients
871 /// elements long.
872 
873 Double_t TMultiDimFit::Eval(const Double_t *x, const Double_t* coeff) const
874 {
875  Double_t returnValue = fMeanQuantity;
876  Double_t term = 0;
877  Int_t i, j;
878 
879  for (i = 0; i < fNCoefficients; i++) {
880  // Evaluate the ith term in the expansion
881  term = (coeff ? coeff[i] : fCoefficients(i));
882  for (j = 0; j < fNVariables; j++) {
883  // Evaluate the factor (polynomial) in the j-th variable.
884  Int_t p = fPowers[fPowerIndex[i] * fNVariables + j];
885  Double_t y = 1 + 2. / (fMaxVariables(j) - fMinVariables(j))
886  * (x[j] - fMaxVariables(j));
887  term *= EvalFactor(p,y);
888  }
889  // Add this term to the final result
890  returnValue += term;
891  }
892  return returnValue;
893 }
894 
895 
896 ////////////////////////////////////////////////////////////////////////////////
897 /// Evaluate parameterization error at point x. Optional argument coeff is
898 /// a vector of coefficients for the parameterisation, fNCoefficients
899 /// elements long.
900 
902 {
903  Double_t returnValue = 0;
904  Double_t term = 0;
905  Int_t i, j;
906 
907  for (i = 0; i < fNCoefficients; i++) {
908  // std::cout << "Error coef " << i << " -> " << fCoefficientsRMS(i) << std::endl;
909  }
910  for (i = 0; i < fNCoefficients; i++) {
911  // Evaluate the ith term in the expansion
912  term = (coeff ? coeff[i] : fCoefficientsRMS(i));
913  for (j = 0; j < fNVariables; j++) {
914  // Evaluate the factor (polynomial) in the j-th variable.
915  Int_t p = fPowers[fPowerIndex[i] * fNVariables + j];
916  Double_t y = 1 + 2. / (fMaxVariables(j) - fMinVariables(j))
917  * (x[j] - fMaxVariables(j));
918  term *= EvalFactor(p,y);
919  // std::cout << "i,j " << i << ", " << j << " " << p << " " << y << " " << EvalFactor(p,y) << " " << term << std::endl;
920  }
921  // Add this term to the final result
922  returnValue += term*term;
923  // std::cout << " i = " << i << " value = " << returnValue << std::endl;
924  }
925  returnValue = sqrt(returnValue);
926  return returnValue;
927 }
928 
929 
930 ////////////////////////////////////////////////////////////////////////////////
931 /// PRIVATE METHOD:
932 /// Calculate the control parameter from the passed powers
933 
935 {
936  Double_t s = 0;
937  Double_t epsilon = 1e-6; // a small number
938  for (Int_t i = 0; i < fNVariables; i++) {
939  if (fMaxPowers[i] != 1)
940  s += (epsilon + iv[i] - 1) / (epsilon + fMaxPowers[i] - 1);
941  }
942  return s;
943 }
944 
945 ////////////////////////////////////////////////////////////////////////////////
946 /// PRIVATE METHOD:
947 /// Evaluate function with power p at variable value x
948 
950 {
951  Int_t i = 0;
952  Double_t p1 = 1;
953  Double_t p2 = 0;
954  Double_t p3 = 0;
955  Double_t r = 0;
956 
957  switch(p) {
958  case 1:
959  r = 1;
960  break;
961  case 2:
962  r = x;
963  break;
964  default:
965  p2 = x;
966  for (i = 3; i <= p; i++) {
967  p3 = p2 * x;
968  if (fPolyType == kLegendre)
969  p3 = ((2 * i - 3) * p2 * x - (i - 2) * p1) / (i - 1);
970  else if (fPolyType == kChebyshev)
971  p3 = 2 * x * p2 - p1;
972  p1 = p2;
973  p2 = p3;
974  }
975  r = p3;
976  }
977 
978  return r;
979 }
980 
981 
982 ////////////////////////////////////////////////////////////////////////////////
983 /// Find the parameterization
984 ///
985 /// Options:
986 /// None so far
987 ///
988 /// For detailed description of what this entails, please refer to the
989 /// Begin_Html<a href="#TMultiDimFit:description">class description</a>End_Html
990 
992 {
993  MakeNormalized();
994  MakeCandidates();
998  MakeCorrelation();
999 }
1000 
1001 ////////////////////////////////////////////////////////////////////////////////
1002 /// Try to fit the found parameterisation to the test sample.
1003 ///
1004 /// Options
1005 /// M use Minuit to improve coefficients
1006 ///
1007 /// Also, refer to
1008 /// Begin_Html<a href="#TMultiDimFit:description">class description</a>End_Html
1009 
1011 {
1012  Int_t i, j;
1013  Double_t* x = new Double_t[fNVariables];
1014  Double_t sumSqD = 0;
1015  Double_t sumD = 0;
1016  Double_t sumSqR = 0;
1017  Double_t sumR = 0;
1018 
1019  // Calculate the residuals over the test sample
1020  for (i = 0; i < fTestSampleSize; i++) {
1021  for (j = 0; j < fNVariables; j++)
1022  x[j] = fTestVariables(i * fNVariables + j);
1023  Double_t res = fTestQuantity(i) - Eval(x);
1024  sumD += fTestQuantity(i);
1025  sumSqD += fTestQuantity(i) * fTestQuantity(i);
1026  sumR += res;
1027  sumSqR += res * res;
1029  ((TH1D*)fHistograms->FindObject("res_test"))->Fill(res);
1030  }
1031  Double_t dAvg = sumSqD - (sumD * sumD) / fTestSampleSize;
1032  Double_t rAvg = sumSqR - (sumR * sumR) / fTestSampleSize;
1033  fTestCorrelationCoeff = (dAvg - rAvg) / dAvg;
1034  fTestError = sumSqR;
1035  fTestPrecision = sumSqR / sumSqD;
1036 
1037  TString opt(option);
1038  opt.ToLower();
1039 
1040  if (!opt.Contains("m"))
1041  MakeChi2();
1042 
1043  if (fNCoefficients * 50 > fTestSampleSize)
1044  Warning("Fit", "test sample is very small");
1045 
1046  if (!opt.Contains("m")) {
1047  Error("Fit", "invalid option");
1048  delete [] x;
1049  return;
1050  }
1051 
1053  if (!fFitter) {
1054  Error("Fit", "Vannot create Fitter");
1055  delete [] x;
1056  return;
1057  }
1059 
1060  const Int_t maxArgs = 16;
1061  Int_t args = 1;
1062  Double_t* arglist = new Double_t[maxArgs];
1063  arglist[0] = -1;
1064  fFitter->ExecuteCommand("SET PRINT",arglist,args);
1065 
1066  for (i = 0; i < fNCoefficients; i++) {
1067  Double_t startVal = fCoefficients(i);
1068  Double_t startErr = fCoefficientsRMS(i);
1069  fFitter->SetParameter(i, Form("coeff%02d",i),
1070  startVal, startErr, 0, 0);
1071  }
1072 
1073  // arglist[0] = 0;
1074  args = 1;
1075  // fFitter->ExecuteCommand("SET PRINT",arglist,args);
1076  fFitter->ExecuteCommand("MIGRAD",arglist,args);
1077 
1078  for (i = 0; i < fNCoefficients; i++) {
1079  Double_t val = 0, err = 0, low = 0, high = 0;
1080  fFitter->GetParameter(i, Form("coeff%02d",i),
1081  val, err, low, high);
1082  fCoefficients(i) = val;
1083  fCoefficientsRMS(i) = err;
1084  }
1085  delete [] x;
1086 }
1087 
1088 ////////////////////////////////////////////////////////////////////////////////
1089 /// Return the static instance.
1090 
1092 {
1093  return fgInstance;
1094 }
1095 
1096 ////////////////////////////////////////////////////////////////////////////////
1097 /// PRIVATE METHOD:
1098 /// Create list of candidate functions for the parameterisation. See
1099 /// also
1100 /// Begin_Html<a href="#TMultiDimFit:description">class description</a>End_Html
1101 
1103 {
1104  Int_t i = 0;
1105  Int_t j = 0;
1106  Int_t k = 0;
1107 
1108  // The temporary array to store the powers in. We don't need to
1109  // initialize this array however.
1111  Int_t *powers = new Int_t[fMaxFuncNV];
1112 
1113  // store of `control variables'
1114  Double_t* control = new Double_t[fMaxFunctions];
1115 
1116  // We've better initialize the variables
1117  Int_t *iv = new Int_t[fNVariables];
1118  for (i = 0; i < fNVariables; i++)
1119  iv[i] = 1;
1120 
1121  if (!fIsUserFunction) {
1122 
1123  // Number of funcs selected
1124  Int_t numberFunctions = 0;
1125 
1126  // Absolute max number of functions
1127  Int_t maxNumberFunctions = 1;
1128  for (i = 0; i < fNVariables; i++)
1129  maxNumberFunctions *= fMaxPowers[i];
1130 
1131  while (kTRUE) {
1132  // Get the control value for this function
1133  Double_t s = EvalControl(iv);
1134 
1135  if (s <= fPowerLimit) {
1136 
1137  // Call over-loadable method Select, as to allow the user to
1138  // interfere with the selection of functions.
1139  if (Select(iv)) {
1140  numberFunctions++;
1141 
1142  // If we've reached the user defined limit of how many
1143  // functions we can consider, break out of the loop
1144  if (numberFunctions > fMaxFunctions)
1145  break;
1146 
1147  // Store the control value, so we can sort array of powers
1148  // later on
1149  control[numberFunctions-1] = Int_t(1.0e+6*s);
1150 
1151  // Store the powers in powers array.
1152  for (i = 0; i < fNVariables; i++) {
1153  j = (numberFunctions - 1) * fNVariables + i;
1154  powers[j] = iv[i];
1155  }
1156  } // if (Select())
1157  } // if (s <= fPowerLimit)
1158 
1159  for (i = 0; i < fNVariables; i++)
1160  if (iv[i] < fMaxPowers[i])
1161  break;
1162 
1163  // If all variables have reached their maximum power, then we
1164  // break out of the loop
1165  if (i == fNVariables) {
1166  fMaxFunctions = numberFunctions;
1167  break;
1168  }
1169 
1170  // Next power in variable i
1171  if (i < fNVariables) iv[i]++;
1172 
1173  for (j = 0; j < i; j++)
1174  iv[j] = 1;
1175  } // while (kTRUE)
1176  }
1177  else {
1178  // In case the user gave an explicit function
1179  for (i = 0; i < fMaxFunctions; i++) {
1180  // Copy the powers to working arrays
1181  for (j = 0; j < fNVariables; j++) {
1182  powers[i * fNVariables + j] = fPowers[i * fNVariables + j];
1183  iv[j] = fPowers[i * fNVariables + j];
1184  }
1185 
1186  control[i] = Int_t(1.0e+6*EvalControl(iv));
1187  }
1188  }
1189 
1190  // Now we need to sort the powers according to least `control
1191  // variable'
1192  Int_t *order = new Int_t[fMaxFunctions];
1193  for (i = 0; i < fMaxFunctions; i++)
1194  order[i] = i;
1195  fMaxFuncNV = fMaxFunctions * fNVariables;
1196  fPowers = new Int_t[fMaxFuncNV];
1197 
1198  for (i = 0; i < fMaxFunctions; i++) {
1199  Double_t x = control[i];
1200  Int_t l = order[i];
1201  k = i;
1202 
1203  for (j = i; j < fMaxFunctions; j++) {
1204  if (control[j] <= x) {
1205  x = control[j];
1206  l = order[j];
1207  k = j;
1208  }
1209  }
1210 
1211  if (k != i) {
1212  control[k] = control[i];
1213  control[i] = x;
1214  order[k] = order[i];
1215  order[i] = l;
1216  }
1217  }
1218 
1219  for (i = 0; i < fMaxFunctions; i++)
1220  for (j = 0; j < fNVariables; j++)
1221  fPowers[i * fNVariables + j] = powers[order[i] * fNVariables + j];
1222 
1223  delete [] control;
1224  delete [] powers;
1225  delete [] order;
1226  delete [] iv;
1227 }
1228 
1229 
1230 ////////////////////////////////////////////////////////////////////////////////
1231 /// Calculate Chi square over either the test sample. The optional
1232 /// argument coeff is a vector of coefficients to use in the
1233 /// evaluation of the parameterisation. If coeff == 0, then the found
1234 /// coefficients is used.
1235 /// Used my MINUIT for fit (see TMultDimFit::Fit)
1236 
1238 {
1239  fChi2 = 0;
1240  Int_t i, j;
1241  Double_t* x = new Double_t[fNVariables];
1242  for (i = 0; i < fTestSampleSize; i++) {
1243  // Get the stored point
1244  for (j = 0; j < fNVariables; j++)
1245  x[j] = fTestVariables(i * fNVariables + j);
1246 
1247  // Evaluate function. Scale to shifted values
1248  Double_t f = Eval(x,coeff);
1249 
1250  // Calculate contribution to Chic square
1251  fChi2 += 1. / TMath::Max(fTestSqError(i),1e-20)
1252  * (fTestQuantity(i) - f) * (fTestQuantity(i) - f);
1253  }
1254 
1255  // Clean up
1256  delete [] x;
1257 
1258  return fChi2;
1259 }
1260 
1261 
1262 ////////////////////////////////////////////////////////////////////////////////
1263 /// Generate the file <filename> with .C appended if argument doesn't
1264 /// end in .cxx or .C. The contains the implementation of the
1265 /// function:
1266 ///
1267 /// Double_t <funcname>(Double_t *x)
1268 ///
1269 /// which does the same as TMultiDimFit::Eval. Please refer to this
1270 /// method.
1271 ///
1272 /// Further, the static variables:
1273 ///
1274 /// Int_t gNVariables
1275 /// Int_t gNCoefficients
1276 /// Double_t gDMean
1277 /// Double_t gXMean[]
1278 /// Double_t gXMin[]
1279 /// Double_t gXMax[]
1280 /// Double_t gCoefficient[]
1281 /// Int_t gPower[]
1282 ///
1283 /// are initialized. The only ROOT header file needed is Rtypes.h
1284 ///
1285 /// See TMultiDimFit::MakeRealCode for a list of options
1286 
1287 void TMultiDimFit::MakeCode(const char* filename, Option_t *option)
1288 {
1289 
1290  TString outName(filename);
1291  if (!outName.EndsWith(".C") && !outName.EndsWith(".cxx"))
1292  outName += ".C";
1293 
1294  MakeRealCode(outName.Data(),"",option);
1295 }
1296 
1297 
1298 
1299 ////////////////////////////////////////////////////////////////////////////////
1300 /// PRIVATE METHOD:
1301 /// Compute the errors on the coefficients. For this to be done, the
1302 /// curvature matrix of the non-orthogonal functions, is computed.
1303 
1305 {
1306  Int_t i = 0;
1307  Int_t j = 0;
1308  Int_t k = 0;
1309  TVectorD iF(fSampleSize);
1310  TVectorD jF(fSampleSize);
1312 
1313  TMatrixDSym curvatureMatrix(fNCoefficients);
1314 
1315  // Build the curvature matrix
1316  for (i = 0; i < fNCoefficients; i++) {
1317  iF = TMatrixDRow(fFunctions,i);
1318  for (j = 0; j <= i; j++) {
1319  jF = TMatrixDRow(fFunctions,j);
1320  for (k = 0; k < fSampleSize; k++)
1321  curvatureMatrix(i,j) +=
1322  1 / TMath::Max(fSqError(k), 1e-20) * iF(k) * jF(k);
1323  curvatureMatrix(j,i) = curvatureMatrix(i,j);
1324  }
1325  }
1326 
1327  // Calculate Chi Square
1328  fChi2 = 0;
1329  for (i = 0; i < fSampleSize; i++) {
1330  Double_t f = 0;
1331  for (j = 0; j < fNCoefficients; j++)
1332  f += fCoefficients(j) * fFunctions(j,i);
1333  fChi2 += 1. / TMath::Max(fSqError(i),1e-20) * (fQuantity(i) - f)
1334  * (fQuantity(i) - f);
1335  }
1336 
1337  // Invert the curvature matrix
1338  const TVectorD diag = TMatrixDDiag_const(curvatureMatrix);
1339  curvatureMatrix.NormByDiag(diag);
1340 
1341  TDecompChol chol(curvatureMatrix);
1342  if (!chol.Decompose())
1343  Error("MakeCoefficientErrors", "curvature matrix is singular");
1344  chol.Invert(curvatureMatrix);
1345 
1346  curvatureMatrix.NormByDiag(diag);
1347 
1348  for (i = 0; i < fNCoefficients; i++)
1349  fCoefficientsRMS(i) = TMath::Sqrt(curvatureMatrix(i,i));
1350 }
1351 
1352 
1353 ////////////////////////////////////////////////////////////////////////////////
1354 /// PRIVATE METHOD:
1355 /// Invert the model matrix B, and compute final coefficients. For a
1356 /// more thorough discussion of what this means, please refer to the
1357 /// Begin_Html<a href="#TMultiDimFit:description">class description</a>End_Html
1358 ///
1359 /// First we invert the lower triangle matrix fOrthCurvatureMatrix
1360 /// and store the inverted matrix in the upper triangle.
1361 
1363 {
1364  Int_t i = 0, j = 0;
1365  Int_t col = 0, row = 0;
1366 
1367  // Invert the B matrix
1368  for (col = 1; col < fNCoefficients; col++) {
1369  for (row = col - 1; row > -1; row--) {
1370  fOrthCurvatureMatrix(row,col) = 0;
1371  for (i = row; i <= col ; i++)
1372  fOrthCurvatureMatrix(row,col) -=
1373  fOrthCurvatureMatrix(i,row)
1374  * fOrthCurvatureMatrix(i,col);
1375  }
1376  }
1377 
1378  // Compute the final coefficients
1379  fCoefficients.ResizeTo(fNCoefficients);
1380 
1381  for (i = 0; i < fNCoefficients; i++) {
1382  Double_t sum = 0;
1383  for (j = i; j < fNCoefficients; j++)
1384  sum += fOrthCurvatureMatrix(i,j) * fOrthCoefficients(j);
1385  fCoefficients(i) = sum;
1386  }
1387 
1388  // Compute the final residuals
1390  for (i = 0; i < fSampleSize; i++)
1391  fResiduals(i) = fQuantity(i);
1392 
1393  for (i = 0; i < fNCoefficients; i++)
1394  for (j = 0; j < fSampleSize; j++)
1395  fResiduals(j) -= fCoefficients(i) * fFunctions(i,j);
1396 
1397  // Compute the max and minimum, and squared sum of the evaluated
1398  // residuals
1399  fMinResidual = 10e10;
1400  fMaxResidual = -10e10;
1401  Double_t sqRes = 0;
1402  for (i = 0; i < fSampleSize; i++){
1403  sqRes += fResiduals(i) * fResiduals(i);
1404  if (fResiduals(i) <= fMinResidual) {
1405  fMinResidual = fResiduals(i);
1406  fMinResidualRow = i;
1407  }
1408  if (fResiduals(i) >= fMaxResidual) {
1409  fMaxResidual = fResiduals(i);
1410  fMaxResidualRow = i;
1411  }
1412  }
1413 
1416 
1417  // If we use histograms, fill some more
1421  for (i = 0; i < fSampleSize; i++) {
1423  ((TH2D*)fHistograms->FindObject("res_d"))->Fill(fQuantity(i),
1424  fResiduals(i));
1426  ((TH1D*)fHistograms->FindObject("res_train"))->Fill(fResiduals(i));
1427 
1429  for (j = 0; j < fNVariables; j++)
1430  ((TH2D*)fHistograms->FindObject(Form("res_x_%d",j)))
1431  ->Fill(fVariables(i * fNVariables + j),fResiduals(i));
1432  }
1433  } // If histograms
1434 
1435 }
1436 
1437 
1438 ////////////////////////////////////////////////////////////////////////////////
1439 /// PRIVATE METHOD:
1440 /// Compute the correlation matrix
1441 
1443 {
1444  if (!fShowCorrelation)
1445  return;
1446 
1448 
1449  Double_t d2 = 0;
1450  Double_t ddotXi = 0; // G.Q. needs to be reinitialized in the loop over i fNVariables
1451  Double_t xiNorm = 0; // G.Q. needs to be reinitialized in the loop over i fNVariables
1452  Double_t xidotXj = 0; // G.Q. needs to be reinitialized in the loop over j fNVariables
1453  Double_t xjNorm = 0; // G.Q. needs to be reinitialized in the loop over j fNVariables
1454 
1455  Int_t i, j, k, l, m; // G.Q. added m variable
1456  for (i = 0; i < fSampleSize; i++)
1457  d2 += fQuantity(i) * fQuantity(i);
1458 
1459  for (i = 0; i < fNVariables; i++) {
1460  ddotXi = 0.; // G.Q. reinitialisation
1461  xiNorm = 0.; // G.Q. reinitialisation
1462  for (j = 0; j< fSampleSize; j++) {
1463  // Index of sample j of variable i
1464  k = j * fNVariables + i;
1465  ddotXi += fQuantity(j) * (fVariables(k) - fMeanVariables(i));
1466  xiNorm += (fVariables(k) - fMeanVariables(i))
1467  * (fVariables(k) - fMeanVariables(i));
1468  }
1469  fCorrelationMatrix(i,0) = ddotXi / TMath::Sqrt(d2 * xiNorm);
1470 
1471  for (j = 0; j < i; j++) {
1472  xidotXj = 0.; // G.Q. reinitialisation
1473  xjNorm = 0.; // G.Q. reinitialisation
1474  for (k = 0; k < fSampleSize; k++) {
1475  // Index of sample j of variable i
1476  // l = j * fNVariables + k; // G.Q.
1477  l = k * fNVariables + j; // G.Q.
1478  m = k * fNVariables + i; // G.Q.
1479  // G.Q. xidotXj += (fVariables(i) - fMeanVariables(i))
1480  // G.Q. * (fVariables(l) - fMeanVariables(j));
1481  xidotXj += (fVariables(m) - fMeanVariables(i))
1482  * (fVariables(l) - fMeanVariables(j)); // G.Q. modified index for Xi
1483  xjNorm += (fVariables(l) - fMeanVariables(j))
1484  * (fVariables(l) - fMeanVariables(j));
1485  }
1486  //fCorrelationMatrix(i+1,j) = xidotXj / TMath::Sqrt(xiNorm * xjNorm);
1487  fCorrelationMatrix(i,j+1) = xidotXj / TMath::Sqrt(xiNorm * xjNorm);
1488  }
1489  }
1490 }
1491 
1492 
1493 
1494 ////////////////////////////////////////////////////////////////////////////////
1495 /// PRIVATE METHOD:
1496 /// Make Gram-Schmidt orthogonalisation. The class description gives
1497 /// a thorough account of this algorithm, as well as
1498 /// references. Please refer to the
1499 /// Begin_Html<a href="#TMultiDimFit:description">class description</a>End_Html
1500 
1502 {
1503 
1504  // calculate w_i, that is, evaluate the current function at data
1505  // point i
1506  Double_t f2 = 0;
1509  Int_t j = 0;
1510  Int_t k = 0;
1511 
1512  for (j = 0; j < fSampleSize; j++) {
1513  fFunctions(fNCoefficients, j) = 1;
1515  // First, however, we need to calculate f_fNCoefficients
1516  for (k = 0; k < fNVariables; k++) {
1517  Int_t p = fPowers[function * fNVariables + k];
1518  Double_t x = fVariables(j * fNVariables + k);
1519  fFunctions(fNCoefficients, j) *= EvalFactor(p,x);
1520  }
1521 
1522  // Calculate f dot f in f2
1524  // Assign to w_fNCoefficients f_fNCoefficients
1526  }
1527 
1528  // the first column of w is equal to f
1529  for (j = 0; j < fNCoefficients; j++) {
1530  Double_t fdw = 0;
1531  // Calculate (f_fNCoefficients dot w_j) / w_j^2
1532  for (k = 0; k < fSampleSize; k++) {
1533  fdw += fFunctions(fNCoefficients, k) * fOrthFunctions(j,k)
1534  / fOrthFunctionNorms(j);
1535  }
1536 
1537  fOrthCurvatureMatrix(fNCoefficients,j) = fdw;
1538  // and subtract it from the current value of w_ij
1539  for (k = 0; k < fSampleSize; k++)
1540  fOrthFunctions(fNCoefficients,k) -= fdw * fOrthFunctions(j,k);
1541  }
1542 
1543  for (j = 0; j < fSampleSize; j++) {
1544  // calculate squared length of w_fNCoefficients
1545  fOrthFunctionNorms(fNCoefficients) +=
1546  fOrthFunctions(fNCoefficients,j)
1547  * fOrthFunctions(fNCoefficients,j);
1548 
1549  // calculate D dot w_fNCoefficients in A
1550  fOrthCoefficients(fNCoefficients) += fQuantity(j)
1551  * fOrthFunctions(fNCoefficients, j);
1552  }
1553 
1554  // First test, but only if didn't user specify
1555  if (!fIsUserFunction)
1556  if (TMath::Sqrt(fOrthFunctionNorms(fNCoefficients) / (f2 + 1e-10))
1558  return 0;
1559 
1560  // The result found by this code for the first residual is always
1561  // much less then the one found be MUDIFI. That's because it's
1562  // supposed to be. The cause is the improved precision of Double_t
1563  // over DOUBLE PRECISION!
1564  fOrthCurvatureMatrix(fNCoefficients,fNCoefficients) = 1;
1565  Double_t b = fOrthCoefficients(fNCoefficients);
1566  fOrthCoefficients(fNCoefficients) /= fOrthFunctionNorms(fNCoefficients);
1567 
1568  // Calculate the residual from including this fNCoefficients.
1569  Double_t dResidur = fOrthCoefficients(fNCoefficients) * b;
1570 
1571  return dResidur;
1572 }
1573 
1574 
1575 ////////////////////////////////////////////////////////////////////////////////
1576 /// Make histograms of the result of the analysis. This message
1577 /// should be sent after having read all data points, but before
1578 /// finding the parameterization
1579 ///
1580 /// Options:
1581 /// A All the below
1582 /// X Original independent variables
1583 /// D Original dependent variables
1584 /// N Normalised independent variables
1585 /// S Shifted dependent variables
1586 /// R1 Residuals versus normalised independent variables
1587 /// R2 Residuals versus dependent variable
1588 /// R3 Residuals computed on training sample
1589 /// R4 Residuals computed on test sample
1590 ///
1591 /// For a description of these quantities, refer to
1592 /// Begin_Html<a href="#TMultiDimFit:description">class description</a>End_Html
1593 
1595 {
1596  TString opt(option);
1597  opt.ToLower();
1598 
1599  if (opt.Length() < 1)
1600  return;
1601 
1602  if (!fHistograms)
1603  fHistograms = new TList;
1604 
1605  // Counter variable
1606  Int_t i = 0;
1607 
1608  // Histogram of original variables
1609  if (opt.Contains("x") || opt.Contains("a")) {
1611  for (i = 0; i < fNVariables; i++)
1612  if (!fHistograms->FindObject(Form("x_%d_orig",i)))
1613  fHistograms->Add(new TH1D(Form("x_%d_orig",i),
1614  Form("Original variable # %d",i),
1615  fBinVarX, fMinVariables(i),
1616  fMaxVariables(i)));
1617  }
1618 
1619  // Histogram of original dependent variable
1620  if (opt.Contains("d") || opt.Contains("a")) {
1622  if (!fHistograms->FindObject("d_orig"))
1623  fHistograms->Add(new TH1D("d_orig", "Original Quantity",
1625  }
1626 
1627  // Histograms of normalized variables
1628  if (opt.Contains("n") || opt.Contains("a")) {
1630  for (i = 0; i < fNVariables; i++)
1631  if (!fHistograms->FindObject(Form("x_%d_norm",i)))
1632  fHistograms->Add(new TH1D(Form("x_%d_norm",i),
1633  Form("Normalized variable # %d",i),
1634  fBinVarX, -1,1));
1635  }
1636 
1637  // Histogram of shifted dependent variable
1638  if (opt.Contains("s") || opt.Contains("a")) {
1640  if (!fHistograms->FindObject("d_shifted"))
1641  fHistograms->Add(new TH1D("d_shifted", "Shifted Quantity",
1644  }
1645 
1646  // Residual from training sample versus independent variables
1647  if (opt.Contains("r1") || opt.Contains("a")) {
1649  for (i = 0; i < fNVariables; i++)
1650  if (!fHistograms->FindObject(Form("res_x_%d",i)))
1651  fHistograms->Add(new TH2D(Form("res_x_%d",i),
1652  Form("Computed residual versus x_%d", i),
1653  fBinVarX, -1, 1,
1654  fBinVarY,
1657  }
1658 
1659  // Residual from training sample versus. dependent variable
1660  if (opt.Contains("r2") || opt.Contains("a")) {
1662  if (!fHistograms->FindObject("res_d"))
1663  fHistograms->Add(new TH2D("res_d",
1664  "Computed residuals vs Quantity",
1665  fBinVarX,
1668  fBinVarY,
1671  }
1672 
1673  // Residual from training sample
1674  if (opt.Contains("r3") || opt.Contains("a")) {
1676  if (!fHistograms->FindObject("res_train"))
1677  fHistograms->Add(new TH1D("res_train",
1678  "Computed residuals over training sample",
1681 
1682  }
1683  if (opt.Contains("r4") || opt.Contains("a")) {
1685  if (!fHistograms->FindObject("res_test"))
1686  fHistograms->Add(new TH1D("res_test",
1687  "Distribution of residuals from test",
1690  }
1691 }
1692 
1693 
1694 ////////////////////////////////////////////////////////////////////////////////
1695 /// Generate the file <classname>MDF.cxx which contains the
1696 /// implementation of the method:
1697 ///
1698 /// Double_t <classname>::MDF(Double_t *x)
1699 ///
1700 /// which does the same as TMultiDimFit::Eval. Please refer to this
1701 /// method.
1702 ///
1703 /// Further, the public static members:
1704 ///
1705 /// Int_t <classname>::fgNVariables
1706 /// Int_t <classname>::fgNCoefficients
1707 /// Double_t <classname>::fgDMean
1708 /// Double_t <classname>::fgXMean[] //[fgNVariables]
1709 /// Double_t <classname>::fgXMin[] //[fgNVariables]
1710 /// Double_t <classname>::fgXMax[] //[fgNVariables]
1711 /// Double_t <classname>::fgCoefficient[] //[fgNCoeffficents]
1712 /// Int_t <classname>::fgPower[] //[fgNCoeffficents*fgNVariables]
1713 ///
1714 /// are initialized, and assumed to exist. The class declaration is
1715 /// assumed to be in <classname>.h and assumed to be provided by the
1716 /// user.
1717 ///
1718 /// See TMultiDimFit::MakeRealCode for a list of options
1719 ///
1720 /// The minimal class definition is:
1721 ///
1722 /// class <classname> {
1723 /// public:
1724 /// Int_t <classname>::fgNVariables; // Number of variables
1725 /// Int_t <classname>::fgNCoefficients; // Number of terms
1726 /// Double_t <classname>::fgDMean; // Mean from training sample
1727 /// Double_t <classname>::fgXMean[]; // Mean from training sample
1728 /// Double_t <classname>::fgXMin[]; // Min from training sample
1729 /// Double_t <classname>::fgXMax[]; // Max from training sample
1730 /// Double_t <classname>::fgCoefficient[]; // Coefficients
1731 /// Int_t <classname>::fgPower[]; // Function powers
1732 ///
1733 /// Double_t Eval(Double_t *x);
1734 /// };
1735 ///
1736 /// Whether the method <classname>::Eval should be static or not, is
1737 /// up to the user.
1738 
1739 void TMultiDimFit::MakeMethod(const Char_t* classname, Option_t* option)
1740 {
1741  MakeRealCode(Form("%sMDF.cxx", classname), classname, option);
1742 }
1743 
1744 
1745 
1746 ////////////////////////////////////////////////////////////////////////////////
1747 /// PRIVATE METHOD:
1748 /// Normalize data to the interval [-1;1]. This is needed for the
1749 /// classes method to work.
1750 
1752 {
1753  Int_t i = 0;
1754  Int_t j = 0;
1755  Int_t k = 0;
1756 
1757  for (i = 0; i < fSampleSize; i++) {
1759  ((TH1D*)fHistograms->FindObject("d_orig"))->Fill(fQuantity(i));
1760 
1761  fQuantity(i) -= fMeanQuantity;
1763 
1765  ((TH1D*)fHistograms->FindObject("d_shifted"))->Fill(fQuantity(i));
1766 
1767  for (j = 0; j < fNVariables; j++) {
1768  Double_t range = 1. / (fMaxVariables(j) - fMinVariables(j));
1769  k = i * fNVariables + j;
1770 
1771  // Fill histograms of original independent variables
1773  ((TH1D*)fHistograms->FindObject(Form("x_%d_orig",j)))
1774  ->Fill(fVariables(k));
1775 
1776  // Normalise independent variables
1777  fVariables(k) = 1 + 2 * range * (fVariables(k) - fMaxVariables(j));
1778 
1779  // Fill histograms of normalised independent variables
1781  ((TH1D*)fHistograms->FindObject(Form("x_%d_norm",j)))
1782  ->Fill(fVariables(k));
1783 
1784  }
1785  }
1786  // Shift min and max of dependent variable
1789 
1790  // Shift mean of independent variables
1791  for (i = 0; i < fNVariables; i++) {
1792  Double_t range = 1. / (fMaxVariables(i) - fMinVariables(i));
1793  fMeanVariables(i) = 1 + 2 * range * (fMeanVariables(i)
1794  - fMaxVariables(i));
1795  }
1796 }
1797 
1798 
1799 ////////////////////////////////////////////////////////////////////////////////
1800 /// PRIVATE METHOD:
1801 /// Find the parameterization over the training sample. A full account
1802 /// of the algorithm is given in the
1803 /// Begin_Html<a href="#TMultiDimFit:description">class description</a>End_Html
1804 
1806 {
1807  Int_t i = -1;
1808  Int_t j = 0;
1809  Int_t k = 0;
1810  Int_t maxPass = 3;
1811  Int_t studied = 0;
1812  Double_t squareResidual = fSumSqAvgQuantity;
1813  fNCoefficients = 0;
1820  fFunctions = 1;
1821 
1823  fPowerIndex = new Int_t[fMaxTerms];
1824  Int_t l;
1825  for (l=0;l<fMaxFunctions;l++) fFunctionCodes[l] = 0;
1826  for (l=0;l<fMaxTerms;l++) fPowerIndex[l] = 0;
1827 
1828  if (fMaxAngle != 0) maxPass = 100;
1829  if (fIsUserFunction) maxPass = 1;
1830 
1831  // Loop over the number of functions we want to study.
1832  // increment inspection counter
1833  while(kTRUE) {
1834 
1835  // Reach user defined limit of studies
1836  if (studied++ >= fMaxStudy) {
1838  break;
1839  }
1840 
1841  // Considered all functions several times
1842  if (k >= maxPass) {
1844  break;
1845  }
1846 
1847  // increment function counter
1848  i++;
1849 
1850  // If we've reached the end of the functions, restart pass
1851  if (i == fMaxFunctions) {
1852  if (fMaxAngle != 0)
1853  fMaxAngle += (90 - fMaxAngle) / 2;
1854  i = 0;
1855  studied--;
1856  k++;
1857  continue;
1858  }
1859  if (studied == 1)
1860  fFunctionCodes[i] = 0;
1861  else if (fFunctionCodes[i] >= 2)
1862  continue;
1863 
1864  // Print a happy message
1865  if (fIsVerbose && studied == 1)
1866  std::cout << "Coeff SumSqRes Contrib Angle QM Func"
1867  << " Value W^2 Powers" << std::endl;
1868 
1869  // Make the Gram-Schmidt
1870  Double_t dResidur = MakeGramSchmidt(i);
1871 
1872  if (dResidur == 0) {
1873  // This function is no good!
1874  // First test is in MakeGramSchmidt
1875  fFunctionCodes[i] = 1;
1876  continue;
1877  }
1878 
1879  // If user specified function, assume they know what they are doing
1880  if (!fIsUserFunction) {
1881  // Flag this function as considered
1882  fFunctionCodes[i] = 2;
1883 
1884  // Test if this function contributes to the fit
1885  if (!TestFunction(squareResidual, dResidur)) {
1886  fFunctionCodes[i] = 1;
1887  continue;
1888  }
1889  }
1890 
1891  // If we get to here, the function currently considered is
1892  // fNCoefficients, so we increment the counter
1893  // Flag this function as OK, and store and the number in the
1894  // index.
1895  fFunctionCodes[i] = 3;
1897  fNCoefficients++;
1898 
1899  // We add the current contribution to the sum of square of
1900  // residuals;
1901  squareResidual -= dResidur;
1902 
1903 
1904  // Calculate control parameter from this function
1905  for (j = 0; j < fNVariables; j++) {
1906  if (fNCoefficients == 1
1907  || fMaxPowersFinal[j] <= fPowers[i * fNVariables + j] - 1)
1908  fMaxPowersFinal[j] = fPowers[i * fNVariables + j] - 1;
1909  }
1910  Double_t s = EvalControl(&fPowers[i * fNVariables]);
1911 
1912  // Print the statistics about this function
1913  if (fIsVerbose) {
1914  std::cout << std::setw(5) << fNCoefficients << " "
1915  << std::setw(10) << std::setprecision(4) << squareResidual << " "
1916  << std::setw(10) << std::setprecision(4) << dResidur << " "
1917  << std::setw(7) << std::setprecision(3) << fMaxAngle << " "
1918  << std::setw(7) << std::setprecision(3) << s << " "
1919  << std::setw(5) << i << " "
1920  << std::setw(10) << std::setprecision(4)
1921  << fOrthCoefficients(fNCoefficients-1) << " "
1922  << std::setw(10) << std::setprecision(4)
1923  << fOrthFunctionNorms(fNCoefficients-1) << " "
1924  << std::flush;
1925  for (j = 0; j < fNVariables; j++)
1926  std::cout << " " << fPowers[i * fNVariables + j] - 1 << std::flush;
1927  std::cout << std::endl;
1928  }
1929 
1930  if (fNCoefficients >= fMaxTerms /* && fIsVerbose */) {
1932  break;
1933  }
1934 
1935  Double_t err = TMath::Sqrt(TMath::Max(1e-20,squareResidual) /
1937  if (err < fMinRelativeError) {
1939  break;
1940  }
1941 
1942  }
1943 
1944  fError = TMath::Max(1e-20,squareResidual);
1947 }
1948 
1949 
1950 ////////////////////////////////////////////////////////////////////////////////
1951 /// PRIVATE METHOD:
1952 /// This is the method that actually generates the code for the
1953 /// evaluation the parameterization on some point.
1954 /// It's called by TMultiDimFit::MakeCode and TMultiDimFit::MakeMethod.
1955 ///
1956 /// The options are: NONE so far
1957 
1959  const char *classname,
1960  Option_t *)
1961 {
1962  Int_t i, j;
1963 
1964  Bool_t isMethod = (classname[0] == '\0' ? kFALSE : kTRUE);
1965  const char *prefix = (isMethod ? Form("%s::", classname) : "");
1966  const char *cv_qual = (isMethod ? "" : "static ");
1967 
1968  std::ofstream outFile(filename,std::ios::out|std::ios::trunc);
1969  if (!outFile) {
1970  Error("MakeRealCode","couldn't open output file '%s'",filename);
1971  return;
1972  }
1973 
1974  if (fIsVerbose)
1975  std::cout << "Writing on file \"" << filename << "\" ... " << std::flush;
1976  //
1977  // Write header of file
1978  //
1979  // Emacs mode line ;-)
1980  outFile << "// -*- mode: c++ -*-" << std::endl;
1981  // Info about creator
1982  outFile << "// " << std::endl
1983  << "// File " << filename
1984  << " generated by TMultiDimFit::MakeRealCode" << std::endl;
1985  // Time stamp
1986  TDatime date;
1987  outFile << "// on " << date.AsString() << std::endl;
1988  // ROOT version info
1989  outFile << "// ROOT version " << gROOT->GetVersion()
1990  << std::endl << "//" << std::endl;
1991  // General information on the code
1992  outFile << "// This file contains the function " << std::endl
1993  << "//" << std::endl
1994  << "// double " << prefix << "MDF(double *x); " << std::endl
1995  << "//" << std::endl
1996  << "// For evaluating the parameterization obtained" << std::endl
1997  << "// from TMultiDimFit and the point x" << std::endl
1998  << "// " << std::endl
1999  << "// See TMultiDimFit class documentation for more "
2000  << "information " << std::endl << "// " << std::endl;
2001  // Header files
2002  if (isMethod)
2003  // If these are methods, we need the class header
2004  outFile << "#include \"" << classname << ".h\"" << std::endl;
2005 
2006  //
2007  // Now for the data
2008  //
2009  outFile << "//" << std::endl
2010  << "// Static data variables" << std::endl
2011  << "//" << std::endl;
2012  outFile << cv_qual << "int " << prefix << "gNVariables = "
2013  << fNVariables << ";" << std::endl;
2014  outFile << cv_qual << "int " << prefix << "gNCoefficients = "
2015  << fNCoefficients << ";" << std::endl;
2016  outFile << cv_qual << "double " << prefix << "gDMean = "
2017  << fMeanQuantity << ";" << std::endl;
2018 
2019  // Assignment to mean vector.
2020  outFile << "// Assignment to mean vector." << std::endl;
2021  outFile << cv_qual << "double " << prefix
2022  << "gXMean[] = {" << std::endl;
2023  for (i = 0; i < fNVariables; i++)
2024  outFile << (i != 0 ? ", " : " ") << fMeanVariables(i) << std::flush;
2025  outFile << " };" << std::endl << std::endl;
2026 
2027  // Assignment to minimum vector.
2028  outFile << "// Assignment to minimum vector." << std::endl;
2029  outFile << cv_qual << "double " << prefix
2030  << "gXMin[] = {" << std::endl;
2031  for (i = 0; i < fNVariables; i++)
2032  outFile << (i != 0 ? ", " : " ") << fMinVariables(i) << std::flush;
2033  outFile << " };" << std::endl << std::endl;
2034 
2035  // Assignment to maximum vector.
2036  outFile << "// Assignment to maximum vector." << std::endl;
2037  outFile << cv_qual << "double " << prefix
2038  << "gXMax[] = {" << std::endl;
2039  for (i = 0; i < fNVariables; i++)
2040  outFile << (i != 0 ? ", " : " ") << fMaxVariables(i) << std::flush;
2041  outFile << " };" << std::endl << std::endl;
2042 
2043  // Assignment to coefficients vector.
2044  outFile << "// Assignment to coefficients vector." << std::endl;
2045  outFile << cv_qual << "double " << prefix
2046  << "gCoefficient[] = {" << std::flush;
2047  for (i = 0; i < fNCoefficients; i++)
2048  outFile << (i != 0 ? "," : "") << std::endl
2049  << " " << fCoefficients(i) << std::flush;
2050  outFile << std::endl << " };" << std::endl << std::endl;
2051 
2052  // Assignment to error coefficients vector.
2053  outFile << "// Assignment to error coefficients vector." << std::endl;
2054  outFile << cv_qual << "double " << prefix
2055  << "gCoefficientRMS[] = {" << std::flush;
2056  for (i = 0; i < fNCoefficients; i++)
2057  outFile << (i != 0 ? "," : "") << std::endl
2058  << " " << fCoefficientsRMS(i) << std::flush;
2059  outFile << std::endl << " };" << std::endl << std::endl;
2060 
2061  // Assignment to powers vector.
2062  outFile << "// Assignment to powers vector." << std::endl
2063  << "// The powers are stored row-wise, that is" << std::endl
2064  << "// p_ij = " << prefix
2065  << "gPower[i * NVariables + j];" << std::endl;
2066  outFile << cv_qual << "int " << prefix
2067  << "gPower[] = {" << std::flush;
2068  for (i = 0; i < fNCoefficients; i++) {
2069  for (j = 0; j < fNVariables; j++) {
2070  if (j != 0) outFile << std::flush << " ";
2071  else outFile << std::endl << " ";
2072  outFile << fPowers[fPowerIndex[i] * fNVariables + j]
2073  << (i == fNCoefficients - 1 && j == fNVariables - 1 ? "" : ",")
2074  << std::flush;
2075  }
2076  }
2077  outFile << std::endl << "};" << std::endl << std::endl;
2078 
2079 
2080  //
2081  // Finally we reach the function itself
2082  //
2083  outFile << "// " << std::endl
2084  << "// The "
2085  << (isMethod ? "method " : "function ")
2086  << " double " << prefix
2087  << "MDF(double *x)"
2088  << std::endl << "// " << std::endl;
2089  outFile << "double " << prefix
2090  << "MDF(double *x) {" << std::endl
2091  << " double returnValue = " << prefix << "gDMean;" << std::endl
2092  << " int i = 0, j = 0, k = 0;" << std::endl
2093  << " for (i = 0; i < " << prefix << "gNCoefficients ; i++) {"
2094  << std::endl
2095  << " // Evaluate the ith term in the expansion" << std::endl
2096  << " double term = " << prefix << "gCoefficient[i];"
2097  << std::endl
2098  << " for (j = 0; j < " << prefix << "gNVariables; j++) {"
2099  << std::endl
2100  << " // Evaluate the polynomial in the jth variable." << std::endl
2101  << " int power = "<< prefix << "gPower["
2102  << prefix << "gNVariables * i + j]; " << std::endl
2103  << " double p1 = 1, p2 = 0, p3 = 0, r = 0;" << std::endl
2104  << " double v = 1 + 2. / ("
2105  << prefix << "gXMax[j] - " << prefix
2106  << "gXMin[j]) * (x[j] - " << prefix << "gXMax[j]);" << std::endl
2107  << " // what is the power to use!" << std::endl
2108  << " switch(power) {" << std::endl
2109  << " case 1: r = 1; break; " << std::endl
2110  << " case 2: r = v; break; " << std::endl
2111  << " default: " << std::endl
2112  << " p2 = v; " << std::endl
2113  << " for (k = 3; k <= power; k++) { " << std::endl
2114  << " p3 = p2 * v;" << std::endl;
2115  if (fPolyType == kLegendre)
2116  outFile << " p3 = ((2 * i - 3) * p2 * v - (i - 2) * p1)"
2117  << " / (i - 1);" << std::endl;
2118  if (fPolyType == kChebyshev)
2119  outFile << " p3 = 2 * v * p2 - p1; " << std::endl;
2120  outFile << " p1 = p2; p2 = p3; " << std::endl << " }" << std::endl
2121  << " r = p3;" << std::endl << " }" << std::endl
2122  << " // multiply this term by the poly in the jth var" << std::endl
2123  << " term *= r; " << std::endl << " }" << std::endl
2124  << " // Add this term to the final result" << std::endl
2125  << " returnValue += term;" << std::endl << " }" << std::endl
2126  << " return returnValue;" << std::endl << "}" << std::endl << std::endl;
2127 
2128  // EOF
2129  outFile << "// EOF for " << filename << std::endl;
2130 
2131  // Close the file
2132  outFile.close();
2133 
2134  if (fIsVerbose)
2135  std::cout << "done" << std::endl;
2136 }
2137 
2138 
2139 ////////////////////////////////////////////////////////////////////////////////
2140 /// Print statistics etc.
2141 /// Options are
2142 /// P Parameters
2143 /// S Statistics
2144 /// C Coefficients
2145 /// R Result of parameterisation
2146 /// F Result of fit
2147 /// K Correlation Matrix
2148 /// M Pretty print formula
2149 ///
2150 
2151 void TMultiDimFit::Print(Option_t *option) const
2152 {
2153  Int_t i = 0;
2154  Int_t j = 0;
2155 
2156  TString opt(option);
2157  opt.ToLower();
2158 
2159  if (opt.Contains("p")) {
2160  // Print basic parameters for this object
2161  std::cout << "User parameters:" << std::endl
2162  << "----------------" << std::endl
2163  << " Variables: " << fNVariables << std::endl
2164  << " Data points: " << fSampleSize << std::endl
2165  << " Max Terms: " << fMaxTerms << std::endl
2166  << " Power Limit Parameter: " << fPowerLimit << std::endl
2167  << " Max functions: " << fMaxFunctions << std::endl
2168  << " Max functions to study: " << fMaxStudy << std::endl
2169  << " Max angle (optional): " << fMaxAngle << std::endl
2170  << " Min angle: " << fMinAngle << std::endl
2171  << " Relative Error accepted: " << fMinRelativeError << std::endl
2172  << " Maximum Powers: " << std::flush;
2173  for (i = 0; i < fNVariables; i++)
2174  std::cout << " " << fMaxPowers[i] - 1 << std::flush;
2175  std::cout << std::endl << std::endl
2176  << " Parameterisation will be done using " << std::flush;
2177  if (fPolyType == kChebyshev)
2178  std::cout << "Chebyshev polynomials" << std::endl;
2179  else if (fPolyType == kLegendre)
2180  std::cout << "Legendre polynomials" << std::endl;
2181  else
2182  std::cout << "Monomials" << std::endl;
2183  std::cout << std::endl;
2184  }
2185 
2186  if (opt.Contains("s")) {
2187  // Print statistics for read data
2188  std::cout << "Sample statistics:" << std::endl
2189  << "------------------" << std::endl
2190  << " D" << std::flush;
2191  for (i = 0; i < fNVariables; i++)
2192  std::cout << " " << std::setw(10) << i+1 << std::flush;
2193  std::cout << std::endl << " Max: " << std::setw(10) << std::setprecision(7)
2194  << fMaxQuantity << std::flush;
2195  for (i = 0; i < fNVariables; i++)
2196  std::cout << " " << std::setw(10) << std::setprecision(4)
2197  << fMaxVariables(i) << std::flush;
2198  std::cout << std::endl << " Min: " << std::setw(10) << std::setprecision(7)
2199  << fMinQuantity << std::flush;
2200  for (i = 0; i < fNVariables; i++)
2201  std::cout << " " << std::setw(10) << std::setprecision(4)
2202  << fMinVariables(i) << std::flush;
2203  std::cout << std::endl << " Mean: " << std::setw(10) << std::setprecision(7)
2204  << fMeanQuantity << std::flush;
2205  for (i = 0; i < fNVariables; i++)
2206  std::cout << " " << std::setw(10) << std::setprecision(4)
2207  << fMeanVariables(i) << std::flush;
2208  std::cout << std::endl << " Function Sum Squares: " << fSumSqQuantity
2209  << std::endl << std::endl;
2210  }
2211 
2212  if (opt.Contains("r")) {
2213  std::cout << "Results of Parameterisation:" << std::endl
2214  << "----------------------------" << std::endl
2215  << " Total reduction of square residuals "
2216  << fSumSqResidual << std::endl
2217  << " Relative precision obtained: "
2218  << fPrecision << std::endl
2219  << " Error obtained: "
2220  << fError << std::endl
2221  << " Multiple correlation coefficient: "
2222  << fCorrelationCoeff << std::endl
2223  << " Reduced Chi square over sample: "
2224  << fChi2 / (fSampleSize - fNCoefficients) << std::endl
2225  << " Maximum residual value: "
2226  << fMaxResidual << std::endl
2227  << " Minimum residual value: "
2228  << fMinResidual << std::endl
2229  << " Estimated root mean square: "
2230  << fRMS << std::endl
2231  << " Maximum powers used: " << std::flush;
2232  for (j = 0; j < fNVariables; j++)
2233  std::cout << fMaxPowersFinal[j] << " " << std::flush;
2234  std::cout << std::endl
2235  << " Function codes of candidate functions." << std::endl
2236  << " 1: considered,"
2237  << " 2: too little contribution,"
2238  << " 3: accepted." << std::flush;
2239  for (i = 0; i < fMaxFunctions; i++) {
2240  if (i % 60 == 0)
2241  std::cout << std::endl << " " << std::flush;
2242  else if (i % 10 == 0)
2243  std::cout << " " << std::flush;
2244  std::cout << fFunctionCodes[i];
2245  }
2246  std::cout << std::endl << " Loop over candidates stopped because " << std::flush;
2247  switch(fParameterisationCode){
2248  case PARAM_MAXSTUDY:
2249  std::cout << "max allowed studies reached" << std::endl; break;
2250  case PARAM_SEVERAL:
2251  std::cout << "all candidates considered several times" << std::endl; break;
2252  case PARAM_RELERR:
2253  std::cout << "wanted relative error obtained" << std::endl; break;
2254  case PARAM_MAXTERMS:
2255  std::cout << "max number of terms reached" << std::endl; break;
2256  default:
2257  std::cout << "some unknown reason" << std::endl;
2258  break;
2259  }
2260  std::cout << std::endl;
2261  }
2262 
2263  if (opt.Contains("f")) {
2264  std::cout << "Results of Fit:" << std::endl
2265  << "---------------" << std::endl
2266  << " Test sample size: "
2267  << fTestSampleSize << std::endl
2268  << " Multiple correlation coefficient: "
2269  << fTestCorrelationCoeff << std::endl
2270  << " Relative precision obtained: "
2271  << fTestPrecision << std::endl
2272  << " Error obtained: "
2273  << fTestError << std::endl
2274  << " Reduced Chi square over sample: "
2275  << fChi2 / (fSampleSize - fNCoefficients) << std::endl
2276  << std::endl;
2277  if (fFitter) {
2278  fFitter->PrintResults(1,1);
2279  std::cout << std::endl;
2280  }
2281  }
2282 
2283  if (opt.Contains("c")){
2284  std::cout << "Coefficients:" << std::endl
2285  << "-------------" << std::endl
2286  << " # Value Error Powers" << std::endl
2287  << " ---------------------------------------" << std::endl;
2288  for (i = 0; i < fNCoefficients; i++) {
2289  std::cout << " " << std::setw(3) << i << " "
2290  << std::setw(12) << fCoefficients(i) << " "
2291  << std::setw(12) << fCoefficientsRMS(i) << " " << std::flush;
2292  for (j = 0; j < fNVariables; j++)
2293  std::cout << " " << std::setw(3)
2294  << fPowers[fPowerIndex[i] * fNVariables + j] - 1 << std::flush;
2295  std::cout << std::endl;
2296  }
2297  std::cout << std::endl;
2298  }
2299  if (opt.Contains("k") && fCorrelationMatrix.IsValid()) {
2300  std::cout << "Correlation Matrix:" << std::endl
2301  << "-------------------";
2303  }
2304 
2305  if (opt.Contains("m")) {
2306  std::cout << "Parameterization:" << std::endl
2307  << "-----------------" << std::endl
2308  << " Normalised variables: " << std::endl;
2309  for (i = 0; i < fNVariables; i++)
2310  std::cout << "\ty_" << i << "\t= 1 + 2 * (x_" << i << " - "
2311  << fMaxVariables(i) << ") / ("
2312  << fMaxVariables(i) << " - " << fMinVariables(i) << ")"
2313  << std::endl;
2314  std::cout << std::endl
2315  << " f(";
2316  for (i = 0; i < fNVariables; i++) {
2317  std::cout << "y_" << i;
2318  if (i != fNVariables-1) std::cout << ", ";
2319  }
2320  std::cout << ") = ";
2321  for (i = 0; i < fNCoefficients; i++) {
2322  if (i != 0)
2323  std::cout << std::endl << "\t" << (fCoefficients(i) < 0 ? "- " : "+ ")
2324  << TMath::Abs(fCoefficients(i));
2325  else
2326  std::cout << fCoefficients(i);
2327  for (j = 0; j < fNVariables; j++) {
2328  Int_t p = fPowers[fPowerIndex[i] * fNVariables + j];
2329  switch (p) {
2330  case 1: break;
2331  case 2: std::cout << " * y_" << j; break;
2332  default:
2333  switch(fPolyType) {
2334  case kLegendre: std::cout << " * L_" << p-1 << "(y_" << j << ")"; break;
2335  case kChebyshev: std::cout << " * C_" << p-1 << "(y_" << j << ")"; break;
2336  default: std::cout << " * y_" << j << "^" << p-1; break;
2337  }
2338  }
2339 
2340  }
2341  }
2342  std::cout << std::endl;
2343  }
2344 }
2345 
2346 
2347 ////////////////////////////////////////////////////////////////////////////////
2348 /// Selection method. User can override this method for specialized
2349 /// selection of acceptable functions in fit. Default is to select
2350 /// all. This message is sent during the build-up of the function
2351 /// candidates table once for each set of powers in
2352 /// variables. Notice, that the argument array contains the powers
2353 /// PLUS ONE. For example, to De select the function
2354 /// f = x1^2 * x2^4 * x3^5,
2355 /// this method should return kFALSE if given the argument
2356 /// { 3, 4, 6 }
2357 
2359 {
2360  return kTRUE;
2361 }
2362 
2363 ////////////////////////////////////////////////////////////////////////////////
2364 /// Set the max angle (in degrees) between the initial data vector to
2365 /// be fitted, and the new candidate function to be included in the
2366 /// fit. By default it is 0, which automatically chooses another
2367 /// selection criteria. See also
2368 /// Begin_Html<a href="#TMultiDimFit:description">class description</a>End_Html
2369 
2371 {
2372  if (ang >= 90 || ang < 0) {
2373  Warning("SetMaxAngle", "angle must be in [0,90)");
2374  return;
2375  }
2376 
2377  fMaxAngle = ang;
2378 }
2379 
2380 ////////////////////////////////////////////////////////////////////////////////
2381 /// Set the min angle (in degrees) between a new candidate function
2382 /// and the subspace spanned by the previously accepted
2383 /// functions. See also
2384 /// Begin_Html<a href="#TMultiDimFit:description">class description</a>End_Html
2385 
2387 {
2388  if (ang > 90 || ang <= 0) {
2389  Warning("SetMinAngle", "angle must be in [0,90)");
2390  return;
2391  }
2392 
2393  fMinAngle = ang;
2394 
2395 }
2396 
2397 
2398 ////////////////////////////////////////////////////////////////////////////////
2399 /// Define a user function. The input array must be of the form
2400 /// (p11, ..., p1N, ... ,pL1, ..., pLN)
2401 /// Where N is the dimension of the data sample, L is the number of
2402 /// terms (given in terms) and the first number, labels the term, the
2403 /// second the variable. More information is given in the
2404 /// Begin_Html<a href="#TMultiDimFit:description">class description</a>End_Html
2405 
2406 void TMultiDimFit::SetPowers(const Int_t* powers, Int_t terms)
2407 {
2409  fMaxFunctions = terms;
2410  fMaxTerms = terms;
2411  fMaxStudy = terms;
2413  fPowers = new Int_t[fMaxFuncNV];
2414  Int_t i, j;
2415  for (i = 0; i < fMaxFunctions; i++)
2416  for(j = 0; j < fNVariables; j++)
2417  fPowers[i * fNVariables + j] = powers[i * fNVariables + j] + 1;
2418 }
2419 
2420 ////////////////////////////////////////////////////////////////////////////////
2421 /// Set the user parameter for the function selection. The bigger the
2422 /// limit, the more functions are used. The meaning of this variable
2423 /// is defined in the
2424 /// Begin_Html<a href="#TMultiDimFit:description">class description</a>End_Html
2425 
2427 {
2428  fPowerLimit = limit;
2429 }
2430 
2431 ////////////////////////////////////////////////////////////////////////////////
2432 /// Set the maximum power to be considered in the fit for each
2433 /// variable. See also
2434 /// Begin_Html<a href="#TMultiDimFit:description">class description</a>End_Html
2435 
2437 {
2438  if (!powers)
2439  return;
2440 
2441  for (Int_t i = 0; i < fNVariables; i++)
2442  fMaxPowers[i] = powers[i]+1;
2443 }
2444 
2445 ////////////////////////////////////////////////////////////////////////////////
2446 /// Set the acceptable relative error for when sum of square
2447 /// residuals is considered minimized. For a full account, refer to
2448 /// the
2449 /// Begin_Html<a href="#TMultiDimFit:description">class description</a>End_Html
2450 
2452 {
2453  fMinRelativeError = error;
2454 }
2455 
2456 
2457 ////////////////////////////////////////////////////////////////////////////////
2458 /// PRIVATE METHOD:
2459 /// Test whether the currently considered function contributes to the
2460 /// fit. See also
2461 /// Begin_Html<a href="#TMultiDimFit:description">class description</a>End_Html
2462 
2464  Double_t dResidur)
2465 {
2466  if (fNCoefficients != 0) {
2467  // Now for the second test:
2468  if (fMaxAngle == 0) {
2469  // If the user hasn't supplied a max angle do the test as,
2470  if (dResidur <
2471  squareResidual / (fMaxTerms - fNCoefficients + 1 + 1E-10)) {
2472  return kFALSE;
2473  }
2474  }
2475  else {
2476  // If the user has provided a max angle, test if the calculated
2477  // angle is less then the max angle.
2478  if (TMath::Sqrt(dResidur/fSumSqAvgQuantity) <
2480  return kFALSE;
2481  }
2482  }
2483  }
2484  // If we get here, the function is OK
2485  return kTRUE;
2486 }
2487 
2488 
2489 ////////////////////////////////////////////////////////////////////////////////
2490 /// Helper function for doing the minimisation of Chi2 using Minuit
2491 
2492 void mdfHelper(int& /*npar*/, double* /*divs*/, double& chi2,
2493  double* coeffs, int /*flag*/)
2494 {
2495  // Get pointer to current TMultiDimFit object.
2497  chi2 = mdf->MakeChi2(coeffs);
2498 }
void Add(TObject *obj, const char *name=0, Int_t check=-1)
Add object with name to browser.
Definition: TBrowser.cxx:259
Int_t * fMaxPowersFinal
Definition: TMultiDimFit.h:78
virtual Double_t MakeChi2(const Double_t *coeff=0)
Calculate Chi square over either the test sample.
TVirtualFitter * fFitter
Definition: TMultiDimFit.h:111
TVectorT< Element > & ResizeTo(Int_t lwb, Int_t upb)
Resize the vector to [lwb:upb] .
Definition: TVectorT.cxx:291
static void mdfHelper(int &, double *, double &, double *, int)
Helper function for doing the minimisation of Chi2 using Minuit.
#define PARAM_RELERR
Double_t fError
Definition: TMultiDimFit.h:98
TVectorD fResiduals
Definition: TMultiDimFit.h:82
Int_t * fPowerIndex
Definition: TMultiDimFit.h:80
static double p3(double t, double a, double b, double c, double d)
TMatrixTDiag_const< Double_t > TMatrixDDiag_const
#define HIST_XORIG
static Vc_ALWAYS_INLINE float_v trunc(float_v::AsArg v)
Definition: math.h:95
ClassImp(TSeqCollection) Int_t TSeqCollection TIter next(this)
Return index of object in collection.
Ssiz_t Length() const
Definition: TString.h:390
virtual TMatrixTBase< Element > & Zero()
Set matrix elements to zero.
Bool_t IsValid() const
Definition: TVectorT.h:89
Double_t fPrecision
Definition: TMultiDimFit.h:100
#define HIST_DORIG
const char Option_t
Definition: RtypesCore.h:62
static TMultiDimFit * fgInstance
Definition: TMultiDimFit.h:35
virtual Int_t SetParameter(Int_t ipar, const char *parname, Double_t value, Double_t verr, Double_t vlow, Double_t vhigh)=0
Bool_t fShowCorrelation
Definition: TMultiDimFit.h:114
virtual void Fit(Option_t *option="")
Try to fit the found parameterisation to the test sample.
virtual ~TMultiDimFit()
Destructor.
#define PARAM_MAXSTUDY
Int_t fMaxFuncNV
Definition: TMultiDimFit.h:72
Int_t GetNrows() const
Definition: TVectorT.h:81
Int_t fTestSampleSize
Definition: TMultiDimFit.h:58
#define DEGRAD
TH1 * h
Definition: legend2.C:5
TVectorD fCoefficients
Definition: TMultiDimFit.h:92
virtual void Print(Option_t *option="ps") const
Print statistics etc.
Double_t fMinRelativeError
Definition: TMultiDimFit.h:63
void SetMaxPowers(const Int_t *powers)
Set the maximum power to be considered in the fit for each variable.
virtual Double_t GetParameter(Int_t ipar) const =0
TMatrixD fOrthFunctions
Definition: TMultiDimFit.h:74
TMatrixD fOrthCurvatureMatrix
Definition: TMultiDimFit.h:91
static const char * filename()
#define gROOT
Definition: TROOT.h:340
Double_t fTestError
Definition: TMultiDimFit.h:99
Basic string class.
Definition: TString.h:137
Bool_t fIsUserFunction
Definition: TMultiDimFit.h:115
void ToLower()
Change string to lower-case.
Definition: TString.cxx:1088
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
const Bool_t kFALSE
Definition: Rtypes.h:92
virtual TObject * FindObject(const char *name) const
Find an object in this list using its name.
Definition: TList.cxx:496
virtual void SetFCN(void *fcn)
To set the address of the minimization objective function.
virtual Double_t EvalFactor(Int_t p, Double_t x) const
PRIVATE METHOD: Evaluate function with power p at variable value x.
virtual void Clear(Option_t *option="")
Clear internal structures and variables.
Double_t fMaxAngle
Definition: TMultiDimFit.h:61
virtual void MakeNormalized()
PRIVATE METHOD: Normalize data to the interval [-1;1].
virtual TMatrixTBase< Element > & ResizeTo(Int_t nrows, Int_t ncols, Int_t=-1)
Set size of the matrix to nrows x ncols New dynamic elements are created, the overlapping part of the...
Definition: TMatrixT.cxx:1201
virtual void MakeRealCode(const char *filename, const char *classname, Option_t *option="")
PRIVATE METHOD: This is the method that actually generates the code for the evaluation the parameteri...
virtual Bool_t TestFunction(Double_t squareResidual, Double_t dResidur)
PRIVATE METHOD: Test whether the currently considered function contributes to the fit...
Short_t Abs(Short_t d)
Definition: TMathBase.h:110
virtual void MakeParameterization()
PRIVATE METHOD: Find the parameterization over the training sample.
TVectorD fTestVariables
Definition: TMultiDimFit.h:56
Double_t fSumSqQuantity
Definition: TMultiDimFit.h:43
virtual void MakeCoefficientErrors()
PRIVATE METHOD: Compute the errors on the coefficients.
Double_t fMinAngle
Definition: TMultiDimFit.h:60
virtual Double_t Eval(const Double_t *x, const Double_t *coeff=0) const
Evaluate parameterization at point x.
Bool_t IsValid() const
Definition: TMatrixTBase.h:157
const char * Data() const
Definition: TString.h:349
TVectorD fTestQuantity
Definition: TMultiDimFit.h:54
virtual void PrintResults(Int_t level, Double_t amin) const =0
double sqrt(double)
TVectorD fTestSqError
Definition: TMultiDimFit.h:55
Int_t * fFunctionCodes
Definition: TMultiDimFit.h:70
Double_t x[n]
Definition: legend1.C:17
void SetPowerLimit(Double_t limit=1e-3)
Set the user parameter for the function selection.
Double_t fSumSqAvgQuantity
Definition: TMultiDimFit.h:44
#define HIST_RD
virtual Double_t EvalControl(const Int_t *powers) const
PRIVATE METHOD: Calculate the control parameter from the passed powers.
The TNamed class is the base class for all named ROOT classes.
Definition: TNamed.h:33
static TMultiDimFit * Instance()
Return the static instance.
TVectorD fSqError
Definition: TMultiDimFit.h:39
TMatrixD fCorrelationMatrix
Definition: TMultiDimFit.h:103
static double p2(double t, double a, double b, double c)
TVectorD fOrthFunctionNorms
Definition: TMultiDimFit.h:75
void SetMinAngle(Double_t angle=1)
Set the min angle (in degrees) between a new candidate function and the subspace spanned by the previ...
Double_t fMaxQuantity
Definition: TMultiDimFit.h:41
void SetMaxAngle(Double_t angle=0)
Set the max angle (in degrees) between the initial data vector to be fitted, and the new candidate fu...
Byte_t fHistogramMask
Definition: TMultiDimFit.h:107
TList * fHistograms
Definition: TMultiDimFit.h:106
virtual void AddRow(const Double_t *x, Double_t D, Double_t E=0)
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.
Bool_t Invert(TMatrixDSym &inv)
For a symmetric matrix A(m,m), its inverse A_inv(m,m) is returned .
virtual Bool_t Select(const Int_t *iv)
Selection method.
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:918
Int_t fNVariables
Definition: TMultiDimFit.h:47
TVectorD fQuantity
Definition: TMultiDimFit.h:38
TVectorD fMeanVariables
Definition: TMultiDimFit.h:48
char * out
Definition: TBase64.cxx:29
TMatrixTRow< Double_t > TMatrixDRow
Double_t fTestPrecision
Definition: TMultiDimFit.h:101
A doubly linked list.
Definition: TList.h:47
Double_t fMeanQuantity
Definition: TMultiDimFit.h:40
void Print(Option_t *name="") const
Print the matrix as a table of elements.
Double_t fMaxResidual
Definition: TMultiDimFit.h:83
Int_t * fPowers
Definition: TMultiDimFit.h:79
Int_t * fMaxPowers
Definition: TMultiDimFit.h:64
Using a TBrowser one can browse all ROOT objects.
Definition: TBrowser.h:41
virtual TMatrixTBase< Element > & NormByDiag(const TVectorT< Element > &v, Option_t *option="D")
option: "D" : b(i,j) = a(i,j)/sqrt(abs*(v(i)*v(j))) (default) else : b(i,j) = a(i,j)*sqrt(abs*(v(i)*v(j))) (default)
Int_t fMaxFunctions
Definition: TMultiDimFit.h:69
ROOT::R::TRInterface & r
Definition: Object.C:4
Bool_t EndsWith(const char *pat, ECaseCompare cmp=kExact) const
Return true if string ends with the specified string.
Definition: TString.cxx:2220
virtual void MakeCandidates()
PRIVATE METHOD: Create list of candidate functions for the parameterisation.
Int_t fSampleSize
Definition: TMultiDimFit.h:52
Double_t fSumSqResidual
Definition: TMultiDimFit.h:87
#define HIST_DSHIF
Int_t fMaxTerms
Definition: TMultiDimFit.h:62
Double_t fChi2
Definition: TMultiDimFit.h:95
TVectorT< Element > & Zero()
Set vector elements to zero.
Definition: TVectorT.cxx:450
Multidimensional Fits in ROOT.
Definition: TMultiDimFit.h:25
TMarker * m
Definition: textangle.C:8
char * Form(const char *fmt,...)
Int_t fParameterisationCode
Definition: TMultiDimFit.h:96
Double_t E()
Definition: TMath.h:54
virtual void FindParameterization(Option_t *option="")
Find the parameterization.
TLine * l
Definition: textangle.C:4
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:51
virtual Int_t ExecuteCommand(const char *command, Double_t *args, Int_t nargs)=0
Bool_t fIsVerbose
Definition: TMultiDimFit.h:116
Double_t fMinQuantity
Definition: TMultiDimFit.h:42
TVectorD fVariables
Definition: TMultiDimFit.h:46
virtual void SetPowers(const Int_t *powers, Int_t terms)
Define a user function.
#define HIST_RX
static double p1(double t, double a, double b)
virtual void AddTestRow(const Double_t *x, Double_t D, Double_t E=0)
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.
1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:613
TMatrixD fFunctions
Definition: TMultiDimFit.h:68
virtual void MakeMethod(const Char_t *className="MDF", Option_t *option="")
Generate the file MDF.cxx which contains the implementation of the method: ...
Double_t fMinResidual
Definition: TMultiDimFit.h:84
REAL epsilon
Definition: triangle.c:617
Double_t Cos(Double_t)
Definition: TMath.h:424
virtual void MakeHistograms(Option_t *option="A")
Make histograms of the result of the analysis.
#define HIST_RTRAI
EMDFPolyType fPolyType
Fit object (MINUIT)
Definition: TMultiDimFit.h:113
double f(double x)
double Double_t
Definition: RtypesCore.h:55
int type
Definition: TGX11.cxx:120
void SetMinRelativeError(Double_t error)
Set the acceptable relative error for when sum of square residuals is considered minimized.
Int_t fMinResidualRow
Definition: TMultiDimFit.h:86
Double_t y[n]
Definition: legend1.C:17
TVectorD fCoefficientsRMS
Definition: TMultiDimFit.h:93
The TH1 histogram class.
Definition: TH1.h:80
Int_t fNCoefficients
Definition: TMultiDimFit.h:89
virtual void MakeCode(const char *functionName="MDF", Option_t *option="")
Generate the file with .C appended if argument doesn't end in .cxx or .C.
TMultiDimFit()
Empty CTOR. Do not use.
virtual void Clear(Option_t *option="")
Remove all objects from the list.
Definition: TList.cxx:348
Int_t fMaxStudy
Definition: TMultiDimFit.h:71
virtual Bool_t Decompose()
Matrix A is decomposed in component U so that A = U^T * U If the decomposition succeeds, bit kDecomposed is set , otherwise kSingular.
virtual Double_t EvalError(const Double_t *x, const Double_t *coeff=0) const
Evaluate parameterization error at point x.
Double_t fPowerLimit
Definition: TMultiDimFit.h:65
#define HIST_XNORM
#define SETBIT(n, i)
Definition: Rtypes.h:121
Int_t fMaxResidualRow
Definition: TMultiDimFit.h:85
#define HIST_RTEST
char Char_t
Definition: RtypesCore.h:29
static TVirtualFitter * Fitter(TObject *obj, Int_t maxpar=25)
Static function returning a pointer to the current fitter.
virtual void Add(TObject *obj)
Definition: TList.h:81
double f2(const double *x)
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:202
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:567
Double_t fRMS
Definition: TMultiDimFit.h:94
Double_t Sin(Double_t)
Definition: TMath.h:421
#define TESTBIT(n, i)
Definition: Rtypes.h:123
virtual void MakeCorrelation()
PRIVATE METHOD: Compute the correlation matrix.
Double_t Sqrt(Double_t x)
Definition: TMath.h:464
virtual Double_t MakeGramSchmidt(Int_t function)
PRIVATE METHOD: Make Gram-Schmidt orthogonalisation.
Double_t fCorrelationCoeff
Definition: TMultiDimFit.h:102
const char * AsString() const
Return the date & time as a string (ctime() format).
Definition: TDatime.cxx:99
const Bool_t kTRUE
Definition: Rtypes.h:91
TVectorD fMinVariables
Definition: TMultiDimFit.h:50
TVectorD fOrthCoefficients
Definition: TMultiDimFit.h:90
const Int_t n
Definition: legend1.C:16
virtual void MakeCoefficients()
PRIVATE METHOD: Invert the model matrix B, and compute final coefficients.
#define PARAM_MAXTERMS
#define PARAM_SEVERAL
ClassImp(TMultiDimFit)
This class stores the date and time with a precision of one second in an unsigned 32 bit word (950130...
Definition: TDatime.h:39
TVectorD fMaxVariables
Definition: TMultiDimFit.h:49
virtual void Browse(TBrowser *b)
Browse the TMultiDimFit object in the TBrowser.
2-D histogram with a double per channel (see TH1 documentation)}
Definition: TH2.h:297
Double_t fTestCorrelationCoeff
Definition: TMultiDimFit.h:104
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:904