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