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