multidimfit

Multi-Dimensional Parametrisation and Fitting

Author: Rene Brun, Christian Holm Christensen
This notebook tutorial was automatically generated with ROOTBOOK-izer from the macro found in the ROOT repository on Sunday, October 02, 2022 at 09:30 AM.

In [1]:
%%cpp -d
#include "Riostream.h"
#include "TROOT.h"
#include "TApplication.h"
#include "TCanvas.h"
#include "TH1.h"
#include "TSystem.h"
#include "TBrowser.h"
#include "TFile.h"
#include "TRandom.h"
#include "TMultiDimFit.h"
#include "TVectorD.h"
#include "TMath.h"

In [2]:
%%cpp -d
void makeData(double* x, double& d, double& e)
{
  // Make data points
  double upp[5] = { 10, 10, 10, 10,  1 };
  double low[5] = {  0,  0,  0,  0, .1 };
  for (int i = 0; i < 4; i++)
    x[i] = (upp[i] - low[i]) * gRandom->Rndm() + low[i];

  d = x[0] * TMath::Sqrt(x[1] * x[1] + x[2] * x[2] + x[3] * x[3]);

  e = gRandom->Gaus(upp[4],low[4]);
}

In [3]:
%%cpp -d
int CompareResults(TMultiDimFit *fit, bool doFit)
{
      //Compare results with reference run


      // the right coefficients (before fit)
   double GoodCoeffsNoFit[] = {
   -4.37056,
   43.1468,
   13.432,
   13.4632,
   13.3964,
   13.328,
   13.3016,
   13.3519,
   4.49724,
   4.63876,
   4.89036,
   -3.69982,
   -3.98618,
   -3.86195,
   4.36054,
   -4.02597,
   4.57037,
   4.69845,
   2.83819,
   -3.48855,
   -3.97612
   };

      // the right coefficients (after fit)
   double GoodCoeffs[] = {
      -4.399,
      43.15,
      13.41,
      13.49,
      13.4,
      13.23,
      13.34,
      13.29,
      4.523,
      4.659,
      4.948,
      -4.026,
      -4.045,
      -3.939,
      4.421,
      -4.006,
      4.626,
      4.378,
      3.516,
      -4.111,
      -3.823,
   };

   // Good Powers
   int GoodPower[] = {
   1,  1,  1,  1,
   2,  1,  1,  1,
   1,  1,  1,  2,
   1,  1,  2,  1,
   1,  2,  1,  1,
   2,  2,  1,  1,
   2,  1,  1,  2,
   2,  1,  2,  1,
   1,  1,  1,  3,
   1,  3,  1,  1,
   1,  1,  5,  1,
   1,  1,  2,  2,
   1,  2,  1,  2,
   1,  2,  2,  1,
   2,  1,  1,  3,
   2,  2,  1,  2,
   2,  1,  3,  1,
   2,  3,  1,  1,
   1,  2,  2,  2,
   2,  1,  2,  2,
   2,  2,  2,  1
   };

   int nc = fit->GetNCoefficients();
   int nv = fit->GetNVariables();
   const int *powers = fit->GetPowers();
   const int *pindex = fit->GetPowerIndex();
   if (nc != 21) return 1;
   const TVectorD *coeffs = fit->GetCoefficients();
   int k = 0;
   for (int i=0;i<nc;i++) {
      if (doFit) {
         if (!TMath::AreEqualRel((*coeffs)[i],GoodCoeffs[i],1e-3)) return 2;
      }
      else {
         if (TMath::Abs((*coeffs)[i] - GoodCoeffsNoFit[i]) > 5e-5) return 2;
      }
      for (int j=0;j<nv;j++) {
         if (powers[pindex[i]*nv+j] != GoodPower[k]) return 3;
         k++;
      }
   }

   // now test the result of the generated function
   gROOT->ProcessLine(".L MDF.C");

   double refMDF = (doFit) ? 43.95 : 43.98;
   // this does not work in CLing since the function is not defined
   //double x[]    = {5,5,5,5};
   //double rMDF   = MDF(x);
   //LM:  need to return the address of the result since it is casted to a long (this should not be in a tutorial !)
   std::intptr_t iret = gROOT->ProcessLine(" double xvalues[] = {5,5,5,5}; double result=MDF(xvalues); &result;");
   double rMDF = * ( (double*)iret);
   //printf("%f\n",rMDF);
   if (TMath::Abs(rMDF -refMDF) > 1e-2) return 4;
   return 0;
}

