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