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