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   //
 636   //
 637   // spline fit
 638   // each knot 4 parameters
 639   //
 640   TMatrixD       *pmatrix = 0;
 641   TVectorD       *pvalues = 0;
 642   if (nder>1){
 643     pmatrix = new TMatrixD(4*(nknots-1)+3*(nknots-2), 4*(nknots-1)+3*(nknots-2));
 644     pvalues = new TVectorD(4*(nknots-1)+3*(nknots-2));
 645   }
 646   if (nder==1){
 647     pmatrix = new TMatrixD(4*(nknots-1)+2*(nknots-2), 4*(nknots-1)+2*(nknots-2));
 648     pvalues = new TVectorD(4*(nknots-1)+2*(nknots-2));
 649   }
 650   if (nder==0){
 651     pmatrix = new TMatrixD(4*(nknots-1)+1*(nknots-2), 4*(nknots-1)+1*(nknots-2));
 652     pvalues = new TVectorD(4*(nknots-1)+1*(nknots-2));
 653   }
 654   if (nder<0){
 655     pmatrix = new TMatrixD(4*(nknots-1)+0*(nknots-2), 4*(nknots-1)+0*(nknots-2));
 656     pvalues = new TVectorD(4*(nknots-1)+0*(nknots-2));
 657   }
 658
 659
 660   TMatrixD &matrix = *pmatrix;
 661   TVectorD &values = *pvalues;
 662   Int_t    current = 0;
 663 //
 664 //  defined extra variables (current4 etc.) to save processing time.
 665 //  fill normal matrices, then copy to sparse matrix.
 666 //
 667   Double_t *graphX = graph->GetX();
 668   Double_t *graphY = graph->GetY();
 669   for (Int_t ip=0;ip<npoints;ip++){
 670     if (current<nknots-2&&graphX[ip]>fX[current+1]) current++;
 671     Double_t xmiddle = (fX[current+1]+fX[current])*0.5;
 672     Double_t x1 = graphX[ip]- xmiddle;
 673     Double_t x2 = x1*x1;
 674     Double_t x3 = x2*x1;
 675     Double_t x4 = x2*x2;
 676     Double_t x5 = x3*x2;
 677     Double_t x6 = x3*x3;
 678     Double_t y  = graphY[ip];
 679     Int_t current4 = 4*current;
 680
 681     matrix(current4  , current4  )+=1;
 682     matrix(current4  , current4+1)+=x1;
 683     matrix(current4  , current4+2)+=x2;
 684     matrix(current4  , current4+3)+=x3;
 685     //
 686     matrix(current4+1, current4  )+=x1;
 687     matrix(current4+1, current4+1)+=x2;
 688     matrix(current4+1, current4+2)+=x3;
 689     matrix(current4+1, current4+3)+=x4;
 690     //
 691     matrix(current4+2, current4  )+=x2;
 692     matrix(current4+2, current4+1)+=x3;
 693     matrix(current4+2, current4+2)+=x4;
 694     matrix(current4+2, current4+3)+=x5;
 695     //
 696     matrix(current4+3, current4  )+=x3;
 697     matrix(current4+3, current4+1)+=x4;
 698     matrix(current4+3, current4+2)+=x5;
 699     matrix(current4+3, current4+3)+=x6;
 700     //
 701     values(current4  ) += y;
 702     values(current4+1) += y*x1;
 703     values(current4+2) += y*x2;
 704     values(current4+3) += y*x3;
 705   }
 706   //
 707   // constraint 0
 708   //
 709   Int_t offset =4*(nknots-1)-1;
 710   if (nder>=0) for (Int_t iknot = 1; iknot<nknots-1; iknot++){
 711
 712     Double_t dxm  =  (fX[iknot]-fX[iknot-1])*0.5;
 713     Double_t dxp  = -(fX[iknot+1]-fX[iknot])*0.5;
 714     Double_t dxm2 = dxm*dxm;
 715     Double_t dxp2 = dxp*dxp;
 716     Double_t dxm3 = dxm2*dxm;
 717     Double_t dxp3 = dxp2*dxp;
 718     Int_t iknot4  = 4*iknot;
 719     Int_t iknot41 = 4*(iknot-1);
 720     Int_t offsKnot = offset+iknot;
 721     //
 722     // condition on knot
 723     //
 724     // a0[i] = a0m[i-1]  + a1m[i-1]*dxm + a2m[i-1]*dxm^2 + a3m[i-1]*dxm^3
 725     // a0[i] = a0m[i-0]  + a1m[i-0]*dxp + a2m[i-0]*dxp^2 + a3m[i-0]*dxp^3
 726     // (a0m[i-1]  + a1m[i-1]*dxm + a2m[i-1]*dxm^2 + a3m[i-1]*dxm^3) -
 727     // (a0m[i-0]  + a1m[i-0]*dxp + a2m[i-0]*dxp^2 + a3m[i-0]*dxp^3)  = 0
 728
 729     matrix(offsKnot, iknot41  )=1;
 730     matrix(offsKnot, iknot4   )=-1;
 731
 732     matrix(offsKnot, iknot41+1)=dxm;
 733     matrix(offsKnot, iknot4 +1)=-dxp;
 734
 735     matrix(offsKnot, iknot41+2)=dxm2;
 736     matrix(offsKnot, iknot4 +2)=-dxp2;
 737
 738     matrix(offsKnot, iknot41+3)=dxm3;
 739     matrix(offsKnot, iknot4 +3)=-dxp3;
 740
 741     matrix(iknot41  , offsKnot)=1;
 742     matrix(iknot41+1, offsKnot)=dxm;
 743     matrix(iknot41+2, offsKnot)=dxm2;
 744     matrix(iknot41+3, offsKnot)=dxm3;
 745     matrix(iknot4  , offsKnot)=-1;
 746     matrix(iknot4+1, offsKnot)=-dxp;
 747     matrix(iknot4+2, offsKnot)=-dxp2;
 748     matrix(iknot4+3, offsKnot)=-dxp3;
 749   }
 750   //
 751   // constraint 1
 752   //
 753   offset =4*(nknots-1)-1+(nknots-2);
 754   if (nder>=1)for (Int_t iknot = 1; iknot<nknots-1; iknot++){
 755
 756     Double_t dxm  =  (fX[iknot]-fX[iknot-1])*0.5;
 757     Double_t dxp  = -(fX[iknot+1]-fX[iknot])*0.5;
 758     Double_t dxm2 = dxm*dxm;
 759     Double_t dxp2 = dxp*dxp;
 760     Int_t iknot4  = 4*iknot;
 761     Int_t iknot41 = 4*(iknot-1);
 762     Int_t offsKnot = offset+iknot;
 763     //
 764     // condition on knot derivation
 765     //
 766     // a0d[i] =  a1m[i-1] + 2*a2m[i-1]*dxm + 3*a3m[i-1]*dxm^2
 767     // a0d[i] =  a1m[i-0] + 2*a2m[i-0]*dxp + 3*a3m[i-0]*dxp^2
 768
 769     //
 770     matrix(offsKnot, iknot41+1)= 1;
 771     matrix(offsKnot, iknot4 +1)=-1;
 772
 773     matrix(offsKnot, iknot41+2)= 2.*dxm;
 774     matrix(offsKnot, iknot4 +2)=-2.*dxp;
 775
 776     matrix(offsKnot, iknot41+3)= 3.*dxm2;
 777     matrix(offsKnot, iknot4 +3)=-3.*dxp2;
 778
 779     matrix(iknot41+1, offsKnot)=1;
 780     matrix(iknot41+2, offsKnot)=2.*dxm;
 781     matrix(iknot41+3, offsKnot)=3.*dxm2;
 782
 783     matrix(iknot4+1, offsKnot)=-1.;
 784     matrix(iknot4+2, offsKnot)=-2.*dxp;
 785     matrix(iknot4+3, offsKnot)=-3.*dxp2;
 786   }
 787   //
 788   // constraint 2
 789   //
 790   offset =4*(nknots-1)-1+2*(nknots-2);
 791   if (nder>=2) for (Int_t iknot = 1; iknot<nknots-1; iknot++){
 792
 793     Double_t dxm  =  (fX[iknot]-fX[iknot-1])*0.5;
 794     Double_t dxp  = -(fX[iknot+1]-fX[iknot])*0.5;
 795     Int_t iknot4  = 4*iknot;
 796     Int_t iknot41 = 4*(iknot-1);
 797     Int_t offsKnot = offset+iknot;
 798     //
 799     // condition on knot second derivative
 800     //
 801     // a0dd[i] =  2*a2m[i-1] + 6*a3m[i-1]*dxm
 802     // a0dd[i] =  2*a2m[i-0] + 6*a3m[i-0]*dxp
 803     //
 804     //
 805     matrix(offsKnot, iknot41+2)= 2.;
 806     matrix(offsKnot, iknot4 +2)=-2.;
 807
 808     matrix(offsKnot, iknot41+3)= 6.*dxm;
 809     matrix(offsKnot, iknot4 +3)=-6.*dxp;
 810
 811     matrix(iknot41+2, offsKnot)=2.;
 812     matrix(iknot41+3, offsKnot)=6.*dxm;
 813
 814     matrix(iknot4+2, offsKnot)=-2.;
 815     matrix(iknot4+3, offsKnot)=-6.*dxp;
 816   }
 817
 818 // sparse matrix to do fit
 819
 820   TMatrixDSparse smatrix(matrix);
 821   TDecompSparse svd(smatrix,0);
 822   Bool_t ok;
 823   const TVectorD results = svd.Solve(values,ok);
 824
 825   for (Int_t iknot = 0; iknot<nknots-1; iknot++){
 826
 827     Double_t dxm  =  -(fX[iknot+1]-fX[iknot])*0.5;
 828
 829     fY0[iknot] = results(4*iknot)+ results(4*iknot+1)*dxm+results(4*iknot+2)*dxm*dxm+
 830       results(4*iknot+3)*dxm*dxm*dxm;
 831
 832     fY1[iknot] = results(4*iknot+1)+2.*results(4*iknot+2)*dxm+
 833       3*results(4*iknot+3)*dxm*dxm;
 834   }
 835   Int_t   iknot2= nknots-1;
 836   Int_t   iknot = nknots-2;
 837   Double_t dxm   =  (fX[iknot2]-fX[iknot2-1])*0.5;
 838
 839   fY0[iknot2] = results(4*iknot)+ results(4*iknot+1)*dxm+results(4*iknot+2)*dxm*dxm+
 840     results(4*iknot+3)*dxm*dxm*dxm;
 841
 842   fY1[iknot2] = results(4*iknot+1)+2.*results(4*iknot+2)*dxm+
 843       3*results(4*iknot+3)*dxm*dxm;
 844
 845   delete  pmatrix;
 846   delete  pvalues;
 847
 848 }
 849
 850
 851
 852
 853
 854 void AliSplineFit::MakeKnots0(TGraph * graph, Double_t maxdelta, Int_t minpoints){
 855   //
 856   // make knots  - restriction max distance and minimum points
 857   //
 858
 859   Int_t npoints  = graph->GetN();
 860   Double_t *xknots = new Double_t[npoints];
 861   Int_t nknots =0;
 862   Int_t ipoints =0;
 863   //
 864   // generate knots
 865   //
 866   for (Int_t ip=0;ip<npoints;ip++){
 867     if (graph->GetX()[ip]-xknots[nknots-1]>maxdelta && ipoints>minpoints){
 868       xknots[nknots] = graph->GetX()[ip];
 869       ipoints=1;
 870       nknots++;
 871     }
 872     ipoints++;
 873   }
 874   if (npoints-ipoints>minpoints){
 875     xknots[nknots] = graph->GetX()[npoints-1];
 876     nknots++;
 877   }else{
 878     xknots[nknots-1] = graph->GetX()[npoints-1];
 879   }
 880
 881   fN = nknots;
 882   fX = new Double_t[nknots];
 883   fY0 = new Double_t[nknots];
 884   fY1 = new Double_t[nknots];
 885   fChi2I= new Double_t[nknots];
 886   for (Int_t i=0; i<nknots; i++) fX[i]= xknots[i];
 887   delete [] xknots;
 888 }
 889
 890
 891
 892
 893 void AliSplineFit::MakeSmooth(TGraph * graph, Float_t ratio, char * type){
 894   //
 895   // Interface to GraphSmooth
 896   //
 897
 898   TGraphSmooth smooth;
 899   Int_t    npoints2 = TMath::Nint(graph->GetN()*ratio);
 900   TGraph * graphT0 = smooth.SmoothKern(graph,type,ratio);
 901   if (!graphT0) return;
 902   TGraph  graphT1(npoints2);
 903   for (Int_t ipoint=0; ipoint<npoints2; ipoint++){
 904     Int_t pointS = TMath::Nint(ipoint/ratio);
 905     if (ipoint==npoints2-1) pointS=graph->GetN()-1;
 906     graphT1.SetPoint(ipoint, graphT0->GetX()[pointS] , graphT0->GetY()[pointS]);
 907   }
 908   TSpline3 spline2("spline", &graphT1);
 909   Update(&spline2, npoints2);
 910 }
 911
 912
 913 void AliSplineFit::Update(TSpline3 *spline, Int_t nknots){
 914   //
 915   //
 916   //
 917
 918   fN = nknots;
 919   fX = new Double_t[nknots];
 920   fY0 = new Double_t[nknots];
 921   fY1 = new Double_t[nknots];
 922   Double_t d0, d1;
 923   fChi2I= 0;
 924   for (Int_t i=0; i<nknots; i++) {
 925     spline->GetCoeff(i,fX[i],fY0[i], fY1[i],d0,d1);
 926   }
 927 }
 928
 929
 930
 931
 932 void AliSplineFit::Test(Int_t npoints, Int_t ntracks, Float_t snoise){
 933   //
 934   // test function
 935   //
 936
 937   AliSplineFit fit;
 938   AliSplineFit fitS;
 939   TGraph * graph0=0;
 940   TGraph * graph1=0;
 941
 942   TTreeSRedirector *pcstream = new TTreeSRedirector("TestSmooth.root");
 943   for (Int_t i=0; i<ntracks; i++){
 944     graph0 = AliSplineFit::GenerGraph(npoints,0.05,0,0,1,0);
 945     graph1 = AliSplineFit::GenerNoise(graph0,snoise);
 946     fit.InitKnots(graph1, 10,10, 0.00);
 947     TGraph *d0 = fit.MakeDiff(graph0);
 948     TGraph *g0 = fit.MakeGraph(0,1,1000,0);
 949     fit.SplineFit(2);
 950     TH1F * h2 = fit.MakeDiffHisto(graph0);
 951     TGraph *d2 = fit.MakeDiff(graph0);
 952     TGraph *g2 = fit.MakeGraph(0,1,1000,0);
 953     fit.SplineFit(1);
 954     TH1F * h1 = fit.MakeDiffHisto(graph0);
 955     TGraph *d1 = fit.MakeDiff(graph0);
 956     TGraph *g1 = fit.MakeGraph(0,1,1000,0);
 957
 958     Float_t ratio = Float_t(fit.fN)/Float_t(npoints);
 959     fitS.MakeSmooth(graph1,ratio,"box");
 960     TGraph *dS = fitS.MakeDiff(graph0);
 961     TGraph *gS = fit.MakeGraph(0,1,1000,0);
 962
 963     TH1F * hS = fitS.MakeDiffHisto(graph0);
 964     Double_t mean2  = h2->GetMean();
 965     Double_t sigma2 = h2->GetRMS();
 966     Double_t mean1  = h1->GetMean();
 967     Double_t sigma1 = h1->GetRMS();
 968     Double_t meanS  = hS->GetMean();
 969     Double_t sigmaS = hS->GetRMS();
 970     char fname[100];
 971     if (fit.fN<20){
 972       sprintf(fname,"pol%d",fit.fN);
 973     }else{
 974       sprintf(fname,"pol%d",19);
 975     }
 976     TF1 fpol("fpol",fname);
 977     graph1->Fit(&fpol);
 978     TGraph dpol(*graph1);
 979     TGraph gpol(*graph1);
 980     for (Int_t ipoint=0; ipoint<graph1->GetN(); ipoint++){
 981       dpol.GetY()[ipoint]= graph0->GetY()[ipoint]-
 982         fpol.Eval(graph0->GetX()[ipoint]);
 983       gpol.GetY()[ipoint]= fpol.Eval(graph0->GetX()[ipoint]);
 984     }
 985     (*pcstream)<<"Test"<<
 986       "Event="<<i<<
 987       "Graph0.="<<graph0<<
 988       "Graph1.="<<graph1<<
 989       "G0.="<<g0<<
 990       "G1.="<<g1<<
 991       "G2.="<<g2<<
 992       "GS.="<<gS<<
 993       "GP.="<<&gpol<<
 994       "D0.="<<d0<<
 995       "D1.="<<d1<<
 996       "D2.="<<d2<<
 997       "DS.="<<dS<<
 998       "DP.="<<&dpol<<
 999       "Npoints="<<fit.fN<<
1000       "Mean1="<<mean1<<
1001       "Mean2="<<mean2<<
1002       "MeanS="<<meanS<<
1003       "Sigma1="<<sigma1<<
1004       "Sigma2="<<sigma2<<
1005       "SigmaS="<<sigmaS<<
1006       "\n";
1007
1008     delete graph0;
1009     delete graph1;
1010     delete g1;
1011     delete g2;
1012     delete gS;
1013     delete h1;
1014     delete h2;
1015     delete hS;
1016   }
1017   delete pcstream;
1018 }
1019
1020 void AliSplineFit::Cleanup(){
1021  //
1022  // deletes extra information to reduce amount of information stored on the data
1023  // base
1024
1025      delete fGraph;  fGraph=0;
1026      delete fParams; fParams=0;
1027      delete fCovars; fCovars=0;
1028      delete [] fIndex; fIndex=0;
1029      delete [] fChi2I; fChi2I=0;
1030 }
1031
1032
1033 void AliSplineFit::CopyGraph() {
1034   //
1035   // enter graph points directly to fit parameters
1036   // (to bee used when too few points are available)
1037   //
1038   fN = fGraph->GetN();
1039   fX = new Double_t[fN];
1040   fY0 = new Double_t[fN];
1041   fY1 = new Double_t[fN];
1042   for (Int_t i=0; i<fN; i++ ) {
1043     fX[i]  = fGraph->GetX()[i];
1044     fY0[i] = fGraph->GetY()[i];
1045     fY1[i] = 0;
1046   }
1047 }