STEER/CDB/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)
 269       return new TGraph(fN,fX,fY0);
 270     else
 271       if (deriv==1)
 272         return new TGraph(fN,fX,fY1);
 273       else
 274         return new TGraph(fN-1,fX,fChi2I);
 275   }
 276   Double_t * x = new Double_t[npoints+1];
 277   Double_t * y = new Double_t[npoints+1];
 278   for (Int_t ip=0; ip<=npoints; ip++){
 279     x[ip] = xmin+ (xmax-xmin)*(Double_t(ip)/Double_t(npoints));
 280     y[ip] = Eval(x[ip],deriv);
 281   }
 282
 283   graph = new TGraph(npoints,x,y);
 284   delete [] x;
 285   delete [] y;
 286   return graph;
 287 }
 288
 289 TGraph * AliSplineFit::MakeDiff(TGraph * graph0) const {
 290   //
 291   // Make graph of difference to reference graph
 292   //
 293
 294   Int_t npoints=graph0->GetN();
 295   TGraph *graph =0;
 296   Double_t * x = new Double_t[npoints];
 297   Double_t * y = new Double_t[npoints];
 298   for (Int_t ip=0; ip<npoints; ip++){
 299     x[ip] = graph0->GetX()[ip];
 300     y[ip] = Eval(x[ip],0)-graph0->GetY()[ip];
 301   }
 302   graph = new TGraph(npoints,x,y);
 303   delete [] x;
 304   delete [] y;
 305   return graph;
 306 }
 307
 308
 309 TH1F * AliSplineFit::MakeDiffHisto(TGraph * graph0) const {
 310   //
 311   // Make histogram of difference to reference graph
 312   //
 313
 314   Int_t npoints=graph0->GetN();
 315   Float_t min=1e38,max=-1e38;
 316   for (Int_t ip=0; ip<npoints; ip++){
 317     Double_t x = graph0->GetX()[ip];
 318     Double_t y = Eval(x,0)-graph0->GetY()[ip];
 319     if (ip==0) {
 320       min = y;
 321       max = y;
 322     }else{
 323       if (y<min) min=y;
 324       if (y>max) max=y;
 325     }
 326   }
 327
 328   TH1F *his = new TH1F("hdiff","hdiff", 100, min, max);
 329   for (Int_t ip=0; ip<npoints; ip++){
 330     Double_t x = graph0->GetX()[ip];
 331     Double_t y = Eval(x,0)-graph0->GetY()[ip];
 332     his->Fill(y);
 333   }
 334
 335   return his;
 336 }
 337
 338
 339
 340 void AliSplineFit::InitKnots(TGraph * graph, Int_t min, Int_t iter, Double_t maxDelta){
 341   //
 342   // initialize knots + estimate sigma of noise + make initial parameters
 343   //
 344   //
 345
 346   const Double_t kEpsilon = 1.e-7;
 347   fGraph  = graph;
 348   fNmin   = min;
 349   fMaxDelta = maxDelta;
 350   Int_t npoints = fGraph->GetN();
 351
 352   // Make simple spline if too few points in graph
 353
 354   if (npoints < fMinPoints ) {
 355     CopyGraph();
 356     return;
 357   }
 358
 359   fN0           = (npoints/fNmin)+1;
 360   Float_t delta = Double_t(npoints)/Double_t(fN0-1);
 361
 362   fParams = new TClonesArray("TVectorD",fN0);
 363   fCovars = new TClonesArray("TMatrixD",fN0);
 364   fIndex  = new Int_t[fN0];
 365   TLinearFitter fitterLocal(4,"pol3");  // local fitter
 366   Double_t sigma2 =0;
 367
 368
 369   Double_t yMin=graph->GetY()[0];
 370   Double_t yMax=graph->GetY()[0];
 371
 372   for (Int_t iKnot=0; iKnot<fN0; iKnot++){
 373     Int_t index0 = TMath::Nint(Double_t(iKnot)*Double_t(delta));
 374     Int_t index1 = TMath::Min(TMath::Nint(Double_t(iKnot+1)*Double_t(delta)),npoints-1);
 375     Int_t indexM = (iKnot>0) ? fIndex[iKnot-1]:index0;
 376     fIndex[iKnot]=TMath::Min(index0, npoints-1);
 377     Float_t startX =graph->GetX()[fIndex[iKnot]];
 378
 379     for (Int_t ipoint=indexM; ipoint<index1; ipoint++){
 380       Double_t dxl   =graph->GetX()[ipoint]-startX;
 381       Double_t  y    = graph->GetY()[ipoint];
 382       if (y<yMin) yMin=y;
 383       if (y>yMax) yMax=y;
 384       fitterLocal.AddPoint(&dxl,y,1);
 385     }
 386
 387     fitterLocal.Eval();
 388     sigma2 += fitterLocal.GetChisquare()/Double_t((index1-indexM)-4.);
 389     TMatrixD   * covar = new ((*fCovars)[iKnot]) TMatrixD(4,4);
 390     TVectorD   * param = new ((*fParams)[iKnot]) TVectorD(4);
 391     fitterLocal.GetParameters(*param);
 392     fitterLocal.GetCovarianceMatrix(*covar);
 393     fitterLocal.ClearPoints();
 394   }
 395   fSigma  =TMath::Sqrt(sigma2/Double_t(fN0));   // mean sigma
 396   Double_t tDiff = ((yMax-yMin)+TMath::Abs(yMax)+TMath::Abs(yMin))*kEpsilon;
 397   fSigma += tDiff+fMaxDelta/TMath::Sqrt(npoints);
 398   fMaxDelta +=tDiff;
 399   for (Int_t iKnot=0; iKnot<fN0; iKnot++){
 400     TMatrixD & cov = *((TMatrixD*)fCovars->At(iKnot));
 401     cov*=fSigma*fSigma;
 402   }
 403   OptimizeKnots(iter);
 404
 405   fN = 0;
 406   for (Int_t iKnot=0; iKnot<fN0; iKnot++) if (fIndex[iKnot]>=0) fN++;
 407   fX  = new Double_t[fN];
 408   fY0 = new Double_t[fN];
 409   fY1 = new Double_t[fN];
 410   fChi2I = new Double_t[fN];
 411   Int_t iKnot=0;
 412   for (Int_t i=0; i<fN0; i++){
 413     if (fIndex[i]<0) continue;
 414     if (iKnot>=fN) {
 415       printf("AliSplineFit::InitKnots: Knot number > Max knot number\n");
 416       break;
 417     }
 418     TVectorD   * param = (TVectorD*) fParams->At(i);
 419     fX[iKnot]  = fGraph->GetX()[fIndex[i]];
 420     fY0[iKnot] = (*param)(0);
 421     fY1[iKnot] = (*param)(1);
 422     fChi2I[iKnot] = 0;
 423     iKnot++;
 424   }
 425 }
 426
 427
 428 Int_t AliSplineFit::OptimizeKnots(Int_t nIter){
 429   //
 430   //
 431   //
 432   const Double_t kMaxChi2= 5;
 433   Int_t nKnots=0;
 434   TTreeSRedirector cstream("SplineIter.root");
 435   for (Int_t iIter=0; iIter<nIter; iIter++){
 436     if (fBDump) cstream<<"Fit"<<
 437       "iIter="<<iIter<<
 438       "fit.="<<this<<
 439       "\n";
 440     nKnots=2;
 441     for (Int_t iKnot=1; iKnot<fN0-1; iKnot++){
 442       if (fIndex[iKnot]<0) continue;   //disabled knot
 443       Double_t chi2 = CheckKnot(iKnot);
 444       Double_t startX = fGraph->GetX()[fIndex[iKnot]];
 445       if (fBDump) {
 446         TMatrixD   * covar = (TMatrixD*)fCovars->At(iKnot);
 447         TVectorD   * param = (TVectorD*)fParams->At(iKnot);
 448         cstream<<"Chi2"<<
 449          "iIter="<<iIter<<
 450          "iKnot="<<iKnot<<
 451          "chi2="<<chi2<<
 452          "x="<<startX<<
 453          "param="<<param<<
 454          "covar="<<covar<<
 455          "\n";
 456       }
 457       if (chi2>kMaxChi2) { nKnots++;continue;}
 458       fIndex[iKnot]*=-1;
 459       Int_t iPrevious=iKnot-1;
 460       Int_t iNext    =iKnot+1;
 461       while (fIndex[iPrevious]<0) iPrevious--;
 462       while (fIndex[iNext]<0) iNext++;
 463       RefitKnot(iPrevious);
 464       RefitKnot(iNext);
 465       iKnot++;
 466       while (iKnot<fN0-1&& fIndex[iKnot]<0) iKnot++;
 467     }
 468   }
 469   return nKnots;
 470 }
 471
 472
 473 Bool_t   AliSplineFit::RefitKnot(Int_t iKnot){
 474   //
 475   //
 476   //
 477
 478   Int_t iPrevious=(iKnot>0)  ?iKnot-1: 0;
 479   Int_t iNext    =(iKnot<fN0)?iKnot+1: fN0-1;
 480   while (iPrevious>0&&fIndex[iPrevious]<0) iPrevious--;
 481   while (iNext<fN0&&fIndex[iNext]<0) iNext++;
 482   if (iPrevious<0) iPrevious=0;
 483   if (iNext>=fN0) iNext=fN0-1;
 484
 485   Double_t startX = fGraph->GetX()[fIndex[iKnot]];
 486   AliSplineFit::fitterStatic()->ClearPoints();
 487   Int_t indPrev = fIndex[iPrevious];
 488   Int_t indNext = fIndex[iNext];
 489   Double_t *graphX = fGraph->GetX();
 490   Double_t *graphY = fGraph->GetY();
 491
 492   // make arrays for points to fit (to save time)
 493
 494   Int_t nPoints = indNext-indPrev;
 495   Double_t *xPoint = new Double_t[3*nPoints];
 496   Double_t *yPoint = &xPoint[nPoints];
 497   Double_t *ePoint = &xPoint[2*nPoints];
 498   Int_t indVec=0;
 499   for (Int_t iPoint=indPrev; iPoint<indNext; iPoint++, indVec++){
 500     Double_t dxl   = graphX[iPoint]-startX;
 501     Double_t  y    = graphY[iPoint];
 502     xPoint[indVec] = dxl;
 503     yPoint[indVec] = y;
 504     ePoint[indVec] =  fSigma;
 505 //    ePoint[indVec] =  fSigma+TMath::Abs(y)*kEpsilon;
 506 //    AliSplineFit::fitterStatic.AddPoint(&dxl,y,fSigma+TMath::Abs(y)*kEpsilon);
 507   }
 508   AliSplineFit::fitterStatic()->AssignData(nPoints,1,xPoint,yPoint,ePoint);
 509   AliSplineFit::fitterStatic()->Eval();
 510
 511 //  delete temporary arrays
 512
 513   delete [] xPoint;
 514
 515   TMatrixD   * covar = (TMatrixD*)fCovars->At(iKnot);
 516   TVectorD   * param = (TVectorD*)fParams->At(iKnot);
 517   AliSplineFit::fitterStatic()->GetParameters(*param);
 518   AliSplineFit::fitterStatic()->GetCovarianceMatrix(*covar);
 519   return 0;
 520 }
 521
 522
 523 Float_t AliSplineFit::CheckKnot(Int_t iKnot){
 524   //
 525   //
 526   //
 527
 528   Int_t iPrevious=iKnot-1;
 529   Int_t iNext    =iKnot+1;
 530   while (fIndex[iPrevious]<0) iPrevious--;
 531   while (fIndex[iNext]<0) iNext++;
 532   TVectorD &pPrevious = *((TVectorD*)fParams->At(iPrevious));
 533   TVectorD &pNext     = *((TVectorD*)fParams->At(iNext));
 534   TVectorD &pKnot     = *((TVectorD*)fParams->At(iKnot));
 535   TMatrixD &cPrevious = *((TMatrixD*)fCovars->At(iPrevious));
 536   TMatrixD &cNext     = *((TMatrixD*)fCovars->At(iNext));
 537   TMatrixD &cKnot     = *((TMatrixD*)fCovars->At(iKnot));
 538   Double_t xPrevious = fGraph->GetX()[fIndex[iPrevious]];
 539   Double_t xNext     = fGraph->GetX()[fIndex[iNext]];
 540   Double_t xKnot     = fGraph->GetX()[fIndex[iKnot]];
 541
 542   // extra variables introduced to save processing time
 543
 544   Double_t dxc  = xNext-xPrevious;
 545   Double_t invDxc = 1./dxc;
 546   Double_t invDxc2 = invDxc*invDxc;
 547   TMatrixD  tPrevious(4,4);
 548   TMatrixD  tNext(4,4);
 549
 550   tPrevious(0,0) = 1;    tPrevious(1,1) = 1;
 551   tPrevious(2,0) = -3.*invDxc2;
 552   tPrevious(2,1) = -2.*invDxc;
 553   tPrevious(3,0) =  2.*invDxc2*invDxc;
 554   tPrevious(3,1) =  1.*invDxc2;
 555   tNext(2,0)     =  3.*invDxc2;      tNext(2,1)     = -1*invDxc;
 556   tNext(3,0)     = -2.*invDxc2*invDxc;  tNext(3,1)     =  1.*invDxc2;
 557   TMatrixD  tpKnot(4,4);
 558   TMatrixD  tpNext(4,4);
 559   Double_t dx = xKnot-xPrevious;
 560   tpKnot(0,0) = 1;      tpKnot(1,1) = 1;        tpKnot(2,2) = 1;        tpKnot(3,3) = 1;
 561   tpKnot(0,1) = dx;     tpKnot(0,2) = dx*dx;    tpKnot(0,3) = dx*dx*dx;
 562   tpKnot(1,2) = 2.*dx;  tpKnot(1,3) = 3.*dx*dx;
 563   tpKnot(2,3) = 3.*dx;
 564   Double_t dxn = xNext-xPrevious;
 565   tpNext(0,0) = 1;       tpNext(1,1) = 1;        tpNext(2,2) = 1;        tpNext(3,3) = 1;
 566   tpNext(0,1) = dxn;     tpNext(0,2) = dxn*dxn;    tpNext(0,3) = dxn*dxn*dxn;
 567   tpNext(1,2) = 2.*dxn;  tpNext(1,3) = 3.*dxn*dxn;
 568   tpNext(2,3) = 3.*dxn;
 569
 570   //
 571   // matrix and vector at previous
 572   //
 573
 574   TVectorD  sPrevious = tPrevious*pPrevious+tNext*pNext;
 575   TVectorD  sKnot     = tpKnot*sPrevious;
 576   TVectorD  sNext     = tpNext*sPrevious;
 577
 578   TMatrixD csPrevious00(tPrevious, TMatrixD::kMult,cPrevious);
 579   csPrevious00 *= tPrevious.T();
 580   TMatrixD csPrevious01(tNext,TMatrixD::kMult,cNext);
 581   csPrevious01*=tNext.T();
 582   TMatrixD  csPrevious(csPrevious00,TMatrixD::kPlus,csPrevious01);
 583   TMatrixD  csKnot(tpKnot,TMatrixD::kMult,csPrevious);
 584   csKnot*=tpKnot.T();
 585   TMatrixD  csNext(tpNext,TMatrixD::kMult,csPrevious);
 586   csNext*=tpNext.T();
 587
 588   TVectorD dPrevious = pPrevious-sPrevious;
 589   TVectorD dKnot     = pKnot-sKnot;
 590   TVectorD dNext     = pNext-sNext;
 591   //
 592   //
 593   TMatrixD prec(4,4);
 594   prec(0,0) = (fMaxDelta*fMaxDelta);
 595   prec(1,1) = prec(0,0)*invDxc2;
 596   prec(2,2) = prec(1,1)*invDxc2;
 597   prec(3,3) = prec(2,2)*invDxc2;
 598
 599 //   prec(0,0) = (fMaxDelta*fMaxDelta);
 600 //   prec(1,1) = (fMaxDelta*fMaxDelta)/(dxc*dxc);
 601 //   prec(2,2) = (fMaxDelta*fMaxDelta)/(dxc*dxc*dxc*dxc);
 602 //   prec(3,3) = (fMaxDelta*fMaxDelta)/(dxc*dxc*dxc*dxc*dxc*dxc);
 603
 604   csPrevious+=cPrevious;
 605   csPrevious+=prec;
 606   csPrevious.Invert();
 607   Double_t  chi2P     = dPrevious*(csPrevious*dPrevious);
 608
 609   csKnot+=cKnot;
 610   csKnot+=prec;
 611   csKnot.Invert();
 612   Double_t  chi2K     = dKnot*(csKnot*dKnot);
 613
 614   csNext+=cNext;
 615   csNext+=prec;
 616   csNext.Invert();
 617   Double_t  chi2N     = dNext*(csNext*dNext);
 618
 619   return (chi2P+chi2K+chi2N)/8.;
 620
 621
 622 }
 623
 624 void AliSplineFit::SplineFit(Int_t nder){
 625   //
 626   // Cubic spline fit of graph
 627   //
 628   // nder
 629   // nder<0  - no continuity requirement
 630   //     =0  - continous  0 derivative
 631   //     =1  - continous  1 derivative
 632   //     >1  - continous  2 derivative
 633   //
 634   if (!fGraph) return;
 635   TGraph * graph = fGraph;
 636   if (nder>1) nder=2;
 637   Int_t nknots  = fN;
 638   if (nknots < 2 ) return;
 639   Int_t npoints = graph->GetN();
 640   if (npoints < fMinPoints ) return;
 641   //
 642   //
 643   // spline fit
 644   // each knot 4 parameters
 645   //
 646   TMatrixD       *pmatrix = 0;
 647   TVectorD       *pvalues = 0;
 648   if (nder>1){
 649     pmatrix = new TMatrixD(4*(nknots-1)+3*(nknots-2), 4*(nknots-1)+3*(nknots-2));
 650     pvalues = new TVectorD(4*(nknots-1)+3*(nknots-2));
 651   }
 652   if (nder==1){
 653     pmatrix = new TMatrixD(4*(nknots-1)+2*(nknots-2), 4*(nknots-1)+2*(nknots-2));
 654     pvalues = new TVectorD(4*(nknots-1)+2*(nknots-2));
 655   }
 656   if (nder==0){
 657     pmatrix = new TMatrixD(4*(nknots-1)+1*(nknots-2), 4*(nknots-1)+1*(nknots-2));
 658     pvalues = new TVectorD(4*(nknots-1)+1*(nknots-2));
 659   }
 660   if (nder<0){
 661     pmatrix = new TMatrixD(4*(nknots-1)+0*(nknots-2), 4*(nknots-1)+0*(nknots-2));
 662     pvalues = new TVectorD(4*(nknots-1)+0*(nknots-2));
 663   }
 664
 665
 666   TMatrixD &matrix = *pmatrix;
 667   TVectorD &values = *pvalues;
 668   Int_t    current = 0;
 669 //
 670 //  defined extra variables (current4 etc.) to save processing time.
 671 //  fill normal matrices, then copy to sparse matrix.
 672 //
 673   Double_t *graphX = graph->GetX();
 674   Double_t *graphY = graph->GetY();
 675   for (Int_t ip=0;ip<npoints;ip++){
 676     if (current<nknots-2&&graphX[ip]>fX[current+1]) current++;
 677     Double_t xmiddle = (fX[current+1]+fX[current])*0.5;
 678     Double_t x1 = graphX[ip]- xmiddle;
 679     Double_t x2 = x1*x1;
 680     Double_t x3 = x2*x1;
 681     Double_t x4 = x2*x2;
 682     Double_t x5 = x3*x2;
 683     Double_t x6 = x3*x3;
 684     Double_t y  = graphY[ip];
 685     Int_t current4 = 4*current;
 686
 687     matrix(current4  , current4  )+=1;
 688     matrix(current4  , current4+1)+=x1;
 689     matrix(current4  , current4+2)+=x2;
 690     matrix(current4  , current4+3)+=x3;
 691     //
 692     matrix(current4+1, current4  )+=x1;
 693     matrix(current4+1, current4+1)+=x2;
 694     matrix(current4+1, current4+2)+=x3;
 695     matrix(current4+1, current4+3)+=x4;
 696     //
 697     matrix(current4+2, current4  )+=x2;
 698     matrix(current4+2, current4+1)+=x3;
 699     matrix(current4+2, current4+2)+=x4;
 700     matrix(current4+2, current4+3)+=x5;
 701     //
 702     matrix(current4+3, current4  )+=x3;
 703     matrix(current4+3, current4+1)+=x4;
 704     matrix(current4+3, current4+2)+=x5;
 705     matrix(current4+3, current4+3)+=x6;
 706     //
 707     values(current4  ) += y;
 708     values(current4+1) += y*x1;
 709     values(current4+2) += y*x2;
 710     values(current4+3) += y*x3;
 711   }
 712   //
 713   // constraint 0
 714   //
 715   Int_t offset =4*(nknots-1)-1;
 716   if (nder>=0) for (Int_t iknot = 1; iknot<nknots-1; iknot++){
 717
 718     Double_t dxm  =  (fX[iknot]-fX[iknot-1])*0.5;
 719     Double_t dxp  = -(fX[iknot+1]-fX[iknot])*0.5;
 720     Double_t dxm2 = dxm*dxm;
 721     Double_t dxp2 = dxp*dxp;
 722     Double_t dxm3 = dxm2*dxm;
 723     Double_t dxp3 = dxp2*dxp;
 724     Int_t iknot4  = 4*iknot;
 725     Int_t iknot41 = 4*(iknot-1);
 726     Int_t offsKnot = offset+iknot;
 727     //
 728     // condition on knot
 729     //
 730     // a0[i] = a0m[i-1]  + a1m[i-1]*dxm + a2m[i-1]*dxm^2 + a3m[i-1]*dxm^3
 731     // a0[i] = a0m[i-0]  + a1m[i-0]*dxp + a2m[i-0]*dxp^2 + a3m[i-0]*dxp^3
 732     // (a0m[i-1]  + a1m[i-1]*dxm + a2m[i-1]*dxm^2 + a3m[i-1]*dxm^3) -
 733     // (a0m[i-0]  + a1m[i-0]*dxp + a2m[i-0]*dxp^2 + a3m[i-0]*dxp^3)  = 0
 734
 735     matrix(offsKnot, iknot41  )=1;
 736     matrix(offsKnot, iknot4   )=-1;
 737
 738     matrix(offsKnot, iknot41+1)=dxm;
 739     matrix(offsKnot, iknot4 +1)=-dxp;
 740
 741     matrix(offsKnot, iknot41+2)=dxm2;
 742     matrix(offsKnot, iknot4 +2)=-dxp2;
 743
 744     matrix(offsKnot, iknot41+3)=dxm3;
 745     matrix(offsKnot, iknot4 +3)=-dxp3;
 746
 747     matrix(iknot41  , offsKnot)=1;
 748     matrix(iknot41+1, offsKnot)=dxm;
 749     matrix(iknot41+2, offsKnot)=dxm2;
 750     matrix(iknot41+3, offsKnot)=dxm3;
 751     matrix(iknot4  , offsKnot)=-1;
 752     matrix(iknot4+1, offsKnot)=-dxp;
 753     matrix(iknot4+2, offsKnot)=-dxp2;
 754     matrix(iknot4+3, offsKnot)=-dxp3;
 755   }
 756   //
 757   // constraint 1
 758   //
 759   offset =4*(nknots-1)-1+(nknots-2);
 760   if (nder>=1)for (Int_t iknot = 1; iknot<nknots-1; iknot++){
 761
 762     Double_t dxm  =  (fX[iknot]-fX[iknot-1])*0.5;
 763     Double_t dxp  = -(fX[iknot+1]-fX[iknot])*0.5;
 764     Double_t dxm2 = dxm*dxm;
 765     Double_t dxp2 = dxp*dxp;
 766     Int_t iknot4  = 4*iknot;
 767     Int_t iknot41 = 4*(iknot-1);
 768     Int_t offsKnot = offset+iknot;
 769     //
 770     // condition on knot derivation
 771     //
 772     // a0d[i] =  a1m[i-1] + 2*a2m[i-1]*dxm + 3*a3m[i-1]*dxm^2
 773     // a0d[i] =  a1m[i-0] + 2*a2m[i-0]*dxp + 3*a3m[i-0]*dxp^2
 774
 775     //
 776     matrix(offsKnot, iknot41+1)= 1;
 777     matrix(offsKnot, iknot4 +1)=-1;
 778
 779     matrix(offsKnot, iknot41+2)= 2.*dxm;
 780     matrix(offsKnot, iknot4 +2)=-2.*dxp;
 781
 782     matrix(offsKnot, iknot41+3)= 3.*dxm2;
 783     matrix(offsKnot, iknot4 +3)=-3.*dxp2;
 784
 785     matrix(iknot41+1, offsKnot)=1;
 786     matrix(iknot41+2, offsKnot)=2.*dxm;
 787     matrix(iknot41+3, offsKnot)=3.*dxm2;
 788
 789     matrix(iknot4+1, offsKnot)=-1.;
 790     matrix(iknot4+2, offsKnot)=-2.*dxp;
 791     matrix(iknot4+3, offsKnot)=-3.*dxp2;
 792   }
 793   //
 794   // constraint 2
 795   //
 796   offset =4*(nknots-1)-1+2*(nknots-2);
 797   if (nder>=2) for (Int_t iknot = 1; iknot<nknots-1; iknot++){
 798
 799     Double_t dxm  =  (fX[iknot]-fX[iknot-1])*0.5;
 800     Double_t dxp  = -(fX[iknot+1]-fX[iknot])*0.5;
 801     Int_t iknot4  = 4*iknot;
 802     Int_t iknot41 = 4*(iknot-1);
 803     Int_t offsKnot = offset+iknot;
 804     //
 805     // condition on knot second derivative
 806     //
 807     // a0dd[i] =  2*a2m[i-1] + 6*a3m[i-1]*dxm
 808     // a0dd[i] =  2*a2m[i-0] + 6*a3m[i-0]*dxp
 809     //
 810     //
 811     matrix(offsKnot, iknot41+2)= 2.;
 812     matrix(offsKnot, iknot4 +2)=-2.;
 813
 814     matrix(offsKnot, iknot41+3)= 6.*dxm;
 815     matrix(offsKnot, iknot4 +3)=-6.*dxp;
 816
 817     matrix(iknot41+2, offsKnot)=2.;
 818     matrix(iknot41+3, offsKnot)=6.*dxm;
 819
 820     matrix(iknot4+2, offsKnot)=-2.;
 821     matrix(iknot4+3, offsKnot)=-6.*dxp;
 822   }
 823
 824 // sparse matrix to do fit
 825
 826   TMatrixDSparse smatrix(matrix);
 827   TDecompSparse svd(smatrix,0);
 828   Bool_t ok;
 829   const TVectorD results = svd.Solve(values,ok);
 830
 831   for (Int_t iknot = 0; iknot<nknots-1; iknot++){
 832
 833     Double_t dxm  =  -(fX[iknot+1]-fX[iknot])*0.5;
 834
 835     fY0[iknot] = results(4*iknot)+ results(4*iknot+1)*dxm+results(4*iknot+2)*dxm*dxm+
 836       results(4*iknot+3)*dxm*dxm*dxm;
 837
 838     fY1[iknot] = results(4*iknot+1)+2.*results(4*iknot+2)*dxm+
 839       3*results(4*iknot+3)*dxm*dxm;
 840   }
 841   Int_t   iknot2= nknots-1;
 842   Int_t   iknot = nknots-2;
 843   Double_t dxm   =  (fX[iknot2]-fX[iknot2-1])*0.5;
 844
 845   fY0[iknot2] = results(4*iknot)+ results(4*iknot+1)*dxm+results(4*iknot+2)*dxm*dxm+
 846     results(4*iknot+3)*dxm*dxm*dxm;
 847
 848   fY1[iknot2] = results(4*iknot+1)+2.*results(4*iknot+2)*dxm+
 849       3*results(4*iknot+3)*dxm*dxm;
 850
 851   delete  pmatrix;
 852   delete  pvalues;
 853
 854 }
 855
 856
 857
 858
 859
 860 void AliSplineFit::MakeKnots0(TGraph * graph, Double_t maxdelta, Int_t minpoints){
 861   //
 862   // make knots  - restriction max distance and minimum points
 863   //
 864
 865   Int_t npoints  = graph->GetN();
 866   Double_t *xknots = new Double_t[npoints];
 867   for (Int_t ip=0;ip<npoints;ip++) xknots[ip]=0;
 868   //
 869   Int_t nknots =0;
 870   Int_t ipoints =0;
 871   //
 872   // generate knots
 873   //
 874   for (Int_t ip=0;ip<npoints;ip++){
 875     if (graph->GetX()[ip]-xknots[nknots-1]>maxdelta && ipoints>minpoints){
 876       xknots[nknots] = graph->GetX()[ip];
 877       ipoints=1;
 878       nknots++;
 879     }
 880     ipoints++;
 881   }
 882   if (npoints-ipoints>minpoints){
 883     xknots[nknots] = graph->GetX()[npoints-1];
 884     nknots++;
 885   }else{
 886     xknots[nknots-1] = graph->GetX()[npoints-1];
 887   }
 888
 889   fN = nknots;
 890   fX = new Double_t[nknots];
 891   fY0 = new Double_t[nknots];
 892   fY1 = new Double_t[nknots];
 893   fChi2I= new Double_t[nknots];
 894   for (Int_t i=0; i<nknots; i++) fX[i]= xknots[i];
 895   delete [] xknots;
 896 }
 897
 898
 899
 900
 901 void AliSplineFit::MakeSmooth(TGraph * graph, Float_t ratio, Option_t * type){
 902   //
 903   // Interface to GraphSmooth
 904   //
 905
 906   TGraphSmooth smooth;
 907   Int_t    npoints2 = TMath::Nint(graph->GetN()*ratio);
 908   TGraph * graphT0 = smooth.SmoothKern(graph,type,ratio);
 909   if (!graphT0) return;
 910   TGraph  graphT1(npoints2);
 911   for (Int_t ipoint=0; ipoint<npoints2; ipoint++){
 912     Int_t pointS = TMath::Nint(ipoint/ratio);
 913     if (ipoint==npoints2-1) pointS=graph->GetN()-1;
 914     graphT1.SetPoint(ipoint, graphT0->GetX()[pointS] , graphT0->GetY()[pointS]);
 915   }
 916   TSpline3 spline2("spline", &graphT1);
 917   Update(&spline2, npoints2);
 918 }
 919
 920
 921 void AliSplineFit::Update(TSpline3 *spline, Int_t nknots){
 922   //
 923   //
 924   //
 925
 926   fN = nknots;
 927   fX = new Double_t[nknots];
 928   fY0 = new Double_t[nknots];
 929   fY1 = new Double_t[nknots];
 930   Double_t d0, d1;
 931   fChi2I= 0;
 932   for (Int_t i=0; i<nknots; i++) {
 933     spline->GetCoeff(i,fX[i],fY0[i], fY1[i],d0,d1);
 934   }
 935 }
 936
 937
 938
 939
 940 void AliSplineFit::Test(Int_t npoints, Int_t ntracks, Float_t snoise){
 941   //
 942   // test function
 943   //
 944
 945   AliSplineFit fit;
 946   AliSplineFit fitS;
 947   TGraph * graph0=0;
 948   TGraph * graph1=0;
 949
 950   TTreeSRedirector *pcstream = new TTreeSRedirector("TestSmooth.root");
 951   for (Int_t i=0; i<ntracks; i++){
 952     graph0 = AliSplineFit::GenerGraph(npoints,0.05,0,0,1,0);
 953     graph1 = AliSplineFit::GenerNoise(graph0,snoise);
 954     fit.InitKnots(graph1, 10,10, 0.00);
 955     TGraph *d0 = fit.MakeDiff(graph0);
 956     TGraph *g0 = fit.MakeGraph(0,1,1000,0);
 957     fit.SplineFit(2);
 958     TH1F * h2 = fit.MakeDiffHisto(graph0);
 959     TGraph *d2 = fit.MakeDiff(graph0);
 960     TGraph *g2 = fit.MakeGraph(0,1,1000,0);
 961     fit.SplineFit(1);
 962     TH1F * h1 = fit.MakeDiffHisto(graph0);
 963     TGraph *d1 = fit.MakeDiff(graph0);
 964     TGraph *g1 = fit.MakeGraph(0,1,1000,0);
 965
 966     Float_t ratio = Float_t(fit.fN)/Float_t(npoints);
 967     fitS.MakeSmooth(graph1,ratio,"box");
 968     TGraph *dS = fitS.MakeDiff(graph0);
 969     TGraph *gS = fit.MakeGraph(0,1,1000,0);
 970
 971     TH1F * hS = fitS.MakeDiffHisto(graph0);
 972     Double_t mean2  = h2->GetMean();
 973     Double_t sigma2 = h2->GetRMS();
 974     Double_t mean1  = h1->GetMean();
 975     Double_t sigma1 = h1->GetRMS();
 976     Double_t meanS  = hS->GetMean();
 977     Double_t sigmaS = hS->GetRMS();
 978     char fname[100];
 979     if (fit.fN<20){
 980       snprintf(fname,100,"pol%d",fit.fN);
 981     }else{
 982       snprintf(fname,100,"pol%d",19);
 983     }
 984     TF1 fpol("fpol",fname);
 985     graph1->Fit(&fpol);
 986     TGraph dpol(*graph1);
 987     TGraph gpol(*graph1);
 988     for (Int_t ipoint=0; ipoint<graph1->GetN(); ipoint++){
 989       dpol.GetY()[ipoint]= graph0->GetY()[ipoint]-
 990         fpol.Eval(graph0->GetX()[ipoint]);
 991       gpol.GetY()[ipoint]= fpol.Eval(graph0->GetX()[ipoint]);
 992     }
 993     (*pcstream)<<"Test"<<
 994       "Event="<<i<<
 995       "Graph0.="<<graph0<<
 996       "Graph1.="<<graph1<<
 997       "G0.="<<g0<<
 998       "G1.="<<g1<<
 999       "G2.="<<g2<<
1000       "GS.="<<gS<<
1001       "GP.="<<&gpol<<
1002       "D0.="<<d0<<
1003       "D1.="<<d1<<
1004       "D2.="<<d2<<
1005       "DS.="<<dS<<
1006       "DP.="<<&dpol<<
1007       "Npoints="<<fit.fN<<
1008       "Mean1="<<mean1<<
1009       "Mean2="<<mean2<<
1010       "MeanS="<<meanS<<
1011       "Sigma1="<<sigma1<<
1012       "Sigma2="<<sigma2<<
1013       "SigmaS="<<sigmaS<<
1014       "\n";
1015
1016     delete graph0;
1017     delete graph1;
1018     delete g1;
1019     delete g2;
1020     delete gS;
1021     delete h1;
1022     delete h2;
1023     delete hS;
1024   }
1025   delete pcstream;
1026 }
1027
1028 void AliSplineFit::Cleanup(){
1029  //
1030  // deletes extra information to reduce amount of information stored on the data
1031  // base
1032
1033  //    delete fGraph;  fGraph=0;   // Don't delete fGraph -- input parameter
1034      delete fParams; fParams=0;
1035      delete fCovars; fCovars=0;
1036      delete [] fIndex; fIndex=0;
1037      delete [] fChi2I; fChi2I=0;
1038 }
1039
1040
1041 void AliSplineFit::CopyGraph() {
1042   //
1043   // enter graph points directly to fit parameters
1044   // (to bee used when too few points are available)
1045   //
1046   fN = fGraph->GetN();
1047   fX = new Double_t[fN];
1048   fY0 = new Double_t[fN];
1049   fY1 = new Double_t[fN];
1050   for (Int_t i=0; i<fN; i++ ) {
1051     fX[i]  = fGraph->GetX()[i];
1052     fY0[i] = fGraph->GetY()[i];
1053     fY1[i] = 0;
1054   }
1055 }