Logo ROOT  
Reference Guide
multidimfit.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_fit
3/// \notebook -nodraw
4/// Multi-Dimensional Parametrisation and Fitting
5///
6/// \macro_output
7/// \macro_code
8///
9/// \authors Rene Brun, Christian Holm Christensen
10
11#include "Riostream.h"
12#include "TROOT.h"
13#include "TApplication.h"
14#include "TCanvas.h"
15#include "TH1.h"
16#include "TSystem.h"
17#include "TBrowser.h"
18#include "TFile.h"
19#include "TRandom.h"
20#include "TMultiDimFit.h"
21#include "TVectorD.h"
22#include "TMath.h"
23
24
25//____________________________________________________________________
26void makeData(Double_t* x, Double_t& d, Double_t& e)
27{
28 // Make data points
29 Double_t upp[5] = { 10, 10, 10, 10, 1 };
30 Double_t low[5] = { 0, 0, 0, 0, .1 };
31 for (int i = 0; i < 4; i++)
32 x[i] = (upp[i] - low[i]) * gRandom->Rndm() + low[i];
33
34 d = x[0] * TMath::Sqrt(x[1] * x[1] + x[2] * x[2] + x[3] * x[3]);
35
36 e = gRandom->Gaus(upp[4],low[4]);
37}
38
39//____________________________________________________________________
40int CompareResults(TMultiDimFit *fit, bool doFit)
41{
42 //Compare results with reference run
43
44
45 // the right coefficients (before fit)
46 double GoodCoeffsNoFit[] = {
47 -4.37056,
48 43.1468,
49 13.432,
50 13.4632,
51 13.3964,
52 13.328,
53 13.3016,
54 13.3519,
55 4.49724,
56 4.63876,
57 4.89036,
58 -3.69982,
59 -3.98618,
60 -3.86195,
61 4.36054,
62 -4.02597,
63 4.57037,
64 4.69845,
65 2.83819,
66 -3.48855,
67 -3.97612
68 };
69
70 // the right coefficients (after fit)
71 double GoodCoeffs[] = {
72 -4.399,
73 43.15,
74 13.41,
75 13.49,
76 13.4,
77 13.23,
78 13.34,
79 13.29,
80 4.523,
81 4.659,
82 4.948,
83 -4.026,
84 -4.045,
85 -3.939,
86 4.421,
87 -4.006,
88 4.626,
89 4.378,
90 3.516,
91 -4.111,
92 -3.823,
93 };
94
95 // Good Powers
96 int GoodPower[] = {
97 1, 1, 1, 1,
98 2, 1, 1, 1,
99 1, 1, 1, 2,
100 1, 1, 2, 1,
101 1, 2, 1, 1,
102 2, 2, 1, 1,
103 2, 1, 1, 2,
104 2, 1, 2, 1,
105 1, 1, 1, 3,
106 1, 3, 1, 1,
107 1, 1, 5, 1,
108 1, 1, 2, 2,
109 1, 2, 1, 2,
110 1, 2, 2, 1,
111 2, 1, 1, 3,
112 2, 2, 1, 2,
113 2, 1, 3, 1,
114 2, 3, 1, 1,
115 1, 2, 2, 2,
116 2, 1, 2, 2,
117 2, 2, 2, 1
118 };
119
120 Int_t nc = fit->GetNCoefficients();
121 Int_t nv = fit->GetNVariables();
122 const Int_t *powers = fit->GetPowers();
123 const Int_t *pindex = fit->GetPowerIndex();
124 if (nc != 21) return 1;
125 const TVectorD *coeffs = fit->GetCoefficients();
126 int k = 0;
127 for (Int_t i=0;i<nc;i++) {
128 if (doFit) {
129 if (!TMath::AreEqualRel((*coeffs)[i],GoodCoeffs[i],1e-3)) return 2;
130 }
131 else {
132 if (TMath::Abs((*coeffs)[i] - GoodCoeffsNoFit[i]) > 5e-5) return 2;
133 }
134 for (Int_t j=0;j<nv;j++) {
135 if (powers[pindex[i]*nv+j] != GoodPower[k]) return 3;
136 k++;
137 }
138 }
139
140 // now test the result of the generated function
141 gROOT->ProcessLine(".L MDF.C");
142
143 Double_t refMDF = (doFit) ? 43.95 : 43.98;
144 // this does not work in CLing since the function is not defined
145 //Double_t x[] = {5,5,5,5};
146 //Double_t rMDF = MDF(x);
147 //LM: need to return the address of the result since it is casted to a long (this should not be in a tutorial !)
148 Long_t iret = gROOT->ProcessLine(" Double_t xvalues[] = {5,5,5,5}; double result=MDF(xvalues); &result;");
149 Double_t rMDF = * ( (Double_t*)iret);
150 //printf("%f\n",rMDF);
151 if (TMath::Abs(rMDF -refMDF) > 1e-2) return 4;
152 return 0;
153}
154
155//____________________________________________________________________
156Int_t multidimfit(bool doFit = true)
157{
158
159 cout << "*************************************************" << endl;
160 cout << "* Multidimensional Fit *" << endl;
161 cout << "* *" << endl;
162 cout << "* By Christian Holm <cholm@nbi.dk> 14/10/00 *" << endl;
163 cout << "*************************************************" << endl;
164 cout << endl;
165
166 // Initialize global TRannom object.
167 gRandom = new TRandom();
168
169 // Open output file
170 TFile* output = new TFile("mdf.root", "RECREATE");
171
172 // Global data parameters
173 Int_t nVars = 4;
174 Int_t nData = 500;
175 Double_t x[4];
176
177 // make fit object and set parameters on it.
179
180 Int_t mPowers[] = { 6 , 6, 6, 6 };
181 fit->SetMaxPowers(mPowers);
182 fit->SetMaxFunctions(1000);
183 fit->SetMaxStudy(1000);
184 fit->SetMaxTerms(30);
185 fit->SetPowerLimit(1);
186 fit->SetMinAngle(10);
187 fit->SetMaxAngle(10);
188 fit->SetMinRelativeError(.01);
189
190 // variables to hold the temporary input data
191 Double_t d;
192 Double_t e;
193
194 // Print out the start parameters
195 fit->Print("p");
196
197 printf("======================================\n");
198
199 // Create training sample
200 Int_t i;
201 for (i = 0; i < nData ; i++) {
202
203 // Make some data
204 makeData(x,d,e);
205
206 // Add the row to the fit object
207 fit->AddRow(x,d,e);
208 }
209
210 // Print out the statistics
211 fit->Print("s");
212
213 // Book histograms
214 fit->MakeHistograms();
215
216 // Find the parameterization
218
219 // Print coefficents
220 fit->Print("rc");
221
222 // Get the min and max of variables from the training sample, used
223 // for cuts in test sample.
224 Double_t *xMax = new Double_t[nVars];
225 Double_t *xMin = new Double_t[nVars];
226 for (i = 0; i < nVars; i++) {
227 xMax[i] = (*fit->GetMaxVariables())(i);
228 xMin[i] = (*fit->GetMinVariables())(i);
229 }
230
231 nData = fit->GetNCoefficients() * 100;
232 Int_t j;
233
234 // Create test sample
235 for (i = 0; i < nData ; i++) {
236 // Make some data
237 makeData(x,d,e);
238
239 for (j = 0; j < nVars; j++)
240 if (x[j] < xMin[j] || x[j] > xMax[j])
241 break;
242
243 // If we get through the loop above, all variables are in range
244 if (j == nVars)
245 // Add the row to the fit object
246 fit->AddTestRow(x,d,e);
247 else
248 i--;
249 }
250 //delete gRandom;
251
252 // Test the parameterizatio and coefficents using the test sample.
253 if (doFit)
254 fit->Fit("M");
255
256 // Print result
257 fit->Print("fc v");
258
259 // Write code to file
260 fit->MakeCode();
261
262 // Write histograms to disk, and close file
263 output->Write();
264 output->Close();
265 delete output;
266
267 // Compare results with reference run
268 Int_t compare = CompareResults(fit, doFit);
269 if (!compare) {
270 printf("\nmultidimfit .............................................. OK\n");
271 } else {
272 printf("\nmultidimfit .............................................. fails case %d\n",compare);
273 }
274
275 // We're done
276 delete fit;
277 delete [] xMin;
278 delete [] xMax;
279 return compare;
280}
#define d(i)
Definition: RSha256.hxx:102
#define e(i)
Definition: RSha256.hxx:103
int Int_t
Definition: RtypesCore.h:43
long Long_t
Definition: RtypesCore.h:52
double Double_t
Definition: RtypesCore.h:57
#define gROOT
Definition: TROOT.h:406
R__EXTERN TRandom * gRandom
Definition: TRandom.h:62
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format.
Definition: TFile.h:53
Multidimensional Fits in ROOT.
Definition: TMultiDimFit.h:15
Int_t GetNCoefficients() const
Definition: TMultiDimFit.h:162
void SetMaxStudy(Int_t n)
Definition: TMultiDimFit.h:200
void SetMaxFunctions(Int_t n)
Definition: TMultiDimFit.h:198
void SetMaxTerms(Int_t terms)
Definition: TMultiDimFit.h:201
Int_t GetNVariables() const
Definition: TMultiDimFit.h:161
virtual void MakeCode(const char *functionName="MDF", Option_t *option="")
Generate the file <filename> with .C appended if argument doesn't end in .cxx or ....
const TVectorD * GetMaxVariables() const
Definition: TMultiDimFit.h:154
void SetPowerLimit(Double_t limit=1e-3)
Set the user parameter for the function selection.
virtual void FindParameterization(Option_t *option="")
Find the parameterization.
const TVectorD * GetMinVariables() const
Definition: TMultiDimFit.h:160
virtual void AddRow(const Double_t *x, Double_t D, Double_t E=0)
Add a row consisting of fNVariables independent variables, the known, dependent quantity,...
const Int_t * GetPowers() const
Definition: TMultiDimFit.h:166
void SetMaxPowers(const Int_t *powers)
Set the maximum power to be considered in the fit for each variable.
virtual void AddTestRow(const Double_t *x, Double_t D, Double_t E=0)
Add a row consisting of fNVariables independent variables, the known, dependent quantity,...
virtual void Fit(Option_t *option="")
Try to fit the found parameterisation to the test sample.
void SetMinRelativeError(Double_t error)
Set the acceptable relative error for when sum of square residuals is considered minimized.
void SetMaxAngle(Double_t angle=0)
Set the max angle (in degrees) between the initial data vector to be fitted, and the new candidate fu...
Int_t * GetPowerIndex() const
Definition: TMultiDimFit.h:164
void SetMinAngle(Double_t angle=1)
Set the min angle (in degrees) between a new candidate function and the subspace spanned by the previ...
virtual void Print(Option_t *option="ps") const
Print statistics etc.
const TVectorD * GetCoefficients() const
Definition: TMultiDimFit.h:142
virtual void MakeHistograms(Option_t *option="A")
Make histograms of the result of the analysis.
This is the base class for the ROOT Random number generators.
Definition: TRandom.h:27
virtual Double_t Gaus(Double_t mean=0, Double_t sigma=1)
Samples a random number from the standard Normal (Gaussian) Distribution with the given mean and sigm...
Definition: TRandom.cxx:263
virtual Double_t Rndm()
Machine independent random number generator.
Definition: TRandom.cxx:541
Double_t x[n]
Definition: legend1.C:17
Double_t Sqrt(Double_t x)
Definition: TMath.h:681
Bool_t AreEqualRel(Double_t af, Double_t bf, Double_t relPrec)
Definition: TMath.h:418
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
static void output(int code)
Definition: gifencode.c:226