[u/mrichter/AliRoot.git] / STEER / AliTrackFitterRieman.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */

///////////////////////////////////////////////////////////////////////////////
//
// Class to the track points on the Riemann sphere. Inputs are
// the set of id's (volids) of the volumes in which residuals are
// calculated to construct a chi2 function to be minimized during 
// the alignment procedures. For the moment the track extrapolation is 
// taken at the space-point reference plane. The reference plane is
// found using the covariance matrix of the point
// (assuming sigma(x)=0 at the reference coordinate system.
//
// Internal usage of AliRieman class for minimization
//
//////////////////////////////////////////////////////////////////////////////

#include <TArrayI.h>
#include <TLinearFitter.h>
#include <TMath.h>
#include <TMatrixD.h>
#include <TMatrixDSym.h>
#include <TRandom.h>
#include <TTreeStream.h>

#include "AliLog.h"
#include "AliLog.h"
#include "AliRieman.h"
#include "AliTrackFitterRieman.h"

ClassImp(AliTrackFitterRieman)


AliTrackFitterRieman::AliTrackFitterRieman():
  AliTrackFitter(),
  fBCorrection(kFALSE),
  fAlpha(0.),
  fNUsed(0),
  fConv(kFALSE),
  fMaxDelta(3),
  fRieman(new AliRieman(10000)),  // allocate rieman
  fMinPointRadius(2.),
  fMaxPointRadius(500.),
  fDebugStream(new TTreeSRedirector("RiemanAlignDebug.root"))
{
  //
  // default constructor
  //
  fCorrY[0] = fCorrY[1] = fCorrY[2] = fCorrY[3] = 0;
  fCorrZ[0] = fCorrZ[1] = fCorrZ[2] = fCorrZ[3] = 0;
}


AliTrackFitterRieman::AliTrackFitterRieman(AliTrackPointArray *array, Bool_t owner):
  AliTrackFitter(array,owner),
  fBCorrection(kFALSE),
  fAlpha(0.),
  fNUsed(0),
  fConv(kFALSE),
  fMaxDelta(3),
  fRieman(new AliRieman(10000)),  //allocate rieman
  fMinPointRadius(2.),
  fMaxPointRadius(500.),
  fDebugStream(0)
{
  //
  // Constructor
  //
  fCorrY[0] = fCorrY[1] = fCorrY[2] = fCorrY[3] = 0;
  fCorrZ[0] = fCorrZ[1] = fCorrZ[2] = fCorrZ[3] = 0;
  if (AliLog::GetDebugLevel("","AliTrackFitterRieman")) fDebugStream = new TTreeSRedirector("RiemanAlignDebug.root");
}

AliTrackFitterRieman::AliTrackFitterRieman(const AliTrackFitterRieman &rieman):
  AliTrackFitter(rieman),
  fBCorrection(rieman.fBCorrection),
  fAlpha(rieman.fAlpha),
  fNUsed(rieman.fNUsed),
  fConv(rieman.fConv),
  fMaxDelta(rieman.fMaxDelta),
  fRieman(new AliRieman(*(rieman.fRieman))),
  fMinPointRadius(rieman.fMinPointRadius),
  fMaxPointRadius(rieman.fMaxPointRadius),
  fDebugStream(0)
{
  //
  // copy constructor
  //
  fCorrY[0] = fCorrY[1] = fCorrY[2] = fCorrY[3] = 0;
  fCorrZ[0] = fCorrZ[1] = fCorrZ[2] = fCorrZ[3] = 0;
  if (AliLog::GetDebugLevel("","AliTrackFitterRieman")) fDebugStream = new TTreeSRedirector("RiemanAlignDebug.root");
}

//_____________________________________________________________________________
AliTrackFitterRieman &AliTrackFitterRieman::operator =(const AliTrackFitterRieman& rieman)
{
  //
  // Assignment operator
  //
  if(this==&rieman) return *this;
  ((AliTrackFitter *)this)->operator=(rieman);

  fBCorrection = rieman.fBCorrection;
  fAlpha  = rieman.fAlpha;
  fNUsed  = rieman.fNUsed;
  fConv   = rieman.fConv;
  fMaxDelta = rieman.fMaxDelta;
  fRieman = new AliRieman(*(rieman.fRieman));
  fMinPointRadius = rieman.fMinPointRadius;
  fMaxPointRadius = rieman.fMaxPointRadius;
  fDebugStream = 0;
  if (AliLog::GetDebugLevel("","AliTrackFitterRieman")) fDebugStream = new TTreeSRedirector("RiemanAlignDebug.root");
  return *this;
}


