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