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 0;
 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=1e+39,max=-1e+39;
 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   Int_t nknots =0;
 864   Int_t ipoints =0;
 865   //
 866   // generate knots
 867   //
 868   for (Int_t ip=0;ip<npoints;ip++){
 869     if (graph->GetX()[ip]-xknots[nknots-1]>maxdelta && ipoints>minpoints){
 870       xknots[nknots] = graph->GetX()[ip];
 871       ipoints=1;
 872       nknots++;
 873     }
 874     ipoints++;
 875   }
 876   if (npoints-ipoints>minpoints){
 877     xknots[nknots] = graph->GetX()[npoints-1];
 878     nknots++;
 879   }else{
 880     xknots[nknots-1] = graph->GetX()[npoints-1];
 881   }
 882
 883   fN = nknots;
 884   fX = new Double_t[nknots];
 885   fY0 = new Double_t[nknots];
 886   fY1 = new Double_t[nknots];
 887   fChi2I= new Double_t[nknots];
 888   for (Int_t i=0; i<nknots; i++) fX[i]= xknots[i];
 889   delete [] xknots;
 890 }
 891
 892
 893
 894
 895 void AliSplineFit::MakeSmooth(TGraph * graph, Float_t ratio, Option_t * type){
 896   //
 897   // Interface to GraphSmooth
 898   //
 899
 900   TGraphSmooth smooth;
 901   Int_t    npoints2 = TMath::Nint(graph->GetN()*ratio);
 902   TGraph * graphT0 = smooth.SmoothKern(graph,type,ratio);
 903   if (!graphT0) return;
 904   TGraph  graphT1(npoints2);
 905   for (Int_t ipoint=0; ipoint<npoints2; ipoint++){
 906     Int_t pointS = TMath::Nint(ipoint/ratio);
 907     if (ipoint==npoints2-1) pointS=graph->GetN()-1;
 908     graphT1.SetPoint(ipoint, graphT0->GetX()[pointS] , graphT0->GetY()[pointS]);
 909   }
 910   TSpline3 spline2("spline", &graphT1);
 911   Update(&spline2, npoints2);
 912 }
 913
 914
 915 void AliSplineFit::Update(TSpline3 *spline, Int_t nknots){
 916   //
 917   //
 918   //
 919
 920   fN = nknots;
 921   fX = new Double_t[nknots];
 922   fY0 = new Double_t[nknots];
 923   fY1 = new Double_t[nknots];
 924   Double_t d0, d1;
 925   fChi2I= 0;
 926   for (Int_t i=0; i<nknots; i++) {
 927     spline->GetCoeff(i,fX[i],fY0[i], fY1[i],d0,d1);
 928   }
 929 }
 930
 931
 932
 933
 934 void AliSplineFit::Test(Int_t npoints, Int_t ntracks, Float_t snoise){
 935   //
 936   // test function
 937   //
 938
 939   AliSplineFit fit;
 940   AliSplineFit fitS;
 941   TGraph * graph0=0;
 942   TGraph * graph1=0;
 943
 944   TTreeSRedirector *pcstream = new TTreeSRedirector("TestSmooth.root");
 945   for (Int_t i=0; i<ntracks; i++){
 946     graph0 = AliSplineFit::GenerGraph(npoints,0.05,0,0,1,0);
 947     graph1 = AliSplineFit::GenerNoise(graph0,snoise);
 948     fit.InitKnots(graph1, 10,10, 0.00);
 949     TGraph *d0 = fit.MakeDiff(graph0);
 950     TGraph *g0 = fit.MakeGraph(0,1,1000,0);
 951     fit.SplineFit(2);
 952     TH1F * h2 = fit.MakeDiffHisto(graph0);
 953     TGraph *d2 = fit.MakeDiff(graph0);
 954     TGraph *g2 = fit.MakeGraph(0,1,1000,0);
 955     fit.SplineFit(1);
 956     TH1F * h1 = fit.MakeDiffHisto(graph0);
 957     TGraph *d1 = fit.MakeDiff(graph0);
 958     TGraph *g1 = fit.MakeGraph(0,1,1000,0);
 959
 960     Float_t ratio = Float_t(fit.fN)/Float_t(npoints);
 961     fitS.MakeSmooth(graph1,ratio,"box");
 962     TGraph *dS = fitS.MakeDiff(graph0);
 963     TGraph *gS = fit.MakeGraph(0,1,1000,0);
 964
 965     TH1F * hS = fitS.MakeDiffHisto(graph0);
 966     Double_t mean2  = h2->GetMean();
 967     Double_t sigma2 = h2->GetRMS();
 968     Double_t mean1  = h1->GetMean();
 969     Double_t sigma1 = h1->GetRMS();
 970     Double_t meanS  = hS->GetMean();
 971     Double_t sigmaS = hS->GetRMS();
 972     char fname[100];
 973     if (fit.fN<20){
 974       sprintf(fname,"pol%d",fit.fN);
 975     }else{
 976       sprintf(fname,"pol%d",19);
 977     }
 978     TF1 fpol("fpol",fname);
 979     graph1->Fit(&fpol);
 980     TGraph dpol(*graph1);
 981     TGraph gpol(*graph1);
 982     for (Int_t ipoint=0; ipoint<graph1->GetN(); ipoint++){
 983       dpol.GetY()[ipoint]= graph0->GetY()[ipoint]-
 984         fpol.Eval(graph0->GetX()[ipoint]);
 985       gpol.GetY()[ipoint]= fpol.Eval(graph0->GetX()[ipoint]);
 986     }
 987     (*pcstream)<<"Test"<<
 988       "Event="<<i<<
 989       "Graph0.="<<graph0<<
 990       "Graph1.="<<graph1<<
 991       "G0.="<<g0<<
 992       "G1.="<<g1<<
 993       "G2.="<<g2<<
 994       "GS.="<<gS<<
 995       "GP.="<<&gpol<<
 996       "D0.="<<d0<<
 997       "D1.="<<d1<<
 998       "D2.="<<d2<<
 999       "DS.="<<dS<<
1000       "DP.="<<&dpol<<
1001       "Npoints="<<fit.fN<<
1002       "Mean1="<<mean1<<
1003       "Mean2="<<mean2<<
1004       "MeanS="<<meanS<<
1005       "Sigma1="<<sigma1<<
1006       "Sigma2="<<sigma2<<
1007       "SigmaS="<<sigmaS<<
1008       "\n";
1009
1010     delete graph0;
1011     delete graph1;
1012     delete g1;
1013     delete g2;
1014     delete gS;
1015     delete h1;
1016     delete h2;
1017     delete hS;
1018   }
1019   delete pcstream;
1020 }
1021
1022 void AliSplineFit::Cleanup(){
1023  //
1024  // deletes extra information to reduce amount of information stored on the data
1025  // base
1026
1027  //    delete fGraph;  fGraph=0;   // Don't delete fGraph -- input parameter
1028      delete fParams; fParams=0;
1029      delete fCovars; fCovars=0;
1030      delete [] fIndex; fIndex=0;
1031      delete [] fChi2I; fChi2I=0;
1032 }
1033
1034
1035 void AliSplineFit::CopyGraph() {
1036   //
1037   // enter graph points directly to fit parameters
1038   // (to bee used when too few points are available)
1039   //
1040   fN = fGraph->GetN();
1041   fX = new Double_t[fN];
1042   fY0 = new Double_t[fN];
1043   fY1 = new Double_t[fN];
1044   for (Int_t i=0; i<fN; i++ ) {
1045     fX[i]  = fGraph->GetX()[i];
1046     fY0[i] = fGraph->GetY()[i];
1047     fY1[i] = 0;
1048   }
1049 }