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