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
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 = nullptr;
450 fPowerLimit = 0;
451
452 fMaxFunctions = 0;
453 fFunctionCodes = nullptr;
454 fMaxStudy = 0;
455 fMaxFuncNV = 0;
456
457 fMaxPowersFinal = nullptr;
458 fPowers = nullptr;
459 fPowerIndex = nullptr;
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 = nullptr;
480 fHistogramMask = 0;
481 fBinVarX = 100;
482 fBinVarY = 100;
483
484 fFitter = nullptr;
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
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 = nullptr;
540 fMaxStudy = 0;
541 fMaxFuncNV = 0;
542
543 fMaxPowersFinal = new Int_t[dimension];
544 fPowers = nullptr;
545 fPowerIndex = nullptr;
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 = nullptr;
566 fHistogramMask = 0;
567 fBinVarX = 100;
568 fBinVarY = 100;
569
570 fFitter = nullptr;
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.
635 if (fSampleSize > size) {
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
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) {
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
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 = nullptr;
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
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.
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.
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, TString::Format("coeff%02d",i).Data(),
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
1084 // use big enough string buffer to get variable name which is not used
1085 char namebuf[512];
1086 fFitter->GetParameter(i, namebuf,
1087 val, err, low, high);
1088 (void) namebuf;
1089 fCoefficients(i) = val;
1090 fCoefficientsRMS(i) = err;
1091 }
1092 delete [] x;
1093 delete [] arglist;
1094}
1095
1096////////////////////////////////////////////////////////////////////////////////
1097/// Return the static instance.
1098
1100{
1101 return fgInstance;
1102}
1103
1104////////////////////////////////////////////////////////////////////////////////
1105/// PRIVATE METHOD:
1106/// Create list of candidate functions for the parameterisation. See
1107/// also
1108/// <a href="#TMultiDimFit:description">class description</a>
1109
1111{
1112 Int_t i = 0;
1113 Int_t j = 0;
1114 Int_t k = 0;
1115
1116 // The temporary array to store the powers in. We don't need to
1117 // initialize this array however.
1119 Int_t *powers = new Int_t[fMaxFuncNV];
1120
1121 // store of `control variables'
1122 Double_t* control = new Double_t[fMaxFunctions];
1123
1124 // We've better initialize the variables
1125 Int_t *iv = new Int_t[fNVariables];
1126 for (i = 0; i < fNVariables; i++)
1127 iv[i] = 1;
1128
1129 if (!fIsUserFunction) {
1130
1131 // Number of funcs selected
1132 Int_t numberFunctions = 0;
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
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(TString::Format("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++) {
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(TString::Format("x_%d_orig",i)))
1616 fHistograms->Add(new TH1D(TString::Format("x_%d_orig",i),
1617 TString::Format("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(TString::Format("x_%d_norm",i)))
1635 fHistograms->Add(new TH1D(TString::Format("x_%d_norm",i),
1636 TString::Format("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(TString::Format("res_x_%d",i)))
1654 fHistograms->Add(new TH2D(TString::Format("res_x_%d",i),
1655 TString::Format("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/// \code{.cpp}
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/// \endcode
1717///
1718/// are initialized, and assumed to exist. The class declaration is
1719/// assumed to be in `<classname>.h` and assumed to be provided by the
1720/// user.
1721///
1722/// \see TMultiDimFit::MakeRealCode for a list of options
1723///
1724/// The minimal class definition is:
1725/// \code{.cpp}
1726/// class <classname> {
1727/// public:
1728/// Int_t <classname>::fgNVariables; // Number of variables
1729/// Int_t <classname>::fgNCoefficients; // Number of terms
1730/// Double_t <classname>::fgDMean; // Mean from training sample
1731/// Double_t <classname>::fgXMean[]; // Mean from training sample
1732/// Double_t <classname>::fgXMin[]; // Min from training sample
1733/// Double_t <classname>::fgXMax[]; // Max from training sample
1734/// Double_t <classname>::fgCoefficient[]; // Coefficients
1735/// Int_t <classname>::fgPower[]; // Function powers
1736///
1737/// Double_t Eval(Double_t *x);
1738/// };
1739/// \endcode
1740///
1741/// Whether the method `<classname>::%Eval` should be static or not, is
1742/// up to the user.
1743
1745{
1746 MakeRealCode(TString::Format("%sMDF.cxx", classname), classname, option);
1747}
1748
1749
1750
1751////////////////////////////////////////////////////////////////////////////////
1752/// PRIVATE METHOD:
1753/// Normalize data to the interval [-1;1]. This is needed for the
1754/// classes method to work.
1755
1757{
1758 Int_t i = 0;
1759 Int_t j = 0;
1760 Int_t k = 0;
1761
1762 for (i = 0; i < fSampleSize; i++) {
1764 ((TH1D*)fHistograms->FindObject("d_orig"))->Fill(fQuantity(i));
1765
1768
1770 ((TH1D*)fHistograms->FindObject("d_shifted"))->Fill(fQuantity(i));
1771
1772 for (j = 0; j < fNVariables; j++) {
1773 Double_t range = 1. / (fMaxVariables(j) - fMinVariables(j));
1774 k = i * fNVariables + j;
1775
1776 // Fill histograms of original independent variables
1778 ((TH1D*)fHistograms->FindObject(TString::Format("x_%d_orig",j)))
1779 ->Fill(fVariables(k));
1780
1781 // Normalise independent variables
1782 fVariables(k) = 1 + 2 * range * (fVariables(k) - fMaxVariables(j));
1783
1784 // Fill histograms of normalised independent variables
1786 ((TH1D*)fHistograms->FindObject(TString::Format("x_%d_norm",j)))
1787 ->Fill(fVariables(k));
1788
1789 }
1790 }
1791 // Shift min and max of dependent variable
1794
1795 // Shift mean of independent variables
1796 for (i = 0; i < fNVariables; i++) {
1797 Double_t range = 1. / (fMaxVariables(i) - fMinVariables(i));
1798 fMeanVariables(i) = 1 + 2 * range * (fMeanVariables(i)
1799 - fMaxVariables(i));
1800 }
1801}
1802
1803
1804////////////////////////////////////////////////////////////////////////////////
1805/// PRIVATE METHOD:
1806/// Find the parameterization over the training sample. A full account
1807/// of the algorithm is given in the
1808/// <a href="#TMultiDimFit:description">class description</a>
1809
1811{
1812 Int_t i = -1;
1813 Int_t j = 0;
1814 Int_t k = 0;
1815 Int_t maxPass = 3;
1816 Int_t studied = 0;
1817 Double_t squareResidual = fSumSqAvgQuantity;
1818 fNCoefficients = 0;
1825 fFunctions = 1;
1826
1829 Int_t l;
1830 for (l=0;l<fMaxFunctions;l++) fFunctionCodes[l] = 0;
1831 for (l=0;l<fMaxTerms;l++) fPowerIndex[l] = 0;
1832
1833 if (fMaxAngle != 0) maxPass = 100;
1834 if (fIsUserFunction) maxPass = 1;
1835
1836 // Loop over the number of functions we want to study.
1837 // increment inspection counter
1838 while(kTRUE) {
1839
1840 // Reach user defined limit of studies
1841 if (studied++ >= fMaxStudy) {
1843 break;
1844 }
1845
1846 // Considered all functions several times
1847 if (k >= maxPass) {
1849 break;
1850 }
1851
1852 // increment function counter
1853 i++;
1854
1855 // If we've reached the end of the functions, restart pass
1856 if (i == fMaxFunctions) {
1857 if (fMaxAngle != 0)
1858 fMaxAngle += (90 - fMaxAngle) / 2;
1859 i = 0;
1860 studied--;
1861 k++;
1862 continue;
1863 }
1864 if (studied == 1)
1865 fFunctionCodes[i] = 0;
1866 else if (fFunctionCodes[i] >= 2)
1867 continue;
1868
1869 // Print a happy message
1870 if (fIsVerbose && studied == 1)
1871 std::cout << "Coeff SumSqRes Contrib Angle QM Func"
1872 << " Value W^2 Powers" << std::endl;
1873
1874 // Make the Gram-Schmidt
1875 Double_t dResidur = MakeGramSchmidt(i);
1876
1877 if (dResidur == 0) {
1878 // This function is no good!
1879 // First test is in MakeGramSchmidt
1880 fFunctionCodes[i] = 1;
1881 continue;
1882 }
1883
1884 // If user specified function, assume they know what they are doing
1885 if (!fIsUserFunction) {
1886 // Flag this function as considered
1887 fFunctionCodes[i] = 2;
1888
1889 // Test if this function contributes to the fit
1890 if (!TestFunction(squareResidual, dResidur)) {
1891 fFunctionCodes[i] = 1;
1892 continue;
1893 }
1894 }
1895
1896 // If we get to here, the function currently considered is
1897 // fNCoefficients, so we increment the counter
1898 // Flag this function as OK, and store and the number in the
1899 // index.
1900 fFunctionCodes[i] = 3;
1903
1904 // We add the current contribution to the sum of square of
1905 // residuals;
1906 squareResidual -= dResidur;
1907
1908
1909 // Calculate control parameter from this function
1910 for (j = 0; j < fNVariables; j++) {
1911 if (fNCoefficients == 1
1912 || fMaxPowersFinal[j] <= fPowers[i * fNVariables + j] - 1)
1913 fMaxPowersFinal[j] = fPowers[i * fNVariables + j] - 1;
1914 }
1916
1917 // Print the statistics about this function
1918 if (fIsVerbose) {
1919 std::cout << std::setw(5) << fNCoefficients << " "
1920 << std::setw(10) << std::setprecision(4) << squareResidual << " "
1921 << std::setw(10) << std::setprecision(4) << dResidur << " "
1922 << std::setw(7) << std::setprecision(3) << fMaxAngle << " "
1923 << std::setw(7) << std::setprecision(3) << s << " "
1924 << std::setw(5) << i << " "
1925 << std::setw(10) << std::setprecision(4)
1927 << std::setw(10) << std::setprecision(4)
1929 << std::flush;
1930 for (j = 0; j < fNVariables; j++)
1931 std::cout << " " << fPowers[i * fNVariables + j] - 1 << std::flush;
1932 std::cout << std::endl;
1933 }
1934
1935 if (fNCoefficients >= fMaxTerms /* && fIsVerbose */) {
1937 break;
1938 }
1939
1940 Double_t err = TMath::Sqrt(TMath::Max(1e-20,squareResidual) /
1942 if (err < fMinRelativeError) {
1944 break;
1945 }
1946
1947 }
1948
1949 fError = TMath::Max(1e-20,squareResidual);
1952}
1953
1954
1955////////////////////////////////////////////////////////////////////////////////
1956/// PRIVATE METHOD:
1957/// This is the method that actually generates the code for the
1958/// evaluation the parameterization on some point.
1959/// It's called by TMultiDimFit::MakeCode and TMultiDimFit::MakeMethod.
1960///
1961/// The options are: NONE so far
1962
1964 const char *classname,
1965 Option_t *)
1966{
1967 Int_t i, j;
1968
1969 Bool_t isMethod = (classname[0] == '\0' ? kFALSE : kTRUE);
1970 TString prefix;
1971 const char *cv_qual = isMethod ? "" : "static ";
1972 if (isMethod)
1973 prefix.Form("%s::", classname);
1974
1975 std::ofstream outFile(filename,std::ios::out|std::ios::trunc);
1976 if (!outFile) {
1977 Error("MakeRealCode","couldn't open output file '%s'",filename);
1978 return;
1979 }
1980
1981 if (fIsVerbose)
1982 std::cout << "Writing on file \"" << filename << "\" ... " << std::flush;
1983 //
1984 // Write header of file
1985 //
1986 // Emacs mode line ;-)
1987 outFile << "// -*- mode: c++ -*-" << std::endl;
1988 // Info about creator
1989 outFile << "// " << std::endl
1990 << "// File " << filename
1991 << " generated by TMultiDimFit::MakeRealCode" << std::endl;
1992 // Time stamp
1993 TDatime date;
1994 outFile << "// on " << date.AsString() << std::endl;
1995 // ROOT version info
1996 outFile << "// ROOT version " << gROOT->GetVersion()
1997 << std::endl << "//" << std::endl;
1998 // General information on the code
1999 outFile << "// This file contains the function " << std::endl
2000 << "//" << std::endl
2001 << "// double " << prefix << "MDF(double *x); " << std::endl
2002 << "//" << std::endl
2003 << "// For evaluating the parameterization obtained" << std::endl
2004 << "// from TMultiDimFit and the point x" << std::endl
2005 << "// " << std::endl
2006 << "// See TMultiDimFit class documentation for more "
2007 << "information " << std::endl << "// " << std::endl;
2008 // Header files
2009 if (isMethod)
2010 // If these are methods, we need the class header
2011 outFile << "#include \"" << classname << ".h\"" << std::endl;
2012
2013 //
2014 // Now for the data
2015 //
2016 outFile << "//" << std::endl
2017 << "// Static data variables" << std::endl
2018 << "//" << std::endl;
2019 outFile << cv_qual << "int " << prefix << "gNVariables = "
2020 << fNVariables << ";" << std::endl;
2021 outFile << cv_qual << "int " << prefix << "gNCoefficients = "
2022 << fNCoefficients << ";" << std::endl;
2023 outFile << cv_qual << "double " << prefix << "gDMean = "
2024 << fMeanQuantity << ";" << std::endl;
2025
2026 // Assignment to mean vector.
2027 outFile << "// Assignment to mean vector." << std::endl;
2028 outFile << cv_qual << "double " << prefix
2029 << "gXMean[] = {" << std::endl;
2030 for (i = 0; i < fNVariables; i++)
2031 outFile << (i != 0 ? ", " : " ") << fMeanVariables(i) << std::flush;
2032 outFile << " };" << std::endl << std::endl;
2033
2034 // Assignment to minimum vector.
2035 outFile << "// Assignment to minimum vector." << std::endl;
2036 outFile << cv_qual << "double " << prefix
2037 << "gXMin[] = {" << std::endl;
2038 for (i = 0; i < fNVariables; i++)
2039 outFile << (i != 0 ? ", " : " ") << fMinVariables(i) << std::flush;
2040 outFile << " };" << std::endl << std::endl;
2041
2042 // Assignment to maximum vector.
2043 outFile << "// Assignment to maximum vector." << std::endl;
2044 outFile << cv_qual << "double " << prefix
2045 << "gXMax[] = {" << std::endl;
2046 for (i = 0; i < fNVariables; i++)
2047 outFile << (i != 0 ? ", " : " ") << fMaxVariables(i) << std::flush;
2048 outFile << " };" << std::endl << std::endl;
2049
2050 // Assignment to coefficients vector.
2051 outFile << "// Assignment to coefficients vector." << std::endl;
2052 outFile << cv_qual << "double " << prefix
2053 << "gCoefficient[] = {" << std::flush;
2054 for (i = 0; i < fNCoefficients; i++)
2055 outFile << (i != 0 ? "," : "") << std::endl
2056 << " " << fCoefficients(i) << std::flush;
2057 outFile << std::endl << " };" << std::endl << std::endl;
2058
2059 // Assignment to error coefficients vector.
2060 outFile << "// Assignment to error coefficients vector." << std::endl;
2061 outFile << cv_qual << "double " << prefix
2062 << "gCoefficientRMS[] = {" << std::flush;
2063 for (i = 0; i < fNCoefficients; i++)
2064 outFile << (i != 0 ? "," : "") << std::endl
2065 << " " << fCoefficientsRMS(i) << std::flush;
2066 outFile << std::endl << " };" << std::endl << std::endl;
2067
2068 // Assignment to powers vector.
2069 outFile << "// Assignment to powers vector." << std::endl
2070 << "// The powers are stored row-wise, that is" << std::endl
2071 << "// p_ij = " << prefix
2072 << "gPower[i * NVariables + j];" << std::endl;
2073 outFile << cv_qual << "int " << prefix
2074 << "gPower[] = {" << std::flush;
2075 for (i = 0; i < fNCoefficients; i++) {
2076 for (j = 0; j < fNVariables; j++) {
2077 if (j != 0) outFile << std::flush << " ";
2078 else outFile << std::endl << " ";
2079 outFile << fPowers[fPowerIndex[i] * fNVariables + j]
2080 << (i == fNCoefficients - 1 && j == fNVariables - 1 ? "" : ",")
2081 << std::flush;
2082 }
2083 }
2084 outFile << std::endl << "};" << std::endl << std::endl;
2085
2086
2087 //
2088 // Finally we reach the function itself
2089 //
2090 outFile << "// " << std::endl
2091 << "// The "
2092 << (isMethod ? "method " : "function ")
2093 << " double " << prefix
2094 << "MDF(double *x)"
2095 << std::endl << "// " << std::endl;
2096 outFile << "double " << prefix
2097 << "MDF(double *x) {" << std::endl
2098 << " double returnValue = " << prefix << "gDMean;" << std::endl
2099 << " int i = 0, j = 0, k = 0;" << std::endl
2100 << " for (i = 0; i < " << prefix << "gNCoefficients ; i++) {"
2101 << std::endl
2102 << " // Evaluate the ith term in the expansion" << std::endl
2103 << " double term = " << prefix << "gCoefficient[i];"
2104 << std::endl
2105 << " for (j = 0; j < " << prefix << "gNVariables; j++) {"
2106 << std::endl
2107 << " // Evaluate the polynomial in the jth variable." << std::endl
2108 << " int power = "<< prefix << "gPower["
2109 << prefix << "gNVariables * i + j]; " << std::endl
2110 << " double p1 = 1, p2 = 0, p3 = 0, r = 0;" << std::endl
2111 << " double v = 1 + 2. / ("
2112 << prefix << "gXMax[j] - " << prefix
2113 << "gXMin[j]) * (x[j] - " << prefix << "gXMax[j]);" << std::endl
2114 << " // what is the power to use!" << std::endl
2115 << " switch(power) {" << std::endl
2116 << " case 1: r = 1; break; " << std::endl
2117 << " case 2: r = v; break; " << std::endl
2118 << " default: " << std::endl
2119 << " p2 = v; " << std::endl
2120 << " for (k = 3; k <= power; k++) { " << std::endl
2121 << " p3 = p2 * v;" << std::endl;
2122 if (fPolyType == kLegendre)
2123 outFile << " p3 = ((2 * k - 3) * p2 * v - (k - 2) * p1)"
2124 << " / (k - 1);" << std::endl;
2125 if (fPolyType == kChebyshev)
2126 outFile << " p3 = 2 * v * p2 - p1; " << std::endl;
2127 outFile << " p1 = p2; p2 = p3; " << std::endl << " }" << std::endl
2128 << " r = p3;" << std::endl << " }" << std::endl
2129 << " // multiply this term by the poly in the jth var" << std::endl
2130 << " term *= r; " << std::endl << " }" << std::endl
2131 << " // Add this term to the final result" << std::endl
2132 << " returnValue += term;" << std::endl << " }" << std::endl
2133 << " return returnValue;" << std::endl << "}" << std::endl << std::endl;
2134
2135 // EOF
2136 outFile << "// EOF for " << filename << std::endl;
2137
2138 // Close the file
2139 outFile.close();
2140
2141 if (fIsVerbose)
2142 std::cout << "done" << std::endl;
2143}
2144
2145
2146////////////////////////////////////////////////////////////////////////////////
2147/// Print statistics etc.
2148/// Options are
2149/// P Parameters
2150/// S Statistics
2151/// C Coefficients
2152/// R Result of parameterisation
2153/// F Result of fit
2154/// K Correlation Matrix
2155/// M Pretty print formula
2156///
2157
2159{
2160 Int_t i = 0;
2161 Int_t j = 0;
2162
2163 TString opt(option);
2164 opt.ToLower();
2165
2166 if (opt.Contains("p")) {
2167 // Print basic parameters for this object
2168 std::cout << "User parameters:" << std::endl
2169 << "----------------" << std::endl
2170 << " Variables: " << fNVariables << std::endl
2171 << " Data points: " << fSampleSize << std::endl
2172 << " Max Terms: " << fMaxTerms << std::endl
2173 << " Power Limit Parameter: " << fPowerLimit << std::endl
2174 << " Max functions: " << fMaxFunctions << std::endl
2175 << " Max functions to study: " << fMaxStudy << std::endl
2176 << " Max angle (optional): " << fMaxAngle << std::endl
2177 << " Min angle: " << fMinAngle << std::endl
2178 << " Relative Error accepted: " << fMinRelativeError << std::endl
2179 << " Maximum Powers: " << std::flush;
2180 for (i = 0; i < fNVariables; i++)
2181 std::cout << " " << fMaxPowers[i] - 1 << std::flush;
2182 std::cout << std::endl << std::endl
2183 << " Parameterisation will be done using " << std::flush;
2184 if (fPolyType == kChebyshev)
2185 std::cout << "Chebyshev polynomials" << std::endl;
2186 else if (fPolyType == kLegendre)
2187 std::cout << "Legendre polynomials" << std::endl;
2188 else
2189 std::cout << "Monomials" << std::endl;
2190 std::cout << std::endl;
2191 }
2192
2193 if (opt.Contains("s")) {
2194 // Print statistics for read data
2195 std::cout << "Sample statistics:" << std::endl
2196 << "------------------" << std::endl
2197 << " D" << std::flush;
2198 for (i = 0; i < fNVariables; i++)
2199 std::cout << " " << std::setw(10) << i+1 << std::flush;
2200 std::cout << std::endl << " Max: " << std::setw(10) << std::setprecision(7)
2201 << fMaxQuantity << std::flush;
2202 for (i = 0; i < fNVariables; i++)
2203 std::cout << " " << std::setw(10) << std::setprecision(4)
2204 << fMaxVariables(i) << std::flush;
2205 std::cout << std::endl << " Min: " << std::setw(10) << std::setprecision(7)
2206 << fMinQuantity << std::flush;
2207 for (i = 0; i < fNVariables; i++)
2208 std::cout << " " << std::setw(10) << std::setprecision(4)
2209 << fMinVariables(i) << std::flush;
2210 std::cout << std::endl << " Mean: " << std::setw(10) << std::setprecision(7)
2211 << fMeanQuantity << std::flush;
2212 for (i = 0; i < fNVariables; i++)
2213 std::cout << " " << std::setw(10) << std::setprecision(4)
2214 << fMeanVariables(i) << std::flush;
2215 std::cout << std::endl << " Function Sum Squares: " << fSumSqQuantity
2216 << std::endl << std::endl;
2217 }
2218
2219 if (opt.Contains("r")) {
2220 std::cout << "Results of Parameterisation:" << std::endl
2221 << "----------------------------" << std::endl
2222 << " Total reduction of square residuals "
2223 << fSumSqResidual << std::endl
2224 << " Relative precision obtained: "
2225 << fPrecision << std::endl
2226 << " Error obtained: "
2227 << fError << std::endl
2228 << " Multiple correlation coefficient: "
2229 << fCorrelationCoeff << std::endl
2230 << " Reduced Chi square over sample: "
2231 << fChi2 / (fSampleSize - fNCoefficients) << std::endl
2232 << " Maximum residual value: "
2233 << fMaxResidual << std::endl
2234 << " Minimum residual value: "
2235 << fMinResidual << std::endl
2236 << " Estimated root mean square: "
2237 << fRMS << std::endl
2238 << " Maximum powers used: " << std::flush;
2239 for (j = 0; j < fNVariables; j++)
2240 std::cout << fMaxPowersFinal[j] << " " << std::flush;
2241 std::cout << std::endl
2242 << " Function codes of candidate functions." << std::endl
2243 << " 1: considered,"
2244 << " 2: too little contribution,"
2245 << " 3: accepted." << std::flush;
2246 for (i = 0; i < fMaxFunctions; i++) {
2247 if (i % 60 == 0)
2248 std::cout << std::endl << " " << std::flush;
2249 else if (i % 10 == 0)
2250 std::cout << " " << std::flush;
2251 std::cout << fFunctionCodes[i];
2252 }
2253 std::cout << std::endl << " Loop over candidates stopped because " << std::flush;
2254 switch(fParameterisationCode){
2255 case PARAM_MAXSTUDY:
2256 std::cout << "max allowed studies reached" << std::endl; break;
2257 case PARAM_SEVERAL:
2258 std::cout << "all candidates considered several times" << std::endl; break;
2259 case PARAM_RELERR:
2260 std::cout << "wanted relative error obtained" << std::endl; break;
2261 case PARAM_MAXTERMS:
2262 std::cout << "max number of terms reached" << std::endl; break;
2263 default:
2264 std::cout << "some unknown reason" << std::endl;
2265 break;
2266 }
2267 std::cout << std::endl;
2268 }
2269
2270 if (opt.Contains("f")) {
2271 std::cout << "Results of Fit:" << std::endl
2272 << "---------------" << std::endl
2273 << " Test sample size: "
2274 << fTestSampleSize << std::endl
2275 << " Multiple correlation coefficient: "
2276 << fTestCorrelationCoeff << std::endl
2277 << " Relative precision obtained: "
2278 << fTestPrecision << std::endl
2279 << " Error obtained: "
2280 << fTestError << std::endl
2281 << " Reduced Chi square over sample: "
2282 << fChi2 / (fSampleSize - fNCoefficients) << std::endl
2283 << std::endl;
2284 if (fFitter) {
2285 fFitter->PrintResults(1,1);
2286 std::cout << std::endl;
2287 }
2288 }
2289
2290 if (opt.Contains("c")){
2291 std::cout << "Coefficients:" << std::endl
2292 << "-------------" << std::endl
2293 << " # Value Error Powers" << std::endl
2294 << " ---------------------------------------" << std::endl;
2295 for (i = 0; i < fNCoefficients; i++) {
2296 std::cout << " " << std::setw(3) << i << " "
2297 << std::setw(12) << fCoefficients(i) << " "
2298 << std::setw(12) << fCoefficientsRMS(i) << " " << std::flush;
2299 for (j = 0; j < fNVariables; j++)
2300 std::cout << " " << std::setw(3)
2301 << fPowers[fPowerIndex[i] * fNVariables + j] - 1 << std::flush;
2302 std::cout << std::endl;
2303 }
2304 std::cout << std::endl;
2305 }
2306 if (opt.Contains("k") && fCorrelationMatrix.IsValid()) {
2307 std::cout << "Correlation Matrix:" << std::endl
2308 << "-------------------";
2310 }
2311
2312 if (opt.Contains("m")) {
2313 std::cout << "Parameterization:" << std::endl
2314 << "-----------------" << std::endl
2315 << " Normalised variables: " << std::endl;
2316 for (i = 0; i < fNVariables; i++)
2317 std::cout << "\ty_" << i << "\t= 1 + 2 * (x_" << i << " - "
2318 << fMaxVariables(i) << ") / ("
2319 << fMaxVariables(i) << " - " << fMinVariables(i) << ")"
2320 << std::endl;
2321 std::cout << std::endl
2322 << " f(";
2323 for (i = 0; i < fNVariables; i++) {
2324 std::cout << "y_" << i;
2325 if (i != fNVariables-1) std::cout << ", ";
2326 }
2327 std::cout << ") = ";
2328 for (i = 0; i < fNCoefficients; i++) {
2329 if (i != 0)
2330 std::cout << std::endl << "\t" << (fCoefficients(i) < 0 ? "- " : "+ ")
2332 else
2333 std::cout << fCoefficients(i);
2334 for (j = 0; j < fNVariables; j++) {
2336 switch (p) {
2337 case 1: break;
2338 case 2: std::cout << " * y_" << j; break;
2339 default:
2340 switch(fPolyType) {
2341 case kLegendre: std::cout << " * L_" << p-1 << "(y_" << j << ")"; break;
2342 case kChebyshev: std::cout << " * C_" << p-1 << "(y_" << j << ")"; break;
2343 default: std::cout << " * y_" << j << "^" << p-1; break;
2344 }
2345 }
2346
2347 }
2348 }
2349 std::cout << std::endl;
2350 }
2351}
2352
2353
2354////////////////////////////////////////////////////////////////////////////////
2355/// Selection method. User can override this method for specialized
2356/// selection of acceptable functions in fit. Default is to select
2357/// all. This message is sent during the build-up of the function
2358/// candidates table once for each set of powers in
2359/// variables. Notice, that the argument array contains the powers
2360/// PLUS ONE. For example, to De select the function
2361/// f = x1^2 * x2^4 * x3^5,
2362/// this method should return kFALSE if given the argument
2363/// { 3, 4, 6 }
2364
2366{
2367 return kTRUE;
2368}
2369
2370////////////////////////////////////////////////////////////////////////////////
2371/// Set the max angle (in degrees) between the initial data vector to
2372/// be fitted, and the new candidate function to be included in the
2373/// fit. By default it is 0, which automatically chooses another
2374/// selection criteria. See also
2375/// <a href="#TMultiDimFit:description">class description</a>
2376
2378{
2379 if (ang >= 90 || ang < 0) {
2380 Warning("SetMaxAngle", "angle must be in [0,90)");
2381 return;
2382 }
2383
2384 fMaxAngle = ang;
2385}
2386
2387////////////////////////////////////////////////////////////////////////////////
2388/// Set the min angle (in degrees) between a new candidate function
2389/// and the subspace spanned by the previously accepted
2390/// functions. See also
2391/// <a href="#TMultiDimFit:description">class description</a>
2392
2394{
2395 if (ang > 90 || ang <= 0) {
2396 Warning("SetMinAngle", "angle must be in [0,90)");
2397 return;
2398 }
2399
2400 fMinAngle = ang;
2401
2402}
2403
2404
2405////////////////////////////////////////////////////////////////////////////////
2406/// Define a user function. The input array must be of the form
2407/// (p11, ..., p1N, ... ,pL1, ..., pLN)
2408/// Where N is the dimension of the data sample, L is the number of
2409/// terms (given in terms) and the first number, labels the term, the
2410/// second the variable. More information is given in the
2411/// <a href="#TMultiDimFit:description">class description</a>
2412
2413void TMultiDimFit::SetPowers(const Int_t* powers, Int_t terms)
2414{
2416 fMaxFunctions = terms;
2417 fMaxTerms = terms;
2418 fMaxStudy = terms;
2420 fPowers = new Int_t[fMaxFuncNV];
2421 Int_t i, j;
2422 for (i = 0; i < fMaxFunctions; i++)
2423 for(j = 0; j < fNVariables; j++)
2424 fPowers[i * fNVariables + j] = powers[i * fNVariables + j] + 1;
2425}
2426
2427////////////////////////////////////////////////////////////////////////////////
2428/// Set the user parameter for the function selection. The bigger the
2429/// limit, the more functions are used. The meaning of this variable
2430/// is defined in the
2431/// <a href="#TMultiDimFit:description">class description</a>
2432
2434{
2435 fPowerLimit = limit;
2436}
2437
2438////////////////////////////////////////////////////////////////////////////////
2439/// Set the maximum power to be considered in the fit for each
2440/// variable. See also
2441/// <a href="#TMultiDimFit:description">class description</a>
2442
2444{
2445 if (!powers)
2446 return;
2447
2448 for (Int_t i = 0; i < fNVariables; i++)
2449 fMaxPowers[i] = powers[i]+1;
2450}
2451
2452////////////////////////////////////////////////////////////////////////////////
2453/// Set the acceptable relative error for when sum of square
2454/// residuals is considered minimized. For a full account, refer to
2455/// the
2456/// <a href="#TMultiDimFit:description">class description</a>
2457
2459{
2460 fMinRelativeError = error;
2461}
2462
2463
2464////////////////////////////////////////////////////////////////////////////////
2465/// PRIVATE METHOD:
2466/// Test whether the currently considered function contributes to the
2467/// fit. See also
2468/// <a href="#TMultiDimFit:description">class description</a>
2469
2471 Double_t dResidur)
2472{
2473 if (fNCoefficients != 0) {
2474 // Now for the second test:
2475 if (fMaxAngle == 0) {
2476 // If the user hasn't supplied a max angle do the test as,
2477 if (dResidur <
2478 squareResidual / (fMaxTerms - fNCoefficients + 1 + 1E-10)) {
2479 return kFALSE;
2480 }
2481 }
2482 else {
2483 // If the user has provided a max angle, test if the calculated
2484 // angle is less then the max angle.
2485 if (TMath::Sqrt(dResidur/fSumSqAvgQuantity) <
2487 return kFALSE;
2488 }
2489 }
2490 }
2491 // If we get here, the function is OK
2492 return kTRUE;
2493}
2494
2495
2496////////////////////////////////////////////////////////////////////////////////
2497/// Helper function for doing the minimisation of Chi2 using Minuit
2498
2499void mdfHelper(int& /*npar*/, double* /*divs*/, double& chi2,
2500 double* coeffs, int /*flag*/)
2501{
2502 // Get pointer to current TMultiDimFit object.
2504 chi2 = mdf->MakeChi2(coeffs);
2505}
#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
RooAbsReal & function()
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
int Int_t
Definition RtypesCore.h:45
char Char_t
Definition RtypesCore.h:37
constexpr Bool_t kFALSE
Definition RtypesCore.h:101
double Double_t
Definition RtypesCore.h:59
constexpr Bool_t kTRUE
Definition RtypesCore.h:100
const char Option_t
Definition RtypesCore.h:66
#define ClassImp(name)
Definition Rtypes.h:377
#define SETBIT(n, i)
Definition Rtypes.h:86
#define TESTBIT(n, i)
Definition Rtypes.h:88
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
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
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 Decompose() override
Matrix A is decomposed in component U so that A = U^T * U If the decomposition succeeds,...
Bool_t Invert(TMatrixDSym &inv)
For a symmetric matrix A(m,m), its inverse A_inv(m,m) is returned .
1-D histogram with a double per channel (see TH1 documentation)
Definition TH1.h:664
TH1 is the base class of all histogram classes in ROOT.
Definition TH1.h:59
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:824
2-D histogram with a double per channel (see TH1 documentation)
Definition TH2.h:338
A doubly linked list.
Definition TList.h:38
void Clear(Option_t *option="") override
Remove all objects from the list.
Definition TList.cxx:400
TObject * FindObject(const char *name) const override
Find an object in this list using its name.
Definition TList.cxx:576
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
virtual TMatrixTBase< Element > & NormByDiag(const TVectorT< Element > &v, Option_t *option="D")
option:
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:47
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition TObject.cxx:962
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition TObject.cxx:976
Basic string class.
Definition TString.h:139
Ssiz_t Length() const
Definition TString.h:417
void ToLower()
Change string to lower-case.
Definition TString.cxx:1182
Bool_t EndsWith(const char *pat, ECaseCompare cmp=kExact) const
Return true if string ends with the specified string.
Definition TString.cxx:2244
const char * Data() const
Definition TString.h:376
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:2378
void Form(const char *fmt,...)
Formats a string using a printf style format descriptor.
Definition TString.cxx:2356
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition TString.h:632
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:73
Bool_t IsValid() const
Definition TVectorT.h:81
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:250
Double_t Sqrt(Double_t x)
Returns the square root of x.
Definition TMath.h:662
Double_t Cos(Double_t)
Returns the cosine of an angle of x radians.
Definition TMath.h:594
Double_t Sin(Double_t)
Returns the sine of an angle of x radians.
Definition TMath.h:588
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition TMathBase.h:123
TMarker m
Definition textangle.C:8
TLine l
Definition textangle.C:4
static uint64_t sum(uint64_t i)
Definition Factory.cxx:2345