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