Commit 0c10b6c3 authored by Vadim Kostyukhin's avatar Vadim Kostyukhin 😠 Committed by Graeme Stewart
Browse files

Better protection against small eigenvalue in track error matrix in CFit.cxx...

Better protection against small eigenvalue in track error matrix in CFit.cxx (TrkVKalVrtCore-00-03-63)

	* New return codes
	* Add to diagonal if eigenvalue of track error matrix is <10^-15 in CFit.cxx
	* TrkVKalVrtCore-00-03-63

2016-09-20  Vadim Kostyukhin <Vadim.Kostyukhin@cern.ch>

	* Correct error in cfTrkCovarCorr (correlations were always downscaled!)
	* Replace 3.14159.. by M_PI
	* TrkVKalVrtCore-00-03-62

2016-09-19  scott snyder  <snyder@bnl.gov>

	* Tagging TrkVKalVrtCore-00-03-61.
	* Fix most gcc6 misleading indentation warnings.


Former-commit-id: d378f733
parent 1487638d
......@@ -269,7 +269,7 @@ long int fitVertex(VKVertex * vk, long int iflag)
savedExtrapVertices[it-1].Y=dxyzst[1];
savedExtrapVertices[it-1].Z=dxyzst[2];
if(it==1){
if(!myPropagator.checkTarget(dxyzst)) { IERR=-6; return IERR; } // First guess is definitely outside working volume
if(!myPropagator.checkTarget(dxyzst)) { IERR=-10; return IERR; } // First guess is definitely outside working volume
}
/* ---------------------------------------------------------------- */
/* Propagate parameters and errors at point dXYZST */
......@@ -291,19 +291,21 @@ long int fitVertex(VKVertex * vk, long int iflag)
double ddx=savedExtrapVertices[it-1].X - oldX;
double ddy=savedExtrapVertices[it-1].Y - oldY;
double ddz=savedExtrapVertices[it-1].Z - oldZ;
if( sqrt(ddx*ddx+ddy*ddy+ddz*ddz)<5.) { IERR=-6; return IERR; } // Impossible to extrapolate
if( sqrt(ddx*ddx+ddy*ddy+ddz*ddz)<5.) { IERR=-11; return IERR; } // Impossible to extrapolate
}
double targV[3]={dxyzst[0],dxyzst[1],dxyzst[2]};
for (tk = 0; tk < NTRK; ++tk) {
myPropagator.Propagate(vk->TrackList[tk], vk->refV, targV, tmpPer, tmpCov);
cfTrkCovarCorr(tmpCov);
double eig5=cfSmallEigenvalue(tmpCov,5 );
if(eig5<1.e-15 || tmpCov[0]>1.e9) { //Bad propagation with material. Try without it.
if(eig5<1.e-15 ){
tmpCov[0]+=1.e-15; tmpCov[2]+=1.e-15; tmpCov[5]+=1.e-15; tmpCov[9]+=1.e-15; tmpCov[14]+=1.e-15;
}else if(tmpCov[0]>1.e9) { //Bad propagation with material. Try without it.
myPropagator.Propagate(-999, vk->TrackList[tk]->Charge,
vk->TrackList[tk]->refPerig,vk->TrackList[tk]->refCovar,
vk->refV, targV, tmpPer, tmpCov);
if(cfSmallEigenvalue(tmpCov,5 )<1.e-15){ //Final protection
tmpCov[1]=0.;tmpCov[3]=0.;tmpCov[6]=0.;tmpCov[10]=0.;
tmpCov[1]=0.;tmpCov[3]=0.;tmpCov[6]=0.;tmpCov[10]=0.;
tmpCov[4]=0.;tmpCov[7]=0.;tmpCov[11]=0.;
tmpCov[8]=0.;tmpCov[12]=0.;
tmpCov[13]=0.;
......
......@@ -400,7 +400,8 @@ int fitVertexCascade( VKVertex * vk, int Pointing)
}
if( (Chi2Cur>1.e8 && Iter>0) || IERR){
delete[] fullMatrix; delete[] fullLSide; delete[] tmpLSide; delete[] iniCovMatrix;
if(Chi2Cur>1.e8) return -1; return IERR;
if(Chi2Cur>1.e8) return -1;
return IERR;
}
//
//--- Attempt to dump oscillations by returning to good position
......@@ -786,7 +787,8 @@ void getFittedCascade( std::vector< Vect3DF > & cVertices,
cascadeCovar.push_back( transformCovar(DIM, Deriv, cascadeCovarFit[iv] ) ); //Transform covariance and save it
for(int kk=0; kk<6; kk++) tmpCov[kk]=cascadeCovar[iv][kk];
covVertices.push_back(tmpCov);
for(it=0; it<DIM; it++) delete[]Deriv[it]; delete[]Deriv;
for(it=0; it<DIM; it++) delete[]Deriv[it];
delete[]Deriv;
}
//
// Get full physics covariance
......@@ -799,7 +801,8 @@ void getFittedCascade( std::vector< Vect3DF > & cVertices,
fullCovar[it*(it+1)/2+jt]=tmp; // symmetrical index formula works ONLY if it>=jt!!!
}
}
for(it=0; it<PDIM; it++) delete[]DPhys[it]; delete[]DPhys;
for(it=0; it<PDIM; it++) delete[]DPhys[it];
delete[]DPhys;
//
cleanCascade();
return;
......
......@@ -51,7 +51,8 @@ void calcPointConstraint( VKPointConstraint * cnst )
double mom= sqrt(ptot[2]*ptot[2] + Pt*Pt) ;
double phi=atan2(ptot[1],ptot[0]);
double theta=acos( ptot[2]/mom );
if(theta<1.e-5) theta=1.e-5; if(theta>3.141592653-1.e-5) theta=3.141592653-1.e-5;
if(theta<1.e-5) theta=1.e-5;
if(theta>M_PI-1.e-5) theta=M_PI-1.e-5;
double invR = magConst / Pt; if(Charge) invR *= Charge;
//
// Fitted vertex position in global reference frame
......
......@@ -23,7 +23,6 @@ void calcPhiConstraint( VKPhiConstraint * cnst)
VKVertex * vk=cnst->getOriginVertex();
int NTRK = vk->TrackList.size();
int i,j;
double cPi=3.141592653589;
double Scale=1.; // Scaling for better precision VK 28.03.2011 Wrong for error matrix!!! Should always be 1.!!!
double diff=0., aa=0;
......@@ -31,8 +30,8 @@ void calcPhiConstraint( VKPhiConstraint * cnst)
for(i=0; i<NTRK-1; i++){
for(j=i+1; j<NTRK; j++){
diff = vk->TrackList[i]->cnstP[1] - vk->TrackList[j]->cnstP[1];
while(diff > cPi) diff -= 2.*cPi;
while(diff < -cPi) diff += 2.*cPi;
while(diff > M_PI) diff -= 2.*M_PI;
while(diff < -M_PI) diff += 2.*M_PI;
aa += diff*Scale;
deriv[i] += Scale;
deriv[j] -= Scale;
......
......@@ -66,7 +66,7 @@ void robtest(VKVertex * vk, long int ifl)
return;
}
/* -- */
double halfPi=3.141592653589/2.;
double halfPi=M_PI/2.;
for (it = 0; it < NTRK; ++it) {
VKTrack *trk=vk->TrackList[it];
if(trk->Id < 0) continue; // Not a real track
......
......@@ -39,8 +39,8 @@ double cfchi2(double *xyzt, long int *ich, double *part,
d3 = par0[2] - part[1];
d4 = par0[3] - phif;
d5 = par0[4] - part[3];
while(d4 > 3.141592653589)d4-=2.*3.141592653589;
while(d4 < -3.141592653589)d4+=2.*3.141592653589;
while(d4 > M_PI)d4-=2.*M_PI;
while(d4 < -M_PI)d4+=2.*M_PI;
// -----------------------Check of propagation
// double paro[5],parn[5],s,ref[3],peri[3];
// paro[0]=0.; paro[1]=0.; paro[2]=part[1]; paro[3]=part[2]; paro[4]=part[3];
......@@ -90,8 +90,8 @@ double cfchi2(double *vrtt, double * part, VKTrack * trk)
d3 = theta_f - trk->Perig[2];
d4 = phif - trk->Perig[3];
d5 = invR_f - trk->Perig[4];
if(d4 > 3.141592653589)d4-=2.*3.141592653589;
if(d4 < -3.141592653589)d4+=2.*3.141592653589;
if(d4 > M_PI)d4-=2.*M_PI;
if(d4 < -M_PI)d4+=2.*M_PI;
res = trk->WgtM[0] * (d1*d1)
+ trk->WgtM[2] * (d2*d2)
+ trk->WgtM[5] * (d3*d3)
......@@ -311,22 +311,22 @@ void cfTrkCovarCorr(double *cov){
// ", "<<cov[13]*cov[13]/cov[9]/cov[14]<<'\n';
double Lim=0.99; double Lim2=Lim*Lim;
bool i0,i1,i2,i3,i4; i0=i1=i2=i3=i4=false;
if( fabs(cov[1]*cov[1])/cov[0]/cov[2] > Lim2 ) i0=true; i1=true;
if( fabs(cov[3]*cov[3])/cov[0]/cov[5] > Lim2 ) i0=true; i2=true;
if( fabs(cov[4]*cov[4])/cov[2]/cov[5] > Lim2 ) i1=true; i2=true;
if( fabs(cov[6]*cov[6])/cov[0]/cov[9] > Lim2 ) i0=true; i3=true;
if( fabs(cov[7]*cov[7])/cov[2]/cov[9] > Lim2 ) i1=true; i3=true;
if( fabs(cov[8]*cov[8])/cov[5]/cov[9] > Lim2 ) i2=true; i3=true;
if( fabs(cov[10]*cov[10])/cov[0]/cov[14] > Lim2 ) i0=true; i4=true;
if( fabs(cov[11]*cov[11])/cov[2]/cov[14] > Lim2 ) i1=true; i4=true;
if( fabs(cov[12]*cov[12])/cov[5]/cov[14] > Lim2 ) i2=true; i4=true;
if( fabs(cov[13]*cov[13])/cov[9]/cov[14] > Lim2 ) i3=true; i4=true;
if(i0) cov[1]*=Lim; cov[3]*=Lim; cov[6]*=Lim; cov[10]*=Lim;
if(i1) cov[1]*=Lim; cov[4]*=Lim; cov[7]*=Lim; cov[11]*=Lim;
if(i2) cov[3]*=Lim; cov[4]*=Lim; cov[8]*=Lim; cov[12]*=Lim;
if(i3) cov[6]*=Lim; cov[7]*=Lim; cov[8]*=Lim; cov[13]*=Lim;
if(i4) cov[10]*=Lim; cov[11]*=Lim; cov[12]*=Lim; cov[13]*=Lim;
if( fabs(cov[1]*cov[1])/cov[0]/cov[2] > Lim2 ){ i0=true; i1=true;}
if( fabs(cov[3]*cov[3])/cov[0]/cov[5] > Lim2 ){ i0=true; i2=true;}
if( fabs(cov[4]*cov[4])/cov[2]/cov[5] > Lim2 ){ i1=true; i2=true;}
if( fabs(cov[6]*cov[6])/cov[0]/cov[9] > Lim2 ){ i0=true; i3=true;}
if( fabs(cov[7]*cov[7])/cov[2]/cov[9] > Lim2 ){ i1=true; i3=true;}
if( fabs(cov[8]*cov[8])/cov[5]/cov[9] > Lim2 ){ i2=true; i3=true;}
if( fabs(cov[10]*cov[10])/cov[0]/cov[14] > Lim2 ){ i0=true; i4=true;}
if( fabs(cov[11]*cov[11])/cov[2]/cov[14] > Lim2 ){ i1=true; i4=true;}
if( fabs(cov[12]*cov[12])/cov[5]/cov[14] > Lim2 ){ i2=true; i4=true;}
if( fabs(cov[13]*cov[13])/cov[9]/cov[14] > Lim2 ){ i3=true; i4=true;}
if(i0){ cov[1]*=Lim; cov[3]*=Lim; cov[6]*=Lim; cov[10]*=Lim; }
if(i1){ cov[1]*=Lim; cov[4]*=Lim; cov[7]*=Lim; cov[11]*=Lim; }
if(i2){ cov[3]*=Lim; cov[4]*=Lim; cov[8]*=Lim; cov[12]*=Lim; }
if(i3){ cov[6]*=Lim; cov[7]*=Lim; cov[8]*=Lim; cov[13]*=Lim; }
if(i4){ cov[10]*=Lim; cov[11]*=Lim; cov[12]*=Lim; cov[13]*=Lim; }
}
} /*End of namespace */
......
......@@ -145,7 +145,8 @@ void vkvFastV(double *p1, double *p2, double* vRef, double dbmag, double *out)
double dir1=ptmp[0]*(cross_ref(1,1)-vRef[0])+ptmp[1]*(cross_ref(2,1)-vRef[1]);
double dir2=ptmp[0]*(cross_ref(1,2)-vRef[0])+ptmp[1]*(cross_ref(2,2)-vRef[1]);
if(dir1*dir2<0){ // points are on different sides of ref. vertex
if(dir1<0) dz1+=1000.; if(dir2<0) dz2+=1000.;
if(dir1<0) dz1+=1000.;
if(dir2<0) dz2+=1000.;
}
}
/* - choice of best dz */
......
......@@ -236,8 +236,8 @@ extern DerivT derivt_;
dzp[kt] -= zp; dphi[kt] -= phip;
zp += xyz[2]; //To gain precision
phip += phi_ini; //To gain precision
while(dphi[kt] > 3.141592653589)dphi[kt]-=2.*3.141592653589;
while(dphi[kt] < -3.141592653589)dphi[kt]+=2.*3.141592653589;
while(dphi[kt] > M_PI)dphi[kt]-=2.*M_PI;
while(dphi[kt] < -M_PI)dphi[kt]+=2.*M_PI;
//std::cout.precision(11);
//std::cout<<" newpar="<<deps[kt]<<", "<<dzp[kt]<<", "<<dtet[kt]<<", "<<dphi[kt]<<", "<<drho[kt]<<", "<<trk->Id<<'\n';
......@@ -525,7 +525,7 @@ extern DerivT derivt_;
for (kt = 0; kt < NTRK; ++kt) { /* variation on PAR is WCI * TT - (WBCI)t * DXYZ */
for( int tpn=0; tpn<9; tpn++){
if(tpn<6)twci[tpn] = vk->tmpArr[kt]->wci[tpn];
twbci[tpn] = vk->tmpArr[kt]->wbci[tpn];
twbci[tpn] = vk->tmpArr[kt]->wbci[tpn];
}
double tt1=vk->tmpArr[kt]->tt[0];
double tt2=vk->tmpArr[kt]->tt[1];
......@@ -792,7 +792,8 @@ extern DerivT derivt_;
//Having 3 points (0,-0.02,0.02) find a pabolic minimum
if (j == 6) { bet = finter(chi2t[4], chi2t[3], chi2t[5], -0.02, 0., 0.02);
if (chi2t[3] == chi2t[4] && chi2t[4] == chi2t[5]) bet = 0.;
if (bet >0.3)bet = 0.3;if (bet <-0.3)bet = -0.3;
if (bet >0.3)bet = 0.3;
if (bet <-0.3)bet = -0.3;
}
//Check if minimum is really a minimum. Otherwise use bet=0. point
if (j == 7 && ContribC[6]>ContribC[3]) { bet = 0.; j--; /*repeat step 7*/}
......
......@@ -203,8 +203,10 @@ int getFullVrtCov(VKVertex * vk, double *ader, double *dcv, double *verr)
delete[] tv; delete[] tr; delete[] tw;
IERR=0; //return IERR;
}
for (i=1; i<=NVar; ++i) for (j = 1; j<=NVar; ++j) ader_ref(i,j)*=Scale[i-1]*Scale[j-1]; delete[] Scale; //Restore scale
for(i=0; i<NVar+1; i++) delete[] ta[i]; delete[] ta; //Clean memory
for (i=1; i<=NVar; ++i) for (j = 1; j<=NVar; ++j) ader_ref(i,j)*=Scale[i-1]*Scale[j-1];
delete[] Scale; //Restore scale
for(i=0; i<NVar+1; i++) delete[] ta[i];
delete[] ta; //Clean memory
/* ---------------------------------------- */
} else {
/* ---------------------------------------- */
......@@ -356,8 +358,10 @@ int getFullVrtCov(VKVertex * vk, double *ader, double *dcv, double *verr)
}
}
// Delete temporary matrices
for(ic=0; ic<totNC; ic++) delete[] R[ic]; delete[] R;
for(ic=0; ic<totNC; ic++) delete[] RC[ic]; delete[] RC;
for(ic=0; ic<totNC; ic++) delete[] R[ic];
delete[] R;
for(ic=0; ic<totNC; ic++) delete[] RC[ic];
delete[] RC;
delete[] RCRt;
//for(int ii=1; ii<=9; ii++)std::cout<<ader_ref(ii,ii)<<", "; std::cout<<" avery full m NEW"<<'\n';
} //end of Avery matrix
......
......@@ -96,7 +96,8 @@ void vpderiv(bool UseTrackErr, long int Charge, double *pari0, double *covi, dou
par[0] = 0.;
par[1] = 0.;
par[2] = acos(pari[5] / sqrt(pp));
if(par[2]<1.e-5) par[2]=1.e-5; if(par[2]>3.141592653-1.e-5) par[2]=3.141592653-1.e-5;
if(par[2]<1.e-5) par[2]=1.e-5;
if(par[2]>M_PI-1.e-5) par[2]=M_PI-1.e-5;
par[3] = atan2(pari[4], pari[3]);
par[4] = rho;
if (Charge == 0) par[4] = constB / pt;
......
......@@ -56,7 +56,8 @@ void xyztrp(long int *ich, double *vrt0, double *pv0, double *covi, double *par
par[0] = 0.; /* (-Yv*cos + Xv*sin) */
par[1] = 0.; /* Zv - cotth*(Xv*cos + Yv*sin) */
par[2] = acos(pv0[2] / sqrt(pp));
if(par[2]<1.e-5)par[2]=1.e-5; if(par[2]>3.141592653-1.e-5) par[2]=3.141592653-1.e-5;
if(par[2]<1.e-5)par[2]=1.e-5;
if(par[2]>M_PI-1.e-5) par[2]=M_PI-1.e-5;
par[3] = atan2(pv0[1], pv0[0]);
par[4] = rho;
if ((*ich) == 0)par[4] = const__ / pt;
......@@ -152,7 +153,8 @@ void combinedTrack(long int ICH, double *vrt0, double *pv0, double *covi, doubl
par[0] = 0.; /* (-Yv*cos + Xv*sin) */
par[1] = 0.; /* Zv - cotth*(Xv*cos + Yv*sin) */
par[2] = acos(pv0[2] / sqrt(pp));
if(par[2]<1.e-5)par[2]=1.e-5; if(par[2]>3.141592653-1.e-5) par[2]=3.141592653-1.e-5;
if(par[2]<1.e-5)par[2]=1.e-5;
if(par[2]>M_PI-1.e-5) par[2]=M_PI-1.e-5;
par[3] = atan2(pv0[1], pv0[0]);
par[4] = rho;
if ( ICH==0 )par[4] = const__ / pt;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment