    // @(#)root/hist:$Name:  $:$Id: TMultiDimFit.cxx,v 1.11 2004/06/01 11:15:57 brun Exp $
    // Author: Christian Holm Christensen 07/11/2000

    //____________________________________________________________________
    //
    //
    //

/*

# Multidimensional Fits in ROOT

#
Overview

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

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

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

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

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

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

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

#
The Method

Let by the dependent quantity of interest, which depends smoothly on the observable quantities , which we'll denote by . Given a training sample of tuples of the form, (`TMultiDimFit::AddRow`)

where are independent variables, is the known, quantity dependent at , and is the square error in , the class `TMultiDimFit` will try to find the parameterization

| (1)

such that

| (2)

is minimal. Here are monomials, or Chebyshev or Legendre polynomials, labelled , in each variable , .

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

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

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

##
Limiting the Number of Terms

As always when dealing with fits, there's a real chance of over fitting. As is well-known, it's always possible to fit an polynomial in to points with , but the polynomial is not likely to fit new data at all [1]. Therefore, the user is asked to provide an upper limit, to the number of terms in (`TMultiDimFit::SetMaxTerms`).

However, since there's an infinite number of to choose from, the user is asked to give the maximum power. , of each variable to be considered in the minimization of (`TMultiDimFit::SetMaxPowers`).

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

There are still a huge amount of possible choices for ; in fact there are possible choices. Obviously we need to limit this. To this end, the user is asked to set a power control limit, (`TMultiDimFit::SetPowerLimit`), and a function is only accepted if

where is the leading power of variable in function . (`TMultiDimFit::MakeCandidates`). So the number of functions increase with (1, 2 is fine, 5 is way out).

## Gram-Schmidt Orthogonalisation

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

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

So we are let to consider an matrix , an element of which is given by

with | (3)

where labels the rows in the training sample and labels functions of variables, and . That is, is the term (or function) numbered evaluated at the data point . We have to normalise to for this to succeed [5] (`TMultiDimFit::MakeNormalized`). We then define a matrix of which the columns are given by

|  |  | (4)
|  |  | (5)

and is the component of orthogonal to . Hence we obtain [3],

if | (6)

We now take as a new model . We thus want to minimize

| (7)

where is a vector of the dependent quantity in the sample. Differentiation with respect to gives, using (6),

| (8)

or

| (9)

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

| (10)

Using (9), we see that

|  |  | (11)

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

##
Function Selection Based on Residual

Supposing that steps of the procedure have been performed, the problem now is to consider the function.

The sum of squares of residuals can be written as

| (12)

where the relation (9) have been taken into account. The contribution of the function to the reduction of S, is given by

| (13)

Two test are now applied to decide whether this function is to be included in the final expression, or not.

###
Test 1

Denoting by the subspace spanned by the function is by construction (see (4)) the projection of the function onto the direction perpendicular to . Now, if the length of (given by ) is very small compared to the length of this new function can not contribute much to the reduction of the sum of squares of residuals. The test consists then in calculating the angle between the two vectors and (see also figure 1) and requiring that it's greater then a threshold value which the user must set (`TMultiDimFit::SetMinAngle`).

Figure 1: (a) Angle between and , (b) angle between and

###
Test 2

Let be the data vector to be fitted. As illustrated in figure 1, the function will contribute significantly to the reduction of , if the angle between and is smaller than an upper limit , defined by the user (`TMultiDimFit::SetMaxAngle`)

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

In case isn't defined, an alternative method of performing this second test is used: The function is accepted if (refer also to equation (13))

| (14)

where is the sum of the first residuals from the functions previously accepted; and is the total number of functions allowed in the final expression of the fit (defined by user).

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

The more coefficients we evaluate, the more the sum of squares of residuals will be reduced. We can evaluate before inverting as shown below.

## Coefficients and Coefficient Errors

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

| (15)

where

| (16)

Consequently, is an upper triangle matrix, which can be readily inverted. So we now evaluate

| (17)

The model can therefore be written as

The original model is therefore identical with this if

| (18)

The reason we use rather then is to save storage, since can be stored in the same matrix as (`TMultiDimFit::MakeCoefficients`). The errors in the coefficients is calculated by inverting the curvature matrix of the non-orthogonal functions [1] (`TMultiDimFit::MakeCoefficientErrors`).

##
Considerations

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

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

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

1. Define are the 5 dependent quantities that define a track.
2. Compute, for different values of , the tracks through the magnetic field, and determine the corresponding .
3. Use the simulated observations to determine, with a simple approximation, the values of . We call these values .
4. Determine from a set of at least five relevant coordinates , using contrains, or alternative:
5. Perform a Principal Component Analysis (using `TPrincipal`), and use to get a linear transformation , so that are constrained and linear independent.
6. Perform a Principal Component Analysis on , to get linear indenpendent (among themselves, but not independent of ) quantities
7. For each component make a mutlidimensional fit, using as the variables, thus determing a set of coefficents .

To process data, using this parameterisation, do

1. Test wether the observation within the domain of the parameterization, using the result from the Principal Component Analysis.
2. Determine as before.
3. Detetmine as before.
4. Use the result of the fit to determind .
5. Transform back to from , using the result from the Principal Component Analysis.

##
Testing the parameterization

The class also provides functionality for testing the, over the training sample, found parameterization (`TMultiDimFit::Fit`). This is done by passing the class a test sample of tuples of the form , where are the independent variables, the known, dependent quantity, and is the square error in (`TMultiDimFit::AddTestRow`).

The parameterization is then evaluated at every in the test sample, and

is evaluated. The relative error over the test sample

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

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

Christian Holm
November 2000, NBI

## Bibliography

1
  Philip R. Bevington and D. Keith Robinson.
  Data Reduction and Error Analysis for the Physical Sciences.
  McGraw-Hill, 2 edition, 1992.

2
  René Brun et al.
  Mudifi.
  Long writeup DD/75-23, CERN, 1980.

3
  Gene H. Golub and Charles F. van Loan.
  Matrix Computations.
  John Hopkins Univeristy Press, Baltimore, 3 edition, 1996.

4
  F. James.
  Minuit.
  Long writeup D506, CERN, 1998.

5
  H. Wind.
  Function parameterization.
  In Proceedings of the 1972 CERN Computing and Data Processing School, volume 72-21 of Yellow report. CERN, 1972.

6
  H. Wind.
  1. principal component analysis, 2. pattern recognition for track finding, 3. interpolation and functional representation.
  Yellow report EP/81-12, CERN, 1981.

*/ //

