Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages

MatrixCalculator Class Reference

#include <MatrixCalculator.h>

List of all members.

Public Member Functions

 MatrixCalculator ()
 MatrixCalculator (const AlgConfig &ac, const TrackContext &tc)
 ~MatrixCalculator ()
Int_t Solve (const DataFT &data)
TMatrixD GetFitErrorMatrix () const
double GetFitErrorMatrix (int i, int j) const
TVectorD GetTrackOut () const
double GetTrackOut (int i) const
double GetChi2 () const
double GetDChi2 () const
FitResult GetFitResult () const

Private Member Functions

Int_t MakePropagatorMatrix (const DataFT &data)
Int_t MakeCovarianceMatrix (const DataFT &data)
Double_t DiagonalElement (Int_t i, Int_t n, const TVectorD &theta_i, const DataFT &data) const
Double_t NonDiagonalElement (Int_t i, Int_t k, Int_t n, const TVectorD &theta_i, const DataFT &data) const

Private Attributes

TMatrixD fC
TVectorD fRes
TMatrixD fA
TMatrixD fV
TMatrixDSym fW
TMatrixD fFitCovM
TMatrixD fFitErrM
double fChi2
double fDChi2
TVectorD fTrackIn
TVectorD fTrackOut
Int_t fNPlanesUsed

Constructor & Destructor Documentation

MatrixCalculator::MatrixCalculator (   ) 
 

default constructor

Definition at line 31 of file MatrixCalculator.cxx.

00031                                    : 
00032         fFitCovM(NTrackParams, NTrackParams), 
00033         fFitErrM(NTrackParams, NTrackParams),
00034         fTrackIn(NTrackParams), fTrackOut(NTrackParams)
00035 {}

MatrixCalculator::MatrixCalculator (  const AlgConfig &  ac,
const TrackContext &  tc
) 
 

constructor

Definition at line 41 of file MatrixCalculator.cxx.

00041                                                                                       :
00042         fFitCovM(NTrackParams, NTrackParams), 
00043         fFitErrM(NTrackParams, NTrackParams),
00044         fTrackIn(NTrackParams), fTrackOut(NTrackParams)
00045 {}

MatrixCalculator::~MatrixCalculator (   ) 
 

destructor

Definition at line 51 of file MatrixCalculator.cxx.

00052 {}

Member Function Documentation

Double_t MatrixCalculator::DiagonalElement (  Int_t  i,
Int_t  n,
const TVectorD &  theta_i,
const DataFT &  data
)  const [private]
 

calcuate a diagonal element of the covariance matrix

Definition at line 321 of file MatrixCalculator.cxx.

References DataFT::GetCos(), DataFT::GetdZSteel(), and DataFT::T().

Referenced by MakeCovarianceMatrix().

00323 {
00324     return  pow(theta_i(i),2) *
00325         ( pow(data.GetdZSteel(i)/data.GetCos(i),2)/3.+
00326           data.GetdZSteel(i)/data.GetCos(i)*data.T(i,n) +
00327           pow(data.T(i,n),2)                              );
00328 }

double MatrixCalculator::GetChi2 (   )  const [inline]
 

get the fit chi-square

Definition at line 59 of file MatrixCalculator.h.

00059 { return fChi2; };

double MatrixCalculator::GetDChi2 (   )  const [inline]
 

get the "delta chi-square" between the input and output vectors of the fit parameters

Definition at line 65 of file MatrixCalculator.h.

Referenced by FitStateIterating::Iterate().

00065 { return fDChi2; };

double MatrixCalculator::GetFitErrorMatrix (  int  i,
int  j
)  const [inline]
 

get an [i,j] element of the error matrix

Definition at line 44 of file MatrixCalculator.h.

00044 { return fFitErrM(i,j); };

TMatrixD MatrixCalculator::GetFitErrorMatrix (   )  const [inline]
 

get the error matrix of fit parameters

Definition at line 39 of file MatrixCalculator.h.

