Index: external/TrackCovariance/TrkUtil.cc
===================================================================
--- external/TrackCovariance/TrkUtil.cc	(revision 82db14533dc34d87c1c2389909a6bb0f351b6854)
+++ external/TrackCovariance/TrkUtil.cc	(revision c5696dd9091f908085bb41c6daa33633af86b66c)
@@ -1,4 +1,5 @@
 #include "TrkUtil.h"
 #include <iostream>
+#include <algorithm>
 
 // Constructor
@@ -46,6 +47,6 @@
 	Double_t r2 = x.Perp2();
 	Double_t cross = x(0) * p(1) - x(1) * p(0);
-	Double_t T = TMath::Sqrt(pt * pt - 2 * a * cross + a * a * r2);
-	Double_t phi0 = TMath::ATan2((p(1) - a * x(0)) / T, (p(0) + a * x(1)) / T);	// Phi0
+	Double_t T = sqrt(pt * pt - 2 * a * cross + a * a * r2);
+	Double_t phi0 = atan2((p(1) - a * x(0)) / T, (p(0) + a * x(1)) / T);	// Phi0
 	Double_t D;							// Impact parameter D
 	if (pt < 10.0) D = (T - pt) / a;
@@ -56,6 +57,6 @@
 	Par(2) = C;		// Store C
 	//Longitudinal parameters
-	Double_t B = C * TMath::Sqrt(TMath::Max(r2 - D * D, 0.0) / (1 + 2 * C * D));
-	Double_t st = TMath::ASin(B) / C;
+	Double_t B = C * sqrt(TMath::Max(r2 - D * D, 0.0) / (1 + 2 * C * D));
+	Double_t st = asin(B) / C;
 	Double_t ct = p(2) / pt;
 	Double_t z0 = x(2) - ct * st;
@@ -84,6 +85,6 @@
 	//
 	TVector3 Xval;
-	Xval(0) = -D * TMath::Sin(phi0);
-	Xval(1) = D * TMath::Cos(phi0);
+	Xval(0) = -D * sin(phi0);
+	Xval(1) = D * cos(phi0);
 	Xval(2) = z0;
 	//
@@ -106,6 +107,6 @@
 	TVector3 Pval;
 	Double_t pt = Bz * cSpeed() / TMath::Abs(2 * C);
-	Pval(0) = pt * TMath::Cos(phi0);
-	Pval(1) = pt * TMath::Sin(phi0);
+	Pval(0) = pt * cos(phi0);
+	Pval(1) = pt * sin(phi0);
 	Pval(2) = pt * ct;
 	//
@@ -128,6 +129,6 @@
 	pACTS(1) = 1000 * Par(3);	// z0 from m to mm
 	pACTS(2) = Par(1);			// Phi0 is unchanged
-	pACTS(3) = TMath::ATan2(1.0, Par(4));		// Theta in [0, pi] range
-	pACTS(4) = Par(2) / (b * TMath::Sqrt(1 + Par(4) * Par(4)));		// q/p in GeV
+	pACTS(3) = atan2(1.0, Par(4));		// Theta in [0, pi] range
+	pACTS(4) = Par(2) / (b * sqrt(1 + Par(4) * Par(4)));		// q/p in GeV
 	pACTS(5) = 0.0;				// Time: currently undefined
 	//
@@ -146,5 +147,5 @@
 	A(0, 0) = 1000.;		// D-D	conversion to mm
 	A(1, 2) = 1.0;		// phi0-phi0
-	A(2, 4) = 1.0 / (TMath::Sqrt(1.0 + ct * ct) * b);	// q/p-C
+	A(2, 4) = 1.0 / (sqrt(1.0 + ct * ct) * b);	// q/p-C
 	A(3, 1) = 1000.;		// z0-z0 conversion to mm
 	A(4, 3) = -1.0 / (1.0 + ct * ct); // theta - cot(theta)
@@ -267,5 +268,5 @@
 		//
 		// Track length per unit phase change 
-		Double_t Scale = TMath::Sqrt(1.0 + ct*ct) / (2.0*TMath::Abs(C));
+		Double_t Scale = sqrt(1.0 + ct*ct) / (2.0*TMath::Abs(C));
 		//
 		// Find intersections with chamber boundaries
@@ -281,5 +282,5 @@
 		if (Rtop > fRmin && TMath::Abs(D) < fRmin) // *** don't treat large D tracks for the moment ***
 		{
-			Double_t ph = 2 * TMath::ASin(C*TMath::Sqrt((fRmin*fRmin - D*D) / (1.0 + 2.0*C*D)));
+			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);
 
@@ -289,5 +290,5 @@
 			//
 			// Include second branch of loopers
-			Double_t ph2 = TMath::TwoPi() - TMath::Abs(ph);
+			Double_t ph2 = 2*TMath::Pi() - TMath::Abs(ph);
 			if (ph < 0)ph2 = -ph2;
 			z = z0 + ct * ph2 / (2.0 * C);
@@ -297,5 +298,5 @@
 		if (Rtop > fRmax && TMath::Abs(D) < fRmax) // *** don't treat large D tracks for the moment ***
 		{
-			Double_t ph = 2 * TMath::ASin(C*TMath::Sqrt((fRmax*fRmax - D*D) / (1.0 + 2.0*C*D)));
+			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	
@@ -306,5 +307,5 @@
 		{
 			Double_t ph = 2.0*C*Zdir;
-			Double_t Rint = TMath::Sqrt(D*D + (1.0 + 2.0*C*D)*pow(TMath::Sin(ph / 2), 2) / (C*C));
+			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	
 		}
@@ -314,5 +315,5 @@
 		{
 			Double_t ph = 2.0*C*Zdir;
-			Double_t Rint = TMath::Sqrt(D*D + (1.0 + 2.0*C*D)*pow(TMath::Sin(ph / 2), 2) / (C*C));
+			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	
 		}
@@ -324,5 +325,7 @@
 		Double_t ph_arr[Nint] = { phRin, phRin2, phRhi, phZmn, phZmx };
 		Int_t srtind[Nint];
-		TMath::Sort(Nint, ph_arr, srtind, kFALSE);
+		//TMath::Sort(Nint, ph_arr, srtind, kFALSE);
+		for (Int_t i = 0; i < Nint; i++) { srtind[i] = i; }
+		std::sort(srtind, srtind + Nint, CompareAsc<const Double_t*>(ph_arr));
 		Int_t iPos = -1;		// First element > 0
 		for (Int_t i = 0; i < Nint; i++)
@@ -462,5 +465,5 @@
 	//	<< "), result= (" <<begam<<", "<< interp<<")" << std::endl;
 	//
-	if (TMath::IsNaN(interp))std::cout << "NaN found: bg= " << begam << ", Opt= " << Opt << 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;
