STEER/CDB/AliSplineFit.cxx

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