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