Arguments are defined.

In [4]:
bool doFit = true;
In [5]:
std::cout << "*************************************************" << std::endl;
std::cout << "*             Multidimensional Fit              *" << std::endl;
std::cout << "*                                               *" << std::endl;
std::cout << "* By Christian Holm <[email protected]> 14/10/00     *" << std::endl;
std::cout << "*************************************************" << std::endl;
std::cout << std::endl;
*************************************************
*             Multidimensional Fit              *
*                                               *
* By Christian Holm <[email protected]> 14/10/00     *
*************************************************

Initialize global TRannom object.

In [6]:
gRandom = new TRandom();

Open output file

In [7]:
TFile* output = new TFile("mdf.root", "RECREATE");

Global data parameters

In [8]:
int nVars       = 4;
int nData       = 500;
double x[4];

make fit object and set parameters on it.

In [9]:
TMultiDimFit* fit = new TMultiDimFit(nVars, TMultiDimFit::kMonomials,"v");

int mPowers[]   = { 6 , 6, 6, 6 };
fit->SetMaxPowers(mPowers);
fit->SetMaxFunctions(1000);
fit->SetMaxStudy(1000);
fit->SetMaxTerms(30);
fit->SetPowerLimit(1);
fit->SetMinAngle(10);
fit->SetMaxAngle(10);
fit->SetMinRelativeError(.01);

variables to hold the temporary input data

In [10]:
double d;
double e;

Print out the start parameters

In [11]:
fit->Print("p");

printf("======================================\n");
User parameters:
----------------
 Variables:                    4
 Data points:                  0
 Max Terms:                    30
 Power Limit Parameter:        1
 Max functions:                1000
 Max functions to study:       1000
 Max angle (optional):         10
 Min angle:                    10
 Relative Error accepted:      0.01
 Maximum Powers:                6 6 6 6

 Parameterisation will be done using Monomials

======================================

Create training sample

In [12]:
int i;
for (i = 0; i < nData ; i++) {

   // Make some data
   makeData(x,d,e);

   // Add the row to the fit object
   fit->AddRow(x,d,e);
}

Print out the statistics

In [13]:
fit->Print("s");
Sample statistics:
------------------
                 D          1          2          3          4
 Max:     141.6264      9.954       9.99      9.998      9.995
 Min:     0.149448     0.0455    0.01523    0.04109   0.003819
 Mean:    48.40441      5.033      5.044          5      5.002
 Function Sum Squares:         1.678e+06

Book histograms

In [14]:
fit->MakeHistograms();

Find the parameterization

