Index: external/TrackCovariance/TrkUtil.cc
===================================================================
--- external/TrackCovariance/TrkUtil.cc	(revision a95da74cda6d94acac311904287f73ef752d400f)
+++ external/TrackCovariance/TrkUtil.cc	(revision 65776c0b380ce077906fc756046dfe0b553a9da2)
@@ -2,4 +2,5 @@
 #include <iostream>
 #include <algorithm>
+#include <TSpline.h>
 
 // Constructor
@@ -45,5 +46,5 @@
 	Double_t C = a / (2 * pt);			// Half curvature
 	//std::cout << "ObsTrk::XPtoPar: fB = " << fB << ", a = " << a << ", pt = " << pt << ", C = " << C << std::endl;
-	Double_t r2 = x.Perp2();
+	Double_t r2 = x(0) * x(0) + x(1) * x(1);
 	Double_t cross = x(0) * p(1) - x(1) * p(0);
 	Double_t T = sqrt(pt * pt - 2 * a * cross + a * a * r2);
@@ -60,5 +61,8 @@
 	Double_t st = asin(B) / C;
 	Double_t ct = p(2) / pt;
-	Double_t z0 = x(2) - ct * st;
+	Double_t z0;
+	Double_t dot = x(0) * p(0) + x(1) * p(1);
+	if (dot > 0.0) z0 = x(2) - ct * st;
+	else z0 = x(2) + ct * st;
 	//
 	Par(3) = z0;		// Store z0
@@ -96,5 +100,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 +237,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 +271,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 +303,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 +311,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 +319,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 +339,5 @@
 		{
 			dPhase = ph_arr[iPos + 2] - ph_arr[iPos + 1];
-			tLength = dPhase*Scale;
+			tLength = dPhase * Scale;
 		}
 	}
@@ -349,5 +346,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 +383,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 +389,5 @@
 
 	}
-//
+	//
 	return Signal;
 }
@@ -413,39 +410,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 +431,23 @@
 	//
 	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);
+	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 65776c0b380ce077906fc756046dfe0b553a9da2)
@@ -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);
@@ -53,5 +54,5 @@
 		Double_t c = 2.99792458e8;	// speed of light m/sec
 		//return TMath::C()*1.0e-9;	// Incompatible with root5
-		return c*1.0e-9; 		// Reduced speed of light
+		return c*1.0e-9; 		// Reduced speed of light	
 	}
 	//
@@ -70,6 +71,6 @@
 	// Cluster counting in gas
 	//
-	void SetBfield(Double_t Bz);
-	// Define gas volume (units = meters)
+	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);
 	// Gas mixture selection