00039 { return fFitErrM; };

FitResult MatrixCalculator::GetFitResult (   )  const
 

get fit results

Definition at line 57 of file MatrixCalculator.cxx.

References fChi2, fDChi2, fFitErrM, fNPlanesUsed, fTrackOut, and fV.

Referenced by FitStateConverged::Iterate().

00058 {
00059     return FitResult(   fFitErrM, fTrackOut, fChi2, 
00060                         fDChi2, fNPlanesUsed, fV.GetNcols() );
00061 }

double MatrixCalculator::GetTrackOut (  int  i  )  const [inline]
 

get the i-th element of the result vector

Definition at line 54 of file MatrixCalculator.h.

00054 { return fTrackOut(i); };

TVectorD MatrixCalculator::GetTrackOut (   )  const [inline]
 

get the result vector of fit parameters

Definition at line 49 of file MatrixCalculator.h.

Referenced by FitStateIterating::Iterate(), and FitStateInitial::Iterate().

00049 { return fTrackOut; };

Int_t MatrixCalculator::MakeCovarianceMatrix (  const DataFT &  data  )  [private]
 

calculate the covariance matrix

Definition at line 213 of file MatrixCalculator.cxx.

References DiagonalElement(), fV, DataFT::GetNPlanesUsed(), DataFT::GetNUHitsUsed(), DataFT::GetNVHitsUsed(), DataFT::GetSigmaU(), DataFT::GetSigmaV(), NonDiagonalElement(), DataFT::ThetaMCS(), DataFT::UHitUse(), and DataFT::VHitUse().

Referenced by Solve().

00214 {
00215     Int_t nplanes;
00216     nplanes = data.GetNPlanesUsed();
00217 
00218     Int_t nuhits, nvhits, nhits;
00219     nuhits = data.GetNUHitsUsed();
00220     nvhits = data.GetNVHitsUsed();
00221     nhits = nuhits+nvhits;
00222     fV.ResizeTo(nhits, nhits);
00223     fV.Zero();
00224 
00225     // Calculate MCS scattering angles for each plane
00226     TVectorD theta_i(nplanes);
00227     for (Int_t i = 0; i<nplanes; i++) {
00228         theta_i(i) = data.ThetaMCS(i);
00229     }
00230 
00231     // Calculate MCS covariance matrix
00232     Double_t diag, non_diag;
00233     Int_t index1, index2;
00234     // U part
00235     index1 = 0;
00236     for (Int_t n = 0; n<nplanes; n++) {
00237         if ( data.UHitUse(n) ) {
00238             // calculate diagonal elements
00239             diag = 0.0;
00240             for (Int_t i = 0; i<=n; i++) {
00241                 diag += DiagonalElement(i, n, theta_i, data);
00242             }
00243             fV(index1,index1) = diag;
00244 
00245             // non-diagonal elements
00246             index2 = index1+1;
00247             for (Int_t k = n+1; k<nplanes; k++) {
00248                 if ( data.UHitUse(k) ) {
00249                     non_diag = 0.0;
00250                     for (Int_t i = 0; i<=n; i++) {
00251                         non_diag += NonDiagonalElement(i, k, n, theta_i, data);
00252                     }
00253                     fV(index1,index2) = non_diag;
00254                     fV(index2,index1) = non_diag;
00255                     index2++;
00256                 }
00257             }
00258             index1++;
00259         }
00260     }
00261     // V part
00262     index1 = 0;
00263     for (Int_t n = 0; n<nplanes; n++) {
00264         if ( data.VHitUse(n) ) {
00265             // calculate diagonal elements
00266             diag = 0.0;
00267             for (Int_t i = 0; i<=n; i++) {
00268                 diag += DiagonalElement(i, n, theta_i, data);
00269             }
00270             fV(index1+nuhits,index1+nuhits) = diag;
00271 
00272             // non-diagonal elements
00273             index2 = index1+1;
00274             for (Int_t k = n+1; k<nplanes; k++) {
00275                 if ( data.VHitUse(k) ) {
00276                     non_diag = 0.0;
00277                     for (Int_t i = 0; i<=n; i++) {
00278                         non_diag += NonDiagonalElement(i, k, n, theta_i, data);
00279                     }
00280                     fV(index1+nuhits,index2+nuhits) = non_diag;
00281                     fV(index2+nuhits,index1+nuhits) = non_diag;
00282                     index2++;
00283                 }
00284             }
00285             index1++;
00286         }
00287     }
00288     
00289     // Add resolutions
00290     Int_t i_hit = 0;
00291     for (Int_t i = 0; i<nplanes; i++) {
00292         if ( data.UHitUse(i) ) {
00293             fV(i_hit,i_hit) += pow(data.GetSigmaU(i),2);
00294             i_hit++;
00295         }
00296     }
00297     for (Int_t i = 0; i<nplanes; i++) {
00298         if ( data.VHitUse(i) ) {
00299             fV(i_hit,i_hit) += pow(data.GetSigmaV(i),2);
00300             i_hit++;
00301         }
00302     }
00303 
00304 //     MsgStream *mftsa = &MSGSTREAM("FitTrackSA", Msg::kVerbose);
00305 //     (*mftsa) << "Covariance Matrix: \n";
00306 //     MsgFormat ffmt("%10.3e");
00307 // 
00308 //     for (Int_t i=0; i<nhits; i++) {
00309 //         for (Int_t j=0; j<nhits; j++) {
00310 //             (*mftsa) <<  ffmt(fV(i,j));
00311 //         }
00312 //         (*mftsa) << "\n";
00313 //     }
00314 
00315     return 0;
00316 }

