Check compilation warnings

Tracking/TrkVertexFitter/TrkVKalVrtCore/src/CFit.cxx

@@ -171,7 +171,7 @@ long int fitVertex(VKVertex * vk, long int iflag)
 {
     long int  i, jerr, tk, it=0;
 
-    double chi2df, dparst[6];
+    double dparst[6];
     double chi2min, dxyzst[3],  chi21s=11., chi22s=10., vShift;
     double aermd[30],tmpd[30];  // temporary array
     double tmpPer[5],tmpCov[15], tmpWgt[15];
@@ -206,7 +206,6 @@ long int fitVertex(VKVertex * vk, long int iflag)
     VKTrack * trk=0;
 
     chi2min = 1e15;
-    chi2df = 0.;
 
     if( forcft_1.nmcnst && forcft_1.useMassCnst )  {       //mass constraints are present
       std::vector<int> index;
@@ -389,6 +388,7 @@ long int fitVertex(VKVertex * vk, long int iflag)
 
 	  if ( vk->passNearVertex && it>1 ) {  //No necessary information at first iteration
             jerr = afterFit(vk, workarray_1.ader, vk->FVC.dcv, PartMom, VrtMomCov);
+            if(jerr!=0) return -17;  // Non-invertable matrix for combined track
 	    cfdcopy( PartMom, &dparst[3], 3);  //vertex part of it is filled above
             cfdcopy(VrtMomCov,vk->FVC.dcovf,21);  //Used in chi2 caclulation later...
 	    cfdcopy(  PartMom, vk->fitMom, 3);          //save Momentum
@@ -425,7 +425,7 @@ long int fitVertex(VKVertex * vk, long int iflag)
 	//std::cout<<"NNVertex="<<dxyzst[0]<<", "<<dxyzst[1]<<", "<<dxyzst[2]<<'\n';
 	//std::cout<<"-----------------------------------------------"<<'\n';
 /*  Test of convergence */
-	chi2df = fabs(chi21s - chi22s);
+	double chi2df = fabs(chi21s - chi22s);
   /*---------------------Normal convergence--------------------*/
         double PrecLimit = min(chi22s*1.e-4, forcft_1.IterationPrecision);
 	//std::cout<<"Convergence="<< chi2df <<"<"<<PrecLimit<<" cnst="<<cnstRemnants<<"<"<<ConstraintAccuracy<<'\n';
@@ -447,6 +447,7 @@ long int fitVertex(VKVertex * vk, long int iflag)
 	if ( vk->passNearVertex ) {
             if(it==1)jerr = afterFit(vk,                0, vk->FVC.dcv, PartMom, VrtMomCov);
             else     jerr = afterFit(vk, workarray_1.ader, vk->FVC.dcv, PartMom, VrtMomCov);
+            if(jerr!=0) return -17;  // Non-invertable matrix for combined track
             for( i=0; i<3; i++) dparst[i] = vk->refIterV[i]+vk->fitV[i]; // fitted vertex at global frame
             cfdcopy( PartMom, &dparst[3], 3);
             cfdcopy(VrtMomCov,vk->FVC.dcovf,21);  //Used in chi2 caclulation later...

Tracking/TrkVertexFitter/TrkVKalVrtCore/src/CascadeUtils.cxx

@@ -87,9 +87,9 @@ int fixPseudoTrackPt(long int NPar, double * fullMtx, double * LSide)
           DerivT[iniPosTrk+it*3+2] =  (sinth2sum*pt*cth)/(curv*ptsum);                 //  dTheta/dC_i      
 	  tpx+=pp[0]; tpy+=pp[1];
         }
-        double iniV0Curv=myMagFld.getCnvCst()*vMagFld[iv]/sqrt(tpx*tpx+tpy*tpy);    //initial PseudoTrack Curvature
-        if(csum<0)iniV0Curv *= -1.;
-        iniV0Curv *= vMagFld[ivnext]/vMagFld[iv];  //magnetic field correction
+//        double iniV0Curv=myMagFld.getCnvCst()*vMagFld[iv]/sqrt(tpx*tpx+tpy*tpy);    //initial PseudoTrack Curvature
+//        if(csum<0)iniV0Curv *= -1.;
+//        iniV0Curv *= vMagFld[ivnext]/vMagFld[iv];  //magnetic field correction
 //
 //fill Full Matrix and left side vector
 //

Tracking/TrkVertexFitter/TrkVKalVrtCore/src/Derclc1.cxx

@@ -2,103 +2,103 @@
   Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
 */
 
-#include <math.h> 
-#include "TrkVKalVrtCore/Derivt.h"
-#include "TrkVKalVrtCore/CommonPars.h"
-#include "TrkVKalVrtCore/TrkVKalVrtCore.h"
-#include <iostream>
-
-namespace Trk {
-
-
-//
-//   Mass constraint calculation in new data model.
-//
-//      cnstV and cnstP values are used!!!
-//-----------------------------------------------
-
-extern std::vector<double> getCnstParticleMom( VKTrack * );
-
-void  calcMassConstraint( VKMassConstraint * cnst )
-{
-    int it,itc;
-    double ptot[4]={0.,0.,0.,0.};  
-    double cth, invR, pp2, pt; 
-    VKConstraintBase * base_cnst = (VKConstraintBase*) cnst;
-    VKVertex * vk=cnst->getOriginVertex();
-    int usedNTRK = cnst->usedParticles.size();
-    VKTrack * trk;
-    std::vector< std::vector<double> > pp(usedNTRK);
-    for( itc=0; itc<usedNTRK; itc++){
-      it = cnst->usedParticles[itc];
-      trk = vk->TrackList.at(it);
-      pp[itc]=getCnstParticleMom( trk );
-      ptot[0] += pp[itc][0];    
-      ptot[1] += pp[itc][1];    
-      ptot[2] += pp[itc][2];    
-      ptot[3] += pp[itc][3];    
-   }
-   
-   double temp = 1.e-10;
-   int ip=0; if( fabs(ptot[1]) > fabs(ptot[ip]) )ip=1; if( fabs(ptot[2]) > fabs(ptot[ip]) )ip=2;
-   int im=0; if( fabs(ptot[1]) < fabs(ptot[im]) )im=1; if( fabs(ptot[2]) < fabs(ptot[im]) )im=2;
-   int id=4; for(int i=0;i<3;i++) if(i!=ip && i!=im)id=i;
-   if(id==4){
-     std::cout<<"ERROR in mass constraint!!!"<<'\n';
-     temp=ptot[3]*ptot[3]-ptot[0]*ptot[0]-ptot[1]*ptot[1]-ptot[2]*ptot[2];
-   } else {
-     temp = sqrt( (ptot[3]-ptot[ip])*(ptot[3]+ptot[ip]) );
-     temp = sqrt( (  temp -ptot[id])*(  temp +ptot[id]) );
-     temp =       (  temp -ptot[im])*(  temp +ptot[im]);
-   }
-   if(temp<1.e-10)temp=1.e-10;  //protection
-   temp=sqrt(temp);             //system mass
-//
-//
-    int numCNST=0;   //constraint number. Single constraint in this case
-//
-//Difference   
-   base_cnst->aa[numCNST] = ( temp - cnst->targetMass ) * ( temp + cnst->targetMass ); 
-//
-//Derivatives               Here pp[][3] - particle energy, pp[][4] - squared particle mom
-    for( itc=0; itc<usedNTRK; itc++){
-      it = cnst->usedParticles[itc];
-      trk  = vk->TrackList.at(it);
-      invR = trk->cnstP[2];
-      cth  = 1. / tan( trk->cnstP[0] );
-      pt   = sqrt(pp[itc][0]*pp[itc][0] + pp[itc][1]*pp[itc][1]);
-      pp2  =      pp[itc][0]*pp[itc][0] + pp[itc][1]*pp[itc][1] + pp[itc][2]*pp[itc][2];
-      base_cnst->f0t[it][numCNST].X =  ptot[3] * (-pp2* cth / pp[itc][3]) 
-			             - ptot[2] * (-pt   * (cth*cth + 1.));
-      base_cnst->f0t[it][numCNST].Y = -ptot[0] * (-pp[itc][1]) 
-			             - ptot[1] *   pp[itc][0];
-      base_cnst->f0t[it][numCNST].Z =  ptot[3] * (-pp2/ (invR * pp[itc][3])) 
-	                             - ptot[0] * (-pp[itc][0] / invR) 
-			             - ptot[1] * (-pp[itc][1] / invR) 
-			             - ptot[2] * (-pp[itc][2] / invR);
-
-    }
-    base_cnst->h0t[numCNST].X = 0.;
-    base_cnst->h0t[numCNST].Y = 0.;
-    base_cnst->h0t[numCNST].Z = 0.;
-/* -- Relative normalisation. Now it is target mass (squared?) to make performance uniform */
-    double Scale=0.025/(2.*cnst->targetMass+1.); //28.03.2011 wrong idea for cascade. VK 06.05.2011 actually correct!
-    //Scale=0.01;
-    for (it = 0; it < (int)vk->TrackList.size(); ++it) {
-	base_cnst->f0t[it][numCNST].X *=  Scale * 2;
-	base_cnst->f0t[it][numCNST].Y *=  Scale * 2;
-	base_cnst->f0t[it][numCNST].Z *=  Scale * 2;
-    }
-    base_cnst->aa[numCNST] *= Scale;
-//std::cout.precision(11);
-//std::cout<<" mass="<<temp<<", "<<cnst->targetMass<<", "<<base_cnst->aa[numCNST]<<'\n';
-//std::cout<<base_cnst->f0t[0][numCNST].X<<", "<<base_cnst->f0t[0][numCNST].Y<<", "
-//         <<base_cnst->f0t[0][numCNST].Z<<'\n';
-//std::cout<<base_cnst->f0t[1][numCNST].X<<", "<<base_cnst->f0t[1][numCNST].Y<<", "
-//         <<base_cnst->f0t[1][numCNST].Z<<'\n';
-}
-
-
-
-
-} /* End of namespace */
+#include <math.h> 
+#include "TrkVKalVrtCore/Derivt.h"
+#include "TrkVKalVrtCore/CommonPars.h"
+#include "TrkVKalVrtCore/TrkVKalVrtCore.h"
+#include <iostream>
+
+namespace Trk {
+
+
+//
+//   Mass constraint calculation in new data model.
+//
+//      cnstV and cnstP values are used!!!
+//-----------------------------------------------
+
+extern std::vector<double> getCnstParticleMom( VKTrack * );
+
+void  calcMassConstraint( VKMassConstraint * cnst )
+{
+    int it,itc;
+    double ptot[4]={0.,0.,0.,0.};  
+    double cth, invR, pp2, pt; 
+    VKConstraintBase * base_cnst = (VKConstraintBase*) cnst;
+    VKVertex * vk=cnst->getOriginVertex();
+    int usedNTRK = cnst->usedParticles.size();
+    VKTrack * trk;
+    std::vector< std::vector<double> > pp(usedNTRK);
+    for( itc=0; itc<usedNTRK; itc++){
+      it = cnst->usedParticles[itc];
+      trk = vk->TrackList.at(it);
+      pp[itc]=getCnstParticleMom( trk );
+      ptot[0] += pp[itc][0];    
+      ptot[1] += pp[itc][1];    
+      ptot[2] += pp[itc][2];    
+      ptot[3] += pp[itc][3];    
+   }
+   
+   double temp = 1.e-10;
+   int ip=0; if( fabs(ptot[1]) > fabs(ptot[ip]) )ip=1; if( fabs(ptot[2]) > fabs(ptot[ip]) )ip=2;
+   int im=0; if( fabs(ptot[1]) < fabs(ptot[im]) )im=1; if( fabs(ptot[2]) < fabs(ptot[im]) )im=2;
+   int id=4; for(int i=0;i<3;i++) if(i!=ip && i!=im)id=i;
+   if(id==4){
+     std::cout<<"ERROR in mass constraint!!!"<<'\n';
+     temp=ptot[3]*ptot[3]-ptot[0]*ptot[0]-ptot[1]*ptot[1]-ptot[2]*ptot[2];
+   } else {
+     temp = sqrt( (ptot[3]-ptot[ip])*(ptot[3]+ptot[ip]) );
+     temp = sqrt( (  temp -ptot[id])*(  temp +ptot[id]) );
+     temp =       (  temp -ptot[im])*(  temp +ptot[im]);
+   }
+   if(temp<1.e-10)temp=1.e-10;  //protection
+   temp=sqrt(temp);             //system mass
+//
+//
+    int numCNST=0;   //constraint number. Single constraint in this case
+//
+//Difference   
+   base_cnst->aa[numCNST] = ( temp - cnst->targetMass ) * ( temp + cnst->targetMass ); 
+//
+//Derivatives               Here pp[][3] - particle energy, pp[][4] - squared particle mom
+    for( itc=0; itc<usedNTRK; itc++){
+      it = cnst->usedParticles[itc];
+      trk  = vk->TrackList.at(it);
+      invR = trk->cnstP[2];
+      cth  = 1. / tan( trk->cnstP[0] );
+      pt   = sqrt(pp[itc][0]*pp[itc][0] + pp[itc][1]*pp[itc][1]);
+      pp2  =      pp[itc][0]*pp[itc][0] + pp[itc][1]*pp[itc][1] + pp[itc][2]*pp[itc][2];
+      base_cnst->f0t[it][numCNST].X =  ptot[3] * (-pp2* cth / pp[itc][3]) 
+			             - ptot[2] * (-pt   * (cth*cth + 1.));
+      base_cnst->f0t[it][numCNST].Y = -ptot[0] * (-pp[itc][1]) 
+			             - ptot[1] *   pp[itc][0];
+      base_cnst->f0t[it][numCNST].Z =  ptot[3] * (-pp2/ (invR * pp[itc][3])) 
+	                             - ptot[0] * (-pp[itc][0] / invR) 
+			             - ptot[1] * (-pp[itc][1] / invR) 
+			             - ptot[2] * (-pp[itc][2] / invR);
+
+    }
+    base_cnst->h0t[numCNST].X = 0.;
+    base_cnst->h0t[numCNST].Y = 0.;
+    base_cnst->h0t[numCNST].Z = 0.;
+/* -- Relative normalisation. Now it is target mass (squared?) to make performance uniform */
+    double Scale=0.025/(2.*cnst->targetMass+1.); //28.03.2011 wrong idea for cascade. VK 06.05.2011 actually correct!
+    //Scale=0.01;
+    for (it = 0; it < (int)vk->TrackList.size(); ++it) {
+	base_cnst->f0t[it][numCNST].X *=  Scale * 2;
+	base_cnst->f0t[it][numCNST].Y *=  Scale * 2;
+	base_cnst->f0t[it][numCNST].Z *=  Scale * 2;
+    }
+    base_cnst->aa[numCNST] *= Scale;
+//std::cout.precision(11);
+//std::cout<<" mass="<<temp<<", "<<cnst->targetMass<<", "<<base_cnst->aa[numCNST]<<'\n';
+//std::cout<<base_cnst->f0t[0][numCNST].X<<", "<<base_cnst->f0t[0][numCNST].Y<<", "
+//         <<base_cnst->f0t[0][numCNST].Z<<'\n';
+//std::cout<<base_cnst->f0t[1][numCNST].X<<", "<<base_cnst->f0t[1][numCNST].Y<<", "
+//         <<base_cnst->f0t[1][numCNST].Z<<'\n';
+}
+
+
+
+
+} /* End of namespace */

Tracking/TrkVertexFitter/TrkVKalVrtCore/src/VtCFitC.cxx

@@ -62,6 +62,7 @@ int vtcfitc( VKVertex * vk )
 	 tf0t.push_back( tmpVec );
        }
     }
+    if(totNC==0)return 0;
     tmpVec.clear();
 //
     std::vector< std::vector<double> > denom;

Tracking/TrkVertexFitter/TrkVKalVrtCore/src/VtDeriv.cxx

@@ -27,7 +27,7 @@ void vpderiv(bool UseTrackErr, long int Charge, double *pari0, double *covi, dou
     /* Local variables */
     double pari[6], covd[15], dcov[3], rvec[50]; /* was [2][6*4+1] */
     double  paro[5];
-    long int jerr, j, ij, ip, ipp, id=0;
+    long int j, ij, ip, ipp, id=0;
     double dwgt0[3]={0.,0.,0.}, constB;
     //double deriv[6], dchi2[4*6+1];                     //VK for debugging
     double cs, pp, sn, pt, rho;
@@ -66,7 +66,6 @@ void vpderiv(bool UseTrackErr, long int Charge, double *pari0, double *covi, dou
 
     /* Function Body */
 /* --------------------- */
-    jerr = 0;
 //VK    constB = *localbmag * .0029979246;
     constB =vkalvrtbmag.bmag  * vkalMagCnvCst;
     
@@ -181,7 +180,7 @@ void vpderiv(bool UseTrackErr, long int Charge, double *pari0, double *covi, dou
 /*---------------------------------------------------------------- */	
 //     Weight matrix and chi2 for given shift	
 //
-	jerr=cfdinv(dcov, &dwgt[0], -2); if(jerr){jerr=cfdinv(dcov, &dwgt[0], 2); if(jerr){dwgt[0]=dwgt[2]=1.e6; dwgt[1]=0.; jerr=0;}};
+	int jerr=cfdinv(dcov, &dwgt[0], -2); if(jerr){jerr=cfdinv(dcov, &dwgt[0], 2); if(jerr){dwgt[0]=dwgt[2]=1.e6; dwgt[1]=0.;}};
 	//dchi2[ip] = sqrt(fabs(dwgt[0]*paro[0]*paro[0] + 2.*dwgt[1]*paro[0]*paro[1] + dwgt[2]*paro[1]*paro[1]));
 	rvec_ref(1, ip) = paro[0];
 	rvec_ref(2, ip) = paro[1];

Tracking/TrkVertexFitter/TrkVKalVrtCore/src/cfImp.cxx

@@ -2,185 +2,185 @@
   Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
 */
 
-#include <math.h>
-#include "TrkVKalVrtCore/Propagator.h"
-#include <iostream>
-
-namespace Trk {
-
-extern vkalPropagator  myPropagator;
-extern int cfdinv(double *, double *, long int); 
-
-#define min(a,b) ((a) <= (b) ? (a) : (b))
-#define max(a,b) ((a) >= (b) ? (a) : (b))
-
-
- void   cfimp(long int TrkID, long int ich, long int IFL, double *par, 
-	double *err, double *vrt, double *vcov, double *rimp, 
-	double *rcov, double *sign )
-{
-    double dcov[3], errd[15], paro[5];
-    double dwgt[3], errn[15];
-    long int jerr, i__, j, ij;
-
-
-    double cs, sn;
-
-    double cnv[6]	/* was [2][3] */;
-/* --------------------------------------------------------- */
-/*  Impact parameter calculation                             */
-/*  Input:                                                   */
-/*      ICH -charge(-1,0,1)                                  */
-/*      IFL = 1 contribution from VRT is added               */
-/*      IFL =-1 contribution from VRT is subtructed          */
-/*      PAR(5),ERR(15) - peregee parameters and error matrix */
-/*      VRT(3),VCOV(6) - vertex and its error matrix         */
-/*                                                           */
-/*    SIGNIFICANCE IS CALCULATED FOR PERIGEE POINT BY DEF.   */
-/*           (NOT FOR THE CLOSEST POINT!!!)                  */
-/*                                                           */
-/*  Output:                                                  */
-/*     RIMP(1) - impact in XY                                */
-/*     RIMP(2) - impact in Z                                 */
-/*     RIMP(3) - Theta at new vertex                         */
-/*     RIMP(4) - Phi at new vertex                           */
-/*     RIMP(5) - curvature at vertex                         */
-/*     RCOV(3) - error matrix for RIMP                       */
-/*     SIGN    - impact significance                         */
-/* Author: V.Kostyukhin                                      */
-/* --------------------------------------------------------- */
-    /* Parameter adjustments */
-    --vcov;
-
-    /* Function Body */
-    for (int ii = 0; ii < 15; ++ii) {errd[ii] = err[ii];}
-
-
-    double Ref0[3]={0.,0.,0.};  //base reference point for standard perigee
-    myPropagator.Propagate( TrkID, ich, par, errd, Ref0, vrt, paro, errn);
-
-//std::cout <<" CFImp new par R,Z="<<paro[0]<<", "<<paro[1]<<'\n';
-/* ---------- */
-    rimp[0] = paro[0];
-    rimp[1] = paro[1];
-    rimp[2] = paro[2];
-    rimp[3] = paro[3];
-    rimp[4] = paro[4];
-/*  X=paro(1)*sn, Y=-paro(1)*cs */
-    sn = sin(paro[3]);
-    cs = cos(paro[3]);
-/* -- New error version */
-    cnv[0] = -sn;
-    cnv[2] = cs;
-    cnv[4] = 0.;
-    cnv[1] = sn / tan(paro[2]);
-    cnv[3] = cs / tan(paro[2]);
-    cnv[5] = -1.;
-    dcov[0] = 0.;
-    dcov[1] = 0.;
-    dcov[2] = 0.;
-    for (i__ = 1; i__ <= 3; ++i__) {
-	for (j = 1; j <= 3; ++j) {
-	    ij = max(i__,j);
-	    ij = ij * (ij - 1) / 2 + min(i__,j);
-	    dcov[0] += cnv[(i__ << 1) - 2] * cnv[(j << 1) - 2] * vcov[ij];
-	    dcov[2] += cnv[(i__ << 1) - 1] * cnv[(j << 1) - 1] * vcov[ij];
-	    dcov[1] += cnv[(i__ << 1) - 2] * cnv[(j << 1) - 1] * vcov[ij];
-	}
-    }
-/* --------------------------------------------------------------- */
-    if (IFL == 1) {
-	rcov[0] = errn[0] + dcov[0];
-	rcov[1] = errn[1] + dcov[1];
-	rcov[2] = errn[2] + dcov[2];
-    } else if (IFL == -1) {
-	rcov[0] = errn[0] - dcov[0];
-	rcov[1] = errn[1] - dcov[1];
-	rcov[2] = errn[2] - dcov[2];
-    } else {
-	rcov[0] = errn[0];
-	rcov[1] = errn[1];
-	rcov[2] = errn[2];
-    }
-    jerr=cfdinv(rcov, dwgt, -2);
-    if (jerr) {jerr=cfdinv(rcov, dwgt, 2); if(jerr){dwgt[0]=dwgt[2]=1.e6; dwgt[1]=0.; jerr=0;}}
-    (*sign) = sqrt(fabs(dwgt[0] * rimp[0] * rimp[0] + dwgt[1] * 2. * rimp[0] * rimp[1] + 
-	    dwgt[2] * rimp[1] * rimp[1]));
-} 
-
-
- void   cfimpc(long int TrkID, long int ich, long int IFL, double *par, 
-	double *err, double *vrt, double *vcov, double *rimp, 
-	double *rcov, double *sign )
-{
-    double dcov[3], errd[15], paro[5];
-    double dwgt[3], errn[15], cnv[6];	/* was [2][3] */
-    long int jerr, i__, j, ij;
-
-
-    double cs, sn;
-    extern void cfClstPnt(double *p, double *, double *);
-/* --------------------------------------------------------- */
-/*    SIGNIFICANCE IS CALCULATED FOR THE CLOSEST POINT NOW!!!*/
-/* Author: V.Kostyukhin                                      */
-/* --------------------------------------------------------- */
-    /* Parameter adjustments */
-    --vcov;
-    for (int ii = 0; ii < 15; ++ii) {errd[ii] = err[ii];}
-
-    double Ref0[3]={0.,0.,0.};  //base reference point for standard perigee
-    myPropagator.Propagate( TrkID, ich, par, errd, Ref0, vrt, paro, errn);
-
-    double tmpVrt[3]={0.,0.,0.}; double ClosestPnt[3];
-    cfClstPnt( paro, tmpVrt, ClosestPnt);
-    paro[0]=sqrt(ClosestPnt[0]*ClosestPnt[0] + ClosestPnt[1]*ClosestPnt[1]);
-    paro[1]=ClosestPnt[2];
-/* ---------- */
-    rimp[0] = paro[0];
-    rimp[1] = paro[1];
-    rimp[2] = paro[2];
-    rimp[3] = paro[3];
-    rimp[4] = paro[4];
-    sn = sin(paro[3]);
-    cs = cos(paro[3]);
-/* -- New error version */
-    cnv[0] = -sn;
-    cnv[2] = cs;
-    cnv[4] = 0.;
-    cnv[1] = sn / tan(paro[2]);
-    cnv[3] = cs / tan(paro[2]);
-    cnv[5] = -1.;
-    dcov[0] = 0.;
-    dcov[1] = 0.;
-    dcov[2] = 0.;
-    for (i__ = 1; i__ <= 3; ++i__) {
-	for (j = 1; j <= 3; ++j) {
-	    ij = max(i__,j);
-	    ij = ij * (ij - 1) / 2 + min(i__,j);
-	    dcov[0] += cnv[(i__ << 1) - 2] * cnv[(j << 1) - 2] * vcov[ij];
-	    dcov[2] += cnv[(i__ << 1) - 1] * cnv[(j << 1) - 1] * vcov[ij];
-	    dcov[1] += cnv[(i__ << 1) - 2] * cnv[(j << 1) - 1] * vcov[ij];
-	}
-    }
-/* --------------------------------------------------------------- */
-    if (IFL == 1) {
-	rcov[0] = errn[0] + dcov[0];
-	rcov[1] = errn[1] + dcov[1];
-	rcov[2] = errn[2] + dcov[2];
-    } else if (IFL == -1) {
-	rcov[0] = errn[0] - dcov[0];
-	rcov[1] = errn[1] - dcov[1];
-	rcov[2] = errn[2] - dcov[2];
-    } else {
-	rcov[0] = errn[0];
-	rcov[1] = errn[1];
-	rcov[2] = errn[2];
-    }
-    jerr=cfdinv(rcov, dwgt, -2);
-    if (jerr) {jerr=cfdinv(rcov, dwgt, 2);if(jerr){dwgt[0]=dwgt[2]=1.e6; dwgt[1]=0.; jerr=0;}}
-    (*sign) = sqrt(fabs(dwgt[0] * rimp[0] * rimp[0] + dwgt[1] * 2. * rimp[0] * rimp[1] + 
-	    dwgt[2] * rimp[1] * rimp[1]));
-} 
-
-} /* end of namespace */
-
+#include <math.h>
+#include "TrkVKalVrtCore/Propagator.h"
+#include <iostream>
+
+namespace Trk {
+
+extern vkalPropagator  myPropagator;
+extern int cfdinv(double *, double *, long int); 
+
+#define min(a,b) ((a) <= (b) ? (a) : (b))
+#define max(a,b) ((a) >= (b) ? (a) : (b))
+
+
+ void   cfimp(long int TrkID, long int ich, long int IFL, double *par, 
+	double *err, double *vrt, double *vcov, double *rimp, 
+	double *rcov, double *sign )
+{
+    double dcov[3], errd[15], paro[5];
+    double dwgt[3], errn[15];
+    int i__, j, ij;
+
+
+    double cs, sn;
+
+    double cnv[6]	/* was [2][3] */;
+/* --------------------------------------------------------- */
+/*  Impact parameter calculation                             */
+/*  Input:                                                   */
+/*      ICH -charge(-1,0,1)                                  */
+/*      IFL = 1 contribution from VRT is added               */
+/*      IFL =-1 contribution from VRT is subtructed          */
+/*      PAR(5),ERR(15) - peregee parameters and error matrix */
+/*      VRT(3),VCOV(6) - vertex and its error matrix         */
+/*                                                           */
+/*    SIGNIFICANCE IS CALCULATED FOR PERIGEE POINT BY DEF.   */
+/*           (NOT FOR THE CLOSEST POINT!!!)                  */
+/*                                                           */
+/*  Output:                                                  */
+/*     RIMP(1) - impact in XY                                */
+/*     RIMP(2) - impact in Z                                 */
+/*     RIMP(3) - Theta at new vertex                         */
+/*     RIMP(4) - Phi at new vertex                           */
+/*     RIMP(5) - curvature at vertex                         */
+/*     RCOV(3) - error matrix for RIMP                       */
+/*     SIGN    - impact significance                         */
+/* Author: V.Kostyukhin                                      */
+/* --------------------------------------------------------- */
+    /* Parameter adjustments */
+    --vcov;
+
+    /* Function Body */
+    for (int ii = 0; ii < 15; ++ii) {errd[ii] = err[ii];}
+
+
+    double Ref0[3]={0.,0.,0.};  //base reference point for standard perigee
+    myPropagator.Propagate( TrkID, ich, par, errd, Ref0, vrt, paro, errn);
+
+//std::cout <<" CFImp new par R,Z="<<paro[0]<<", "<<paro[1]<<'\n';
+/* ---------- */
+    rimp[0] = paro[0];
+    rimp[1] = paro[1];
+    rimp[2] = paro[2];
+    rimp[3] = paro[3];
+    rimp[4] = paro[4];
+/*  X=paro(1)*sn, Y=-paro(1)*cs */
+    sn = sin(paro[3]);
+    cs = cos(paro[3]);
+/* -- New error version */
+    cnv[0] = -sn;
+    cnv[2] = cs;
+    cnv[4] = 0.;
+    cnv[1] = sn / tan(paro[2]);
+    cnv[3] = cs / tan(paro[2]);
+    cnv[5] = -1.;
+    dcov[0] = 0.;
+    dcov[1] = 0.;
+    dcov[2] = 0.;
+    for (i__ = 1; i__ <= 3; ++i__) {
+	for (j = 1; j <= 3; ++j) {
+	    ij = max(i__,j);
+	    ij = ij * (ij - 1) / 2 + min(i__,j);
+	    dcov[0] += cnv[(i__ << 1) - 2] * cnv[(j << 1) - 2] * vcov[ij];
+	    dcov[2] += cnv[(i__ << 1) - 1] * cnv[(j << 1) - 1] * vcov[ij];
+	    dcov[1] += cnv[(i__ << 1) - 2] * cnv[(j << 1) - 1] * vcov[ij];
+	}
+    }
+/* --------------------------------------------------------------- */
+    if (IFL == 1) {
+	rcov[0] = errn[0] + dcov[0];
+	rcov[1] = errn[1] + dcov[1];
+	rcov[2] = errn[2] + dcov[2];
+    } else if (IFL == -1) {
+	rcov[0] = errn[0] - dcov[0];
+	rcov[1] = errn[1] - dcov[1];
+	rcov[2] = errn[2] - dcov[2];
+    } else {
+	rcov[0] = errn[0];
+	rcov[1] = errn[1];
+	rcov[2] = errn[2];
+    }
+    int jerr=cfdinv(rcov, dwgt, -2);
+    if (jerr) {jerr=cfdinv(rcov, dwgt, 2); if(jerr){dwgt[0]=dwgt[2]=1.e6; dwgt[1]=0.;}}
+    (*sign) = sqrt(fabs(dwgt[0] * rimp[0] * rimp[0] + dwgt[1] * 2. * rimp[0] * rimp[1] + 
+	    dwgt[2] * rimp[1] * rimp[1]));
+} 
+
+
+ void   cfimpc(long int TrkID, long int ich, long int IFL, double *par, 
+	double *err, double *vrt, double *vcov, double *rimp, 
+	double *rcov, double *sign )
+{
+    double dcov[3], errd[15], paro[5];
+    double dwgt[3], errn[15], cnv[6];	/* was [2][3] */
+    int i__, j, ij;
+
+
+    double cs, sn;
+    extern void cfClstPnt(double *p, double *, double *);
+/* --------------------------------------------------------- */
+/*    SIGNIFICANCE IS CALCULATED FOR THE CLOSEST POINT NOW!!!*/
+/* Author: V.Kostyukhin                                      */
+/* --------------------------------------------------------- */
+    /* Parameter adjustments */
+    --vcov;
+    for (int ii = 0; ii < 15; ++ii) {errd[ii] = err[ii];}
+
+    double Ref0[3]={0.,0.,0.};  //base reference point for standard perigee
+    myPropagator.Propagate( TrkID, ich, par, errd, Ref0, vrt, paro, errn);
+
+    double tmpVrt[3]={0.,0.,0.}; double ClosestPnt[3];
+    cfClstPnt( paro, tmpVrt, ClosestPnt);
+    paro[0]=sqrt(ClosestPnt[0]*ClosestPnt[0] + ClosestPnt[1]*ClosestPnt[1]);
+    paro[1]=ClosestPnt[2];
+/* ---------- */
+    rimp[0] = paro[0];
+    rimp[1] = paro[1];
+    rimp[2] = paro[2];
+    rimp[3] = paro[3];
+    rimp[4] = paro[4];
+    sn = sin(paro[3]);
+    cs = cos(paro[3]);
+/* -- New error version */
+    cnv[0] = -sn;
+    cnv[2] = cs;
+    cnv[4] = 0.;
+    cnv[1] = sn / tan(paro[2]);
+    cnv[3] = cs / tan(paro[2]);
+    cnv[5] = -1.;
+    dcov[0] = 0.;
+    dcov[1] = 0.;
+    dcov[2] = 0.;
+    for (i__ = 1; i__ <= 3; ++i__) {
+	for (j = 1; j <= 3; ++j) {
+	    ij = max(i__,j);
+	    ij = ij * (ij - 1) / 2 + min(i__,j);
+	    dcov[0] += cnv[(i__ << 1) - 2] * cnv[(j << 1) - 2] * vcov[ij];
+	    dcov[2] += cnv[(i__ << 1) - 1] * cnv[(j << 1) - 1] * vcov[ij];
+	    dcov[1] += cnv[(i__ << 1) - 2] * cnv[(j << 1) - 1] * vcov[ij];
+	}
+    }
+/* --------------------------------------------------------------- */
+    if (IFL == 1) {
+	rcov[0] = errn[0] + dcov[0];
+	rcov[1] = errn[1] + dcov[1];
+	rcov[2] = errn[2] + dcov[2];
+    } else if (IFL == -1) {
+	rcov[0] = errn[0] - dcov[0];
+	rcov[1] = errn[1] - dcov[1];
+	rcov[2] = errn[2] - dcov[2];
+    } else {
+	rcov[0] = errn[0];
+	rcov[1] = errn[1];
+	rcov[2] = errn[2];
+    }
+    int jerr=cfdinv(rcov, dwgt, -2);
+    if (jerr) {jerr=cfdinv(rcov, dwgt, 2);if(jerr){dwgt[0]=dwgt[2]=1.e6; dwgt[1]=0.;}}
+    (*sign) = sqrt(fabs(dwgt[0] * rimp[0] * rimp[0] + dwgt[1] * 2. * rimp[0] * rimp[1] + 
+	    dwgt[2] * rimp[1] * rimp[1]));
+} 
+
+} /* end of namespace */
+

Tracking/TrkVertexFitter/TrkVKalVrtCore/src/cfNewPm.cxx

@@ -2,133 +2,134 @@
   Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
 */
 
-#include <math.h>
-#include <iostream>
-#include "TrkVKalVrtCore/VKalVrtBMag.h"
-#include "TrkVKalVrtCore/Propagator.h"
-
-namespace Trk {
-
-extern vkalPropagator  myPropagator;
-extern VKalVrtBMag vkalvrtbmag;
-extern vkalMagFld  myMagFld;
-
-
-extern double d_sign(double, double);
-extern void cfnewp(long int *, double *, double *, double *, double *, double *);
-extern void vkgrkuta_(double *, double *, double *, double *);
-
-
-void cfnewpm(double *par, double *xyzStart, double *xyzEnd, double *ustep, double *parn, double *closePoint)
-{
-    double d__1, d__2,dist_left;
-    double vect[7], stmg, vout[7], dpar0[5];
-    double p, perig[3], dstep, constB, xyzst[3], charge, fx, fy, fz, pt, px, py, pz;
-    double posold, poscur, totway, dp;
-    double dX, dY, dZ;
-    long int ich;
-
-/* --------------------------------------------------------- */
-/*  The same as CFNEWP but for the case on nonuniform        */
-/*   magnetic field                                          */
-/*                                                           */
-/*  Propagation from xyzStart reference point to xyzEnd point*/
-/*     PAR  - perigee parameters wrt xyzStart                */
-/*     PARN - perigee parameters wrt xyzEnd                  */
-/*  Author: V.Kostyukhin                                     */
-/* --------------------------------------------------------- */
-    /* Parameter adjustments */
-    --par;
-
-    d__1 = tan(par[3]);
-    totway = (*ustep) * sqrt(1. / (d__1 * d__1) + 1.);
-
-    if (fabs(*ustep) < 10. && fabs(totway) < 20.)  return;   // Distance(mm) is small. Simplest propagation is used
-
-    stmg = 40.;  //Propagation step in mm for nonuniform field
-
-    vect[0] =  sin(par[4]) * par[1]   +xyzStart[0];
-    vect[1] = -cos(par[4]) * par[1]   +xyzStart[1];
-    vect[2] = par[2]                  +xyzStart[2];
-
-    myMagFld.getMagFld( vect[0],vect[1],vect[2],fx,fy,fz);
-    constB = fz * myMagFld.getCnvCst();
-
-    pt = constB / fabs(par[5]);
-    px = pt * cos(par[4]);
-    py = pt * sin(par[4]);
-    pz = pt / tan(par[3]);
-    p = sqrt(pt * pt + pz * pz);
-
-    vect[3] = px / p;
-    vect[4] = py / p;
-    vect[5] = pz / p;
-    vect[6] = p;
-    charge = -d_sign( 1., par[5]);
-    poscur = 0.;
-//std::cout <<"VkGrkuta vect="<<vect[0]<<", "<<vect[1]<<", "<<vect[2]<<", "<<vect[3]<<", "
-//                            <<vect[4]<<", "<<vect[5]<<", "<<vect[6]<<'\n';
-    while(fabs(poscur) < fabs(totway)) {
-	posold = poscur;
-	d__1 = fabs(poscur) + stmg;
-        d__2 = fabs(totway);
-	poscur = d__1 < d__2 ? d__1 : d__2;
-	poscur = d_sign(poscur, totway);
-	dstep = poscur - posold;
-	vkgrkuta_(&charge, &dstep, vect, vout);
-	vect[0] = vout[0];
-	vect[1] = vout[1];
-	vect[2] = vout[2];
-	vect[3] = vout[3];
-	vect[4] = vout[4];
-	vect[5] = vout[5];
-// safety for strongly nonuniform field
-	dX = vect[0]-xyzEnd[0];
-	dY = vect[1]-xyzEnd[1];
-	dZ = vect[2]-xyzEnd[2];
-	dist_left = sqrt( dX*dX + dY*dY + dZ*dZ);
-	if(dist_left < stmg) break;    // arrived
-    }
-
-/* --- Now we are in the new point. Calculate track parameters */
-/*   at new point with new mag.field */
-
-    myMagFld.getMagFld( xyzEnd[0],xyzEnd[1],xyzEnd[2],fx,fy,fz);
-    //myMagFld.getMagFld( vect[0], vect[1], vect[2], fx, fy, fz);
-    constB = fz * myMagFld.getCnvCst();
-
-    dpar0[0] = 0.;
-    dpar0[1] = 0.;
-    dpar0[2] = acos(vout[5]);
-    dpar0[3] = atan2(vout[4], vout[3]);
-//    if (dpar0[3] < 0.) dpar0[3] += 6.2831853071794;
-
-    px = vect[3] * vect[6];
-    py = vect[4] * vect[6];
-    dpar0[4] = d_sign(  (constB/sqrt(px*px+py*py)), par[5]);    /* new value of curvature */
-    ich = (long int) charge;
-    xyzst[0] = xyzEnd[0] - vout[0];
-    xyzst[1] = xyzEnd[1] - vout[1];
-    xyzst[2] = xyzEnd[2] - vout[2];
-
-    cfnewp(&ich, dpar0, xyzst, &dp, parn, perig); // Last step of propagation 
-    closePoint[0] = perig[0] + vout[0];                  // with simple program
-    closePoint[1] = perig[1] + vout[1];
-    closePoint[2] = perig[2] + vout[2];
-
-//std::cout <<" Exit from cfNewPm="<<dp<<'\n';
-//std::cout <<parn[0]<<'\n';
-//   if(parn[0] == 0.) {
-//std::cout <<dp<<", "<<(*ustep)<<'\n';
-//       for (int jj=0; jj<6;jj++) std::cout <<" vout="<<vout[jj]<<'\n';
-//       for (int jj=0; jj<3;jj++) std::cout <<" xyzst="<<xyzst[jj]<<'\n';
-//       for (int jj=0; jj<5;jj++) std::cout <<" dpar0="<<dpar0[jj]<<'\n';
-//       for (int jj=0; jj<5;jj++) std::cout <<" parn="<<parn[jj]<<'\n';
-//   }
-    return;
-} 
-
-
-
-
-} /* End of namespace */
+#include <math.h>
+#include <iostream>
+#include "TrkVKalVrtCore/VKalVrtBMag.h"
+#include "TrkVKalVrtCore/Propagator.h"
+
+namespace Trk {
+
+extern vkalPropagator  myPropagator;
+extern VKalVrtBMag vkalvrtbmag;
+extern vkalMagFld  myMagFld;
+
+
+extern double d_sign(double, double);
+extern void cfnewp(long int *, double *, double *, double *, double *, double *);
+extern void vkgrkuta_(double *, double *, double *, double *);
+
+
+void cfnewpm(double *par, double *xyzStart, double *xyzEnd, double *ustep, double *parn, double *closePoint)
+{
+    double d__1, d__2,dist_left;
+    double vect[7], stmg, dpar0[5];
+    double vout[7]={0.};
+    double p, perig[3], dstep, constB, xyzst[3], charge, fx, fy, fz, pt, px, py, pz;
+    double posold, poscur, totway, dp;
+    double dX, dY, dZ;
+    long int ich;
+
+/* --------------------------------------------------------- */
+/*  The same as CFNEWP but for the case on nonuniform        */
+/*   magnetic field                                          */
+/*                                                           */
+/*  Propagation from xyzStart reference point to xyzEnd point*/
+/*     PAR  - perigee parameters wrt xyzStart                */
+/*     PARN - perigee parameters wrt xyzEnd                  */
+/*  Author: V.Kostyukhin                                     */
+/* --------------------------------------------------------- */
+    /* Parameter adjustments */
+    --par;
+
+    d__1 = tan(par[3]);
+    totway = (*ustep) * sqrt(1. / (d__1 * d__1) + 1.);
+
+    if (fabs(*ustep) < 10. && fabs(totway) < 20.)  return;   // Distance(mm) is small. Simplest propagation is used
+
+    stmg = 40.;  //Propagation step in mm for nonuniform field
+
+    vect[0] =  sin(par[4]) * par[1]   +xyzStart[0];
+    vect[1] = -cos(par[4]) * par[1]   +xyzStart[1];
+    vect[2] = par[2]                  +xyzStart[2];
+
+    myMagFld.getMagFld( vect[0],vect[1],vect[2],fx,fy,fz);
+    constB = fz * myMagFld.getCnvCst();
+
+    pt = constB / fabs(par[5]);
+    px = pt * cos(par[4]);
+    py = pt * sin(par[4]);
+    pz = pt / tan(par[3]);
+    p = sqrt(pt * pt + pz * pz);
+
+    vect[3] = px / p;
+    vect[4] = py / p;
+    vect[5] = pz / p;
+    vect[6] = p;
+    charge = -d_sign( 1., par[5]);
+    poscur = 0.;
+//std::cout <<"VkGrkuta vect="<<vect[0]<<", "<<vect[1]<<", "<<vect[2]<<", "<<vect[3]<<", "
+//                            <<vect[4]<<", "<<vect[5]<<", "<<vect[6]<<'\n';
+    while(fabs(poscur) < fabs(totway)) {
+	posold = poscur;
+	d__1 = fabs(poscur) + stmg;
+        d__2 = fabs(totway);
+	poscur = d__1 < d__2 ? d__1 : d__2;
+	poscur = d_sign(poscur, totway);
+	dstep = poscur - posold;
+	vkgrkuta_(&charge, &dstep, vect, vout);
+	vect[0] = vout[0];
+	vect[1] = vout[1];
+	vect[2] = vout[2];
+	vect[3] = vout[3];
+	vect[4] = vout[4];
+	vect[5] = vout[5];
+// safety for strongly nonuniform field
+	dX = vect[0]-xyzEnd[0];
+	dY = vect[1]-xyzEnd[1];
+	dZ = vect[2]-xyzEnd[2];
+	dist_left = sqrt( dX*dX + dY*dY + dZ*dZ);
+	if(dist_left < stmg) break;    // arrived
+    }
+
+/* --- Now we are in the new point. Calculate track parameters */
+/*   at new point with new mag.field */
+
+    myMagFld.getMagFld( xyzEnd[0],xyzEnd[1],xyzEnd[2],fx,fy,fz);
+    //myMagFld.getMagFld( vect[0], vect[1], vect[2], fx, fy, fz);
+    constB = fz * myMagFld.getCnvCst();
+
+    dpar0[0] = 0.;
+    dpar0[1] = 0.;
+    dpar0[2] = acos(vout[5]);
+    dpar0[3] = atan2(vout[4], vout[3]);
+//    if (dpar0[3] < 0.) dpar0[3] += 6.2831853071794;
+
+    px = vect[3] * vect[6];
+    py = vect[4] * vect[6];
+    dpar0[4] = d_sign(  (constB/sqrt(px*px+py*py)), par[5]);    /* new value of curvature */
+    ich = (long int) charge;
+    xyzst[0] = xyzEnd[0] - vout[0];
+    xyzst[1] = xyzEnd[1] - vout[1];
+    xyzst[2] = xyzEnd[2] - vout[2];
+
+    cfnewp(&ich, dpar0, xyzst, &dp, parn, perig); // Last step of propagation 
+    closePoint[0] = perig[0] + vout[0];                  // with simple program
+    closePoint[1] = perig[1] + vout[1];
+    closePoint[2] = perig[2] + vout[2];
+
+//std::cout <<" Exit from cfNewPm="<<dp<<'\n';
+//std::cout <<parn[0]<<'\n';
+//   if(parn[0] == 0.) {
+//std::cout <<dp<<", "<<(*ustep)<<'\n';
+//       for (int jj=0; jj<6;jj++) std::cout <<" vout="<<vout[jj]<<'\n';
+//       for (int jj=0; jj<3;jj++) std::cout <<" xyzst="<<xyzst[jj]<<'\n';
+//       for (int jj=0; jj<5;jj++) std::cout <<" dpar0="<<dpar0[jj]<<'\n';
+//       for (int jj=0; jj<5;jj++) std::cout <<" parn="<<parn[jj]<<'\n';
+//   }
+    return;
+} 
+
+
+
+
+} /* End of namespace */

Tracking/TrkVertexFitter/TrkVKalVrtCore/src/cfVrtDst.cxx

@@ -32,7 +32,7 @@ extern vkalPropagator  myPropagator;
 //--------------------------------------------------------------------------------------------------
 double cfVrtDstSig( VKVertex * vk, bool UseTrkErr)
 {
-    long int it,IERR; int i,j,ij;
+    int i,j,ij,it;
     double parV0[5], covParV0[15];
     double Signif;
 
@@ -96,7 +96,7 @@ double cfVrtDstSig( VKVertex * vk, bool UseTrkErr)
     if ( UseTrkErr){ covImp[0] += nCov[0]; covImp[1] += nCov[1]; covImp[2] += nCov[2];}
 
     double dwgt[3];
-    IERR=cfdinv(covImp, dwgt, -2); if(IERR){ IERR=cfdinv(covImp, dwgt, 2); if(IERR){dwgt[0]=dwgt[2]=1.e6; dwgt[1]=0.; IERR=0; }}
+    int IERR=cfdinv(covImp, dwgt, -2); if(IERR){ IERR=cfdinv(covImp, dwgt, 2); if(IERR){dwgt[0]=dwgt[2]=1.e6; dwgt[1]=0.;}}
     Signif = sqrt(dwgt[0] * nPar[0] * nPar[0] +  2. * dwgt[1] * nPar[0] * nPar[1] + 
                   dwgt[2] * nPar[1] * nPar[1]);