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