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* x, double& d, double& e)
27{
28 // Make data points
29 double upp[5] = { 10, 10, 10, 10, 1 };
30 double 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 nc = fit->GetNCoefficients();
121 int nv = fit->GetNVariables();
122 const int *powers = fit->GetPowers();
123 const int *pindex = fit->GetPowerIndex();
124 if (nc != 21) return 1;
125 const TVectorD *coeffs = fit->GetCoefficients();
126 int k = 0;
127 for (int 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 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 refMDF = (doFit) ? 43.95 : 43.98;
144 // this does not work in CLing since the function is not defined
145 //double x[] = {5,5,5,5};
146 //double 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 std::intptr_t iret = gROOT->ProcessLine(" double xvalues[] = {5,5,5,5}; double result=MDF(xvalues); &result;");
149 double rMDF = * ( (double*)iret);
150 //printf("%f\n",rMDF);
151 if (TMath::Abs(rMDF -refMDF) > 1e-2) return 4;
152 return 0;
153}
154
155//____________________________________________________________________
156int multidimfit(bool doFit = true)
157{
158
159 std::cout << "*************************************************" << std::endl;
160 std::cout << "* Multidimensional Fit *" << std::endl;
161 std::cout << "* *" << std::endl;
162 std::cout << "* By Christian Holm <cholm@nbi.dk> 14/10/00 *" << std::endl;
163 std::cout << "*************************************************" << std::endl;
164 std::cout << std::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 nVars = 4;
174 int nData = 500;
175 double x[4];
176
177 // make fit object and set parameters on it.
179
180 int 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 d;
192 double e;
193
194 // Print out the start parameters
195 fit->Print("p");
196
197 printf("======================================\n");
198
199 // Create training sample
200 int 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 *xMax = new double[nVars];
225 double *xMin = new double[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 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 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
#define gROOT
Definition TROOT.h:407
R__EXTERN TRandom * gRandom
Definition TRandom.h:62
A ROOT file is composed of a header, followed by consecutive data records (TKey instances) with a wel...
Definition TFile.h:53
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
Double_t Rndm() override
Machine independent random number generator.
Definition TRandom.cxx:552
Double_t x[n]
Definition legend1.C:17