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