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