ROOT logo

From $ROOTSYS/tutorials/fit/multidimfit.C

//Multi-Dimensional Parametrisation and Fitting
//Authors: Rene Brun, Christian Holm Christensen
   
#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"


//____________________________________________________________________
void makeData(Double_t* x, Double_t& d, Double_t& e) 
{
  // Make data points 
  Double_t upp[5] = { 10, 10, 10, 10,  1 };
  Double_t 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]);
}

//____________________________________________________________________
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_t nc = fit->GetNCoefficients();
  Int_t nv = fit->GetNVariables();
  const Int_t *powers = fit->GetPowers();
  const Int_t *pindex = fit->GetPowerIndex();
  if (nc != 21) return 1;
  const TVectorD *coeffs = fit->GetCoefficients();
  int k = 0;
  for (Int_t 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_t 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_t refMDF = (doFit) ? 43.95 : 43.98;
  // this does not work in CLing or ACLIC since the function is not defined
#ifdef __CINT__
  Double_t x[]    = {5,5,5,5};
  Double_t rMDF   = MDF(x);
#else
  //LM:  need to return the address of the result since it is casted to a long (this should not be in a tutorial !)
  Long_t iret = gROOT->ProcessLine(" Double_t x[] = {5,5,5,5}; double result=MDF(x); &result;");
  Double_t rMDF = * ( (Double_t*)iret);
#endif
  //printf("%f\n",rMDF);
  if (TMath::Abs(rMDF -refMDF) > 1e-2) return 4;
  return 0;     
}

