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