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