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