In [15]:
fit->FindParameterization();
Coeff   SumSqRes    Contrib   Angle      QM   Func     Value        W^2  Powers
    1  5.065e+05  1.367e-26      10 6.67e-07     0  -5.23e-15        500  0 0 0 0
    2   1.15e+05  3.915e+05      50   0.167     2      47.33      174.7  1 0 0 0
    3  8.755e+04  2.749e+04      80   0.167     1      13.26      156.3  0 0 0 1
    4  6.188e+04  2.568e+04      80   0.167     3      12.39      167.3  0 0 1 0
    5  3.708e+04   2.48e+04      80   0.167     4       12.6      156.3  0 1 0 0
    6  2.596e+04  1.112e+04      85   0.333     8      14.91      50.03  1 1 0 0
    7  1.667e+04       9290      85   0.333     9      13.02      54.78  1 0 0 1
    8       7382       9287      85   0.333    14      12.64      58.13  1 0 1 0
    9       6235       1147    87.5   0.333     5      5.095      44.16  0 0 0 2
   10       5218       1018    87.5   0.333    12      4.983      40.99  0 2 0 0
   11       4193       1025    87.5   0.667    53      5.229       37.5  0 0 4 0
   12       3299      893.8    88.8   0.333     6     -4.058      54.27  0 0 1 1
   13       2458      841.2    88.8   0.333     7     -4.155      48.73  0 1 0 1
   14       1933      524.7    88.8   0.333    13     -3.291      48.45  0 1 1 0
   15       1675      258.1    88.8     0.5    19      4.211      14.56  1 0 0 2
   16       1334      340.6    88.8     0.5    26     -4.731      15.22  1 1 0 1
   17       1079      255.5    88.8     0.5    33      3.953      16.35  1 0 2 0
   18      788.2      290.4    88.8     0.5    34      4.687      13.22  1 2 0 0
   19      709.2      78.94    89.4     0.5    21       2.23      15.88  0 1 1 1
   20      473.4      235.8    89.4     0.5    23     -3.543      18.78  1 0 1 1
   21      235.4        238    89.4     0.5    28     -3.976      15.06  1 1 1 0

Print coefficents

In [16]:
fit->Print("rc");
Results of Parameterisation:
----------------------------
 Total reduction of square residuals    5.063e+05
 Relative precision obtained:           0.01185
 Error obtained:                        235.4
 Multiple correlation coefficient:      0.9995
 Reduced Chi square over sample:        0.4975
 Maximum residual value:                3.243
 Minimum residual value:                -2.59
 Estimated root mean square:            0.6862
 Maximum powers used:                   1 2 4 2 
 Function codes of candidate functions.
  1: considered,  2: too little contribution,  3: accepted.
 3333333333 1133311113 1313113131 1113311111 1111111111 1113111111
 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111
 111111
 Loop over candidates stopped because max allowed studies reached

Coefficients:
-------------
   #         Value        Error   Powers
 ---------------------------------------
   0        -4.371       0.08798     0   0   0   0
   1         43.15        0.1601     1   0   0   0
   2         13.43       0.08032     0   0   0   1
   3         13.46       0.07805     0   0   1   0
   4          13.4       0.08054     0   1   0   0
   5         13.33        0.1423     1   1   0   0
   6          13.3        0.1367     1   0   0   1
   7         13.35        0.1331     1   0   1   0
   8         4.497        0.1511     0   0   0   2
   9         4.639        0.1585     0   2   0   0
  10          4.89         0.164     0   0   4   0
  11          -3.7        0.1364     0   0   1   1
  12        -3.986        0.1438     0   1   0   1
  13        -3.862        0.1458     0   1   1   0
  14         4.361        0.2614     1   0   0   2
  15        -4.026        0.2555     1   1   0   1
  16          4.57        0.2477     1   0   2   0
  17         4.698        0.2729     1   2   0   0
  18         2.838        0.2525     0   1   1   1
  19        -3.489        0.2292     1   0   1   1
  20        -3.976        0.2566     1   1   1   0

Get the min and max of variables from the training sample, used for cuts in test sample.

In [17]:
double *xMax = new double[nVars];
double *xMin = new double[nVars];
for (i = 0; i < nVars; i++) {
   xMax[i] = (*fit->GetMaxVariables())(i);
   xMin[i] = (*fit->GetMinVariables())(i);
}

nData = fit->GetNCoefficients() * 100;
int j;

Create test sample

In [18]:
for (i = 0; i < nData ; i++) {
   // Make some data
   makeData(x,d,e);

   for (j = 0; j < nVars; j++)
      if (x[j] < xMin[j] || x[j] > xMax[j])
   break;

   // If we get through the loop above, all variables are in range
   if (j == nVars)
      // Add the row to the fit object
      fit->AddTestRow(x,d,e);
   else
      i--;
}