// #include "Riostream.h" #include "TMultiDimFit.h" #include "TMath.h" #include "TH1.h" #include "TH2.h" #include "TROOT.h" #include "TBrowser.h" #define RADDEG (180. / TMath::Pi()) #define DEGRAD (TMath::Pi() / 180.) #define HIST_XORIG 0 #define HIST_DORIG 1 #define HIST_XNORM 2 #define HIST_DSHIF 3 #define HIST_RX 4 #define HIST_RD 5 #define HIST_RTRAI 6 #define HIST_RTEST 7 #define PARAM_MAXSTUDY 1 #define PARAM_SEVERAL 2 #define PARAM_RELERR 3 #define PARAM_MAXTERMS 4 //____________________________________________________________________ static void mdfHelper(int&, double*, double&, double*, int); //____________________________________________________________________ ClassImp(TMultiDimFit); //____________________________________________________________________ // Static instance. Used with mdfHelper and TMinuit TMultiDimFit* TMultiDimFit::fgInstance = 0; //____________________________________________________________________ TMultiDimFit::TMultiDimFit() { // Empty CTOR. Do not use fMeanQuantity = 0; fMaxQuantity = 0; fMinQuantity = 0; fSumSqQuantity = 0; fSumSqAvgQuantity = 0; fPowerLimit = 1; fMaxAngle = 0; fMinAngle = 1; fNVariables = 0; fMaxVariables = 0; fMinVariables = 0; fSampleSize = 0; fMaxAngle = 0; fMinAngle = 0; fPolyType = kMonomials; fShowCorrelation = kFALSE; fIsUserFunction = kFALSE; fPowers = 0; fMaxPowers = 0; fMaxPowersFinal = 0; fHistograms = 0; fHistogramMask = 0; fQuantity(0); fVariables(0); fMaxVariables(0); fMinVariables(0); fMeanVariables(0); fFitter = 0; } //____________________________________________________________________ TMultiDimFit::TMultiDimFit(Int_t dimension, EMDFPolyType type, Option_t *option) : TNamed("multidimfit","Multi-dimensional fit object"), fQuantity(dimension), fSqError(dimension), fVariables(dimension*100), fMeanVariables(dimension), fMaxVariables(dimension), fMinVariables(dimension) { // Constructor // Second argument is the type of polynomials to use in // parameterisation, one of: // TMultiDimFit::kMonomials // TMultiDimFit::kChebyshev // TMultiDimFit::kLegendre // // Options: // K Compute (k)correlation matrix // V Be verbose // // Default is no options. // fgInstance = this; fMeanQuantity = 0; fMaxQuantity = 0; fMinQuantity = 0; fSumSqQuantity = 0; fSumSqAvgQuantity = 0; fPowerLimit = 1; fMaxAngle = 0; fMinAngle = 1; fNVariables = dimension; fMaxVariables = 0; fMinVariables = 0; fSampleSize = 0; fTestSampleSize = 0; fMinRelativeError = 0.01; fError = 0; fTestError = 0; fPrecision = 0; fTestPrecision = 0; fParameterisationCode = 0; fPolyType = type; fShowCorrelation = kFALSE; fIsVerbose = kFALSE; TString opt = option; opt.ToLower(); if (opt.Contains("k")) fShowCorrelation = kTRUE; if (opt.Contains("v")) fIsVerbose = kTRUE; fIsUserFunction = kFALSE; fHistograms = 0; fHistogramMask = 0; fPowers = 0; fMaxPowers = new Int_t[dimension]; fMaxPowersFinal = new Int_t[dimension]; fFitter = 0; } //____________________________________________________________________ TMultiDimFit::~TMultiDimFit() { // DTOR delete [] fPowers; delete [] fMaxPowers; delete [] fMaxPowersFinal; delete [] fPowerIndex; delete [] fFunctionCodes; fHistograms->Clear("nodelete"); delete fHistograms; } //____________________________________________________________________ void TMultiDimFit::AddRow(const Double_t *x, Double_t D, Double_t E) { // Add a row consisting of fNVariables independent variables, the // known, dependent quantity, and optionally, the square error in // the dependent quantity, to the training sample to be used for the // parameterization. // The mean of the variables and quantity is calculated on the fly, // as outlined in TPrincipal::AddRow. // This sample should be representive of the problem at hand. // Please note, that if no error is given Poisson statistics is // assumed and the square error is set to the value of dependent // quantity. See also the // class description

if (!x) return; if (++fSampleSize == 1) { fMeanQuantity = D; fMaxQuantity = D; fMinQuantity = D; } else { fMeanQuantity *= 1 - 1./Double_t(fSampleSize); fMeanQuantity += D / Double_t(fSampleSize); fSumSqQuantity += D * D; if (D >= fMaxQuantity) fMaxQuantity = D; if (D <= fMinQuantity) fMinQuantity = D; } // If the vector isn't big enough to hold the new data, then // expand the vector by half it's size. Int_t size = fQuantity.GetNrows(); if (fSampleSize > size) { fQuantity.ResizeTo(size + size/2); fSqError.ResizeTo(size + size/2); } // Store the value fQuantity(fSampleSize-1) = D; fSqError(fSampleSize-1) = (E == 0 ? D : E); // Store data point in internal vector // If the vector isn't big enough to hold the new data, then // expand the vector by half it's size size = fVariables.GetNrows(); if (fSampleSize * fNVariables > size) fVariables.ResizeTo(size + size/2); // Increment the data point counter Int_t i,j; for (i = 0; i < fNVariables; i++) { if (fSampleSize == 1) { fMeanVariables(i) = x[i]; fMaxVariables(i) = x[i]; fMinVariables(i) = x[i]; } else { fMeanVariables(i) *= 1 - 1./Double_t(fSampleSize); fMeanVariables(i) += x[i] / Double_t(fSampleSize); // Update the maximum value for this component if (x[i] >= fMaxVariables(i)) fMaxVariables(i) = x[i]; // Update the minimum value for this component if (x[i] <= fMinVariables(i)) fMinVariables(i) = x[i]; } // Store the data. j = (fSampleSize-1) * fNVariables + i; fVariables(j) = x[i]; } } //____________________________________________________________________ void TMultiDimFit::AddTestRow(const Double_t *x, Double_t D, Double_t E) { // Add a row consisting of fNVariables independent variables, the // known, dependent quantity, and optionally, the square error in // the dependent quantity, to the test sample to be used for the // test of the parameterization. // This sample needn't be representive of the problem at hand. // Please note, that if no error is given Poisson statistics is // assumed and the square error is set to the value of dependent // quantity. See also the // class description

if (fTestSampleSize++ == 0) { fTestQuantity.ResizeTo(fNVariables); fTestSqError.ResizeTo(fNVariables); fTestVariables.ResizeTo(fNVariables * 100); } // If the vector isn't big enough to hold the new data, then // expand the vector by half it's size. Int_t size = fTestQuantity.GetNrows(); if (fTestSampleSize > size) { fTestQuantity.ResizeTo(size + size/2); fTestSqError.ResizeTo(size + size/2); } // Store the value fTestQuantity(fTestSampleSize-1) = D; fTestSqError(fTestSampleSize-1) = (E == 0 ? D : E); // Store data point in internal vector // If the vector isn't big enough to hold the new data, then // expand the vector by half it's size size = fTestVariables.GetNrows(); if (fTestSampleSize * fNVariables > size) fTestVariables.ResizeTo(size + size/2); // Increment the data point counter Int_t i,j; for (i = 0; i < fNVariables; i++) { j = fNVariables * (fTestSampleSize - 1) + i; fTestVariables(j) = x[i]; if (x[i] > fMaxVariables(i)) Warning("AddTestRow", "variable %d (row: %d) too large: %f > %f", i, fTestSampleSize, x[i], fMaxVariables(i)); if (x[i] < fMinVariables(i)) Warning("AddTestRow", "variable %d (row: %d) too small: %f < %f", i, fTestSampleSize, x[i], fMinVariables(i)); } } //____________________________________________________________________ void TMultiDimFit::Browse(TBrowser* b) { // Browse the TMultiDimFit object in the TBrowser. if (fHistograms) { TIter next(fHistograms); TH1* h = 0; while ((h = (TH1*)next())) b->Add(h,h->GetName()); } if (fVariables.IsValid()) b->Add(&fVariables, "Variables (Training)"); if (fQuantity.IsValid()) b->Add(&fQuantity, "Quantity (Training)"); if (fSqError.IsValid()) b->Add(&fSqError, "Error (Training)"); if (fMeanVariables.IsValid()) b->Add(&fMeanVariables, "Mean of Variables (Training)"); if (fMaxVariables.IsValid()) b->Add(&fMaxVariables, "Mean of Variables (Training)"); if (fMinVariables.IsValid()) b->Add(&fMinVariables, "Min of Variables (Training)"); if (fTestVariables.IsValid()) b->Add(&fTestVariables, "Variables (Test)"); if (fTestQuantity.IsValid()) b->Add(&fTestQuantity, "Quantity (Test)"); if (fTestSqError.IsValid()) b->Add(&fTestSqError, "Error (Test)"); if (fFunctions.IsValid()) b->Add(&fFunctions, "Functions"); if(fCoefficients.IsValid()) b->Add(&fCoefficients,"Coefficients"); if(fCoefficientsRMS.IsValid()) b->Add(&fCoefficientsRMS,"Coefficients Errors"); if (fOrthFunctions.IsValid()) b->Add(&fOrthFunctions, "Orthogonal Functions"); if (fOrthFunctionNorms.IsValid()) b->Add(&fOrthFunctionNorms, "Orthogonal Functions Norms"); if (fResiduals.IsValid()) b->Add(&fResiduals, "Residuals"); if(fOrthCoefficients.IsValid()) b->Add(&fOrthCoefficients,"Orthogonal Coefficients"); if (fOrthCurvatureMatrix.IsValid()) b->Add(&fOrthCurvatureMatrix,"Orthogonal curvature matrix"); if(fCorrelationMatrix.IsValid()) b->Add(&fCorrelationMatrix,"Correlation Matrix"); if (fFitter) b->Add(fFitter, fFitter->GetName()); } //____________________________________________________________________ void TMultiDimFit::Clear(Option_t *option) { // Clear internal structures and variables Int_t i, j, n = fNVariables, m = fMaxFunctions; // Training sample, dependent quantity fQuantity.Zero(); fSqError.Zero(); fMeanQuantity = 0; fMaxQuantity = 0; fMinQuantity = 0; fSumSqQuantity = 0; fSumSqAvgQuantity = 0; // Training sample, independent variables fVariables.Zero(); fNVariables = 0; fSampleSize = 0; fMeanVariables.Zero(); fMaxVariables.Zero(); fMinVariables.Zero(); // Test sample fTestQuantity.Zero(); fTestSqError.Zero(); fTestVariables.Zero(); fTestSampleSize = 0; // Functions fFunctions.Zero(); //for (i = 0; i < fMaxTerms; i++) fPowerIndex[i] = 0; //for (i = 0; i < fMaxTerms; i++) fFunctionCodes[i] = 0; fMaxFunctions = 0; fMaxStudy = 0; fOrthFunctions.Zero(); fOrthFunctionNorms.Zero(); // Control parameters fMinRelativeError = 0; fMinAngle = 0; fMaxAngle = 0; fMaxTerms = 0; // Powers for (i = 0; i < n; i++) { fMaxPowers[i] = 0; fMaxPowersFinal[i] = 0; for (j = 0; j < m; j++) fPowers[i * n + j] = 0; } fPowerLimit = 0; // Residuals fMaxResidual = 0; fMinResidual = 0; fMaxResidualRow = 0; fMinResidualRow = 0; fSumSqResidual = 0; // Fit fNCoefficients = 0; fOrthCoefficients = 0; fOrthCurvatureMatrix = 0; fRMS = 0; fCorrelationMatrix.Zero(); fError = 0; fTestError = 0; fPrecision = 0; fTestPrecision = 0; // Coefficients fCoefficients.Zero(); fCoefficientsRMS.Zero(); fResiduals.Zero(); fHistograms->Clear(option); // Options fPolyType = kMonomials; fShowCorrelation = kFALSE; fIsUserFunction = kFALSE; } //____________________________________________________________________ Double_t TMultiDimFit::Eval(const Double_t *x, const Double_t* coeff) { // Evaluate parameterization at point x. Optional argument coeff is // a vector of coefficients for the parameterisation, fNCoefficients // elements long. Double_t returnValue = fMeanQuantity; Double_t term = 0; Int_t i, j; for (i = 0; i < fNCoefficients; i++) { // Evaluate the ith term in the expansion term = (coeff ? coeff[i] : fCoefficients(i)); for (j = 0; j < fNVariables; j++) { // Evaluate the factor (polynomial) in the j-th variable. Int_t p = fPowers[fPowerIndex[i] * fNVariables + j]; Double_t X = 1 + 2. / (fMaxVariables(j) - fMinVariables(j)) * (x[j] - fMaxVariables(j)); term *= EvalFactor(p,X); } // Add this term to the final result returnValue += term; } return returnValue; } //____________________________________________________________________ Double_t TMultiDimFit::EvalControl(const Int_t *iv) { // PRIVATE METHOD: // Calculate the control parameter from the passed powers Double_t s = 0; Double_t epsilon = 1e-6; // a small number for (Int_t i = 0; i < fNVariables; i++) { if (fMaxPowers[i] != 1) s += (epsilon + iv[i] - 1) / (epsilon + fMaxPowers[i] - 1); } return s; } //____________________________________________________________________ Double_t TMultiDimFit::EvalFactor(Int_t p, Double_t x) { // PRIVATE METHOD: // Evaluate function with power p at variable value x Int_t i = 0; Double_t p1 = 1; Double_t p2 = 0; Double_t p3 = 0; Double_t r = 0; switch(p) { case 1: r = 1; break; case 2: r = x; break; default: p2 = x; for (i = 3; i <= p; i++) { p3 = p2 * x; if (fPolyType == kLegendre) p3 = ((2 * i - 3) * p2 * x - (i - 2) * p1) / (i - 1); else if (fPolyType == kChebyshev) p3 = 2 * x * p2 - p1; p1 = p2; p2 = p3; } r = p3; } return r; } //____________________________________________________________________ void TMultiDimFit::FindParameterization(Option_t *) { // Find the parameterization // // Options: // None so far // // For detailed description of what this entails, please refer to the // class description

MakeNormalized(); MakeCandidates(); MakeParameterization(); MakeCoefficients(); MakeCoefficientErrors(); MakeCorrelation(); } //____________________________________________________________________ void TMultiDimFit::Fit(Option_t *option) { // Try to fit the found parameterisation to the test sample. // // Options // M use Minuit to improve coefficients // // Also, refer to // class description

Int_t i, j; Double_t* x = new Double_t[fNVariables]; Double_t sumSqD = 0; Double_t sumD = 0; Double_t sumSqR = 0; Double_t sumR = 0; // Calculate the residuals over the test sample for (i = 0; i < fTestSampleSize; i++) { for (j = 0; j < fNVariables; j++) x[j] = fTestVariables(i * fNVariables + j); Double_t res = fTestQuantity(i) - Eval(x); sumD += fTestQuantity(i); sumSqD += fTestQuantity(i) * fTestQuantity(i); sumR += res; sumSqR += res * res; if (TESTBIT(fHistogramMask,HIST_RTEST)) ((TH1F*)fHistograms->FindObject("res_test"))->Fill(res); } Double_t dAvg = sumSqD - (sumD * sumD) / fTestSampleSize; Double_t rAvg = sumSqR - (sumR * sumR) / fTestSampleSize; fTestCorrelationCoeff = (dAvg - rAvg) / dAvg; fTestError = sumSqR; fTestPrecision = sumSqR / sumSqD; TString opt(option); opt.ToLower(); if (!opt.Contains("m")) MakeChi2(); if (fNCoefficients * 50 > fTestSampleSize) Warning("Fit", "test sample is very small"); if (!opt.Contains("m")) return; fFitter = TVirtualFitter::Fitter(0,fNCoefficients); fFitter->SetFCN(mdfHelper); const Int_t maxArgs = 16; Int_t args = 1; Double_t* arglist = new Double_t[maxArgs]; arglist[0] = -1; fFitter->ExecuteCommand("SET PRINT",arglist,args); for (i = 0; i < fNCoefficients; i++) { Double_t startVal = fCoefficients(i); Double_t startErr = fCoefficientsRMS(i); fFitter->SetParameter(i, Form("coeff%02d",i), startVal, startErr, 0, 0); } // arglist[0] = 0; args = 1; // fFitter->ExecuteCommand("SET PRINT",arglist,args); fFitter->ExecuteCommand("MIGRAD",arglist,args); for (i = 0; i < fNCoefficients; i++) { Double_t val = 0, err = 0, low = 0, high = 0; fFitter->GetParameter(i, Form("coeff%02d",i), val, err, low, high); fCoefficients(i) = val; fCoefficientsRMS(i) = err; } } //____________________________________________________________________ void TMultiDimFit::MakeCandidates() { // PRIVATE METHOD: // Create list of candidate functions for the parameterisation. See // also // class description

