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