AliTrackFitterRieman::~AliTrackFitterRieman(){
  //
  //
  //
  delete fRieman;
  delete fDebugStream;
}

void AliTrackFitterRieman::Reset()
{
  // Reset the track parameters and
  // rieman sums
  //
  AliTrackFitter::Reset();
  fRieman->Reset();
  fAlpha = 0.;
  fNUsed = 0;
  fConv =kFALSE;
}

Bool_t AliTrackFitterRieman::Fit(const TArrayI *volIds,const TArrayI *volIdsFit,
				 AliGeomManager::ELayerID layerRangeMin,
			 AliGeomManager::ELayerID layerRangeMax)
{
  // Fit the track points. The method takes as an input
  // the set of id's (volids) of the volumes in which
  // one wants to calculate the residuals.
  // The following parameters are used to define the
  // range of volumes to be used in the fitting
  // As a result two AliTrackPointArray's obects are filled.
  // The first one contains the space points with
  // volume id's from volids list. The second array of points represents
  // the track extrapolations corresponding to the space points
  // in the first array. The two arrays can be used to find
  // the residuals in the volids and consequently construct a
  // chi2 function to be minimized during the alignment
  // procedures. For the moment the track extrapolation is taken
  // at the space-point reference plane. The reference plane is
  // found using the covariance matrix of the point
  // (assuming sigma(x)=0 at the reference coordinate system.
  Int_t debugLevel = AliLog::GetDebugLevel("","AliTrackFitterRieman");
  
  //  Float_t debugRatio = 1./(1.+debugLevel);
  Float_t debugRatio = debugLevel? 1.0/debugLevel : 1.0;

  if (!fPoints) {
    AliError("Track points array not available! Exiting...");
    return kFALSE;
  }
  Int_t npoints = fPoints->GetNPoints();
  if ( npoints<fMinNPoints) return kFALSE;
  //
  // fast count points
  if (volIdsFit != 0x0) {
    Int_t countFit = 0;
    Int_t countPoint = 0;
    for (Int_t ipoint = 0; ipoint < npoints; ipoint++) {
      if (FindVolId(volIds,fPoints->GetVolumeID()[ipoint]))
	countPoint++;
      if (volIdsFit != 0x0) {
	if (FindVolId(volIdsFit,fPoints->GetVolumeID()[ipoint]))
	  countFit++;
      }
    }
    if (countPoint==0) return kFALSE;
    if ((countFit<fMinNPoints) && (volIdsFit != 0x0)) return kFALSE;
  }
  //

  Reset();

  if (fPoints && volIdsFit && AliLog::GetDebugLevel("","AliTrackFitterRieman")>1&& gRandom->Rndm()<debugRatio){
    Int_t nVol    = volIds->GetSize();
    Int_t nVolFit = volIdsFit->GetSize();
    Int_t volId  = volIds->At(0);    
    (*fDebugStream)<<"PInput"<<
      "NPoints="<<npoints<<    // number of points
      "VolId="<<volId<<        // first vol ID
      "NVol="<<nVol<<          // number of volumes 
      "NvolFit="<<nVolFit<<    // number of volumes to fit
      "fPoints.="<<fPoints<<   // input points
      "\n";
  }

  Bool_t isAlphaCalc = kFALSE;
  AliTrackPoint p,plocal;
//   fPoints->GetPoint(p,0);
//   fAlpha = TMath::ATan2(p.GetY(),p.GetX());

  Int_t npVolId = 0;
  fNUsed = 0;
  Int_t *pindex = new Int_t[npoints];
  for (Int_t ipoint = 0; ipoint < npoints; ipoint++)
    {
      fPoints->GetPoint(p,ipoint);
      UShort_t iVolId = p.GetVolumeID();
      if (FindVolId(volIds,iVolId)) {
	pindex[npVolId] = ipoint;
	npVolId++;
      }
      if (volIdsFit != 0x0) {
	if (!FindVolId(volIdsFit,iVolId)) continue;
      }
      else {
	if (iVolId < AliGeomManager::LayerToVolUID(layerRangeMin,0) ||
	    iVolId > AliGeomManager::LayerToVolUID(layerRangeMax,
						AliGeomManager::LayerSize(layerRangeMax))) continue;
      }
      if (!isAlphaCalc) {
	fAlpha = p.GetAngle();
	isAlphaCalc = kTRUE;
      }
      plocal = p.Rotate(fAlpha);
      if (TMath::Abs(plocal.GetX())>fMaxPointRadius || TMath::Abs(plocal.GetX())<fMinPointRadius || plocal.GetCov()[3]<=0 ||plocal.GetCov()[5]<=0 ){
	printf("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<</n");
	p.Dump();
	plocal.Dump();
	printf("Problematic point\n");	
	printf("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<</n");
      }else{
	AddPoint(plocal.GetX(),plocal.GetY(),plocal.GetZ(),
		 TMath::Sqrt(plocal.GetCov()[3]),TMath::Sqrt(plocal.GetCov()[5]));
      }
      // fNUsed++;  AddPoint should be responsible
    }

  if (npVolId == 0 || fNUsed < fMinNPoints) {
    delete [] pindex;
    return kFALSE;
  }

  Update();
 
  if (!fConv) {
    delete [] pindex;
    return kFALSE;
  }

  if ((fParams[0] == 0) ||
      ((-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1) <= 0)) {
    delete [] pindex;
    return kFALSE;
  }


  if (fNUsed < fMinNPoints) {
    delete [] pindex;
    return kFALSE;
  }

  fPVolId = new AliTrackPointArray(npVolId);
  fPTrack = new AliTrackPointArray(npVolId);
  AliTrackPoint p2;
  for (Int_t ipoint = 0; ipoint < npVolId; ipoint++)
    {
      Int_t index = pindex[ipoint];
      fPoints->GetPoint(p,index);
      if (GetPCA(p,p2) && (
	  TMath::Abs(p.GetX()-p2.GetX())<fMaxDelta && 
	  TMath::Abs(p.GetY()-p2.GetY())<fMaxDelta &&
	  TMath::Abs(p.GetZ()-p2.GetZ())<fMaxDelta
	  )) {
	Float_t xyz[3],xyz2[3];
	p.GetXYZ(xyz); p2.GetXYZ(xyz2);
	//	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]);	
	fPVolId->AddPoint(ipoint,&p);
	fPTrack->AddPoint(ipoint,&p2);
      }else{
	// what should be default bahavior -  
	delete [] pindex;
	delete fPVolId;
	delete fPTrack;
	fPVolId =0;
	fPTrack =0;
	return kFALSE;
      }
    }  
  
  if (AliLog::GetDebugLevel("","AliTrackFitterRieman")>0 && gRandom->Rndm()<debugRatio){
    //
    // Debug Info
    //
    AliTrackPointArray *lPVolId = new AliTrackPointArray(npVolId);
    AliTrackPointArray *lPTrack = new AliTrackPointArray(npVolId);
    AliTrackPointArray *lPTrackE = new AliTrackPointArray(npVolId);
    AliRieman * residual = fRieman->MakeResiduals();
    for (Int_t ipoint = 0; ipoint < npVolId; ipoint++){
      AliTrackPoint p0, p0local;      
      AliTrackPoint pFit, pFitlocal, pFitLocalE;      
      fPVolId->GetPoint(p0,ipoint);
      Float_t lAngle = p0.GetAngle();
      p0local= p0.MasterToLocal();
      fPTrack->GetPoint(pFit,ipoint);
      pFitlocal= pFit.Rotate(lAngle);
      //
      Float_t xyz[3], cov[6];
      xyz[0] = pFitlocal.GetX();
      xyz[1] = pFitlocal.GetY();
      xyz[2] = pFitlocal.GetZ();
      for (Int_t icov=0; icov<6; icov++) cov[icov]=0;
      cov[3] = GetErrY2at(xyz[0]);
      cov[5] = GetErrZ2at(xyz[0]);      
      pFitLocalE.SetXYZ(xyz,cov);
      //
      lPVolId->AddPoint(ipoint,&p0local);
      lPTrack->AddPoint(ipoint,&pFitlocal);
      lPTrackE->AddPoint(ipoint,&pFitLocalE);
    }      
    //
    // debug info
    //
    Int_t nVol    = volIds->GetSize();
    Int_t nVolFit = volIdsFit->GetSize();
    Int_t volId   = volIds->At(0);   
    Int_t modId   =0;
    Int_t layer   =  AliGeomManager::VolUIDToLayer(volId,modId);
    Int_t volIdFit   = volIdsFit->At(0);   
    Int_t modIdFit   =0;
    Int_t layerFit   =  AliGeomManager::VolUIDToLayer(volIdFit,modIdFit);
    
    (*fDebugStream)<<"Fit"<<
      "VolId="<<volId<<        // volume ID
      "Layer="<<layer<<        // layer ID
      "Module="<<modId<<       // module ID
      "LayerFit="<<layerFit<<        // layer ID fit
      "ModuleFit="<<modIdFit<<       // module ID fit
      "NVol="<<nVol<<          // number of volumes 
      "NvolFit="<<nVolFit<<    // number of volumes to fit
      "Points0.="<<fPVolId<<    // original points
      "Points1.="<<fPTrack<<    // fitted points 
      "LPoints0.="<<lPVolId<<   // original points - local frame
      "LPoints1.="<<lPTrack<<   // fitted points   - local frame      
      "LPointsE.="<<lPTrackE<<   // fitted points with ext error - local frame      
      "Rieman.="<<this<<        // original rieman fit
      "Res.="<<residual<<       // residuals of rieman fit  
      "\n";
    delete lPVolId;
    delete lPTrack;
    delete residual;
  }

  delete [] pindex;
  return kTRUE;
}


