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