Int_t MatrixCalculator::MakePropagatorMatrix (  const DataFT &  data  )  [private]
 

calculate the propagator matrix

Definition at line 167 of file MatrixCalculator.cxx.

References fA, fTrackIn, DataFT::GetNPlanesUsed(), DataFT::GetNUHitsUsed(), DataFT::GetNVHitsUsed(), DataFT::GetUf(), DataFT::GetVf(), DataFT::GetZ(), DataFT::UHitUse(), and DataFT::VHitUse().

Referenced by Solve().

00168 {
00169     Int_t nuhits, nvhits;
00170     nuhits = data.GetNUHitsUsed();
00171     nvhits = data.GetNVHitsUsed();
00172 
00173     fA.ResizeTo(nuhits+nvhits, NTrackParams);
00174 
00175     Int_t nplanes = data.GetNPlanesUsed();
00176     Int_t ihits_u = 0;
00177     //Protect against rare arithmetic exception
00178     Double_t epsilon = 0.00001;
00179     if (TMath::Abs(fTrackIn(kQoverP))<1e-6) 
00180       fTrackIn(kQoverP) = epsilon;
00181 
00182     for (Int_t i=0; i<nplanes; i++) {
00183         if ( data.UHitUse(i) ) {
00184             fA(ihits_u,kU) = 1.;
00185             fA(ihits_u,kdUdZ) = data.GetZ(i) - data.GetZ(0);
00186             fA(ihits_u,kV) = 0.;
00187             fA(ihits_u,kdVdZ) = 0.;
00188             fA(ihits_u,kQoverP) = (data.GetUf(i) - fTrackIn(kU) -
00189                fTrackIn(kdUdZ)*(data.GetZ(i)-data.GetZ(0)))/fTrackIn(kQoverP);
00190             ihits_u++;
00191         }
00192     }
00193 
00194     Int_t ihits_v = 0;
00195     for (Int_t i = 0; i<nplanes; i++) {
00196         if (data.VHitUse(i) ) {
00197             fA(ihits_v+nuhits,kV) = 1.;
00198             fA(ihits_v+nuhits,kdVdZ) = data.GetZ(i) - data.GetZ(0);
00199             fA(ihits_v+nuhits,kU) = 0.;
00200             fA(ihits_v+nuhits,kdUdZ) = 0.;
00201             fA(ihits_v+nuhits,kQoverP) = (data.GetVf(i) - fTrackIn(kV) -
00202                fTrackIn(kdVdZ)*(data.GetZ(i)-data.GetZ(0)))/fTrackIn(kQoverP);
00203             ihits_v++;
00204         }
00205     }
00206     
00207     return 0;
00208 }

