Index: external/TrackCovariance/ObsTrk.cc
===================================================================
--- external/TrackCovariance/ObsTrk.cc	(revision a95da74cda6d94acac311904287f73ef752d400f)
+++ external/TrackCovariance/ObsTrk.cc	(revision 53f474642cbcbb6dae79428731080cbdabba2e1d)
@@ -14,5 +14,5 @@
 ObsTrk::ObsTrk(TVector3 x, TVector3 p, Double_t Q, Double_t B, SolGridCov *GC)
 {
-	SetBfield(B);
+	SetB(B);
 	fGC = GC;
 	fGenX = x;
@@ -56,5 +56,5 @@
 ObsTrk::ObsTrk(Double_t *x, Double_t *p, Double_t Q, Double_t B, SolGridCov* GC)
 {
-	SetBfield(B);
+	SetB(B);
 	fGC = GC;
 	fGenX.SetXYZ(x[0],x[1],x[2]);
Index: external/TrackCovariance/TrkUtil.cc
===================================================================
--- external/TrackCovariance/TrkUtil.cc	(revision a95da74cda6d94acac311904287f73ef752d400f)
+++ external/TrackCovariance/TrkUtil.cc	(revision 53f474642cbcbb6dae79428731080cbdabba2e1d)
@@ -2,4 +2,5 @@
 #include <iostream>
 #include <algorithm>
+#include <TSpline.h>
 
 // Constructor
@@ -96,5 +97,5 @@
 	if (fBz == 0.0)std::cout << "TrkUtil::ParToP: Warning Bz not set" << std::endl;
 	//
-	return ParToP(Par,fBz);
+	return ParToP(Par, fBz);
 }
 //
@@ -233,11 +234,4 @@
 }
 //
-
-//
-void TrkUtil::SetBfield(Double_t Bz)
-{
-	fBz = Bz;
-}
-
 // Setup chamber volume
 void TrkUtil::SetDchBoundaries(Double_t Rmin, Double_t Rmax, Double_t Zmin, Double_t Zmax)
@@ -274,29 +268,29 @@
 		//	<< ", C= " << C << ", z0= " << z0 << ", ct= " << ct << std::endl;
 		//
-		// Track length per unit phase change
-		Double_t Scale = sqrt(1.0 + ct*ct) / (2.0*TMath::Abs(C));
+		// Track length per unit phase change 
+		Double_t Scale = sqrt(1.0 + ct * ct) / (2.0 * TMath::Abs(C));
 		//
 		// Find intersections with chamber boundaries
 		//
-		Double_t phRin = 0.0;			// phase of inner cylinder
-		Double_t phRin2= 0.0;			// phase of inner cylinder intersection (2nd branch)
+		Double_t phRin = 0.0;			// phase of inner cylinder 
+		Double_t phRin2 = 0.0;			// phase of inner cylinder intersection (2nd branch)
 		Double_t phRhi = 0.0;			// phase of outer cylinder intersection
 		Double_t phZmn = 0.0;			// phase of left wall intersection
 		Double_t phZmx = 0.0;			// phase of right wall intersection
 		//  ... with inner cylinder
