STEER/AliSplineFit.cxx

   1 /**************************************************************************
   2  * Copyright(c) 2006-07, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16
  17 #include "AliSplineFit.h"
  18
  19 ClassImp(AliSplineFit);
  20
  21 TLinearFitter*
  22 AliSplineFit::fitterStatic()
  23 {
  24   static TLinearFitter* fit = new TLinearFitter(4,"pol3","");
  25   return fit;
  26 }
  27
  28 AliSplineFit::AliSplineFit() :
  29   fBDump(kFALSE),
  30   fGraph (0),
  31   fNmin (0),
  32   fSigma (0),
  33   fMaxDelta (0),
  34   fN0 (0),
  35   fParams (0),
  36   fCovars (0),
  37   fIndex (0),
  38   fN    (0),
  39   fChi2  (0.0),
  40   fX   (0),
  41   fY0  (0),
  42   fY1  (0),
  43   fChi2I (0)
  44   //
  45   // Default constructor
  46   //
  47 { }
  48
  49
  50
  51 AliSplineFit::AliSplineFit(const AliSplineFit& source) :
  52   TObject(source),
  53   fBDump (source.fBDump),
  54   fGraph (source.fGraph),
  55   fNmin  (source.fNmin),
  56   fSigma (source.fSigma),
  57   fMaxDelta (source.fMaxDelta),
  58   fN0    (source.fN0),
  59   fN     (source.fN),
  60   fChi2  (source.fChi2)
  61 {
  62 //
  63 //  Copy constructor
  64 //
  65   fIndex = new Int_t[fN0];
  66   fParams = new TClonesArray("TVectorD",fN0);
  67   fCovars = new TClonesArray("TMatrixD",fN0);
  68   fParams = (TClonesArray*)source.fParams->Clone();
  69   fCovars = (TClonesArray*)source.fCovars->Clone();
  70   for (Int_t i=0; i<fN0; i++) fIndex[i] = source.fIndex[i];
  71
  72   fX     = new Double_t[fN];
  73   fY0    = new Double_t[fN];
  74   fY1    = new Double_t[fN];
  75   fChi2I = new Double_t[fN];
  76   for (Int_t i=0; i<fN; i++){
  77     fX[i]  = source.fX[i];
  78     fY0[i] = source.fY0[i];
  79     fY1[i] = source.fY1[i];
  80   }
  81 }
  82 AliSplineFit& AliSplineFit::operator=(const AliSplineFit& source){
  83 //
  84 //  assignment operator
  85 //
  86   if (&source == this) return *this;
  87
  88 //
  89 // reassign memory as previous fit could have a different size
  90 //
  91
  92   if ( fN0 != source.fN0) {
  93
  94     delete fParams;
  95     delete fCovars;
  96     delete []fIndex;
  97
  98     fN0 = source.fN0;
  99     fIndex = new Int_t[fN0];
 100     fParams = new TClonesArray("TVectorD",fN0);
 101     fCovars = new TClonesArray("TMatrixD",fN0);
 102   }
 103   if ( fN != source.fN) {
 104
 105     delete []fX;
 106     delete []fY0;
 107     delete []fY1;
 108     delete []fChi2I;
 109     fN = source.fN;
 110     fX     = new Double_t[fN];
 111     fY0    = new Double_t[fN];
 112     fY1    = new Double_t[fN];
 113     fChi2I = new Double_t[fN];
 114   }
 115
 116 // use copy constructor (without reassigning memory) to copy values
 117
 118   new (this) AliSplineFit(source);
 119
 120   return *this;
 121 }
 122
 123
 124 AliSplineFit::~AliSplineFit(){
 125   //
 126   // destructor. Don't delete fGraph, as this normally comes as input parameter
 127   //
 128   delete []fX;
 129   delete []fY0;
 130   delete []fY1;
 131   delete []fChi2I;
 132   delete fParams;
 133   delete fCovars;
 134   delete []fIndex;
 135 }
 136
 137 Double_t   AliSplineFit::Eval(Double_t x, Int_t deriv) const{
 138   //
 139   // evaluate value at x
 140   //   deriv = 0: function value
 141   //         = 1: first derivative
 142   //         = 2: 2nd derivative
 143   //         = 3: 3rd derivative
 144   //
 145   //  a2 = -(3*a0 -3*b0 + 2*a1*dx +b1*dx)/(dx*dx)
 146   //  a3 = -(-2*a0+2*b0 - a1*dx - b1*dx)/(dx*dx*dx)
 147
 148   Int_t index = TMath::BinarySearch(fN,fX,x);
 149   if (index<0) index =0;
 150   if (index>fN-2) index =fN-2;
 151   //
 152   Double_t dx   = x-fX[index];
 153   Double_t dxc  = fX[index+1]-fX[index];
 154   Double_t y0   = fY0[index];
 155   Double_t y1   = fY1[index];
 156   Double_t y01  = fY0[index+1];
 157   Double_t y11  = fY1[index+1];
 158   Double_t y2   = -(3.*y0-3.*y01+2*y1*dxc+y11*dxc)/(dxc*dxc);
 159   Double_t y3   = -(-2.* y0 + 2*y01 -  y1*dxc - y11*dxc) /(dxc*dxc*dxc);
 160   Double_t val = y0+y1*dx+y2*dx*dx+y3*dx*dx*dx;
 161   if (deriv==1) val = y1+2.*y2*dx+3.*y3*dx*dx;
 162   if (deriv==2) val = 2.*y2+6.*y3*dx;
 163   if (deriv==3) val = 6*y3;
 164   return val;
 165 }
 166
 167
 168 TGraph * AliSplineFit::GenerGraph(Int_t npoints, Double_t fraction, Double_t s1, Double_t s2, Double_t s3, Int_t der){
 169   //
 170   // generate random graph
 171   // xrange 0,1
 172   // yrange 0,1
 173   // s1, s2, s3 -  sigma of derivative
 174   // fraction   -
 175
 176   Double_t *value = new Double_t[npoints];
 177   Double_t *time  = new Double_t[npoints];
 178   Double_t d0=0, d1=0,d2=0,d3=0;
 179   value[0] = d0;
 180   time[0]  = 0;
 181   for(Int_t i=1; i<npoints; i++){
 182     Double_t dtime = 1./npoints;
 183     Double_t dd1 = dtime;
 184     Double_t dd2 = dd1*dd1;
 185     Double_t dd3 = dd2*dd1;
 186     d0 += d1*dd1 + d2*dd2/2. + d3*dd3/6.;
 187     d1 += d2*dd1 +d3*dd2/2;
 188     d2 += d3*dd1;
 189     value[i] = d0;
 190     time[i]  = time[i-1]+dtime;
 191     d1 =(1.-fraction)*d1+fraction*(gRandom->Exp(s1))*(gRandom->Rndm()-0.5);
 192     d2 =(1.-fraction)*d2+fraction*(gRandom->Exp(s2))*(gRandom->Rndm()-0.5);
 193     d3 =(1.-fraction)*d3+fraction*(gRandom->Exp(s3))*(gRandom->Rndm()-0.5);
 194     if (gRandom->Rndm()<fraction) d3 =(1.-fraction)*d3+fraction*(gRandom->BreitWigner(0,s3));
 195   }
 196   Double_t dmean = (value[npoints-1]-value[0])/(time[npoints-1]-time[0]);
 197   Double_t min = value[0];
 198   Double_t max = value[0];
 199   for (Int_t i=0; i<npoints; i++){
 200     value[i]  = value[i]-dmean*(time[i]-time[0]);
 201     if (value[i]<min) min=value[i];
 202     if (value[i]>max) max=value[i];
 203   }
 204
 205   for (Int_t i=0; i<npoints; i++){
 206     value[i]  = (value[i]-min)/(max-min);
 207   }
 208   if (der==1) for (Int_t i=1; i<npoints; i++){
 209     value[i-1]  =  (value[i]-value[i-1])/(time[i]-time[i-1]);
 210   }
 211
 212   TGraph * graph = new TGraph(npoints,time,value);
 213
 214   delete [] value;
 215   delete [] time;
 216   return graph;
 217 }
 218
 219
 220 TGraph * AliSplineFit::GenerNoise(TGraph * graph0, Double_t sigma0){
 221   //
 222   // add noise to graph
 223   //
 224
 225   Int_t npoints=graph0->GetN();
 226   Double_t *value = new Double_t[npoints];
 227   Double_t *time  = new Double_t[npoints];
 228   for(Int_t i=0; i<npoints; i++){
 229     time[i]  = graph0->GetX()[i];
 230     value[i] = graph0->GetY()[i]+gRandom->Gaus(0,sigma0);
 231   }
 232   TGraph * graph = new TGraph(npoints,time,value);
 233
 234   delete [] value;
 235   delete [] time;
 236   return graph;
 237 }
 238
 239
 240 TGraph * AliSplineFit::MakeGraph(Double_t xmin, Double_t xmax, Int_t npoints, Int_t deriv) const {
 241   //
 242   // if npoints<=0 draw derivative
 243   //
 244
 245   TGraph *graph =0;
 246   if (npoints<=0) {
 247     if (deriv<=0) return new TGraph(fN,fX,fY0);
 248     if (deriv==1) return new TGraph(fN,fX,fY1);
 249     if (deriv>2) return new TGraph(fN-1,fX,fChi2I);
 250   }
 251   Double_t * x = new Double_t[npoints+1];
 252   Double_t * y = new Double_t[npoints+1];
 253   for (Int_t ip=0; ip<=npoints; ip++){
 254     x[ip] = xmin+ (xmax-xmin)*(Double_t(ip)/Double_t(npoints));
 255     y[ip] = Eval(x[ip],deriv);
 256   }
 257
 258   graph = new TGraph(npoints,x,y);
 259   delete [] x;
 260   delete [] y;
 261   return graph;
 262 }
 263
 264 TGraph * AliSplineFit::MakeDiff(TGraph * graph0) const {
 265   //
 266   // Make graph of difference to reference graph
 267   //
 268
 269   Int_t npoints=graph0->GetN();
 270   TGraph *graph =0;
 271   Double_t * x = new Double_t[npoints];
 272   Double_t * y = new Double_t[npoints];
 273   for (Int_t ip=0; ip<npoints; ip++){
 274     x[ip] = graph0->GetX()[ip];
 275     y[ip] = Eval(x[ip],0)-graph0->GetY()[ip];
 276   }
 277   graph = new TGraph(npoints,x,y);
 278   delete [] x;
 279   delete [] y;
 280   return graph;
 281 }
 282
 283
 284 TH1F * AliSplineFit::MakeDiffHisto(TGraph * graph0) const {
 285   //
 286   // Make histogram of difference to reference graph
 287   //
 288
 289   Int_t npoints=graph0->GetN();
 290   Float_t min=1e+39,max=-1e+39;
 291   for (Int_t ip=0; ip<npoints; ip++){
 292     Double_t x = graph0->GetX()[ip];
 293     Double_t y = Eval(x,0)-graph0->GetY()[ip];
 294     if (ip==0) {
 295       min = y;
 296       max = y;
 297     }else{
 298       if (y<min) min=y;
 299       if (y>max) max=y;
 300     }
 301   }
 302
 303   TH1F *his = new TH1F("hdiff","hdiff", 100, min, max);
 304   for (Int_t ip=0; ip<npoints; ip++){
 305     Double_t x = graph0->GetX()[ip];
 306     Double_t y = Eval(x,0)-graph0->GetY()[ip];
 307     his->Fill(y);
 308   }
 309
 310   return his;
 311 }
 312
 313
 314
 315 void AliSplineFit::InitKnots(TGraph * graph, Int_t min, Int_t iter, Double_t maxDelta){
 316   //
 317   // initialize knots + estimate sigma of noise + make initial parameters
 318   //
 319   //
 320
 321   const Double_t kEpsilon = 1.e-7;
 322   fGraph  = graph;
 323   fNmin   = min;
 324   fMaxDelta = maxDelta;
 325   Int_t npoints = fGraph->GetN();
 326   fN0           = (npoints/fNmin)+1;
 327   Float_t delta = Double_t(npoints)/Double_t(fN0-1);
 328
 329   fParams = new TClonesArray("TVectorD",fN0);
 330   fCovars = new TClonesArray("TMatrixD",fN0);
 331   fIndex  = new Int_t[fN0];
 332   TLinearFitter fitterLocal(4,"pol3");  // local fitter
 333   Double_t sigma2 =0;
 334
 335
 336   Double_t yMin=graph->GetY()[0];
 337   Double_t yMax=graph->GetY()[0];
 338
 339   for (Int_t iKnot=0; iKnot<fN0; iKnot++){
 340     Int_t index0 = TMath::Nint(Double_t(iKnot)*Double_t(delta));
 341     Int_t index1 = TMath::Min(TMath::Nint(Double_t(iKnot+1)*Double_t(delta)),npoints-1);
 342     Int_t indexM = (iKnot>0) ? fIndex[iKnot-1]:index0;
 343     fIndex[iKnot]=TMath::Min(index0, npoints-1);
 344     Float_t startX =graph->GetX()[fIndex[iKnot]];
 345
 346     for (Int_t ipoint=indexM; ipoint<index1; ipoint++){
 347       Double_t dxl   =graph->GetX()[ipoint]-startX;
 348       Double_t  y    = graph->GetY()[ipoint];
 349       if (y<yMin) yMin=y;
 350       if (y>yMax) yMax=y;
 351       fitterLocal.AddPoint(&dxl,y,1);
 352     }
 353
 354     fitterLocal.Eval();
 355     sigma2 += fitterLocal.GetChisquare()/Double_t((index1-indexM)-4.);
 356     TMatrixD   * covar = new ((*fCovars)[iKnot]) TMatrixD(4,4);
 357     TVectorD   * param = new ((*fParams)[iKnot]) TVectorD(4);
 358     fitterLocal.GetParameters(*param);
 359     fitterLocal.GetCovarianceMatrix(*covar);
 360     fitterLocal.ClearPoints();
 361   }
 362   fSigma  =TMath::Sqrt(sigma2/Double_t(fN0));   // mean sigma
 363   Double_t tDiff = ((yMax-yMin)+TMath::Abs(yMax)+TMath::Abs(yMin))*kEpsilon;
 364   fSigma += tDiff+fMaxDelta/TMath::Sqrt(npoints);
 365   fMaxDelta +=tDiff;
 366   for (Int_t iKnot=0; iKnot<fN0; iKnot++){
 367     TMatrixD & cov = *((TMatrixD*)fCovars->At(iKnot));
 368     cov*=fSigma*fSigma;
 369   }
 370   OptimizeKnots(iter);
 371
 372   fN = 0;
 373   for (Int_t iKnot=0; iKnot<fN0; iKnot++) if (fIndex[iKnot]>=0) fN++;
 374   fX  = new Double_t[fN];
 375   fY0 = new Double_t[fN];
 376   fY1 = new Double_t[fN];
 377   fChi2I = new Double_t[fN];
 378   Int_t iKnot=0;
 379   for (Int_t i=0; i<fN0; i++){
 380     if (fIndex[i]<0) continue;
 381     if (iKnot>=fN) {
 382       printf("AliSplineFit::InitKnots: Knot number > Max knot number\n");
 383       break;
 384     }
 385     TVectorD   * param = (TVectorD*) fParams->At(i);
 386     fX[iKnot]  = fGraph->GetX()[fIndex[i]];
 387     fY0[iKnot] = (*param)(0);
 388     fY1[iKnot] = (*param)(1);
 389     fChi2I[iKnot] = 0;
 390     iKnot++;
 391   }
 392 }
 393
 394
 395 Int_t AliSplineFit::OptimizeKnots(Int_t nIter){
 396   //
 397   //
 398   //
 399   const Double_t kMaxChi2= 5;
 400   Int_t nKnots=0;
 401   TTreeSRedirector cstream("SplineIter.root");
 402   for (Int_t iIter=0; iIter<nIter; iIter++){
 403     if (fBDump) cstream<<"Fit"<<
 404       "iIter="<<iIter<<
 405       "fit.="<<this<<
 406       "\n";
 407     nKnots=2;
 408     for (Int_t iKnot=1; iKnot<fN0-1; iKnot++){
 409       if (fIndex[iKnot]<0) continue;   //disabled knot
 410       Double_t chi2 = CheckKnot(iKnot);
 411       Double_t startX = fGraph->GetX()[fIndex[iKnot]];
 412       if (fBDump) {
 413         TMatrixD   * covar = (TMatrixD*)fCovars->At(iKnot);
 414         TVectorD   * param = (TVectorD*)fParams->At(iKnot);
 415         cstream<<"Chi2"<<
 416          "iIter="<<iIter<<
 417          "iKnot="<<iKnot<<
 418          "chi2="<<chi2<<
 419          "x="<<startX<<
 420          "param="<<param<<
 421          "covar="<<covar<<
 422          "\n";
 423       }
 424       if (chi2>kMaxChi2) { nKnots++;continue;}
 425       fIndex[iKnot]*=-1;
 426       Int_t iPrevious=iKnot-1;
 427       Int_t iNext    =iKnot+1;
 428       while (fIndex[iPrevious]<0) iPrevious--;
 429       while (fIndex[iNext]<0) iNext++;
 430       RefitKnot(iPrevious);
 431       RefitKnot(iNext);
 432       iKnot++;
 433       while (iKnot<fN0-1&& fIndex[iKnot]<0) iKnot++;
 434     }
 435   }
 436   return nKnots;
 437 }
 438
 439
 440 Bool_t   AliSplineFit::RefitKnot(Int_t iKnot){
 441   //
 442   //
 443   //
 444   const Double_t kEpsilon = 1.e-7;
 445
 446   Int_t iPrevious=(iKnot>0)  ?iKnot-1: 0;
 447   Int_t iNext    =(iKnot<fN0)?iKnot+1: fN0-1;
 448   while (iPrevious>0&&fIndex[iPrevious]<0) iPrevious--;
 449   while (iNext<fN0&&fIndex[iNext]<0) iNext++;
 450   if (iPrevious<0) iPrevious=0;
 451   if (iNext>=fN0) iNext=fN0-1;
 452
 453   Double_t startX = fGraph->GetX()[fIndex[iKnot]];
 454   AliSplineFit::fitterStatic()->ClearPoints();
 455   Int_t indPrev = fIndex[iPrevious];
 456   Int_t indNext = fIndex[iNext];
 457   Double_t *graphX = fGraph->GetX();
 458   Double_t *graphY = fGraph->GetY();
 459
 460   // make arrays for points to fit (to save time)
 461
 462   Int_t nPoints = indNext-indPrev;
 463   Double_t *xPoint = new Double_t[3*nPoints];
 464   Double_t *yPoint = &xPoint[nPoints];
 465   Double_t *ePoint = &xPoint[2*nPoints];
 466   Int_t indVec=0;
 467   for (Int_t iPoint=indPrev; iPoint<indNext; iPoint++, indVec++){
 468     Double_t dxl   = graphX[iPoint]-startX;
 469     Double_t  y    = graphY[iPoint];
 470     xPoint[indVec] = dxl;
 471     yPoint[indVec] = y;
 472     ePoint[indVec] =  fSigma;
 473 //    ePoint[indVec] =  fSigma+TMath::Abs(y)*kEpsilon;
 474 //    AliSplineFit::fitterStatic.AddPoint(&dxl,y,fSigma+TMath::Abs(y)*kEpsilon);
 475   }
 476   AliSplineFit::fitterStatic()->AssignData(nPoints,1,xPoint,yPoint,ePoint);
 477   AliSplineFit::fitterStatic()->Eval();
 478
 479 //  delete temporary arrays
 480
 481   delete [] xPoint;
 482
 483   TMatrixD   * covar = (TMatrixD*)fCovars->At(iKnot);
 484   TVectorD   * param = (TVectorD*)fParams->At(iKnot);
 485   AliSplineFit::fitterStatic()->GetParameters(*param);
 486   AliSplineFit::fitterStatic()->GetCovarianceMatrix(*covar);
 487   return 0;
 488 }
 489
 490
 491 Float_t AliSplineFit::CheckKnot(Int_t iKnot){
 492   //
 493   //
 494   //
 495
 496   Int_t iPrevious=iKnot-1;
 497   Int_t iNext    =iKnot+1;
 498   while (fIndex[iPrevious]<0) iPrevious--;
 499   while (fIndex[iNext]<0) iNext++;
 500   TVectorD &pPrevious = *((TVectorD*)fParams->At(iPrevious));
 501   TVectorD &pNext     = *((TVectorD*)fParams->At(iNext));
 502   TVectorD &pKnot     = *((TVectorD*)fParams->At(iKnot));
 503   TMatrixD &cPrevious = *((TMatrixD*)fCovars->At(iPrevious));
 504   TMatrixD &cNext     = *((TMatrixD*)fCovars->At(iNext));
 505   TMatrixD &cKnot     = *((TMatrixD*)fCovars->At(iKnot));
 506   Double_t xPrevious = fGraph->GetX()[fIndex[iPrevious]];
 507   Double_t xNext     = fGraph->GetX()[fIndex[iNext]];
 508   Double_t xKnot     = fGraph->GetX()[fIndex[iKnot]];
 509
 510   // extra variables introduced to save processing time
 511
 512   Double_t dxc  = xNext-xPrevious;
 513   Double_t invDxc = 1./dxc;
 514   Double_t invDxc2 = invDxc*invDxc;
 515   TMatrixD  tPrevious(4,4);
 516   TMatrixD  tNext(4,4);
 517
 518   tPrevious(0,0) = 1;    tPrevious(1,1) = 1;
 519   tPrevious(2,0) = -3.*invDxc2;
 520   tPrevious(2,1) = -2.*invDxc;
 521   tPrevious(3,0) =  2.*invDxc2*invDxc;
 522   tPrevious(3,1) =  1.*invDxc2;
 523   tNext(2,0)     =  3.*invDxc2;      tNext(2,1)     = -1*invDxc;
 524   tNext(3,0)     = -2.*invDxc2*invDxc;  tNext(3,1)     =  1.*invDxc2;
 525   TMatrixD  tpKnot(4,4);
 526   TMatrixD  tpNext(4,4);
 527   Double_t dx = xKnot-xPrevious;
 528   tpKnot(0,0) = 1;      tpKnot(1,1) = 1;        tpKnot(2,2) = 1;        tpKnot(3,3) = 1;
 529   tpKnot(0,1) = dx;     tpKnot(0,2) = dx*dx;    tpKnot(0,3) = dx*dx*dx;
 530   tpKnot(1,2) = 2.*dx;  tpKnot(1,3) = 3.*dx*dx;
 531   tpKnot(2,3) = 3.*dx;
 532   Double_t dxn = xNext-xPrevious;
 533   tpNext(0,0) = 1;       tpNext(1,1) = 1;        tpNext(2,2) = 1;        tpNext(3,3) = 1;
 534   tpNext(0,1) = dxn;     tpNext(0,2) = dxn*dxn;    tpNext(0,3) = dxn*dxn*dxn;
 535   tpNext(1,2) = 2.*dxn;  tpNext(1,3) = 3.*dxn*dxn;
 536   tpNext(2,3) = 3.*dxn;
 537
 538   //
 539   // matrix and vector at previous
 540   //
 541
 542   TVectorD  sPrevious = tPrevious*pPrevious+tNext*pNext;
 543   TVectorD  sKnot     = tpKnot*sPrevious;
 544   TVectorD  sNext     = tpNext*sPrevious;
 545
 546   TMatrixD csPrevious00(tPrevious, TMatrixD::kMult,cPrevious);
 547   csPrevious00 *= tPrevious.T();
 548   TMatrixD csPrevious01(tNext,TMatrixD::kMult,cNext);
 549   csPrevious01*=tNext.T();
 550   TMatrixD  csPrevious(csPrevious00,TMatrixD::kPlus,csPrevious01);
 551   TMatrixD  csKnot(tpKnot,TMatrixD::kMult,csPrevious);
 552   csKnot*=tpKnot.T();
 553   TMatrixD  csNext(tpNext,TMatrixD::kMult,csPrevious);
 554   csNext*=tpNext.T();
 555
 556   TVectorD dPrevious = pPrevious-sPrevious;
 557   TVectorD dKnot     = pKnot-sKnot;
 558   TVectorD dNext     = pNext-sNext;
 559   //
 560   //
 561   TMatrixD prec(4,4);
 562   prec(0,0) = (fMaxDelta*fMaxDelta);
 563   prec(1,1) = prec(0,0)*invDxc2;
 564   prec(2,2) = prec(1,1)*invDxc2;
 565   prec(3,3) = prec(2,2)*invDxc2;
 566
 567 //   prec(0,0) = (fMaxDelta*fMaxDelta);
 568 //   prec(1,1) = (fMaxDelta*fMaxDelta)/(dxc*dxc);
 569 //   prec(2,2) = (fMaxDelta*fMaxDelta)/(dxc*dxc*dxc*dxc);
 570 //   prec(3,3) = (fMaxDelta*fMaxDelta)/(dxc*dxc*dxc*dxc*dxc*dxc);
 571
 572   csPrevious+=cPrevious;
 573   csPrevious+=prec;
 574   csPrevious.Invert();
 575   Double_t  chi2P     = dPrevious*(csPrevious*dPrevious);
 576
 577   csKnot+=cKnot;
 578   csKnot+=prec;
 579   csKnot.Invert();
 580   Double_t  chi2K     = dKnot*(csKnot*dKnot);
 581
 582   csNext+=cNext;
 583   csNext+=prec;
 584   csNext.Invert();
 585   Double_t  chi2N     = dNext*(csNext*dNext);
 586
 587   return (chi2P+chi2K+chi2N)/8.;
 588
 589
 590 }
 591
 592 void AliSplineFit::SplineFit(Int_t nder){
 593   //
 594   // Cubic spline fit of graph
 595   //
 596   // nder
 597   // nder<0  - no continuity requirement
 598   //     =0  - continous  0 derivative
 599   //     =1  - continous  1 derivative
 600   //     >1  - continous  2 derivative
 601   //
 602   if (!fGraph) return;
 603   TGraph * graph = fGraph;
 604   if (nder>1) nder=2;
 605   Int_t nknots  = fN;
 606   Int_t npoints = graph->GetN();
 607   //
 608   //
 609   // spline fit
 610   // each knot 4 parameters
 611   //
 612   TMatrixD       *pmatrix = 0;
 613   TVectorD       *pvalues = 0;
 614   if (nder>1){
 615     pmatrix = new TMatrixD(4*(nknots-1)+3*(nknots-2), 4*(nknots-1)+3*(nknots-2));
 616     pvalues = new TVectorD(4*(nknots-1)+3*(nknots-2));
 617   }
 618   if (nder==1){
 619     pmatrix = new TMatrixD(4*(nknots-1)+2*(nknots-2), 4*(nknots-1)+2*(nknots-2));
 620     pvalues = new TVectorD(4*(nknots-1)+2*(nknots-2));
 621   }
 622   if (nder==0){
 623     pmatrix = new TMatrixD(4*(nknots-1)+1*(nknots-2), 4*(nknots-1)+1*(nknots-2));
 624     pvalues = new TVectorD(4*(nknots-1)+1*(nknots-2));
 625   }
 626   if (nder<0){
 627     pmatrix = new TMatrixD(4*(nknots-1)+0*(nknots-2), 4*(nknots-1)+0*(nknots-2));
 628     pvalues = new TVectorD(4*(nknots-1)+0*(nknots-2));
 629   }
 630
 631
 632   TMatrixD &matrix = *pmatrix;
 633   TVectorD &values = *pvalues;
 634   Int_t    current = 0;
 635 //
 636 //  defined extra variables (current4 etc.) to save processing time.
 637 //  fill normal matrices, then copy to sparse matrix.
 638 //
 639   Double_t *graphX = graph->GetX();
 640   Double_t *graphY = graph->GetY();
 641   for (Int_t ip=0;ip<npoints;ip++){
 642     if (current<nknots-2&&graphX[ip]>fX[current+1]) current++;
 643     Double_t xmiddle = (fX[current+1]+fX[current])*0.5;
 644     Double_t x1 = graphX[ip]- xmiddle;
 645     Double_t x2 = x1*x1;
 646     Double_t x3 = x2*x1;
 647     Double_t x4 = x2*x2;
 648     Double_t x5 = x3*x2;
 649     Double_t x6 = x3*x3;
 650     Double_t y  = graphY[ip];
 651     Int_t current4 = 4*current;
 652
 653     matrix(current4  , current4  )+=1;
 654     matrix(current4  , current4+1)+=x1;
 655     matrix(current4  , current4+2)+=x2;
 656     matrix(current4  , current4+3)+=x3;
 657     //
 658     matrix(current4+1, current4  )+=x1;
 659     matrix(current4+1, current4+1)+=x2;
 660     matrix(current4+1, current4+2)+=x3;
 661     matrix(current4+1, current4+3)+=x4;
 662     //
 663     matrix(current4+2, current4  )+=x2;
 664     matrix(current4+2, current4+1)+=x3;
 665     matrix(current4+2, current4+2)+=x4;
 666     matrix(current4+2, current4+3)+=x5;
 667     //
 668     matrix(current4+3, current4  )+=x3;
 669     matrix(current4+3, current4+1)+=x4;
 670     matrix(current4+3, current4+2)+=x5;
 671     matrix(current4+3, current4+3)+=x6;
 672     //
 673     values(current4  ) += y;
 674     values(current4+1) += y*x1;
 675     values(current4+2) += y*x2;
 676     values(current4+3) += y*x3;
 677   }
 678   //
 679   // constraint 0
 680   //
 681   Int_t offset =4*(nknots-1)-1;
 682   if (nder>=0) for (Int_t iknot = 1; iknot<nknots-1; iknot++){
 683
 684     Double_t dxm  =  (fX[iknot]-fX[iknot-1])*0.5;
 685     Double_t dxp  = -(fX[iknot+1]-fX[iknot])*0.5;
 686     Double_t dxm2 = dxm*dxm;
 687     Double_t dxp2 = dxp*dxp;
 688     Double_t dxm3 = dxm2*dxm;
 689     Double_t dxp3 = dxp2*dxp;
 690     Int_t iknot4  = 4*iknot;
 691     Int_t iknot41 = 4*(iknot-1);
 692     Int_t offsKnot = offset+iknot;
 693     //
 694     // condition on knot
 695     //
 696     // a0[i] = a0m[i-1]  + a1m[i-1]*dxm + a2m[i-1]*dxm^2 + a3m[i-1]*dxm^3
 697     // a0[i] = a0m[i-0]  + a1m[i-0]*dxp + a2m[i-0]*dxp^2 + a3m[i-0]*dxp^3
 698     // (a0m[i-1]  + a1m[i-1]*dxm + a2m[i-1]*dxm^2 + a3m[i-1]*dxm^3) -
 699     // (a0m[i-0]  + a1m[i-0]*dxp + a2m[i-0]*dxp^2 + a3m[i-0]*dxp^3)  = 0
 700
 701     matrix(offsKnot, iknot41  )=1;
 702     matrix(offsKnot, iknot4   )=-1;
 703
 704     matrix(offsKnot, iknot41+1)=dxm;
 705     matrix(offsKnot, iknot4 +1)=-dxp;
 706
 707     matrix(offsKnot, iknot41+2)=dxm2;
 708     matrix(offsKnot, iknot4 +2)=-dxp2;
 709
 710     matrix(offsKnot, iknot41+3)=dxm3;
 711     matrix(offsKnot, iknot4 +3)=-dxp3;
 712
 713     matrix(iknot41  , offsKnot)=1;
 714     matrix(iknot41+1, offsKnot)=dxm;
 715     matrix(iknot41+2, offsKnot)=dxm2;
 716     matrix(iknot41+3, offsKnot)=dxm3;
 717     matrix(iknot4  , offsKnot)=-1;
 718     matrix(iknot4+1, offsKnot)=-dxp;
 719     matrix(iknot4+2, offsKnot)=-dxp2;
 720     matrix(iknot4+3, offsKnot)=-dxp3;
 721   }
 722   //
 723   // constraint 1
 724   //
 725   offset =4*(nknots-1)-1+(nknots-2);
 726   if (nder>=1)for (Int_t iknot = 1; iknot<nknots-1; iknot++){
 727
 728     Double_t dxm  =  (fX[iknot]-fX[iknot-1])*0.5;
 729     Double_t dxp  = -(fX[iknot+1]-fX[iknot])*0.5;
 730     Double_t dxm2 = dxm*dxm;
 731     Double_t dxp2 = dxp*dxp;
 732     Int_t iknot4  = 4*iknot;
 733     Int_t iknot41 = 4*(iknot-1);
 734     Int_t offsKnot = offset+iknot;
 735     //
 736     // condition on knot derivation
 737     //
 738     // a0d[i] =  a1m[i-1] + 2*a2m[i-1]*dxm + 3*a3m[i-1]*dxm^2
 739     // a0d[i] =  a1m[i-0] + 2*a2m[i-0]*dxp + 3*a3m[i-0]*dxp^2
 740
 741     //
 742     matrix(offsKnot, iknot41+1)= 1;
 743     matrix(offsKnot, iknot4 +1)=-1;
 744
 745     matrix(offsKnot, iknot41+2)= 2.*dxm;
 746     matrix(offsKnot, iknot4 +2)=-2.*dxp;
 747
 748     matrix(offsKnot, iknot41+3)= 3.*dxm2;
 749     matrix(offsKnot, iknot4 +3)=-3.*dxp2;
 750
 751     matrix(iknot41+1, offsKnot)=1;
 752     matrix(iknot41+2, offsKnot)=2.*dxm;
 753     matrix(iknot41+3, offsKnot)=3.*dxm2;
 754
 755     matrix(iknot4+1, offsKnot)=-1.;
 756     matrix(iknot4+2, offsKnot)=-2.*dxp;
 757     matrix(iknot4+3, offsKnot)=-3.*dxp2;
 758   }
 759   //
 760   // constraint 2
 761   //
 762   offset =4*(nknots-1)-1+2*(nknots-2);
 763   if (nder>=2) for (Int_t iknot = 1; iknot<nknots-1; iknot++){
 764
 765     Double_t dxm  =  (fX[iknot]-fX[iknot-1])*0.5;
 766     Double_t dxp  = -(fX[iknot+1]-fX[iknot])*0.5;
 767     Int_t iknot4  = 4*iknot;
 768     Int_t iknot41 = 4*(iknot-1);
 769     Int_t offsKnot = offset+iknot;
 770     //
 771     // condition on knot second derivative
 772     //
 773     // a0dd[i] =  2*a2m[i-1] + 6*a3m[i-1]*dxm
 774     // a0dd[i] =  2*a2m[i-0] + 6*a3m[i-0]*dxp
 775     //
 776     //
 777     matrix(offsKnot, iknot41+2)= 2.;
 778     matrix(offsKnot, iknot4 +2)=-2.;
 779
 780     matrix(offsKnot, iknot41+3)= 6.*dxm;
 781     matrix(offsKnot, iknot4 +3)=-6.*dxp;
 782
 783     matrix(iknot41+2, offsKnot)=2.;
 784     matrix(iknot41+3, offsKnot)=6.*dxm;
 785
 786     matrix(iknot4+2, offsKnot)=-2.;
 787     matrix(iknot4+3, offsKnot)=-6.*dxp;
 788   }
 789
 790 // sparse matrix to do fit
 791
 792   TMatrixDSparse smatrix(matrix);
 793   TDecompSparse svd(smatrix,0);
 794   Bool_t ok;
 795   const TVectorD results = svd.Solve(values,ok);
 796
 797   for (Int_t iknot = 0; iknot<nknots-1; iknot++){
 798
 799     Double_t dxm  =  -(fX[iknot+1]-fX[iknot])*0.5;
 800
 801     fY0[iknot] = results(4*iknot)+ results(4*iknot+1)*dxm+results(4*iknot+2)*dxm*dxm+
 802       results(4*iknot+3)*dxm*dxm*dxm;
 803
 804     fY1[iknot] = results(4*iknot+1)+2.*results(4*iknot+2)*dxm+
 805       3*results(4*iknot+3)*dxm*dxm;
 806   }
 807   Int_t   iknot2= nknots-1;
 808   Int_t   iknot = nknots-2;
 809   Double_t dxm   =  (fX[iknot2]-fX[iknot2-1])*0.5;
 810
 811   fY0[iknot2] = results(4*iknot)+ results(4*iknot+1)*dxm+results(4*iknot+2)*dxm*dxm+
 812     results(4*iknot+3)*dxm*dxm*dxm;
 813
 814   fY1[iknot2] = results(4*iknot+1)+2.*results(4*iknot+2)*dxm+
 815       3*results(4*iknot+3)*dxm*dxm;
 816
 817   delete  pmatrix;
 818   delete  pvalues;
 819
 820 }
 821
 822
 823
 824
 825
 826 void AliSplineFit::MakeKnots0(TGraph * graph, Double_t maxdelta, Int_t minpoints){
 827   //
 828   // make knots  - restriction max distance and minimum points
 829   //
 830
 831   Int_t npoints  = graph->GetN();
 832   Double_t *xknots = new Double_t[npoints];
 833   Int_t nknots =0;
 834   Int_t ipoints =0;
 835   //
 836   // generate knots
 837   //
 838   for (Int_t ip=0;ip<npoints;ip++){
 839     if (graph->GetX()[ip]-xknots[nknots-1]>maxdelta && ipoints>minpoints){
 840       xknots[nknots] = graph->GetX()[ip];
 841       ipoints=1;
 842       nknots++;
 843     }
 844     ipoints++;
 845   }
 846   if (npoints-ipoints>minpoints){
 847     xknots[nknots] = graph->GetX()[npoints-1];
 848     nknots++;
 849   }else{
 850     xknots[nknots-1] = graph->GetX()[npoints-1];
 851   }
 852
 853   fN = nknots;
 854   fX = new Double_t[nknots];
 855   fY0 = new Double_t[nknots];
 856   fY1 = new Double_t[nknots];
 857   fChi2I= new Double_t[nknots];
 858   for (Int_t i=0; i<nknots; i++) fX[i]= xknots[i];
 859   delete [] xknots;
 860 }
 861
 862
 863
 864
 865 void AliSplineFit::MakeSmooth(TGraph * graph, Float_t ratio, char * type){
 866   //
 867   // Interface to GraphSmooth
 868   //
 869
 870   TGraphSmooth smooth;
 871   Int_t    npoints2 = TMath::Nint(graph->GetN()*ratio);
 872   TGraph * graphT0 = smooth.SmoothKern(graph,type,ratio);
 873   if (!graphT0) return;
 874   TGraph  graphT1(npoints2);
 875   for (Int_t ipoint=0; ipoint<npoints2; ipoint++){
 876     Int_t pointS = TMath::Nint(ipoint/ratio);
 877     if (ipoint==npoints2-1) pointS=graph->GetN()-1;
 878     graphT1.SetPoint(ipoint, graphT0->GetX()[pointS] , graphT0->GetY()[pointS]);
 879   }
 880   TSpline3 spline2("spline", &graphT1);
 881   Update(&spline2, npoints2);
 882 }
 883
 884
 885 void AliSplineFit::Update(TSpline3 *spline, Int_t nknots){
 886   //
 887   //
 888   //
 889
 890   fN = nknots;
 891   fX = new Double_t[nknots];
 892   fY0 = new Double_t[nknots];
 893   fY1 = new Double_t[nknots];
 894   Double_t d0, d1;
 895   fChi2I= 0;
 896   for (Int_t i=0; i<nknots; i++) {
 897     spline->GetCoeff(i,fX[i],fY0[i], fY1[i],d0,d1);
 898   }
 899 }
 900
 901
 902
 903
 904 void AliSplineFit::Test(Int_t npoints, Int_t ntracks, Float_t snoise){
 905   //
 906   // test function
 907   //
 908
 909   AliSplineFit fit;
 910   AliSplineFit fitS;
 911   TGraph * graph0=0;
 912   TGraph * graph1=0;
 913
 914   TTreeSRedirector *pcstream = new TTreeSRedirector("TestSmooth.root");
 915   for (Int_t i=0; i<ntracks; i++){
 916     graph0 = AliSplineFit::GenerGraph(npoints,0.05,0,0,1,0);
 917     graph1 = AliSplineFit::GenerNoise(graph0,snoise);
 918     fit.InitKnots(graph1, 10,10, 0.00);
 919     TGraph *d0 = fit.MakeDiff(graph0);
 920     TGraph *g0 = fit.MakeGraph(0,1,1000,0);
 921     fit.SplineFit(2);
 922     TH1F * h2 = fit.MakeDiffHisto(graph0);
 923     TGraph *d2 = fit.MakeDiff(graph0);
 924     TGraph *g2 = fit.MakeGraph(0,1,1000,0);
 925     fit.SplineFit(1);
 926     TH1F * h1 = fit.MakeDiffHisto(graph0);
 927     TGraph *d1 = fit.MakeDiff(graph0);
 928     TGraph *g1 = fit.MakeGraph(0,1,1000,0);
 929
 930     Float_t ratio = Float_t(fit.fN)/Float_t(npoints);
 931     fitS.MakeSmooth(graph1,ratio,"box");
 932     TGraph *dS = fitS.MakeDiff(graph0);
 933     TGraph *gS = fit.MakeGraph(0,1,1000,0);
 934
 935     TH1F * hS = fitS.MakeDiffHisto(graph0);
 936     Double_t mean2  = h2->GetMean();
 937     Double_t sigma2 = h2->GetRMS();
 938     Double_t mean1  = h1->GetMean();
 939     Double_t sigma1 = h1->GetRMS();
 940     Double_t meanS  = hS->GetMean();
 941     Double_t sigmaS = hS->GetRMS();
 942     char fname[100];
 943     if (fit.fN<20){
 944       sprintf(fname,"pol%d",fit.fN);
 945     }else{
 946       sprintf(fname,"pol%d",19);
 947     }
 948     TF1 fpol("fpol",fname);
 949     graph1->Fit(&fpol);
 950     TGraph dpol(*graph1);
 951     TGraph gpol(*graph1);
 952     for (Int_t ipoint=0; ipoint<graph1->GetN(); ipoint++){
 953       dpol.GetY()[ipoint]= graph0->GetY()[ipoint]-
 954         fpol.Eval(graph0->GetX()[ipoint]);
 955       gpol.GetY()[ipoint]= fpol.Eval(graph0->GetX()[ipoint]);
 956     }
 957     (*pcstream)<<"Test"<<
 958       "Event="<<i<<
 959       "Graph0.="<<graph0<<
 960       "Graph1.="<<graph1<<
 961       "G0.="<<g0<<
 962       "G1.="<<g1<<
 963       "G2.="<<g2<<
 964       "GS.="<<gS<<
 965       "GP.="<<&gpol<<
 966       "D0.="<<d0<<
 967       "D1.="<<d1<<
 968       "D2.="<<d2<<
 969       "DS.="<<dS<<
 970       "DP.="<<&dpol<<
 971       "Npoints="<<fit.fN<<
 972       "Mean1="<<mean1<<
 973       "Mean2="<<mean2<<
 974       "MeanS="<<meanS<<
 975       "Sigma1="<<sigma1<<
 976       "Sigma2="<<sigma2<<
 977       "SigmaS="<<sigmaS<<
 978       "\n";
 979
 980     delete graph0;
 981     delete graph1;
 982     delete g1;
 983     delete g2;
 984     delete gS;
 985     delete h1;
 986     delete h2;
 987     delete hS;
 988   }
 989   delete pcstream;
 990 }