]> 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
Switch off branches by hand.
[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;
633 Int_t npoints = graph->GetN();
634 //
635 //
636 // spline fit
637 // each knot 4 parameters
52fdcd41 638 //
0dd3a2ac 639 TMatrixD *pmatrix = 0;
640 TVectorD *pvalues = 0;
641 if (nder>1){
642 pmatrix = new TMatrixD(4*(nknots-1)+3*(nknots-2), 4*(nknots-1)+3*(nknots-2));
643 pvalues = new TVectorD(4*(nknots-1)+3*(nknots-2));
644 }
645 if (nder==1){
646 pmatrix = new TMatrixD(4*(nknots-1)+2*(nknots-2), 4*(nknots-1)+2*(nknots-2));
647 pvalues = new TVectorD(4*(nknots-1)+2*(nknots-2));
648 }
649 if (nder==0){
650 pmatrix = new TMatrixD(4*(nknots-1)+1*(nknots-2), 4*(nknots-1)+1*(nknots-2));
651 pvalues = new TVectorD(4*(nknots-1)+1*(nknots-2));
652 }
653 if (nder<0){
654 pmatrix = new TMatrixD(4*(nknots-1)+0*(nknots-2), 4*(nknots-1)+0*(nknots-2));
655 pvalues = new TVectorD(4*(nknots-1)+0*(nknots-2));
656 }
657
658
659 TMatrixD &matrix = *pmatrix;
660 TVectorD &values = *pvalues;
661 Int_t current = 0;
662//
663// defined extra variables (current4 etc.) to save processing time.
664// fill normal matrices, then copy to sparse matrix.
665//
666 Double_t *graphX = graph->GetX();
667 Double_t *graphY = graph->GetY();
668 for (Int_t ip=0;ip<npoints;ip++){
669 if (current<nknots-2&&graphX[ip]>fX[current+1]) current++;
670 Double_t xmiddle = (fX[current+1]+fX[current])*0.5;
671 Double_t x1 = graphX[ip]- xmiddle;
672 Double_t x2 = x1*x1;
673 Double_t x3 = x2*x1;
674 Double_t x4 = x2*x2;
675 Double_t x5 = x3*x2;
676 Double_t x6 = x3*x3;
677 Double_t y = graphY[ip];
678 Int_t current4 = 4*current;
679
680 matrix(current4 , current4 )+=1;
681 matrix(current4 , current4+1)+=x1;
682 matrix(current4 , current4+2)+=x2;
683 matrix(current4 , current4+3)+=x3;
684 //
685 matrix(current4+1, current4 )+=x1;
686 matrix(current4+1, current4+1)+=x2;
687 matrix(current4+1, current4+2)+=x3;
688 matrix(current4+1, current4+3)+=x4;
689 //
690 matrix(current4+2, current4 )+=x2;
691 matrix(current4+2, current4+1)+=x3;
692 matrix(current4+2, current4+2)+=x4;
693 matrix(current4+2, current4+3)+=x5;
694 //
695 matrix(current4+3, current4 )+=x3;
696 matrix(current4+3, current4+1)+=x4;
697 matrix(current4+3, current4+2)+=x5;
698 matrix(current4+3, current4+3)+=x6;
699 //
700 values(current4 ) += y;
701 values(current4+1) += y*x1;
702 values(current4+2) += y*x2;
703 values(current4+3) += y*x3;
704 }
705 //
706 // constraint 0
707 //
708 Int_t offset =4*(nknots-1)-1;
709 if (nder>=0) for (Int_t iknot = 1; iknot<nknots-1; iknot++){
710
711 Double_t dxm = (fX[iknot]-fX[iknot-1])*0.5;
712 Double_t dxp = -(fX[iknot+1]-fX[iknot])*0.5;
713 Double_t dxm2 = dxm*dxm;
714 Double_t dxp2 = dxp*dxp;
715 Double_t dxm3 = dxm2*dxm;
716 Double_t dxp3 = dxp2*dxp;
717 Int_t iknot4 = 4*iknot;
718 Int_t iknot41 = 4*(iknot-1);
719 Int_t offsKnot = offset+iknot;
720 //
721 // condition on knot
722 //
723 // a0[i] = a0m[i-1] + a1m[i-1]*dxm + a2m[i-1]*dxm^2 + a3m[i-1]*dxm^3
724 // a0[i] = a0m[i-0] + a1m[i-0]*dxp + a2m[i-0]*dxp^2 + a3m[i-0]*dxp^3
725 // (a0m[i-1] + a1m[i-1]*dxm + a2m[i-1]*dxm^2 + a3m[i-1]*dxm^3) -
726 // (a0m[i-0] + a1m[i-0]*dxp + a2m[i-0]*dxp^2 + a3m[i-0]*dxp^3) = 0
727
728 matrix(offsKnot, iknot41 )=1;
729 matrix(offsKnot, iknot4 )=-1;
730
731 matrix(offsKnot, iknot41+1)=dxm;
732 matrix(offsKnot, iknot4 +1)=-dxp;
733
734 matrix(offsKnot, iknot41+2)=dxm2;
735 matrix(offsKnot, iknot4 +2)=-dxp2;
736
737 matrix(offsKnot, iknot41+3)=dxm3;
738 matrix(offsKnot, iknot4 +3)=-dxp3;
739
740 matrix(iknot41 , offsKnot)=1;
741 matrix(iknot41+1, offsKnot)=dxm;
742 matrix(iknot41+2, offsKnot)=dxm2;
743 matrix(iknot41+3, offsKnot)=dxm3;
744 matrix(iknot4 , offsKnot)=-1;
745 matrix(iknot4+1, offsKnot)=-dxp;
746 matrix(iknot4+2, offsKnot)=-dxp2;
747 matrix(iknot4+3, offsKnot)=-dxp3;
748 }
749 //
750 // constraint 1
751 //
752 offset =4*(nknots-1)-1+(nknots-2);
753 if (nder>=1)for (Int_t iknot = 1; iknot<nknots-1; iknot++){
754
755 Double_t dxm = (fX[iknot]-fX[iknot-1])*0.5;
756 Double_t dxp = -(fX[iknot+1]-fX[iknot])*0.5;
757 Double_t dxm2 = dxm*dxm;
758 Double_t dxp2 = dxp*dxp;
759 Int_t iknot4 = 4*iknot;
760 Int_t iknot41 = 4*(iknot-1);
761 Int_t offsKnot = offset+iknot;
762 //
763 // condition on knot derivation
764 //
765 // a0d[i] = a1m[i-1] + 2*a2m[i-1]*dxm + 3*a3m[i-1]*dxm^2
766 // a0d[i] = a1m[i-0] + 2*a2m[i-0]*dxp + 3*a3m[i-0]*dxp^2
767
768 //
769 matrix(offsKnot, iknot41+1)= 1;
770 matrix(offsKnot, iknot4 +1)=-1;
771
772 matrix(offsKnot, iknot41+2)= 2.*dxm;
773 matrix(offsKnot, iknot4 +2)=-2.*dxp;
774
775 matrix(offsKnot, iknot41+3)= 3.*dxm2;
776 matrix(offsKnot, iknot4 +3)=-3.*dxp2;
777
778 matrix(iknot41+1, offsKnot)=1;
779 matrix(iknot41+2, offsKnot)=2.*dxm;
780 matrix(iknot41+3, offsKnot)=3.*dxm2;
781
782 matrix(iknot4+1, offsKnot)=-1.;
783 matrix(iknot4+2, offsKnot)=-2.*dxp;
784 matrix(iknot4+3, offsKnot)=-3.*dxp2;
785 }
786 //
787 // constraint 2
788 //
789 offset =4*(nknots-1)-1+2*(nknots-2);
790 if (nder>=2) for (Int_t iknot = 1; iknot<nknots-1; iknot++){
791
792 Double_t dxm = (fX[iknot]-fX[iknot-1])*0.5;
793 Double_t dxp = -(fX[iknot+1]-fX[iknot])*0.5;
794 Int_t iknot4 = 4*iknot;
795 Int_t iknot41 = 4*(iknot-1);
796 Int_t offsKnot = offset+iknot;
797 //
798 // condition on knot second derivative
799 //
800 // a0dd[i] = 2*a2m[i-1] + 6*a3m[i-1]*dxm
801 // a0dd[i] = 2*a2m[i-0] + 6*a3m[i-0]*dxp
802 //
803 //
804 matrix(offsKnot, iknot41+2)= 2.;
805 matrix(offsKnot, iknot4 +2)=-2.;
806
807 matrix(offsKnot, iknot41+3)= 6.*dxm;
808 matrix(offsKnot, iknot4 +3)=-6.*dxp;
809
810 matrix(iknot41+2, offsKnot)=2.;
811 matrix(iknot41+3, offsKnot)=6.*dxm;
812
813 matrix(iknot4+2, offsKnot)=-2.;
814 matrix(iknot4+3, offsKnot)=-6.*dxp;
815 }
816
817// sparse matrix to do fit
818
819 TMatrixDSparse smatrix(matrix);
820 TDecompSparse svd(smatrix,0);
821 Bool_t ok;
822 const TVectorD results = svd.Solve(values,ok);
823
824 for (Int_t iknot = 0; iknot<nknots-1; iknot++){
825
826 Double_t dxm = -(fX[iknot+1]-fX[iknot])*0.5;
827
828 fY0[iknot] = results(4*iknot)+ results(4*iknot+1)*dxm+results(4*iknot+2)*dxm*dxm+
829 results(4*iknot+3)*dxm*dxm*dxm;
830
831 fY1[iknot] = results(4*iknot+1)+2.*results(4*iknot+2)*dxm+
832 3*results(4*iknot+3)*dxm*dxm;
833 }
834 Int_t iknot2= nknots-1;
835 Int_t iknot = nknots-2;
836 Double_t dxm = (fX[iknot2]-fX[iknot2-1])*0.5;
837
838 fY0[iknot2] = results(4*iknot)+ results(4*iknot+1)*dxm+results(4*iknot+2)*dxm*dxm+
839 results(4*iknot+3)*dxm*dxm*dxm;
840
841 fY1[iknot2] = results(4*iknot+1)+2.*results(4*iknot+2)*dxm+
842 3*results(4*iknot+3)*dxm*dxm;
843
844 delete pmatrix;
845 delete pvalues;
846
847}
848
849
850
851
852
853void AliSplineFit::MakeKnots0(TGraph * graph, Double_t maxdelta, Int_t minpoints){
854 //
855 // make knots - restriction max distance and minimum points
856 //
857
858 Int_t npoints = graph->GetN();
859 Double_t *xknots = new Double_t[npoints];
860 Int_t nknots =0;
861 Int_t ipoints =0;
862 //
863 // generate knots
864 //
865 for (Int_t ip=0;ip<npoints;ip++){
866 if (graph->GetX()[ip]-xknots[nknots-1]>maxdelta && ipoints>minpoints){
867 xknots[nknots] = graph->GetX()[ip];
868 ipoints=1;
869 nknots++;
870 }
871 ipoints++;
872 }
873 if (npoints-ipoints>minpoints){
874 xknots[nknots] = graph->GetX()[npoints-1];
875 nknots++;
876 }else{
877 xknots[nknots-1] = graph->GetX()[npoints-1];
878 }
879
880 fN = nknots;
881 fX = new Double_t[nknots];
882 fY0 = new Double_t[nknots];
883 fY1 = new Double_t[nknots];
884 fChi2I= new Double_t[nknots];
885 for (Int_t i=0; i<nknots; i++) fX[i]= xknots[i];
886 delete [] xknots;
887}
888
889
890
891
892void AliSplineFit::MakeSmooth(TGraph * graph, Float_t ratio, char * type){
893 //
894 // Interface to GraphSmooth
895 //
896
897 TGraphSmooth smooth;
898 Int_t npoints2 = TMath::Nint(graph->GetN()*ratio);
899 TGraph * graphT0 = smooth.SmoothKern(graph,type,ratio);
900 if (!graphT0) return;
901 TGraph graphT1(npoints2);
902 for (Int_t ipoint=0; ipoint<npoints2; ipoint++){
903 Int_t pointS = TMath::Nint(ipoint/ratio);
904 if (ipoint==npoints2-1) pointS=graph->GetN()-1;
905 graphT1.SetPoint(ipoint, graphT0->GetX()[pointS] , graphT0->GetY()[pointS]);
906 }
907 TSpline3 spline2("spline", &graphT1);
908 Update(&spline2, npoints2);
909}
910
911
912void AliSplineFit::Update(TSpline3 *spline, Int_t nknots){
913 //
914 //
915 //
916
917 fN = nknots;
918 fX = new Double_t[nknots];
919 fY0 = new Double_t[nknots];
920 fY1 = new Double_t[nknots];
921 Double_t d0, d1;
922 fChi2I= 0;
923 for (Int_t i=0; i<nknots; i++) {
924 spline->GetCoeff(i,fX[i],fY0[i], fY1[i],d0,d1);
925 }
926}
927
928
929
930
931void AliSplineFit::Test(Int_t npoints, Int_t ntracks, Float_t snoise){
932 //
933 // test function
934 //
935
936 AliSplineFit fit;
937 AliSplineFit fitS;
938 TGraph * graph0=0;
939 TGraph * graph1=0;
940
941 TTreeSRedirector *pcstream = new TTreeSRedirector("TestSmooth.root");
942 for (Int_t i=0; i<ntracks; i++){
943 graph0 = AliSplineFit::GenerGraph(npoints,0.05,0,0,1,0);
944 graph1 = AliSplineFit::GenerNoise(graph0,snoise);
945 fit.InitKnots(graph1, 10,10, 0.00);
946 TGraph *d0 = fit.MakeDiff(graph0);
947 TGraph *g0 = fit.MakeGraph(0,1,1000,0);
948 fit.SplineFit(2);
949 TH1F * h2 = fit.MakeDiffHisto(graph0);
950 TGraph *d2 = fit.MakeDiff(graph0);
951 TGraph *g2 = fit.MakeGraph(0,1,1000,0);
952 fit.SplineFit(1);
953 TH1F * h1 = fit.MakeDiffHisto(graph0);
954 TGraph *d1 = fit.MakeDiff(graph0);
955 TGraph *g1 = fit.MakeGraph(0,1,1000,0);
956
957 Float_t ratio = Float_t(fit.fN)/Float_t(npoints);
958 fitS.MakeSmooth(graph1,ratio,"box");
959 TGraph *dS = fitS.MakeDiff(graph0);
960 TGraph *gS = fit.MakeGraph(0,1,1000,0);
961
962 TH1F * hS = fitS.MakeDiffHisto(graph0);
963 Double_t mean2 = h2->GetMean();
964 Double_t sigma2 = h2->GetRMS();
965 Double_t mean1 = h1->GetMean();
966 Double_t sigma1 = h1->GetRMS();
967 Double_t meanS = hS->GetMean();
968 Double_t sigmaS = hS->GetRMS();
969 char fname[100];
970 if (fit.fN<20){
971 sprintf(fname,"pol%d",fit.fN);
972 }else{
973 sprintf(fname,"pol%d",19);
974 }
975 TF1 fpol("fpol",fname);
976 graph1->Fit(&fpol);
977 TGraph dpol(*graph1);
978 TGraph gpol(*graph1);
979 for (Int_t ipoint=0; ipoint<graph1->GetN(); ipoint++){
980 dpol.GetY()[ipoint]= graph0->GetY()[ipoint]-
981 fpol.Eval(graph0->GetX()[ipoint]);
982 gpol.GetY()[ipoint]= fpol.Eval(graph0->GetX()[ipoint]);
983 }
984 (*pcstream)<<"Test"<<
985 "Event="<<i<<
986 "Graph0.="<<graph0<<
987 "Graph1.="<<graph1<<
988 "G0.="<<g0<<
989 "G1.="<<g1<<
990 "G2.="<<g2<<
991 "GS.="<<gS<<
992 "GP.="<<&gpol<<
993 "D0.="<<d0<<
994 "D1.="<<d1<<
995 "D2.="<<d2<<
996 "DS.="<<dS<<
997 "DP.="<<&dpol<<
998 "Npoints="<<fit.fN<<
999 "Mean1="<<mean1<<
1000 "Mean2="<<mean2<<
1001 "MeanS="<<meanS<<
1002 "Sigma1="<<sigma1<<
1003 "Sigma2="<<sigma2<<
1004 "SigmaS="<<sigmaS<<
1005 "\n";
1006
1007 delete graph0;
1008 delete graph1;
1009 delete g1;
1010 delete g2;
1011 delete gS;
1012 delete h1;
1013 delete h2;
1014 delete hS;
1015 }
52fdcd41 1016 delete pcstream;
0dd3a2ac 1017}
8625679a 1018
1019void AliSplineFit::Cleanup(){
1020 //
1021 // deletes extra information to reduce amount of information stored on the data
1022 // base
1023
1024 delete fGraph; fGraph=0;
1025 delete fParams; fParams=0;
1026 delete fCovars; fCovars=0;
1027 delete [] fIndex; fIndex=0;
1028 delete [] fChi2I; fChi2I=0;
1029}
52fdcd41 1030
1031
1032void AliSplineFit::CopyGraph() {
1033 //
1034 // enter graph points directly to fit parameters
1035 // (to bee used when too few points are available)
1036 //
1037 fN = fGraph->GetN();
1038 fX = new Double_t[fN];
1039 fY0 = new Double_t[fN];
1040 fY1 = new Double_t[fN];
1041 for (Int_t i=0; i<fN; i++ ) {
1042 fX[i] = fGraph->GetX()[i];
1043 fY0[i] = fGraph->GetY()[i];
1044 fY1[i] = 0;
1045 }
1046}