-		Double_t Rtop = TMath::Abs((1.0 + C*D) / C);
+		Double_t Rtop = TMath::Abs((1.0 + C * D) / C);
 
 		if (Rtop > fRmin && TMath::Abs(D) < fRmin) // *** don't treat large D tracks for the moment ***
 		{
-			Double_t ph = 2 * asin(C*sqrt((fRmin*fRmin - D*D) / (1.0 + 2.0*C*D)));
-			Double_t z = z0 + ct*ph / (2.0*C);
+			Double_t ph = 2 * asin(C * sqrt((fRmin * fRmin - D * D) / (1.0 + 2.0 * C * D)));
+			Double_t z = z0 + ct * ph / (2.0 * C);
 
 			//std::cout << "Rin intersection: ph = " << ph<<", z= "<<z << std::endl;
 
-			if (z < fZmax && z > fZmin)	phRin = TMath::Abs(ph);	// Intersection inside chamber volume
+			if (z < fZmax && z > fZmin)	phRin = TMath::Abs(ph);	// Intersection inside chamber volume	
 			//
 			// Include second branch of loopers
 			Double_t Pi = 3.14159265358979323846;
-			Double_t ph2 = 2*Pi - TMath::Abs(ph);
+			Double_t ph2 = 2 * Pi - TMath::Abs(ph);
 			if (ph < 0)ph2 = -ph2;
 			z = z0 + ct * ph2 / (2.0 * C);
@@ -306,7 +300,7 @@
 		if (Rtop > fRmax && TMath::Abs(D) < fRmax) // *** don't treat large D tracks for the moment ***
 		{
-			Double_t ph = 2 * asin(C*sqrt((fRmax*fRmax - D*D) / (1.0 + 2.0*C*D)));
-			Double_t z = z0 + ct*ph / (2.0*C);
-			if (z < fZmax && z > fZmin)	phRhi = TMath::Abs(ph);	// Intersection inside chamber volume
+			Double_t ph = 2 * asin(C * sqrt((fRmax * fRmax - D * D) / (1.0 + 2.0 * C * D)));
+			Double_t z = z0 + ct * ph / (2.0 * C);
+			if (z < fZmax && z > fZmin)	phRhi = TMath::Abs(ph);	// Intersection inside chamber volume	
 		}
 		//  ... with left wall
@@ -314,7 +308,7 @@
 		if (Zdir > 0.0)
 		{
-			Double_t ph = 2.0*C*Zdir;
-			Double_t Rint = sqrt(D*D + (1.0 + 2.0*C*D)*pow(sin(ph / 2), 2) / (C*C));
-			if (Rint < fRmax && Rint > fRmin)	phZmn = TMath::Abs(ph);	// Intersection inside chamber volume
+			Double_t ph = 2.0 * C * Zdir;
+			Double_t Rint = sqrt(D * D + (1.0 + 2.0 * C * D) * pow(sin(ph / 2), 2) / (C * C));
+			if (Rint < fRmax && Rint > fRmin)	phZmn = TMath::Abs(ph);	// Intersection inside chamber volume	
 		}
 		//  ... with right wall
@@ -322,7 +316,7 @@
 		if (Zdir > 0.0)
 		{
-			Double_t ph = 2.0*C*Zdir;
-			Double_t Rint = sqrt(D*D + (1.0 + 2.0*C*D)*pow(sin(ph / 2), 2) / (C*C));
-			if (Rint < fRmax && Rint > fRmin)	phZmx = TMath::Abs(ph);	// Intersection inside chamber volume
+			Double_t ph = 2.0 * C * Zdir;
+			Double_t Rint = sqrt(D * D + (1.0 + 2.0 * C * D) * pow(sin(ph / 2), 2) / (C * C));
+			if (Rint < fRmax && Rint > fRmin)	phZmx = TMath::Abs(ph);	// Intersection inside chamber volume	
 		}
 		//
@@ -342,5 +336,5 @@
 		{
 			dPhase = ph_arr[iPos + 2] - ph_arr[iPos + 1];
-			tLength = dPhase*Scale;
+			tLength = dPhase * Scale;
 		}
 	}
@@ -349,5 +343,5 @@
 //
 // Return number of ionization clusters
-Bool_t TrkUtil::IonClusters(Double_t &Ncl, Double_t mass, TVectorD Par)
+Bool_t TrkUtil::IonClusters(Double_t& Ncl, Double_t mass, TVectorD Par)
 {
 	//
@@ -386,5 +380,5 @@
 			TVector3 p = ParToP(Par);
 			bg = p.Mag() / mass;
-			muClu = Nclusters(bg)*tLen;				// Avg. number of clusters
+			muClu = Nclusters(bg) * tLen;				// Avg. number of clusters
 
 			Ncl = gRandom->PoissonD(muClu);			// Actual number of clusters
@@ -392,5 +386,5 @@
 
 	}
-//
+	//
 	return Signal;
 }
