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