From 8919d7600ce2f11d1767272ac53d8c56b668d80f Mon Sep 17 00:00:00 2001
From: Stephen Nicholas Swatman <stephen.nicholas.swatman@cern.ch>
Date: Tue, 20 Oct 2020 13:12:18 +0200
Subject: [PATCH] Fix aliasing bug in newCovarianceMatrix method
The `PatternTrackParameters` class conains a method by the name of
`newCovarianceMatrix` which updates the covariance matrix of the object
it is called on according to the covariance matrix in a second set of
parameters and a Jacobian matrix. The calculation performed is J*C*J^T.
In many cases where this method is called, the object on which it is
called is the same as the object given as the right hand side of the
operation. This means that the covariance matrices alias, leading to
possibly incorrect results. It also invites behaviour where the
calculations are done using invalid covariance matrices.
This commit fixes these issues by changing the signature of the method
to accept a covariance matrix directly instead of a second set of track
parameters. This makes it harder, albeit not impossible, to alias the
covariance matrices as described earlier. It also solves the invalid
covariance matrix problem described in !37151.
---
.../src/SiTrajectoryElement_xk.cxx | 148 +++++++++---------
.../PatternTrackParameters.h | 17 +-
.../src/PatternTrackParameters.cxx | 109 ++++++-------
.../src/RungeKuttaPropagator.cxx | 6 +-
4 files changed, 151 insertions(+), 129 deletions(-)
diff --git a/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/src/SiTrajectoryElement_xk.cxx b/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/src/SiTrajectoryElement_xk.cxx
index 5b20bf15667..8bd9a196ff2 100644
--- a/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/src/SiTrajectoryElement_xk.cxx
+++ b/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/src/SiTrajectoryElement_xk.cxx
@@ -1261,82 +1261,86 @@ bool InDet::SiTrajectoryElement_xk::transformGlobalToPlane
acos(globalPars[5]), /// theta (from global cosTheta)
globalPars[6] /// qoverp
};
- /// write into output parameters. Assign our surface to them
- outputParameters.setParameters(m_surface,p);
- /// update local direction vector
- m_localDir[0] = globalPars[3];
- m_localDir[1] = globalPars[4];
- m_localDir[2] = globalPars[5];
-
- /// if we don't need the cov, we are done
- if(!useJac) return true;
- /// Condition trajectory on surface
- double A = Az[0]*globalPars[3]+Az[1]*globalPars[4]+Az[2]*globalPars[5];
- if(A!=0.) A=1./A;
- double s0 = Az[0]*globalPars[ 7]+Az[1]*globalPars[ 8]+Az[2]*globalPars[ 9];
- double s1 = Az[0]*globalPars[14]+Az[1]*globalPars[15]+Az[2]*globalPars[16];
- double s2 = Az[0]*globalPars[21]+Az[1]*globalPars[22]+Az[2]*globalPars[23];
- double s3 = Az[0]*globalPars[28]+Az[1]*globalPars[29]+Az[2]*globalPars[30];
- double s4 = Az[0]*globalPars[35]+Az[1]*globalPars[36]+Az[2]*globalPars[37];
- double T0 =(Ax[0]*globalPars[ 3]+Ax[1]*globalPars[ 4]+Ax[2]*globalPars[ 5])*A;
- double T1 =(Ay[0]*globalPars[ 3]+Ay[1]*globalPars[ 4]+Ay[2]*globalPars[ 5])*A;
- double n = 1./globalPars[6];
+ /// update local direction vector
+ m_localDir[0] = globalPars[3];
+ m_localDir[1] = globalPars[4];
+ m_localDir[2] = globalPars[5];
- double Jac[21];
+ if (useJac) {
+ /// Condition trajectory on surface
+ double A = Az[0]*globalPars[3]+Az[1]*globalPars[4]+Az[2]*globalPars[5];
+ if(A!=0.) A=1./A;
+ double s0 = Az[0]*globalPars[ 7]+Az[1]*globalPars[ 8]+Az[2]*globalPars[ 9];
+ double s1 = Az[0]*globalPars[14]+Az[1]*globalPars[15]+Az[2]*globalPars[16];
+ double s2 = Az[0]*globalPars[21]+Az[1]*globalPars[22]+Az[2]*globalPars[23];
+ double s3 = Az[0]*globalPars[28]+Az[1]*globalPars[29]+Az[2]*globalPars[30];
+ double s4 = Az[0]*globalPars[35]+Az[1]*globalPars[36]+Az[2]*globalPars[37];
+ double T0 =(Ax[0]*globalPars[ 3]+Ax[1]*globalPars[ 4]+Ax[2]*globalPars[ 5])*A;
+ double T1 =(Ay[0]*globalPars[ 3]+Ay[1]*globalPars[ 4]+Ay[2]*globalPars[ 5])*A;
+ double n = 1./globalPars[6];
+
+ double Jac[21];
+
+ // Jacobian production
+ //
+ Jac[ 0] = (Ax[0]*globalPars[ 7]+Ax[1]*globalPars[ 8])+(Ax[2]*globalPars[ 9]-s0*T0); // dL0/dL0
+ Jac[ 1] = (Ax[0]*globalPars[14]+Ax[1]*globalPars[15])+(Ax[2]*globalPars[16]-s1*T0); // dL0/dL1
+ Jac[ 2] = (Ax[0]*globalPars[21]+Ax[1]*globalPars[22])+(Ax[2]*globalPars[23]-s2*T0); // dL0/dPhi
+ Jac[ 3] = (Ax[0]*globalPars[28]+Ax[1]*globalPars[29])+(Ax[2]*globalPars[30]-s3*T0); // dL0/dThe
+ Jac[ 4] =((Ax[0]*globalPars[35]+Ax[1]*globalPars[36])+(Ax[2]*globalPars[37]-s4*T0))*n; // dL0/dCM
+
+ Jac[ 5] = (Ay[0]*globalPars[ 7]+Ay[1]*globalPars[ 8])+(Ay[2]*globalPars[ 9]-s0*T1); // dL1/dL0
+ Jac[ 6] = (Ay[0]*globalPars[14]+Ay[1]*globalPars[15])+(Ay[2]*globalPars[16]-s1*T1); // dL1/dL1
+ Jac[ 7] = (Ay[0]*globalPars[21]+Ay[1]*globalPars[22])+(Ay[2]*globalPars[23]-s2*T1); // dL1/dPhi
+ Jac[ 8] = (Ay[0]*globalPars[28]+Ay[1]*globalPars[29])+(Ay[2]*globalPars[30]-s3*T1); // dL1/dThe
+ Jac[ 9] =((Ay[0]*globalPars[35]+Ay[1]*globalPars[36])+(Ay[2]*globalPars[37]-s4*T1))*n; // dL1/dCM
+
+ double P3=0;
+ double P4=0;
+ /// transverse direction component
+ double C = globalPars[3]*globalPars[3]+globalPars[4]*globalPars[4];
+ if(C > 1.e-20) {
+ C= 1./C ;
+ /// unit direction vector in the X,Y plane
+ P3 = globalPars[3]*C;
+ P4 =globalPars[4]*C;
+ C =-sqrt(C);
+ }
+ else{
+ C=-1.e10;
+ P3 = 1.;
+ P4 =0.;
+ }
- // Jacobian production
- //
- Jac[ 0] = (Ax[0]*globalPars[ 7]+Ax[1]*globalPars[ 8])+(Ax[2]*globalPars[ 9]-s0*T0); // dL0/dL0
- Jac[ 1] = (Ax[0]*globalPars[14]+Ax[1]*globalPars[15])+(Ax[2]*globalPars[16]-s1*T0); // dL0/dL1
- Jac[ 2] = (Ax[0]*globalPars[21]+Ax[1]*globalPars[22])+(Ax[2]*globalPars[23]-s2*T0); // dL0/dPhi
- Jac[ 3] = (Ax[0]*globalPars[28]+Ax[1]*globalPars[29])+(Ax[2]*globalPars[30]-s3*T0); // dL0/dThe
- Jac[ 4] =((Ax[0]*globalPars[35]+Ax[1]*globalPars[36])+(Ax[2]*globalPars[37]-s4*T0))*n; // dL0/dCM
-
- Jac[ 5] = (Ay[0]*globalPars[ 7]+Ay[1]*globalPars[ 8])+(Ay[2]*globalPars[ 9]-s0*T1); // dL1/dL0
- Jac[ 6] = (Ay[0]*globalPars[14]+Ay[1]*globalPars[15])+(Ay[2]*globalPars[16]-s1*T1); // dL1/dL1
- Jac[ 7] = (Ay[0]*globalPars[21]+Ay[1]*globalPars[22])+(Ay[2]*globalPars[23]-s2*T1); // dL1/dPhi
- Jac[ 8] = (Ay[0]*globalPars[28]+Ay[1]*globalPars[29])+(Ay[2]*globalPars[30]-s3*T1); // dL1/dThe
- Jac[ 9] =((Ay[0]*globalPars[35]+Ay[1]*globalPars[36])+(Ay[2]*globalPars[37]-s4*T1))*n; // dL1/dCM
-
- double P3=0;
- double P4=0;
- /// transverse direction component
- double C = globalPars[3]*globalPars[3]+globalPars[4]*globalPars[4];
- if(C > 1.e-20) {
- C= 1./C ;
- /// unit direction vector in the X,Y plane
- P3 = globalPars[3]*C;
- P4 =globalPars[4]*C;
- C =-sqrt(C);
+ double T2 =(P3*globalPars[43]-P4*globalPars[42])*A;
+ double C44 = C*globalPars[44] *A;
+
+ Jac[10] = P3*globalPars[11]-P4*globalPars[10]-s0*T2; // dPhi/dL0
+ Jac[11] = P3*globalPars[18]-P4*globalPars[17]-s1*T2; // dPhi/dL1
+ Jac[12] = P3*globalPars[25]-P4*globalPars[24]-s2*T2; // dPhi/dPhi
+ Jac[13] = P3*globalPars[32]-P4*globalPars[31]-s3*T2; // dPhi/dThe
+ Jac[14] =(P3*globalPars[39]-P4*globalPars[38]-s4*T2)*n; // dPhi/dCM
+
+ Jac[15] = C*globalPars[12]-s0*C44; // dThe/dL0
+ Jac[16] = C*globalPars[19]-s1*C44; // dThe/dL1
+ Jac[17] = C*globalPars[26]-s2*C44; // dThe/dPhi
+ Jac[18] = C*globalPars[33]-s3*C44; // dThe/dThe
+ Jac[19] =(C*globalPars[40]-s4*C44)*n; // dThe/dCM
+ Jac[20] = 1.; // dCM /dCM
+
+ /// covariance matrix production using jacobian - CovNEW = J*CovOLD*Jt
+ AmgSymMatrix(5) newCov = Trk::PatternTrackParameters::newCovarianceMatrix(startingParameters.covariance(), Jac);
+ outputParameters.setParametersWithCovariance(m_surface, p, newCov);
+
+ /// check for negative diagonals in the cov
+ const double* t = &outputParameters.cov()[0];
+ if(t[0]<=0. || t[2]<=0. || t[5]<=0. || t[9]<=0. || t[14]<=0.) return false;
+ } else {
+ /// write into output parameters. Assign our surface to them
+ outputParameters.setParameters(m_surface,p);
}
- else{
- C=-1.e10;
- P3 = 1.;
- P4 =0.;
- }
-
- double T2 =(P3*globalPars[43]-P4*globalPars[42])*A;
- double C44 = C*globalPars[44] *A;
-
- Jac[10] = P3*globalPars[11]-P4*globalPars[10]-s0*T2; // dPhi/dL0
- Jac[11] = P3*globalPars[18]-P4*globalPars[17]-s1*T2; // dPhi/dL1
- Jac[12] = P3*globalPars[25]-P4*globalPars[24]-s2*T2; // dPhi/dPhi
- Jac[13] = P3*globalPars[32]-P4*globalPars[31]-s3*T2; // dPhi/dThe
- Jac[14] =(P3*globalPars[39]-P4*globalPars[38]-s4*T2)*n; // dPhi/dCM
-
- Jac[15] = C*globalPars[12]-s0*C44; // dThe/dL0
- Jac[16] = C*globalPars[19]-s1*C44; // dThe/dL1
- Jac[17] = C*globalPars[26]-s2*C44; // dThe/dPhi
- Jac[18] = C*globalPars[33]-s3*C44; // dThe/dThe
- Jac[19] =(C*globalPars[40]-s4*C44)*n; // dThe/dCM
- Jac[20] = 1.; // dCM /dCM
-
- /// covariance matrix production using jacobian - CovNEW = J*CovOLD*Jt
- outputParameters.newCovarianceMatrix(startingParameters,Jac);
- /// check for negative diagonals in the cov
- const double* t = &outputParameters.cov()[0];
- if(t[0]<=0. || t[2]<=0. || t[5]<=0. || t[9]<=0. || t[14]<=0.) return false;
+
return true;
}
diff --git a/Tracking/TrkEvent/TrkPatternParameters/TrkPatternParameters/PatternTrackParameters.h b/Tracking/TrkEvent/TrkPatternParameters/TrkPatternParameters/PatternTrackParameters.h
index fecd3e65279..b98f0883c85 100755
--- a/Tracking/TrkEvent/TrkPatternParameters/TrkPatternParameters/PatternTrackParameters.h
+++ b/Tracking/TrkEvent/TrkPatternParameters/TrkPatternParameters/PatternTrackParameters.h
@@ -74,6 +74,7 @@ namespace Trk {
void setParameters (const Surface*,const double* );
void setCovariance ( const double*);
void setParametersWithCovariance(const Surface*,const double*,const double*);
+ void setParametersWithCovariance(const Surface*,const double*,const AmgSymMatrix(5)&);
///////////////////////////////////////////////////////////////////
// Convertors
@@ -101,8 +102,7 @@ namespace Trk {
// Covariance matrix production using jacobian CovNEW = J*CovOLD*Jt
///////////////////////////////////////////////////////////////////
- void newCovarianceMatrix
- (PatternTrackParameters&,double*);
+ static AmgSymMatrix(5) newCovarianceMatrix(const AmgSymMatrix(5) &, const double *);
///////////////////////////////////////////////////////////////////
// Print
@@ -261,6 +261,19 @@ namespace Trk {
setParameters(s,p);
setCovariance(c );
}
+
+ inline void PatternTrackParameters::setParametersWithCovariance
+ (const Surface* s,const double* p,const AmgSymMatrix(5)& c)
+ {
+ double C[15] = {
+ c(0, 0),
+ c(1, 0), c(1, 1),
+ c(2, 0), c(2, 1), c(2, 2),
+ c(3, 0), c(3, 1), c(3, 2), c(3, 3),
+ c(4, 0), c(4, 1), c(4, 2), c(4, 3), c(4, 4)
+ };
+ setParametersWithCovariance(s, p, C);
+ }
///////////////////////////////////////////////////////////////////
// Diagonal symetric matrix production
diff --git a/Tracking/TrkEvent/TrkPatternParameters/src/PatternTrackParameters.cxx b/Tracking/TrkEvent/TrkPatternParameters/src/PatternTrackParameters.cxx
index ef8e2b005d0..5d38864bf1c 100755
--- a/Tracking/TrkEvent/TrkPatternParameters/src/PatternTrackParameters.cxx
+++ b/Tracking/TrkEvent/TrkPatternParameters/src/PatternTrackParameters.cxx
@@ -81,60 +81,63 @@ bool Trk::PatternTrackParameters::production(const Trk::ParametersBase<5,Trk::Ch
// Covariance matrix production using jacobian cov(this) = J*( Tp cov)*Jt
///////////////////////////////////////////////////////////////////
-void Trk::PatternTrackParameters::newCovarianceMatrix
-(Trk::PatternTrackParameters& Tp,double* Jac)
+
+AmgSymMatrix(5) Trk::PatternTrackParameters::newCovarianceMatrix
+(const AmgSymMatrix(5) & V, const double * Jac)
{
- double* V = &Tp.m_covariance[0];
- double a11 = (Jac[ 0]*V[ 0]+Jac[ 1]*V[ 1]+Jac[ 2]*V[ 3])+(Jac[ 3]*V[ 6]+Jac[ 4]*V[10]);
- double a12 = (Jac[ 0]*V[ 1]+Jac[ 1]*V[ 2]+Jac[ 2]*V[ 4])+(Jac[ 3]*V[ 7]+Jac[ 4]*V[11]);
- double a13 = (Jac[ 0]*V[ 3]+Jac[ 1]*V[ 4]+Jac[ 2]*V[ 5])+(Jac[ 3]*V[ 8]+Jac[ 4]*V[12]);
- double a14 = (Jac[ 0]*V[ 6]+Jac[ 1]*V[ 7]+Jac[ 2]*V[ 8])+(Jac[ 3]*V[ 9]+Jac[ 4]*V[13]);
- double a15 = (Jac[ 0]*V[10]+Jac[ 1]*V[11]+Jac[ 2]*V[12])+(Jac[ 3]*V[13]+Jac[ 4]*V[14]);
-
- m_covariance[ 0] = (a11*Jac[ 0]+a12*Jac[ 1]+a13*Jac[ 2])+(a14*Jac[ 3]+a15*Jac[ 4]);
-
- double a21 = (Jac[ 5]*V[ 0]+Jac[ 6]*V[ 1]+Jac[ 7]*V[ 3])+(Jac[ 8]*V[ 6]+Jac[ 9]*V[10]);
- double a22 = (Jac[ 5]*V[ 1]+Jac[ 6]*V[ 2]+Jac[ 7]*V[ 4])+(Jac[ 8]*V[ 7]+Jac[ 9]*V[11]);
- double a23 = (Jac[ 5]*V[ 3]+Jac[ 6]*V[ 4]+Jac[ 7]*V[ 5])+(Jac[ 8]*V[ 8]+Jac[ 9]*V[12]);
- double a24 = (Jac[ 5]*V[ 6]+Jac[ 6]*V[ 7]+Jac[ 7]*V[ 8])+(Jac[ 8]*V[ 9]+Jac[ 9]*V[13]);
- double a25 = (Jac[ 5]*V[10]+Jac[ 6]*V[11]+Jac[ 7]*V[12])+(Jac[ 8]*V[13]+Jac[ 9]*V[14]);
-
- m_covariance[ 1] = (a21*Jac[ 0]+a22*Jac[ 1]+a23*Jac[ 2])+(a24*Jac[ 3]+a25*Jac[ 4]);
- m_covariance[ 2] = (a21*Jac[ 5]+a22*Jac[ 6]+a23*Jac[ 7])+(a24*Jac[ 8]+a25*Jac[ 9]);
-
- double a31 = (Jac[10]*V[ 0]+Jac[11]*V[ 1]+Jac[12]*V[ 3])+(Jac[13]*V[ 6]+Jac[14]*V[10]);
- double a32 = (Jac[10]*V[ 1]+Jac[11]*V[ 2]+Jac[12]*V[ 4])+(Jac[13]*V[ 7]+Jac[14]*V[11]);
- double a33 = (Jac[10]*V[ 3]+Jac[11]*V[ 4]+Jac[12]*V[ 5])+(Jac[13]*V[ 8]+Jac[14]*V[12]);
- double a34 = (Jac[10]*V[ 6]+Jac[11]*V[ 7]+Jac[12]*V[ 8])+(Jac[13]*V[ 9]+Jac[14]*V[13]);
- double a35 = (Jac[10]*V[10]+Jac[11]*V[11]+Jac[12]*V[12])+(Jac[13]*V[13]+Jac[14]*V[14]);
-
- m_covariance[ 3] = (a31*Jac[ 0]+a32*Jac[ 1]+a33*Jac[ 2])+(a34*Jac[ 3]+a35*Jac[ 4]);
- m_covariance[ 4] = (a31*Jac[ 5]+a32*Jac[ 6]+a33*Jac[ 7])+(a34*Jac[ 8]+a35*Jac[ 9]);
- m_covariance[ 5] = (a31*Jac[10]+a32*Jac[11]+a33*Jac[12])+(a34*Jac[13]+a35*Jac[14]);
-
- double a41 = (Jac[15]*V[ 0]+Jac[16]*V[ 1]+Jac[17]*V[ 3])+(Jac[18]*V[ 6]+Jac[19]*V[10]);
- double a42 = (Jac[15]*V[ 1]+Jac[16]*V[ 2]+Jac[17]*V[ 4])+(Jac[18]*V[ 7]+Jac[19]*V[11]);
- double a43 = (Jac[15]*V[ 3]+Jac[16]*V[ 4]+Jac[17]*V[ 5])+(Jac[18]*V[ 8]+Jac[19]*V[12]);
- double a44 = (Jac[15]*V[ 6]+Jac[16]*V[ 7]+Jac[17]*V[ 8])+(Jac[18]*V[ 9]+Jac[19]*V[13]);
- double a45 = (Jac[15]*V[10]+Jac[16]*V[11]+Jac[17]*V[12])+(Jac[18]*V[13]+Jac[19]*V[14]);
-
- m_covariance[ 6] = (a41*Jac[ 0]+a42*Jac[ 1]+a43*Jac[ 2])+(a44*Jac[ 3]+a45*Jac[ 4]);
- m_covariance[ 7] = (a41*Jac[ 5]+a42*Jac[ 6]+a43*Jac[ 7])+(a44*Jac[ 8]+a45*Jac[ 9]);
- m_covariance[ 8] = (a41*Jac[10]+a42*Jac[11]+a43*Jac[12])+(a44*Jac[13]+a45*Jac[14]);
- m_covariance[ 9] = (a41*Jac[15]+a42*Jac[16]+a43*Jac[17])+(a44*Jac[18]+a45*Jac[19]);
-
- double a51 = Jac[20]*V[10];
- double a52 = Jac[20]*V[11];
- double a53 = Jac[20]*V[12];
- double a54 = Jac[20]*V[13];
- double a55 = Jac[20]*V[14];
-
- m_covariance[10] = (a51*Jac[ 0]+a52*Jac[ 1]+a53*Jac[ 2])+(a54*Jac[ 3]+a55*Jac[ 4]);
- m_covariance[11] = (a51*Jac[ 5]+a52*Jac[ 6]+a53*Jac[ 7])+(a54*Jac[ 8]+a55*Jac[ 9]);
- m_covariance[12] = (a51*Jac[10]+a52*Jac[11]+a53*Jac[12])+(a54*Jac[13]+a55*Jac[14]);
- m_covariance[13] = (a51*Jac[15]+a52*Jac[16]+a53*Jac[17])+(a54*Jac[18]+a55*Jac[19]);
- m_covariance[14] = a55*Jac[20];
- m_iscovariance = true;
+ AmgSymMatrix(5) rv;
+
+ double a11 = (Jac[ 0]*V(0, 0)+Jac[ 1]*V(0, 1)+Jac[ 2]*V(0, 2))+(Jac[ 3]*V(0, 3)+Jac[ 4]*V(0, 4));
+ double a12 = (Jac[ 0]*V(0, 1)+Jac[ 1]*V(1, 1)+Jac[ 2]*V(1, 2))+(Jac[ 3]*V(1, 3)+Jac[ 4]*V(1, 4));
+ double a13 = (Jac[ 0]*V(0, 2)+Jac[ 1]*V(1, 2)+Jac[ 2]*V(2, 2))+(Jac[ 3]*V(2, 3)+Jac[ 4]*V(2, 4));
+ double a14 = (Jac[ 0]*V(0, 3)+Jac[ 1]*V(1, 3)+Jac[ 2]*V(2, 3))+(Jac[ 3]*V(3, 3)+Jac[ 4]*V(3, 4));
+ double a15 = (Jac[ 0]*V(0, 4)+Jac[ 1]*V(1, 4)+Jac[ 2]*V(2, 4))+(Jac[ 3]*V(3, 4)+Jac[ 4]*V(4, 4));
+
+ rv.fillSymmetric(0, 0, (a11*Jac[ 0]+a12*Jac[ 1]+a13*Jac[ 2])+(a14*Jac[ 3]+a15*Jac[ 4]));
+
+ double a21 = (Jac[ 5]*V(0, 0)+Jac[ 6]*V(0, 1)+Jac[ 7]*V(0, 2))+(Jac[ 8]*V(0, 3)+Jac[ 9]*V(0, 4));
+ double a22 = (Jac[ 5]*V(0, 1)+Jac[ 6]*V(1, 1)+Jac[ 7]*V(1, 2))+(Jac[ 8]*V(1, 3)+Jac[ 9]*V(1, 4));
+ double a23 = (Jac[ 5]*V(0, 2)+Jac[ 6]*V(1, 2)+Jac[ 7]*V(2, 2))+(Jac[ 8]*V(2, 3)+Jac[ 9]*V(2, 4));
+ double a24 = (Jac[ 5]*V(0, 3)+Jac[ 6]*V(1, 3)+Jac[ 7]*V(2, 3))+(Jac[ 8]*V(3, 3)+Jac[ 9]*V(3, 4));
+ double a25 = (Jac[ 5]*V(0, 4)+Jac[ 6]*V(1, 4)+Jac[ 7]*V(2, 4))+(Jac[ 8]*V(3, 4)+Jac[ 9]*V(4, 4));
+
+ rv.fillSymmetric(1, 0, (a21*Jac[ 0]+a22*Jac[ 1]+a23*Jac[ 2])+(a24*Jac[ 3]+a25*Jac[ 4]));
+ rv.fillSymmetric(1, 1, (a21*Jac[ 5]+a22*Jac[ 6]+a23*Jac[ 7])+(a24*Jac[ 8]+a25*Jac[ 9]));
+
+ double a31 = (Jac[10]*V(0, 0)+Jac[11]*V(0, 1)+Jac[12]*V(0, 2))+(Jac[13]*V(0, 3)+Jac[14]*V(0, 4));
+ double a32 = (Jac[10]*V(0, 1)+Jac[11]*V(1, 1)+Jac[12]*V(1, 2))+(Jac[13]*V(1, 3)+Jac[14]*V(1, 4));
+ double a33 = (Jac[10]*V(0, 2)+Jac[11]*V(1, 2)+Jac[12]*V(2, 2))+(Jac[13]*V(2, 3)+Jac[14]*V(2, 4));
+ double a34 = (Jac[10]*V(0, 3)+Jac[11]*V(1, 3)+Jac[12]*V(2, 3))+(Jac[13]*V(3, 3)+Jac[14]*V(3, 4));
+ double a35 = (Jac[10]*V(0, 4)+Jac[11]*V(1, 4)+Jac[12]*V(2, 4))+(Jac[13]*V(3, 4)+Jac[14]*V(4, 4));
+
+ rv.fillSymmetric(2, 0, (a31*Jac[ 0]+a32*Jac[ 1]+a33*Jac[ 2])+(a34*Jac[ 3]+a35*Jac[ 4]));
+ rv.fillSymmetric(2, 1, (a31*Jac[ 5]+a32*Jac[ 6]+a33*Jac[ 7])+(a34*Jac[ 8]+a35*Jac[ 9]));
+ rv.fillSymmetric(2, 2, (a31*Jac[10]+a32*Jac[11]+a33*Jac[12])+(a34*Jac[13]+a35*Jac[14]));
+
+ double a41 = (Jac[15]*V(0, 0)+Jac[16]*V(0, 1)+Jac[17]*V(0, 2))+(Jac[18]*V(0, 3)+Jac[19]*V(0, 4));
+ double a42 = (Jac[15]*V(0, 1)+Jac[16]*V(1, 1)+Jac[17]*V(1, 2))+(Jac[18]*V(1, 3)+Jac[19]*V(1, 4));
+ double a43 = (Jac[15]*V(0, 2)+Jac[16]*V(1, 2)+Jac[17]*V(2, 2))+(Jac[18]*V(2, 3)+Jac[19]*V(2, 4));
+ double a44 = (Jac[15]*V(0, 3)+Jac[16]*V(1, 3)+Jac[17]*V(2, 3))+(Jac[18]*V(3, 3)+Jac[19]*V(3, 4));
+ double a45 = (Jac[15]*V(0, 4)+Jac[16]*V(1, 4)+Jac[17]*V(2, 4))+(Jac[18]*V(3, 4)+Jac[19]*V(4, 4));
+
+ rv.fillSymmetric(3, 0, (a41*Jac[ 0]+a42*Jac[ 1]+a43*Jac[ 2])+(a44*Jac[ 3]+a45*Jac[ 4]));
+ rv.fillSymmetric(3, 1, (a41*Jac[ 5]+a42*Jac[ 6]+a43*Jac[ 7])+(a44*Jac[ 8]+a45*Jac[ 9]));
+ rv.fillSymmetric(3, 2, (a41*Jac[10]+a42*Jac[11]+a43*Jac[12])+(a44*Jac[13]+a45*Jac[14]));
+ rv.fillSymmetric(3, 3, (a41*Jac[15]+a42*Jac[16]+a43*Jac[17])+(a44*Jac[18]+a45*Jac[19]));
+
+ double a51 = Jac[20]*V(0, 4);
+ double a52 = Jac[20]*V(1, 4);
+ double a53 = Jac[20]*V(2, 4);
+ double a54 = Jac[20]*V(3, 4);
+ double a55 = Jac[20]*V(4, 4);
+
+ rv.fillSymmetric(4, 0, (a51*Jac[ 0]+a52*Jac[ 1]+a53*Jac[ 2])+(a54*Jac[ 3]+a55*Jac[ 4]));
+ rv.fillSymmetric(4, 1, (a51*Jac[ 5]+a52*Jac[ 6]+a53*Jac[ 7])+(a54*Jac[ 8]+a55*Jac[ 9]));
+ rv.fillSymmetric(4, 2, (a51*Jac[10]+a52*Jac[11]+a53*Jac[12])+(a54*Jac[13]+a55*Jac[14]));
+ rv.fillSymmetric(4, 3, (a51*Jac[15]+a52*Jac[16]+a53*Jac[17])+(a54*Jac[18]+a55*Jac[19]));
+ rv.fillSymmetric(4, 4, a55*Jac[20]);
+
+ return rv;
}
///////////////////////////////////////////////////////////////////
diff --git a/Tracking/TrkExtrapolation/TrkExRungeKuttaPropagator/src/RungeKuttaPropagator.cxx b/Tracking/TrkExtrapolation/TrkExRungeKuttaPropagator/src/RungeKuttaPropagator.cxx
index 969de3df8ab..4ce1808b1b2 100755
--- a/Tracking/TrkExtrapolation/TrkExRungeKuttaPropagator/src/RungeKuttaPropagator.cxx
+++ b/Tracking/TrkExtrapolation/TrkExRungeKuttaPropagator/src/RungeKuttaPropagator.cxx
@@ -1499,11 +1499,13 @@ bool Trk::RungeKuttaPropagator::propagateRungeKutta
// New simple track parameters production
//
- Tb.setParameters(&Su,p);
if(useJac) {
- Tb.newCovarianceMatrix(Ta,Jac);
+ AmgSymMatrix(5) newCov = Trk::PatternTrackParameters::newCovarianceMatrix(Ta.covariance(), Jac);
+ Tb.setParametersWithCovariance(&Su, p, newCov);
const double* cv = Tb.cov();
if( cv[0]<=0. || cv[2]<=0. || cv[5]<=0. || cv[9]<=0. || cv[14]<=0.) return false;
+ } else {
+ Tb.setParameters(&Su,p);
}
return true;
}
--
GitLab