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