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