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