//____________________________________________________________________
Int_t multidimfit(bool doFit = true) 
{

  cout << "*************************************************" << endl; 
  cout << "*             Multidimensional Fit              *" << endl;
  cout << "*                                               *" << endl;
  cout << "* By Christian Holm <cholm@nbi.dk> 14/10/00     *" << endl;
  cout << "*************************************************" << endl; 
  cout << endl;

  // Initialize global TRannom object. 
  gRandom = new TRandom();

  // Open output file 
  TFile* output = new TFile("mdf.root", "RECREATE");
  
  // Global data parameters 
  Int_t nVars       = 4;
  Int_t nData       = 500;
  Double_t x[4];

  // make fit object and set parameters on it. 
  TMultiDimFit* fit = new TMultiDimFit(nVars, TMultiDimFit::kMonomials,"v");

  Int_t 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 
  Double_t d;
  Double_t e;
  
  // Print out the start parameters
  fit->Print("p");

  printf("======================================\n");

  // Create training sample 
  Int_t 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
  fit->Print("s");

  // Book histograms 
  fit->MakeHistograms();

  // Find the parameterization 
  fit->FindParameterization();

  // Print coefficents 
  fit->Print("rc");

  // Get the min and max of variables from the training sample, used
  // for cuts in test sample. 
  Double_t *xMax = new Double_t[nVars];
  Double_t *xMin = new Double_t[nVars];
  for (i = 0; i < nVars; i++) {
    xMax[i] = (*fit->GetMaxVariables())(i);
    xMin[i] = (*fit->GetMinVariables())(i);
  }

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

  // Create test sample 
  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. 
  if (doFit) 
     fit->Fit("M");

  // Print result 
  fit->Print("fc v");

  // Write code to file 
  fit->MakeCode();

  // Write histograms to disk, and close file 
  output->Write();
  output->Close();
  delete output;
  
  // Compare results with reference run
  Int_t compare = CompareResults(fit, doFit);
  if (!compare) {
     printf("\nmultidimfit ..............................................  OK\n");
  } else {
     printf("\nmultidimfit ..............................................  fails case %d\n",compare);
  }

  // We're done 
  delete fit;
  return compare;
}
 multidimfit.C:1
 multidimfit.C:2
 multidimfit.C:3
 multidimfit.C:4
 multidimfit.C:5
 multidimfit.C:6
 multidimfit.C:7
 multidimfit.C:8
 multidimfit.C:9
 multidimfit.C:10
 multidimfit.C:11
 multidimfit.C:12
 multidimfit.C:13
 multidimfit.C:14
 multidimfit.C:15
 multidimfit.C:16
 multidimfit.C:17
 multidimfit.C:18
 multidimfit.C:19
 multidimfit.C:20
 multidimfit.C:21
 multidimfit.C:22
 multidimfit.C:23
 multidimfit.C:24
 multidimfit.C:25
 multidimfit.C:26
 multidimfit.C:27
 multidimfit.C:28
 multidimfit.C:29
 multidimfit.C:30
 multidimfit.C:31
 multidimfit.C:32
 multidimfit.C:33
 multidimfit.C:34
 multidimfit.C:35
 multidimfit.C:36
 multidimfit.C:37
 multidimfit.C:38
 multidimfit.C:39
 multidimfit.C:40
 multidimfit.C:41
 multidimfit.C:42
 multidimfit.C:43
 multidimfit.C:44
 multidimfit.C:45
 multidimfit.C:46
 multidimfit.C:47
 multidimfit.C:48
 multidimfit.C:49
 multidimfit.C:50
 multidimfit.C:51
 multidimfit.C:52
 multidimfit.C:53
 multidimfit.C:54
 multidimfit.C:55
 multidimfit.C:56
 multidimfit.C:57
 multidimfit.C:58
 multidimfit.C:59
 multidimfit.C:60
 multidimfit.C:61
 multidimfit.C:62
 multidimfit.C:63
 multidimfit.C:64
 multidimfit.C:65
 multidimfit.C:66
 multidimfit.C:67
 multidimfit.C:68
 multidimfit.C:69
 multidimfit.C:70
 multidimfit.C:71
 multidimfit.C:72
 multidimfit.C:73
 multidimfit.C:74
 multidimfit.C:75
 multidimfit.C:76
 multidimfit.C:77
 multidimfit.C:78
 multidimfit.C:79
 multidimfit.C:80
 multidimfit.C:81
 multidimfit.C:82
 multidimfit.C:83
 multidimfit.C:84
 multidimfit.C:85
 multidimfit.C:86
 multidimfit.C:87
 multidimfit.C:88
 multidimfit.C:89
 multidimfit.C:90
 multidimfit.C:91
 multidimfit.C:92
 multidimfit.C:93
 multidimfit.C:94
 multidimfit.C:95
 multidimfit.C:96
 multidimfit.C:97
 multidimfit.C:98
 multidimfit.C:99
 multidimfit.C:100
 multidimfit.C:101
 multidimfit.C:102
 multidimfit.C:103
 multidimfit.C:104
 multidimfit.C:105
 multidimfit.C:106
 multidimfit.C:107
 multidimfit.C:108
 multidimfit.C:109
 multidimfit.C:110
 multidimfit.C:111
 multidimfit.C:112
 multidimfit.C:113
 multidimfit.C:114
 multidimfit.C:115
 multidimfit.C:116
 multidimfit.C:117
 multidimfit.C:118
 multidimfit.C:119
 multidimfit.C:120
 multidimfit.C:121
 multidimfit.C:122
 multidimfit.C:123
 multidimfit.C:124
 multidimfit.C:125
 multidimfit.C:126
 multidimfit.C:127
 multidimfit.C:128
 multidimfit.C:129
 multidimfit.C:130
 multidimfit.C:131
 multidimfit.C:132
 multidimfit.C:133
 multidimfit.C:134
 multidimfit.C:135
 multidimfit.C:136
 multidimfit.C:137
 multidimfit.C:138
 multidimfit.C:139
 multidimfit.C:140
 multidimfit.C:141
 multidimfit.C:142
 multidimfit.C:143
 multidimfit.C:144
 multidimfit.C:145
 multidimfit.C:146
 multidimfit.C:147
 multidimfit.C:148
 multidimfit.C:149
 multidimfit.C:150
 multidimfit.C:151
 multidimfit.C:152
 multidimfit.C:153
 multidimfit.C:154
 multidimfit.C:155
 multidimfit.C:156
 multidimfit.C:157
 multidimfit.C:158
 multidimfit.C:159
 multidimfit.C:160
 multidimfit.C:161
 multidimfit.C:162
 multidimfit.C:163
 multidimfit.C:164
 multidimfit.C:165
 multidimfit.C:166
 multidimfit.C:167
 multidimfit.C:168
 multidimfit.C:169
 multidimfit.C:170
 multidimfit.C:171
 multidimfit.C:172
 multidimfit.C:173
 multidimfit.C:174
 multidimfit.C:175
 multidimfit.C:176
 multidimfit.C:177
 multidimfit.C:178
 multidimfit.C:179
 multidimfit.C:180
 multidimfit.C:181
 multidimfit.C:182
 multidimfit.C:183
 multidimfit.C:184
 multidimfit.C:185
 multidimfit.C:186
 multidimfit.C:187
 multidimfit.C:188
 multidimfit.C:189
 multidimfit.C:190
 multidimfit.C:191
 multidimfit.C:192
 multidimfit.C:193
 multidimfit.C:194
 multidimfit.C:195
 multidimfit.C:196
 multidimfit.C:197
 multidimfit.C:198
 multidimfit.C:199
 multidimfit.C:200
 multidimfit.C:201
 multidimfit.C:202
 multidimfit.C:203
 multidimfit.C:204
 multidimfit.C:205
 multidimfit.C:206
 multidimfit.C:207
 multidimfit.C:208
 multidimfit.C:209
 multidimfit.C:210
 multidimfit.C:211
 multidimfit.C:212
 multidimfit.C:213
 multidimfit.C:214
 multidimfit.C:215
 multidimfit.C:216
 multidimfit.C:217
 multidimfit.C:218
 multidimfit.C:219
 multidimfit.C:220
 multidimfit.C:221
 multidimfit.C:222
 multidimfit.C:223
 multidimfit.C:224
 multidimfit.C:225
 multidimfit.C:226
 multidimfit.C:227
 multidimfit.C:228
 multidimfit.C:229
 multidimfit.C:230
 multidimfit.C:231
 multidimfit.C:232
 multidimfit.C:233
 multidimfit.C:234
 multidimfit.C:235
 multidimfit.C:236
 multidimfit.C:237
 multidimfit.C:238
 multidimfit.C:239
 multidimfit.C:240
 multidimfit.C:241
 multidimfit.C:242
 multidimfit.C:243
 multidimfit.C:244
 multidimfit.C:245
 multidimfit.C:246
 multidimfit.C:247
 multidimfit.C:248
 multidimfit.C:249
 multidimfit.C:250
 multidimfit.C:251
 multidimfit.C:252
 multidimfit.C:253
 multidimfit.C:254
 multidimfit.C:255
 multidimfit.C:256
 multidimfit.C:257
 multidimfit.C:258
 multidimfit.C:259
 multidimfit.C:260
 multidimfit.C:261
 multidimfit.C:262
 multidimfit.C:263
 multidimfit.C:264
 multidimfit.C:265
 multidimfit.C:266
 multidimfit.C:267
 multidimfit.C:268
 multidimfit.C:269
 multidimfit.C:270
 multidimfit.C:271
 multidimfit.C:272
 multidimfit.C:273
 multidimfit.C:274
 multidimfit.C:275