ROOT   6.10/09 Reference Guide
line3Dfit.C
Go to the documentation of this file.
1 /// \file
2 /// \ingroup tutorial_fit
3 /// \notebook
4 /// Fitting of a TGraph2D with a 3D straight line
5 ///
6 /// run this macro by doing:
7 ///
8 /// ~~~{.cpp}
9 /// root>.x line3Dfit.C+
10 /// ~~~
11 ///
12 /// \macro_image
13 /// \macro_output
14 /// \macro_code
15 ///
16 /// \author Lorenzo Moneta
17
18 #include <TMath.h>
19 #include <TGraph2D.h>
20 #include <TRandom2.h>
21 #include <TStyle.h>
22 #include <TCanvas.h>
23 #include <TF2.h>
24 #include <TH1.h>
25 #include <Math/Functor.h>
26 #include <TPolyLine3D.h>
27 #include <Math/Vector3D.h>
28 #include <Fit/Fitter.h>
29
30 using namespace ROOT::Math;
31
32
33 // define the parametric line equation
34 void line(double t, const double *p, double &x, double &y, double &z) {
35  // a parametric line is define from 6 parameters but 4 are independent
36  // x0,y0,z0,z1,y1,z1 which are the coordinates of two points on the line
37  // can choose z0 = 0 if line not parallel to x-y plane and z1 = 1;
38  x = p[0] + p[1]*t;
39  y = p[2] + p[3]*t;
40  z = t;
41 }
42
43
44 bool first = true;
45
46 // function Object to be minimized
47 struct SumDistance2 {
48  // the TGraph is a data member of the object
49  TGraph2D *fGraph;
50
51  SumDistance2(TGraph2D *g) : fGraph(g) {}
52
53  // calculate distance line-point
54  double distance2(double x,double y,double z, const double *p) {
55  // distance line point is D= | (xp-x0) cross ux |
56  // where ux is direction of line and x0 is a point in the line (like t = 0)
57  XYZVector xp(x,y,z);
58  XYZVector x0(p[0], p[2], 0. );
59  XYZVector x1(p[0] + p[1], p[2] + p[3], 1. );
60  XYZVector u = (x1-x0).Unit();
61  double d2 = ((xp-x0).Cross(u)).Mag2();
62  return d2;
63  }
64
65  // implementation of the function to be minimized
66  double operator() (const double *par) {
67  assert(fGraph != 0);
68  double * x = fGraph->GetX();
69  double * y = fGraph->GetY();
70  double * z = fGraph->GetZ();
71  int npoints = fGraph->GetN();
72  double sum = 0;
73  for (int i = 0; i < npoints; ++i) {
74  double d = distance2(x[i],y[i],z[i],par);
75  sum += d;
76  }
77  if (first) {
78  std::cout << "Total Initial distance square = " << sum << std::endl;
79  }
80  first = false;
81  return sum;
82  }
83
84 };
85
86 Int_t line3Dfit()
87 {
88  gStyle->SetOptStat(0);
89  gStyle->SetOptFit();
90
91
92  //double e = 0.1;
93  Int_t nd = 10000;
94
95
96  // double xmin = 0; double ymin = 0;
97  // double xmax = 10; double ymax = 10;
98
99  TGraph2D * gr = new TGraph2D();
100
101  // Fill the 2D graph
102  double p0[4] = {10,20,1,2};
103
104  // generate graph with the 3d points
105  for (Int_t N=0; N<nd; N++) {
106  double x,y,z = 0;
107  // Generate a random number
108  double t = gRandom->Uniform(0,10);
109  line(t,p0,x,y,z);
110  double err = 1;
111  // do a gaussian smearing around the points in all coordinates
112  x += gRandom->Gaus(0,err);
113  y += gRandom->Gaus(0,err);
114  z += gRandom->Gaus(0,err);
115  gr->SetPoint(N,x,y,z);
116  //dt->SetPointError(N,0,0,err);
117  }
118  // fit the graph now
119
120  ROOT::Fit::Fitter fitter;
121
122
123  // make the functor objet
124  SumDistance2 sdist(gr);
125  ROOT::Math::Functor fcn(sdist,4);
126  // set the function and the initial parameter values
127  double pStart[4] = {1,1,1,1};
128  fitter.SetFCN(fcn,pStart);
129  // set step sizes different than default ones (0.3 times parameter values)
130  for (int i = 0; i < 4; ++i) fitter.Config().ParSettings(i).SetStepSize(0.01);
131
132  bool ok = fitter.FitFCN();
133  if (!ok) {
134  Error("line3Dfit","Line3D Fit failed");
135  return 1;
136  }
137
138  const ROOT::Fit::FitResult & result = fitter.Result();
139
140  std::cout << "Total final distance square " << result.MinFcnValue() << std::endl;
141  result.Print(std::cout);
142
143
144  gr->Draw("p0");
145
146  // get fit parameters
147  const double * parFit = result.GetParams();
148
149  // draw the fitted line
150  int n = 1000;
151  double t0 = 0;
152  double dt = 10;
153  TPolyLine3D *l = new TPolyLine3D(n);
154  for (int i = 0; i <n;++i) {
155  double t = t0+ dt*i/n;
156  double x,y,z;
157  line(t,parFit,x,y,z);
158  l->SetPoint(i,x,y,z);
159  }
160  l->SetLineColor(kRed);
161  l->Draw("same");
162
163  // draw original line
164  TPolyLine3D *l0 = new TPolyLine3D(n);
165  for (int i = 0; i <n;++i) {
166  double t = t0+ dt*i/n;
167  double x,y,z;
168  line(t,p0,x,y,z);
169  l0->SetPoint(i,x,y,z);
170  }
171  l0->SetLineColor(kBlue);
172  l0->Draw("same");
173  return 0;
174 }
175
176 int main() {
177  return line3Dfit();
178 }
double par[1]
Definition: unuranDistr.cxx:38
static long int sum(long int i)
Definition: Factory.cxx:2162
virtual void Draw(Option_t *option="")
Specific drawing options can be used to paint a TGraph2D:
Definition: TGraph2D.cxx:704
TLine * line
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:235
Definition: Rtypes.h:56
Documentation for class Functor class.
Definition: Functor.h:392
R__EXTERN TStyle * gStyle
Definition: TStyle.h:402
bool FitFCN(unsigned int npar, Function &fcn, const double *params=0, unsigned int dataSize=0, bool chi2fit=false)
Fit using the a generic FCN function as a C++ callable object implementing double () (const double *)...
Definition: Fitter.h:518
#define N
int Int_t
Definition: RtypesCore.h:41
A 3-dimensional polyline.
Definition: TPolyLine3D.h:31
double MinFcnValue() const
Return value of the objective function (chi2 or likelihood) used in the fit.
Definition: FitResult.h:121
TRObject operator()(const T1 &t1) const
const ParameterSettings & ParSettings(unsigned int i) const
get the parameter settings for the i-th parameter (const method)
Definition: FitConfig.h:74
const FitResult & Result() const
get fit result
Definition: Fitter.h:337
Double_t x[n]
Definition: legend1.C:17
const FitConfig & Config() const
access to the fit configuration (const method)
Definition: Fitter.h:365
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition: TAttLine.h:40
Class describing a generic displacement vector in 3 dimensions.
void SetStepSize(double err)
set the step size
Fitter class, entry point for performing all type of fits.
Definition: Fitter.h:77
void SetOptFit(Int_t fit=1)
The type of information about fit parameters printed in the histogram statistics box can be selected ...
Definition: TStyle.cxx:1219
const double * GetParams() const
parameter values (return const pointer)
Definition: FitResult.h:177
TLine * l
Definition: textangle.C:4
T Mag2(const SVector< T, D > &rhs)
Vector magnitude square Template to compute .
Definition: Functions.h:229
R__EXTERN TRandom * gRandom
Definition: TRandom.h:62
TGraphErrors * gr
Definition: legend1.C:25
Double_t * GetY() const
Definition: TGraph2D.h:119
virtual void Draw(Option_t *option="")
Draw this 3-D polyline with its current attributes.
SVector< T, 3 > Cross(const SVector< T, 3 > &lhs, const SVector< T, 3 > &rhs)
Vector Cross Product (only for 3-dim vectors) .
Definition: Functions.h:322
class containg the result of the fit and all the related information (fitted parameter values...
Definition: FitResult.h:48
virtual void SetPoint(Int_t point, Double_t x, Double_t y, Double_t z)
Set point n to x, y, z.
static const double x1[5]
virtual void SetPoint(Int_t point, Double_t x, Double_t y, Double_t z)
Sets point number n.
Definition: TGraph2D.cxx:1692
Double_t y[n]
Definition: legend1.C:17
Double_t * GetZ() const
Definition: TGraph2D.h:120
virtual Double_t Uniform(Double_t x1=1)
Returns a uniform deviate on the interval (0, x1).
Definition: TRandom.cxx:606
you should not use this method at all Int_t Int_t z
Definition: TRolke.cxx:630
bool SetFCN(unsigned int npar, Function &fcn, const double *params=0, unsigned int dataSize=0, bool chi2fit=false)
Set a generic FCN function as a C++ callable object implementing double () (const double *) Note that...
Definition: Fitter.h:523
void Print(std::ostream &os, bool covmat=false) const
print the result and optionaly covariance matrix and correlations
Definition: FitResult.cxx:446
Int_t GetN() const
Definition: TGraph2D.h:117
void SetOptStat(Int_t stat=1)
The type of information printed in the histogram statistics box can be selected via the parameter mod...
Definition: TStyle.cxx:1267
double result[121]
Definition: Rtypes.h:56
Definition: first.py:1
SVector< T, D > Unit(const SVector< T, D > &rhs)
Unit.
Definition: Functions.h:381
Graphics object made of three arrays X, Y and Z with the same number of points each.
Definition: TGraph2D.h:40
const Int_t n
Definition: legend1.C:16
int main(int argc, char **argv)
void Error(ErrorHandler_t func, int code, const char *va_(fmt),...)
Write error message and call a handler, if required.
Double_t * GetX() const
Definition: TGraph2D.h:118