~maus-scifi/maus/tracker_devel

« back to all changes in this revision

Viewing changes to src/common_cpp/Recon/SciFi/LeastSquaresFitter.cc

  • Committer: Christopher Hunt
  • Date: 2015-01-29 16:17:20 UTC
  • mto: This revision was merged to the branch mainline in revision 1143.
  • Revision ID: christopher.hunt08@imperial.ac.uk-20150129161720-n8pjzvdoffnsqj3o
Updatest to kalman. Use PatRec for state vector and covariance seeds

Show diffs side-by-side

added added

removed removed

Lines of Context:
102
102
  At.T();                          // Transpose At (leaving A unchanged)
103
103
  TMatrixD V_p(At * V_m * A);      // The covariance matrix of the parameters of model (inv)
104
104
  V_p.Invert(det);                 // Invert in place
105
 
  covariance = V_p;
106
105
  TMatrixD P(V_p * At * V_m * K);  // The least sqaures estimate of the parameters
107
106
 
108
107
  // Extract the fit parameters
120
119
  else
121
120
    R = sqrt((4 * alpha) + (beta * beta) + (gamma * gamma)) / (2 * alpha);
122
121
 
 
122
  // Transform the covariance matrix to the same basis
 
123
  TMatrixD jacobian( 3, 3 );
 
124
  jacobian(0,0) = beta / (2.0*alpha*alpha);
 
125
  jacobian(0,1) = -1.0 / (2.0*alpha);
 
126
  jacobian(1,0) = gamma / (2.0*alpha*alpha);
 
127
  jacobian(1,2) = -1.0 / (2.0*alpha);
 
128
  jacobian(2,0) = ( -1.0/(2.0*alpha) ) * ( ( (beta*beta + gamma*gamma) / (2.0*alpha) ) + 1 ) /
 
129
                                                sqrt( ( (beta*beta + gamma*gamma) / 4.0 ) + alpha );
 
130
  jacobian(2,1) = ( beta/(4.0*alpha*alpha) ) /
 
131
                            sqrt( ( (beta*beta + gamma*gamma)/(4.0*alpha*alpha) ) + ( 1.0/alpha ) );
 
132
  jacobian(2,2) = ( gamma/(4.0*alpha*alpha) ) /
 
133
                            sqrt( ( (beta*beta + gamma*gamma)/(4.0*alpha*alpha) ) + ( 1.0/alpha ) );
 
134
  TMatrixD jacobianT(3,3);
 
135
  jacobianT.Transpose( jacobian );
 
136
 
 
137
  covariance = jacobian * V_p * jacobianT;
 
138
 
123
139
  // if ( R < 0. )
124
140
  //  std::cout << "R was < 0 but taking abs_val for physical correctness\n";
125
141
  R = fabs(R);