Int_t i = 0; Int_t j = 0; Int_t k = 0; // The temporary array to store the powers in. We don't need to // initialize this array however. Int_t *powers = new Int_t[fNVariables * fMaxFunctions]; // store of `control variables' Double_t* control = new Double_t[fMaxFunctions]; // We've better initialize the variables Int_t *iv = new Int_t[fNVariables]; for (i = 0; i < fNVariables; i++) iv[i] = 1; if (!fIsUserFunction) { // Number of funcs selected Int_t numberFunctions = 0; // Absolute max number of functions Int_t maxNumberFunctions = 1; for (i = 0; i < fNVariables; i++) maxNumberFunctions *= fMaxPowers[i]; while (kTRUE) { // Get the control value for this function Double_t s = EvalControl(iv); if (s <= fPowerLimit) { // Call over-loadable method Select, as to allow the user to // interfere with the selection of functions. if (Select(iv)) { numberFunctions++; // If we've reached the user defined limit of how many // functions we can consider, break out of the loop if (numberFunctions > fMaxFunctions) break; // Store the control value, so we can sort array of powers // later on control[numberFunctions-1] = s; // Store the powers in powers array. for (i = 0; i < fNVariables; i++) { j = (numberFunctions - 1) * fNVariables + i; powers[j] = iv[i]; } } // if (Select()) } // if (s <= fPowerLimit) for (i = 0; i < fNVariables; i++) if (iv[i] < fMaxPowers[i]) break; // If all variables have reached their maximum power, then we // break out of the loop if (i == fNVariables) { fMaxFunctions = numberFunctions; break; } // Next power in variable i iv[i]++; for (j = 0; j < i; j++) iv[j] = 1; } // while (kTRUE) } else { // In case the user gave an explicit function for (i = 0; i < fMaxFunctions; i++) { // Copy the powers to working arrays for (j = 0; j < fNVariables; j++) { powers[i * fNVariables + j] = fPowers[i * fNVariables + j]; iv[j] = fPowers[i * fNVariables + j]; } control[i] = EvalControl(iv); } } // Now we need to sort the powers according to least `control // variable' Int_t *order = new Int_t[fMaxFunctions]; for (i = 0; i < fMaxFunctions; i++) order[i] = i; fPowers = new Int_t[fMaxFunctions * fNVariables]; for (i = 0; i < fMaxFunctions; i++) { Double_t x = control[i]; Int_t l = order[i]; k = i; for (j = i; j < fMaxFunctions; j++) { if (control[j] <= x) { x = control[j]; l = order[j]; k = j; } } if (k != i) { control[k] = control[i]; control[i] = x; order[k] = order[i]; order[i] = l; } } for (i = 0; i < fMaxFunctions; i++) for (j = 0; j < fNVariables; j++) fPowers[i * fNVariables + j] = powers[order[i] * fNVariables + j]; delete [] control; delete [] powers; delete [] order; delete [] iv; } //____________________________________________________________________ Double_t TMultiDimFit::MakeChi2(const Double_t* coeff) { // Calculate Chi square over either the test sample. The optional // argument coeff is a vector of coefficients to use in the // evaluation of the parameterisation. If coeff == 0, then the found // coefficients is used. // Used my MINUIT for fit (see TMultDimFit::Fit) fChi2 = 0; Int_t i, j; Double_t* x = new Double_t[fNVariables]; for (i = 0; i < fTestSampleSize; i++) { // Get the stored point for (j = 0; j < fNVariables; j++) x[j] = fTestVariables(i * fNVariables + j); // Evaluate function. Scale to shifted values Double_t f = Eval(x,coeff); // Calculate contribution to Chic square fChi2 += 1. / TMath::Max(fTestSqError(i),1e-20) * (fTestQuantity(i) - f) * (fTestQuantity(i) - f); } // Clean up delete [] x; return fChi2; } //____________________________________________________________________ void TMultiDimFit::MakeCode(const char* filename, Option_t *option) { // Generate the file <filename> with .C appended if argument doesn't // end in .cxx or .C. The contains the implementation of the // function: // // Double_t <funcname>(Double_t *x) // // which does the same as TMultiDimFit::Eval. Please refer to this // method. // // Further, the static variables: // // Int_t gNVariables // Int_t gNCoefficients // Double_t gDMean // Double_t gXMean[] // Double_t gXMin[] // Double_t gXMax[] // Double_t gCoefficient[] // Int_t gPower[] // // are initialized. The only ROOT header file needed is Rtypes.h // // See TMultiDimFit::MakeRealCode for a list of options TString outName(filename); if (!outName.EndsWith(".C") && !outName.EndsWith(".cxx")) outName += ".C"; MakeRealCode(outName.Data(),"",option); } //____________________________________________________________________ void TMultiDimFit::MakeCoefficientErrors() { // PRIVATE METHOD: // Compute the errors on the coefficients. For this to be done, the // curvature matrix of the non-orthogonal functions, is computed. Int_t i = 0; Int_t j = 0; Int_t k = 0; TVectorD iF(fSampleSize); TVectorD jF(fSampleSize); fCoefficientsRMS.ResizeTo(fNCoefficients); TMatrixD curvatureMatrix(fNCoefficients,fNCoefficients); // Build the curvature matrix for (i = 0; i < fNCoefficients; i++) { iF = TMatrixDRow(fFunctions,i); for (j = 0; j <= i; j++) { jF = TMatrixDRow(fFunctions,j); for (k = 0; k < fSampleSize; k++) curvatureMatrix(i,j) += 1 / TMath::Max(fSqError(k), 1e-20) * iF(k) * jF(k); curvatureMatrix(j,i) = curvatureMatrix(i,j); } } // Calculate Chi Square fChi2 = 0; for (i = 0; i < fSampleSize; i++) { Double_t f = 0; for (j = 0; j < fNCoefficients; j++) f += fCoefficients(j) * fFunctions(j,i); fChi2 += 1. / TMath::Max(fSqError(i),1e-20) * (fQuantity(i) - f) * (fQuantity(i) - f); } // Invert the curvature matrix Double_t det = 1; curvatureMatrix.Invert(&det); if (det == 0) Warning("MakeCoefficientErrors", "curvature matrix is singular"); if (fIsVerbose) cout << "Determinant of curvature matrix: " << det << endl; for (i = 0; i < fNCoefficients; i++) fCoefficientsRMS(i) = TMath::Sqrt(curvatureMatrix(i,i)); } //____________________________________________________________________ void TMultiDimFit::MakeCoefficients() { // PRIVATE METHOD: // Invert the model matrix B, and compute final coefficients. For a // more thorough discussion of what this means, please refer to the // class description

// // First we invert the lower triangle matrix fOrthCurvatureMatrix // and store the inverted matrix in the upper triangle. Int_t i = 0, j = 0; Int_t col = 0, row = 0; // Invert the B matrix for (col = 1; col < fNCoefficients; col++) { for (row = col - 1; row > -1; row--) { fOrthCurvatureMatrix(row,col) = 0; for (i = row; i <= col ; i++) fOrthCurvatureMatrix(row,col) -= fOrthCurvatureMatrix(i,row) * fOrthCurvatureMatrix(i,col); } } // Compute the final coefficients fCoefficients.ResizeTo(fNCoefficients); for (i = 0; i < fNCoefficients; i++) { Double_t sum = 0; for (j = i; j < fNCoefficients; j++) sum += fOrthCurvatureMatrix(i,j) * fOrthCoefficients(j); fCoefficients(i) = sum; } // Compute the final residuals fResiduals.ResizeTo(fSampleSize); for (i = 0; i < fSampleSize; i++) fResiduals(i) = fQuantity(i); for (i = 0; i < fNCoefficients; i++) for (j = 0; j < fSampleSize; j++) fResiduals(j) -= fCoefficients(i) * fFunctions(i,j); // Compute the max and minimum, and squared sum of the evaluated // residuals fMinResidual = 10e10; fMaxResidual = -10e10; Double_t sqRes = 0; for (i = 0; i < fSampleSize; i++){ sqRes += fResiduals(i) * fResiduals(i); if (fResiduals(i) <= fMinResidual) { fMinResidual = fResiduals(i); fMinResidualRow = i; } if (fResiduals(i) >= fMaxResidual) { fMaxResidual = fResiduals(i); fMaxResidualRow = i; } } fCorrelationCoeff = fSumSqResidual / fSumSqAvgQuantity; fPrecision = TMath::Sqrt(sqRes / fSumSqQuantity); // If we use histograms, fill some more if (TESTBIT(fHistogramMask,HIST_RD) || TESTBIT(fHistogramMask,HIST_RTRAI) || TESTBIT(fHistogramMask,HIST_RX)) { for (i = 0; i < fSampleSize; i++) { if (TESTBIT(fHistogramMask,HIST_RD)) ((TH2D*)fHistograms->FindObject("res_d"))->Fill(fQuantity(i), fResiduals(i)); if (TESTBIT(fHistogramMask,HIST_RTRAI)) ((TH1D*)fHistograms->FindObject("res_train"))->Fill(fResiduals(i)); if (TESTBIT(fHistogramMask,HIST_RX)) for (j = 0; j < fNVariables; j++) ((TH2D*)fHistograms->FindObject(Form("res_x_%d",j))) ->Fill(fVariables(i * fNVariables + j),fResiduals(i)); } } // If histograms } //____________________________________________________________________ void TMultiDimFit::MakeCorrelation() { // PRIVATE METHOD: // Compute the correlation matrix if (!fShowCorrelation) return; fCorrelationMatrix.ResizeTo(fNVariables,fNVariables+1); Double_t d2 = 0; Double_t DdotXi = 0; Double_t XiNorm = 0; Double_t XidotXj = 0; Double_t XjNorm = 0; Int_t i, j, k, l; for (i = 0; i < fSampleSize; i++) d2 += fQuantity(i) * fQuantity(i); for (i = 0; i < fNVariables; i++) { for (j = 0; j< fSampleSize; j++) { // Index of sample j of variable i k = j * fNVariables + i; DdotXi += fQuantity(j) * (fVariables(k) - fMeanVariables(i)); XiNorm += (fVariables(k) - fMeanVariables(i)) * (fVariables(k) - fMeanVariables(i)); } fCorrelationMatrix(i,0) = DdotXi / TMath::Sqrt(d2 * XiNorm); for (j = 0; j < i; j++) { for (k = 0; k < fSampleSize; k++) { // Index of sample j of variable i l = j * fNVariables + k; XidotXj += (fVariables(i) - fMeanVariables(i)) * (fVariables(l) - fMeanVariables(j)); XjNorm += (fVariables(l) - fMeanVariables(j)) * (fVariables(l) - fMeanVariables(j)); } fCorrelationMatrix(i+1,j) = XidotXj / TMath::Sqrt(XiNorm * XjNorm); } } } //____________________________________________________________________ Double_t TMultiDimFit::MakeGramSchmidt(Int_t function) { // PRIVATE METHOD: // Make Gram-Schmidt orthogonalisation. The class description gives // a thorough account of this algorithm, as well as // references. Please refer to the // class description

// calculate w_i, that is, evaluate the current function at data // point i Double_t f2 = 0; fOrthCoefficients(fNCoefficients) = 0; fOrthFunctionNorms(fNCoefficients) = 0; Int_t j = 0; Int_t k = 0; for (j = 0; j < fSampleSize; j++) { fFunctions(fNCoefficients, j) = 1; fOrthFunctions(fNCoefficients, j) = 0; // First, however, we need to calculate f_fNCoefficients for (k = 0; k < fNVariables; k++) { Int_t p = fPowers[function * fNVariables + k]; Double_t x = fVariables(j * fNVariables + k); fFunctions(fNCoefficients, j) *= EvalFactor(p,x); } // Calculate f dot f in f2 f2 += fFunctions(fNCoefficients,j) * fFunctions(fNCoefficients,j); // Assign to w_fNCoefficients f_fNCoefficients fOrthFunctions(fNCoefficients, j) = fFunctions(fNCoefficients, j); } // the first column of w is equal to f for (j = 0; j < fNCoefficients; j++) { Double_t fdw = 0; // Calculate (f_fNCoefficients dot w_j) / w_j^2 for (k = 0; k < fSampleSize; k++) { fdw += fFunctions(fNCoefficients, k) * fOrthFunctions(j,k) / fOrthFunctionNorms(j); } fOrthCurvatureMatrix(fNCoefficients,j) = fdw; // and subtract it from the current value of w_ij for (k = 0; k < fSampleSize; k++) fOrthFunctions(fNCoefficients,k) -= fdw * fOrthFunctions(j,k); } for (j = 0; j < fSampleSize; j++) { // calculate squared length of w_fNCoefficients fOrthFunctionNorms(fNCoefficients) += fOrthFunctions(fNCoefficients,j) * fOrthFunctions(fNCoefficients,j); // calculate D dot w_fNCoefficients in A fOrthCoefficients(fNCoefficients) += fQuantity(j) * fOrthFunctions(fNCoefficients, j); } // First test, but only if didn't user specify if (!fIsUserFunction) if (TMath::Sqrt(fOrthFunctionNorms(fNCoefficients) / (f2 + 1e-10)) < TMath::Sin(fMinAngle*DEGRAD)) return 0; // The result found by this code for the first residual is always // much less then the one found be MUDIFI. That's because it's // supposed to be. The cause is the improved precision of Double_t // over DOUBLE PRECISION! fOrthCurvatureMatrix(fNCoefficients,fNCoefficients) = 1; Double_t b = fOrthCoefficients(fNCoefficients); fOrthCoefficients(fNCoefficients) /= fOrthFunctionNorms(fNCoefficients); // Calculate the residual from including this fNCoefficients. Double_t dResidur = fOrthCoefficients(fNCoefficients) * b; return dResidur; } //____________________________________________________________________ void TMultiDimFit::MakeHistograms(Option_t *option) { // Make histograms of the result of the analysis. This message // should be sent after having read all data points, but before // finding the parameterization // // Options: // A All the below // X Original independent variables // D Original dependent variables // N Normalised independent variables // S Shifted dependent variables // R1 Residuals versus normalised independent variables // R2 Residuals versus dependent variable // R3 Residuals computed on training sample // R4 Residuals computed on test sample // // For a description of these quantities, refer to // class description