void AliTrackFitterRieman::AddPoint(Float_t x, Float_t y, Float_t z, Float_t sy, Float_t sz)
{
  //
  // add point to rieman fitter
  //
  fRieman->AddPoint(x,y,z,sy,sz);
  fNUsed = fRieman->GetN();
}


Bool_t AliTrackFitterRieman::Update(){
  //
  // 
  //
  fRieman->Update();
  fConv = kFALSE;
  if (fRieman->IsValid()){
    for (Int_t ipar=0; ipar<6; ipar++){
      fParams[ipar] = fRieman->GetParam()[ipar];
    }
    fChi2  =  fRieman->GetChi2();
    fNdf   =  fRieman->GetN()- 2;
    fNUsed =  fRieman->GetN();
    fConv = kTRUE;
  }
  //
  //
  //
  TLinearFitter fitY(3,"pol2");
  TLinearFitter fitZ(3,"pol2");
  for (Int_t ip=0; ip<fRieman->GetN();ip++){
    Double_t x = fRieman->GetX()[ip]; 
    fitY.AddPoint(&x,fRieman->GetY()[ip]-fRieman->GetYat(x),1);
    fitZ.AddPoint(&x,fRieman->GetZ()[ip]-fRieman->GetZat(x),1);    
  }
  fitY.Eval();
  fitZ.Eval();
  for (Int_t iparam=0; iparam<3; iparam++){
    fCorrY[iparam]=fitY.GetParameter(iparam);
    fCorrZ[iparam]=fitZ.GetParameter(iparam);
  }
  fCorrY[3]=fitY.GetChisquare()/Float_t(fRieman->GetN()-3);   
  fCorrZ[3]=fitZ.GetChisquare()/Float_t(fRieman->GetN()-3);

  return kTRUE;
}