Double_t MatrixCalculator::NonDiagonalElement (  Int_t  i,
Int_t  k,
Int_t  n,
const TVectorD &  theta_i,
const DataFT &  data
)  const [private]
 

calcuate a non-diagonal element of the covariance matrix

Definition at line 333 of file MatrixCalculator.cxx.

References DataFT::GetCos(), DataFT::GetdZSteel(), DataFT::GetZ(), and DataFT::T().

Referenced by MakeCovarianceMatrix().

00335 {
00336     return pow(theta_i(i),2) *
00337         (  pow(data.GetdZSteel(i)/data.GetCos(i),2)/3.+
00338            data.GetdZSteel(i)/data.GetCos(i)*data.T(i,n) +
00339            pow(data.T(i,n),2) +
00340            TMath::Abs(data.GetZ(k)-data.GetZ(n))*
00341            (data.GetdZSteel(i)/data.GetCos(i)/2.+data.T(i,n))  );
00342 }

Int_t MatrixCalculator::Solve (  const DataFT &  data  ) 
 

The main method called by AlgFitTrackSA - based on the track data from DataFT object, calculate covariance and propagator matrices and solve the matrix equation; return error status - 0 is succes, otherwise failure

Definition at line 91 of file MatrixCalculator.cxx.

References fA, fC, fChi2, fDChi2, fFitCovM, fFitErrM, DataFT::FillVectorC(), DataFT::FillVectorRes(), fNPlanesUsed, fTrackIn, fTrackOut, fV, DataFT::GetNPlanesUsed(), DataFT::GetTrack(), MakeCovarianceMatrix(), MakePropagatorMatrix(), and MSGSTREAM.

Referenced by FitStateIterating::Iterate(), and FitStateInitial::Iterate().

00092 {
00093     fTrackIn = data.GetTrack();
00094     
00095     fNPlanesUsed = data.GetNPlanesUsed();
00096         
00097     MakePropagatorMatrix(data);
00098     
00099     TMatrixD At(TMatrixD::kTransposed, fA);
00100 
00101     MakeCovarianceMatrix(data);
00102 
00103     //fW.ResizeTo(fV.GetNrows(), fV.GetNcols());
00104 
00105 #if ROOT_VERSION_CODE >= ROOT_VERSION(4,0,0)
00106     // using Cholesky decomposition here - works on symmetric
00107     // matrices TMatrixDSym only and should be faster than
00108     // regular inversion
00109     // TDecompChol chol(fV); 
00110     // chol.Invert(fW);
00111     
00112     // fW = TMatrixD(TMatrixD::kInverted, fV);
00113     
00114     TMatrixD W(TMatrixD::kInverted, fV);
00115 #else
00116     // fW = TMatrixD(TMatrixD::kInvertedPosDef, fV);
00117     
00118     TMatrixD W(TMatrixD::kInverted, fV);
00119 #endif
00120         
00121     fFitCovM = TMatrixD( TMatrixD(At,TMatrixD::kMult,W), TMatrixD::kMult, fA);
00122 
00123 #if ROOT_VERSION_CODE >= ROOT_VERSION(4,0,0)
00124     fFitErrM = TMatrixD(TMatrixD::kInverted, fFitCovM);
00125 #else    
00126     fFitErrM = TMatrixD(TMatrixD::kInvertedPosDef, fFitCovM);
00127 #endif
00128 
00129     MsgStream *mftsa = &MSGSTREAM("FitTrackSA", Msg::kVerbose);
00130     (*mftsa) << "Solution Error Matrix:\n";
00131     MsgFormat efmt("%10.2e");
00132     for (Int_t i = 0; i<fFitErrM.GetNrows(); i++) {
00133         for (Int_t j = 0; j<fFitErrM.GetNcols(); j++) {
00134             (*mftsa) << efmt(fFitErrM(i,j));
00135         }
00136         (*mftsa) << "\n";
00137     }
00138 
00139     data.FillVectorC(fC);
00140 
00141     TMatrixD solution( fFitErrM, TMatrixD::kMult, 
00142         TMatrixD(At, TMatrixD::kMult, TMatrixD(W,TMatrixD::kMult,fC)) );
00143 
00144     fTrackOut = TMatrixDColumn(solution, 0);
00145 
00146     TMatrixD  vRes(1,1);
00147     data.FillVectorRes(vRes);
00148     TMatrixD chi2(TMatrixD(TMatrixD::kTransposed,vRes), 
00149                     TMatrixD::kMult, TMatrixD(W, TMatrixD::kMult, vRes));
00150     fChi2 = chi2(0,0);
00151 
00152     TVectorD vres(fTrackOut);
00153     for (Int_t i = 0; i<NTrackParams; i++) {
00154         vres(i) -= fTrackIn(i);
00155     }
00156     TVectorD vresT(vres);
00157     vresT *= fFitCovM;
00158     
00159     fDChi2 = vresT*vres;
00160     
00161     return 0;
00162 }

Member Data Documentation

TMatrixD MatrixCalculator::fA [private]
 

propagation matrix

Definition at line 109 of file MatrixCalculator.h.

Referenced by MakePropagatorMatrix(), and Solve().

TMatrixD MatrixCalculator::fC [private]
 

vector of hit coordinates (single column TMatrixD)

Definition at line 99 of file MatrixCalculator.h.

Referenced by Solve().

double MatrixCalculator::fChi2 [private]
 

fit chi-square

Definition at line 134 of file MatrixCalculator.h.

Referenced by GetFitResult(), and Solve().

double MatrixCalculator::fDChi2 [private]
 

"delta chi-square" between the input and output vectors of the fit parameters

Definition at line 140 of file MatrixCalculator.h.

Referenced by GetFitResult(), and Solve().

TMatrixD MatrixCalculator::fFitCovM [private]
 

fit parameters covariance matrix

Definition at line 124 of file MatrixCalculator.h.

Referenced by Solve().

TMatrixD MatrixCalculator::fFitErrM [private]
 

fit parameters error matrix

Definition at line 129 of file MatrixCalculator.h.

Referenced by GetFitResult(), and Solve().

Int_t MatrixCalculator::fNPlanesUsed [private]
 

number of planes used in the fit

Definition at line 155 of file MatrixCalculator.h.

Referenced by GetFitResult(), and Solve().

TVectorD MatrixCalculator::fRes [private]
 

vector of hit residuals

Definition at line 104 of file MatrixCalculator.h.

TVectorD MatrixCalculator::fTrackIn [private]
 

initial vector of the fit parameters (u, du/dz, v, dv/dz, q/p)

Definition at line 145 of file MatrixCalculator.h.

Referenced by MakePropagatorMatrix(), and Solve().

TVectorD MatrixCalculator::fTrackOut [private]
 

solution vector of the fit parameters (u, du/dz, v, dv/dz, q/p)

Definition at line 150 of file MatrixCalculator.h.

Referenced by GetFitResult(), and Solve().

TMatrixD MatrixCalculator::fV [private]
 

covariance matrix

Definition at line 114 of file MatrixCalculator.h.

Referenced by GetFitResult(), MakeCovarianceMatrix(), and Solve().

TMatrixDSym MatrixCalculator::fW [private]
 

weight matrix (inverse of covariance matrix)

Definition at line 119 of file MatrixCalculator.h.

The documentation for this class was generated from the following files:
Generated on Mon Feb 15 11:09:30 2010 for loon by  doxygen 1.3.9.1