StRoot: StarRoot/TPolinom.cxx Source File

00001 #include <stdlib.h>
00002 #include <math.h>
00003 #include "TError.h"
00004 #include "TArrayD.h"
00005 #if ROOT_VERSION_CODE < 331013
00006 #include "TCL.h"
00007 #else
00008 #include "TCernLib.h"
00009 #endif
00010 #include "TPolinom.h"
00011 #include <cassert>
00012 
00013 ClassImp(TPolinom)
00014 //_____________________________________________________________________________
00015 TPolinom::TPolinom(int npw,const double *coefs)
00016 {
00017   fNP=0; fCoe=0,fEmx=0;
00018   SetCoefs(npw,coefs);
00019 }
00020 //_____________________________________________________________________________
00021 TPolinom::TPolinom(const TPolinom& from )
00022 {
00023   fNP=0; fCoe=0;fEmx=0;
00024   *this = from;
00025 }
00026 //_____________________________________________________________________________
00027 TPolinom &TPolinom::operator=(const TPolinom& from )
00028 {
00029   Clear();
00030   SetCoefs(from.fNP,from.fCoe);
00031   if (!from.fEmx) return *this;
00032   int n = (fNP+1)*(fNP+2)/2;
00033   fEmx= new double[n];
00034   TCL::ucopy(from.fEmx,fEmx,n);
00035   return *this;
00036 }
00037 //_____________________________________________________________________________
00038 TPolinom::TPolinom(double c0)
00039 {
00040   fNP=0; fCoe=0,fEmx=0;
00041   SetCoefs(0,&c0);
00042 }
00043 //_____________________________________________________________________________
00044 TPolinom::TPolinom(double c0,double c1)
00045 {
00046   fNP=0; fCoe=0,fEmx=0;
00047   SetCoefs(1,0);
00048   SetCoeff(0,c0);
00049   SetCoeff(1,c1);
00050 }
00051 //_____________________________________________________________________________
00052 TPolinom::TPolinom(double c0,double c1,double c2)
00053 {
00054   fNP=0; fCoe=0,fEmx=0;
00055   SetCoefs(2,0);
00056   SetCoeff(0,c0);
00057   SetCoeff(1,c1);
00058   SetCoeff(2,c2);
00059 }
00060 //_____________________________________________________________________________
00061 TPolinom::~TPolinom()
00062 {
00063   Clear();
00064 }
00065 //_____________________________________________________________________________
00066 void TPolinom::Clear(const char *)
00067 {
00068   if (fCoe!=f2Coe)  delete [] fCoe; fCoe=0;
00069   delete [] fEmx; fEmx=0;
00070   fNP=-1;
00071 }
00072 //_____________________________________________________________________________
00073 void TPolinom::Print(const char*) const
00074 {
00075   Info("Print","Power %d ",fNP);
00076   if (!fCoe) return;
00077   double *e = fEmx;
00078   for (int i=0,n=1;i<=fNP;i++,e+=n++) {
00079     double err = (fEmx)? sqrt(e[i]):0; 
00080     Info("Print","Coef[%d] = %g +- %g",i,fCoe[i],err);
00081   }
00082 }
00083 
00084 
00085 //_____________________________________________________________________________
00086 void TPolinom::SetCoefs(int npw,const double *coefs)
00087 {
00088   if (fNP!=npw || !fCoe) {
00089     fNP = npw;
00090     if (fCoe !=f2Coe) delete [] fCoe;
00091     fCoe=f2Coe;
00092     if(fNP>2) {
00093       fCoe = new double[fNP+1];
00094       memset(fCoe,0,(fNP+1)*sizeof(double));
00095     }
00096   }
00097   if (coefs) {memcpy(fCoe,coefs,(fNP+1)*sizeof(double));}
00098   else       {memset(fCoe,0    ,(fNP+1)*sizeof(double));}
00099 }
00100 //_____________________________________________________________________________
00101 double TPolinom::Eval(double x,int n,double *coe) 
00102 {
00103   double res = 0;
00104   for (int i = n; i>=0; i--) {res = coe[i]+x*res;}
00105   return res;
00106 }  
00107 //_____________________________________________________________________________
00108 double TPolinom::Eval(double x) const
00109 {
00110   return Eval(x,fNP,fCoe);
00111 }  
00112 //_____________________________________________________________________________
00113 double TPolinom::Deriv(double x) const
00114 {
00115   double res = 0;
00116   for (int i = fNP; i>0; i--) {res = fCoe[i]*i+x*res;}
00117   return res;
00118 }  
00119 //_____________________________________________________________________________
00120 void TPolinom::Move(double x)
00121 {
00122   if (!fNP) return;
00123   int il = fNP+1;
00124   if (fNP==1) {
00125     fCoe[0]+=fCoe[1]*x;
00126     if (!fEmx) return;
00127     fEmx[0]+=x*(2*fEmx[1]+x*fEmx[2]);
00128     fEmx[1]+=x*fEmx[2];
00129     return;
00130   }
00131 
00132 
00133   TArrayD arr((il*(il*3+1))/2);
00134   double *tra = arr.GetArray();
00135   double *tmp = tra+il*il;
00136   tra[0]=1; for (int i=1;i<il;i++) tra[i]=tra[i-1]*x;
00137   double *t=tra;
00138   int div=1,bak=il+1;
00139   for (int irow=1;irow<il;irow++) {
00140     div *=irow; t+=il;
00141     for (int icol=irow;icol<il;icol++) {
00142        t[icol] = (t[icol-bak]*icol)/div;
00143   } }
00144 
00145   TCL::vmatl(tra,fCoe,tmp ,il);
00146   TCL::ucopy(tmp,     fCoe,il);
00147   if (!fEmx) return;
00148   TCL::trasat(tra,fEmx,tmp,il,il);
00149   TCL::ucopy(tmp, fEmx ,il*(il+1)/2);
00150 }
00151 //_____________________________________________________________________________
00152 void TPolinom::Backward() 
00153 {
00154   for (int i=1;i<=fNP;i+=2) { fCoe[i]=-fCoe[i];}
00155   if (!fEmx) return;
00156   for (int i=0,li=0;i<=fNP;li+=++i) {
00157     for (int j=0   ;j<=i  ;j++)     {if ((i^j)&1) fEmx[li+j]*=-1;}}
00158 }
00159 //_____________________________________________________________________________
00160 double TPolinom::GetEmx(int i,int j) const
00161 {
00162   if (!fEmx) return 0;
00163   int ii = i,jj=j;
00164   if (i<j) {ii=j;jj=i;}
00165   return fEmx[((ii+1)*ii)/2+jj];
00166 }
00167 //_____________________________________________________________________________
00168 double TPolinom::Evrr(double x) const
00169 {
00170   if (!fEmx) return 0;
00171   TArrayD arr(fNP+1);
00172   double *xx = arr.GetArray();
00173   xx[0]=1;
00174   for (int i=1;i<=fNP;i++) xx[i]=xx[i-1]*x;
00175   double res;
00176   TCL::trasat(xx,fEmx,&res,1,fNP+1);
00177   return res;
00178 }  
00179 //_____________________________________________________________________________
00180 //_____________________________________________________________________________
00181 //_____________________________________________________________________________
00182 //_____________________________________________________________________________
00183 ClassImp(TPoliFitter)
00184 //_____________________________________________________________________________
00185 enum EPoliFitter {kX,kY,kW,kXYW=3};
00186 
00187 //_____________________________________________________________________________
00188 TPoliFitter::TPoliFitter(int np):TPolinom(np),fArr(100)
00189 {
00190   Clear();
00191 }
00192 //_____________________________________________________________________________
00193 TPoliFitter::TPoliFitter(const TPoliFitter &fr):TPolinom(fr),fArr(fr.fArr)
00194 {
00195   memcpy(fBeg,fr.fBeg,fEnd-fBeg+1);
00196   fDat = fArr.GetArray();
00197   fC   = fDat+fN; fP = fC+fNP+1;
00198 }
00199 //_____________________________________________________________________________
00200 void TPoliFitter::Add(double x, double y,double err2)
00201 {
00202   if (fArr.GetSize()<=fN+kXYW) { fArr.Set(fN*2+kXYW); fDat = fArr.GetArray();}
00203   fArr[fN+kX]=x;
00204   fArr[fN+kY]=y;
00205   fArr[fN+kW]=1./err2;
00206   fN+=kXYW;fNuse++;
00207 }
00208 //_____________________________________________________________________________
00209 void TPoliFitter::AddErr(double err2)
00210 {
00211   fArr[fN-kXYW+kW]=1./err2;
00212 }
00213 //_____________________________________________________________________________
00214 const double *TPoliFitter::GetX(int i) const    
00215 {return fDat+i*kXYW;}
00216 //_____________________________________________________________________________
00217       double *TPoliFitter::GetX(int i)
00218 {return fDat+i*kXYW;}
00219 //_____________________________________________________________________________
00220 void TPoliFitter::Clear(const char *)
00221 {
00222   memset(fBeg,0,fEnd-fBeg+1);
00223   fArr.Reset();
00224   fDat = fArr.GetArray();
00225 }
00226 //_____________________________________________________________________________
00227 void TPoliFitter::Print(const char* tit) const
00228 {
00229   if (!tit || !*tit) tit = "Print";
00230 
00231   Info(tit,"NPoints %d Wtot=%g",fN/kXYW,fWtot);
00232   for (int l=0;l<fN;l+=kXYW) {
00233     double dy = fDat[l+kY]-Eval(fDat[l+kX]);
00234     double dXi2 = dy*dy*fDat[l+kW]*fWtot;
00235     printf("%d - \tX=%g \tY=%g \tW=%g \tdY=%g \tdXi2=%g\n"
00236           ,l/kXYW,fDat[l+kX],fDat[l+kY],fDat[l+kW]
00237           ,dy,dXi2);
00238   }
00239   TPolinom::Print();
00240 }
00241 //_____________________________________________________________________________
00242 void TPoliFitter::Prepare()
00243 {
00244   int need = fN+(fNP+1) + (fNP+1)*(fNP+2)/2;
00245   if (need > fArr.GetSize()) fArr.Set(need);
00246   fDat = fArr.GetArray();
00247   fC   = fDat+fN; fP = fC+fNP+1;
00248   TCL::vzero(fC,(fNP+1) + (fNP+1)*(fNP+2)/2);
00249   if (!fWtot) {
00250     for (int l=0;l<fN;l+=kXYW) {fWtot+=fDat[l+kW];}
00251     for (int l=0;l<fN;l+=kXYW) {fDat[l+kW]/=fWtot;}
00252   }
00253   fP[0]=1;
00254   int lNew = 1,lLst=0,pLst=0,lOld,pOld;
00255   for (int pNew=1;pNew<=fNP;) {
00256     TCL::vzero(fC,pNew+1);
00257     for (int l=0;l<fN;l+=kXYW) {//Data loop
00258       double x = fDat[l+kX],wt=fDat[l+kW];
00259       if (wt<=0) continue;
00260       double yNew = Eval(x,pLst,fP+lLst)*x;
00261       fC[pNew]+=yNew*yNew*wt;
00262       yNew*=wt;
00263       for (lOld = lLst,pOld=pLst;pOld>=0;lOld-=pOld,pOld--) {//Old loop
00264         fC[pOld]+= yNew * TPolinom::Eval(x,pOld,fP+lOld); 
00265       }//end old loop
00266     }//end Data loop
00267 //      New ort polinom         
00268     fP[lNew]=0;
00269     TCL::ucopy(fP+lLst,fP+lNew+1,pLst+1);
00270     for (lOld = lLst,pOld=pLst;pOld>=0;lOld-=pOld,pOld--) {//Old loop
00271       TCL::vlinco(fP+lNew,1.,fP+lOld,-fC[pOld],fP+lNew,pOld+1);
00272     }
00273     double nor = fC[pNew]-TCL::vdot(fC,fC,pNew);
00274     nor = sqrt(nor);
00275     TCL::vscale(fP+lNew,1./nor,fP+lNew,pNew+1);
00276     lLst = lNew; pLst=pNew;
00277     lNew += pNew+1; pNew++;
00278   }
00279 }  
00280 //_____________________________________________________________________________
00281 double TPoliFitter::Fit()   
00282 {
00283   Prepare();
00284   int lp,np;
00285 
00286   TCL::vzero(fC,fNP+1);
00287   fChi2=0;
00288   for (int l=0;l<fN;l+=kXYW) {//Data loop
00289     double x = fDat[l+kX];
00290     double y = fDat[l+kY];
00291     double w = fDat[l+kW];
00292     if (w<=0) continue;
00293     fChi2 += y*y*w;
00294     for (lp=0,np=0;np<=fNP; lp+=++np ) {
00295       fC[np]+=Eval(x,np,fP+lp)*y*w;
00296     }
00297   }
00298   TCL::vzero(fCoe,fNP+1);
00299   for (lp=0,np=0;np<=fNP; np++, lp+=np ) {
00300     fChi2-=fC[np]*fC[np];
00301     TCL::vlinco(fCoe,1.,fP+lp,fC[np],fCoe,np+1);
00302   }
00303   fChi2*=fWtot;
00304   fNdf = fNuse-(fNP+1);
00305   if (fNdf>0) fChi2/=fNdf;
00306   return fChi2;
00307 }
00308 //_____________________________________________________________________________
00309 void   TPoliFitter::MakeErrs()
00310 {
00311   delete fEmx;
00312   int n = (fNP+1)*(fNP+2)/2;
00313   fEmx = new double[n];
00314   TCL::trsmul(fP,fEmx,fNP+1);
00315 //TCL::trsmlu(fP,fEmx,fNP+1);
00316   TCL::vscale(fEmx,1./fWtot,fEmx,n);
00317 }
00318 //_____________________________________________________________________________
00319 double TPoliFitter::EvalChi2() 
00320 {
00321   double Xi2=0; int n=0;
00322   for (int l=0;l<fN;l+=kXYW) {//Data loop
00323     double x = fDat[l+kX];
00324     double y = fDat[l+kY];
00325     double w = fDat[l+kW];
00326     if (w<=0) continue;
00327     double ye = Eval(x);
00328     Xi2+= (y-ye)*(y-ye)*w;
00329     n++;
00330   }
00331   Xi2*=fWtot;
00332   if (fNdf) Xi2/=fNdf;
00333   return Xi2;
00334 }
00335 //_____________________________________________________________________________
00336 double TPoliFitter::FixAt(double x, double y) 
00337 {
00338    assert(fEmx);
00339    if (fabs(x)>1e-8) Move(x);
00340 
00341    double emx0 = fEmx[0];
00342    double h    = y-fCoe[0];
00343    double lamda = h/emx0;
00344    double *e = fEmx;
00345    for (int i=0,n=1;i<=fNP;i++,e+=n++) {fCoe[i]+= e[0]*lamda;}
00346    TCL::trsinv(fEmx,fEmx,fNP+1);
00347    e = fEmx;
00348    for (int i=0,n=1;i<=fNP;i++,e+=n++) {e[0]=0;}
00349    fEmx[0]=1;
00350    TCL::trsinv(fEmx,fEmx,fNP+1);
00351    fEmx[0]=0;
00352    fNdf++;
00353    double addXi2 = h*h/emx0;
00354    fChi2 += (addXi2-fChi2)/fNdf;
00355    if (fabs(x)>1e-8) Move(-x);
00356    return fChi2;
00357 }
00358 //_____________________________________________________________________________
00359 void TPoliFitter::Skip(int idx) 
00360 {
00361   fDat[kXYW*idx+kW]=-1;
00362   SetNdf(fNdf-1);
00363 }
00364 //_____________________________________________________________________________
00365 void TPoliFitter::SetNdf(int ndf) 
00366 {
00367   fChi2*=fNdf;
00368   if (ndf) fChi2/=ndf;
00369   fNdf=ndf;
00370 }
00371 //_____________________________________________________________________________
00372 void TPoliFitter::DCoeDy(int iy,double *dcoe) const
00373 {
00374 //      derivative dCoe/dY[idat]
00375   TCL::vzero(dcoe,fNP+1);
00376   TArrayD ac(fNP+1); double *c=ac.GetArray();
00377   int l = iy*kXYW;
00378   double x = fDat[l+kX];
00379   double w = fDat[l+kW];
00380   if (w<=0) return;
00381   for (int lp=0,np=0;np<=fNP; lp+=++np ) {
00382     c[np]=Eval(x,np,fP+lp)*w;
00383   }
00384   for (int lp=0,np=0;np<=fNP; np++, lp+=np ) {
00385     TCL::vlinco(dcoe,1.,fP+lp,c[np],dcoe,np+1);
00386   }
00387 }
00388 //_____________________________________________________________________________
00389 double TPoliFitter::EvalOrt(int idx,double x) const
00390 {
00391   int lp = (idx*(idx+1))/2;
00392   return Eval(x,idx,fP+lp);
00393 }  
00394 //_____________________________________________________________________________
00395 double TPoliFitter::MakeMatrix(TMatrixD &Akj) const
00396 {
00397 //  Create matrix such that sum(Akj*Yk*Yj) == Chi2
00398 
00399   int N = NPts();
00400   Akj.ResizeTo(N,N); 
00401   for (int k=0;k<N;k++) {
00402     double Xk = GetX(k)[0];
00403     double Wk = GetX(k)[2];
00404     for (int j=0;j<=k;j++) {
00405       double Xj = GetX(j)[0];
00406       double Wj = GetX(j)[2];
00407       double Fkj=0;
00408       for (int i=0;i<=fNP;i++) {//loop over ort polinoms       
00409         double Fik = EvalOrt(i,Xk);
00410         double Fij = EvalOrt(i,Xj);
00411         Fkj+= Fik*Fij;
00412       }//end i
00413       Akj[k][j] = -Fkj*Wk*Wj*fWtot;
00414       Akj[j][k] = Akj[k][j];
00415      }//end j
00416      Akj[k][k]+= Wk*fWtot;
00417    }//end k
00418    return fChi2*fNdf;
00419 }
00420 
00421 //_____________________________________________________________________________
00422 //_____________________________________________________________________________
00423 #include "TCanvas.h"
00424 #include "TH1F.h"
00425 #include "TGraph.h"
00426 #include "TSystem.h"
00427 #include "TRandom.h"
00428 #include "TVectorD.h"
00429 //_____________________________________________________________________________
00430 //______________________________________________________________________________
00431 void TPoliFitter::Show() const
00432 {
00433 static TCanvas *myCanvas = 0;
00434 static TGraph  *ptGraph  = 0;
00435 static TGraph  *ciGraph  = 0;
00436   double x[100],y[100];
00437   int nPts = Size();
00438   if (nPts>100) nPts=100;
00439   for (int i=0;i<nPts;i++) {
00440     x[i]=fDat[i*3+0];  y[i]=fDat[i*3+1];
00441   }
00442 
00443 
00444   if(!myCanvas) myCanvas = new TCanvas("TPoliFitter","",600,800);
00445   myCanvas->Clear();
00446 
00447   delete ptGraph; delete ciGraph;
00448   ptGraph  = new TGraph(nPts  , x, y);
00449   ptGraph->SetMarkerColor(kRed);
00450   ptGraph->Draw("A*");
00451 
00452   double x0=x[0];
00453   double dx = (x[nPts-1]-x[0])/99;
00454   for (int i=0;i<100;i++) {
00455     x[i]=x0+dx*i;
00456     y[i]=Eval(x[i]);
00457   }
00458   
00459   ciGraph  = new TGraph(100  , x, y);
00460   ciGraph->Draw("Same CP");
00461 
00462   myCanvas->Modified();
00463   myCanvas->Update();
00464   while(!gSystem->ProcessEvents()){}; 
00465 
00466 }
00467 //_____________________________________________________________________________
00468 void TPoliFitter::Test(int kase)   
00469 {
00470 static TCanvas* myCanvas=0;
00471 static TH1F *hh[6]={0};
00472 static const char *hNams[]={"dA0","dA1","dA2","dY","Xi2","Xi2E",0};
00473   if(!myCanvas)  myCanvas=new TCanvas("C1","",600,800);
00474   myCanvas->Clear();
00475   myCanvas->Divide(1,6);
00476 
00477   for (int i=0;i<6;i++) { 
00478     int low = (i>=4)? 0:-5;
00479     delete hh[i]; hh[i]= new TH1F(hNams[i],hNams[i],100,low,5);
00480     myCanvas->cd(i+1); hh[i]->Draw();
00481   }
00482 
00483 
00484 
00485   double A[3]={1,2,3},B[3]={1,2,3};
00486   if (kase==1) {
00487     B[0] = A[0]+ 5*(A[1]+5*A[2]);
00488     B[1] = A[1] + 2*A[2]*5;
00489     B[2] = A[2];
00490   }
00491 
00492   double y5 = B[0] + 5.5*(B[1]+5.5*B[2]);
00493   for (int ievt=0;ievt <10000; ievt++) {
00494   
00495     TPoliFitter pf(2);
00496     int nPts=20;
00497     for (double x=0;x<nPts;x++) {
00498       double y = A[0]+x*(A[1]+x*A[2]);
00499       double dy = gRandom->Gaus(0,0.1);
00500       pf.Add(x,y+dy,0.1*0.1);
00501     }
00502     double Xi2 = pf.Fit();
00503     pf.MakeErrs();
00504     if (kase==1) { pf.Move(5)    ;}
00505     if (kase==2) { pf.FixAt(0.,A[0]);}
00506     double Xi2E = pf.EvalChi2();
00507     hh[4]->Fill(Xi2);
00508     hh[5]->Fill(Xi2E);
00509 
00510 
00511     const double *c = pf.Coe();
00512 //    printf("Xi2 = %g Coe= %g %g %g\n",Xi2,c[0],c[1],c[2]);
00513     double del = pf.Eval(5.5)-y5;
00514     del /= sqrt(pf.Evrr(5.5));
00515     hh[3]->Fill(del);
00516     for (int i=0;i<3;i++) {
00517       del = (c[i]-B[i])/sqrt(pf.GetEmx(i,i)+1e-10);  
00518       hh[i]->Fill(del);
00519     }
00520   }
00521 
00522   myCanvas->Modified();
00523   myCanvas->Update();
00524   while(!gSystem->ProcessEvents()){}; 
00525 }
00526 //_____________________________________________________________________________
00527 void TPoliFitter::TestCorr()   
00528 {
00529 static TCanvas* myCanvas=0;
00530 static TH1F *hh[6]={0,0,0,0,0,0};
00531 static const char *hNams[]={"C01","C01-","C02","C02-","C12","C12-",0};
00532   if(!myCanvas)  myCanvas=new TCanvas("C1","",600,800);
00533   myCanvas->Clear();
00534   myCanvas->Divide(1,6);
00535 
00536   for (int i=0;i<6;i++) { 
00537     delete hh[i]; hh[i]= new TH1F(hNams[i],hNams[i],100,-1,1);
00538     myCanvas->cd(i+1); hh[i]->Draw();
00539   }
00540 
00541 
00542 
00543   double A[3]={1,2,3},B[3]={1,2,3};
00544 
00545   for (int ievt=0;ievt <1000; ievt++) {
00546   
00547     TPoliFitter pf(2);
00548     for (double x=0;x<10;x++) {
00549       double y = A[0]+x*(A[1]+x*A[2]);
00550       double dy = gRandom->Gaus(0,0.1);
00551       pf.Add(x,y+dy,0.1*0.1);
00552     }
00553     double Xi2 = pf.Fit(); if (Xi2){};
00554     pf.MakeErrs();
00555 
00556     const double *c = pf.Coe();
00557     int ih = 0;
00558     for (int i=0;i<3;i++) {
00559       double deli = (c[i]-B[i]);
00560 //      double diai = pf.GetEmx(i,i);
00561     for (int j=i+1;j<3;j++) {
00562       double delj = (c[j]-B[j]);
00563 //      double diaj = pf.GetEmx(j,j);
00564 //       hh[ih+0]->Fill(deli*delj/sqrt(diai*diaj));
00565 //       hh[ih+1]->Fill((pf.GetEmx(i,j))/sqrt(diai*diaj));
00566       hh[ih+0]->Fill(deli*delj*100);
00567       hh[ih+1]->Fill((deli*delj-pf.GetEmx(i,j))*100);
00568       ih+=2;
00569     } }
00570   }
00571 
00572   myCanvas->Modified();
00573   myCanvas->Update();
00574   while(!gSystem->ProcessEvents()){}; 
00575 }
00576 //_____________________________________________________________________________
00577 void TPoliFitter::Dest(int kase)   
00578 {
00579 if (kase){};
00580 static TCanvas* myCanvas=0;
00581 static TH1F *hh[5]={0,0,0,0,0};
00582 static const char *hNams[]={"Der0","Der1","Der2",0};
00583   if(!myCanvas)  myCanvas=new TCanvas("C1","",600,800);
00584   myCanvas->Clear();
00585   myCanvas->Divide(1,3);
00586 
00587   for (int i=0;i<3;i++) { 
00588     delete hh[i]; hh[i]= new TH1F(hNams[i],hNams[i],100,-1,1);
00589     myCanvas->cd(i+1); hh[i]->Draw();
00590   }
00591 
00592   double A[3]={1,2,3};
00593 
00594   for (int ievt=0;ievt <1000; ievt++) {
00595   
00596     TPoliFitter pf(2);
00597     int nPts=20;
00598     for (double x=0;x<nPts;x++) {
00599       double y = A[0]+x*(A[1]+x*A[2]);
00600       double dy = gRandom->Gaus(0,0.1);
00601       pf.Add(x,y+dy,0.1*0.1);
00602     }
00603     double Xi2 = pf.Fit();
00604     pf.MakeErrs();
00605     double coe[3],doe[3];
00606     TCL::ucopy(pf.Coe(),coe,3);
00607     
00608     for (int ip=0;ip<nPts;ip++) {
00609       TPoliFitter pff(pf);
00610       pf.DCoeDy(ip,doe);
00611       double y=pf.GetX(ip)[1];
00612       double dy = 0.1*y;
00613       pff.GetX(ip)[1]=y+dy;
00614       Xi2 = pff.Fit();
00615       for (int ih=0;ih<3;ih++) {
00616         double tst = (pff.Coe()[ih]-coe[ih])/dy;
00617         tst = (tst-doe[ih])/(fabs(doe[ih])+1e-10);
00618         hh[ih]->Fill(tst);
00619       }
00620     }  
00621   }
00622 
00623   myCanvas->Modified();
00624   myCanvas->Update();
00625   while(!gSystem->ProcessEvents()){}; 
00626 }
00627 //_____________________________________________________________________________
00628 void TPoliFitter::Test2()   
00629 {
00630   enum {nPts=10};
00631   double A[3]={1,2,3};
00632   TMatrixD G(nPts,nPts);
00633   TVectorD Y(nPts),D(nPts),Y1(nPts),D1(nPts);
00634   int np = 2;
00635   
00636   TPoliFitter pf(np);
00637   for (int ix=0;ix<nPts;ix++) {
00638     double x = ix;
00639     double y = A[0]+x*(A[1]+x*A[2]);
00640     double err = 0.1*sqrt(ix+1.);
00641     double dy = gRandom->Gaus(0,err);
00642     Y[ix]=y+dy;
00643     pf.Add(x,y+dy,err*err);
00644   }
00645   double Xi2 = pf.Fit();
00646   pf.MakeErrs();
00647   printf("Make Akj[][] matrix\n");
00648   Xi2= pf.MakeMatrix(G);
00649   double myXi2 = (Y*(G*Y));
00650   D = 2.*(G*Y);
00651   Y1 = Y;
00652   printf("TPoliFitter::Test2() Xi2=%g myXi2=%g\n",Xi2,myXi2);
00653   for (int ix=0;ix<nPts;ix++) {
00654     double sav = pf.GetX(ix)[1];
00655     double delta = sav*1e-3;
00656     pf.GetX(ix)[1] = sav + delta;
00657     Y1[ix]         = sav + delta;
00658     D1 = 2.*(G*Y1);
00659     double Xi21 = pf.Fit()*pf.Ndf();
00660     pf.GetX(ix)[1] = sav;
00661     Y1[ix]         = sav;
00662     double num = (Xi21-Xi2)/delta;
00663     double ana = 0.5*(D[ix]+D1[ix]);
00664     double dif = 200*fabs(num-ana)/(fabs(num)+fabs(ana)+1e-10);
00665     printf("TPoliFitter::Test2() d/d[%d] \tAna=%g \tNum=%g \tDif=%g\n"
00666             ,ix,ana,num,dif);
00667   }
00668 
00669 }
00670 
00671 
00672 
00673 
00674 
00675   
00676   
00677   
00678