StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
TPolinom.cxx
1 #include <stdlib.h>
2 #include <math.h>
3 #include "TError.h"
4 #include "TArrayD.h"
5 #if ROOT_VERSION_CODE < 331013
6 #include "TCL.h"
7 #else
8 #include "TCernLib.h"
9 #endif
10 #include "TPolinom.h"
11 #include <cassert>
12 
13 ClassImp(TPolinom)
14 //_____________________________________________________________________________
15 TPolinom::TPolinom(int npw,const double *coefs)
16 {
17  fNP=0; fCoe=0,fEmx=0;
18  SetCoefs(npw,coefs);
19 }
20 //_____________________________________________________________________________
21 TPolinom::TPolinom(const TPolinom& from )
22 {
23  fNP=0; fCoe=0;fEmx=0;
24  *this = from;
25 }
26 //_____________________________________________________________________________
27 TPolinom &TPolinom::operator=(const TPolinom& from )
28 {
29  Clear();
30  SetCoefs(from.fNP,from.fCoe);
31  if (!from.fEmx) return *this;
32  int n = (fNP+1)*(fNP+2)/2;
33  fEmx= new double[n];
34  TCL::ucopy(from.fEmx,fEmx,n);
35  return *this;
36 }
37 //_____________________________________________________________________________
38 TPolinom::TPolinom(double c0)
39 {
40  fNP=0; fCoe=0,fEmx=0;
41  SetCoefs(0,&c0);
42 }
43 //_____________________________________________________________________________
44 TPolinom::TPolinom(double c0,double c1)
45 {
46  fNP=0; fCoe=0,fEmx=0;
47  SetCoefs(1,0);
48  SetCoeff(0,c0);
49  SetCoeff(1,c1);
50 }
51 //_____________________________________________________________________________
52 TPolinom::TPolinom(double c0,double c1,double c2)
53 {
54  fNP=0; fCoe=0,fEmx=0;
55  SetCoefs(2,0);
56  SetCoeff(0,c0);
57  SetCoeff(1,c1);
58  SetCoeff(2,c2);
59 }
60 //_____________________________________________________________________________
61 TPolinom::~TPolinom()
62 {
63  Clear();
64 }
65 //_____________________________________________________________________________
66 void TPolinom::Clear(const char *)
67 {
68  if (fCoe!=f2Coe) delete [] fCoe; fCoe=0;
69  delete [] fEmx; fEmx=0;
70  fNP=-1;
71 }
72 //_____________________________________________________________________________
73 void TPolinom::Print(const char*) const
74 {
75  Info("Print","Power %d ",fNP);
76  if (!fCoe) return;
77  double *e = fEmx;
78  for (int i=0,n=1;i<=fNP;i++,e+=n++) {
79  double err = (fEmx)? sqrt(e[i]):0;
80  Info("Print","Coef[%d] = %g +- %g",i,fCoe[i],err);
81  }
82 }
83 
84 
85 //_____________________________________________________________________________
86 void TPolinom::SetCoefs(int npw,const double *coefs)
87 {
88  if (fNP!=npw || !fCoe) {
89  fNP = npw;
90  if (fCoe !=f2Coe) delete [] fCoe;
91  fCoe=f2Coe;
92  if(fNP>2) {
93  fCoe = new double[fNP+1];
94  memset(fCoe,0,(fNP+1)*sizeof(double));
95  }
96  }
97  if (coefs) {memcpy(fCoe,coefs,(fNP+1)*sizeof(double));}
98  else {memset(fCoe,0 ,(fNP+1)*sizeof(double));}
99 }
100 //_____________________________________________________________________________
101 double TPolinom::Eval(double x,int n,double *coe)
102 {
103  double res = 0;
104  for (int i = n; i>=0; i--) {res = coe[i]+x*res;}
105  return res;
106 }
107 //_____________________________________________________________________________
108 double TPolinom::Eval(double x) const
109 {
110  return Eval(x,fNP,fCoe);
111 }
112 //_____________________________________________________________________________
113 double TPolinom::Deriv(double x) const
114 {
115  double res = 0;
116  for (int i = fNP; i>0; i--) {res = fCoe[i]*i+x*res;}
117  return res;
118 }
119 //_____________________________________________________________________________
120 void TPolinom::Move(double x)
121 {
122  if (!fNP) return;
123  int il = fNP+1;
124  if (fNP==1) {
125  fCoe[0]+=fCoe[1]*x;
126  if (!fEmx) return;
127  fEmx[0]+=x*(2*fEmx[1]+x*fEmx[2]);
128  fEmx[1]+=x*fEmx[2];
129  return;
130  }
131 
132 
133  TArrayD arr((il*(il*3+1))/2);
134  double *tra = arr.GetArray();
135  double *tmp = tra+il*il;
136  tra[0]=1; for (int i=1;i<il;i++) tra[i]=tra[i-1]*x;
137  double *t=tra;
138  int div=1,bak=il+1;
139  for (int irow=1;irow<il;irow++) {
140  div *=irow; t+=il;
141  for (int icol=irow;icol<il;icol++) {
142  t[icol] = (t[icol-bak]*icol)/div;
143  } }
144 
145  TCL::vmatl(tra,fCoe,tmp ,il);
146  TCL::ucopy(tmp, fCoe,il);
147  if (!fEmx) return;
148  TCL::trasat(tra,fEmx,tmp,il,il);
149  TCL::ucopy(tmp, fEmx ,il*(il+1)/2);
150 }
151 //_____________________________________________________________________________
152 void TPolinom::Backward()
153 {
154  for (int i=1;i<=fNP;i+=2) { fCoe[i]=-fCoe[i];}
155  if (!fEmx) return;
156  for (int i=0,li=0;i<=fNP;li+=++i) {
157  for (int j=0 ;j<=i ;j++) {if ((i^j)&1) fEmx[li+j]*=-1;}}
158 }
159 //_____________________________________________________________________________
160 double TPolinom::GetEmx(int i,int j) const
161 {
162  if (!fEmx) return 0;
163  int ii = i,jj=j;
164  if (i<j) {ii=j;jj=i;}
165  return fEmx[((ii+1)*ii)/2+jj];
166 }
167 //_____________________________________________________________________________
168 double TPolinom::Evrr(double x) const
169 {
170  if (!fEmx) return 0;
171  TArrayD arr(fNP+1);
172  double *xx = arr.GetArray();
173  xx[0]=1;
174  for (int i=1;i<=fNP;i++) xx[i]=xx[i-1]*x;
175  double res;
176  TCL::trasat(xx,fEmx,&res,1,fNP+1);
177  return res;
178 }
179 //_____________________________________________________________________________
180 //_____________________________________________________________________________
181 //_____________________________________________________________________________
182 //_____________________________________________________________________________
183 ClassImp(TPoliFitter)
184 //_____________________________________________________________________________
185 enum EPoliFitter {kX,kY,kW,kXYW=3};
186 
187 //_____________________________________________________________________________
188 TPoliFitter::TPoliFitter(int np):TPolinom(np),fArr(100)
189 {
190  Clear();
191 }
192 //_____________________________________________________________________________
193 TPoliFitter::TPoliFitter(const TPoliFitter &fr):TPolinom(fr),fArr(fr.fArr)
194 {
195  memcpy(fBeg,fr.fBeg,fEnd-fBeg+1);
196  fDat = fArr.GetArray();
197  fC = fDat+fN; fP = fC+fNP+1;
198 }
199 //_____________________________________________________________________________
200 void TPoliFitter::Add(double x, double y,double err2)
201 {
202  if (fArr.GetSize()<=fN+kXYW) { fArr.Set(fN*2+kXYW); fDat = fArr.GetArray();}
203  fArr[fN+kX]=x;
204  fArr[fN+kY]=y;
205  fArr[fN+kW]=1./err2;
206  fN+=kXYW;fNuse++;
207 }
208 //_____________________________________________________________________________
209 void TPoliFitter::AddErr(double err2)
210 {
211  fArr[fN-kXYW+kW]=1./err2;
212 }
213 //_____________________________________________________________________________
214 const double *TPoliFitter::GetX(int i) const
215 {return fDat+i*kXYW;}
216 //_____________________________________________________________________________
217  double *TPoliFitter::GetX(int i)
218 {return fDat+i*kXYW;}
219 //_____________________________________________________________________________
220 void TPoliFitter::Clear(const char *)
221 {
222  memset(fBeg,0,fEnd-fBeg+1);
223  fArr.Reset();
224  fDat = fArr.GetArray();
225 }
226 //_____________________________________________________________________________
227 void TPoliFitter::Print(const char* tit) const
228 {
229  if (!tit || !*tit) tit = "Print";
230 
231  Info(tit,"NPoints %d Wtot=%g",fN/kXYW,fWtot);
232  for (int l=0;l<fN;l+=kXYW) {
233  double dy = fDat[l+kY]-Eval(fDat[l+kX]);
234  double dXi2 = dy*dy*fDat[l+kW]*fWtot;
235  printf("%d - \tX=%g \tY=%g \tW=%g \tdY=%g \tdXi2=%g\n"
236  ,l/kXYW,fDat[l+kX],fDat[l+kY],fDat[l+kW]
237  ,dy,dXi2);
238  }
239  TPolinom::Print();
240 }
241 //_____________________________________________________________________________
242 void TPoliFitter::Prepare()
243 {
244  int need = fN+(fNP+1) + (fNP+1)*(fNP+2)/2;
245  if (need > fArr.GetSize()) fArr.Set(need);
246  fDat = fArr.GetArray();
247  fC = fDat+fN; fP = fC+fNP+1;
248  TCL::vzero(fC,(fNP+1) + (fNP+1)*(fNP+2)/2);
249  if (!fWtot) {
250  for (int l=0;l<fN;l+=kXYW) {fWtot+=fDat[l+kW];}
251  for (int l=0;l<fN;l+=kXYW) {fDat[l+kW]/=fWtot;}
252  }
253  fP[0]=1;
254  int lNew = 1,lLst=0,pLst=0,lOld,pOld;
255  for (int pNew=1;pNew<=fNP;) {
256  TCL::vzero(fC,pNew+1);
257  for (int l=0;l<fN;l+=kXYW) {//Data loop
258  double x = fDat[l+kX],wt=fDat[l+kW];
259  if (wt<=0) continue;
260  double yNew = Eval(x,pLst,fP+lLst)*x;
261  fC[pNew]+=yNew*yNew*wt;
262  yNew*=wt;
263  for (lOld = lLst,pOld=pLst;pOld>=0;lOld-=pOld,pOld--) {//Old loop
264  fC[pOld]+= yNew * TPolinom::Eval(x,pOld,fP+lOld);
265  }//end old loop
266  }//end Data loop
267 // New ort polinom
268  fP[lNew]=0;
269  TCL::ucopy(fP+lLst,fP+lNew+1,pLst+1);
270  for (lOld = lLst,pOld=pLst;pOld>=0;lOld-=pOld,pOld--) {//Old loop
271  TCL::vlinco(fP+lNew,1.,fP+lOld,-fC[pOld],fP+lNew,pOld+1);
272  }
273  double nor = fC[pNew]-TCL::vdot(fC,fC,pNew);
274  nor = sqrt(nor);
275  TCL::vscale(fP+lNew,1./nor,fP+lNew,pNew+1);
276  lLst = lNew; pLst=pNew;
277  lNew += pNew+1; pNew++;
278  }
279 }
280 //_____________________________________________________________________________
281 double TPoliFitter::Fit()
282 {
283  Prepare();
284  int lp,np;
285 
286  TCL::vzero(fC,fNP+1);
287  fChi2=0;
288  for (int l=0;l<fN;l+=kXYW) {//Data loop
289  double x = fDat[l+kX];
290  double y = fDat[l+kY];
291  double w = fDat[l+kW];
292  if (w<=0) continue;
293  fChi2 += y*y*w;
294  for (lp=0,np=0;np<=fNP; lp+=++np ) {
295  fC[np]+=Eval(x,np,fP+lp)*y*w;
296  }
297  }
298  TCL::vzero(fCoe,fNP+1);
299  for (lp=0,np=0;np<=fNP; np++, lp+=np ) {
300  fChi2-=fC[np]*fC[np];
301  TCL::vlinco(fCoe,1.,fP+lp,fC[np],fCoe,np+1);
302  }
303  fChi2*=fWtot;
304  fNdf = fNuse-(fNP+1);
305  if (fNdf>0) fChi2/=fNdf;
306  return fChi2;
307 }
308 //_____________________________________________________________________________
309 void TPoliFitter::MakeErrs()
310 {
311  delete [] fEmx;
312  int n = (fNP+1)*(fNP+2)/2;
313  fEmx = new double[n];
314  TCL::trsmul(fP,fEmx,fNP+1);
315  TCL::vscale(fEmx,1./fWtot,fEmx,n);
316 }
317 //_____________________________________________________________________________
318 double TPoliFitter::EvalChi2()
319 {
320  double Xi2=0; int n=0;
321  for (int l=0;l<fN;l+=kXYW) {//Data loop
322  double x = fDat[l+kX];
323  double y = fDat[l+kY];
324  double w = fDat[l+kW];
325  if (w<=0) continue;
326  double ye = Eval(x);
327  Xi2+= (y-ye)*(y-ye)*w;
328  n++;
329  }
330  Xi2*=fWtot;
331  if (fNdf) Xi2/=fNdf;
332  return Xi2;
333 }
334 //_____________________________________________________________________________
335 double TPoliFitter::FixAt(double x, double y)
336 {
337  assert(fEmx);
338  if (fabs(x)>1e-8) Move(x);
339 
340  double emx0 = fEmx[0];
341  double h = y-fCoe[0];
342  double lamda = h/emx0;
343  double *e = fEmx;
344  for (int i=0,n=1;i<=fNP;i++,e+=n++) {fCoe[i]+= e[0]*lamda;}
345  TCL::trsinv(fEmx,fEmx,fNP+1);
346  e = fEmx;
347  for (int i=0,n=1;i<=fNP;i++,e+=n++) {e[0]=0;}
348  fEmx[0]=1;
349  TCL::trsinv(fEmx,fEmx,fNP+1);
350  fEmx[0]=0;
351  fNdf++;
352  double addXi2 = h*h/emx0;
353  fChi2 += (addXi2-fChi2)/fNdf;
354  if (fabs(x)>1e-8) Move(-x);
355  return fChi2;
356 }
357 //_____________________________________________________________________________
358 void TPoliFitter::Skip(int idx)
359 {
360  fDat[kXYW*idx+kW]=-1;
361  SetNdf(fNdf-1);
362 }
363 //_____________________________________________________________________________
364 void TPoliFitter::SetNdf(int ndf)
365 {
366  fChi2*=fNdf;
367  if (ndf) fChi2/=ndf;
368  fNdf=ndf;
369 }
370 //_____________________________________________________________________________
371 void TPoliFitter::DCoeDy(int iy,double *dcoe) const
372 {
373 // derivative dCoe/dY[idat]
374  TCL::vzero(dcoe,fNP+1);
375  TArrayD ac(fNP+1); double *c=ac.GetArray();
376  int l = iy*kXYW;
377  double x = fDat[l+kX];
378  double w = fDat[l+kW];
379  if (w<=0) return;
380  for (int lp=0,np=0;np<=fNP; lp+=++np ) {
381  c[np]=Eval(x,np,fP+lp)*w;
382  }
383  for (int lp=0,np=0;np<=fNP; np++, lp+=np ) {
384  TCL::vlinco(dcoe,1.,fP+lp,c[np],dcoe,np+1);
385  }
386 }
387 //_____________________________________________________________________________
388 double TPoliFitter::EvalOrt(int idx,double x) const
389 {
390  int lp = (idx*(idx+1))/2;
391  return Eval(x,idx,fP+lp);
392 }
393 //_____________________________________________________________________________
394 double TPoliFitter::MakeMatrix(TMatrixD &Akj) const
395 {
396 // Create matrix such that sum(Akj*Yk*Yj) == Chi2
397 
398  int N = NPts();
399  Akj.ResizeTo(N,N);
400  for (int k=0;k<N;k++) {
401  double Xk = GetX(k)[0];
402  double Wk = GetX(k)[2];
403  for (int j=0;j<=k;j++) {
404  double Xj = GetX(j)[0];
405  double Wj = GetX(j)[2];
406  double Fkj=0;
407  for (int i=0;i<=fNP;i++) {//loop over ort polinoms
408  double Fik = EvalOrt(i,Xk);
409  double Fij = EvalOrt(i,Xj);
410  Fkj+= Fik*Fij;
411  }//end i
412  Akj[k][j] = -Fkj*Wk*Wj*fWtot;
413  Akj[j][k] = Akj[k][j];
414  }//end j
415  Akj[k][k]+= Wk*fWtot;
416  }//end k
417  return fChi2*fNdf;
418 }
419 
420 //_____________________________________________________________________________
421 //_____________________________________________________________________________
422 #include "TCanvas.h"
423 #include "TH1F.h"
424 #include "TGraph.h"
425 #include "TSystem.h"
426 #include "TRandom.h"
427 #include "TVectorD.h"
428 //_____________________________________________________________________________
429 //______________________________________________________________________________
430 void TPoliFitter::Show() const
431 {
432 static TCanvas *myCanvas = 0;
433 static TGraph *ptGraph = 0;
434 static TGraph *ciGraph = 0;
435  double x[100],y[100];
436  int nPts = Size();
437  if (nPts>100) nPts=100;
438  for (int i=0;i<nPts;i++) {
439  x[i]=fDat[i*3+0]; y[i]=fDat[i*3+1];
440  }
441 
442 
443  if(!myCanvas) myCanvas = new TCanvas("TPoliFitter","",600,800);
444  myCanvas->Clear();
445 
446  delete ptGraph; delete ciGraph;
447  ptGraph = new TGraph(nPts , x, y);
448  ptGraph->SetMarkerColor(kRed);
449  ptGraph->Draw("A*");
450 
451  double x0=x[0];
452  double dx = (x[nPts-1]-x[0])/99;
453  for (int i=0;i<100;i++) {
454  x[i]=x0+dx*i;
455  y[i]=Eval(x[i]);
456  }
457 
458  ciGraph = new TGraph(100 , x, y);
459  ciGraph->Draw("Same CP");
460 
461  myCanvas->Modified();
462  myCanvas->Update();
463  while(!gSystem->ProcessEvents()){};
464 
465 }
466 //_____________________________________________________________________________
467 void TPoliFitter::Test(int kase)
468 {
469 static TCanvas* myCanvas=0;
470 static TH1F *hh[6]={0};
471 static const char *hNams[]={"dA0","dA1","dA2","dY","Xi2","Xi2E",0};
472  if(!myCanvas) myCanvas=new TCanvas("C1","",600,800);
473  myCanvas->Clear();
474  myCanvas->Divide(1,6);
475 
476  for (int i=0;i<6;i++) {
477  int low = (i>=4)? 0:-5;
478  delete hh[i]; hh[i]= new TH1F(hNams[i],hNams[i],100,low,5);
479  myCanvas->cd(i+1); hh[i]->Draw();
480  }
481 
482 
483 
484  double A[3]={1,2,3},B[3]={1,2,3};
485  if (kase==1) {
486  B[0] = A[0]+ 5*(A[1]+5*A[2]);
487  B[1] = A[1] + 2*A[2]*5;
488  B[2] = A[2];
489  }
490 
491  double y5 = B[0] + 5.5*(B[1]+5.5*B[2]);
492  for (int ievt=0;ievt <10000; ievt++) {
493 
494  TPoliFitter pf(2);
495  int nPts=20;
496  for (double x=0;x<nPts;x++) {
497  double y = A[0]+x*(A[1]+x*A[2]);
498  double dy = gRandom->Gaus(0,0.1);
499  pf.Add(x,y+dy,0.1*0.1);
500  }
501  double Xi2 = pf.Fit();
502  pf.MakeErrs();
503  if (kase==1) { pf.Move(5) ;}
504  if (kase==2) { pf.FixAt(0.,A[0]);}
505  double Xi2E = pf.EvalChi2();
506  hh[4]->Fill(Xi2);
507  hh[5]->Fill(Xi2E);
508 
509 
510  const double *c = pf.Coe();
511 // printf("Xi2 = %g Coe= %g %g %g\n",Xi2,c[0],c[1],c[2]);
512  double del = pf.Eval(5.5)-y5;
513  del /= sqrt(pf.Evrr(5.5));
514  hh[3]->Fill(del);
515  for (int i=0;i<3;i++) {
516  del = (c[i]-B[i])/sqrt(pf.GetEmx(i,i)+1e-10);
517  hh[i]->Fill(del);
518  }
519  }
520 
521  myCanvas->Modified();
522  myCanvas->Update();
523  while(!gSystem->ProcessEvents()){};
524 }
525 //_____________________________________________________________________________
526 void TPoliFitter::TestCorr()
527 {
528 static TCanvas* myCanvas=0;
529 static TH1F *hh[6]={0,0,0,0,0,0};
530 static const char *hNams[]={"C01","C01-","C02","C02-","C12","C12-",0};
531  if(!myCanvas) myCanvas=new TCanvas("C1","",600,800);
532  myCanvas->Clear();
533  myCanvas->Divide(1,6);
534 
535  for (int i=0;i<6;i++) {
536  delete hh[i]; hh[i]= new TH1F(hNams[i],hNams[i],100,-1,1);
537  myCanvas->cd(i+1); hh[i]->Draw();
538  }
539 
540 
541 
542  double A[3]={1,2,3},B[3]={1,2,3};
543 
544  for (int ievt=0;ievt <1000; ievt++) {
545 
546  TPoliFitter pf(2);
547  for (double x=0;x<10;x++) {
548  double y = A[0]+x*(A[1]+x*A[2]);
549  double dy = gRandom->Gaus(0,0.1);
550  pf.Add(x,y+dy,0.1*0.1);
551  }
552  double Xi2 = pf.Fit(); if (Xi2){};
553  pf.MakeErrs();
554 
555  const double *c = pf.Coe();
556  int ih = 0;
557  for (int i=0;i<3;i++) {
558  double deli = (c[i]-B[i]);
559 // double diai = pf.GetEmx(i,i);
560  for (int j=i+1;j<3;j++) {
561  double delj = (c[j]-B[j]);
562 // double diaj = pf.GetEmx(j,j);
563 // hh[ih+0]->Fill(deli*delj/sqrt(diai*diaj));
564 // hh[ih+1]->Fill((pf.GetEmx(i,j))/sqrt(diai*diaj));
565  hh[ih+0]->Fill(deli*delj*100);
566  hh[ih+1]->Fill((deli*delj-pf.GetEmx(i,j))*100);
567  ih+=2;
568  } }
569  }
570 
571  myCanvas->Modified();
572  myCanvas->Update();
573  while(!gSystem->ProcessEvents()){};
574 }
575 //_____________________________________________________________________________
576 void TPoliFitter::Dest(int kase)
577 {
578 if (kase){};
579 static TCanvas* myCanvas=0;
580 static TH1F *hh[5]={0,0,0,0,0};
581 static const char *hNams[]={"Der0","Der1","Der2",0};
582  if(!myCanvas) myCanvas=new TCanvas("C1","",600,800);
583  myCanvas->Clear();
584  myCanvas->Divide(1,3);
585 
586  for (int i=0;i<3;i++) {
587  delete hh[i]; hh[i]= new TH1F(hNams[i],hNams[i],100,-1,1);
588  myCanvas->cd(i+1); hh[i]->Draw();
589  }
590 
591  double A[3]={1,2,3};
592 
593  for (int ievt=0;ievt <1000; ievt++) {
594 
595  TPoliFitter pf(2);
596  int nPts=20;
597  for (double x=0;x<nPts;x++) {
598  double y = A[0]+x*(A[1]+x*A[2]);
599  double dy = gRandom->Gaus(0,0.1);
600  pf.Add(x,y+dy,0.1*0.1);
601  }
602  double Xi2 = pf.Fit();
603  pf.MakeErrs();
604  double coe[3],doe[3];
605  TCL::ucopy(pf.Coe(),coe,3);
606 
607  for (int ip=0;ip<nPts;ip++) {
608  TPoliFitter pff(pf);
609  pf.DCoeDy(ip,doe);
610  double y=pf.GetX(ip)[1];
611  double dy = 0.1*y;
612  pff.GetX(ip)[1]=y+dy;
613  Xi2 = pff.Fit(); if(Xi2){};
614  for (int ih=0;ih<3;ih++) {
615  double tst = (pff.Coe()[ih]-coe[ih])/dy;
616  tst = (tst-doe[ih])/(fabs(doe[ih])+1e-10);
617  hh[ih]->Fill(tst);
618  }
619  }
620  }
621 
622  myCanvas->Modified();
623  myCanvas->Update();
624  while(!gSystem->ProcessEvents()){};
625 }
626 //_____________________________________________________________________________
627 void TPoliFitter::Test2()
628 {
629  enum {nPts=10};
630  double A[3]={1,2,3};
631  TMatrixD G(nPts,nPts);
632  TVectorD Y(nPts),D(nPts),Y1(nPts),D1(nPts);
633  int np = 2;
634 
635  TPoliFitter pf(np);
636  for (int ix=0;ix<nPts;ix++) {
637  double x = ix;
638  double y = A[0]+x*(A[1]+x*A[2]);
639  double err = 0.1*sqrt(ix+1.);
640  double dy = gRandom->Gaus(0,err);
641  Y[ix]=y+dy;
642  pf.Add(x,y+dy,err*err);
643  }
644  double Xi2 = pf.Fit();
645  pf.MakeErrs();
646  printf("Make Akj[][] matrix\n");
647  Xi2= pf.MakeMatrix(G);
648  double myXi2 = (Y*(G*Y));
649  D = 2.*(G*Y);
650  Y1 = Y;
651  printf("TPoliFitter::Test2() Xi2=%g myXi2=%g\n",Xi2,myXi2);
652  for (int ix=0;ix<nPts;ix++) {
653  double sav = pf.GetX(ix)[1];
654  double delta = sav*1e-3;
655  pf.GetX(ix)[1] = sav + delta;
656  Y1[ix] = sav + delta;
657  D1 = 2.*(G*Y1);
658  double Xi21 = pf.Fit()*pf.Ndf();
659  pf.GetX(ix)[1] = sav;
660  Y1[ix] = sav;
661  double num = (Xi21-Xi2)/delta;
662  double ana = 0.5*(D[ix]+D1[ix]);
663  double dif = 200*fabs(num-ana)/(fabs(num)+fabs(ana)+1e-10);
664  printf("TPoliFitter::Test2() d/d[%d] \tAna=%g \tNum=%g \tDif=%g\n"
665  ,ix,ana,num,dif);
666  }
667 
668 }
669 
670 
671 
672 
673 
674 
675 
676 
677 
TCL::trsmul
static float * trsmul(const float *g, float *gi, int n)
Definition: TCernLib.cxx:1212
TCL::trsinv
static float * trsinv(const float *g, float *gi, int n)
Definition: TCernLib.cxx:1160
TCL::trasat
static float * trasat(const float *a, const float *s, float *r, int m, int n)
Definition: TCernLib.cxx:540
Generated by   doxygen 1.8.5