//_____________________________________________________________________________
Bool_t AliTrackFitterRieman::GetPCA(const AliTrackPoint &p, AliTrackPoint &p2) const
{
  //
  // Get the closest to a given spacepoint track trajectory point
  // Look for details in the description of the Fit() method
  //
  if (!fConv) return kFALSE;

  // First X and Y coordinates
  Double_t sin = TMath::Sin(fAlpha);
  Double_t cos = TMath::Cos(fAlpha);
  //  fParam[0]  = 1/y0
  //  fParam[1]  = -x0/y0
  //  fParam[2]  = - (R^2 - x0^2 - y0^2)/y0
  if (fParams[0] == 0) return kFALSE;
  // Track parameters in the global coordinate system
  Double_t x0 = -fParams[1]/fParams[0]*cos -         1./fParams[0]*sin;
  Double_t y0 =          1./fParams[0]*cos - fParams[1]/fParams[0]*sin;
  if ((-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1) <= 0) return kFALSE;
  Double_t r  = TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1)/
                fParams[0];

  // Define space-point refence plane
  Double_t alphap = p.GetAngle();
  Double_t sinp = TMath::Sin(alphap);
  Double_t cosp = TMath::Cos(alphap);
  Double_t x  = p.GetX()*cosp + p.GetY()*sinp;
  Double_t y  = p.GetY()*cosp - p.GetX()*sinp;
  Double_t x0p= x0*cosp + y0*sinp;
  Double_t y0p= y0*cosp - x0*sinp;
  if ((r*r - (x-x0p)*(x-x0p))<0) {
    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));
    return kFALSE;
  }
  Double_t temp = TMath::Sqrt(r*r - (x-x0p)*(x-x0p));
  Double_t y1 = y0p + temp;
  Double_t y2 = y0p - temp;
  Double_t yprime = y1;
  if(TMath::Abs(y2-y) < TMath::Abs(y1-y)) yprime = y2;

  // Back to the global coordinate system
  Double_t xsecond = x*cosp - yprime*sinp;
  Double_t ysecond = yprime*cosp + x*sinp;

  // Now Z coordinate and track angles
  Double_t x2 = xsecond*cos + ysecond*sin;
  Double_t zsecond = GetZat(x2);
  Double_t dydx = GetDYat(x2);
  Double_t dzdx = GetDZat(x2);

  // Fill the cov matrix of the track extrapolation point
  Double_t cov[6] = {0,0,0,0,0,0};
  Double_t sigmax = 100*100.;
  cov[0] = sigmax;           cov[1] = sigmax*dydx;      cov[2] = sigmax*dzdx;
  cov[3] = sigmax*dydx*dydx; cov[4] = sigmax*dydx*dzdx;
  cov[5] = sigmax*dzdx*dzdx;

  Double_t sigmay2 = GetErrY2at(x2);
  Double_t sigmaz2 = GetErrZ2at(x2);
  cov[3] += sigmay2;
  cov[5] += sigmaz2;
   

  Float_t  newcov[6];
  newcov[0] = cov[0]*cos*cos-
            2*cov[1]*sin*cos+
              cov[3]*sin*sin;
  newcov[1] = cov[1]*(cos*cos-sin*sin)-
             (cov[3]-cov[0])*sin*cos;
  newcov[2] = cov[2]*cos-
              cov[4]*sin;
  newcov[3] = cov[0]*sin*sin+
            2*cov[1]*sin*cos+
              cov[3]*cos*cos;
  newcov[4] = cov[4]*cos+
              cov[2]*sin;
  newcov[5] = cov[5];

  p2.SetXYZ(xsecond,ysecond,zsecond,newcov);
  Int_t debugLevel = AliLog::GetDebugLevel("","AliTrackFitterRieman");
  Float_t debugRatio = 1./(1.+debugLevel);
  if (AliLog::GetDebugLevel("","AliTrackFitterRieman")>0 && gRandom->Rndm()<debugRatio){
    AliTrackPoint lp0(p);
    AliTrackPoint lp2(p2);
    AliTrackPoint localp0(p);
    AliTrackPoint localp2(p2);
    Float_t lAngle = lp0.GetAngle();
    localp0 = localp0.Rotate(lAngle);
    localp2 = localp2.Rotate(lAngle);

    (*fDebugStream)<<"PCA"<<
      "P0.="<<&lp0<<  //global position
      "P2.="<<&lp2<<
      "LP0.="<<&localp0<<  //local position
      "LP2.="<<&localp2<<
      "\n";
  }
  return kTRUE;
}

