1f69321795e5212c6b9b199737532787336b789c
[u/mrichter/AliRoot.git] / STEER / AliTrackFitterRieman.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, 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 /* $Id$ */
17
18 ///////////////////////////////////////////////////////////////////////////////
19 //
20 // Class to the track points on the Riemann sphere. Inputs are
21 // the set of id's (volids) of the volumes in which residuals are
22 // calculated to construct a chi2 function to be minimized during 
23 // the alignment procedures. For the moment the track extrapolation is 
24 // taken at the space-point reference plane. The reference plane is
25 // found using the covariance matrix of the point
26 // (assuming sigma(x)=0 at the reference coordinate system.
27 //
28 // Internal usage of AliRieman class for minimization
29 //
30 //////////////////////////////////////////////////////////////////////////////
31
32 #include <TArrayI.h>
33 #include <TLinearFitter.h>
34 #include <TMath.h>
35 #include <TMatrixD.h>
36 #include <TMatrixDSym.h>
37 #include <TRandom.h>
38 #include <TTreeStream.h>
39
40 #include "AliLog.h"
41 #include "AliLog.h"
42 #include "AliRieman.h"
43 #include "AliTrackFitterRieman.h"
44
45 ClassImp(AliTrackFitterRieman)
46
47
48 AliTrackFitterRieman::AliTrackFitterRieman():
49   AliTrackFitter(),
50   fBCorrection(kFALSE),
51   fAlpha(0.),
52   fNUsed(0),
53   fConv(kFALSE),
54   fMaxDelta(3),
55   fRieman(new AliRieman(10000)),  // allocate rieman
56   fMinPointRadius(2.),
57   fMaxPointRadius(500.),
58   fDebugStream(new TTreeSRedirector("RiemanAlignDebug.root"))
59 {
60   //
61   // default constructor
62   //
63   fCorrY[0] = fCorrY[1] = fCorrY[2] = fCorrY[3] = 0;
64   fCorrZ[0] = fCorrZ[1] = fCorrZ[2] = fCorrZ[3] = 0;
65 }
66
67
68 AliTrackFitterRieman::AliTrackFitterRieman(AliTrackPointArray *array, Bool_t owner):
69   AliTrackFitter(array,owner),
70   fBCorrection(kFALSE),
71   fAlpha(0.),
72   fNUsed(0),
73   fConv(kFALSE),
74   fMaxDelta(3),
75   fRieman(new AliRieman(10000)),  //allocate rieman
76   fMinPointRadius(2.),
77   fMaxPointRadius(500.),
78   fDebugStream(0)
79 {
80   //
81   // Constructor
82   //
83   fCorrY[0] = fCorrY[1] = fCorrY[2] = fCorrY[3] = 0;
84   fCorrZ[0] = fCorrZ[1] = fCorrZ[2] = fCorrZ[3] = 0;
85   if (AliLog::GetDebugLevel("","AliTrackFitterRieman")) fDebugStream = new TTreeSRedirector("RiemanAlignDebug.root");
86 }
87
88 AliTrackFitterRieman::AliTrackFitterRieman(const AliTrackFitterRieman &rieman):
89   AliTrackFitter(rieman),
90   fBCorrection(rieman.fBCorrection),
91   fAlpha(rieman.fAlpha),
92   fNUsed(rieman.fNUsed),
93   fConv(rieman.fConv),
94   fMaxDelta(rieman.fMaxDelta),
95   fRieman(new AliRieman(*(rieman.fRieman))),
96   fMinPointRadius(rieman.fMinPointRadius),
97   fMaxPointRadius(rieman.fMaxPointRadius),
98   fDebugStream(0)
99 {
100   //
101   // copy constructor
102   //
103   fCorrY[0] = fCorrY[1] = fCorrY[2] = fCorrY[3] = 0;
104   fCorrZ[0] = fCorrZ[1] = fCorrZ[2] = fCorrZ[3] = 0;
105   if (AliLog::GetDebugLevel("","AliTrackFitterRieman")) fDebugStream = new TTreeSRedirector("RiemanAlignDebug.root");
106 }
107
108 //_____________________________________________________________________________
109 AliTrackFitterRieman &AliTrackFitterRieman::operator =(const AliTrackFitterRieman& rieman)
110 {
111   //
112   // Assignment operator
113   //
114   if(this==&rieman) return *this;
115   ((AliTrackFitter *)this)->operator=(rieman);
116
117   fBCorrection = rieman.fBCorrection;
118   fAlpha  = rieman.fAlpha;
119   fNUsed  = rieman.fNUsed;
120   fConv   = rieman.fConv;
121   fMaxDelta = rieman.fMaxDelta;
122   fRieman = new AliRieman(*(rieman.fRieman));
123   fMinPointRadius = rieman.fMinPointRadius;
124   fMaxPointRadius = rieman.fMaxPointRadius;
125   fDebugStream = 0;
126   if (AliLog::GetDebugLevel("","AliTrackFitterRieman")) fDebugStream = new TTreeSRedirector("RiemanAlignDebug.root");
127   return *this;
128 }
129
130
131 AliTrackFitterRieman::~AliTrackFitterRieman(){
132   //
133   //
134   //
135   delete fRieman;
136   delete fDebugStream;
137 }
138
139 void AliTrackFitterRieman::Reset()
140 {
141   // Reset the track parameters and
142   // rieman sums
143   //
144   AliTrackFitter::Reset();
145   fRieman->Reset();
146   fAlpha = 0.;
147   fNUsed = 0;
148   fConv =kFALSE;
149 }
150
151 Bool_t AliTrackFitterRieman::Fit(const TArrayI *volIds,const TArrayI *volIdsFit,
152                                  AliGeomManager::ELayerID layerRangeMin,
153                          AliGeomManager::ELayerID layerRangeMax)
154 {
155   // Fit the track points. The method takes as an input
156   // the set of id's (volids) of the volumes in which
157   // one wants to calculate the residuals.
158   // The following parameters are used to define the
159   // range of volumes to be used in the fitting
160   // As a result two AliTrackPointArray's obects are filled.
161   // The first one contains the space points with
162   // volume id's from volids list. The second array of points represents
163   // the track extrapolations corresponding to the space points
164   // in the first array. The two arrays can be used to find
165   // the residuals in the volids and consequently construct a
166   // chi2 function to be minimized during the alignment
167   // procedures. For the moment the track extrapolation is taken
168   // at the space-point reference plane. The reference plane is
169   // found using the covariance matrix of the point
170   // (assuming sigma(x)=0 at the reference coordinate system.
171   Int_t debugLevel = AliLog::GetDebugLevel("","AliTrackFitterRieman");
172   
173   //  Float_t debugRatio = 1./(1.+debugLevel);
174   Float_t debugRatio = debugLevel? 1.0/debugLevel : 1.0;
175
176   Int_t npoints = fPoints->GetNPoints();
177   if ( npoints<fMinNPoints) return kFALSE;
178   //
179   // fast count points
180   if (volIdsFit != 0x0) {
181     Int_t countFit = 0;
182     Int_t countPoint = 0;
183     for (Int_t ipoint = 0; ipoint < npoints; ipoint++) {
184       if (FindVolId(volIds,fPoints->GetVolumeID()[ipoint]))
185         countPoint++;
186       if (volIdsFit != 0x0) {
187         if (FindVolId(volIdsFit,fPoints->GetVolumeID()[ipoint]))
188           countFit++;
189       }
190     }
191     if (countPoint==0) return kFALSE;
192     if ((countFit<fMinNPoints) && (volIdsFit != 0x0)) return kFALSE;
193   }
194   //
195
196   Reset();
197
198   if (fPoints && AliLog::GetDebugLevel("","AliTrackFitterRieman")>1&& gRandom->Rndm()<debugRatio){
199     Int_t nVol    = volIds->GetSize();
200     Int_t nVolFit = volIdsFit->GetSize();
201     Int_t volId  = volIds->At(0);    
202     (*fDebugStream)<<"PInput"<<
203       "NPoints="<<npoints<<    // number of points
204       "VolId="<<volId<<        // first vol ID
205       "NVol="<<nVol<<          // number of volumes 
206       "NvolFit="<<nVolFit<<    // number of volumes to fit
207       "fPoints.="<<fPoints<<   // input points
208       "\n";
209   }
210
211   Bool_t isAlphaCalc = kFALSE;
212   AliTrackPoint p,plocal;
213 //   fPoints->GetPoint(p,0);
214 //   fAlpha = TMath::ATan2(p.GetY(),p.GetX());
215
216   Int_t npVolId = 0;
217   fNUsed = 0;
218   Int_t *pindex = new Int_t[npoints];
219   for (Int_t ipoint = 0; ipoint < npoints; ipoint++)
220     {
221       fPoints->GetPoint(p,ipoint);
222       UShort_t iVolId = p.GetVolumeID();
223       if (FindVolId(volIds,iVolId)) {
224         pindex[npVolId] = ipoint;
225         npVolId++;
226       }
227       if (volIdsFit != 0x0) {
228         if (!FindVolId(volIdsFit,iVolId)) continue;
229       }
230       else {
231         if (iVolId < AliGeomManager::LayerToVolUID(layerRangeMin,0) ||
232             iVolId > AliGeomManager::LayerToVolUID(layerRangeMax,
233                                                 AliGeomManager::LayerSize(layerRangeMax))) continue;
234       }
235       if (!isAlphaCalc) {
236         fAlpha = p.GetAngle();
237         isAlphaCalc = kTRUE;
238       }
239       plocal = p.Rotate(fAlpha);
240       if (TMath::Abs(plocal.GetX())>fMaxPointRadius || TMath::Abs(plocal.GetX())<fMinPointRadius || plocal.GetCov()[3]<=0 ||plocal.GetCov()[5]<=0 ){
241         printf("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<</n");
242         p.Dump();
243         plocal.Dump();
244         printf("Problematic point\n");  
245         printf("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<</n");
246       }else{
247         AddPoint(plocal.GetX(),plocal.GetY(),plocal.GetZ(),
248                  TMath::Sqrt(plocal.GetCov()[3]),TMath::Sqrt(plocal.GetCov()[5]));
249       }
250       // fNUsed++;  AddPoint should be responsible
251     }
252
253   if (npVolId == 0 || fNUsed < fMinNPoints) {
254     delete [] pindex;
255     return kFALSE;
256   }
257
258   Update();
259  
260   if (!fConv) {
261     delete [] pindex;
262     return kFALSE;
263   }
264
265   if ((fParams[0] == 0) ||
266       ((-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1) <= 0)) {
267     delete [] pindex;
268     return kFALSE;
269   }
270
271
272   if (fNUsed < fMinNPoints) {
273     delete [] pindex;
274     return kFALSE;
275   }
276
277   fPVolId = new AliTrackPointArray(npVolId);
278   fPTrack = new AliTrackPointArray(npVolId);
279   AliTrackPoint p2;
280   for (Int_t ipoint = 0; ipoint < npVolId; ipoint++)
281     {
282       Int_t index = pindex[ipoint];
283       fPoints->GetPoint(p,index);
284       if (GetPCA(p,p2) && (
285           TMath::Abs(p.GetX()-p2.GetX())<fMaxDelta && 
286           TMath::Abs(p.GetY()-p2.GetY())<fMaxDelta &&
287           TMath::Abs(p.GetZ()-p2.GetZ())<fMaxDelta
288           )) {
289         Float_t xyz[3],xyz2[3];
290         p.GetXYZ(xyz); p2.GetXYZ(xyz2);
291         //      printf("residuals %f %d %d %f %f %f %f %f %f\n",fChi2,fNUsed,fConv,xyz[0],xyz[1],xyz[2],xyz2[0]-xyz[0],xyz2[1]-xyz[1],xyz2[2]-xyz[2]);  
292         fPVolId->AddPoint(ipoint,&p);
293         fPTrack->AddPoint(ipoint,&p2);
294       }else{
295         // what should be default bahavior -  
296         delete [] pindex;
297         delete fPVolId;
298         delete fPTrack;
299         fPVolId =0;
300         fPTrack =0;
301         return kFALSE;
302       }
303     }  
304   
305   if (AliLog::GetDebugLevel("","AliTrackFitterRieman")>0 && gRandom->Rndm()<debugRatio){
306     //
307     // Debug Info
308     //
309     AliTrackPointArray *lPVolId = new AliTrackPointArray(npVolId);
310     AliTrackPointArray *lPTrack = new AliTrackPointArray(npVolId);
311     AliTrackPointArray *lPTrackE = new AliTrackPointArray(npVolId);
312     AliRieman * residual = fRieman->MakeResiduals();
313     for (Int_t ipoint = 0; ipoint < npVolId; ipoint++){
314       AliTrackPoint p0, p0local;      
315       AliTrackPoint pFit, pFitlocal, pFitLocalE;      
316       fPVolId->GetPoint(p0,ipoint);
317       Float_t lAngle = p0.GetAngle();
318       p0local= p0.MasterToLocal();
319       fPTrack->GetPoint(pFit,ipoint);
320       pFitlocal= pFit.Rotate(lAngle);
321       //
322       Float_t xyz[3], cov[6];
323       xyz[0] = pFitlocal.GetX();
324       xyz[1] = pFitlocal.GetY();
325       xyz[2] = pFitlocal.GetZ();
326       for (Int_t icov=0; icov<6; icov++) cov[icov]=0;
327       cov[3] = GetErrY2at(xyz[0]);
328       cov[5] = GetErrZ2at(xyz[0]);      
329       pFitLocalE.SetXYZ(xyz,cov);
330       //
331       lPVolId->AddPoint(ipoint,&p0local);
332       lPTrack->AddPoint(ipoint,&pFitlocal);
333       lPTrackE->AddPoint(ipoint,&pFitLocalE);
334     }      
335     //
336     // debug info
337     //
338     Int_t nVol    = volIds->GetSize();
339     Int_t nVolFit = volIdsFit->GetSize();
340     Int_t volId   = volIds->At(0);   
341     Int_t modId   =0;
342     Int_t layer   =  AliGeomManager::VolUIDToLayer(volId,modId);
343     Int_t volIdFit   = volIdsFit->At(0);   
344     Int_t modIdFit   =0;
345     Int_t layerFit   =  AliGeomManager::VolUIDToLayer(volIdFit,modIdFit);
346     
347     (*fDebugStream)<<"Fit"<<
348       "VolId="<<volId<<        // volume ID
349       "Layer="<<layer<<        // layer ID
350       "Module="<<modId<<       // module ID
351       "LayerFit="<<layerFit<<        // layer ID fit
352       "ModuleFit="<<modIdFit<<       // module ID fit
353       "NVol="<<nVol<<          // number of volumes 
354       "NvolFit="<<nVolFit<<    // number of volumes to fit
355       "Points0.="<<fPVolId<<    // original points
356       "Points1.="<<fPTrack<<    // fitted points 
357       "LPoints0.="<<lPVolId<<   // original points - local frame
358       "LPoints1.="<<lPTrack<<   // fitted points   - local frame      
359       "LPointsE.="<<lPTrackE<<   // fitted points with ext error - local frame      
360       "Rieman.="<<this<<        // original rieman fit
361       "Res.="<<residual<<       // residuals of rieman fit  
362       "\n";
363     delete lPVolId;
364     delete lPTrack;
365     delete residual;
366   }
367
368   delete [] pindex;
369   return kTRUE;
370 }
371
372
373 void AliTrackFitterRieman::AddPoint(Float_t x, Float_t y, Float_t z, Float_t sy, Float_t sz)
374 {
375   //
376   // add point to rieman fitter
377   //
378   fRieman->AddPoint(x,y,z,sy,sz);
379   fNUsed = fRieman->GetN();
380 }
381
382
383
384 Bool_t AliTrackFitterRieman::Update(){
385   //
386   // 
387   //
388   fRieman->Update();
389   fConv = kFALSE;
390   if (fRieman->IsValid()){
391     for (Int_t ipar=0; ipar<6; ipar++){
392       fParams[ipar] = fRieman->GetParam()[ipar];
393     }
394     fChi2  =  fRieman->GetChi2();
395     fNdf   =  fRieman->GetN()- 2;
396     fNUsed =  fRieman->GetN();
397     fConv = kTRUE;
398   }
399   //
400   //
401   //
402   TLinearFitter fitY(3,"pol2");
403   TLinearFitter fitZ(3,"pol2");
404   for (Int_t ip=0; ip<fRieman->GetN();ip++){
405     Double_t x = fRieman->GetX()[ip]; 
406     fitY.AddPoint(&x,fRieman->GetY()[ip]-fRieman->GetYat(x),1);
407     fitZ.AddPoint(&x,fRieman->GetZ()[ip]-fRieman->GetZat(x),1);    
408   }
409   fitY.Eval();
410   fitZ.Eval();
411   for (Int_t iparam=0; iparam<3; iparam++){
412     fCorrY[iparam]=fitY.GetParameter(iparam);
413     fCorrZ[iparam]=fitZ.GetParameter(iparam);
414   }
415   fCorrY[3]=fitY.GetChisquare()/Float_t(fRieman->GetN()-3);   
416   fCorrZ[3]=fitZ.GetChisquare()/Float_t(fRieman->GetN()-3);
417
418   return kTRUE;
419 }
420
421
422
423 //_____________________________________________________________________________
424 Bool_t AliTrackFitterRieman::GetPCA(const AliTrackPoint &p, AliTrackPoint &p2) const
425 {
426   //
427   // Get the closest to a given spacepoint track trajectory point
428   // Look for details in the description of the Fit() method
429   //
430   if (!fConv) return kFALSE;
431
432   // First X and Y coordinates
433   Double_t sin = TMath::Sin(fAlpha);
434   Double_t cos = TMath::Cos(fAlpha);
435   //  fParam[0]  = 1/y0
436   //  fParam[1]  = -x0/y0
437   //  fParam[2]  = - (R^2 - x0^2 - y0^2)/y0
438   if (fParams[0] == 0) return kFALSE;
439   // Track parameters in the global coordinate system
440   Double_t x0 = -fParams[1]/fParams[0]*cos -         1./fParams[0]*sin;
441   Double_t y0 =          1./fParams[0]*cos - fParams[1]/fParams[0]*sin;
442   if ((-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1) <= 0) return kFALSE;
443   Double_t r  = TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1)/
444                 fParams[0];
445
446   // Define space-point refence plane
447   Double_t alphap = p.GetAngle();
448   Double_t sinp = TMath::Sin(alphap);
449   Double_t cosp = TMath::Cos(alphap);
450   Double_t x  = p.GetX()*cosp + p.GetY()*sinp;
451   Double_t y  = p.GetY()*cosp - p.GetX()*sinp;
452   Double_t x0p= x0*cosp + y0*sinp;
453   Double_t y0p= y0*cosp - x0*sinp;
454   if ((r*r - (x-x0p)*(x-x0p))<0) {
455     AliWarning(Form("Track extrapolation failed ! (Track radius = %f, track circle x = %f, space-point x = %f, reference plane angle = %f\n",r,x0p,x,alphap));
456     return kFALSE;
457   }
458   Double_t temp = TMath::Sqrt(r*r - (x-x0p)*(x-x0p));
459   Double_t y1 = y0p + temp;
460   Double_t y2 = y0p - temp;
461   Double_t yprime = y1;
462   if(TMath::Abs(y2-y) < TMath::Abs(y1-y)) yprime = y2;
463
464   // Back to the global coordinate system
465   Double_t xsecond = x*cosp - yprime*sinp;
466   Double_t ysecond = yprime*cosp + x*sinp;
467
468   // Now Z coordinate and track angles
469   Double_t x2 = xsecond*cos + ysecond*sin;
470   Double_t zsecond = GetZat(x2);
471   Double_t dydx = GetDYat(x2);
472   Double_t dzdx = GetDZat(x2);
473
474   // Fill the cov matrix of the track extrapolation point
475   Double_t cov[6] = {0,0,0,0,0,0};
476   Double_t sigmax = 100*100.;
477   cov[0] = sigmax;           cov[1] = sigmax*dydx;      cov[2] = sigmax*dzdx;
478   cov[3] = sigmax*dydx*dydx; cov[4] = sigmax*dydx*dzdx;
479   cov[5] = sigmax*dzdx*dzdx;
480
481   Double_t sigmay2 = GetErrY2at(x2);
482   Double_t sigmaz2 = GetErrZ2at(x2);
483   cov[3] += sigmay2;
484   cov[5] += sigmaz2;
485    
486
487   Float_t  newcov[6];
488   newcov[0] = cov[0]*cos*cos-
489             2*cov[1]*sin*cos+
490               cov[3]*sin*sin;
491   newcov[1] = cov[1]*(cos*cos-sin*sin)-
492              (cov[3]-cov[0])*sin*cos;
493   newcov[2] = cov[2]*cos-
494               cov[4]*sin;
495   newcov[3] = cov[0]*sin*sin+
496             2*cov[1]*sin*cos+
497               cov[3]*cos*cos;
498   newcov[4] = cov[4]*cos+
499               cov[2]*sin;
500   newcov[5] = cov[5];
501
502   p2.SetXYZ(xsecond,ysecond,zsecond,newcov);
503   Int_t debugLevel = AliLog::GetDebugLevel("","AliTrackFitterRieman");
504   Float_t debugRatio = 1./(1.+debugLevel);
505   if (AliLog::GetDebugLevel("","AliTrackFitterRieman")>0 && gRandom->Rndm()<debugRatio){
506     AliTrackPoint lp0(p);
507     AliTrackPoint lp2(p2);
508     AliTrackPoint localp0(p);
509     AliTrackPoint localp2(p2);
510     Float_t lAngle = lp0.GetAngle();
511     localp0 = localp0.Rotate(lAngle);
512     localp2 = localp2.Rotate(lAngle);
513
514     (*fDebugStream)<<"PCA"<<
515       "P0.="<<&lp0<<  //global position
516       "P2.="<<&lp2<<
517       "LP0.="<<&localp0<<  //local position
518       "LP2.="<<&localp2<<
519       "\n";
520   }
521   return kTRUE;
522 }
523
524 Double_t AliTrackFitterRieman::GetYat(Double_t x) const  {
525   //
526   // get y position at given point 
527   //
528   Double_t correction=0;
529   if (fBCorrection){ // systematic effect correction
530     correction  =  fCorrY[0]+fCorrY[1]*x +fCorrY[2]*x*x;
531   }
532   return fRieman->GetYat(x)+correction;
533 }
534
535 Double_t AliTrackFitterRieman::GetZat(Double_t x) const  {
536   //
537   // get z position at given point 
538   //
539   Double_t correction=0;
540   if (fBCorrection){ // systematic effect correction
541     correction  =  fCorrZ[0]+fCorrZ[1]*x +fCorrZ[2]*x*x;
542   }
543   return fRieman->GetZat(x)+correction;
544 }
545
546 Double_t AliTrackFitterRieman::GetErrY2at(Double_t x) const {
547   //
548   // get estimate of extrapolation error
549   //
550   Double_t error = fRieman->GetErrY(x);
551   Double_t correction=0;
552   if (fBCorrection){ // everestimate error due systematic effect
553     error      *=fCorrY[3];
554     correction  =  fCorrY[0]+fCorrY[1]*x +fCorrY[2]*x*x;
555     correction *=correction;
556   }
557   return TMath::Sqrt(error+correction);
558 }
559
560 Double_t AliTrackFitterRieman::GetErrZ2at(Double_t x) const {
561   //
562   // get estimate of extrapolation error
563   //
564   Double_t error = fRieman->GetErrZ(x)*fCorrZ[3];
565   Double_t correction=0;
566   if (fBCorrection){  
567      error    *=  fCorrZ[3];
568     correction =  fCorrZ[0]+fCorrZ[1]*x +fCorrZ[2]*x*x;
569     correction*=  correction;
570   }
571   return TMath::Sqrt(error+correction);
572 }
573
574 void AliTrackFitterRieman::SetParam(Int_t i, Double_t par) {
575   if (i<0 || i>5) return;
576   fParams[i]=par;
577   fRieman->GetParam()[i]=par;
578 }