Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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 representative
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 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 representative 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 multidimensional fit,
326 using \f$\mathbf{x}^\prime\f$ as the variables, thus determining a set of
327 coefficients \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. Determine \f$\mathbf{x}^\prime\f$ as before.
335 4. Use the result of the fit to determine \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 representative 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 University 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 "TList.h"
396#include "TBrowser.h"
397#include "TDecompChol.h"
398#include "TDatime.h"
399
400
401#define RADDEG (180. / TMath::Pi())
402#define DEGRAD (TMath::Pi() / 180.)
403#define HIST_XORIG 0
404#define HIST_DORIG 1
405#define HIST_XNORM 2
406#define HIST_DSHIF 3
407#define HIST_RX 4
408#define HIST_RD 5
409#define HIST_RTRAI 6
410#define HIST_RTEST 7
411#define PARAM_MAXSTUDY 1
412#define PARAM_SEVERAL 2
413#define PARAM_RELERR 3
414#define PARAM_MAXTERMS 4
415
416
417////////////////////////////////////////////////////////////////////////////////
418
419static void mdfHelper(int&, double*, double&, double*, int);
420
421////////////////////////////////////////////////////////////////////////////////
422
424
425//____________________________________________________________________
426// Static instance. Used with mdfHelper and TMinuit
428
429
430////////////////////////////////////////////////////////////////////////////////
431/// Empty CTOR. Do not use
432
434{
435 fMeanQuantity = 0;
436 fMaxQuantity = 0;
437 fMinQuantity = 0;
438 fSumSqQuantity = 0;
440
441 fNVariables = 0;
442 fSampleSize = 0;
443 fTestSampleSize = 0;
444
445 fMinAngle = 1;
446 fMaxAngle = 0;
447 fMaxTerms = 0;
449 fMaxPowers = 0;
450 fPowerLimit = 0;
451
452 fMaxFunctions = 0;
453 fFunctionCodes = 0;
454 fMaxStudy = 0;
455 fMaxFuncNV = 0;
456
457 fMaxPowersFinal = 0;
458 fPowers = 0;
459 fPowerIndex = 0;
460
461 fMaxResidual = 0;
462 fMinResidual = 0;
463 fMaxResidualRow = 0;
464 fMinResidualRow = 0;
465 fSumSqResidual = 0;
466
467 fNCoefficients = 0;
468 fRMS = 0;
469 fChi2 = 0;
471
472 fError = 0;
473 fTestError = 0;
474 fPrecision = 0;
475 fTestPrecision = 0;
478
479 fHistograms = 0;
480 fHistogramMask = 0;
481 fBinVarX = 100;
482 fBinVarY = 100;
483
484 fFitter = 0;
489
490}
491
492
493////////////////////////////////////////////////////////////////////////////////
494/// Constructor
495/// Second argument is the type of polynomials to use in
496/// parameterisation, one of:
497/// TMultiDimFit::kMonomials
498/// TMultiDimFit::kChebyshev
499/// TMultiDimFit::kLegendre
500///
501/// Options:
502/// K Compute (k)correlation matrix
503/// V Be verbose
504///
505/// Default is no options.
506///
507
510 Option_t *option)
511: TNamed("multidimfit","Multi-dimensional fit object"),
512fQuantity(dimension),
513fSqError(dimension),
514fVariables(dimension*100),
515fMeanVariables(dimension),
516fMaxVariables(dimension),
517fMinVariables(dimension)
518{
519 fgInstance = this;
520
521 fMeanQuantity = 0;
522 fMaxQuantity = 0;
523 fMinQuantity = 0;
524 fSumSqQuantity = 0;
526
527 fNVariables = dimension;
528 fSampleSize = 0;
529 fTestSampleSize = 0;
530
531 fMinAngle = 1;
532 fMaxAngle = 0;
533 fMaxTerms = 0;
534 fMinRelativeError = 0.01;
535 fMaxPowers = new Int_t[dimension];
536 fPowerLimit = 1;
537
538 fMaxFunctions = 0;
539 fFunctionCodes = 0;
540 fMaxStudy = 0;
541 fMaxFuncNV = 0;
542
543 fMaxPowersFinal = new Int_t[dimension];
544 fPowers = 0;
545 fPowerIndex = 0;
546
547 fMaxResidual = 0;
548 fMinResidual = 0;
549 fMaxResidualRow = 0;
550 fMinResidualRow = 0;
551 fSumSqResidual = 0;
552
553 fNCoefficients = 0;
554 fRMS = 0;
555 fChi2 = 0;
557
558 fError = 0;
559 fTestError = 0;
560 fPrecision = 0;
561 fTestPrecision = 0;
564
565 fHistograms = 0;
566 fHistogramMask = 0;
567 fBinVarX = 100;
568 fBinVarY = 100;
569
570 fFitter = 0;
571 fPolyType = type;
575 TString opt = option;
576 opt.ToLower();
577
578 if (opt.Contains("k")) fShowCorrelation = kTRUE;
579 if (opt.Contains("v")) fIsVerbose = kTRUE;
580}
581
582
583////////////////////////////////////////////////////////////////////////////////
584/// Destructor
585
587{
588 delete [] fPowers;
589 delete [] fMaxPowers;
590 delete [] fMaxPowersFinal;
591 delete [] fPowerIndex;
592 delete [] fFunctionCodes;
593 if (fHistograms) fHistograms->Clear("nodelete");
594 delete fHistograms;
595}
596
597
598////////////////////////////////////////////////////////////////////////////////
599/// Add a row consisting of fNVariables independent variables, the
600/// known, dependent quantity, and optionally, the square error in
601/// the dependent quantity, to the training sample to be used for the
602/// parameterization.
603/// The mean of the variables and quantity is calculated on the fly,
604/// as outlined in TPrincipal::AddRow.
605/// This sample should be representative of the problem at hand.
606/// Please note, that if no error is given Poisson statistics is
607/// assumed and the square error is set to the value of dependent
608/// quantity. See also the
609/// <a href="#TMultiDimFit:description">class description</a>
610
612{
613 if (!x)
614 return;
615
616 if (++fSampleSize == 1) {
617 fMeanQuantity = D;
618 fMaxQuantity = D;
619 fMinQuantity = D;
620 fSumSqQuantity = D * D;// G.Q. erratum on August 15th, 2008
621 }
622 else {
625 fSumSqQuantity += D * D;
626
627 if (D >= fMaxQuantity) fMaxQuantity = D;
628 if (D <= fMinQuantity) fMinQuantity = D;
629 }
630
631
632 // If the vector isn't big enough to hold the new data, then
633 // expand the vector by half it's size.
634 Int_t size = fQuantity.GetNrows();
635 if (fSampleSize > size) {
636 fQuantity.ResizeTo(size + size/2);
637 fSqError.ResizeTo(size + size/2);
638 }
639
640 // Store the value
641 fQuantity(fSampleSize-1) = D;
642 fSqError(fSampleSize-1) = (E == 0 ? D : E);
643
644 // Store data point in internal vector
645 // If the vector isn't big enough to hold the new data, then
646 // expand the vector by half it's size
647 size = fVariables.GetNrows();
648 if (fSampleSize * fNVariables > size)
649 fVariables.ResizeTo(size + size/2);
650
651
652 // Increment the data point counter
653 Int_t i,j;
654 for (i = 0; i < fNVariables; i++) {
655 if (fSampleSize == 1) {
656 fMeanVariables(i) = x[i];
657 fMaxVariables(i) = x[i];
658 fMinVariables(i) = x[i];
659 }
660 else {
661 fMeanVariables(i) *= 1 - 1./Double_t(fSampleSize);
663
664 // Update the maximum value for this component
665 if (x[i] >= fMaxVariables(i)) fMaxVariables(i) = x[i];
666
667 // Update the minimum value for this component
668 if (x[i] <= fMinVariables(i)) fMinVariables(i) = x[i];
669
670 }
671
672 // Store the data.
673 j = (fSampleSize-1) * fNVariables + i;
674 fVariables(j) = x[i];
675 }
676}
677
678
679////////////////////////////////////////////////////////////////////////////////
680/// Add a row consisting of fNVariables independent variables, the
681/// known, dependent quantity, and optionally, the square error in
682/// the dependent quantity, to the test sample to be used for the
683/// test of the parameterization.
684/// This sample needn't be representative of the problem at hand.
685/// Please note, that if no error is given Poisson statistics is
686/// assumed and the square error is set to the value of dependent
687/// quantity. See also the
688/// <a href="#TMultiDimFit:description">class description</a>
689
691{
692 if (fTestSampleSize++ == 0) {
696 }
697
698 // If the vector isn't big enough to hold the new data, then
699 // expand the vector by half it's size.
701 if (fTestSampleSize > size) {
702 fTestQuantity.ResizeTo(size + size/2);
703 fTestSqError.ResizeTo(size + size/2);
704 }
705
706 // Store the value
708 fTestSqError(fTestSampleSize-1) = (E == 0 ? D : E);
709
710 // Store data point in internal vector
711 // If the vector isn't big enough to hold the new data, then
712 // expand the vector by half it's size
713 size = fTestVariables.GetNrows();
714 if (fTestSampleSize * fNVariables > size)
715 fTestVariables.ResizeTo(size + size/2);
716
717
718 // Increment the data point counter
719 Int_t i,j;
720 for (i = 0; i < fNVariables; i++) {
721 j = fNVariables * (fTestSampleSize - 1) + i;
722 fTestVariables(j) = x[i];
723
724 if (x[i] > fMaxVariables(i))
725 Warning("AddTestRow", "variable %d (row: %d) too large: %f > %f",
726 i, fTestSampleSize, x[i], fMaxVariables(i));
727 if (x[i] < fMinVariables(i))
728 Warning("AddTestRow", "variable %d (row: %d) too small: %f < %f",
729 i, fTestSampleSize, x[i], fMinVariables(i));
730 }
731}
732
733
734////////////////////////////////////////////////////////////////////////////////
735/// Browse the TMultiDimFit object in the TBrowser.
736
738{
739 if (fHistograms) {
740 TIter next(fHistograms);
741 TH1* h = 0;
742 while ((h = (TH1*)next()))
743 b->Add(h,h->GetName());
744 }
745 if (fVariables.IsValid())
746 b->Add(&fVariables, "Variables (Training)");
747 if (fQuantity.IsValid())
748 b->Add(&fQuantity, "Quantity (Training)");
749 if (fSqError.IsValid())
750 b->Add(&fSqError, "Error (Training)");
752 b->Add(&fMeanVariables, "Mean of Variables (Training)");
754 b->Add(&fMaxVariables, "Mean of Variables (Training)");
756 b->Add(&fMinVariables, "Min of Variables (Training)");
758 b->Add(&fTestVariables, "Variables (Test)");
760 b->Add(&fTestQuantity, "Quantity (Test)");
761 if (fTestSqError.IsValid())
762 b->Add(&fTestSqError, "Error (Test)");
763 if (fFunctions.IsValid())
764 b->Add(&fFunctions, "Functions");
766 b->Add(&fCoefficients,"Coefficients");
768 b->Add(&fCoefficientsRMS,"Coefficients Errors");
770 b->Add(&fOrthFunctions, "Orthogonal Functions");
772 b->Add(&fOrthFunctionNorms, "Orthogonal Functions Norms");
773 if (fResiduals.IsValid())
774 b->Add(&fResiduals, "Residuals");
776 b->Add(&fOrthCoefficients,"Orthogonal Coefficients");
778 b->Add(&fOrthCurvatureMatrix,"Orthogonal curvature matrix");
780 b->Add(&fCorrelationMatrix,"Correlation Matrix");
781 if (fFitter)
782 b->Add(fFitter, fFitter->GetName());
783}
784
785
786////////////////////////////////////////////////////////////////////////////////
787/// Clear internal structures and variables
788
790{
791 Int_t i, j, n = fNVariables, m = fMaxFunctions;
792
793 // Training sample, dependent quantity
794 fQuantity.Zero();
795 fSqError.Zero();
796 fMeanQuantity = 0;
797 fMaxQuantity = 0;
798 fMinQuantity = 0;
799 fSumSqQuantity = 0;
801
802 // Training sample, independent variables
804 fNVariables = 0;
805 fSampleSize = 0;
809
810 // Test sample
814 fTestSampleSize = 0;
815
816 // Functions
818 //for (i = 0; i < fMaxTerms; i++) fPowerIndex[i] = 0;
819 //for (i = 0; i < fMaxTerms; i++) fFunctionCodes[i] = 0;
820 fMaxFunctions = 0;
821 fMaxStudy = 0;
824
825 // Control parameters
827 fMinAngle = 0;
828 fMaxAngle = 0;
829 fMaxTerms = 0;
830
831 // Powers
832 for (i = 0; i < n; i++) {
833 fMaxPowers[i] = 0;
834 fMaxPowersFinal[i] = 0;
835 for (j = 0; j < m; j++)
836 fPowers[i * n + j] = 0;
837 }
838 fPowerLimit = 0;
839
840 // Residuals
841 fMaxResidual = 0;
842 fMinResidual = 0;
843 fMaxResidualRow = 0;
844 fMinResidualRow = 0;
845 fSumSqResidual = 0;
846
847 // Fit
848 fNCoefficients = 0;
851 fRMS = 0;
853 fError = 0;
854 fTestError = 0;
855 fPrecision = 0;
856 fTestPrecision = 0;
857
858 // Coefficients
862 fHistograms->Clear(option);
863
864 // Options
868}
869
870
871////////////////////////////////////////////////////////////////////////////////
872/// Evaluate parameterization at point x. Optional argument coeff is
873/// a vector of coefficients for the parameterisation, fNCoefficients
874/// elements long.
875
876Double_t TMultiDimFit::Eval(const Double_t *x, const Double_t* coeff) const
877{
878 Double_t returnValue = fMeanQuantity;
879 Double_t term = 0;
880 Int_t i, j;
881
882 for (i = 0; i < fNCoefficients; i++) {
883 // Evaluate the ith term in the expansion
884 term = (coeff ? coeff[i] : fCoefficients(i));
885 for (j = 0; j < fNVariables; j++) {
886 // Evaluate the factor (polynomial) in the j-th variable.
887 Int_t p = fPowers[fPowerIndex[i] * fNVariables + j];
888 Double_t y = 1 + 2. / (fMaxVariables(j) - fMinVariables(j))
889 * (x[j] - fMaxVariables(j));
890 term *= EvalFactor(p,y);
891 }
892 // Add this term to the final result
893 returnValue += term;
894 }
895 return returnValue;
896}
897
898
899////////////////////////////////////////////////////////////////////////////////
900/// Evaluate parameterization error at point x. Optional argument coeff is
901/// a vector of coefficients for the parameterisation, fNCoefficients
902/// elements long.
903
905{
906 Double_t returnValue = 0;
907 Double_t term = 0;
908 Int_t i, j;
909
910 for (i = 0; i < fNCoefficients; i++) {
911 // std::cout << "Error coef " << i << " -> " << fCoefficientsRMS(i) << std::endl;
912 }
913 for (i = 0; i < fNCoefficients; i++) {
914 // Evaluate the ith term in the expansion
915 term = (coeff ? coeff[i] : fCoefficientsRMS(i));
916 for (j = 0; j < fNVariables; j++) {
917 // Evaluate the factor (polynomial) in the j-th variable.
918 Int_t p = fPowers[fPowerIndex[i] * fNVariables + j];
919 Double_t y = 1 + 2. / (fMaxVariables(j) - fMinVariables(j))
920 * (x[j] - fMaxVariables(j));
921 term *= EvalFactor(p,y);
922 // std::cout << "i,j " << i << ", " << j << " " << p << " " << y << " " << EvalFactor(p,y) << " " << term << std::endl;
923 }
924 // Add this term to the final result
925 returnValue += term*term;
926 // std::cout << " i = " << i << " value = " << returnValue << std::endl;
927 }
928 returnValue = sqrt(returnValue);
929 return returnValue;
930}
931
932
933////////////////////////////////////////////////////////////////////////////////
934/// PRIVATE METHOD:
935/// Calculate the control parameter from the passed powers
936
938{
939 Double_t s = 0;
940 Double_t epsilon = 1e-6; // a small number
941 for (Int_t i = 0; i < fNVariables; i++) {
942 if (fMaxPowers[i] != 1)
943 s += (epsilon + iv[i] - 1) / (epsilon + fMaxPowers[i] - 1);
944 }
945 return s;
946}
947
948////////////////////////////////////////////////////////////////////////////////
949/// PRIVATE METHOD:
950/// Evaluate function with power p at variable value x
951
953{
954 Int_t i = 0;
955 Double_t p1 = 1;
956 Double_t p2 = 0;
957 Double_t p3 = 0;
958 Double_t r = 0;
959
960 switch(p) {
961 case 1:
962 r = 1;
963 break;
964 case 2:
965 r = x;
966 break;
967 default:
968 p2 = x;
969 for (i = 3; i <= p; i++) {
970 p3 = p2 * x;
971 if (fPolyType == kLegendre)
972 p3 = ((2 * i - 3) * p2 * x - (i - 2) * p1) / (i - 1);
973 else if (fPolyType == kChebyshev)
974 p3 = 2 * x * p2 - p1;
975 p1 = p2;
976 p2 = p3;
977 }
978 r = p3;
979 }
980
981 return r;
982}
983
984
985////////////////////////////////////////////////////////////////////////////////
986/// Find the parameterization
987///
988/// Options:
989/// None so far
990///
991/// For detailed description of what this entails, please refer to the
992/// <a href="#TMultiDimFit:description">class description</a>
993
995{
1002}
1003
1004////////////////////////////////////////////////////////////////////////////////
1005/// Try to fit the found parameterisation to the test sample.
1006///
1007/// Options
1008/// M use Minuit to improve coefficients
1009///
1010/// Also, refer to
1011/// <a href="#TMultiDimFit:description">class description</a>
1012
1014{
1015 Int_t i, j;
1017 Double_t sumSqD = 0;
1018 Double_t sumD = 0;
1019 Double_t sumSqR = 0;
1020 Double_t sumR = 0;
1021
1022 // Calculate the residuals over the test sample
1023 for (i = 0; i < fTestSampleSize; i++) {
1024 for (j = 0; j < fNVariables; j++)
1025 x[j] = fTestVariables(i * fNVariables + j);
1026 Double_t res = fTestQuantity(i) - Eval(x);
1027 sumD += fTestQuantity(i);
1028 sumSqD += fTestQuantity(i) * fTestQuantity(i);
1029 sumR += res;
1030 sumSqR += res * res;
1032 ((TH1D*)fHistograms->FindObject("res_test"))->Fill(res);
1033 }
1034 Double_t dAvg = sumSqD - (sumD * sumD) / fTestSampleSize;
1035 Double_t rAvg = sumSqR - (sumR * sumR) / fTestSampleSize;
1036 fTestCorrelationCoeff = (dAvg - rAvg) / dAvg;
1037 fTestError = sumSqR;
1038 fTestPrecision = sumSqR / sumSqD;
1039
1040 TString opt(option);
1041 opt.ToLower();
1042
1043 if (!opt.Contains("m"))
1044 MakeChi2();
1045
1047 Warning("Fit", "test sample is very small");
1048
1049 if (!opt.Contains("m")) {
1050 Error("Fit", "invalid option");
1051 delete [] x;
1052 return;
1053 }
1054
1056 if (!fFitter) {
1057 Error("Fit", "Cannot create Fitter");
1058 delete [] x;
1059 return;
1060 }
1062
1063 const Int_t maxArgs = 16;
1064 Int_t args = 1;
1065 Double_t* arglist = new Double_t[maxArgs];
1066 arglist[0] = -1;
1067 fFitter->ExecuteCommand("SET PRINT",arglist,args);
1068
1069 for (i = 0; i < fNCoefficients; i++) {
1070 Double_t startVal = fCoefficients(i);
1071 Double_t startErr = fCoefficientsRMS(i);
1072 fFitter->SetParameter(i, Form("coeff%02d",i),
1073 startVal, startErr, 0, 0);
1074 }
1075
1076 // arglist[0] = 0;
1077 args = 1;
1078 // fFitter->ExecuteCommand("SET PRINT",arglist,args);
1079 fFitter->ExecuteCommand("MIGRAD",arglist,args);
1080
1081 for (i = 0; i < fNCoefficients; i++) {
1082 Double_t val = 0, err = 0, low = 0, high = 0;
1083 fFitter->GetParameter(i, Form("coeff%02d",i),
1084 val, err, low, high);
1085 fCoefficients(i) = val;
1086 fCoefficientsRMS(i) = err;
1087 }
1088 delete [] x;
1089}
1090
1091////////////////////////////////////////////////////////////////////////////////
1092/// Return the static instance.
1093
1095{
1096 return fgInstance;
1097}
1098
1099////////////////////////////////////////////////////////////////////////////////
1100/// PRIVATE METHOD:
1101/// Create list of candidate functions for the parameterisation. See
1102/// also
1103/// <a href="#TMultiDimFit:description">class description</a>
1104
1106{
1107 Int_t i = 0;
1108 Int_t j = 0;
1109 Int_t k = 0;
1110
1111 // The temporary array to store the powers in. We don't need to
1112 // initialize this array however.
1114 Int_t *powers = new Int_t[fMaxFuncNV];
1115
1116 // store of `control variables'
1117 Double_t* control = new Double_t[fMaxFunctions];
1118
1119 // We've better initialize the variables
1120 Int_t *iv = new Int_t[fNVariables];
1121 for (i = 0; i < fNVariables; i++)
1122 iv[i] = 1;
1123
1124 if (!fIsUserFunction) {
1125
1126 // Number of funcs selected
1127 Int_t numberFunctions = 0;
1128
1129 // Absolute max number of functions
1130 Int_t maxNumberFunctions = 1;
1131 for (i = 0; i < fNVariables; i++)
1132 maxNumberFunctions *= fMaxPowers[i];
1133
1134 while (kTRUE) {
1135 // Get the control value for this function
1136 Double_t s = EvalControl(iv);
1137
1138 if (s <= fPowerLimit) {
1139
1140 // Call over-loadable method Select, as to allow the user to
1141 // interfere with the selection of functions.
1142 if (Select(iv)) {
1143 numberFunctions++;
1144
1145 // If we've reached the user defined limit of how many
1146 // functions we can consider, break out of the loop
1147 if (numberFunctions > fMaxFunctions)
1148 break;
1149
1150 // Store the control value, so we can sort array of powers
1151 // later on
1152 control[numberFunctions-1] = Int_t(1.0e+6*s);
1153
1154 // Store the powers in powers array.
1155 for (i = 0; i < fNVariables; i++) {
1156 j = (numberFunctions - 1) * fNVariables + i;
1157 powers[j] = iv[i];
1158 }
1159 } // if (Select())
1160 } // if (s <= fPowerLimit)
1161
1162 for (i = 0; i < fNVariables; i++)
1163 if (iv[i] < fMaxPowers[i])
1164 break;
1165
1166 // If all variables have reached their maximum power, then we
1167 // break out of the loop
1168 if (i == fNVariables) {
1169 fMaxFunctions = numberFunctions;
1170 break;
1171 }
1172
1173 // Next power in variable i
1174 if (i < fNVariables) iv[i]++;
1175
1176 for (j = 0; j < i; j++)
1177 iv[j] = 1;
1178 } // while (kTRUE)
1179 }
1180 else {
1181 // In case the user gave an explicit function
1182 for (i = 0; i < fMaxFunctions; i++) {
1183 // Copy the powers to working arrays
1184 for (j = 0; j < fNVariables; j++) {
1185 powers[i * fNVariables + j] = fPowers[i * fNVariables + j];
1186 iv[j] = fPowers[i * fNVariables + j];
1187 }
1188
1189 control[i] = Int_t(1.0e+6*EvalControl(iv));
1190 }
1191 }
1192
1193 // Now we need to sort the powers according to least `control
1194 // variable'
1195 Int_t *order = new Int_t[fMaxFunctions];
1196 for (i = 0; i < fMaxFunctions; i++)
1197 order[i] = i;
1199 fPowers = new Int_t[fMaxFuncNV];
1200
1201 for (i = 0; i < fMaxFunctions; i++) {
1202 Double_t x = control[i];
1203 Int_t l = order[i];
1204 k = i;
1205
1206 for (j = i; j < fMaxFunctions; j++) {
1207 if (control[j] <= x) {
1208 x = control[j];
1209 l = order[j];
1210 k = j;
1211 }
1212 }
1213
1214 if (k != i) {
1215 control[k] = control[i];
1216 control[i] = x;
1217 order[k] = order[i];
1218 order[i] = l;
1219 }
1220 }
1221
1222 for (i = 0; i < fMaxFunctions; i++)
1223 for (j = 0; j < fNVariables; j++)
1224 fPowers[i * fNVariables + j] = powers[order[i] * fNVariables + j];
1225
1226 delete [] control;
1227 delete [] powers;
1228 delete [] order;
1229 delete [] iv;
1230}
1231
1232
1233////////////////////////////////////////////////////////////////////////////////
1234/// Calculate Chi square over either the test sample. The optional
1235/// argument coeff is a vector of coefficients to use in the
1236/// evaluation of the parameterisation. If coeff == 0, then the found
1237/// coefficients is used.
1238/// Used my MINUIT for fit (see TMultDimFit::Fit)
1239
1241{
1242 fChi2 = 0;
1243 Int_t i, j;
1245 for (i = 0; i < fTestSampleSize; i++) {
1246 // Get the stored point
1247 for (j = 0; j < fNVariables; j++)
1248 x[j] = fTestVariables(i * fNVariables + j);
1249
1250 // Evaluate function. Scale to shifted values
1251 Double_t f = Eval(x,coeff);
1252
1253 // Calculate contribution to Chic square
1254 fChi2 += 1. / TMath::Max(fTestSqError(i),1e-20)
1255 * (fTestQuantity(i) - f) * (fTestQuantity(i) - f);
1256 }
1257
1258 // Clean up
1259 delete [] x;
1260
1261 return fChi2;
1262}
1263
1264
1265////////////////////////////////////////////////////////////////////////////////
1266/// Generate the file <filename> with .C appended if argument doesn't
1267/// end in .cxx or .C. The contains the implementation of the
1268/// function:
1269///
1270/// Double_t <funcname>(Double_t *x)
1271///
1272/// which does the same as TMultiDimFit::Eval. Please refer to this
1273/// method.
1274///
1275/// Further, the static variables:
1276///
1277/// Int_t gNVariables
1278/// Int_t gNCoefficients
1279/// Double_t gDMean
1280/// Double_t gXMean[]
1281/// Double_t gXMin[]
1282/// Double_t gXMax[]
1283/// Double_t gCoefficient[]
1284/// Int_t gPower[]
1285///
1286/// are initialized. The only ROOT header file needed is Rtypes.h
1287///
1288/// See TMultiDimFit::MakeRealCode for a list of options
1289
1290void TMultiDimFit::MakeCode(const char* filename, Option_t *option)
1291{
1292
1293 TString outName(filename);
1294 if (!outName.EndsWith(".C") && !outName.EndsWith(".cxx"))
1295 outName += ".C";
1296
1297 MakeRealCode(outName.Data(),"",option);
1298}
1299
1300
1301
1302////////////////////////////////////////////////////////////////////////////////
1303/// PRIVATE METHOD:
1304/// Compute the errors on the coefficients. For this to be done, the
1305/// curvature matrix of the non-orthogonal functions, is computed.
1306
1308{
1309 Int_t i = 0;
1310 Int_t j = 0;
1311 Int_t k = 0;
1315
1316 TMatrixDSym curvatureMatrix(fNCoefficients);
1317
1318 // Build the curvature matrix
1319 for (i = 0; i < fNCoefficients; i++) {
1320 iF = TMatrixDRow(fFunctions,i);
1321 for (j = 0; j <= i; j++) {
1322 jF = TMatrixDRow(fFunctions,j);
1323 for (k = 0; k < fSampleSize; k++)
1324 curvatureMatrix(i,j) +=
1325 1 / TMath::Max(fSqError(k), 1e-20) * iF(k) * jF(k);
1326 curvatureMatrix(j,i) = curvatureMatrix(i,j);
1327 }
1328 }
1329
1330 // Calculate Chi Square
1331 fChi2 = 0;
1332 for (i = 0; i < fSampleSize; i++) {
1333 Double_t f = 0;
1334 for (j = 0; j < fNCoefficients; j++)
1335 f += fCoefficients(j) * fFunctions(j,i);
1336 fChi2 += 1. / TMath::Max(fSqError(i),1e-20) * (fQuantity(i) - f)
1337 * (fQuantity(i) - f);
1338 }
1339
1340 // Invert the curvature matrix
1341 const TVectorD diag = TMatrixDDiag_const(curvatureMatrix);
1342 curvatureMatrix.NormByDiag(diag);
1343
1344 TDecompChol chol(curvatureMatrix);
1345 if (!chol.Decompose())
1346 Error("MakeCoefficientErrors", "curvature matrix is singular");
1347 chol.Invert(curvatureMatrix);
1348
1349 curvatureMatrix.NormByDiag(diag);
1350
1351 for (i = 0; i < fNCoefficients; i++)
1352 fCoefficientsRMS(i) = TMath::Sqrt(curvatureMatrix(i,i));
1353}
1354
1355
1356////////////////////////////////////////////////////////////////////////////////
1357/// PRIVATE METHOD:
1358/// Invert the model matrix B, and compute final coefficients. For a
1359/// more thorough discussion of what this means, please refer to the
1360/// <a href="#TMultiDimFit:description">class description</a>
1361///
1362/// First we invert the lower triangle matrix fOrthCurvatureMatrix
1363/// and store the inverted matrix in the upper triangle.
1364
1366{
1367 Int_t i = 0, j = 0;
1368 Int_t col = 0, row = 0;
1369
1370 // Invert the B matrix
1371 for (col = 1; col < fNCoefficients; col++) {
1372 for (row = col - 1; row > -1; row--) {
1373 fOrthCurvatureMatrix(row,col) = 0;
1374 for (i = row; i <= col ; i++)
1375 fOrthCurvatureMatrix(row,col) -=
1377 * fOrthCurvatureMatrix(i,col);
1378 }
1379 }
1380
1381 // Compute the final coefficients
1383
1384 for (i = 0; i < fNCoefficients; i++) {
1385 Double_t sum = 0;
1386 for (j = i; j < fNCoefficients; j++)
1388 fCoefficients(i) = sum;
1389 }
1390
1391 // Compute the final residuals
1393 for (i = 0; i < fSampleSize; i++)
1394 fResiduals(i) = fQuantity(i);
1395
1396 for (i = 0; i < fNCoefficients; i++)
1397 for (j = 0; j < fSampleSize; j++)
1398 fResiduals(j) -= fCoefficients(i) * fFunctions(i,j);
1399
1400 // Compute the max and minimum, and squared sum of the evaluated
1401 // residuals
1402 fMinResidual = 10e10;
1403 fMaxResidual = -10e10;
1404 Double_t sqRes = 0;
1405 for (i = 0; i < fSampleSize; i++){
1406 sqRes += fResiduals(i) * fResiduals(i);
1407 if (fResiduals(i) <= fMinResidual) {
1409 fMinResidualRow = i;
1410 }
1411 if (fResiduals(i) >= fMaxResidual) {
1413 fMaxResidualRow = i;
1414 }
1415 }
1416
1419
1420 // If we use histograms, fill some more
1424 for (i = 0; i < fSampleSize; i++) {
1426 ((TH2D*)fHistograms->FindObject("res_d"))->Fill(fQuantity(i),
1427 fResiduals(i));
1429 ((TH1D*)fHistograms->FindObject("res_train"))->Fill(fResiduals(i));
1430
1432 for (j = 0; j < fNVariables; j++)
1433 ((TH2D*)fHistograms->FindObject(Form("res_x_%d",j)))
1434 ->Fill(fVariables(i * fNVariables + j),fResiduals(i));
1435 }
1436 } // If histograms
1437
1438}
1439
1440
1441////////////////////////////////////////////////////////////////////////////////
1442/// PRIVATE METHOD:
1443/// Compute the correlation matrix
1444
1446{
1447 if (!fShowCorrelation)
1448 return;
1449
1451
1452 Double_t d2 = 0;
1453 Double_t ddotXi = 0; // G.Q. needs to be reinitialized in the loop over i fNVariables
1454 Double_t xiNorm = 0; // G.Q. needs to be reinitialized in the loop over i fNVariables
1455 Double_t xidotXj = 0; // G.Q. needs to be reinitialized in the loop over j fNVariables
1456 Double_t xjNorm = 0; // G.Q. needs to be reinitialized in the loop over j fNVariables
1457
1458 Int_t i, j, k, l, m; // G.Q. added m variable
1459 for (i = 0; i < fSampleSize; i++)
1460 d2 += fQuantity(i) * fQuantity(i);
1461
1462 for (i = 0; i < fNVariables; i++) {
1463 ddotXi = 0.; // G.Q. reinitialisation
1464 xiNorm = 0.; // G.Q. reinitialisation
1465 for (j = 0; j< fSampleSize; j++) {
1466 // Index of sample j of variable i
1467 k = j * fNVariables + i;
1468 ddotXi += fQuantity(j) * (fVariables(k) - fMeanVariables(i));
1469 xiNorm += (fVariables(k) - fMeanVariables(i))
1470 * (fVariables(k) - fMeanVariables(i));
1471 }
1472 fCorrelationMatrix(i,0) = ddotXi / TMath::Sqrt(d2 * xiNorm);
1473
1474 for (j = 0; j < i; j++) {
1475 xidotXj = 0.; // G.Q. reinitialisation
1476 xjNorm = 0.; // G.Q. reinitialisation
1477 for (k = 0; k < fSampleSize; k++) {
1478 // Index of sample j of variable i
1479 // l = j * fNVariables + k; // G.Q.
1480 l = k * fNVariables + j; // G.Q.
1481 m = k * fNVariables + i; // G.Q.
1482 // G.Q. xidotXj += (fVariables(i) - fMeanVariables(i))
1483 // G.Q. * (fVariables(l) - fMeanVariables(j));
1484 xidotXj += (fVariables(m) - fMeanVariables(i))
1485 * (fVariables(l) - fMeanVariables(j)); // G.Q. modified index for Xi
1486 xjNorm += (fVariables(l) - fMeanVariables(j))
1487 * (fVariables(l) - fMeanVariables(j));
1488 }
1489 //fCorrelationMatrix(i+1,j) = xidotXj / TMath::Sqrt(xiNorm * xjNorm);
1490 fCorrelationMatrix(i,j+1) = xidotXj / TMath::Sqrt(xiNorm * xjNorm);
1491 }
1492 }
1493}
1494
1495
1496
1497////////////////////////////////////////////////////////////////////////////////
1498/// PRIVATE METHOD:
1499/// Make Gram-Schmidt orthogonalisation. The class description gives
1500/// a thorough account of this algorithm, as well as
1501/// references. Please refer to the
1502/// <a href="#TMultiDimFit:description">class description</a>
1503
1505{
1506
1507 // calculate w_i, that is, evaluate the current function at data
1508 // point i
1509 Double_t f2 = 0;
1512 Int_t j = 0;
1513 Int_t k = 0;
1514
1515 for (j = 0; j < fSampleSize; j++) {
1516 fFunctions(fNCoefficients, j) = 1;
1518 // First, however, we need to calculate f_fNCoefficients
1519 for (k = 0; k < fNVariables; k++) {
1520 Int_t p = fPowers[function * fNVariables + k];
1521 Double_t x = fVariables(j * fNVariables + k);
1523 }
1524
1525 // Calculate f dot f in f2
1527 // Assign to w_fNCoefficients f_fNCoefficients
1529 }
1530
1531 // the first column of w is equal to f
1532 for (j = 0; j < fNCoefficients; j++) {
1533 Double_t fdw = 0;
1534 // Calculate (f_fNCoefficients dot w_j) / w_j^2
1535 for (k = 0; k < fSampleSize; k++) {
1536 fdw += fFunctions(fNCoefficients, k) * fOrthFunctions(j,k)
1537 / fOrthFunctionNorms(j);
1538 }
1539
1541 // and subtract it from the current value of w_ij
1542 for (k = 0; k < fSampleSize; k++)
1544 }
1545
1546 for (j = 0; j < fSampleSize; j++) {
1547 // calculate squared length of w_fNCoefficients
1551
1552 // calculate D dot w_fNCoefficients in A
1555 }
1556
1557 // First test, but only if didn't user specify
1558 if (!fIsUserFunction)
1561 return 0;
1562
1563 // The result found by this code for the first residual is always
1564 // much less then the one found be MUDIFI. That's because it's
1565 // supposed to be. The cause is the improved precision of Double_t
1566 // over DOUBLE PRECISION!
1570
1571 // Calculate the residual from including this fNCoefficients.
1573
1574 return dResidur;
1575}
1576
1577
1578////////////////////////////////////////////////////////////////////////////////
1579/// Make histograms of the result of the analysis. This message
1580/// should be sent after having read all data points, but before
1581/// finding the parameterization
1582///
1583/// Options:
1584/// A All the below
1585/// X Original independent variables
1586/// D Original dependent variables
1587/// N Normalised independent variables
1588/// S Shifted dependent variables
1589/// R1 Residuals versus normalised independent variables
1590/// R2 Residuals versus dependent variable
1591/// R3 Residuals computed on training sample
1592/// R4 Residuals computed on test sample
1593///
1594/// For a description of these quantities, refer to
1595/// <a href="#TMultiDimFit:description">class description</a>
1596
1598{
1599 TString opt(option);
1600 opt.ToLower();
1601
1602 if (opt.Length() < 1)
1603 return;
1604
1605 if (!fHistograms)
1606 fHistograms = new TList;
1607
1608 // Counter variable
1609 Int_t i = 0;
1610
1611 // Histogram of original variables
1612 if (opt.Contains("x") || opt.Contains("a")) {
1614 for (i = 0; i < fNVariables; i++)
1615 if (!fHistograms->FindObject(Form("x_%d_orig",i)))
1616 fHistograms->Add(new TH1D(Form("x_%d_orig",i),
1617 Form("Original variable # %d",i),
1619 fMaxVariables(i)));
1620 }
1621
1622 // Histogram of original dependent variable
1623 if (opt.Contains("d") || opt.Contains("a")) {
1625 if (!fHistograms->FindObject("d_orig"))
1626 fHistograms->Add(new TH1D("d_orig", "Original Quantity",
1628 }
1629
1630 // Histograms of normalized variables
1631 if (opt.Contains("n") || opt.Contains("a")) {
1633 for (i = 0; i < fNVariables; i++)
1634 if (!fHistograms->FindObject(Form("x_%d_norm",i)))
1635 fHistograms->Add(new TH1D(Form("x_%d_norm",i),
1636 Form("Normalized variable # %d",i),
1637 fBinVarX, -1,1));
1638 }
1639
1640 // Histogram of shifted dependent variable
1641 if (opt.Contains("s") || opt.Contains("a")) {
1643 if (!fHistograms->FindObject("d_shifted"))
1644 fHistograms->Add(new TH1D("d_shifted", "Shifted Quantity",
1647 }
1648
1649 // Residual from training sample versus independent variables
1650 if (opt.Contains("r1") || opt.Contains("a")) {
1652 for (i = 0; i < fNVariables; i++)
1653 if (!fHistograms->FindObject(Form("res_x_%d",i)))
1654 fHistograms->Add(new TH2D(Form("res_x_%d",i),
1655 Form("Computed residual versus x_%d", i),
1656 fBinVarX, -1, 1,
1657 fBinVarY,
1660 }
1661
1662 // Residual from training sample versus. dependent variable
1663 if (opt.Contains("r2") || opt.Contains("a")) {
1665 if (!fHistograms->FindObject("res_d"))
1666 fHistograms->Add(new TH2D("res_d",
1667 "Computed residuals vs Quantity",
1668 fBinVarX,
1671 fBinVarY,
1674 }
1675
1676 // Residual from training sample
1677 if (opt.Contains("r3") || opt.Contains("a")) {
1679 if (!fHistograms->FindObject("res_train"))
1680 fHistograms->Add(new TH1D("res_train",
1681 "Computed residuals over training sample",
1684
1685 }
1686 if (opt.Contains("r4") || opt.Contains("a")) {
1688 if (!fHistograms->FindObject("res_test"))
1689 fHistograms->Add(new TH1D("res_test",
1690 "Distribution of residuals from test",
1693 }
1694}
1695
1696
1697////////////////////////////////////////////////////////////////////////////////
1698/// Generate the file <classname>MDF.cxx which contains the
1699/// implementation of the method:
1700///
1701/// Double_t <classname>::MDF(Double_t *x)
1702///
1703/// which does the same as TMultiDimFit::Eval. Please refer to this
1704/// method.
1705///
1706/// Further, the public static members:
1707///
1708/// Int_t <classname>::fgNVariables
1709/// Int_t <classname>::fgNCoefficients
1710/// Double_t <classname>::fgDMean
1711/// Double_t <classname>::fgXMean[] //[fgNVariables]
1712/// Double_t <classname>::fgXMin[] //[fgNVariables]
1713/// Double_t <classname>::fgXMax[] //[fgNVariables]
1714/// Double_t <classname>::fgCoefficient[] //[fgNCoeffficents]
1715/// Int_t <classname>::fgPower[] //[fgNCoeffficents*fgNVariables]
1716///
1717/// are initialized, and assumed to exist. The class declaration is
1718/// assumed to be in <classname>.h and assumed to be provided by the
1719/// user.
1720///
1721/// See TMultiDimFit::MakeRealCode for a list of options
1722///
1723/// The minimal class definition is:
1724///
1725/// class <classname> {
1726/// public:
1727/// Int_t <classname>::fgNVariables; // Number of variables
1728/// Int_t <classname>::fgNCoefficients; // Number of terms
1729/// Double_t <classname>::fgDMean; // Mean from training sample
1730/// Double_t <classname>::fgXMean[]; // Mean from training sample
1731/// Double_t <classname>::fgXMin[]; // Min from training sample
1732/// Double_t <classname>::fgXMax[]; // Max from training sample
1733/// Double_t <classname>::fgCoefficient[]; // Coefficients
1734/// Int_t <classname>::fgPower[]; // Function powers
1735///
1736/// Double_t Eval(Double_t *x);
1737/// };
1738///
1739/// Whether the method <classname>::Eval should be static or not, is
1740/// up to the user.
1741
1742void TMultiDimFit::MakeMethod(const Char_t* classname, Option_t* option)
1743{
1744 MakeRealCode(Form("%sMDF.cxx", classname), classname, option);
1745}
1746
1747
1748
1749////////////////////////////////////////////////////////////////////////////////
1750/// PRIVATE METHOD:
1751/// Normalize data to the interval [-1;1]. This is needed for the
1752/// classes method to work.
1753
1755{
1756 Int_t i = 0;
1757 Int_t j = 0;
1758 Int_t k = 0;
1759
1760 for (i = 0; i < fSampleSize; i++) {
1762 ((TH1D*)fHistograms->FindObject("d_orig"))->Fill(fQuantity(i));
1763
1766
1768 ((TH1D*)fHistograms->FindObject("d_shifted"))->Fill(fQuantity(i));
1769
1770 for (j = 0; j < fNVariables; j++) {
1771 Double_t range = 1. / (fMaxVariables(j) - fMinVariables(j));
1772 k = i * fNVariables + j;
1773
1774 // Fill histograms of original independent variables
1776 ((TH1D*)fHistograms->FindObject(Form("x_%d_orig",j)))
1777 ->Fill(fVariables(k));
1778
1779 // Normalise independent variables
1780 fVariables(k) = 1 + 2 * range * (fVariables(k) - fMaxVariables(j));
1781
1782 // Fill histograms of normalised independent variables
1784 ((TH1D*)fHistograms->FindObject(Form("x_%d_norm",j)))
1785 ->Fill(fVariables(k));
1786
1787 }
1788 }
1789 // Shift min and max of dependent variable
1792
1793 // Shift mean of independent variables
1794 for (i = 0; i < fNVariables; i++) {
1795 Double_t range = 1. / (fMaxVariables(i) - fMinVariables(i));
1796 fMeanVariables(i) = 1 + 2 * range * (fMeanVariables(i)
1797 - fMaxVariables(i));
1798 }
1799}
1800
1801
1802////////////////////////////////////////////////////////////////////////////////
1803/// PRIVATE METHOD:
1804/// Find the parameterization over the training sample. A full account
1805/// of the algorithm is given in the
1806/// <a href="#TMultiDimFit:description">class description</a>
1807
1809{
1810 Int_t i = -1;
1811 Int_t j = 0;
1812 Int_t k = 0;
1813 Int_t maxPass = 3;
1814 Int_t studied = 0;
1815 Double_t squareResidual = fSumSqAvgQuantity;
1816 fNCoefficients = 0;
1823 fFunctions = 1;
1824
1827 Int_t l;
1828 for (l=0;l<fMaxFunctions;l++) fFunctionCodes[l] = 0;
1829 for (l=0;l<fMaxTerms;l++) fPowerIndex[l] = 0;
1830
1831 if (fMaxAngle != 0) maxPass = 100;
1832 if (fIsUserFunction) maxPass = 1;
1833
1834 // Loop over the number of functions we want to study.
1835 // increment inspection counter
1836 while(kTRUE) {
1837
1838 // Reach user defined limit of studies
1839 if (studied++ >= fMaxStudy) {
1841 break;
1842 }
1843
1844 // Considered all functions several times
1845 if (k >= maxPass) {
1847 break;
1848 }
1849
1850 // increment function counter
1851 i++;
1852
1853 // If we've reached the end of the functions, restart pass
1854 if (i == fMaxFunctions) {
1855 if (fMaxAngle != 0)
1856 fMaxAngle += (90 - fMaxAngle) / 2;
1857 i = 0;
1858 studied--;
1859 k++;
1860 continue;
1861 }
1862 if (studied == 1)
1863 fFunctionCodes[i] = 0;
1864 else if (fFunctionCodes[i] >= 2)
1865 continue;
1866
1867 // Print a happy message
1868 if (fIsVerbose && studied == 1)
1869 std::cout << "Coeff SumSqRes Contrib Angle QM Func"
1870 << " Value W^2 Powers" << std::endl;
1871
1872 // Make the Gram-Schmidt
1873 Double_t dResidur = MakeGramSchmidt(i);
1874
1875 if (dResidur == 0) {
1876 // This function is no good!
1877 // First test is in MakeGramSchmidt
1878 fFunctionCodes[i] = 1;
1879 continue;
1880 }
1881
1882 // If user specified function, assume they know what they are doing
1883 if (!fIsUserFunction) {
1884 // Flag this function as considered
1885 fFunctionCodes[i] = 2;
1886
1887 // Test if this function contributes to the fit
1888 if (!TestFunction(squareResidual, dResidur)) {
1889 fFunctionCodes[i] = 1;
1890 continue;
1891 }
1892 }
1893
1894 // If we get to here, the function currently considered is
1895 // fNCoefficients, so we increment the counter
1896 // Flag this function as OK, and store and the number in the
1897 // index.
1898 fFunctionCodes[i] = 3;
1901
1902 // We add the current contribution to the sum of square of
1903 // residuals;
1904 squareResidual -= dResidur;
1905
1906
1907 // Calculate control parameter from this function
1908 for (j = 0; j < fNVariables; j++) {
1909 if (fNCoefficients == 1
1910 || fMaxPowersFinal[j] <= fPowers[i * fNVariables + j] - 1)
1911 fMaxPowersFinal[j] = fPowers[i * fNVariables + j] - 1;
1912 }
1914
1915 // Print the statistics about this function
1916 if (fIsVerbose) {
1917 std::cout << std::setw(5) << fNCoefficients << " "
1918 << std::setw(10) << std::setprecision(4) << squareResidual << " "
1919 << std::setw(10) << std::setprecision(4) << dResidur << " "
1920 << std::setw(7) << std::setprecision(3) << fMaxAngle << " "
1921 << std::setw(7) << std::setprecision(3) << s << " "
1922 << std::setw(5) << i << " "
1923 << std::setw(10) << std::setprecision(4)
1925 << std::setw(10) << std::setprecision(4)
1927 << std::flush;
1928 for (j = 0; j < fNVariables; j++)
1929 std::cout << " " << fPowers[i * fNVariables + j] - 1 << std::flush;
1930 std::cout << std::endl;
1931 }
1932
1933 if (fNCoefficients >= fMaxTerms /* && fIsVerbose */) {
1935 break;
1936 }
1937
1938 Double_t err = TMath::Sqrt(TMath::Max(1e-20,squareResidual) /
1940 if (err < fMinRelativeError) {
1942 break;
1943 }
1944
1945 }
1946
1947 fError = TMath::Max(1e-20,squareResidual);
1950}
1951
1952
1953////////////////////////////////////////////////////////////////////////////////
1954/// PRIVATE METHOD:
1955/// This is the method that actually generates the code for the
1956/// evaluation the parameterization on some point.
1957/// It's called by TMultiDimFit::MakeCode and TMultiDimFit::MakeMethod.
1958///
1959/// The options are: NONE so far
1960
1961void TMultiDimFit::MakeRealCode(const char *filename,
1962 const char *classname,
1963 Option_t *)
1964{
1965 Int_t i, j;
1966
1967 Bool_t isMethod = (classname[0] == '\0' ? kFALSE : kTRUE);
1968 const char *prefix = (isMethod ? Form("%s::", classname) : "");
1969 const char *cv_qual = (isMethod ? "" : "static ");
1970
1971 std::ofstream outFile(filename,std::ios::out|std::ios::trunc);
1972 if (!outFile) {
1973 Error("MakeRealCode","couldn't open output file '%s'",filename);
1974 return;
1975 }
1976
1977 if (fIsVerbose)
1978 std::cout << "Writing on file \"" << filename << "\" ... " << std::flush;
1979 //
1980 // Write header of file
1981 //
1982 // Emacs mode line ;-)
1983 outFile << "// -*- mode: c++ -*-" << std::endl;
1984 // Info about creator
1985 outFile << "// " << std::endl
1986 << "// File " << filename
1987 << " generated by TMultiDimFit::MakeRealCode" << std::endl;
1988 // Time stamp
1989 TDatime date;
1990 outFile << "// on " << date.AsString() << std::endl;
1991 // ROOT version info
1992 outFile << "// ROOT version " << gROOT->GetVersion()
1993 << std::endl << "//" << std::endl;
1994 // General information on the code
1995 outFile << "// This file contains the function " << std::endl
1996 << "//" << std::endl
1997 << "// double " << prefix << "MDF(double *x); " << std::endl
1998 << "//" << std::endl
1999 << "// For evaluating the parameterization obtained" << std::endl
2000 << "// from TMultiDimFit and the point x" << std::endl
2001 << "// " << std::endl
2002 << "// See TMultiDimFit class documentation for more "
2003 << "information " << std::endl << "// " << std::endl;
2004 // Header files
2005 if (isMethod)
2006 // If these are methods, we need the class header
2007 outFile << "#include \"" << classname << ".h\"" << std::endl;
2008
2009 //
2010 // Now for the data
2011 //
2012 outFile << "//" << std::endl
2013 << "// Static data variables" << std::endl
2014 << "//" << std::endl;
2015 outFile << cv_qual << "int " << prefix << "gNVariables = "
2016 << fNVariables << ";" << std::endl;
2017 outFile << cv_qual << "int " << prefix << "gNCoefficients = "
2018 << fNCoefficients << ";" << std::endl;
2019 outFile << cv_qual << "double " << prefix << "gDMean = "
2020 << fMeanQuantity << ";" << std::endl;
2021
2022 // Assignment to mean vector.
2023 outFile << "// Assignment to mean vector." << std::endl;
2024 outFile << cv_qual << "double " << prefix
2025 << "gXMean[] = {" << std::endl;
2026 for (i = 0; i < fNVariables; i++)
2027 outFile << (i != 0 ? ", " : " ") << fMeanVariables(i) << std::flush;
2028 outFile << " };" << std::endl << std::endl;
2029
2030 // Assignment to minimum vector.
2031 outFile << "// Assignment to minimum vector." << std::endl;
2032 outFile << cv_qual << "double " << prefix
2033 << "gXMin[] = {" << std::endl;
2034 for (i = 0; i < fNVariables; i++)
2035 outFile << (i != 0 ? ", " : " ") << fMinVariables(i) << std::flush;
2036 outFile << " };" << std::endl << std::endl;
2037
2038 // Assignment to maximum vector.
2039 outFile << "// Assignment to maximum vector." << std::endl;
2040 outFile << cv_qual << "double " << prefix
2041 << "gXMax[] = {" << std::endl;
2042 for (i = 0; i < fNVariables; i++)
2043 outFile << (i != 0 ? ", " : " ") << fMaxVariables(i) << std::flush;
2044 outFile << " };" << std::endl << std::endl;
2045
2046 // Assignment to coefficients vector.
2047 outFile << "// Assignment to coefficients vector." << std::endl;
2048 outFile << cv_qual << "double " << prefix
2049 << "gCoefficient[] = {" << std::flush;
2050 for (i = 0; i < fNCoefficients; i++)
2051 outFile << (i != 0 ? "," : "") << std::endl
2052 << " " << fCoefficients(i) << std::flush;
2053 outFile << std::endl << " };" << std::endl << std::endl;
2054
2055 // Assignment to error coefficients vector.
2056 outFile << "// Assignment to error coefficients vector." << std::endl;
2057 outFile << cv_qual << "double " << prefix
2058 << "gCoefficientRMS[] = {" << std::flush;
2059 for (i = 0; i < fNCoefficients; i++)
2060 outFile << (i != 0 ? "," : "") << std::endl
2061 << " " << fCoefficientsRMS(i) << std::flush;
2062 outFile << std::endl << " };" << std::endl << std::endl;
2063
2064 // Assignment to powers vector.
2065 outFile << "// Assignment to powers vector." << std::endl
2066 << "// The powers are stored row-wise, that is" << std::endl
2067 << "// p_ij = " << prefix
2068 << "gPower[i * NVariables + j];" << std::endl;
2069 outFile << cv_qual << "int " << prefix
2070 << "gPower[] = {" << std::flush;
2071 for (i = 0; i < fNCoefficients; i++) {
2072 for (j = 0; j < fNVariables; j++) {
2073 if (j != 0) outFile << std::flush << " ";
2074 else outFile << std::endl << " ";
2075 outFile << fPowers[fPowerIndex[i] * fNVariables + j]
2076 << (i == fNCoefficients - 1 && j == fNVariables - 1 ? "" : ",")
2077 << std::flush;
2078 }
2079 }
2080 outFile << std::endl << "};" << std::endl << std::endl;
2081
2082
2083 //
2084 // Finally we reach the function itself
2085 //
2086 outFile << "// " << std::endl
2087 << "// The "
2088 << (isMethod ? "method " : "function ")
2089 << " double " << prefix
2090 << "MDF(double *x)"
2091 << std::endl << "// " << std::endl;
2092 outFile << "double " << prefix
2093 << "MDF(double *x) {" << std::endl
2094 << " double returnValue = " << prefix << "gDMean;" << std::endl
2095 << " int i = 0, j = 0, k = 0;" << std::endl
2096 << " for (i = 0; i < " << prefix << "gNCoefficients ; i++) {"
2097 << std::endl
2098 << " // Evaluate the ith term in the expansion" << std::endl
2099 << " double term = " << prefix << "gCoefficient[i];"
2100 << std::endl
2101 << " for (j = 0; j < " << prefix << "gNVariables; j++) {"
2102 << std::endl
2103 << " // Evaluate the polynomial in the jth variable." << std::endl
2104 << " int power = "<< prefix << "gPower["
2105 << prefix << "gNVariables * i + j]; " << std::endl
2106 << " double p1 = 1, p2 = 0, p3 = 0, r = 0;" << std::endl
2107 << " double v = 1 + 2. / ("
2108 << prefix << "gXMax[j] - " << prefix
2109 << "gXMin[j]) * (x[j] - " << prefix << "gXMax[j]);" << std::endl
2110 << " // what is the power to use!" << std::endl
2111 << " switch(power) {" << std::endl
2112 << " case 1: r = 1; break; " << std::endl
2113 << " case 2: r = v; break; " << std::endl
2114 << " default: " << std::endl
2115 << " p2 = v; " << std::endl
2116 << " for (k = 3; k <= power; k++) { " << std::endl
2117 << " p3 = p2 * v;" << std::endl;
2118 if (fPolyType == kLegendre)
2119 outFile << " p3 = ((2 * i - 3) * p2 * v - (i - 2) * p1)"
2120 << " / (i - 1);" << std::endl;
2121 if (fPolyType == kChebyshev)
2122 outFile << " p3 = 2 * v * p2 - p1; " << std::endl;
2123 outFile << " p1 = p2; p2 = p3; " << std::endl << " }" << std::endl
2124 << " r = p3;" << std::endl << " }" << std::endl
2125 << " // multiply this term by the poly in the jth var" << std::endl
2126 << " term *= r; " << std::endl << " }" << std::endl
2127 << " // Add this term to the final result" << std::endl
2128 << " returnValue += term;" << std::endl << " }" << std::endl
2129 << " return returnValue;" << std::endl << "}" << std::endl << std::endl;
2130
2131 // EOF
2132 outFile << "// EOF for " << filename << std::endl;
2133
2134 // Close the file
2135 outFile.close();
2136
2137 if (fIsVerbose)
2138 std::cout << "done" << std::endl;
2139}
2140
2141
2142////////////////////////////////////////////////////////////////////////////////
2143/// Print statistics etc.
2144/// Options are
2145/// P Parameters
2146/// S Statistics
2147/// C Coefficients
2148/// R Result of parameterisation
2149/// F Result of fit
2150/// K Correlation Matrix
2151/// M Pretty print formula
2152///
2153
2155{
2156 Int_t i = 0;
2157 Int_t j = 0;
2158
2159 TString opt(option);
2160 opt.ToLower();
2161
2162 if (opt.Contains("p")) {
2163 // Print basic parameters for this object
2164 std::cout << "User parameters:" << std::endl
2165 << "----------------" << std::endl
2166 << " Variables: " << fNVariables << std::endl
2167 << " Data points: " << fSampleSize << std::endl
2168 << " Max Terms: " << fMaxTerms << std::endl
2169 << " Power Limit Parameter: " << fPowerLimit << std::endl
2170 << " Max functions: " << fMaxFunctions << std::endl
2171 << " Max functions to study: " << fMaxStudy << std::endl
2172 << " Max angle (optional): " << fMaxAngle << std::endl
2173 << " Min angle: " << fMinAngle << std::endl
2174 << " Relative Error accepted: " << fMinRelativeError << std::endl
2175 << " Maximum Powers: " << std::flush;
2176 for (i = 0; i < fNVariables; i++)
2177 std::cout << " " << fMaxPowers[i] - 1 << std::flush;
2178 std::cout << std::endl << std::endl
2179 << " Parameterisation will be done using " << std::flush;
2180 if (fPolyType == kChebyshev)
2181 std::cout << "Chebyshev polynomials" << std::endl;
2182 else if (fPolyType == kLegendre)
2183 std::cout << "Legendre polynomials" << std::endl;
2184 else
2185 std::cout << "Monomials" << std::endl;
2186 std::cout << std::endl;
2187 }
2188
2189 if (opt.Contains("s")) {
2190 // Print statistics for read data
2191 std::cout << "Sample statistics:" << std::endl
2192 << "------------------" << std::endl
2193 << " D" << std::flush;
2194 for (i = 0; i < fNVariables; i++)
2195 std::cout << " " << std::setw(10) << i+1 << std::flush;
2196 std::cout << std::endl << " Max: " << std::setw(10) << std::setprecision(7)
2197 << fMaxQuantity << std::flush;
2198 for (i = 0; i < fNVariables; i++)
2199 std::cout << " " << std::setw(10) << std::setprecision(4)
2200 << fMaxVariables(i) << std::flush;
2201 std::cout << std::endl << " Min: " << std::setw(10) << std::setprecision(7)
2202 << fMinQuantity << std::flush;
2203 for (i = 0; i < fNVariables; i++)
2204 std::cout << " " << std::setw(10) << std::setprecision(4)
2205 << fMinVariables(i) << std::flush;
2206 std::cout << std::endl << " Mean: " << std::setw(10) << std::setprecision(7)
2207 << fMeanQuantity << std::flush;
2208 for (i = 0; i < fNVariables; i++)
2209 std::cout << " " << std::setw(10) << std::setprecision(4)
2210 << fMeanVariables(i) << std::flush;
2211 std::cout << std::endl << " Function Sum Squares: " << fSumSqQuantity
2212 << std::endl << std::endl;
2213 }
2214
2215 if (opt.Contains("r")) {
2216 std::cout << "Results of Parameterisation:" << std::endl
2217 << "----------------------------" << std::endl
2218 << " Total reduction of square residuals "
2219 << fSumSqResidual << std::endl
2220 << " Relative precision obtained: "
2221 << fPrecision << std::endl
2222 << " Error obtained: "
2223 << fError << std::endl
2224 << " Multiple correlation coefficient: "
2225 << fCorrelationCoeff << std::endl
2226 << " Reduced Chi square over sample: "
2227 << fChi2 / (fSampleSize - fNCoefficients) << std::endl
2228 << " Maximum residual value: "
2229 << fMaxResidual << std::endl
2230 << " Minimum residual value: "
2231 << fMinResidual << std::endl
2232 << " Estimated root mean square: "
2233 << fRMS << std::endl
2234 << " Maximum powers used: " << std::flush;
2235 for (j = 0; j < fNVariables; j++)
2236 std::cout << fMaxPowersFinal[j] << " " << std::flush;
2237 std::cout << std::endl
2238 << " Function codes of candidate functions." << std::endl
2239 << " 1: considered,"
2240 << " 2: too little contribution,"
2241 << " 3: accepted." << std::flush;
2242 for (i = 0; i < fMaxFunctions; i++) {
2243 if (i % 60 == 0)
2244 std::cout << std::endl << " " << std::flush;
2245 else if (i % 10 == 0)
2246 std::cout << " " << std::flush;
2247 std::cout << fFunctionCodes[i];
2248 }
2249 std::cout << std::endl << " Loop over candidates stopped because " << std::flush;
2250 switch(fParameterisationCode){
2251 case PARAM_MAXSTUDY:
2252 std::cout << "max allowed studies reached" << std::endl; break;
2253 case PARAM_SEVERAL:
2254 std::cout << "all candidates considered several times" << std::endl; break;
2255 case PARAM_RELERR:
2256 std::cout << "wanted relative error obtained" << std::endl; break;
2257 case PARAM_MAXTERMS:
2258 std::cout << "max number of terms reached" << std::endl; break;
2259 default:
2260 std::cout << "some unknown reason" << std::endl;
2261 break;
2262 }
2263 std::cout << std::endl;
2264 }
2265
2266 if (opt.Contains("f")) {
2267 std::cout << "Results of Fit:" << std::endl
2268 << "---------------" << std::endl
2269 << " Test sample size: "
2270 << fTestSampleSize << std::endl
2271 << " Multiple correlation coefficient: "
2272 << fTestCorrelationCoeff << std::endl
2273 << " Relative precision obtained: "
2274 << fTestPrecision << std::endl
2275 << " Error obtained: "
2276 << fTestError << std::endl
2277 << " Reduced Chi square over sample: "
2278 << fChi2 / (fSampleSize - fNCoefficients) << std::endl
2279 << std::endl;
2280 if (fFitter) {
2281 fFitter->PrintResults(1,1);
2282 std::cout << std::endl;
2283 }
2284 }
2285
2286 if (opt.Contains("c")){
2287 std::cout << "Coefficients:" << std::endl
2288 << "-------------" << std::endl
2289 << " # Value Error Powers" << std::endl
2290 << " ---------------------------------------" << std::endl;
2291 for (i = 0; i < fNCoefficients; i++) {
2292 std::cout << " " << std::setw(3) << i << " "
2293 << std::setw(12) << fCoefficients(i) << " "
2294 << std::setw(12) << fCoefficientsRMS(i) << " " << std::flush;
2295 for (j = 0; j < fNVariables; j++)
2296 std::cout << " " << std::setw(3)
2297 << fPowers[fPowerIndex[i] * fNVariables + j] - 1 << std::flush;
2298 std::cout << std::endl;
2299 }
2300 std::cout << std::endl;
2301 }
2302 if (opt.Contains("k") && fCorrelationMatrix.IsValid()) {
2303 std::cout << "Correlation Matrix:" << std::endl
2304 << "-------------------";
2306 }
2307
2308 if (opt.Contains("m")) {
2309 std::cout << "Parameterization:" << std::endl
2310 << "-----------------" << std::endl
2311 << " Normalised variables: " << std::endl;
2312 for (i = 0; i < fNVariables; i++)
2313 std::cout << "\ty_" << i << "\t= 1 + 2 * (x_" << i << " - "
2314 << fMaxVariables(i) << ") / ("
2315 << fMaxVariables(i) << " - " << fMinVariables(i) << ")"
2316 << std::endl;
2317 std::cout << std::endl
2318 << " f(";
2319 for (i = 0; i < fNVariables; i++) {
2320 std::cout << "y_" << i;
2321 if (i != fNVariables-1) std::cout << ", ";
2322 }
2323 std::cout << ") = ";
2324 for (i = 0; i < fNCoefficients; i++) {
2325 if (i != 0)
2326 std::cout << std::endl << "\t" << (fCoefficients(i) < 0 ? "- " : "+ ")
2328 else
2329 std::cout << fCoefficients(i);
2330 for (j = 0; j < fNVariables; j++) {
2331 Int_t p = fPowers[fPowerIndex[i] * fNVariables + j];
2332 switch (p) {
2333 case 1: break;
2334 case 2: std::cout << " * y_" << j; break;
2335 default:
2336 switch(fPolyType) {
2337 case kLegendre: std::cout << " * L_" << p-1 << "(y_" << j << ")"; break;
2338 case kChebyshev: std::cout << " * C_" << p-1 << "(y_" << j << ")"; break;
2339 default: std::cout << " * y_" << j << "^" << p-1; break;
2340 }
2341 }
2342
2343 }
2344 }
2345 std::cout << std::endl;
2346 }
2347}
2348
2349
2350////////////////////////////////////////////////////////////////////////////////
2351/// Selection method. User can override this method for specialized
2352/// selection of acceptable functions in fit. Default is to select
2353/// all. This message is sent during the build-up of the function
2354/// candidates table once for each set of powers in
2355/// variables. Notice, that the argument array contains the powers
2356/// PLUS ONE. For example, to De select the function
2357/// f = x1^2 * x2^4 * x3^5,
2358/// this method should return kFALSE if given the argument
2359/// { 3, 4, 6 }
2360
2362{
2363 return kTRUE;
2364}
2365
2366////////////////////////////////////////////////////////////////////////////////
2367/// Set the max angle (in degrees) between the initial data vector to
2368/// be fitted, and the new candidate function to be included in the
2369/// fit. By default it is 0, which automatically chooses another
2370/// selection criteria. See also
2371/// <a href="#TMultiDimFit:description">class description</a>
2372
2374{
2375 if (ang >= 90 || ang < 0) {
2376 Warning("SetMaxAngle", "angle must be in [0,90)");
2377 return;
2378 }
2379
2380 fMaxAngle = ang;
2381}
2382
2383////////////////////////////////////////////////////////////////////////////////
2384/// Set the min angle (in degrees) between a new candidate function
2385/// and the subspace spanned by the previously accepted
2386/// functions. See also
2387/// <a href="#TMultiDimFit:description">class description</a>
2388
2390{
2391 if (ang > 90 || ang <= 0) {
2392 Warning("SetMinAngle", "angle must be in [0,90)");
2393 return;
2394 }
2395
2396 fMinAngle = ang;
2397
2398}
2399
2400
2401////////////////////////////////////////////////////////////////////////////////
2402/// Define a user function. The input array must be of the form
2403/// (p11, ..., p1N, ... ,pL1, ..., pLN)
2404/// Where N is the dimension of the data sample, L is the number of
2405/// terms (given in terms) and the first number, labels the term, the
2406/// second the variable. More information is given in the
2407/// <a href="#TMultiDimFit:description">class description</a>
2408
2409void TMultiDimFit::SetPowers(const Int_t* powers, Int_t terms)
2410{
2412 fMaxFunctions = terms;
2413 fMaxTerms = terms;
2414 fMaxStudy = terms;
2416 fPowers = new Int_t[fMaxFuncNV];
2417 Int_t i, j;
2418 for (i = 0; i < fMaxFunctions; i++)
2419 for(j = 0; j < fNVariables; j++)
2420 fPowers[i * fNVariables + j] = powers[i * fNVariables + j] + 1;
2421}
2422
2423////////////////////////////////////////////////////////////////////////////////
2424/// Set the user parameter for the function selection. The bigger the
2425/// limit, the more functions are used. The meaning of this variable
2426/// is defined in the
2427/// <a href="#TMultiDimFit:description">class description</a>
2428
2430{
2431 fPowerLimit = limit;
2432}
2433
2434////////////////////////////////////////////////////////////////////////////////
2435/// Set the maximum power to be considered in the fit for each
2436/// variable. See also
2437/// <a href="#TMultiDimFit:description">class description</a>
2438
2440{
2441 if (!powers)
2442 return;
2443
2444 for (Int_t i = 0; i < fNVariables; i++)
2445 fMaxPowers[i] = powers[i]+1;
2446}
2447
2448////////////////////////////////////////////////////////////////////////////////
2449/// Set the acceptable relative error for when sum of square
2450/// residuals is considered minimized. For a full account, refer to
2451/// the
2452/// <a href="#TMultiDimFit:description">class description</a>
2453
2455{
2456 fMinRelativeError = error;
2457}
2458
2459
2460////////////////////////////////////////////////////////////////////////////////
2461/// PRIVATE METHOD:
2462/// Test whether the currently considered function contributes to the
2463/// fit. See also
2464/// <a href="#TMultiDimFit:description">class description</a>
2465
2467 Double_t dResidur)
2468{
2469 if (fNCoefficients != 0) {
2470 // Now for the second test:
2471 if (fMaxAngle == 0) {
2472 // If the user hasn't supplied a max angle do the test as,
2473 if (dResidur <
2474 squareResidual / (fMaxTerms - fNCoefficients + 1 + 1E-10)) {
2475 return kFALSE;
2476 }
2477 }
2478 else {
2479 // If the user has provided a max angle, test if the calculated
2480 // angle is less then the max angle.
2481 if (TMath::Sqrt(dResidur/fSumSqAvgQuantity) <
2483 return kFALSE;
2484 }
2485 }
2486 }
2487 // If we get here, the function is OK
2488 return kTRUE;
2489}
2490
2491
2492////////////////////////////////////////////////////////////////////////////////
2493/// Helper function for doing the minimisation of Chi2 using Minuit
2494
2495void mdfHelper(int& /*npar*/, double* /*divs*/, double& chi2,
2496 double* coeffs, int /*flag*/)
2497{
2498 // Get pointer to current TMultiDimFit object.
2500 chi2 = mdf->MakeChi2(coeffs);
2501}
ROOT::R::TRInterface & r
Definition Object.C:4
#define b(i)
Definition RSha256.hxx:100
#define f(i)
Definition RSha256.hxx:104
#define h(i)
Definition RSha256.hxx:106
#define e(i)
Definition RSha256.hxx:103
int Int_t
Definition RtypesCore.h:45
char Char_t
Definition RtypesCore.h:37
const Bool_t kFALSE
Definition RtypesCore.h:92
double Double_t
Definition RtypesCore.h:59
const Bool_t kTRUE
Definition RtypesCore.h:91
const char Option_t
Definition RtypesCore.h:66
#define ClassImp(name)
Definition Rtypes.h:364
#define SETBIT(n, i)
Definition Rtypes.h:86
#define TESTBIT(n, i)
Definition Rtypes.h:88
int type
Definition TGX11.cxx:121
double sqrt(double)
TMatrixTRow< Double_t > TMatrixDRow
TMatrixTDiag_const< Double_t > TMatrixDDiag_const
#define HIST_DSHIF
#define PARAM_RELERR
#define HIST_DORIG
#define HIST_XNORM
#define HIST_RX
#define HIST_XORIG
#define PARAM_SEVERAL
#define PARAM_MAXSTUDY
static void mdfHelper(int &, double *, double &, double *, int)
Helper function for doing the minimisation of Chi2 using Minuit.
#define HIST_RTEST
#define DEGRAD
#define HIST_RTRAI
#define HIST_RD
#define PARAM_MAXTERMS
#define gROOT
Definition TROOT.h:406
char * Form(const char *fmt,...)
Using a TBrowser one can browse all ROOT objects.
Definition TBrowser.h:37
This class stores the date and time with a precision of one second in an unsigned 32 bit word (950130...
Definition TDatime.h:37
const char * AsString() const
Return the date & time as a string (ctime() format).
Definition TDatime.cxx:102
Cholesky Decomposition class.
Definition TDecompChol.h:25
Bool_t Invert(TMatrixDSym &inv)
For a symmetric matrix A(m,m), its inverse A_inv(m,m) is returned .
virtual Bool_t Decompose()
Matrix A is decomposed in component U so that A = U^T * U If the decomposition succeeds,...
1-D histogram with a double per channel (see TH1 documentation)}
Definition TH1.h:618
TH1 is the base class of all histogram classes in ROOT.
Definition TH1.h:58
virtual Bool_t Add(TF1 *h1, Double_t c1=1, Option_t *option="")
Performs the operation: this = this + c1*f1 if errors are defined (see TH1::Sumw2),...
Definition TH1.cxx:822
2-D histogram with a double per channel (see TH1 documentation)}
Definition TH2.h:292
A doubly linked list.
Definition TList.h:44
virtual void Add(TObject *obj)
Definition TList.h:87
virtual TObject * FindObject(const char *name) const
Find an object in this list using its name.
Definition TList.cxx:578
virtual void Clear(Option_t *option="")
Remove all objects from the list.
Definition TList.cxx:402
void Print(Option_t *name="") const
Print the matrix as a table of elements.
virtual TMatrixTBase< Element > & Zero()
Set matrix elements to zero.
Bool_t IsValid() const
virtual TMatrixTBase< Element > & NormByDiag(const TVectorT< Element > &v, Option_t *option="D")
option:
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...
Multidimensional Fits in ROOT.
virtual void MakeCorrelation()
PRIVATE METHOD: Compute the correlation matrix.
TVectorD fCoefficients
Double_t fPrecision
Byte_t fHistogramMask
virtual void MakeCoefficientErrors()
PRIVATE METHOD: Compute the errors on the coefficients.
TMatrixD fOrthCurvatureMatrix
Bool_t fShowCorrelation
Bool_t fIsUserFunction
Int_t fNVariables
Int_t fMaxResidualRow
Double_t fMinQuantity
virtual Double_t Eval(const Double_t *x, const Double_t *coeff=0) const
Evaluate parameterization at point x.
Double_t fSumSqQuantity
TVectorD fMaxVariables
Int_t fSampleSize
Double_t fSumSqAvgQuantity
Int_t * fPowers
virtual void MakeNormalized()
PRIVATE METHOD: Normalize data to the interval [-1;1].
Int_t * fMaxPowers
TVectorD fTestQuantity
TVectorD fTestSqError
Int_t fMaxFunctions
static TMultiDimFit * Instance()
Return the static instance.
virtual ~TMultiDimFit()
Destructor.
virtual void MakeCode(const char *functionName="MDF", Option_t *option="")
Generate the file <filename> with .C appended if argument doesn't end in .cxx or ....
virtual void MakeCandidates()
PRIVATE METHOD: Create list of candidate functions for the parameterisation.
Double_t fTestError
Double_t fMinResidual
void SetPowerLimit(Double_t limit=1e-3)
Set the user parameter for the function selection.
virtual void FindParameterization(Option_t *option="")
Find the parameterization.
virtual Double_t MakeGramSchmidt(Int_t function)
PRIVATE METHOD: Make Gram-Schmidt orthogonalisation.
TVectorD fSqError
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,...
TVectorD fQuantity
TVectorD fTestVariables
TVectorD fOrthCoefficients
Int_t fTestSampleSize
TVectorD fResiduals
virtual Bool_t TestFunction(Double_t squareResidual, Double_t dResidur)
PRIVATE METHOD: Test whether the currently considered function contributes to the fit.
TVectorD fVariables
Double_t fCorrelationCoeff
void SetMaxPowers(const Int_t *powers)
Set the maximum power to be considered in the fit for each variable.
Int_t * fMaxPowersFinal
Int_t fMinResidualRow
virtual void SetPowers(const Int_t *powers, Int_t terms)
Define a user function.
Double_t fMaxQuantity
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...
Double_t fRMS
TVirtualFitter * fFitter
TMatrixD fFunctions
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,...
Int_t fParameterisationCode
Double_t fMaxAngle
Int_t fNCoefficients
virtual void Fit(Option_t *option="")
Try to fit the found parameterisation to the test sample.
TMultiDimFit()
Empty CTOR. Do not use.
Double_t fTestCorrelationCoeff
Double_t fMinAngle
void SetMinRelativeError(Double_t error)
Set the acceptable relative error for when sum of square residuals is considered minimized.
Double_t fSumSqResidual
TVectorD fMeanVariables
virtual void MakeMethod(const Char_t *className="MDF", Option_t *option="")
Generate the file <classname>MDF.cxx which contains the implementation of the method:
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...
Double_t fTestPrecision
virtual void MakeCoefficients()
PRIVATE METHOD: Invert the model matrix B, and compute final coefficients.
TMatrixD fCorrelationMatrix
Int_t * fPowerIndex
TMatrixD fOrthFunctions
virtual void Clear(Option_t *option="")
Clear internal structures and variables.
virtual Double_t EvalControl(const Int_t *powers) const
PRIVATE METHOD: Calculate the control parameter from the passed powers.
virtual void Browse(TBrowser *b)
Browse the TMultiDimFit object in the TBrowser.
virtual void MakeParameterization()
PRIVATE METHOD: Find the parameterization over the training sample.
EMDFPolyType fPolyType
Fit object (MINUIT)
virtual Bool_t Select(const Int_t *iv)
Selection method.
virtual Double_t EvalError(const Double_t *x, const Double_t *coeff=0) const
Evaluate parameterization error at point x.
TList * fHistograms
Double_t fChi2
Double_t fMaxResidual
Double_t fError
virtual Double_t EvalFactor(Int_t p, Double_t x) const
PRIVATE METHOD: Evaluate function with power p at variable value x.
void SetMinAngle(Double_t angle=1)
Set the min angle (in degrees) between a new candidate function and the subspace spanned by the previ...
virtual void Print(Option_t *option="ps") const
Print statistics etc.
TVectorD fMinVariables
TVectorD fOrthFunctionNorms
Double_t fMeanQuantity
TVectorD fCoefficientsRMS
Double_t fPowerLimit
static TMultiDimFit * fgInstance
virtual void MakeHistograms(Option_t *option="A")
Make histograms of the result of the analysis.
Bool_t fIsVerbose
Double_t fMinRelativeError
virtual Double_t MakeChi2(const Double_t *coeff=0)
Calculate Chi square over either the test sample.
Int_t * fFunctionCodes
The TNamed class is the base class for all named ROOT classes.
Definition TNamed.h:29
virtual const char * GetName() const
Returns name of object.
Definition TNamed.h:47
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition TObject.cxx:879
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition TObject.cxx:893
Basic string class.
Definition TString.h:136
Ssiz_t Length() const
Definition TString.h:410
void ToLower()
Change string to lower-case.
Definition TString.cxx:1145
Bool_t EndsWith(const char *pat, ECaseCompare cmp=kExact) const
Return true if string ends with the specified string.
Definition TString.cxx:2197
const char * Data() const
Definition TString.h:369
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition TString.h:624
TVectorT< Element > & Zero()
Set vector elements to zero.
Definition TVectorT.cxx:453
TVectorT< Element > & ResizeTo(Int_t lwb, Int_t upb)
Resize the vector to [lwb:upb] .
Definition TVectorT.cxx:294
Int_t GetNrows() const
Definition TVectorT.h:75
Bool_t IsValid() const
Definition TVectorT.h:83
virtual void PrintResults(Int_t level, Double_t amin) const =0
virtual void SetFCN(void(*fcn)(Int_t &, Double_t *, Double_t &f, Double_t *, Int_t))
To set the address of the minimization objective function called by the native compiler (see function...
virtual Int_t SetParameter(Int_t ipar, const char *parname, Double_t value, Double_t verr, Double_t vlow, Double_t vhigh)=0
virtual Int_t ExecuteCommand(const char *command, Double_t *args, Int_t nargs)=0
virtual Double_t GetParameter(Int_t ipar) const =0
static TVirtualFitter * Fitter(TObject *obj, Int_t maxpar=25)
Static function returning a pointer to the current fitter.
Double_t y[n]
Definition legend1.C:17
Double_t x[n]
Definition legend1.C:17
const Int_t n
Definition legend1.C:16
Short_t Max(Short_t a, Short_t b)
Definition TMathBase.h:212
Double_t Sqrt(Double_t x)
Definition TMath.h:691
Double_t Cos(Double_t)
Definition TMath.h:643
Double_t Sin(Double_t)
Definition TMath.h:639
Short_t Abs(Short_t d)
Definition TMathBase.h:120
auto * m
Definition textangle.C:8
auto * l
Definition textangle.C:4
static uint64_t sum(uint64_t i)
Definition Factory.cxx:2345
REAL epsilon
Definition triangle.c:617