Double_t AliTrackFitterRieman::GetYat(Double_t x) const  {
  //
  // get y position at given point 
  //
  Double_t correction=0;
  if (fBCorrection){ // systematic effect correction
    correction  =  fCorrY[0]+fCorrY[1]*x +fCorrY[2]*x*x;
  }
  return fRieman->GetYat(x)+correction;
}

Double_t AliTrackFitterRieman::GetZat(Double_t x) const  {
  //
  // get z position at given point 
  //
  Double_t correction=0;
  if (fBCorrection){ // systematic effect correction
    correction  =  fCorrZ[0]+fCorrZ[1]*x +fCorrZ[2]*x*x;
  }
  return fRieman->GetZat(x)+correction;
}

Double_t AliTrackFitterRieman::GetErrY2at(Double_t x) const {
  //
  // get estimate of extrapolation error
  //
  Double_t error = fRieman->GetErrY(x);
  Double_t correction=0;
  if (fBCorrection){ // everestimate error due systematic effect
    error      *=fCorrY[3];
    correction  =  fCorrY[0]+fCorrY[1]*x +fCorrY[2]*x*x;
    correction *=correction;
  }
  return TMath::Sqrt(error+correction);
}

Double_t AliTrackFitterRieman::GetErrZ2at(Double_t x) const {
  //
  // get estimate of extrapolation error
  //
  Double_t error = fRieman->GetErrZ(x)*fCorrZ[3];
  Double_t correction=0;
  if (fBCorrection){  
     error    *=  fCorrZ[3];
    correction =  fCorrZ[0]+fCorrZ[1]*x +fCorrZ[2]*x*x;
    correction*=  correction;
  }
  return TMath::Sqrt(error+correction);
}

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