TString opt(option); opt.ToLower(); if (opt.Length() < 1) return; if (!fHistograms) fHistograms = new TList; // Counter variable Int_t i = 0; // Histogram of original variables if (opt.Contains("x") || opt.Contains("a")) { SETBIT(fHistogramMask,HIST_XORIG); for (i = 0; i < fNVariables; i++) if (!fHistograms->FindObject(Form("x_%d_orig",i))) fHistograms->Add(new TH1D(Form("x_%d_orig",i), Form("Original variable # %d",i), 100, fMinVariables(i), fMaxVariables(i))); } // Histogram of original dependent variable if (opt.Contains("d") || opt.Contains("a")) { SETBIT(fHistogramMask,HIST_DORIG); if (!fHistograms->FindObject("d_orig")) fHistograms->Add(new TH1D("d_orig", "Original Quantity", 100, fMinQuantity, fMaxQuantity)); } // Histograms of normalized variables if (opt.Contains("n") || opt.Contains("a")) { SETBIT(fHistogramMask,HIST_XNORM); for (i = 0; i < fNVariables; i++) if (!fHistograms->FindObject(Form("x_%d_norm",i))) fHistograms->Add(new TH1D(Form("x_%d_norm",i), Form("Normalized variable # %d",i), 100, -1,1)); } // Histogram of shifted dependent variable if (opt.Contains("s") || opt.Contains("a")) { SETBIT(fHistogramMask,HIST_DSHIF); if (!fHistograms->FindObject("d_shifted")) fHistograms->Add(new TH1D("d_shifted", "Shifted Quantity", 100, fMinQuantity - fMeanQuantity, fMaxQuantity - fMeanQuantity)); } // Residual from training sample versus independent variables if (opt.Contains("r1") || opt.Contains("a")) { SETBIT(fHistogramMask,HIST_RX); for (i = 0; i < fNVariables; i++) if (!fHistograms->FindObject(Form("res_x_%d",i))) fHistograms->Add(new TH2D(Form("res_x_%d",i), Form("Computed residual versus x_%d", i), 100, -1, 1, 35, fMinQuantity - fMeanQuantity, fMaxQuantity - fMeanQuantity)); } // Residual from training sample versus. dependent variable if (opt.Contains("r2") || opt.Contains("a")) { SETBIT(fHistogramMask,HIST_RD); if (!fHistograms->FindObject("res_d")) fHistograms->Add(new TH2D("res_d", "Computed residuals vs Quantity", 100, fMinQuantity - fMeanQuantity, fMaxQuantity - fMeanQuantity, 35, fMinQuantity - fMeanQuantity, fMaxQuantity - fMeanQuantity)); } // Residual from training sample if (opt.Contains("r3") || opt.Contains("a")) { SETBIT(fHistogramMask,HIST_RTRAI); if (!fHistograms->FindObject("res_train")) fHistograms->Add(new TH1D("res_train", "Computed residuals over training sample", 100, fMinQuantity - fMeanQuantity, fMaxQuantity - fMeanQuantity)); } if (opt.Contains("r4") || opt.Contains("a")) { SETBIT(fHistogramMask,HIST_RTEST); if (!fHistograms->FindObject("res_test")) fHistograms->Add(new TH1F("res_test", "Distribution of residuals from test", 100,fMinQuantity - fMeanQuantity, fMaxQuantity - fMeanQuantity)); } } //____________________________________________________________________ void TMultiDimFit::MakeMethod(const Char_t* classname, Option_t* option) { // Generate the file <classname>MDF.cxx which contains the // implementation of the method: // // Double_t <classname>::MDF(Double_t *x) // // which does the same as TMultiDimFit::Eval. Please refer to this // method. // // Further, the public static members: // // Int_t <classname>::fgNVariables // Int_t <classname>::fgNCoefficients // Double_t <classname>::fgDMean // Double_t <classname>::fgXMean[] //[fgNVariables] // Double_t <classname>::fgXMin[] //[fgNVariables] // Double_t <classname>::fgXMax[] //[fgNVariables] // Double_t <classname>::fgCoefficient[] //[fgNCoeffficents] // Int_t <classname>::fgPower[] //[fgNCoeffficents*fgNVariables] // // are initialized, and assumed to exist. The class declaration is // assumed to be in <classname>.h and assumed to be provided by the // user. // // See TMultiDimFit::MakeRealCode for a list of options // // The minimal class definition is: // // class <classname> { // public: // Int_t <classname>::fgNVariables; // Number of variables // Int_t <classname>::fgNCoefficients; // Number of terms // Double_t <classname>::fgDMean; // Mean from training sample // Double_t <classname>::fgXMean[]; // Mean from training sample // Double_t <classname>::fgXMin[]; // Min from training sample // Double_t <classname>::fgXMax[]; // Max from training sample // Double_t <classname>::fgCoefficient[]; // Coefficients // Int_t <classname>::fgPower[]; // Function powers // // Double_t Eval(Double_t *x); // }; // // Whether the method <classname>::Eval should be static or not, is // up to the user. MakeRealCode(Form("%sMDF.cxx", classname), classname, option); } //____________________________________________________________________ void TMultiDimFit::MakeNormalized() { // PRIVATE METHOD: // Normalize data to the interval [-1;1]. This is needed for the // classes method to work. Int_t i = 0; Int_t j = 0; Int_t k = 0; for (i = 0; i < fSampleSize; i++) { if (TESTBIT(fHistogramMask,HIST_DORIG)) ((TH1F*)fHistograms->FindObject("d_orig"))->Fill(fQuantity(i)); fQuantity(i) -= fMeanQuantity; fSumSqAvgQuantity += fQuantity(i) * fQuantity(i); if (TESTBIT(fHistogramMask,HIST_DSHIF)) ((TH1F*)fHistograms->FindObject("d_shifted"))->Fill(fQuantity(i)); for (j = 0; j < fNVariables; j++) { Double_t range = 1. / (fMaxVariables(j) - fMinVariables(j)); k = i * fNVariables + j; // Fill histograms of original independent variables if (TESTBIT(fHistogramMask,HIST_XORIG)) ((TH1F*)fHistograms->FindObject(Form("x_%d_orig",j))) ->Fill(fVariables(k)); // Normalise independent variables fVariables(k) = 1 + 2 * range * (fVariables(k) - fMaxVariables(j)); // Fill histograms of normalised independent variables if (TESTBIT(fHistogramMask,HIST_XNORM)) ((TH1F*)fHistograms->FindObject(Form("x_%d_norm",j))) ->Fill(fVariables(k)); } } // Shift min and max of dependent variable fMaxQuantity -= fMeanQuantity; fMinQuantity -= fMeanQuantity; // Shift mean of independent variables for (i = 0; i < fNVariables; i++) { Double_t range = 1. / (fMaxVariables(i) - fMinVariables(i)); fMeanVariables(i) = 1 + 2 * range * (fMeanVariables(i) - fMaxVariables(i)); } } //____________________________________________________________________ void TMultiDimFit::MakeParameterization() { // PRIVATE METHOD: // Find the parameterization over the training sample. A full account // of the algorithm is given in the // class description

