// neglecting error of xi, and supposing xi>>yi sigmai ~ sigmaVi ~ 2*sigmay*t
-#include "TMatrixDSym.h"
-//#include "TDecompChol.h"
-#include "TMatrixD.h"
+#include <TMatrixDSym.h>
+#include <TMath.h>
+#include <TMatrixD.h>
+
#include "AliRieman.h"
ClassImp(AliRieman)
-AliRieman::AliRieman(){
+AliRieman::AliRieman() :
+ TObject(),
+ fCapacity(0),
+ fN(0),
+ fX(0),
+ fY(0),
+ fZ(0),
+ fSy(0),
+ fSz(0),
+ fCovar(0),
+ fCovarPolY(0),
+ fCovarPolZ(0),
+ fSumZZ(0),
+ fChi2(0),
+ fChi2Y(0),
+ fChi2Z(0),
+ fConv(kFALSE)
+{
//
// default constructor
//
for (Int_t i=0;i<6;i++) fParams[i] = 0;
for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
- fSumZZ = 0;
- fCapacity = 0;
- fN =0;
- fX = 0;
- fY = 0;
- fZ = 0;
- fSy = 0;
- fSz = 0;
- fChi2 = 0;
- fChi2Y = 0;
- fChi2Z = 0;
- fCovar = 0;
- fConv = kFALSE;
+ for (Int_t i=0;i<6;i++) {
+ fSumPolY[i]=0;
+ fSumPolZ[i]=0;
+ }
}
-AliRieman::AliRieman(Int_t capacity){
+AliRieman::AliRieman(Int_t capacity) :
+ TObject(),
+ fCapacity(capacity),
+ fN(0),
+ fX(new Double_t[fCapacity]),
+ fY(new Double_t[fCapacity]),
+ fZ(new Double_t[fCapacity]),
+ fSy(new Double_t[fCapacity]),
+ fSz(new Double_t[fCapacity]),
+ fCovar(new TMatrixDSym(6)),
+ fCovarPolY(new TMatrixDSym(3)),
+ fCovarPolZ(new TMatrixDSym(2)),
+ fSumZZ(0),
+ fChi2(0),
+ fChi2Y(0),
+ fChi2Z(0),
+ fConv(kFALSE)
+{
//
// default constructor
//
for (Int_t i=0;i<6;i++) fParams[i] = 0;
for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
- fSumZZ =0;
- fCapacity = capacity;
- fN =0;
- fX = new Double_t[fCapacity];
- fY = new Double_t[fCapacity];
- fZ = new Double_t[fCapacity];
- fSy = new Double_t[fCapacity];
- fSz = new Double_t[fCapacity];
- fCovar = new TMatrixDSym(6);
- fChi2 = 0;
- fChi2Y = 0;
- fChi2Z = 0;
- fConv = kFALSE;
+ for (Int_t i=0;i<6;i++) {
+ fSumPolY[i]=0;
+ fSumPolZ[i]=0;
+ }
}
-AliRieman::AliRieman(const AliRieman &rieman):TObject(rieman){
+AliRieman::AliRieman(const AliRieman &rieman):
+ TObject(rieman),
+ fCapacity(rieman.fN),
+ fN(rieman.fN),
+ fX(new Double_t[fN]),
+ fY(new Double_t[fN]),
+ fZ(new Double_t[fN]),
+ fSy(new Double_t[fN]),
+ fSz(new Double_t[fN]),
+ fCovar(new TMatrixDSym(*(rieman.fCovar))),
+ fCovarPolY(new TMatrixDSym(*(rieman.fCovarPolY))),
+ fCovarPolZ(new TMatrixDSym(*(rieman.fCovarPolZ))),
+ fSumZZ(rieman.fSumZZ),
+ fChi2(rieman.fChi2),
+ fChi2Y(rieman.fChi2Y),
+ fChi2Z(rieman.fChi2Z),
+ fConv(rieman.fConv)
+
+{
//
// copy constructor
//
for (Int_t i=0;i<6;i++) fParams[i] = rieman.fParams[i];
for (Int_t i=0;i<9;i++) fSumXY[i] = rieman.fSumXY[i];
for (Int_t i=0;i<9;i++) fSumXZ[i] = rieman.fSumXZ[i];
- fSumZZ = rieman.fSumZZ;
- fCapacity = rieman.fN;
- fN =rieman.fN;
- fCovar = new TMatrixDSym(*(rieman.fCovar));
- fConv = rieman.fConv;
- fX = new Double_t[fN];
- fY = new Double_t[fN];
- fZ = new Double_t[fN];
- fSy = new Double_t[fN];
- fSz = new Double_t[fN];
+ for (Int_t i=0;i<6;i++) {
+ fSumPolY[i]=rieman.fSumPolY[i];
+ fSumPolZ[i]=rieman.fSumPolZ[i];
+ }
for (Int_t i=0;i<rieman.fN;i++){
fX[i] = rieman.fX[i];
fY[i] = rieman.fY[i];
delete[]fSy;
delete[]fSz;
delete fCovar;
+ delete fCovarPolY;
+ delete fCovarPolZ;
}
void AliRieman::Reset()
fN=0;
for (Int_t i=0;i<6;i++) fParams[i] = 0;
for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
- for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
+ for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
+ for (Int_t i=0;i<6;i++) {
+ fSumPolY[i]=0;
+ fSumPolZ[i]=0;
+ }
fSumZZ =0;
fConv =kFALSE;
}
fSumXZ[8] +=weight*r*r*r*r*r*r/(24*24);
fSumZZ += z*z*weight;
//
+ // sum accumulation for rough error estimates of the track extrapolation error
+ //
+ Double_t maty = 1./(sy*sy);
+ Double_t matz = 1./(sz*sz);
+ for (Int_t i=0;i<5; i++){
+ fSumPolY[i] += maty;
+ fSumPolZ[i] += matz;
+ maty *= x;
+ matz *= x;
+ }
fN++;
}
-
-
-
void AliRieman::UpdatePol(){
//
// Rieman update
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];
+ smatrix(1,0) = fSumXY[1]; smatrix(2,0)=fSumXY[3]; smatrix(2,1)=fSumXY[7];
sums(0,0) = fSumXY[2]; sums(0,1) =fSumXY[6]; sums(0,2) =fSumXY[8];
smatrix.Invert();
if (smatrix.IsValid()){
//
// XZ part
//
- 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;
+ 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;
+ Double_t sum1 = fSumXZ[1] + fSumXZ[5]*rm1*rm1;
+ Double_t sum2 = fSumXZ[2] + fSumXZ[6]*rm1*rm1 + fSumXZ[8]*rm1*rm1*rm1*rm1;
+ Double_t sum4 = 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(1,0) = fSumXZ[1];
+ // smatrixz(0,0) = fSumXZ[0]; smatrixz(0,1) = fSumXZ[1]; smatrixz(1,1) = fSumXZ[2];
+ //smatrixz(1,0) = fSumXZ[1];
+ smatrixz(0,0) = fSumXZ[0]; smatrixz(0,1) = sum1; smatrixz(1,1) = sum2;
+ smatrixz(1,0) = sum1;
smatrixz.Invert();
TMatrixD sumsxz(1,2);
if (smatrixz.IsValid()){
sumsxz(0,0) = fSumXZ[3];
- sumsxz(0,1) = fSumXZ[4];
+ // sumsxz(0,1) = fSumXZ[4];
+ sumsxz(0,1) = sum4;
TMatrixD res = sumsxz*smatrixz;
fParams[3] = res(0,0);
fParams[4] = res(0,1);
}
conv++;
}
+ UpdateCovariancePol();
// (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; ==>
}
conv++;
}
-
+ UpdateCovariancePol();
// (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; ==>
fChi2 = fChi2Y +fChi2Z;
}
+void AliRieman::UpdateCovariancePol(){
+ //
+ // covariance matrices for rough extrapolation error estimate
+ // covariance matrix - get error estimates at point x = 0
+ // ! NOTE - error estimates is very rough and it is valid only if dl<<R
+ // when dl is the distance between first and last point
+ //
+ //
+ (*fCovarPolY)(0,0) = fSumPolY[0]; (*fCovarPolY)(0,1) = fSumPolY[1]; (*fCovarPolY)(0,2) = fSumPolY[2];
+ (*fCovarPolY)(1,0) = fSumPolY[1]; (*fCovarPolY)(1,1) = fSumPolY[2]; (*fCovarPolY)(1,2) = fSumPolY[3];
+ (*fCovarPolY)(2,0) = fSumPolY[2]; (*fCovarPolY)(2,1) = fSumPolY[3]; (*fCovarPolY)(2,2) = fSumPolY[4];
+ fCovarPolY->Invert();
+ //
+
+ (*fCovarPolZ)(0,0) = fSumPolZ[0]; (*fCovarPolZ)(0,1) = fSumPolZ[1];
+ (*fCovarPolZ)(1,0) = fSumPolZ[1]; (*fCovarPolZ)(1,1) = fSumPolZ[2];
+ fCovarPolZ->Invert();
+ //
+}
+
+Double_t AliRieman::GetErrY(Double_t x) const{
+ //
+ // P0' = P0 + P1 * x + P2 * x^2
+ // P1' = P1 + P2 * x
+ // P2' = + P2
+ TMatrixD trans(3,3);
+ trans(0,0) = 1;
+ trans(0,1) = x;
+ trans(0,2) = x*x;
+ trans(1,1) = 1;
+ trans(1,2) = x;
+ trans(2,2) = 1;
+ //
+ TMatrixD covarX(trans,TMatrixD::kMult,*fCovarPolY);
+ covarX*=trans.T();
+ return covarX(0,0);
+}
+
+Double_t AliRieman::GetErrZ(Double_t x) const{
+ //
+ // assumption error of curvature determination neggligible
+ //
+ // P0' = P0 + P1 * x
+ // P1' = P1
+ TMatrixD trans(2,2);
+ trans(0,0) = 1;
+ trans(0,1) = x;
+ trans(1,1) = 1;
+ //
+ TMatrixD covarX(trans,TMatrixD::kMult,*fCovarPolZ);
+ covarX*=trans.T();
+ return covarX(0,0);
+}
+
+Bool_t AliRieman::GetExternalParameters(Double_t xref, Double_t *params, Double_t * covar){
+ //
+ // Get external parameters
+ // + approximative covariance
+ // 0 1 2 3 4 // y - local y
+ // 5 6 7 8 // z - local z
+ // 9 10 11 // snp - local sin fi
+ // 12 13 // tgl - deep angle
+ // 14 // C - curvature
+ Double_t sign = (GetC()>0) ? 1.:-1.;
+
+ params[0] = GetYat(xref);
+ params[1] = GetZat(xref);
+ params[2] = TMath::Sin(TMath::ATan(GetDYat(xref)));
+ params[3] = sign*fParams[4];
+ params[4] = GetC();
+ //
+ // covariance matrix in y
+ //
+ TMatrixD transY(3,3);
+ transY(0,0) = 1;
+ transY(0,1) = xref;
+ transY(0,2) = xref*xref;
+ transY(1,1) = 1;
+ transY(1,2) = xref;
+ transY(2,2) = 1;
+ TMatrixD covarY(transY,TMatrixD::kMult,*fCovarPolY);
+ covarY*=transY.T();
+ //
+ TMatrixD transZ(2,2);
+ transZ(0,0) = 1;
+ transZ(0,1) = xref;
+ transZ(1,1) = 1;
+ TMatrixD covarZ(transZ,TMatrixD::kMult,*fCovarPolZ);
+ covarZ*=transZ.T();
+ //
+ // C ~ 2*P2 - in rotated frame
+ //
+ covar[0] = covarY(0,0); covar[1] = 0; covar[2] = covarY(0,1); covar[3] = 0; covar[4] = TMath::Sqrt(2.)*covarY(0,2);
+ covar[5] = covarZ(0,0); covar[6] = 0; covar[7] = sign*covarZ(0,1); covar[8] = 0;
+ covar[9] = covarY(1,1); covar[10]= 0; covar[11]= TMath::Sqrt(2.)*covarY(1,2);
+ covar[12] = covarZ(1,1); covar[13]= 0;
+ covar[14] = 2.*covarY(2,2);
+ //
+ return fConv;
+}
+
+
Double_t AliRieman::GetYat(Double_t x) const {
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff