Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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
32using namespace ROOT::Math;
33
34
35// define the parametric line equation
36void 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
46bool first = true;
47
48// function Object to be minimized
49struct 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
88Int_t line3Dfit()
89{
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
178int main() {
179 return line3Dfit();
180}
int main()
Definition Prototype.cxx:12
#define d(i)
Definition RSha256.hxx:102
#define g(i)
Definition RSha256.hxx:105
static const double x1[5]
int Int_t
Definition RtypesCore.h:45
@ kRed
Definition Rtypes.h:66
@ kBlue
Definition Rtypes.h:66
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
Definition TError.cxx:187
#define N
TGLVector3 Cross(const TGLVector3 &v1, const TGLVector3 &v2)
Definition TGLUtil.h:323
TRObject operator()(const T1 &t1) const
R__EXTERN TRandom * gRandom
Definition TRandom.h:62
R__EXTERN TStyle * gStyle
Definition TStyle.h:413
const ParameterSettings & ParSettings(unsigned int i) const
get the parameter settings for the i-th parameter (const method)
Definition FitConfig.h:76
class containg the result of the fit and all the related information (fitted parameter values,...
Definition FitResult.h:47
const double * GetParams() const
parameter values (return const pointer)
Definition FitResult.h:176
void Print(std::ostream &os, bool covmat=false) const
print the result and optionaly covariance matrix and correlations
double MinFcnValue() const
Return value of the objective function (chi2 or likelihood) used in the fit.
Definition FitResult.h:120
Fitter class, entry point for performing all type of fits.
Definition Fitter.h:77
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:610
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:615
const FitResult & Result() const
get fit result
Definition Fitter.h:384
const FitConfig & Config() const
access to the fit configuration (const method)
Definition Fitter.h:412
void SetStepSize(double err)
set the step size
Documentation for class Functor class.
Definition Functor.h:401
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition TAttLine.h:40
Graphics object made of three arrays X, Y and Z with the same number of points each.
Definition TGraph2D.h:41
Double_t * GetY() const
Definition TGraph2D.h:122
Double_t * GetX() const
Definition TGraph2D.h:121
Int_t GetN() const
Definition TGraph2D.h:120
Double_t * GetZ() const
Definition TGraph2D.h:123
virtual void SetPoint(Int_t i, Double_t x, Double_t y)
Set x and y values for point number i.
Definition TGraph.cxx:2298
virtual void Draw(Option_t *chopt="")
Draw this graph with its current attributes.
Definition TGraph.cxx:769
A 3-dimensional polyline.
Definition TPolyLine3D.h:33
virtual void SetPoint(Int_t point, Double_t x, Double_t y, Double_t z)
Set point n to x, y, z.
virtual void Draw(Option_t *option="")
Draw this 3-D polyline with its current attributes.
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 Uniform(Double_t x1=1)
Returns a uniform deviate on the interval (0, x1).
Definition TRandom.cxx:672
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:1589
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:1541
TLine * line
SVector< T, D > Unit(const SVector< T, D > &rhs)
Unit.
Definition Functions.h:382
T Mag2(const SVector< T, D > &rhs)
Vector magnitude square Template to compute .
Definition Functions.h:230
Double_t y[n]
Definition legend1.C:17
Double_t x[n]
Definition legend1.C:17
const Int_t n
Definition legend1.C:16
TGraphErrors * gr
Definition legend1.C:25
Definition first.py:1
auto * l
Definition textangle.C:4
static uint64_t sum(uint64_t i)
Definition Factory.cxx:2345