Logo ROOT   6.12/07
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 #include <cassert>
31 
32 using namespace ROOT::Math;
33 
34 
35 // define the parametric line equation
36 void line(double t, const double *p, double &x, double &y, double &z) {
37  // a parametric line is define from 6 parameters but 4 are independent
38  // x0,y0,z0,z1,y1,z1 which are the coordinates of two points on the line
39  // can choose z0 = 0 if line not parallel to x-y plane and z1 = 1;
40  x = p[0] + p[1]*t;
41  y = p[2] + p[3]*t;
42  z = t;
43 }
44 
45 
46 bool first = true;
47 
48 // function Object to be minimized
49 struct SumDistance2 {
50  // the TGraph is a data member of the object
51  TGraph2D *fGraph;
52 
53  SumDistance2(TGraph2D *g) : fGraph(g) {}
54 
55  // calculate distance line-point
56  double distance2(double x,double y,double z, const double *p) {
57  // distance line point is D= | (xp-x0) cross ux |
58  // where ux is direction of line and x0 is a point in the line (like t = 0)
59  XYZVector xp(x,y,z);
60  XYZVector x0(p[0], p[2], 0. );
61  XYZVector x1(p[0] + p[1], p[2] + p[3], 1. );
62  XYZVector u = (x1-x0).Unit();
63  double d2 = ((xp-x0).Cross(u)).Mag2();
64  return d2;
65  }
66 
67  // implementation of the function to be minimized
68  double operator() (const double *par) {
69  assert(fGraph != 0);
70  double * x = fGraph->GetX();
71  double * y = fGraph->GetY();
72  double * z = fGraph->GetZ();
73  int npoints = fGraph->GetN();
74  double sum = 0;
75  for (int i = 0; i < npoints; ++i) {
76  double d = distance2(x[i],y[i],z[i],par);
77  sum += d;
78  }
79  if (first) {
80  std::cout << "Total Initial distance square = " << sum << std::endl;
81  }
82  first = false;
83  return sum;
84  }
85 
86 };
87 
88 Int_t line3Dfit()
89 {
90  gStyle->SetOptStat(0);
91  gStyle->SetOptFit();
92 
93 
94  //double e = 0.1;
95  Int_t nd = 10000;
96 
97 
98  // double xmin = 0; double ymin = 0;
99  // double xmax = 10; double ymax = 10;
100 
101  TGraph2D * gr = new TGraph2D();
102 
103  // Fill the 2D graph
104  double p0[4] = {10,20,1,2};
105 
106  // generate graph with the 3d points
107  for (Int_t N=0; N<nd; N++) {
108  double x,y,z = 0;
109  // Generate a random number
110  double t = gRandom->Uniform(0,10);
111  line(t,p0,x,y,z);
112  double err = 1;
113  // do a gaussian smearing around the points in all coordinates
114  x += gRandom->Gaus(0,err);
115  y += gRandom->Gaus(0,err);
116  z += gRandom->Gaus(0,err);
117  gr->SetPoint(N,x,y,z);
118  //dt->SetPointError(N,0,0,err);
119  }
120  // fit the graph now
121 
122  ROOT::Fit::Fitter fitter;
123 
124 
125  // make the functor objet
126  SumDistance2 sdist(gr);
127  ROOT::Math::Functor fcn(sdist,4);
128  // set the function and the initial parameter values
129  double pStart[4] = {1,1,1,1};
130  fitter.SetFCN(fcn,pStart);
131  // set step sizes different than default ones (0.3 times parameter values)
132  for (int i = 0; i < 4; ++i) fitter.Config().ParSettings(i).SetStepSize(0.01);
133 
134  bool ok = fitter.FitFCN();
135  if (!ok) {
136  Error("line3Dfit","Line3D Fit failed");
137  return 1;
138  }
139 
140  const ROOT::Fit::FitResult & result = fitter.Result();
141 
142  std::cout << "Total final distance square " << result.MinFcnValue() << std::endl;
143  result.Print(std::cout);
144 
145 
146  gr->Draw("p0");
147 
148  // get fit parameters
149  const double * parFit = result.GetParams();
150 
151  // draw the fitted line
152  int n = 1000;
153  double t0 = 0;
154  double dt = 10;
155  TPolyLine3D *l = new TPolyLine3D(n);
156  for (int i = 0; i <n;++i) {
157  double t = t0+ dt*i/n;
158  double x,y,z;
159  line(t,parFit,x,y,z);
160  l->SetPoint(i,x,y,z);
161  }
162  l->SetLineColor(kRed);
163  l->Draw("same");
164 
165  // draw original line
166  TPolyLine3D *l0 = new TPolyLine3D(n);
167  for (int i = 0; i <n;++i) {
168  double t = t0+ dt*i/n;
169  double x,y,z;
170  line(t,p0,x,y,z);
171  l0->SetPoint(i,x,y,z);
172  }
173  l0->SetLineColor(kBlue);
174  l0->Draw("same");
175  return 0;
176 }
177 
178 int main() {
179  return line3Dfit();
180 }
static long int sum(long int i)
Definition: Factory.cxx:2173
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:256
Definition: Rtypes.h:59
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:590
#define N
virtual void Draw(Option_t *option="P0")
Specific drawing options can be used to paint a TGraph2D:
Definition: TGraph2D.cxx:705
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:75
const FitResult & Result() const
get fit result
Definition: Fitter.h:365
Double_t x[n]
Definition: legend1.C:17
int main(int argc, char **argv)
const FitConfig & Config() const
access to the fit configuration (const method)
Definition: Fitter.h:393
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:1218
const double * GetParams() const
parameter values (return const pointer)
Definition: FitResult.h:177
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:1696
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:627
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:595
void Print(std::ostream &os, bool covmat=false) const
print the result and optionaly covariance matrix and correlations
Definition: FitResult.cxx:446
auto * l
Definition: textangle.C:4
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:1266
Definition: Rtypes.h:59
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
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
static constexpr double g