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