@@ -413,39 +407,4 @@
 	//
 	//
-	/*
-	std::vector<double> bg{ 0.5, 0.8, 1., 2., 3., 4., 5., 8., 10.,
-	12., 15., 20., 50., 100., 200., 500., 1000. };
-	// He 90 - Isobutane 10
-	std::vector<double> ncl_He_Iso{ 42.94, 23.6,18.97,12.98,12.2,12.13,
-	12.24,12.73,13.03,13.29,13.63,14.08,15.56,16.43,16.8,16.95,16.98 };
-	//
-	// pure He
-	std::vector<double> ncl_He{ 11.79,6.5,5.23,3.59,3.38,3.37,3.4,3.54,3.63,
-	3.7,3.8,3.92,4.33,4.61,4.78,4.87,4.89 };
-	//
-	// Argon 50 - Ethane 50
-	std::vector<double> ncl_Ar_Eth{ 130.04,71.55,57.56,39.44,37.08,36.9,
-	37.25,38.76,39.68,40.49,41.53,42.91,46.8,48.09,48.59,48.85,48.93 };
-	//
-	// pure Argon
-	std::vector<double> ncl_Ar{ 88.69,48.93,39.41,27.09,25.51,25.43,25.69,
-	26.78,27.44,28.02,28.77,29.78,32.67,33.75,34.24,34.57,34.68 };
-	//
-	Int_t nPoints = (Int_t)bg.size();
-	bg.push_back(10000.);
-	std::vector<double> ncl;
-	switch (Opt)
-	{
-	case 0: ncl = ncl_He_Iso;			// He-Isobutane
-		break;
-	case 1: ncl = ncl_He;				// pure He
-		break;
-	case 2: ncl = ncl_Ar_Eth;			// Argon - Ethane
-		break;
-	case 3: ncl = ncl_Ar;				// pure Argon
-		break;
-	}
-	ncl.push_back(ncl[nPoints - 1]);
-	*/
 	const Int_t Npt = 18;
 	Double_t bg[Npt] = { 0.5, 0.8, 1., 2., 3., 4., 5., 8., 10.,
@@ -469,47 +428,25 @@
 	//
 	Double_t ncl[Npt];
-    	switch (Opt)
-    	{
-		case 0: std::copy(ncl_He_Iso, ncl_He_Iso + Npt, ncl);	// He-Isobutane
+	switch (Opt)
+	{
+	case 0: std::copy(ncl_He_Iso, ncl_He_Iso + Npt, ncl);	// He-Isobutane
 		break;
-		case 1: std::copy(ncl_He, ncl_He + Npt, ncl);		// pure He
+	case 1: std::copy(ncl_He, ncl_He + Npt, ncl);		// pure He
 		break;
-		case 2: std::copy(ncl_Ar_Eth, ncl_Ar_Eth + Npt, ncl);	// Argon - Ethane
+	case 2: std::copy(ncl_Ar_Eth, ncl_Ar_Eth + Npt, ncl);	// Argon - Ethane
 		break;
-		case 3: std::copy(ncl_Ar, ncl_Ar + Npt, ncl);		// pure Argon
+	case 3: std::copy(ncl_Ar, ncl_Ar + Npt, ncl);		// pure Argon
 		break;
-    	}
-	//
-	Int_t ilow = 0;
-	while (begam > bg[ilow])ilow++;
-	ilow--;
-	//std::cout << "ilow= " << ilow << ", low = " << bg[ilow] << ", val = " << begam
-	//	<< ", high = " << bg[ilow + 1] << std::endl;
-	//
-	Int_t ind[3] = { ilow, ilow + 1, ilow + 2 };
-	TVectorD y(3);
-	for (Int_t i = 0; i < 3; i++)y(i) = ncl[ind[i]];
-	TVectorD x(3);
-	for (Int_t i = 0; i < 3; i++)x(i) = bg[ind[i]];
-	TMatrixD Xval(3, 3);
-	for (Int_t i = 0; i < 3; i++)Xval(i, 0) = 1.0;
-	for (Int_t i = 0; i < 3; i++)Xval(i, 1) = x(i);
-	for (Int_t i = 0; i < 3; i++)Xval(i, 2) = x(i) * x(i);
-	//std::cout << "Xval:" << std::endl; Xval.Print();
-	Xval.Invert();
-	TVectorD coeff = Xval * y;
-	Double_t interp = coeff[0] + coeff[1] * begam + coeff[2] * begam * begam;
-	//std::cout << "val1= (" <<x(0)<<", "<< y(0) << "), val2= ("
-	//	<<x(1)<<", "<< y(1) << "), val3= ("
-	//	<<x(2)<<", "<< y(2)
-	//	<< "), result= (" <<begam<<", "<< interp<<")" << std::endl;
-	//
-	//if (TMath::IsNaN(interp))std::cout << "NaN found: bg= " << begam << ", Opt= " << Opt << std::endl;
-	if (begam < bg[0]) interp = 0.0;
-	//std::cout << "bg= " << begam << ", Opt= " << Opt <<", interp = "<<interp<< std::endl;
-	return 100*interp;
-}
-//
-Double_t TrkUtil::funcNcl(Double_t *xp, Double_t *par){
+	}
+	//
+	Double_t interp = 0.0;
+	TSpline3* sp3 = new TSpline3("sp3", bg, ncl, Npt);
+	if (begam > bg[0] && begam < bg[Npt - 1]) interp = sp3->Eval(begam);
+	if(begam < bg[0]) interp = bg[0];
+	if(begam > bg[Npt-1]) interp = bg[Npt-1];
+	return 100 * interp;
+}
+//
+Double_t TrkUtil::funcNcl(Double_t* xp, Double_t* par) {
 	Double_t bg = xp[0];
 	return Nclusters(bg);
Index: external/TrackCovariance/TrkUtil.h
===================================================================
--- external/TrackCovariance/TrkUtil.h	(revision a95da74cda6d94acac311904287f73ef752d400f)
+++ external/TrackCovariance/TrkUtil.h	(revision 53f474642cbcbb6dae79428731080cbdabba2e1d)
@@ -25,4 +25,5 @@
 	// Service routines
 	//
+	void SetB(Double_t Bz){ fBz = Bz; }
 	TVectorD XPtoPar(TVector3 x, TVector3 p, Double_t Q);
 	TVector3 ParToP(TVectorD Par);
@@ -70,5 +71,5 @@
 	// Cluster counting in gas
 	//
-	void SetBfield(Double_t Bz);
+	void SetBfield(Double_t Bz){ fBz = Bz; }
 	// Define gas volume (units = meters)
 	void SetDchBoundaries(Double_t Rmin, Double_t Rmax, Double_t Zmin, Double_t Zmax);
Index: external/TrackCovariance/VertexFit.cc
===================================================================
--- external/TrackCovariance/VertexFit.cc	(revision a95da74cda6d94acac311904287f73ef752d400f)
+++ external/TrackCovariance/VertexFit.cc	(revision 53f474642cbcbb6dae79428731080cbdabba2e1d)
@@ -18,4 +18,5 @@
 	fxCst.ResizeTo(3);
 	fCovCst.ResizeTo(3, 3);
+	fCovCstInv.ResizeTo(3, 3);
 	fXv.ResizeTo(3);
 	fcovXv.ResizeTo(3, 3);
@@ -31,4 +32,5 @@
 	fxCst.ResizeTo(3);
 	fCovCst.ResizeTo(3, 3);
+	fCovCstInv.ResizeTo(3, 3);
 	fXv.ResizeTo(3);
 	fcovXv.ResizeTo(3, 3);
@@ -38,6 +40,8 @@
 		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));
 	}
 	fChi2List.ResizeTo(fNtr);