Int_t i = -1; Int_t j = 0; Int_t k = 0; Int_t maxPass = 3; Int_t studied = 0; Double_t squareResidual = fSumSqAvgQuantity; fNCoefficients = 0; fSumSqResidual = fSumSqAvgQuantity; fFunctions.ResizeTo(fMaxTerms,fSampleSize); fOrthFunctions.ResizeTo(fMaxTerms,fSampleSize); fOrthFunctionNorms.ResizeTo(fMaxTerms); fOrthCoefficients.ResizeTo(fMaxTerms); fOrthCurvatureMatrix.ResizeTo(fMaxTerms,fMaxTerms); fFunctions = 1; fFunctionCodes = new Int_t[fMaxFunctions]; fPowerIndex = new Int_t[fMaxTerms]; Int_t l; for (l=0;l<fMaxFunctions;l++) fFunctionCodes[l] = 0; for (l=0;l<fMaxTerms;l++) fPowerIndex[l] = 0; if (fMaxAngle != 0) maxPass = 100; if (fIsUserFunction) maxPass = 1; // Loop over the number of functions we want to study. // increment inspection counter while(kTRUE) { // Reach user defined limit of studies if (studied++ >= fMaxStudy) { fParameterisationCode = PARAM_MAXSTUDY; break; } // Considered all functions several times if (k >= maxPass) { fParameterisationCode = PARAM_SEVERAL; break; } // increment function counter i++; // If we've reached the end of the functions, restart pass if (i == fMaxFunctions) { if (fMaxAngle != 0) fMaxAngle += (90 - fMaxAngle) / 2; i = 0; studied--; k++; continue; } if (studied == 1) fFunctionCodes[i] = 0; else if (fFunctionCodes[i] >= 2) continue; // Print a happy message if (fIsVerbose && studied == 1) cout << "Coeff SumSqRes Contrib Angle QM Func" << " Value W^2 Powers" << endl; // Make the Gram-Schmidt Double_t dResidur = MakeGramSchmidt(i); if (dResidur == 0) { // This function is no good! // First test is in MakeGramSchmidt fFunctionCodes[i] = 1; continue; } // If user specified function, assume she/he knows what he's doing if (!fIsUserFunction) { // Flag this function as considered fFunctionCodes[i] = 2; // Test if this function contributes to the fit if (!TestFunction(squareResidual, dResidur)) { fFunctionCodes[i] = 1; continue; } } // If we get to here, the function currently considered is // fNCoefficients, so we increment the counter // Flag this function as OK, and store and the number in the // index. fFunctionCodes[i] = 3; fPowerIndex[fNCoefficients] = i; fNCoefficients++; // We add the current contribution to the sum of square of // residuals; squareResidual -= dResidur; // Calculate control parameter from this function for (j = 0; j < fNVariables; j++) { if (fNCoefficients == 1 || fMaxPowersFinal[j] <= fPowers[i * fNVariables + j] - 1) fMaxPowersFinal[j] = fPowers[i * fNVariables + j] - 1; } Double_t s = EvalControl(&fPowers[i * fNVariables]); // Print the statistics about this function if (fIsVerbose) { #ifndef R__MACOSX cout << setw(5) << fNCoefficients << " " << setw(10) << setprecision(4) << squareResidual << " " << setw(10) << setprecision(4) << dResidur << " " << setw(7) << setprecision(3) << fMaxAngle << " " << setw(7) << setprecision(3) << s << " " << setw(5) << i << " " << setw(10) << setprecision(4) << fOrthCoefficients(fNCoefficients-1) << " " << setw(10) << setprecision(4) << fOrthFunctionNorms(fNCoefficients-1) << " " << flush; #else printf("%d5 %g10.4 %g10.4 %g7.3 %g7.3 %d5 %g10.4 %g10.4 ", fNCoefficients, squareResidual, dResidur, fMaxAngle, s, i, fOrthCoefficients(fNCoefficients-1), fOrthFunctionNorms(fNCoefficients-1)); fflush(stdout); #endif for (j = 0; j < fNVariables; j++) cout << " " << fPowers[i * fNVariables + j] - 1 << flush; cout << endl; } if (fNCoefficients >= fMaxTerms /* && fIsVerbose */) { fParameterisationCode = PARAM_MAXTERMS; break; } Double_t err = TMath::Sqrt(TMath::Max(1e-20,squareResidual) / fSumSqAvgQuantity); if (err < fMinRelativeError) { fParameterisationCode = PARAM_RELERR; break; } } fError = TMath::Max(1e-20,squareResidual); fSumSqResidual -= fError; fRMS = TMath::Sqrt(fError / fSampleSize); } //____________________________________________________________________ void TMultiDimFit::MakeRealCode(const char *filename, const char *classname, Option_t *) { // PRIVATE METHOD: // This is the method that actually generates the code for the // evaluation the parameterization on some point. // It's called by TMultiDimFit::MakeCode and TMultiDimFit::MakeMethod. // // The options are: NONE so far Int_t i, j; Bool_t isMethod = (classname[0] == '0' ? kFALSE : kTRUE); const char *prefix = (isMethod ? Form("%s::", classname) : ""); const char *cv_qual = (isMethod ? "" : "static "); ofstream outFile(filename,ios::out if (!outFile) { Error("MakeRealCode","couldn't open output file '%s'",filename); return; } if (fIsVerbose) cout << "Writing on file "" << filename << "" ... " << flush; // // Write header of file // // Emacs mode line ;-) outFile << "// -*- mode: c++ -*-" << endl; // Info about creator outFile << "// " << endl << "// File " << filename << " generated by TMultiDimFit::MakeRealCode" << endl; // Time stamp TDatime date; outFile << "// on " << date.AsString() << endl; // ROOT version info outFile << "// ROOT version " << gROOT->GetVersion() << endl << "//" << endl; // General information on the code outFile << "// This file contains the function " << endl << "//" << endl << "// double " << prefix << "MDF(double *x); " << endl << "//" << endl << "// For evaluating the parameterization obtained" << endl << "// from TMultiDimFit and the point x" << endl << "// " << endl << "// See TMultiDimFit class documentation for more " << "information " << endl << "// " << endl; // Header files if (isMethod) // If these are methods, we need the class header outFile << "#include "" << classname << ".h"" << endl; // // Now for the data // outFile << "//" << endl << "// Static data variables" << endl << "//" << endl; outFile << cv_qual << "int " << prefix << "gNVariables = " << fNVariables << ";" << endl; outFile << cv_qual << "int " << prefix << "gNCoefficients = " << fNCoefficients << ";" << endl; outFile << cv_qual << "double " << prefix << "gDMean = " << fMeanQuantity << ";" << endl; // Assignment to mean vector. outFile << "// Assignment to mean vector." << endl; outFile << cv_qual << "double " << prefix << "gXMean[] = {" << endl; for (i = 0; i < fNVariables; i++) outFile << (i != 0 ? ", " : " ") << fMeanVariables(i) << flush; outFile << " };" << endl << endl; // Assignment to minimum vector. outFile << "// Assignment to minimum vector." << endl; outFile << cv_qual << "double " << prefix << "gXMin[] = {" << endl; for (i = 0; i < fNVariables; i++) outFile << (i != 0 ? ", " : " ") << fMinVariables(i) << flush; outFile << " };" << endl << endl; // Assignment to maximum vector. outFile << "// Assignment to maximum vector." << endl; outFile << cv_qual << "double " << prefix << "gXMax[] = {" << endl; for (i = 0; i < fNVariables; i++) outFile << (i != 0 ? ", " : " ") << fMaxVariables(i) << flush; outFile << " };" << endl << endl; // Assignment to coefficients vector. outFile << "// Assignment to coefficients vector." << endl; outFile << cv_qual << "double " << prefix << "gCoefficient[] = {" << flush; for (i = 0; i < fNCoefficients; i++) outFile << (i != 0 ? "," : "") << endl << " " << fCoefficients(i) << flush; outFile << endl << " };" << endl << endl; // Assignment to powers vector. outFile << "// Assignment to powers vector." << endl << "// The powers are stored row-wise, that is" << endl << "// p_ij = " << prefix << "gPower[i * NVariables + j];" << endl; outFile << cv_qual << "int " << prefix << "gPower[] = {" << flush; for (i = 0; i < fNCoefficients; i++) { for (j = 0; j < fNVariables; j++) { if (j != 0) outFile << flush << " "; else outFile << endl << " "; outFile << fPowers[fPowerIndex[i] * fNVariables + j] << (i == fNCoefficients - 1 && j == fNVariables - 1 ? "" : ",") << flush; } } outFile << endl << "};" << endl << endl; // // Finally we reach the function itself // outFile << "// " << endl << "// The " << (isMethod ? "method " : "function ") << " double " << prefix << "MDF(double *x)" << endl << "// " << endl; outFile << "double " << prefix << "MDF(double *x) {" << endl << " double returnValue = " << prefix << "gDMean;" << endl << " int i = 0, j = 0;" << endl << " for (i = 0; i < " << prefix << "gNCoefficients ; i++) {" << endl << " // Evaluate the ith term in the expansion" << endl << " double term = " << prefix << "gCoefficient[i]; <<" << endl << " for (j = 0; j < " << prefix << "gNVariables; j++) {" << endl << " // Evaluate the polynomial in the jth variable." << endl << " int power = "<< prefix << "gPower[" << prefix << "gNVariables * i + j]; " << endl << " double p1 = 1, p2 = 0, p3 = 0, r = 0;" << endl << " double v = 1 + 2. / (" << prefix << "gXMax[j] - " << prefix << "gXMin[j]) * (x[j] - " << prefix << "gXMax[j]);" << endl << " // what is the power to use!" << endl << " switch(power) {" << endl << " case 1: r = 1; break; " << endl << " case 2: r = v; break; " << endl << " default: " << endl << " p2 = v; " << endl << " for (i = 3; i <= power; i++) { " << endl << " p3 = p2 * v;" << endl; if (fPolyType == kLegendre) outFile << " p3 = ((2 * i - 3) * p2 * v - (i - 2) * p1)" << " / (i - 1);" << endl; if (fPolyType == kChebyshev) outFile << " p3 = 2 * v * p2 - p1; " << endl; outFile << " p1 = p2; p2 = p3; " << endl << " }" << endl << " r = p3;" << endl << " }" << endl << " // multiply this term by the poly in the jth var" << endl << " term *= r; " << endl << " }" << endl << " // Add this term to the final result" << endl << " returnValue += term;" << endl << " }" << endl << " return returnValue;" << endl << "}" << endl << endl; // EOF outFile << "// EOF for " << filename << endl; // Close the file outFile.close(); if (fIsVerbose) cout << "done" << endl; } //____________________________________________________________________ void TMultiDimFit::Print(Option_t *option) const { // Print statistics etc. // Options are // P Parameters // S Statistics // C Coefficients // R Result of parameterisation // F Result of fit // K Correlation Matrix // M Pretty print formula // Int_t i = 0; Int_t j = 0; TString opt(option); opt.ToLower(); if (opt.Contains("p")) { // Print basic parameters for this object cout << "User parameters:" << endl << "----------------" << endl << " Variables: " << fNVariables << endl << " Data points: " << fSampleSize << endl << " Max Terms: " << fMaxTerms << endl << " Power Limit Parameter: " << fPowerLimit << endl << " Max functions: " << fMaxFunctions << endl << " Max functions to study: " << fMaxStudy << endl << " Max angle (optional): " << fMaxAngle << endl << " Min angle: " << fMinAngle << endl << " Relative Error accepted: " << fMinRelativeError << endl << " Maximum Powers: " << flush; for (i = 0; i < fNVariables; i++) cout << " " << fMaxPowers[i] - 1 << flush; cout << endl << endl << " Parameterisation will be done using " << flush; if (fPolyType == kChebyshev) cout << "Chebyshev polynomials" << endl; else if (fPolyType == kLegendre) cout << "Legendre polynomials" << endl; else cout << "Monomials" << endl; cout << endl; } if (opt.Contains("s")) { // Print statistics for read data cout << "Sample statistics:" << endl << "------------------" << endl << " D" << flush; #ifndef R__MACOSX for (i = 0; i < fNVariables; i++) cout << " " << setw(10) << i+1 << flush; cout << endl << " Max: " << setw(10) << setprecision(7) << fMaxQuantity << flush; for (i = 0; i < fNVariables; i++) cout << " " << setw(10) << setprecision(4) << fMaxVariables(i) << flush; cout << endl << " Min: " << setw(10) << setprecision(7) << fMinQuantity << flush; for (i = 0; i < fNVariables; i++) cout << " " << setw(10) << setprecision(4) << fMinVariables(i) << flush; cout << endl << " Mean: " << setw(10) << setprecision(7) << fMeanQuantity << flush; for (i = 0; i < fNVariables; i++) cout << " " << setw(10) << setprecision(4) << fMeanVariables(i) << flush; cout << endl << " Function Sum Squares: " << fSumSqQuantity << endl << endl; #else for (i = 0; i < fNVariables; i++) { fprintf(stdout," %d10",i+1); fflush(stdout); } fprintf(stdout,"\n Max: %g10.7",fMaxQuantity); fflush(stdout); for (i = 0; i < fNVariables; i++) { fprintf(stdout," %g10.4",fMaxVariables(i)); fflush(stdout); } fprintf(stdout,"\n Min: %g10.7",fMinQuantity); fflush(stdout); for (i = 0; i < fNVariables; i++) { fprintf(stdout," %g10.4",fMinVariables(i)); fflush(stdout); } fprintf(stdout,"\n Mean: %g10.7",fMeanQuantity); fflush(stdout); for (i = 0; i < fNVariables; i++) { fprintf(stdout," %g10.4",fMeanVariables(i)); fflush(stdout); } fprintf(stdout,"\n Function Sum Squares: %g10.4\n\n", fSumSqQuantity); #endif } if (opt.Contains("r")) { cout << "Results of Parameterisation:" << endl << "----------------------------" << endl << " Total reduction of square residuals " << fSumSqResidual << endl << " Relative precision obtained: " << fPrecision << endl << " Error obtained: " << fError << endl << " Multiple correlation coefficient: " << fCorrelationCoeff << endl << " Reduced Chi square over sample: " << fChi2 / (fSampleSize - fNCoefficients) << endl << " Maximum residual value: " << fMaxResidual << endl << " Minimum residual value: " << fMinResidual << endl << " Estimated root mean square: " << fRMS << endl << " Maximum powers used: " << flush; for (j = 0; j < fNVariables; j++) cout << fMaxPowersFinal[j] << " " << flush; cout << endl << " Function codes of candidate functions." << endl << " 1: considered," << " 2: too little contribution," << " 3: accepted." << flush; for (i = 0; i < fMaxFunctions; i++) { if (i % 60 == 0) cout << endl << " " << flush; else if (i % 10 == 0) cout << " " << flush; cout << fFunctionCodes[i]; } cout << endl << " Loop over candidates stopped because " << flush; switch(fParameterisationCode){ case PARAM_MAXSTUDY: cout << "max allowed studies reached" << endl; break; case PARAM_SEVERAL: cout << "all candidates considered several times" << endl; break; case PARAM_RELERR: cout << "wanted relative error obtained" << endl; break; case PARAM_MAXTERMS: cout << "max number of terms reached" << endl; break; default: cout << "some unknown reason" << endl; break; } cout << endl; } if (opt.Contains("f")) { cout << "Results of Fit:" << endl << "---------------" << endl << " Test sample size: " << fTestSampleSize << endl << " Multiple correlation coefficient: " << fTestCorrelationCoeff << endl << " Relative precision obtained: " << fTestPrecision << endl << " Error obtained: " << fTestError << endl << " Reduced Chi square over sample: " << fChi2 / (fSampleSize - fNCoefficients) << endl << endl; if (fFitter) { fFitter->PrintResults(1,1); cout << endl; } } if (opt.Contains("c")){ cout << "Coefficients:" << endl << "-------------" << endl << " # Value Error Powers" << endl << " ---------------------------------------" << endl; for (i = 0; i < fNCoefficients; i++) { #ifndef R__MACOSX cout << " " << setw(3) << i << " " << setw(12) << fCoefficients(i) << " " << setw(12) << fCoefficientsRMS(i) << " " << flush; for (j = 0; j < fNVariables; j++) cout << " " << setw(3) << fPowers[fPowerIndex[i] * fNVariables + j] - 1 << flush; #else fprintf(stdout," %d3 %g12.4 %g12.4 ", i,fCoefficients(i),fCoefficientsRMS(i)); fflush(stdout); for (j = 0; j < fNVariables; j++) { fprintf(stdout," %d3",fPowers[fPowerIndex[i] * fNVariables + j] - 1); fflush(stdout); } #endif cout << endl; } cout << endl; } if (opt.Contains("k") && fCorrelationMatrix.IsValid()) { cout << "Correlation Matrix:" << endl << "-------------------"; fCorrelationMatrix.Print(); } if (opt.Contains("m")) { cout << "Parameterization:" << endl << "-----------------" << endl << " Normalised variables: " << endl; for (i = 0; i < fNVariables; i++) cout << "\ty_" << i << "\t= 1 + 2 * (x_" << i << " - " << fMaxVariables(i) << ") / (" << fMaxVariables(i) << " - " << fMinVariables(i) << ")" << endl; cout << endl << " f("; for (i = 0; i < fNVariables; i++) { cout << "y_" << i; if (i != fNVariables-1) cout << ", "; } cout << ") = "; for (Int_t i = 0; i < fNCoefficients; i++) { if (i != 0) cout << endl << "\t" << (fCoefficients(i) < 0 ? "- " : "+ ") << TMath::Abs(fCoefficients(i)); else cout << fCoefficients(i); for (Int_t j = 0; j < fNVariables; j++) { Int_t p = fPowers[fPowerIndex[i] * fNVariables + j]; switch (p) { case 1: break; case 2: cout << " * y_" << j; break; default: switch(fPolyType) { case kLegendre: cout << " * L_" << p << "(y_" << j << ")"; break; case kChebyshev: cout << " * C_" << p << "(y_" << j << ")"; break; default: cout << " * y_" << j << "^" << p; break; } } } } cout << endl; } } //____________________________________________________________________ Bool_t TMultiDimFit::Select(const Int_t *) { // Selection method. User can override this method for specialized // selection of acceptable functions in fit. Default is to select // all. This message is sent during the build-up of the function // candidates table once for each set of powers in // variables. Notice, that the argument array contains the powers // PLUS ONE. For example, to De select the function // f = x1^2 * x2^4 * x3^5, // this method should return kFALSE if given the argument // { 3, 4, 6 } return kTRUE; } //____________________________________________________________________ void TMultiDimFit::SetMaxAngle(Double_t ang) { // Set the max angle (in degrees) between the initial data vector to // be fitted, and the new candidate function to be included in the // fit. By default it is 0, which automatically chooses another // selection criteria. See also // class description

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

