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