From: hristov Date: Thu, 20 Apr 2006 22:18:32 +0000 (+0000) Subject: Double precision (Marian). Coding conventions (Federico) X-Git-Url: http://git.uio.no/git/?a=commitdiff_plain;h=049179df0b5fdbe965c6067d9bf79fc0b3221232;p=u%2Fmrichter%2FAliRoot.git Double precision (Marian). Coding conventions (Federico) --- diff --git a/STEER/AliRieman.cxx b/STEER/AliRieman.cxx index aabad725094..956c10619c9 100755 --- a/STEER/AliRieman.cxx +++ b/STEER/AliRieman.cxx @@ -1,8 +1,50 @@ -#include "TTree.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 for global helix fit of a track +// Author: M.Ivanov +// The method uses the following idea: +//------------------------------------------------------ +// XY direction +// +// (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 + + #include "TMatrixDSym.h" +//#include "TDecompChol.h" #include "TMatrixD.h" #include "AliRieman.h" -#include "TRandom.h" ClassImp(AliRieman) @@ -23,6 +65,9 @@ AliRieman::AliRieman(){ fZ = 0; fSy = 0; fSz = 0; + fChi2 = 0; + fChi2Y = 0; + fChi2Z = 0; fCovar = 0; fConv = kFALSE; } @@ -38,16 +83,19 @@ AliRieman::AliRieman(Int_t capacity){ fCapacity = capacity; fN =0; - fX = new Float_t[fCapacity]; - fY = new Float_t[fCapacity]; - fZ = new Float_t[fCapacity]; - fSy = new Float_t[fCapacity]; - fSz = new Float_t[fCapacity]; + 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; } -AliRieman::AliRieman(const AliRieman &rieman):TObject(){ +AliRieman::AliRieman(const AliRieman &rieman):TObject(rieman){ // // copy constructor // @@ -58,11 +106,11 @@ AliRieman::AliRieman(const AliRieman &rieman):TObject(){ fN =rieman.fN; fCovar = new TMatrixDSym(*(rieman.fCovar)); fConv = rieman.fConv; - fX = new Float_t[fN]; - fY = new Float_t[fN]; - fZ = new Float_t[fN]; - fSy = new Float_t[fN]; - fSz = new Float_t[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]; for (Int_t i=0;i // - // substitution t = 1/(x^2+y^2), u = 2*x*t, y = 2*y*t, D0 = R^2 - x0^2- y0^2 + // 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 // @@ -246,7 +300,16 @@ void AliRieman::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]; + //TDecompChol decomp(smatrix); + //decomp.SetTol(1.0e-14); + //smatrix = + //decomp.Invert(); smatrix.Invert(); if (smatrix.IsValid()){ for (Int_t i=0;i<3;i++) @@ -263,20 +326,24 @@ void AliRieman::Update(){ fConv = kFALSE; //non converged return; } + if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1.<0){ + fConv = kFALSE; //non converged + return; + } // // XZ part // Double_t x0 = -fParams[1]/fParams[0]; - Double_t Rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); + Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1.); Double_t sumXZ[9]; for (Int_t i=0;i<9;i++) sumXZ[i]=0; for (Int_t i=0;i // // (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 + // 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 // @@ -318,11 +386,17 @@ void AliRieman::Update(){ if (conv>1) fConv =kTRUE; else fConv=kFALSE; + fChi2Y = CalcChi2Y(); + fChi2Z = CalcChi2Z(); + fChi2 = fChi2Y +fChi2Z; } -Double_t AliRieman::GetYat(Double_t x){ +Double_t AliRieman::GetYat(Double_t x) const { + // + // get y at given x position + // if (!fConv) return 0.; Double_t res2 = (x*fParams[0]+fParams[1]); res2*=res2; @@ -333,43 +407,98 @@ Double_t AliRieman::GetYat(Double_t x){ return res2; } -Double_t AliRieman::GetDYat(Double_t x){ +Double_t AliRieman::GetDYat(Double_t x) const { + // + // get dy/dx at given x position + // 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)); + 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 AliRieman::GetZat(Double_t x){ +Double_t AliRieman::GetZat(Double_t x) const { + // + // get z at given x position + // 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 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; + Double_t res = fParams[3]+fParams[4]*dphi/rm1; return res; } -Double_t AliRieman::GetDZat(Double_t x){ +Double_t AliRieman::GetDZat(Double_t x) const { + // + // get dz/dx at given x postion + // 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)); + 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; } -Double_t AliRieman::GetC(){ +Double_t AliRieman::GetC() const{ + // + // get curvature + // return fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); } +Double_t AliRieman::CalcChi2Y() const{ + // + // calculate chi2 for Y + // + Double_t sumchi2 = 0; + for (Int_t i=0;iAddPoint(fX[i],fY[i]-GetYat(fX[i]),fZ[i]-GetZat(fX[i]), fSy[i],fSz[i]); + } + return rieman; +} + diff --git a/STEER/AliRieman.h b/STEER/AliRieman.h index 2b599c19ea9..59026d57f82 100644 --- a/STEER/AliRieman.h +++ b/STEER/AliRieman.h @@ -1,7 +1,16 @@ -#ifndef ALIRIEMANN_H -#define ALIRIEMANN_H +#ifndef ALIRIEMAN_H +#define ALIRIEMAN_H +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ +/* $Id$ */ +// Class for global helix fit of a track +// Author: M.Ivanov +// This class uses decomposition of the chi2 based on the fact that +// one can rotate the coordinate system and provide xi >> yi for each +// space point -#include "TMatrixDSym.h" + +#include class AliRieman : public TObject{ public: @@ -10,40 +19,53 @@ class AliRieman : public TObject{ AliRieman(const AliRieman &rieman); ~AliRieman(); void Reset(); - void AddPoint(Float_t x, Float_t y, Float_t z, Float_t sy, Float_t sz); + void AddPoint(Double_t x, Double_t y, Double_t z, Double_t sy, Double_t sz); Int_t GetN() const {return fN;} Int_t GetCapacity() const {return fCapacity;} - Float_t * GetX(){return fX;} - Float_t * GetY(){return fY;} - Float_t * GetZ(){return fZ;} - Float_t * GetSy(){return fSy;} - Float_t * GetSz(){return fSz;} + Double_t * GetX(){return fX;} + Double_t * GetY(){return fY;} + Double_t * GetZ(){return fZ;} + Double_t * GetSy(){return fSy;} + Double_t * GetSz(){return fSz;} void Update(); void UpdatePol(); Double_t* GetParam(){return fParams;} - const TMatrixDSym & GetCovariance(){return *fCovar;} - Double_t GetC(); - Double_t GetYat(Double_t x); - Double_t GetZat(Double_t x); - Double_t GetDYat(Double_t x); - Double_t GetDZat(Double_t x); + const TMatrixDSym & GetCovariance() const {return *fCovar;} + Double_t GetC() const; + Double_t GetYat(Double_t x) const; + Double_t GetZat(Double_t x) const; + Double_t GetDYat(Double_t x) const; + Double_t GetDZat(Double_t x) const; + // + Double_t GetChi2Y() const { return fChi2Y;} + Double_t GetChi2Z() const { return fChi2Z;} + Double_t GetChi2() const { return fChi2; } + + Double_t CalcChi2Y() const; + Double_t CalcChi2Z() const; + Double_t CalcChi2() const; + AliRieman * MakeResiduals() const; + // protected: // public: Int_t fCapacity; // capacity Int_t fN; // numebr of points - Float_t *fX; //[fN] x coordinate - Float_t *fY; //[fN] y coordinate - Float_t *fZ; //[fN] z coordinate - Float_t *fSy; //[fN] sigma y coordinate - Float_t *fSz; //[fN] sigma z coordinate + Double_t *fX; //[fN] x coordinate + Double_t *fY; //[fN] y coordinate + Double_t *fZ; //[fN] z coordinate + Double_t *fSy; //[fN] sigma y coordinate + Double_t *fSz; //[fN] sigma z coordinate Double_t fParams[6]; //Parameters TMatrixDSym *fCovar; //Covariance Double_t fSumXY[9]; //sums for XY part Double_t fSumXZ[9]; //sums for XZ part + Double_t fChi2; //sums of chi2 + Double_t fChi2Y; //sums of chi2 for y coord + Double_t fChi2Z; //sums of chi2 foz z coord Bool_t fConv; // indicates convergation - protected: - + protected: private: + AliRieman& operator=(const AliRieman &rieman){return *this;} ClassDef(AliRieman,1) // Fast fit of helices on ITS RecPoints };