Index: examples/ExamplePVtxFind.C
===================================================================
--- examples/ExamplePVtxFind.C	(revision bd376e3a1d8131d6e838b21d8cd738f47b325aad)
+++ examples/ExamplePVtxFind.C	(revision b750b0a794794673ed5b2468f10e215307d06a20)
@@ -1,5 +1,5 @@
 /*
 Example of using vertex fitter class to fit primary vertex
-assumed to be generated in (0,0,0)
+assumed to be generated with Pythia8/ee_zh_smr-shf.cmn
 */
 
@@ -44,4 +44,5 @@
 	TH1D* hTrFoundG = new TH1D("hTrFoundG", "Found GOOD primary tracks", 41, -0.5, 40.5);
 	TH1D* hTrFoundB = new TH1D("hTrFoundB", "Found BAD  primary tracks", 41, -0.5, 40.5);
+	TH1D* hTrDiff = new TH1D("hTrDiff", "Found - available tracks", 21, -10.5, 10.5);
 	//
 	// Loop over all events
@@ -87,4 +88,19 @@
 				pr[it] = new TVectorD(obsPar);
 				cv[it] = new TMatrixDSym(trk->CovarianceMatrix());
+				// Check covariance
+				/*
+				TMatrixDSymEigen Eign(*cv[it]);
+				TVectorD lambda = Eign.GetEigenValues();
+				Bool_t good = kTRUE;
+				for(Int_t i=0; i<5; i++) if(lambda(i)<0.0) good = kFALSE;
+				if(!good){
+					std::cout<<"Cov["<<it<<"]"<<" eigenvalues: "
+					<<lambda(0)<<", "<<lambda(1)<<", "<<lambda(2)<<", "
+					<<lambda(3)<<", "<<lambda(4)<<std::endl;
+					TMatrixDSym Ctmp = *cv[it];
+					std::cout<<"Error D (dir - cov) "<<trk->ErrorD0
+					<<" - "<<TMath::Sqrt(Ctmp(0,0))<<std::endl;
+				}
+				*/
 				//std::cout << "Track loaded it= " << it << std::endl;
 			//
@@ -120,4 +136,5 @@
 			}		// End loop on tracks
 		}
+		std::cout<<"PVtxFind true vertex: Nprim= "<<Nprim<<", x,y,z= "<<xpv<<", "<<ypv<<", "<<zpv<<std::endl;
 		//
 		// Find primary vertex
@@ -132,4 +149,9 @@
 		covpvc(1, 1) = 2.55e-05 * 2.55e-05;
 		covpvc(2, 2) = 0.64 * 0.64;
+		if(Nprim == 0){
+			xpv = xpvc(0);
+			ypv = xpvc(1);
+			zpv = xpvc(2);
+		}
 		//
 		//
@@ -137,9 +159,12 @@
 		Int_t nSkim = 0;
 		Int_t* nSkimmed = new Int_t[NtrG];
-		VertexFit* Vskim = new VertexFit();
-		Vskim->AddVtxConstraint(xpvc, covpvc);
-		Double_t MaxChi2 = 2.0+3*sqrt(2.);
+		TVectorD** PrSk = new TVectorD * [1];
+		TMatrixDSym** CvSk = new TMatrixDSym * [1];
+		Double_t MaxChi2 = 2.0+3*2.;
 		for (Int_t n = 0; n < NtrG; n++) {
-			Vskim->AddTrk(pr[n], cv[n]);
+			PrSk[0] = new TVectorD(*pr[n]);
+			CvSk[0] = new TMatrixDSym(*cv[n]);
+			VertexFit* Vskim = new VertexFit(1,PrSk, CvSk);
+			Vskim->AddVtxConstraint(xpvc, covpvc);
 			Double_t Chi2One = Vskim->GetVtxChi2();
 			//std::cout<<"Track "<<n<<", Chi2 = "<<Chi2One<<std::endl;
@@ -149,23 +174,25 @@
 				nSkim++;
 			}
-			Vskim->RemoveTrk(0);
-		}
+			// Cleanup
+			delete Vskim;
+		}
+		delete PrSk[0];
+		delete CvSk[0];
+		delete[] PrSk;
+		delete[] CvSk;
+		//
 		// Load tracks for primary fit
 		Int_t MinTrk = 1;	// Minumum # tracks for vertex fit
-		TVectorD** PrFit = new TVectorD * [nSkim];
-		TMatrixDSym** CvFit = new TMatrixDSym * [nSkim];
 		if (nSkim >= MinTrk) {
+			TVectorD** PrFit = new TVectorD * [nSkim];
+			TMatrixDSym** CvFit = new TMatrixDSym * [nSkim];
 			for (Int_t n = 0; n < nSkim; n++) {
 				PrFit[n] = new TVectorD(*pr[nSkimmed[n]]);
 				CvFit[n] = new TMatrixDSym(*cv[nSkimmed[n]]);
-				TVectorD test = *PrFit[n];
-				//std::cout << "Test *PrFit[" << n << "](0) = " << test(0) << std::endl;
-			}
-		}
-		delete[] nSkimmed;
-		delete Vskim;
-		Int_t Nfound = nSkim;
-		if(entry%500 ==0)std::cout << "Found tracks "<<Nfound << std::endl;
-		if (nSkim >= MinTrk) {
+			}
+			delete[] nSkimmed;
+			Int_t Nfound = nSkim;
+			if(entry%500 ==0)std::cout << "Found tracks "<<Nfound << std::endl;
+			const Int_t MaxFound = 100; Double_t Chi2LL[MaxFound]; Double_t *Chi2L = Chi2LL;
 			VertexFit* Vtx = new VertexFit(nSkim, PrFit, CvFit);
 			//std::cout << "Vertex fit created " << std::endl;
@@ -173,5 +200,5 @@
 			//
 			// Remove tracks with large chi2
-			Double_t MaxChi2Fit = 9.;
+			Double_t MaxChi2Fit = 4.5;
 			Bool_t Done = kFALSE;
 			while (!Done) {
@@ -180,5 +207,5 @@
 				TVectorD Chi2List = Vtx->GetVtxChi2List();	// Get contributions to Chi2
 				//std::cout << "After Chi2List.  " << std::endl; Chi2List.Print();
-				Double_t* Chi2L = new Double_t[Nfound];
+				//Double_t* Chi2L = new Double_t[Nfound];
 				Chi2L = Chi2List.GetMatrixArray();
 				Int_t iMax = TMath::LocMax(Nfound, Chi2L);
@@ -186,5 +213,5 @@
 				Double_t Chi2Mx = Chi2L[iMax];
 				//std::cout << "Chi2Mx "<<Chi2Mx << std::endl;
-				if (Chi2Mx > MaxChi2Fit && Nfound > 2) {
+				if (Chi2Mx > MaxChi2Fit && Nfound > 1) {
 					//std::cout << "Before remove.  Nfound = "<<Nfound << std::endl;
 					Vtx->RemoveTrk(iMax);
@@ -195,7 +222,4 @@
 					Done = kTRUE;
 				}
-				//std::cout << "Done =  " << Done << std::endl;
-				//delete[] Chi2L;
-				//std::cout << "Array Chi2L removed " << std::endl;
 			}
 			//
@@ -203,8 +227,8 @@
 			//
 			// Require minimum number of tracks in vertex
-			Int_t Nmin = 2;
+			Int_t Nmin = 1;
 			if (Nfound >= Nmin) {
 				TVectorD xvtx = Vtx->GetVtx();
-				//std::cout << "Found vertex " << xvtx(0)<<", "<<xvtx(1)<<", "<<xvtx(2) << std::endl;
+				std::cout << "Found vertex " << xvtx(0)<<", "<<xvtx(1)<<", "<<xvtx(2) << std::endl;
 				TMatrixDSym covX = Vtx->GetVtxCov();
 				Double_t Chi2 = Vtx->GetVtxChi2();
@@ -222,8 +246,14 @@
 				hTrPrim->Fill(Nprim);
 				hTrFound->Fill((Double_t)Nfound);
+				hTrDiff->Fill((Double_t)Nfound-Nprim);
 				//std::cout << "Histograms filled " << std::endl;
 			}
 			//
+			// Clean
 			delete Vtx;
+			for (Int_t i = 0; i < nSkim; i++) delete PrFit[i];
+			for (Int_t i = 0; i < nSkim; i++) delete CvFit[i];
+			delete[] PrFit;
+			delete[] CvFit;
 		}
 
@@ -233,10 +263,6 @@
 		for (Int_t i = 0; i < NtrG; i++) delete pr[i];
 		for (Int_t i = 0; i < NtrG; i++) delete cv[i];
-		for (Int_t i = 0; i < nSkim; i++) delete PrFit[i];
-		for (Int_t i = 0; i < nSkim; i++) delete CvFit[i];
 		delete[] pr;
 		delete[] cv;
-		delete[] PrFit;
-		delete[] CvFit;
 	}
 	//
@@ -254,5 +280,5 @@
 	//
 	TCanvas* CnvN = new TCanvas("CnvN", "Primary tracks found", 100, 100, 900, 500);
-	CnvN->Divide(2, 1);
+	CnvN->Divide(2, 2);
 	CnvN->cd(1);
 	hTrPrim->Draw();
@@ -260,3 +286,5 @@
 	hTrFound->SetLineColor(kRed);
 	hTrFound->Draw();
+	CnvN->cd(3);
+	hTrDiff->Draw();
 }
Index: examples/ExamplePVtxFitEC.C
===================================================================
--- examples/ExamplePVtxFitEC.C	(revision bd376e3a1d8131d6e838b21d8cd738f47b325aad)
+++ examples/ExamplePVtxFitEC.C	(revision b750b0a794794673ed5b2468f10e215307d06a20)
@@ -1,5 +1,5 @@
 /*
 Example of using vertex fitter class to fit primary vertex
-assumed to be generated in (0,0,0)
+assumed to be generated with Pythia8/ee_zh_smr-shf.cmn
 */
 
@@ -21,4 +21,20 @@
 void ExamplePVtxFitEC(const char* inputFile, Int_t Nevent = 5)
 {
+	//
+	// Beam constraint 
+	// (consistent with generation with ee_zh_smr-shf.cmnd)
+	//
+	// Mean beam position
+	TVectorD xpvc(3);
+	xpvc(0) = 1.0;
+	xpvc(1) = -2.0;
+	xpvc(2) = 10.0;
+	// Interaction region covariance
+	TMatrixDSym covpvc(3); covpvc.Zero();
+	covpvc(0, 0) = 0.0097 * 0.0097;
+	covpvc(1, 1) = 2.55e-05 * 2.55e-05;
+	covpvc(2, 2) = 0.64 * 0.64;
+
+	//
 	// Create chain of root trees
 	TChain chain("Delphes");
@@ -35,8 +51,14 @@
 	// Book histograms
 	Int_t Nbin = 100;
+	// Vertex fit pulls
 	TH1D* hXpull = new TH1D("hXpull", "Pull X vertex component", Nbin, -10., 10.);
 	TH1D* hYpull = new TH1D("hYpull", "Pull Y vertex component", Nbin, -10., 10.);
 	TH1D* hZpull = new TH1D("hZpull", "Pull Z vertex component", Nbin, -10., 10.);
 	TH1D* hChi2 = new TH1D("hChi2", "Vertex #chi^{2}/N_{dof}", Nbin, 0., 10.);
+	// Generation check
+	TH1D* hXvpull = new TH1D("hXvpull", "Pull X generated vertex", Nbin, -10., 10.);
+	TH1D* hYvpull = new TH1D("hYvpull", "Pull Y generated vertex", Nbin, -10., 10.);
+	TH1D* hZvpull = new TH1D("hZvpull", "Pull Z generated vertex", Nbin, -10., 10.);
+
 	//
 	// Loop over all events
@@ -65,10 +87,13 @@
 				// Get associated generated particle
 				GenParticle* gp = (GenParticle*)trk->Particle.GetObject();
-				TVector3 xg(1.e-3*gp->X,1.e-3*gp->Y,1.e-3*gp->Z);
+				TVector3 xg(1.e-3*gp->X,1.e-3*gp->Y,1.e-3*gp->Z); 	// mm -> meters
 				TVector3 pg(gp->Px,gp->Py,gp->Pz);
 				Double_t Q = (Double_t)gp->Charge;
 				Double_t Bz = 2.0;
 				TVectorD genParM =TrkUtil:: XPtoPar(xg, pg, Q, Bz);
-				TVectorD genPar = TrkUtil::ParToMm(genParM);
+				TVectorD genPar = TrkUtil::ParToMm(genParM); 		// -> back to mm
+
+				// 
+				// Check if original track is from primary vertex
 				//
 				// Position of origin in mm
@@ -79,20 +104,20 @@
 				Int_t mp = gp->M1;	// Mother
 				while(mp>0){
-					GenParticle* gm = 
-					(GenParticle*)branchGenPart->At(mp);
+					GenParticle* gm = (GenParticle*)branchGenPart->At(mp);
 					Double_t xm = gm->X;
                                 	Double_t ym = gm->Y;
                                 	Double_t zm = gm->Z;
 					if(x!=xm || y!=ym || z!=zm){
-						prim = kFALSE;
+						prim = kFALSE;	// Not primary
 						break;
 					}else mp = gm->M1;
 				}
+
 				//
 				// group tracks originating from the primary vertex
 				if (prim)
 				{
-					//std::cout<<"Event: "<<entry<<", track: "<<it;
-					//std::cout<<", x = "<<x<<", y = "<<y<<", z = "<<z<<std::endl; 
+					//
+					// Store true primary vertex for this event
 					xpv = x;
 					ypv = y;
@@ -104,7 +129,4 @@
 					Double_t obsC = trk->C;
 					Double_t obsZ0 = trk->DZ;
-					//std::cout<<"Z0 track = "<< obsZ0
-					//<<", gen Z0 = "<<genPar(3)
-					//<<", gen cot = "<<genPar(4)<<std::endl;
 					Double_t obsCtg = trk->CtgTheta;
 					Double_t oPar[5] = { obsD0, obsPhi, obsC, obsZ0, obsCtg };
@@ -112,6 +134,4 @@
 					//
 					pr[Ntr] = new TVectorD(obsPar);
-					//std::cout<<"Cov Matrix:"<<std::endl;
-					//trk->CovarianceMatrix().Print();
 					cv[Ntr] = new TMatrixDSym(trk->CovarianceMatrix());
 					Ntr++;
@@ -122,14 +142,5 @@
 		// Fit primary vertex with beam constraint
 		//
-		//Beam constraint
-		TVectorD xpvc(3);
-		xpvc(0) = 1.0;
-		xpvc(1) = -2.0;
-		xpvc(2) = 10.0;
-		TMatrixDSym covpvc(3); covpvc.Zero();
-		covpvc(0, 0) = 0.0097 * 0.0097;
-		covpvc(1, 1) = 2.55e-05 * 2.55e-05;
-		covpvc(2, 2) = 0.64 * 0.64;
-		Int_t MinTrk = 2;	// Minumum # tracks for vertex fit
+		Int_t MinTrk = 3;	// Minumum # tracks for vertex fit
 		if (Ntr >= MinTrk) {
 			VertexFit* Vtx = new VertexFit(Ntr, pr, cv);
@@ -139,11 +150,14 @@
 			Double_t Chi2 = Vtx->GetVtxChi2();
 			Double_t Ndof = 2 * (Double_t)Ntr;
+			delete Vtx;
+			//
 			Double_t PullX = (xvtx(0)-xpv) / TMath::Sqrt(covX(0, 0));
 			Double_t PullY = (xvtx(1)-ypv) / TMath::Sqrt(covX(1, 1));
 			Double_t PullZ = (xvtx(2)-zpv) / TMath::Sqrt(covX(2, 2));
-			//std::cout<<"**** True  vertex (x, y, z) "<<xpv<<", "
-			//<<ypv<<", "<<zpv<<std::endl;
-			//std::cout<<"**** Found vertex (x, y, z) "<<xvtx(0)<<", "
-                        //<<xvtx(1)<<", "<<xvtx(2)<<std::endl;
+			//
+			Double_t PullXv = (xpvc(0)-xpv) / TMath::Sqrt(covpvc(0, 0));
+			Double_t PullYv = (xpvc(1)-ypv) / TMath::Sqrt(covpvc(1, 1));
+			Double_t PullZv = (xpvc(2)-zpv) / TMath::Sqrt(covpvc(2, 2));
+			//
 			//
 			// Fill histograms
@@ -152,7 +166,10 @@
 			hZpull->Fill(PullZ);
 			hChi2->Fill(Chi2 / Ndof);
+			//
+			hXvpull->Fill(PullXv);
+			hYvpull->Fill(PullYv);
+			hZvpull->Fill(PullZv);
 		}
 
-		//std::cout << "Vertex chi2/Ndof = " << Chi2 / Ndof << std::endl;
 		//
 		// Cleanup
@@ -162,4 +179,5 @@
 		delete[] cv;
 	}
+
 	//
 	// Show resulting histograms
@@ -174,3 +192,12 @@
 	hZpull->Fit("gaus"); hZpull->Draw();
 	Cnv->cd(4); hChi2->Draw();
+	//
+	TCanvas* Cnv1 = new TCanvas("Cnv1", "Generated primary vertex pulls", 100, 100, 900, 500);
+	Cnv1->Divide(3, 1);
+	Cnv1->cd(1); gPad->SetLogy(1); gStyle->SetOptFit(1111);
+	hXvpull->Fit("gaus"); hXvpull->Draw();
+	Cnv1->cd(2); gPad->SetLogy(1); gStyle->SetOptFit(1111);
+	hYvpull->Fit("gaus"); hYvpull->Draw();
+	Cnv1->cd(3); gPad->SetLogy(1); gStyle->SetOptFit(1111);
+	hZvpull->Fit("gaus"); hZvpull->Draw();
 }
Index: external/TrackCovariance/TrkUtil.h
===================================================================
--- external/TrackCovariance/TrkUtil.h	(revision bd376e3a1d8131d6e838b21d8cd738f47b325aad)
+++ external/TrackCovariance/TrkUtil.h	(revision b750b0a794794673ed5b2468f10e215307d06a20)
@@ -68,15 +68,15 @@
 	//
 	static TVectorD XPtoPar(TVector3 x, TVector3 p, Double_t Q, Double_t Bz);
-	static TVector3 ParToX(TVectorD Par);				// position of minimum distance from z axis
+	static TVector3 ParToX(TVectorD Par);			// position of minimum distance from z axis
 	static TVector3 ParToP(TVectorD Par, Double_t Bz);	// Get Momentum from track parameters
-	static Double_t ParToQ(TVectorD Par);				// Get track charge
+	static Double_t ParToQ(TVectorD Par);			// Get track charge
 	static void LineDistance(TVector3 x0, TVector3 y0, TVector3 dirx, TVector3 diry, Double_t &sx, Double_t &sy, Double_t &distance);
 	//
 	// Track trajectory
 	//
-	static TVector3 Xtrack(TVectorD par, Double_t s);		// Parametric track trajectory
-	TVectorD derRphi_R(TVectorD par, Double_t R);	// Derivatives of R-phi at constant R
+	static TVector3 Xtrack(TVectorD par, Double_t s);	// Parametric track trajectory
+	TVectorD derRphi_R(TVectorD par, Double_t R);		// Derivatives of R-phi at constant R
 	TVectorD derZ_R(TVectorD par, Double_t R);		// Derivatives of z at constant R
-	TVectorD derRphi_Z(TVectorD par, Double_t z);	// Derivatives of R-phi at constant z
+	TVectorD derRphi_Z(TVectorD par, Double_t z);		// Derivatives of R-phi at constant z
 	TVectorD derR_Z(TVectorD par, Double_t z);		// Derivatives of R at constant z
 	//
@@ -88,5 +88,5 @@
 	//
 	static TVectorD ParToMm(TVectorD Par);			// Parameter conversion
-	static TMatrixDSym CovToMm(TMatrixDSym Cov);	// Covariance conversion
+	static TMatrixDSym CovToMm(TMatrixDSym Cov);		// Covariance conversion
 	//
 	// Inside cylindrical volume
Index: external/TrackCovariance/VertexFit.cc
===================================================================
--- external/TrackCovariance/VertexFit.cc	(revision bd376e3a1d8131d6e838b21d8cd738f47b325aad)
+++ external/TrackCovariance/VertexFit.cc	(revision b750b0a794794673ed5b2468f10e215307d06a20)
@@ -1,7 +1,11 @@
-﻿#include <TMath.h>
+/*
+Vertex fitting code
+*/
+#include <TMath.h>
 #include <TVectorD.h>
 #include <TVector3.h>
 #include <TMatrixD.h>
 #include <TMatrixDSym.h>
+#include <TMatrixDSymEigen.h>
 #include "VertexFit.h"
 //
@@ -38,10 +42,8 @@
 	for (Int_t i = 0; i < fNtr; i++)
 	{
-		TVectorD pr = *trkPar[i];
-		fPar.push_back(new TVectorD(pr));
-		fParNew.push_back(new TVectorD(pr));
-		TMatrixDSym cv = *trkCov[i];
-		fCov.push_back(new TMatrixDSym(cv));
-		fCovNew.push_back(new TMatrixDSym(cv));
+		fPar.push_back(new TVectorD(*trkPar[i]));
+		fParNew.push_back(new TVectorD(*trkPar[i]));
+		fCov.push_back(new TMatrixDSym(*trkCov[i]));
+		fCovNew.push_back(new TMatrixDSym(*trkCov[i]));
 	}
 	fChi2List.ResizeTo(fNtr);
@@ -76,26 +78,28 @@
 void VertexFit::ResetWrkArrays()
 {
-	Int_t N = (Int_t)ffi.size();
-	for (Int_t i = 0; i < N; i++)
-	{
-		if (fx0i[i])  { fx0i[i]->Clear();		delete fx0i[i]; }
-		if (fai[i])   { fai[i]->Clear();		delete fai[i]; }
-		if (fdi[i])   { fdi[i]->Clear();		delete fdi[i]; }
-		if (fAti[i])  { fAti[i]->Clear();		delete fAti[i]; }
-		if (fDi[i])   { fDi[i]->Clear();		delete fDi[i]; }
-		if (fWi[i])   { fWi[i]->Clear();		delete fWi[i]; }
-		if (fWinvi[i]){ fWinvi[i]->Clear();     delete fWinvi[i]; }
-	}
-	fa2i.clear();
-	fx0i.clear();
-	fai.clear();
-	fdi.clear();
-	fAti.clear();
-	fDi.clear();
-	fWi.clear();
-	fWinvi.clear();
+	Int_t N = (Int_t)fdi.size();
+	if(N > 0){
+		for (Int_t i = 0; i < N; i++)
+		{
+			if (fx0i[i])  { fx0i[i]->Clear();	delete fx0i[i]; }
+			if (fai[i])   { fai[i]->Clear();	delete fai[i]; }
+			if (fdi[i])   { fdi[i]->Clear();	delete fdi[i]; }
+			if (fAti[i])  { fAti[i]->Clear();	delete fAti[i]; }
+			if (fDi[i])   { fDi[i]->Clear();	delete fDi[i]; }
+			if (fWi[i])   { fWi[i]->Clear();	delete fWi[i]; }
+			if (fWinvi[i]){ fWinvi[i]->Clear();     delete fWinvi[i]; }
+		}
+		fa2i.clear();
+		fx0i.clear();
+		fai.clear();
+		fdi.clear();
+		fAti.clear();
+		fDi.clear();
+		fWi.clear();
+		fWinvi.clear();
+	}
 }
 VertexFit::~VertexFit()
-{
+{	
 	fxCst.Clear();		
 	fCovCst.Clear();
@@ -103,5 +107,5 @@
 	fXv.Clear();		
 	fcovXv.Clear();		
-	fChi2List.Clear();	
+	fChi2List.Clear();
 	//
 	for (Int_t i = 0; i < fNtr; i++)
@@ -111,5 +115,5 @@
 		fCov[i]->Clear();		delete fCov[i];
 		fCovNew[i]->Clear();	delete fCovNew[i];
-	}
+	}	
 	fPar.clear();
 	fParNew.clear();
@@ -118,218 +122,16 @@
 	//
 	ResetWrkArrays();
-	ffi.clear();
+	ffi.clear();	
 	fNtr = 0;
 }
 //
-Double_t VertexFit::FastRv(TVectorD p1, TVectorD p2)
-{
-	//
-	// Find radius of minimum distance between two tracks
-	//
-	// p = (D,phi, C, z0, ct)
-	//
-	// Define arrays
-	//
-	Double_t C1 = p1(2);
-	Double_t C2 = p2(2);
-	Double_t ph1 = p1(1);
-	Double_t ph2 = p2(1);
-	TVectorD x0 = Fill_x0(p1);
-	TVectorD y0 = Fill_x0(p2);
-	TVectorD n = Fill_a(p1, 0.0);
-	n *= (2*C1);
-	TVectorD k = Fill_a(p2, 0.0);
-	k *= (2*C2);
-	//
-	// Setup and solve linear system
-	//
-	Double_t nn = 0; for (Int_t i = 0; i < 3; i++)nn += n(i) * n(i);
-	Double_t nk = 0; for (Int_t i = 0; i < 3; i++)nk += n(i) * k(i);
-	Double_t kk = 0; for (Int_t i = 0; i < 3; i++)kk += k(i) * k(i);
-	Double_t discr = nn * kk - nk * nk;
-	TMatrixDSym H(2);
-	H(0, 0) = kk;
-	H(0, 1) = nk;
-	H(1, 0) = H(0, 1);
-	H(1, 1) = nn;
-	TVectorD c(2);
-	c(0) = 0; for (Int_t i = 0; i < 3; i++)c(0) += n(i) * (y0(i) - x0(i));
-	c(1) = 0; for (Int_t i = 0; i < 3; i++)c(1) += -k(i) * (y0(i) - x0(i));
-	TVectorD smin = (H * c);
-	smin *= 1.0 / discr;
-	//
-	TVectorD X = x0 + smin(0) * n;
-	TVectorD Y = y0 + smin(1) * k;
-	//
-	// Higher order corrections
-	X(0) += -C1 * smin(0) * smin(0) * TMath::Sin(ph1);
-	X(1) +=  C1 * smin(0) * smin(0) * TMath::Cos(ph1);
-	Y(0) += -C2 * smin(1) * smin(1) * TMath::Sin(ph2);
-	Y(1) +=  C2 * smin(1) * smin(1) * TMath::Cos(ph2);
-	//
-	TVectorD Xavg = 0.5 * (X + Y);
-	//
-	//
-	return TMath::Sqrt(Xavg(0)*Xavg(0)+Xavg(1)*Xavg(1));
-}
-//
-// Starting radius determination
-Double_t VertexFit::StartRadius()
-{
-	//
-	// Maximum impact parameter
-	Double_t Rd = 0;
-	for (Int_t i = 0; i < fNtr; i++)
-	{
-		TVectorD par = *fPar[i];
-		Double_t Dabs = TMath::Abs(par(0));
-		if (Dabs > Rd)Rd = Dabs;
-	}
-	//-----------------------------
-	//
-	// Find track pair with phi difference closest to pi/2
-	Int_t isel = 0; Int_t jsel = 0;		// selected track indices
-	Double_t dSinMax = 0.0;				// Max phi difference
-	for (Int_t i = 0; i < fNtr - 1; i++)
-	{
-		TVectorD pari = *fPar[i];
-		Double_t phi1 = pari(1);
-
-		for (Int_t j = i + 1; j < fNtr; j++)
-		{
-			TVectorD parj = *fPar[j];
-			Double_t phi2 = parj(1);
-			Double_t Sindphi = TMath::Abs(TMath::Sin(phi2 - phi1));
-			if (Sindphi > dSinMax)
-			{
-				isel = i; jsel = j;
-				dSinMax = Sindphi;
-			}
-		}
-	}
-	//
-	//------------------------------------------
-	//
-	// Find radius of minimum distrance between tracks
-	TVectorD p1 = *fPar[isel];
-	TVectorD p2 = *fPar[jsel];
-	Double_t R = FastRv(p1, p2);
-	//
-	R = 0.9 * R + 0.1 * Rd;		// Protect for overshoot
-	//
-	return R;
-}
-//
-// Regularized symmetric matrix inversion
-//
-TMatrixDSym VertexFit::RegInv(TMatrixDSym& Min)
-{
-	TMatrixDSym M = Min;				// Decouple from input
-	Int_t N = M.GetNrows();			// Matrix size
-	TMatrixDSym D(N); D.Zero();		// Normaliztion matrix
-	TMatrixDSym R(N);				// Normarized matrix
-	TMatrixDSym Rinv(N);				// Inverse of R
-	TMatrixDSym Minv(N);				// Inverse of M
-	//
-	// Check for 0's and normalize
-	for (Int_t i = 0; i < N; i++)
-	{
-		if (M(i, i) != 0.0) D(i, i) = 1. / TMath::Sqrt(TMath::Abs(M(i, i)));
-		else D(i, i) = 1.0;
-	}
-	R = M.Similarity(D);
-	//
-	// Recursive algorithms stops when N = 2
-	//
-	//****************
-	// case N = 2  ***
-	//****************
-	if (N == 2)
-	{
-		Double_t det = R(0, 0) * R(1, 1) - R(0, 1) * R(1, 0);
-		if (det == 0)
-		{
-			std::cout << "VertexFit::RegInv: null determinant for N = 2" << std::endl;
-			Rinv.Zero();	// Return null matrix
-		}
-		else
-		{
-			// invert matrix 
-			Rinv(0, 0) = R(1, 1);
-			Rinv(0, 1) = -R(0, 1);
-			Rinv(1, 0) = Rinv(0, 1);
-			Rinv(1, 1) = R(0, 0);
-			Rinv *= 1. / det;
-		}
-	}
-	//****************
-	// case N > 2  ***
-	//****************
-	else
-	{
-		// Break up matrix
-		TMatrixDSym Q = R.GetSub(0, N - 2, 0, N - 2);	// Upper left 
-		TVectorD p(N - 1);
-		for (Int_t i = 0; i < N - 1; i++)p(i) = R(N - 1, i);
-		Double_t q = R(N - 1, N - 1);
-		//Invert pieces and re-assemble
-		TMatrixDSym Ainv(N - 1);
-		TMatrixDSym A(N - 1);
-		if (TMath::Abs(q) > 1.0e-15)
-		{
-			// Case |q| > 0
-			Ainv.Rank1Update(p, -1.0 / q);
-			Ainv += Q;
-			A = RegInv(Ainv);		// Recursive call
-			TMatrixDSub(Rinv, 0, N - 2, 0, N - 2) = A;
-			//
-			TVectorD b = (-1.0 / q) * (A * p);
-			for (Int_t i = 0; i < N - 1; i++)
-			{
-				Rinv(N - 1, i) = b(i);
-				Rinv(i, N - 1) = b(i);
-			}
-			//
-			Double_t pdotb = 0.;
-			for (Int_t i = 0; i < N - 1; i++)pdotb += p(i) * b(i);
-			Double_t c = (1.0 - pdotb) / q;
-			Rinv(N - 1, N - 1) = c;
-		}
-		else
-		{
-			// case q = 0
-			TMatrixDSym Qinv = RegInv(Q);		// Recursive call
-			Double_t a = Qinv.Similarity(p);
-			Double_t c = -1.0 / a;
-			Rinv(N - 1, N - 1) = c;
-			//
-			TVectorD b = (1.0 / a) * (Qinv * p);
-			for (Int_t i = 0; i < N - 1; i++)
-			{
-				Rinv(N - 1, i) = b(i);
-				Rinv(i, N - 1) = b(i);
-			}
-			//
-			A.Rank1Update(p, -1 / a);
-			A += Q;
-			A.Similarity(Qinv);
-			TMatrixDSub(Rinv, 0, N - 2, 0, N - 2) = A;
-		}
-	}
-	Minv = Rinv.Similarity(D);
-	return Minv;
-}
-//
-//
-//
-TMatrixD VertexFit::Fill_A(TVectorD par, Double_t phi)
-{
-	//
-	// Derivative of track 3D position vector with respect to track parameters at constant phase 
-	//
-	// par = vector of track parameters
-	// phi = phase
-	//
-	TMatrixD A(3, 5);
+//
+//
+TVectorD VertexFit::Fill_x0(TVectorD par)
+{
+	//
+	// Calculate track 3D position at R = |D| (minimum approach to z-axis)
+	//
+	TVectorD x0(3);
 	//
 	// Decode input arrays
@@ -341,38 +143,17 @@
 	Double_t ct = par(4);
 	//
-	// Fill derivative matrix dx/d alpha
-	// D
-	A(0, 0) = -TMath::Sin(p0);
-	A(1, 0) = TMath::Cos(p0);
-	A(2, 0) = 0.0;
-	// phi0
-	A(0, 1) = -D * TMath::Cos(p0) + (TMath::Cos(phi + p0) - TMath::Cos(p0)) / (2 * C);
-	A(1, 1) = -D * TMath::Sin(p0) + (TMath::Sin(phi + p0) - TMath::Sin(p0)) / (2 * C);
-	A(2, 1) = 0.0;
-	// C
-	A(0, 2) = -(TMath::Sin(phi + p0) - TMath::Sin(p0)) / (2 * C * C);
-	A(1, 2) = (TMath::Cos(phi + p0) - TMath::Cos(p0)) / (2 * C * C);
-	A(2, 2) = -ct * phi / (2 * C * C);
-	// z0
-	A(0, 3) = 0.0;
-	A(1, 3) = 0.0;
-	A(2, 3) = 1.0;
-	// ct = lambda
-	A(0, 4) = 0.0;
-	A(1, 4) = 0.0;
-	A(2, 4) = phi / (2 * C);
-	//
-	return A;
-}
-//
-TVectorD VertexFit::Fill_a(TVectorD par, Double_t phi)
-{
-	//
-	// Derivative of track 3D position vector with respect to phase at constant track parameters
-	//
-	// par = vector of track parameters
-	// phi = phase
-	//
-	TVectorD a(3);
+	x0(0) = -D * TMath::Sin(p0);
+	x0(1) = D * TMath::Cos(p0);
+	x0(2) = z0;
+	//
+	return x0;
+}
+//
+TVectorD VertexFit::Fill_x(TVectorD par, Double_t phi)
+{
+	//
+	// Calculate track 3D position for a given phase, phi
+	//
+	TVectorD x(3);
 	//
 	// Decode input arrays
@@ -384,48 +165,4 @@
 	Double_t ct = par(4);
 	//
-	a(0) = TMath::Cos(phi + p0) / (2 * C);
-	a(1) = TMath::Sin(phi + p0) / (2 * C);
-	a(2) = ct / (2 * C);
-	//
-	return a;
-}
-//
-TVectorD VertexFit::Fill_x0(TVectorD par)
-{
-	//
-	// Calculate track 3D position at R = |D| (minimum approach to z-axis)
-	//
-	TVectorD x0(3);
-	//
-	// Decode input arrays
-	//
-	Double_t D = par(0);
-	Double_t p0 = par(1);
-	Double_t C = par(2);
-	Double_t z0 = par(3);
-	Double_t ct = par(4);
-	//
-	x0(0) = -D * TMath::Sin(p0);
-	x0(1) = D * TMath::Cos(p0);
-	x0(2) = z0;
-	//
-	return x0;
-}
-//
-TVectorD VertexFit::Fill_x(TVectorD par, Double_t phi)
-{
-	//
-	// Calculate track 3D position for a given phase, phi
-	//
-	TVectorD x(3);
-	//
-	// Decode input arrays
-	//
-	Double_t D = par(0);
-	Double_t p0 = par(1);
-	Double_t C = par(2);
-	Double_t z0 = par(3);
-	Double_t ct = par(4);
-	//
 	TVectorD x0 = Fill_x0(par);
 	x(0) = x0(0) + (TMath::Sin(phi + p0) - TMath::Sin(p0)) / (2 * C);
@@ -439,25 +176,28 @@
 {
 	//
-	// Get track parameters and their covariance
+	// Get track parameters, covariance and phase
+	Double_t fs = ffi[i];			// Get phase
 	TVectorD par = *fParNew[i];
 	TMatrixDSym Cov = *fCov[i];
 	//
 	// Fill all track related work arrays arrays
-	Double_t fs = ffi[i];						// Get phase
+	TMatrixD A = derXdPar(par, fs);		// A = dx/da = derivatives wrt track parameters
+	TMatrixDSym Winv = Cov;				
+	Winv.Similarity(A);			// W^-1 = A*C*A'
+
+	TMatrixD At(TMatrixD::kTransposed, A);		// A transposed
+	fAti.push_back(new TMatrixD(At));		// Store A'
+	fWinvi.push_back(new TMatrixDSym(Winv));	// Store W^-1 matrix
+	//
 	TVectorD xs = Fill_x(par, fs);
 	fx0i.push_back(new TVectorD(xs));			// Start helix position
-	//
-	TMatrixD A = Fill_A(par, fs);				// A = dx/da = derivatives wrt track parameters
-	TMatrixD At(TMatrixD::kTransposed, A);		// A transposed
-	fAti.push_back(new TMatrixD(At));			// Store A' 
-	TVectorD di = A * (par - *fPar[i]);		// x-shift
-	fdi.push_back(new TVectorD(di));			// Store x-shift
-	TMatrixDSym Winv = Cov.Similarity(A);		// W^-1 = A*C*A'
-	fWinvi.push_back(new TMatrixDSym(Winv));	// Store W^-1 matrix
+	// 
+	TVectorD di = A * (par - *fPar[i]);			// x-shift
+	fdi.push_back(new TVectorD(di));			// Store x-shift	
 	//
 	TMatrixDSym W = RegInv(Winv);				// W = (A*C*A')^-1
 	fWi.push_back(new TMatrixDSym(W));			// Store W matrix
 	//
-	TVectorD a = Fill_a(par, fs);				// a = dx/ds = derivatives wrt phase
+	TVectorD a = derXds(par, fs);				// a = dx/ds = derivatives wrt phase
 	fai.push_back(new TVectorD(a));				// Store a
 	//
@@ -473,4 +213,73 @@
 }
 //
+void VertexFit::VtxFitNoSteer()
+{
+	//
+	// Initialize
+	//
+	std::vector<TVectorD*> x0i;				// Tracks at ma
+	std::vector<TVectorD*> ni;				// Track derivative wrt phase
+	std::vector<TMatrixDSym*> Ci;			// Position error matrix at fixed phase
+	std::vector<TVectorD*> wi;				// Ci*ni
+	//
+	// Track loop
+	for (Int_t i = 0; i < fNtr; i++)
+	{
+		Double_t s = 0.;
+		TVectorD par = *fPar[i];
+		TMatrixDSym Cov = *fCov[i];
+		x0i.push_back(new TVectorD(Fill_x0(par)));
+		ni.push_back(new TVectorD(derXds(par, s)));
+		TMatrixD A = derXdPar(par, s);
+		Ci.push_back(new TMatrixDSym(Cov.Similarity(A)));
+		TMatrixDSym Cinv = RegInv(*Ci[i]);
+		wi.push_back(new TVectorD(Cinv * (*ni[i])));
+	}
+	//std::cout << "Vtx init completed. fNtr = "<<fNtr << std::endl;
+	//
+	// Get fit vertex
+	//
+	TMatrixDSym D(3); D.Zero();
+	TVectorD Dx(3); Dx.Zero();
+	for (Int_t i = 0; i < fNtr; i++)
+	{
+		TMatrixDSym Cinv = RegInv(*Ci[i]);
+		TMatrixDSym W(3);
+		W.Rank1Update(*wi[i], 1. / Ci[i]->Similarity(*wi[i]));
+		TMatrixDSym Dd = Cinv - W;
+		D += Dd;
+		Dx += Dd * (*x0i[i]);
+	}
+	if(fVtxCst){
+		D  += fCovCstInv;
+		Dx += fCovCstInv*fxCst;
+	}
+	fXv = RegInv(D) * Dx;
+	//std::cout << "Fast vertex (x, y, z) = "<<fXv(0)<<", "<<fXv(1)<<", "<<fXv(2) << std::endl;
+	//
+	// Get fit phases
+	//
+	for (Int_t i = 0; i < fNtr; i++){
+		Double_t si = Dot(*wi[i], fXv - (*x0i[i])) / Ci[i]->Similarity(*wi[i]);
+		ffi.push_back(si);
+		//TVectorD xvi = Fill_x(*fPar[i],si);
+		//std::cout << "Fast vertex "<<i<<": xvi = "<<xvi(0)<<", "<<xvi(1)<<", "<<xvi(2) 
+		//	<<", si = "<<si<< std::endl;
+	}
+	//
+	// Cleanup
+	for (Int_t i = 0; i < fNtr; i++)
+	{
+		x0i[i]->Clear();	delete x0i[i];
+		ni[i]->Clear();		delete ni[i];
+		Ci[i]->Clear();		delete Ci[i];
+		wi[i]->Clear();		delete wi[i];
+	}
+	x0i.clear();;
+	ni.clear();
+	Ci.clear(); 
+	wi.clear();
+}
+//
 void  VertexFit::VertexFitter()
 {
@@ -487,17 +296,7 @@
 	TMatrixDSym covX(3);
 	Double_t Chi2 = 0;
-	TVectorD x0 = fXv;	// If previous fit done
-	if (fRold < 0.0) {
-		// External constraint
-		if (fVtxCst) fRold = TMath::Sqrt(fxCst(0) * fxCst(0) + fxCst(1) * fxCst(1));
-		// No constraint
-		else for (Int_t i = 0; i < 3; i++)x0(i) = 1000.;	// Set to arbitrary large value 
-	}
-	//
-	// Starting vertex radius approximation
-	//
-	Double_t R = fRold;						// Use previous fit if available
-	if (R < 0.0) R = StartRadius();			// Rough vertex estimate
-	//std::cout << "Start radius: " << R << std::endl;
+	//
+	VtxFitNoSteer();	// Fast vertex finder on first pass (set ffi and fXv)
+	TVectorD x0 = fXv;
 	//
 	// Iteration properties
@@ -515,5 +314,5 @@
 		TVectorD cterm(3); TMatrixDSym H(3); TMatrixDSym DW1D(3);
 		covX.Zero();		// Reset vertex covariance
-		cterm.Zero();	// Reset constant term
+		cterm.Zero();		// Reset constant term
 		H.Zero();		// Reset H matrix
 		DW1D.Zero();
@@ -530,15 +329,4 @@
 			TVectorD par = *fPar[i];
 			TMatrixDSym Cov = *fCov[i];
-			//
-			// For first iteration only
-			Double_t fs;
-			if (Ntry <= 0)	// Initialize all phases on first pass
-			{
-				Double_t D = par(0);
-				Double_t C = par(2);
-				Double_t arg = TMath::Max(1.0e-6, (R * R - D * D) / (1 + 2 * C * D));
-				fs = 2 * TMath::ASin(C * TMath::Sqrt(arg));
-				ffi.push_back(fs);
-			}
 			//
 			// Update track related arrays
@@ -555,10 +343,10 @@
 			// update constant term
 			TVectorD xs = *fx0i[i] - *fdi[i];
-			TVectorD xx0 = *fx0i[i];
-			/*
-			std::cout << "Iter. " << Ntry << ", trk " << i << ", x= "
-				<< xx0(0) << ", " << xx0(1) << ", " << xx0(2)<<
-				", ph0= "<<par(1)<< std::endl;
-			*/
+			//TVectorD xx0 = *fx0i[i];
+			
+			//std::cout << "Iter. " << Ntry << ", trk " << i << ", xs= "
+			//	<< xs(0) << ", " << xs(1) << ", " << xs(2)<<
+			//	", ph0= "<<par(1)<< std::endl;
+			
 			cterm += Ds * xs;
 		}				// End loop on tracks
@@ -584,8 +372,12 @@
 			TMatrixDSym Wm1 = *fWinvi[i];
 			fChi2List(i) = Wm1.Similarity(lambda);
+			if(fChi2List(i) < 0.0){
+				//std::cout<<"# "<<i<<", Chi2= "<<fChi2List(i)<<", Wm1:"<<std::endl; Wm1.Print();
+				//std::cout<<"Lambda= "<<std::endl; lambda.Print();
+			}
 			Chi2 += fChi2List(i);
 			TVectorD a = *fai[i];
-			TVectorD b = (*fWi[i]) * (x - (*fx0i[i]));
-			for (Int_t j = 0; j < 3; j++)ffi[i] += a(j) * b(j) / fa2i[i];
+			TVectorD b = (*fWi[i]) * (x - *fx0i[i] + *fdi[i]);
+			ffi[i] += Dot(a, b) / fa2i[i];
 			TVectorD newPar = *fPar[i] - ((*fCov[i]) * (*fAti[i])) * lambda;
 			fParNew[i] = new TVectorD(newPar);
@@ -610,5 +402,7 @@
 	//
 	fVtxDone = kTRUE;		// Set fit completion flag
-	fRold = TMath::Sqrt(fXv(0)*fXv(0) + fXv(1)*fXv(1));	// Store fit radius
+	fRold = TMath::Sqrt(fXv(0)*fXv(0) + fXv(1)*fXv(1));	// Store fit 
+	//std::cout << "Found vertex " << fXv(0) << ", " << fXv(1) << ", " << fXv(2) 
+	//	<< ", after "<<Ntry<<" iterations"<<std::endl;
 	//
 }
Index: external/TrackCovariance/VertexFit.h
===================================================================
--- external/TrackCovariance/VertexFit.h	(revision bd376e3a1d8131d6e838b21d8cd738f47b325aad)
+++ external/TrackCovariance/VertexFit.h	(revision b750b0a794794673ed5b2468f10e215307d06a20)
@@ -6,4 +6,5 @@
 #include <TVectorD.h>
 #include <TMatrixDSym.h>
+#include "TrkUtil.h"
 #include "ObsTrk.h"
 #include <vector>
@@ -12,5 +13,6 @@
 // Class for vertex fitting
 
-class VertexFit {
+class VertexFit: public TrkUtil
+{
 	//
 	// Vertex fitting with track parameters steering
@@ -21,5 +23,5 @@
 	//
 	// Inputs
-	Int_t fNtr;					// Number of tracks
+	Int_t fNtr;							// Number of tracks
 	std::vector<TVectorD*> fPar;		// Input parameter array
 	std::vector<TVectorD*> fParNew;		// Updated parameter array
@@ -27,11 +29,11 @@
 	std::vector<TMatrixDSym*> fCovNew;	// Updated parameter covariances
 	// Constraints
-	Bool_t fVtxCst;				// Vertex constraint flag
-	TVectorD fxCst;				// Constraint value
+	Bool_t fVtxCst;					// Vertex constraint flag
+	TVectorD fxCst;					// Constraint value
 	TMatrixDSym fCovCst;			// Constraint 
 	TMatrixDSym fCovCstInv;			// Inverse of constraint covariance
 	//
 	// Results
-	Bool_t fVtxDone;				// Flag vertex fit completed
+	Bool_t fVtxDone;			// Flag vertex fit completed
 	Double_t fRold;				// Current value of vertex radius
 	TVectorD fXv;				// Found vertex
@@ -52,15 +54,15 @@
 	//
 	// Service routines
-	//void InitWrkArrays();							// Initializations
 	void ResetWrkArrays();							// Clear work arrays
-	Double_t StartRadius();							// Starting vertex radius determination
-	Double_t FastRv(TVectorD p1, TVectorD p2);		// Fast vertex radius determination
-	TMatrixDSym RegInv(TMatrixDSym& Smat0);			// Regularized 3D matrix inversion 
-	TMatrixD Fill_A(TVectorD par, Double_t phi);		// Derivative of track position wrt track parameters
-	TVectorD Fill_a(TVectorD par, Double_t phi);		// Derivative of track position wrt track phase
+	//Double_t StartRadius();							// Starting vertex radius determination
+	//Double_t FastRv(TVectorD p1, TVectorD p2);		// Fast vertex radius determination
+	//TMatrixDSym RegInv(TMatrixDSym& Smat0);		// Regularized 3D matrix inversion 
+	//TMatrixD Fill_A(TVectorD par, Double_t phi);	// Derivative of track position wrt track parameters
+	//TVectorD Fill_a(TVectorD par, Double_t phi);	// Derivative of track position wrt track phase
 	TVectorD Fill_x0(TVectorD par);					// Track position at dma to z-axis
-	TVectorD Fill_x(TVectorD par, Double_t phi);		// Track position at given phase
+	TVectorD Fill_x(TVectorD par, Double_t phi);	// Track position at given phase
 	void UpdateTrkArrays(Int_t i);					// Fill track realted arrays
-	void VertexFitter();								// Vertex finder routine
+	void VtxFitNoSteer();							// Vertex fitter routine w/o parameter steering
+	void VertexFitter();							// Vertex fitter routine w/  parameter steering
 public:
 	//
Index: modules/TreeWriter.cc
===================================================================
--- modules/TreeWriter.cc	(revision bd376e3a1d8131d6e838b21d8cd738f47b325aad)
+++ modules/TreeWriter.cc	(revision b750b0a794794673ed5b2468f10e215307d06a20)
@@ -369,14 +369,14 @@
 
     // add some offdiagonal covariance matrix elements
-    entry->ErrorD0Phi          = candidate->TrackCovariance(0,1);
+    entry->ErrorD0Phi          = candidate->TrackCovariance(0,1)*1.e3;
     entry->ErrorD0C            = candidate->TrackCovariance(0,2);
-    entry->ErrorD0DZ           = candidate->TrackCovariance(0,3);
-    entry->ErrorD0CtgTheta     = candidate->TrackCovariance(0,4);
-    entry->ErrorPhiC           = candidate->TrackCovariance(1,2);
-    entry->ErrorPhiDZ          = candidate->TrackCovariance(1,3);
+    entry->ErrorD0DZ           = candidate->TrackCovariance(0,3)*1.e6;
+    entry->ErrorD0CtgTheta     = candidate->TrackCovariance(0,4)*1.e3;
+    entry->ErrorPhiC           = candidate->TrackCovariance(1,2)*1.e-3;
+    entry->ErrorPhiDZ          = candidate->TrackCovariance(1,3)*1.e3;
     entry->ErrorPhiCtgTheta    = candidate->TrackCovariance(1,4);
     entry->ErrorCDZ            = candidate->TrackCovariance(2,3);
-    entry->ErrorCCtgTheta      = candidate->TrackCovariance(2,4);
-    entry->ErrorDZCtgTheta     = candidate->TrackCovariance(3,4);
+    entry->ErrorCCtgTheta      = candidate->TrackCovariance(2,4)*1.e-3;
+    entry->ErrorDZCtgTheta     = candidate->TrackCovariance(3,4)*1.e3;
 
     entry->Xd = candidate->Xd;
