From 4af3ff12f89da4472ffdb733c18fb1a77f74f84a Mon Sep 17 00:00:00 2001
From: christos <christos@cern.ch>
Date: Mon, 13 Jul 2020 17:11:20 +0100
Subject: [PATCH] review comments
---
.../src/JetFitterInitializationHelper.cxx | 8 ++--
.../src/KalmanVertexOnJetAxisUpdator.cxx | 41 ++++---------------
.../TrkVKalVrtCore/src/Matrix.cxx | 16 ++++----
3 files changed, 20 insertions(+), 45 deletions(-)
diff --git a/Tracking/TrkVertexFitter/TrkJetVxFitter/src/JetFitterInitializationHelper.cxx b/Tracking/TrkVertexFitter/TrkJetVxFitter/src/JetFitterInitializationHelper.cxx
index 2bad68d03d9..8232277534f 100755
--- a/Tracking/TrkVertexFitter/TrkJetVxFitter/src/JetFitterInitializationHelper.cxx
+++ b/Tracking/TrkVertexFitter/TrkJetVxFitter/src/JetFitterInitializationHelper.cxx
@@ -41,7 +41,7 @@ namespace Trk
Amg::Vector3D getSingleVtxPositionWithSignFlip(const Amg::VectorX & myPosition,
int numVertex,
- bool signfliptreatment) {
+ bool signFlipTreatment) {
int numbRow=numRow(numVertex);
double xv=myPosition[Trk::jet_xv];
@@ -56,7 +56,7 @@ namespace Trk
dist=dist/fabs(dist)*300./sin(theta);
}
if (dist<0) {
- if (signfliptreatment) {
+ if (signFlipTreatment) {
dist=-dist;
} else {
dist=0.;
@@ -341,7 +341,7 @@ namespace Trk
}
void JetFitterInitializationHelper::linearizeAllTracks(VxJetCandidate* myJetCandidate,
- bool signfliptreatment,
+ bool signFlipTreatment,
double maxdistance) const {
const VertexPositions & myLinVertexPosition=myJetCandidate->getLinearizationVertexPositions();
@@ -385,7 +385,7 @@ namespace Trk
for (std::vector<VxVertexOnJetAxis*>::const_iterator VtxIter=VtxBegin;VtxIter!=VtxEnd;++VtxIter) {
int numVertex=(*VtxIter)->getNumVertex();
- Amg::Vector3D secondaryVertexPos(getSingleVtxPositionWithSignFlip(myPosition,numVertex,signfliptreatment));
+ Amg::Vector3D secondaryVertexPos(getSingleVtxPositionWithSignFlip(myPosition,numVertex,signFlipTreatment));
// std::cout << " Considering linearization at n. vertex " << numVertex << " pos " << secondaryVertexPos << std::endl;
diff --git a/Tracking/TrkVertexFitter/TrkJetVxFitter/src/KalmanVertexOnJetAxisUpdator.cxx b/Tracking/TrkVertexFitter/TrkJetVxFitter/src/KalmanVertexOnJetAxisUpdator.cxx
index 06cf0df3a00..e4510cc4815 100755
--- a/Tracking/TrkVertexFitter/TrkJetVxFitter/src/KalmanVertexOnJetAxisUpdator.cxx
+++ b/Tracking/TrkVertexFitter/TrkJetVxFitter/src/KalmanVertexOnJetAxisUpdator.cxx
@@ -256,12 +256,7 @@ namespace Trk{
} else {
ndf+=sign*1.;
}
-
- //add the chi2 to the previous chi2
- //Trk::RecVertexPositions r_vtx(new_vrt_weight_times_position,erm,ndf,0.,true);//ndf, chi2);
- Trk::RecVertexPositions r_vtx(new_vrt_pos,new_vrt_cov,ndf,chi2);
-
- return r_vtx;
+ return Trk::RecVertexPositions(new_vrt_pos,new_vrt_cov,ndf,chi2);
}
@@ -278,8 +273,7 @@ namespace Trk{
if(myPosition.covariancePosition().determinant() ==0.0) {
ATH_MSG_WARNING ("The vertex-positions covariance matrix is not invertible");
ATH_MSG_WARNING ("The copy of initial vertex returned");
- const Trk::RecVertexPositions& r_vtx(myPosition);
- return r_vtx;
+ return Trk::RecVertexPositions(myPosition);
}
#ifdef KalmanVertexUpdate_OLD
const Amg::VectorX & old_vrt_pos = myPosition.position();
@@ -299,19 +293,10 @@ namespace Trk{
ATH_MSG_WARNING(" The determinant of the track covariance matrix is zero or negative: " << trackParametersWeight.determinant());
}
-
-// if (old_full_vrt_weight.determinant()<=0)
-// {
-// ATH_MSG_WARNING(" The determinant of the vertex full weight matrix is zero or negative: " << old_full_vrt_weight.determinant());
-// }
-
-
-// std::cout << " WEIGHT MATRIX of POS " << old_full_vrt_weight << std::endl;
Amg::MatrixX old_vrt_weight(numrow_toupdate+1,numrow_toupdate+1);
old_vrt_weight = old_full_vrt_weight.block(0,0,numrow_toupdate+1,numrow_toupdate+1);
-// std::cout << " Reduced Weight of POS " << old_vrt_weight << std::endl;
#endif
@@ -320,16 +305,14 @@ namespace Trk{
AmgSymMatrix(3) S = B.transpose()*(trackParametersWeight*B);
-// std::cout << " Matrix S " << S << std::endl;
if(S.determinant() ==0.0) {
ATH_MSG_WARNING ("The S matrix is not invertible");
ATH_MSG_WARNING ("A copy of initial vertex returned");
- const Trk::RecVertexPositions& r_vtx(myPosition);
- return r_vtx;
- } S = S.inverse().eval();
+ return Trk::RecVertexPositions(myPosition);
+ }
+ S = S.inverse().eval();
-// std::cout << " Matrix S after inversion " << S << std::endl;
//G_b = G_k - G_k*B_k*W_k*B_k^(T)*G_k
@@ -337,7 +320,6 @@ namespace Trk{
//gB=gB/2.0+gB.transpose().eval()/2.0;
-// std::cout << " gB " << gB << std::endl;
#ifdef Updator_DEBUG
std::cout<<"Gain factor obtained: "<<trackParametersWeight*(B*(S*B.transpose()))*trackParametersWeight.transpose()<<std::endl;
@@ -592,7 +574,8 @@ namespace Trk{
if(S.determinant() ==0.0) {
ATH_MSG_ERROR ("The matrix S is not invertible, return chi2 0");
return -0.;
- } S = S.inverse().eval();
+ }
+ S = S.inverse().eval();
//refitted track momentum
Amg::Vector3D newTrackMomentum = S*B.transpose()*trackParametersWeight*(trackParameters - constantTerm - A*new_vrt_position);
@@ -708,7 +691,7 @@ namespace Trk{
int numRows=new_vrt_weight.rows();
if (numRows<=6) {
if(new_vrt_weight.determinant() ==0)
-// if |A|==0; then A is not invertible;
+ // if |A|==0; then A is not invertible;
{
throw std::string("The reduced weight matrix is not invertible. Previous vertex returned ");
}
@@ -726,13 +709,10 @@ namespace Trk{
B(i,j-5)=new_vrt_weight(i,j);
}
}
-// std::cout << " weight matrix: " << new_vrt_weight << std::endl;
-// std::cout << " non diag upper right elem: " << B << std::endl;
// CLHEP::HepSymMatrix D=new_vrt_weight.sub(6,numRows);
Amg::MatrixX D = new_vrt_weight.block(5,5,numRows-5,numRows-5); // will be of dim (numrows-5)
- //std::cout << " D "<<std::endl << D << std::endl;
// Eigen::DiagonalMatrix<double,Eigen::Dynamic> DdiagINV(D.rows()); // was CLHEP::HepDiagMatrix
// if (numRows>6 && D.isDiagonal(/*precision<scalar>=*/1e-7))
@@ -767,15 +747,12 @@ namespace Trk{
}
// D=D.inverse().eval();
-// std::cout << " D after inversion " << D << std::endl;
-
Amg::MatrixX E = A - B*(D*B.transpose());
if (E.determinant() == 0.) {
throw std::string("Cannot invert E matrix...");
}
E=E.inverse().eval();
-// std::cout << " E " << E << std::endl;
Amg::MatrixX finalWeight(numRows,numRows);
finalWeight.setZero();
@@ -786,7 +763,6 @@ namespace Trk{
finalWeight.block(5,5,D.rows(),D.rows()) =
D+(D*((B.transpose()*(E*B))*D.transpose()));
-// std::cout << " finalWeight " << finalWeight << std::endl;
// Amg::MatrixX normalInversion=new_vrt_weight.inverse().eval();
// Amg::MatrixX smartInversion=finalWeight;
@@ -804,7 +780,6 @@ namespace Trk{
}
}
}
- if (mismatch) std::cout << " mismatch, normalInv: " << normalInversion << " smartInv " << smartInversion << " det 1 " << normalInversion.determinant() << " det 2 " << smartInversion.determinant() << std::endl;
*/
// new_vrt_weight = new_vrt_weight.inverse().eval();//finalWeight;//MODIFIED!!
diff --git a/Tracking/TrkVertexFitter/TrkVKalVrtCore/src/Matrix.cxx b/Tracking/TrkVertexFitter/TrkVKalVrtCore/src/Matrix.cxx
index fc842b46292..59126554234 100644
--- a/Tracking/TrkVertexFitter/TrkVKalVrtCore/src/Matrix.cxx
+++ b/Tracking/TrkVertexFitter/TrkVKalVrtCore/src/Matrix.cxx
@@ -270,9 +270,9 @@ void scaleg(double *g, double *scale, long int N, long int mfirst)
for (i__ = 1; i__ <= N; ++i__) {
scale[i__] = 1.;
if (g_ref(i__, i__) == 0.) continue;
- tmp = sqrt(fabs(g_ref(i__, i__)));
+ tmp = std::sqrt(fabs(g_ref(i__, i__)));
//scale[i__] = 1./tmp; g_ref(i__, i__) = d_sign( 1., g_ref(i__, i__)); //VK old version -> diag==1
- scale[i__] = 1./sqrt(tmp); g_ref(i__, i__) = d_sign( tmp, g_ref(i__, i__)); //VK new version -> diag=sqrt(diag_old)
+ scale[i__] = 1./std::sqrt(tmp); g_ref(i__, i__) = d_sign( tmp, g_ref(i__, i__)); //VK new version -> diag=sqrt(diag_old)
}
if (N <= 1) return;
@@ -414,8 +414,8 @@ double vkPythag(double a, double b)
double absa,absb;
absa=fabs(a);
absb=fabs(b);
- if (absa > absb) return absa*sqrt(1.0+(absb/absa)*absb/absa);
- return (absb == 0.0 ? 0.0 : absb*sqrt(1.0+(absa/absb)*absa/absb));
+ if (absa > absb) return absa*std::sqrt(1.0+(absb/absa)*absb/absa);
+ return (absb == 0.0 ? 0.0 : absb*std::sqrt(1.0+(absa/absb)*absa/absb));
}
void vkSVDCmp(double **a, int m, int n, double w[], double **v)
@@ -437,7 +437,7 @@ void vkSVDCmp(double **a, int m, int n, double w[], double **v)
s += a[k][i]*a[k][i];
}
f=a[i][i];
- g = -SIGN(sqrt(s),f);
+ g = -SIGN(std::sqrt(s),f);
h=f*g-s;
a[i][i]=f-g;
for (j=l;j<=n;j++) {
@@ -458,7 +458,7 @@ void vkSVDCmp(double **a, int m, int n, double w[], double **v)
s += a[i][k]*a[i][k];
}
f=a[i][l];
- g = -SIGN(sqrt(s),f);
+ g = -SIGN(std::sqrt(s),f);
h=f*g-s;
a[i][l]=f-g;
for (k=l;k<=n;k++) rv1[k]=a[i][k]/h;
@@ -742,10 +742,10 @@ int vkjacobi(double **a, int n, double d[], double **v)
t=(a[ip][iq])/h;
else {
theta=0.5*h/(a[ip][iq]);
- t=1.0/(fabs(theta)+sqrt(1.0+theta*theta));
+ t=1.0/(fabs(theta)+std::sqrt(1.0+theta*theta));
if (theta < 0.0) t = -t;
}
- c=1.0/sqrt(1+t*t);
+ c=1.0/std::sqrt(1+t*t);
s=t*c;
tau=s/(1.0+c);
h=t*a[ip][iq];
--
GitLab