delete gRandom;

Test the parameterizatio and coefficents using the test sample.

In [19]:
if (doFit)
   fit->Fit("M");

Print result

In [20]:
fit->Print("fc v");
Results of Fit:
---------------
 Test sample size:                      2100
 Multiple correlation coefficient:      0.9994
 Relative precision obtained:           0.0001753
 Error obtained:                        1275
 Reduced Chi square over sample:        2.47

 FCN=1 FROM MIGRAD    STATUS=CONVERGED     861 CALLS         862 TOTAL
                     EDM=1.67352e-06    STRATEGY= 1      ERROR MATRIX ACCURATE 
  EXT PARAMETER                                  PHYSICAL LIMITS       
  NO.   NAME      VALUE            ERROR       NEGATIVE      POSITIVE  
   1  coeff00     -4.39851e+00   4.44260e-02
   2  coeff01      4.31493e+01   8.56451e-02
   3  coeff02      1.34121e+01   3.78565e-02
   4  coeff03      1.34869e+01   3.80951e-02
   5  coeff04      1.33954e+01   3.74054e-02
   6  coeff05      1.32280e+01   6.57916e-02
   7  coeff06      1.33441e+01   6.75855e-02
   8  coeff07      1.32943e+01   6.66410e-02
   9  coeff08      4.52254e+00   7.39945e-02
  10  coeff09      4.65912e+00   7.21745e-02
  11  coeff10      4.94808e+00   8.14935e-02
  12  coeff11     -4.02586e+00   6.53780e-02
  13  coeff12     -4.04534e+00   6.55396e-02
  14  coeff13     -3.93856e+00   6.51725e-02
  15  coeff14      4.42141e+00   1.30526e-01
  16  coeff15     -4.00581e+00   1.17191e-01
  17  coeff16      4.62595e+00   1.30233e-01
  18  coeff17      4.37782e+00   1.28579e-01
  19  coeff18      3.51629e+00   1.13771e-01
  20  coeff19     -4.11068e+00   1.17446e-01
  21  coeff20     -3.82302e+00   1.16486e-01

Coefficients:
-------------
   #         Value        Error   Powers
 ---------------------------------------
   0        -4.399       0.04443     0   0   0   0
   1         43.15       0.08565     1   0   0   0
   2         13.41       0.03786     0   0   0   1
   3         13.49        0.0381     0   0   1   0
   4          13.4       0.03741     0   1   0   0
   5         13.23       0.06579     1   1   0   0
   6         13.34       0.06759     1   0   0   1
   7         13.29       0.06664     1   0   1   0
   8         4.523       0.07399     0   0   0   2
   9         4.659       0.07217     0   2   0   0
  10         4.948       0.08149     0   0   4   0
  11        -4.026       0.06538     0   0   1   1
  12        -4.045       0.06554     0   1   0   1
  13        -3.939       0.06517     0   1   1   0
  14         4.421        0.1305     1   0   0   2
  15        -4.006        0.1172     1   1   0   1
  16         4.626        0.1302     1   0   2   0
  17         4.378        0.1286     1   2   0   0
  18         3.516        0.1138     0   1   1   1
  19        -4.111        0.1174     1   0   1   1
  20        -3.823        0.1165     1   1   1   0

Write code to file

In [21]:
fit->MakeCode();
Writing on file "MDF.C" ... done

Write histograms to disk, and close file

In [22]:
output->Write();
output->Close();
delete output;

Compare results with reference run

In [23]:
int compare = CompareResults(fit, doFit);
if (!compare) {
   printf("\nmultidimfit ..............................................  OK\n");
} else {
   printf("\nmultidimfit ..............................................  fails case %d\n",compare);
}
multidimfit ..............................................  OK

We're done

In [24]:
delete fit;
delete [] xMin;
delete [] xMax;
return compare;