if (ang > 90 || ang <= 0) { Warning("SetMinAngle", "angle must be in [0,90)"); return; } fMinAngle = ang; } //____________________________________________________________________ void TMultiDimFit::SetPowers(const Int_t* powers, Int_t terms) { // Define a user function. The input array must be of the form // (p11, ..., p1N, ... ,pL1, ..., pLN) // Where N is the dimension of the data sample, L is the number of // terms (given in terms) and the first number, labels the term, the // second the variable. More information is given in the // class description

fIsUserFunction = kTRUE; fMaxFunctions = terms; fMaxTerms = terms; fMaxStudy = terms; fPowers = new Int_t[fMaxFunctions * fNVariables]; Int_t i, j; for (i = 0; i < fMaxFunctions; i++) for(j = 0; j < fNVariables; j++) fPowers[i * fNVariables + j] = powers[i * fNVariables + j] + 1; } //____________________________________________________________________ void TMultiDimFit::SetPowerLimit(Double_t limit) { // Set the user parameter for the function selection. The bigger the // limit, the more functions are used. The meaning of this variable // is defined in the // class description

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

if (!powers) return; for (Int_t i = 0; i < fNVariables; i++) fMaxPowers[i] = powers[i]+1; } //____________________________________________________________________ void TMultiDimFit::SetMinRelativeError(Double_t error) { // Set the acceptable relative error for when sum of square // residuals is considered minimized. For a full account, refer to // the // class description

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

if (fNCoefficients != 0) { // Now for the second test: if (fMaxAngle == 0) { // If the user hasn't supplied a max angle do the test as, if (dResidur < squareResidual / (fMaxTerms - fNCoefficients + 1 + 1E-10)) { return kFALSE; } } else { // If the user has provided a max angle, test if the calculated // angle is less then the max angle. if (TMath::Sqrt(dResidur/fSumSqAvgQuantity) < TMath::Cos(fMaxAngle*DEGRAD)) { return kFALSE; } } } // If we get here, the function is OK return kTRUE; } //____________________________________________________________________ void mdfHelper(int& /*npar*/, double* /*divs*/, double& chi2, double* coeffs, int /*flag*/) { // Helper function for doing the minimisation of Chi2 using Minuit // Get pointer to current TMultiDimFit object. TMultiDimFit* mdf = TMultiDimFit::Instance(); chi2 = mdf->MakeChi2(coeffs); }