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