@@ -54,4 +58,5 @@
 	fxCst.ResizeTo(3);
 	fCovCst.ResizeTo(3, 3);
+	fCovCstInv.ResizeTo(3, 3);
 	fXv.ResizeTo(3);
 	fcovXv.ResizeTo(3, 3);
@@ -61,5 +66,7 @@
 	{
 		fPar.push_back(new TVectorD(track[i]->GetObsPar()));
+		fParNew.push_back(new TVectorD(track[i]->GetObsPar()));
 		fCov.push_back(new TMatrixDSym(track[i]->GetCov()));
+		fCovNew.push_back(new TMatrixDSym(track[i]->GetCov()));
 	}
 }
@@ -73,7 +80,9 @@
 	{
 		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 (fWi[i])   { fWi[i]->Clear();			delete fWi[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]; }
 	}
@@ -81,4 +90,6 @@
 	fx0i.clear();
 	fai.clear();
+	fdi.clear();
+	fAti.clear();
 	fDi.clear();
 	fWi.clear();
@@ -88,5 +99,6 @@
 {
 	fxCst.Clear();		
-	fCovCst.Clear();		
+	fCovCst.Clear();
+	fCovCstInv.Clear();
 	fXv.Clear();		
 	fcovXv.Clear();		
@@ -95,9 +107,13 @@
 	for (Int_t i = 0; i < fNtr; i++)
 	{
-		fPar[i]->Clear();	delete fPar[i];
-		fCov[i]->Clear();	delete fCov[i];
+		fPar[i]->Clear();		delete fPar[i];
+		fParNew[i]->Clear();	delete fParNew[i];
+		fCov[i]->Clear();		delete fCov[i];
+		fCovNew[i]->Clear();	delete fCovNew[i];
 	}
 	fPar.clear();
-	fCov.clear();	
+	fParNew.clear();
+	fCov.clear();
+	fCovNew.clear();		
 	//
 	ResetWrkArrays();
@@ -424,5 +440,5 @@
 	//
 	// Get track parameters and their covariance
-	TVectorD par = *fPar[i];
+	TVectorD par = *fParNew[i];
 	TMatrixDSym Cov = *fCov[i];
 	//
@@ -433,6 +449,10 @@
 	//
 	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
+	fWinvi.push_back(new TMatrixDSym(Winv));	// Store W^-1 matrix
 	//
 	TMatrixDSym W = RegInv(Winv);				// W = (A*C*A')^-1
@@ -456,6 +476,5 @@
 {
 	//std::cout << "VertexFitter: just in" << std::endl;
-	if (fNtr < 2)
-	{
+	if (fNtr < 2 && !fVtxCst){
 		std::cout << "VertexFit::VertexFitter - Method called with less than 2 tracks - Aborting " << std::endl;
 		std::exit(1);
@@ -469,5 +488,10 @@
 	Double_t Chi2 = 0;
 	TVectorD x0 = fXv;	// If previous fit done
-	if (fRold < 0.0)for (Int_t i = 0; i < 3; i++)x0(i) = 1000.;	// Set to arbitrary large value if not
+	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
@@ -475,4 +499,5 @@
 	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;
 	//
 	// Iteration properties
@@ -529,7 +554,19 @@
 			H += Ds;
 			// update constant term
-			TVectorD xs = *fx0i[i];
+			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;
+			*/
 			cterm += Ds * xs;
 		}				// End loop on tracks
+		// Some additions in case of external constraints
+		if (fVtxCst) {
+			H += fCovCstInv;
+			cterm += fCovCstInv * fxCst;
+			DW1D += fCovCstInv;
+		}
 		//
 		// update vertex position
@@ -544,5 +581,5 @@
 		for (Int_t i = 0; i < fNtr; i++)
 		{
-			TVectorD lambda = (*fDi[i]) * (*fx0i[i] - x);
+			TVectorD lambda = (*fDi[i]) * (*fx0i[i] - x - *fdi[i]);
 			TMatrixDSym Wm1 = *fWinvi[i];
 			fChi2List(i) = Wm1.Similarity(lambda);
@@ -551,5 +588,9 @@
 			TVectorD b = (*fWi[i]) * (x - (*fx0i[i]));
 			for (Int_t j = 0; j < 3; j++)ffi[i] += a(j) * b(j) / fa2i[i];
+			TVectorD newPar = *fPar[i] - ((*fCov[i]) * (*fAti[i])) * lambda;
+			fParNew[i] = new TVectorD(newPar);
 		}
+		// Add external constraint to Chi2
+		if (fVtxCst) Chi2 += fCovCstInv.Similarity(x - fxCst);
 		//
 		TVectorD dx = x - x0;
@@ -562,7 +603,8 @@
 		// Store result
 		//
-		fXv = x;				// Vertex position
+		fXv = x;			// Vertex position
 		fcovXv = covX;		// Vertex covariance
 		fChi2 = Chi2;		// Vertex fit Chi2
+		//std::cout << "Found vertex " << fXv(0) << ", " << fXv(1) << ", " << fXv(2) << std::endl;
 	}		// end of iteration loop
 	//
@@ -603,5 +645,14 @@
 void VertexFit::AddVtxConstraint(TVectorD xv, TMatrixDSym cov)	// Add gaussian vertex constraint
 {
-	std::cout << "VertexFit::AddVtxConstraint: Not implemented yet" << std::endl;
+	//std::cout << "VertexFit::AddVtxConstraint: Not implemented yet" << std::endl;
+	fVtxCst = kTRUE;				// Vertex constraint flag
+	fxCst = xv;						// Constraint value
+	fCovCst = cov;					// Constraint covariance
+	fCovCstInv = cov;
+	fCovCstInv.Invert();				// Its inverse
+	//
+	// Set starting vertex as external constraint
+	fXv = fxCst;
+	fcovXv = fCovCst;
 }
 //
@@ -613,4 +664,6 @@
 	fPar.push_back(par);			// add new track
 	fCov.push_back(Cov);
+	fParNew.push_back(par);			// add new track
+	fCovNew.push_back(Cov);
 	//
 	// Reset previous vertex temp arrays
@@ -628,4 +681,6 @@
 	fPar.erase(fPar.begin() + iTrk);		// Remove track
 	fCov.erase(fCov.begin() + iTrk);
+	fParNew.erase(fParNew.begin() + iTrk);		// Remove track
+	fCovNew.erase(fCovNew.begin() + iTrk);
 	//
 	// Reset previous vertex temp arrays
Index: external/TrackCovariance/VertexFit.h
===================================================================
--- external/TrackCovariance/VertexFit.h	(revision a95da74cda6d94acac311904287f73ef752d400f)
+++ external/TrackCovariance/VertexFit.h	(revision 53f474642cbcbb6dae79428731080cbdabba2e1d)
@@ -23,9 +23,12 @@
 	Int_t fNtr;					// Number of tracks
 	std::vector<TVectorD*> fPar;		// Input parameter array
-	std::vector<TMatrixDSym*> fCov;// Input parameter covariances
+	std::vector<TVectorD*> fParNew;		// Updated parameter array
+	std::vector<TMatrixDSym*> fCov;		// Input parameter covariances
+	std::vector<TMatrixDSym*> fCovNew;	// Updated parameter covariances
 	// Constraints
 	Bool_t fVtxCst;				// Vertex constraint flag
 	TVectorD fxCst;				// Constraint value
-	TMatrixDSym fCovCst;			// Constraint covariance
+	TMatrixDSym fCovCst;			// Constraint 
+	TMatrixDSym fCovCstInv;			// Inverse of constraint covariance
 	//
 	// Results
@@ -39,7 +42,9 @@
 	// Work arrays
 	std::vector<Double_t> ffi;				// Fit phases
-	std::vector<TVectorD*> fx0i;				// Track expansion points
+	std::vector<TVectorD*> fx0i;			// Track expansion points
 	std::vector<TVectorD*> fai;				// dx/dphi
+	std::vector<TVectorD*> fdi;				// x-shift
 	std::vector<Double_t> fa2i;				// a'Wa
+	std::vector<TMatrixD*> fAti;			// A transposed
 	std::vector<TMatrixDSym*> fDi;			// W-WBW
 	std::vector<TMatrixDSym*> fWi;			// (ACA')^-1
