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