-#include "TMatrixDSym.h"
-#include "TMatrixD.h"
+/**************************************************************************
+ * 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()
+ 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
//
- fAlpha = 0.;
- for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
- fSumYY = 0;
- for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
- fSumZZ = 0;
- fNUsed = 0;
- fConv = kFALSE;
}
AliTrackFitterRieman::AliTrackFitterRieman(AliTrackPointArray *array, Bool_t owner):
- AliTrackFitter(array,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
//
- fAlpha = 0.;
- for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
- fSumYY = 0;
- for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
- fSumZZ = 0;
- fNUsed = 0;
- fConv = kFALSE;
+ if (AliLog::GetDebugLevel("","AliTrackFitterRieman")) fDebugStream = new TTreeSRedirector("RiemanAlignDebug.root");
}
AliTrackFitterRieman::AliTrackFitterRieman(const AliTrackFitterRieman &rieman):
- AliTrackFitter(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
//
- fAlpha = rieman.fAlpha;
- for (Int_t i=0;i<9;i++) fSumXY[i] = rieman.fSumXY[i];
- fSumYY = rieman.fSumYY;
- for (Int_t i=0;i<9;i++) fSumXZ[i] = rieman.fSumXZ[i];
- fSumZZ = rieman.fSumZZ;
- fNUsed = rieman.fNUsed;
- fConv = rieman.fConv;
+ if (AliLog::GetDebugLevel("","AliTrackFitterRieman")) fDebugStream = new TTreeSRedirector("RiemanAlignDebug.root");
}
//_____________________________________________________________________________
AliTrackFitterRieman &AliTrackFitterRieman::operator =(const AliTrackFitterRieman& rieman)
{
- // assignment operator
+ //
+ // Assignment operator
//
if(this==&rieman) return *this;
((AliTrackFitter *)this)->operator=(rieman);
- fAlpha = rieman.fAlpha;
- for (Int_t i=0;i<9;i++) fSumXY[i] = rieman.fSumXY[i];
- fSumYY = rieman.fSumYY;
- for (Int_t i=0;i<9;i++) fSumXZ[i] = rieman.fSumXZ[i];
- fSumZZ = rieman.fSumZZ;
- fNUsed = rieman.fNUsed;
- fConv = rieman.fConv;
-
+ 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()
-{
- // destructor
+
+AliTrackFitterRieman::~AliTrackFitterRieman(){
+ //
//
+ //
+ delete fRieman;
+ delete fDebugStream;
}
void AliTrackFitterRieman::Reset()
{
// Reset the track parameters and
// rieman sums
+ //
AliTrackFitter::Reset();
+ fRieman->Reset();
fAlpha = 0.;
- for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
- fSumYY = 0;
- for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
- fSumZZ = 0;
fNUsed = 0;
fConv =kFALSE;
}
Bool_t AliTrackFitterRieman::Fit(const TArrayI *volIds,const TArrayI *volIdsFit,
- AliAlignObj::ELayerID layerRangeMin,
- AliAlignObj::ELayerID layerRangeMax)
+ 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
// 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;
+
+ 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();
- Int_t npoints = fPoints->GetNPoints();
- if (npoints < 3) return kFALSE;
+ if (fPoints && 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;
Int_t npVolId = 0;
fNUsed = 0;
Int_t *pindex = new Int_t[npoints];
- fX = new Float_t[npoints];
- fY = new Float_t[npoints];
- fZ = new Float_t[npoints];
- fSy = new Float_t[npoints];
- fSz = new Float_t[npoints];
for (Int_t ipoint = 0; ipoint < npoints; ipoint++)
{
fPoints->GetPoint(p,ipoint);
if (!FindVolId(volIdsFit,iVolId)) continue;
}
else {
- if (iVolId < AliAlignObj::LayerToVolUID(layerRangeMin,0) ||
- iVolId > AliAlignObj::LayerToVolUID(layerRangeMax,
- AliAlignObj::LayerSize(layerRangeMax-
- AliAlignObj::kFirstLayer))) continue;
+ 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);
- AddPoint(plocal.GetX(),plocal.GetY(),plocal.GetZ(),
- TMath::Sqrt(plocal.GetCov()[3]),TMath::Sqrt(plocal.GetCov()[5]));
- fNUsed++;
+ 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 (fNUsed < 3) {
+ if (npVolId == 0 || fNUsed < fMinNPoints) {
delete [] pindex;
- delete [] fX;
- delete [] fY;
- delete [] fZ;
- delete [] fSy;
- delete [] fSz;
return kFALSE;
}
Update();
-
- delete [] fX;
- delete [] fY;
- delete [] fZ;
- delete [] fSy;
- delete [] fSz;
if (!fConv) {
delete [] pindex;
{
Int_t index = pindex[ipoint];
fPoints->GetPoint(p,index);
- if (GetPCA(p,p2)) {
+ 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]);
+ // 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;
-
- // debug info
-// Float_t chi2 = 0, chi22 = 0;
-// for (Int_t ipoint = 0; ipoint < npoints; ipoint++)
-// {
-// fPoints->GetPoint(p,ipoint);
-// UShort_t iVolId = p.GetVolumeID();
-// if (volIdFit != 0) {
-// if (iVolId != volIdFit) continue;
-// }
-// else {
-// if (iVolId < AliAlignObj::LayerToVolUID(layerRangeMin,0) ||
-// iVolId > AliAlignObj::LayerToVolUID(layerRangeMax,AliAlignObj::LayerSize(layerRangeMax-
-// AliAlignObj::kFirstLayer))) continue;
-// }
-// plocal = p.Rotate(fAlpha);
-// Float_t delta = (fParams[0]*(plocal.GetX()*plocal.GetX()+plocal.GetY()*plocal.GetY())+
-// 2.*plocal.GetX()*fParams[1]+
-// fParams[2]-
-// 2.*plocal.GetY())/
-// (2.*TMath::Sqrt(plocal.GetCov()[3]));
-// // Float_t delta2 = (fParams[3]+
-// // plocal.GetX()*fParams[4]+
-// // plocal.GetX()*plocal.GetX()*fParams[5]-
-// // plocal.GetZ())/
-// // (TMath::Sqrt(plocal.GetCov()[5]));
-// Double_t r = TMath::Sqrt(plocal.GetX()*plocal.GetX()+plocal.GetY()*plocal.GetY());
-// Double_t Rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
-// Float_t delta2 = (fParams[3]+
-// r*fParams[4]+r*r*r*fParams[4]*Rm1*Rm1/24-
-// plocal.GetZ())/
-// (TMath::Sqrt(plocal.GetCov()[5]));
-// chi2 += delta*delta;
-// chi22 += delta2*delta2;
-// // printf("pulls %d %d %f %f\n",ipoint,iVolId,delta,delta2);
-
-// }
-// printf("My chi2 = %f + %f = %f\n",chi2,chi22,chi2+chi22);
-
return kTRUE;
}
+
void AliTrackFitterRieman::AddPoint(Float_t x, Float_t y, Float_t z, Float_t sy, Float_t sz)
{
//
- // Rieman update
- //
- //------------------------------------------------------
- // XY direction
+ // add point to rieman fitter
//
- // (x-x0)^2+(y-y0)^2-R^2=0 ===>
- //
- // (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0; ==>
- //
- // substitution t = 1/(x^2+y^2), u = 2*x*t, v = 2*y*t, D0 = R^2 - x0^2- y0^2
- //
- // 1 - u*x0 - v*y0 - t *D0 =0 ; - linear equation
- //
- // next substition a = 1/y0 b = -x0/y0 c = -D0/y0
- //
- // final linear equation : a + u*b +t*c - v =0;
- //
- // Minimization :
- //
- // sum( (a + ui*b +ti*c - vi)^2)/(sigmai)^2 = min;
- //
- // where sigmai is the error of maesurement (a + ui*b +ti*c - vi)
- //
- // neglecting error of xi, and supposing xi>>yi sigmai ~ sigmaVi ~ 2*sigmay*t
- //
- fX[fNUsed] = x; fY[fNUsed]=y; fZ[fNUsed]=z; fSy[fNUsed]=sy; fSz[fNUsed]=sz;
- //
- // XY part
- //
- Double_t t = x*x+y*y;
- if (t<2) return;
- t = 1./t;
- Double_t u = 2.*x*t;
- Double_t v = 2.*y*t;
- Double_t error = 2.*sy*t;
- error *=error;
- Double_t weight = 1./error;
- fSumXY[0] +=weight;
- fSumXY[1] +=u*weight; fSumXY[2]+=v*weight; fSumXY[3]+=t*weight;
- fSumXY[4] +=u*u*weight; fSumXY[5]+=t*t*weight;
- fSumXY[6] +=u*v*weight; fSumXY[7]+=u*t*weight; fSumXY[8]+=v*t*weight;
- fSumYY += v*v*weight;
- //
- // XZ part
- //
- if (1) {
- weight = 1./(sz*sz);
-// fSumXZ[0] +=weight;
-// fSumXZ[1] +=weight*x; fSumXZ[2] +=weight*x*x; fSumXZ[3] +=weight*x*x*x; fSumXZ[4] += weight*x*x*x*x;
-// fSumXZ[5] +=weight*z; fSumXZ[6] +=weight*x*z; fSumXZ[7] +=weight*x*x*z;
- fSumZZ += z*z*weight;
- }
- else {
- weight = 1./(sz*sz);
- fSumXZ[0] +=weight;
- Double_t r = TMath::Sqrt(x*x+y*y);
- fSumXZ[1] +=weight*r; fSumXZ[2] +=weight*r*r; fSumXZ[3] +=weight*z; fSumXZ[4] += weight*r*z;
- // Now the auxulary sums
- fSumXZ[5] +=weight*r*r*r/24; fSumXZ[6] +=weight*r*r*r*r/12; fSumXZ[7] +=weight*r*r*r*z/24;
- fSumXZ[8] +=weight*r*r*r*r*r*r/(24*24);
- fSumZZ += z*z*weight;
- }
+ fRieman->AddPoint(x,y,z,sy,sz);
+ fNUsed = fRieman->GetN();
}
-void AliTrackFitterRieman::Update(){
- //
- // Rieman update
- //
- //
- for (Int_t i=0;i<6;i++)fParams[i]=0;
- fChi2 = 0;
- fNdf = 0;
- Int_t conv=0;
- //
- // XY part
- //
- TMatrixDSym smatrix(3);
- TMatrixD sums(1,3);
- //
- // smatrix(0,0) = s0; smatrix(1,1)=su2; smatrix(2,2)=st2;
- // smatrix(0,1) = su; smatrix(0,2)=st; smatrix(1,2)=sut;
- // sums(0,0) = sv; sums(0,1)=suv; sums(0,2)=svt;
-
- smatrix(0,0) = fSumXY[0]; smatrix(1,1)=fSumXY[4]; smatrix(2,2)=fSumXY[5];
- smatrix(0,1) = fSumXY[1]; smatrix(0,2)=fSumXY[3]; smatrix(1,2)=fSumXY[7];
- sums(0,0) = fSumXY[2]; sums(0,1) =fSumXY[6]; sums(0,2) =fSumXY[8];
- smatrix.Invert();
- if (smatrix.IsValid()){
- for (Int_t i=0;i<3;i++)
- for (Int_t j=0;j<=i;j++){
- (*fCov)(i,j)=smatrix(i,j);
- }
- TMatrixD res = sums*smatrix;
- fParams[0] = res(0,0);
- fParams[1] = res(0,1);
- fParams[2] = res(0,2);
- TMatrixD tmp = res*sums.T();
- fChi2 += fSumYY - tmp(0,0);
- fNdf += fNUsed - 3;
- conv++;
- }
+
+
+Bool_t AliTrackFitterRieman::Update(){
//
- // XZ part
+ //
//
- if (1) {
- Double_t x0 = -fParams[1]/fParams[0];
- Double_t Rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
-
- for (Int_t i=0;i<fNUsed;i++){
- Double_t phi = TMath::ASin((fX[i]-x0)*Rm1);
- Double_t phi0 = TMath::ASin((0.-x0)*Rm1);
- Double_t weight = 1/fSz[i];
- weight *=weight;
- Double_t dphi = (phi-phi0)/Rm1;
- fSumXZ[0] +=weight;
- fSumXZ[1] +=weight*dphi;
- fSumXZ[2] +=weight*dphi*dphi;
- fSumXZ[3] +=weight*fZ[i];
- fSumXZ[4] +=weight*dphi*fZ[i];
- }
-
- TMatrixDSym smatrixz(2);
- smatrixz(0,0) = fSumXZ[0]; smatrixz(0,1) = fSumXZ[1]; smatrixz(1,1) = fSumXZ[2];
- smatrixz.Invert();
- TMatrixD sumsxz(1,2);
- if (smatrixz.IsValid()){
- sumsxz(0,0) = fSumXZ[3];
- sumsxz(0,1) = fSumXZ[4];
- TMatrixD res = sumsxz*smatrixz;
- fParams[3] = res(0,0);
- fParams[4] = res(0,1);
- fParams[5] = 0;
- for (Int_t i=0;i<2;i++)
- for (Int_t j=0;j<=i;j++){
- (*fCov)(i+3,j+3)=smatrixz(i,j);
- }
- TMatrixD tmp = res*sumsxz.T();
- fChi2 += fSumZZ - tmp(0,0);
- fNdf += fNUsed - 2;
- conv++;
- }
- }
- else {
- Double_t Rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
- fSumXZ[1] += fSumXZ[5]*Rm1*Rm1;
- fSumXZ[2] += fSumXZ[6]*Rm1*Rm1 + fSumXZ[8]*Rm1*Rm1*Rm1*Rm1;
- fSumXZ[4] += fSumXZ[7]*Rm1*Rm1;
-
- TMatrixDSym smatrixz(2);
- smatrixz(0,0) = fSumXZ[0]; smatrixz(0,1) = fSumXZ[1]; smatrixz(1,1) = fSumXZ[2];
- smatrixz.Invert();
- TMatrixD sumsxz(1,2);
- if (smatrixz.IsValid()){
- sumsxz(0,0) = fSumXZ[3];
- sumsxz(0,1) = fSumXZ[4];
- TMatrixD res = sumsxz*smatrixz;
- fParams[3] = res(0,0);
- fParams[4] = res(0,1);
- fParams[5] = 0;
- for (Int_t i=0;i<2;i++)
- for (Int_t j=0;j<=i;j++){
- (*fCov)(i+3,j+3)=smatrixz(i,j);
- }
- TMatrixD tmp = res*sumsxz.T();
- fChi2 += fSumZZ - tmp(0,0);
- fNdf += fNUsed - 2;
- conv++;
+ 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;
}
-
- // (x-x0)^2+(y-y0)^2-R^2=0 ===>
//
- // (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0; ==>
- // substitution t = 1/(x^2+y^2), u = 2*x*t, y = 2*y*t, D0 = R^2 - x0^2- y0^2
//
- // 1 - u*x0 - v*y0 - t *D0 =0 ; - linear equation
- //
- // next substition a = 1/y0 b = -x0/y0 c = -D0/y0
- // final linear equation : a + u*b +t*c - v =0;
//
- // fParam[0] = 1/y0
- // fParam[1] = -x0/y0
- // fParam[2] = - (R^2 - x0^2 - y0^2)/y0
- if (conv>1) fConv =kTRUE;
- else
- fConv=kFALSE;
-}
-
-Double_t AliTrackFitterRieman::GetYat(Double_t x){
- if (!fConv) return 0.;
- Double_t res2 = (x*fParams[0]+fParams[1]);
- res2*=res2;
- res2 = 1.-fParams[2]*fParams[0]+fParams[1]*fParams[1]-res2;
- if (res2<0) return 0;
- res2 = TMath::Sqrt(res2);
- res2 = (1-res2)/fParams[0];
- return res2;
-}
-
-Double_t AliTrackFitterRieman::GetDYat(Double_t x) const {
- if (!fConv) return 0.;
- Double_t x0 = -fParams[1]/fParams[0];
- if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<0) return 0;
- Double_t Rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
- if ( 1./(Rm1*Rm1)-(x-x0)*(x-x0)<=0) return 0;
- Double_t res = (x-x0)/TMath::Sqrt(1./(Rm1*Rm1)-(x-x0)*(x-x0));
- if (fParams[0]<0) res*=-1.;
- return res;
-}
-
-
-
-Double_t AliTrackFitterRieman::GetZat(Double_t x) const {
- if (!fConv) return 0.;
- Double_t x0 = -fParams[1]/fParams[0];
- if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
- Double_t Rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
- Double_t phi = TMath::ASin((x-x0)*Rm1);
- Double_t phi0 = TMath::ASin((0.-x0)*Rm1);
- Double_t dphi = (phi-phi0);
- Double_t res = fParams[3]+fParams[4]*dphi/Rm1;
- return res;
-}
-
-Double_t AliTrackFitterRieman::GetDZat(Double_t x) const {
- if (!fConv) return 0.;
- Double_t x0 = -fParams[1]/fParams[0];
- if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
- Double_t Rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
- Double_t res = fParams[4]/TMath::Cos(TMath::ASin((x-x0)*Rm1));
- return res;
-}
-
+ 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);
-Double_t AliTrackFitterRieman::GetC(){
- return fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
+ return kTRUE;
}
-Bool_t AliTrackFitterRieman::GetXYZat(Double_t r, Float_t *xyz){
- if (!fConv) return kFALSE;
- Double_t res2 = (r*fParams[0]+fParams[1]);
- res2*=res2;
- res2 = 1.-fParams[2]*fParams[0]+fParams[1]*fParams[1]-res2;
- if (res2<0) return kFALSE;
- res2 = TMath::Sqrt(res2);
- res2 = (1-res2)/fParams[0];
- if (!fConv) return kFALSE;
- Double_t x0 = -fParams[1]/fParams[0];
- if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
- Double_t Rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
- Double_t phi = TMath::ASin((r-x0)*Rm1);
- Double_t phi0 = TMath::ASin((0.-x0)*Rm1);
- Double_t dphi = (phi-phi0);
- Double_t res = fParams[3]+fParams[4]*dphi/Rm1;
-
- Double_t sin = TMath::Sin(fAlpha);
- Double_t cos = TMath::Cos(fAlpha);
- xyz[0] = r*cos - res2*sin;
- xyz[1] = res2*cos + r*sin;
- xyz[2] = res;
-
- 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 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)/
+ Double_t r = TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1)/
fParams[0];
// Define space-point refence plane
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) return kFALSE;
- Double_t temp = TMath::Sqrt(R*R - (x-x0p)*(x-x0p));
+ 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;
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+
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;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff