ROOT   Reference Guide
Searching...
No Matches
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
208 }
209
210 // Print out the statistics
211 fit->Print("s");
212
213 // Book histograms
214 fit->MakeHistograms();
215
216 // Find the parameterization
217 fit->FindParameterization();
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
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:45
long Long_t
Definition RtypesCore.h:54
double Double_t
Definition RtypesCore.h:59
#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:54
Multidimensional Fits in ROOT.
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:274
virtual Double_t Rndm()
Machine independent random number generator.
Definition TRandom.cxx:552
Double_t x[n]
Definition legend1.C:17