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