Index: external/TrackCovariance/ObsTrk.cc
===================================================================
--- external/TrackCovariance/ObsTrk.cc	(revision 1294bbaccb2cdf7fe3fe092f8fdaaa8ab86bad36)
+++ external/TrackCovariance/ObsTrk.cc	(revision f9aeea610d46ad5e5a508a3152713c4d53df921d)
@@ -1,4 +1,2 @@
-#include <iostream>
-
 #include <TMath.h>
 #include <TVector3.h>
@@ -7,135 +5,206 @@
 #include <TDecompChol.h>
 #include <TRandom.h>
-
+#include <iostream>
 #include "SolGridCov.h"
 #include "ObsTrk.h"
-
-using namespace std;
-
+//
+// Constructors
 // x(3) track origin, p(3) track momentum at origin, Q charge, B magnetic field in Tesla
 ObsTrk::ObsTrk(TVector3 x, TVector3 p, Double_t Q, Double_t B, SolGridCov *GC)
 {
-  fGC = GC;
-  fGenX = x;
-  fGenP = p;
-  fGenQ = Q;
-  fB = B;
-  fGenPar.ResizeTo(5);
-  fObsPar.ResizeTo(5);
-  fCov.ResizeTo(5, 5);
-  fGenPar = XPtoPar(x, p, Q);
-  fObsPar = GenToObsPar(fGenPar, fGC);
-  fObsX = ParToX(fObsPar);
-  fObsP = ParToP(fObsPar);
-  fObsQ = ParToQ(fObsPar);
-}
-
+	fGC = GC;
+	fGenX = x;
+	fGenP = p;
+	fGenQ = Q;
+	fB = B;
+	fGenPar.ResizeTo(5);
+	fGenParACTS.ResizeTo(6);
+	fObsPar.ResizeTo(5);
+	fObsParACTS.ResizeTo(6);
+	fCov.ResizeTo(5, 5);
+	fCovACTS.ResizeTo(6, 6);
+	fGenPar = XPtoPar(x,p,Q);
+	fGenParACTS = ParToACTS(fGenPar);
+	/*
+	std::cout << "ObsTrk::ObsTrk: fGenPar";
+	for (Int_t i = 0; i < 5; i++)std::cout << fGenPar(i) << ", ";
+	std::cout << std::endl;
+	*/
+	fObsPar = GenToObsPar(fGenPar, fGC);
+	fObsParACTS = ParToACTS(fObsPar);
+	fObsX = ParToX(fObsPar);
+	fObsP = ParToP(fObsPar);
+	fObsQ = ParToQ(fObsPar);
+	fCovACTS = CovToACTS(fCov);
+}
+//
+// Destructor
 ObsTrk::~ObsTrk()
 {
-}
-
+	fGenX.Clear();
+	fGenP.Clear();
+	fGenPar.Clear();
+	fGenParACTS.Clear();
+	fObsX.Clear();
+	fObsP.Clear();
+	fObsPar.Clear();
+	fObsParACTS.Clear();
+	fCov.Clear();
+	fCovACTS.Clear();
+}
 TVectorD ObsTrk::XPtoPar(TVector3 x, TVector3 p, Double_t Q)
 {
-  TVectorD Par(5);
-  // Transverse parameters
-  Double_t a = -Q * fB * 0.2998; // Units are Tesla, GeV and m
-  Double_t pt = p.Pt();
-  Double_t C = a / (2 * pt); // Half curvature
-
-  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 D; // Impact parameter D
-  if (pt < 10.0) D = (T - pt) / a;
-  else D = (-2 * cross + a * r2) / (T + pt);
-
-  Par(0) = D; // Store D
-  Par(1) = phi0; // Store phi0
-  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 ct = p(2) / pt;
-  Double_t z0 = x(2) - ct * st;
-
-  Par(3) = z0; // Store z0
-  Par(4) = ct; // Store cot(theta)
-
-  return Par;
-}
-
+	//
+	TVectorD Par(5);
+	// Transverse parameters
+	Double_t a = -Q*fB*0.2998;			// Units are Tesla, GeV and meters
+	Double_t pt = p.Pt();
+	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 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 D;							// Impact parameter D
+	if (pt < 10.0) D = (T - pt) / a;
+	else D = (-2 * cross + a*r2) / (T + pt);
+	//
+	Par(0) = D;		// Store D
+	Par(1) = phi0;	// Store phi0
+	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 ct = p(2) / pt;
+	Double_t z0 = x(2) - ct*st;
+	//
+	Par(3) = z0;		// Store z0
+	Par(4) = ct;		// Store cot(theta)
+	//
+	return Par;
+}
+//
 TVector3 ObsTrk::ParToX(TVectorD Par)
 {
-  Double_t D    = Par(0);
-  Double_t phi0 = Par(1);
-  Double_t z0   = Par(3);
-
-  TVector3 Xval;
-  Xval(0) = -D * TMath::Sin(phi0);
-  Xval(1) =  D * TMath::Cos(phi0);
-  Xval(2) =  z0;
-
-  return Xval;
-}
-
+	Double_t D    = Par(0);
+	Double_t phi0 = Par(1);
+	Double_t z0   = Par(3);
+	//
+	TVector3 Xval;
+	Xval(0) = -D*TMath::Sin(phi0); 
+	Xval(1) =  D*TMath::Cos(phi0);  
+	Xval(2) =  z0;
+	//
+	return Xval;
+}
+//
 TVector3 ObsTrk::ParToP(TVectorD Par)
 {
-  Double_t C    = Par(2);
-  Double_t phi0 = Par(1);
-  Double_t ct   = Par(4);
-  //
-  TVector3 Pval;
-  Double_t pt = fB * 0.2998 / TMath::Abs(2 * C);
-  Pval(0) = pt * TMath::Cos(phi0);
-  Pval(1) = pt * TMath::Sin(phi0);
-  Pval(2) = pt * ct;
-  //
-  return Pval;
-}
+	Double_t C    = Par(2);
+	Double_t phi0 = Par(1);
+	Double_t ct   = Par(4);
+	//
+	TVector3 Pval;
+	Double_t pt = fB*0.2998 / TMath::Abs(2 * C);
+	Pval(0) = pt*TMath::Cos(phi0);
+	Pval(1) = pt*TMath::Sin(phi0);
+	Pval(2) = pt*ct;
+	//
+	return Pval;
+}
+//
 
 Double_t ObsTrk::ParToQ(TVectorD Par)
 {
-  return TMath::Sign(1.0, -Par(2));
-}
-
+	return TMath::Sign(1.0, -Par(2));
+}
+//
 TVectorD ObsTrk::GenToObsPar(TVectorD gPar, SolGridCov *GC)
 {
-  TVector3 p = ParToP(gPar);
-  Double_t pt = p.Pt();
-  Double_t tanTh = 1.0 / TMath::Abs(gPar(4));
-  Double_t angd = TMath::ATan(tanTh) * 180. / TMath::Pi();
-  // Check ranges
-  Double_t minPt = GC->GetMinPt ();
-  if (pt < minPt) cout << "Warning ObsTrk::GenToObsPar: pt " << pt << " is below grid range of " << minPt << endl;
-  Double_t maxPt = GC->GetMaxPt();
-  if (pt > maxPt) cout << "Warning ObsTrk::GenToObsPar: pt " << pt << " is above grid range of " << maxPt << endl;
-  Double_t minAn = GC->GetMinAng();
-  if (angd < minAn) cout << "Warning ObsTrk::GenToObsPar: angle " << angd
-    << " is below grid range of " << minAn << endl;
-  Double_t maxAn = GC->GetMaxAng();
-  if (angd > maxAn) cout << "Warning ObsTrk::GenToObsPar: angle " << angd
-    << " is above grid range of " << maxAn << endl;
-  TMatrixDSym Cov = GC->GetCov(pt, angd);
-  fCov = Cov;
-  // Now do Choleski decomposition and random number extraction, with appropriate stabilization
-  TMatrixDSym CvN = Cov;
-  TMatrixDSym DCv(5); DCv.Zero();
-  TMatrixDSym DCvInv(5); DCvInv.Zero();
-  for (Int_t id = 0; id < 5; id++)
-  {
-    Double_t dVal = TMath::Sqrt(Cov(id, id));
-    DCv   (id, id) = dVal;
-    DCvInv(id, id) = 1.0 / dVal;
-  }
-  CvN.Similarity(DCvInv); // Normalize diagonal to 1
-  TDecompChol Chl(CvN);
-  Bool_t OK = Chl.Decompose(); // Choleski decomposition of normalized matrix
-  TMatrixD U = Chl.GetU(); // Get Upper triangular matrix
-  TMatrixD Ut(TMatrixD::kTransposed, U); // Transposed of U (lower triangular)
-  TVectorD r(5);
-  for (Int_t i = 0; i < 5; i++) r(i) = gRandom->Gaus(0.0, 1.0); // Array of normal random numbers
-  TVectorD oPar = gPar + DCv * (Ut * r); // Observed parameter vector
-
-  return oPar;
-}
+	TVector3 p = ParToP(gPar);
+	Double_t pt = p.Pt();
+	Double_t tanTh = 1.0 / TMath::Abs(gPar(4));
+	Double_t angd = TMath::ATan(tanTh)*180. / TMath::Pi();
+	//
+	// Check ranges
+	Double_t minPt = GC->GetMinPt ();
+	if (pt < minPt) std::cout << "Warning ObsTrk::GenToObsPar: pt " << pt << " is below grid range of " << minPt << std::endl;
+	Double_t maxPt = GC->GetMaxPt();
+	if (pt > maxPt) std::cout << "Warning ObsTrk::GenToObsPar: pt " << pt << " is above grid range of " << maxPt << std::endl;
+	Double_t minAn = GC->GetMinAng();
+	if (angd < minAn) std::cout << "Warning ObsTrk::GenToObsPar: angle " << angd 
+		<< " is below grid range of " << minAn << std::endl;
+	Double_t maxAn = GC->GetMaxAng();
+	if (angd > maxAn) std::cout << "Warning ObsTrk::GenToObsPar: angle " << angd
+		<< " is above grid range of " << maxAn << std::endl;
+	//
+	TMatrixDSym Cov = GC->GetCov(pt, angd);
+	fCov = Cov;
+	//
+	// Now do Choleski decomposition and random number extraction, with appropriate stabilization
+	//
+	TMatrixDSym CvN = Cov;
+	TMatrixDSym DCv(5); DCv.Zero();
+	TMatrixDSym DCvInv(5); DCvInv.Zero();
+	for (Int_t id = 0; id < 5; id++)
+	{
+		Double_t dVal = TMath::Sqrt(Cov(id, id));
+		DCv   (id, id) = dVal;
+		DCvInv(id, id) = 1.0 / dVal;
+	}
+	CvN.Similarity(DCvInv);			// Normalize diagonal to 1
+	TDecompChol Chl(CvN);
+	Bool_t OK = Chl.Decompose();		// Choleski decomposition of normalized matrix
+	TMatrixD U = Chl.GetU();			// Get Upper triangular matrix
+	TMatrixD Ut(TMatrixD::kTransposed, U); // Transposed of U (lower triangular)
+	TVectorD r(5);
+	for (Int_t i = 0; i < 5; i++)r(i) = gRandom->Gaus(0.0, 1.0);		// Array of normal random numbers
+	TVectorD oPar = gPar + DCv*(Ut*r);	// Observed parameter vector
+	//
+	return oPar;
+}
+// Parameter conversion to ACTS format
+TVectorD ObsTrk::ParToACTS(TVectorD Par)
+{
+	TVectorD pACTS(6);	// Return vector
+	//
+	Double_t b = -0.29988*fB / 2.;
+	pACTS(0) = 1000*Par(0);		// D from m to mm
+	pACTS(1) = 1000 * Par(3);	// z0 from m to mm
+	pACTS(2) = Par(1);			// Phi0 is unchanged
+	pACTS(3) = TMath::ATan(1.0 / Par(4)) + TMath::PiOver2();		// Theta in [0, pi] range
+	pACTS(4) = Par(2) / (b*TMath::Sqrt(1 + Par(4)*Par(4)));		// q/p in GeV
+	pACTS(5) = 0.0;				// Time: currently undefined
+	//
+	return pACTS;
+}
+// Covariance conversion to ACTS format
+TMatrixDSym ObsTrk::CovToACTS(TMatrixDSym Cov)
+{
+	TMatrixDSym cACTS(6); cACTS.Zero();
+	Double_t b = -0.29988*fB / 2.;
+	//
+	// Fill derivative matrix
+	TMatrixD A(5, 5);	A.Zero();
+	Double_t ct = fGenPar(4);	// cot(theta)
+	Double_t C = fGenPar(2);		// half curvature
+	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(3, 1) = 1000.;		// z0-z0 conversion to mm
+	A(4, 3) = -1.0 / (1.0 + ct*ct); // theta - cot(theta)
+	A(4, 4) = -C*ct / (b*pow(1.0 + ct*ct,3.0/2.0)); // q/p-cot(theta)
+	//
+	TMatrixDSym Cv = Cov;
+	TMatrixD At(5, 5);
+	At.Transpose(A);
+	Cv.Similarity(At);
+	TMatrixDSub(cACTS, 0, 4, 0, 4) = Cv;
+	cACTS(5, 5) = 0.1;	// Currently undefined: set to arbitrary value to avoid crashes
+	//
+	return cACTS;
+}
+
+
+
+	
Index: external/TrackCovariance/ObsTrk.h
===================================================================
--- external/TrackCovariance/ObsTrk.h	(revision 1294bbaccb2cdf7fe3fe092f8fdaaa8ab86bad36)
+++ external/TrackCovariance/ObsTrk.h	(revision f9aeea610d46ad5e5a508a3152713c4d53df921d)
@@ -1,56 +1,82 @@
+//
 #ifndef G__OBSTRK_H
 #define G__OBSTRK_H
-
 #include <TVector3.h>
 #include <TVectorD.h>
 #include <TMatrixDSym.h>
-
-class SolGridCov;
-
+#include <TDecompChol.h>
+#include "SolGridCov.h"
+//
 // Class to handle smearing of generated charged particle tracks
-
+//
 // Author: F. Bedeschi
 //         INFN - Sezione di Pisa, Italy
 //
 class ObsTrk{
-  // Class to handle simulation of tracking resolution
-  // Prefix Obs marks variables after resolution smearing
-  // Prefix Gen marks variables before resolution smearing
-private:
-  Double_t fB;      // Solenoid magnetic field
-  SolGridCov *fGC;  // Covariance matrix grid
-  Double_t fGenQ;   // Generated track charge
-  Double_t fObsQ;   // Observed  track charge
-  TVector3 fGenX;   // Generated track origin (x,y,z)
-  TVector3 fObsX;   // Observed  track origin (x,y,z) @ track min. approach
-  TVector3 fGenP;   // Generated track momentum at track origin
-  TVector3 fObsP;   // Observed  track momentum @ track minimum approach
-  TVectorD fGenPar; // Generated helix track parameters (D, phi0, C, z0, cot(th))
-  TVectorD fObsPar; // Observed  helix track parameters (D, phi0, C, z0, cot(th))
-  TMatrixDSym fCov; // INterpolated covariance of track parameters
+	//
+	// Class to handle simulation of tracking resolution
+	// Prefix Obs marks variables after resolution smearing
+	// Prefix Gen marks variables before resolution smearing
+	//
+private:	
+	Double_t fB;						// Solenoid magnetic field
+	SolGridCov *fGC;					// Covariance matrix grid
+	Double_t fGenQ;					// Generated track charge
+	Double_t fObsQ;					// Observed  track charge
+	TVector3 fGenX;					// Generated track origin (x,y,z)
+	TVector3 fObsX;					// Observed  track origin (x,y,z) @ track min. approach 
+	TVector3 fGenP;					// Generated track momentum at track origin
+	TVector3 fObsP;					// Observed  track momentum @ track minimum approach
+	TVectorD fGenPar;				// Generated helix track parameters (D, phi0, C, z0, cot(th))
+	TVectorD fGenParACTS;			// Generated helix track parameters (D, z0, phi0, th, q/p, time)
+	TVectorD fObsPar;				// Observed  helix track parameters (D, phi0, C, z0, cot(th))
+	TVectorD fObsParACTS;			// Observed  helix track parameters (D, z0, phi0, th, q/p, time)
+	TMatrixDSym fCov;				// INterpolated covariance of track parameters
+	TMatrixDSym fCovACTS;			// Covariance of track parameters in ACTS format
+									// (D, z0, phi0, theta, q/p, time)
+	//
+	// Conversion to ACTS parametrization
+	//
+	TVectorD ParToACTS(TVectorD Par);		// Parameter conversion
+	TMatrixDSym CovToACTS(TMatrixDSym Cov);	// Covariance conversion
+	//
 public:
-  // x(3) track origin, p(3) track momentum at origin, Q charge, B magnetic field in Tesla
-  ObsTrk(TVector3 x, TVector3 p, Double_t Q, Double_t B, SolGridCov *GC); // Initialize and generate smeared track
-  ~ObsTrk();
-  // Service routines
-  TVectorD XPtoPar(TVector3 x, TVector3 p, Double_t Q);
-  TVectorD GenToObsPar(TVectorD gPar, SolGridCov *GC);
-  TVector3 ParToX(TVectorD Par);
-  TVector3 ParToP(TVectorD Par);
-  Double_t ParToQ(TVectorD Par);
-  // Accessors
-  // Generator level X, P, Q
-  Double_t GetGenQ() { return fGenQ; }
-  TVector3 GetGenX() { return fGenX; }
-  TVector3 GetGenP() { return fGenP; }
-  // D, phi0, C, z0, cot(th)
-  TVectorD GetGenPar() { return fGenPar; }
-  // Observed level X, P, Q
-  Double_t GetObsQ() { return fObsQ; }
-  TVector3 GetObsX() { return fObsX; }
-  TVector3 GetObsP() { return fObsP; }
-  // D, phi0, C, z0, cot(th)
-  TVectorD GetObsPar() { return fObsPar; }
-  TMatrixDSym GetCov() { return fCov; }
+	//
+	// Constructors
+	// x(3) track origin, p(3) track momentum at origin, Q charge, B magnetic field in Tesla
+	ObsTrk(TVector3 x, TVector3 p, Double_t Q, Double_t B, SolGridCov *GC);	// Initialize and generate smeared track
+	// Destructor
+	~ObsTrk();
+	//
+	// Service routines
+	//
+	TVectorD XPtoPar(TVector3 x, TVector3 p, Double_t Q);
+	TVectorD GenToObsPar(TVectorD gPar, SolGridCov *GC);
+	TVector3 ParToX(TVectorD Par);
+	TVector3 ParToP(TVectorD Par);
+	Double_t ParToQ(TVectorD Par);
+	//
+	// Accessors
+	//
+	// Generator level:
+	// X, P, Q
+	Double_t GetGenQ()	{ return fGenQ; }
+	TVector3 GetGenX()	{ return fGenX; }
+	TVector3 GetGenP()	{ return fGenP; }
+	// D, phi0, C, z0, cot(th)
+	TVectorD GetGenPar()	{ return fGenPar; }
+	// D, z0, phi0, theta, q/p, time
+	TVectorD GetGenParACTS()	{ return fGenParACTS; }
+	// Observed level X, P, Q
+	Double_t GetObsQ()	{ return fObsQ; }
+	TVector3 GetObsX()	{ return fObsX; }
+	TVector3 GetObsP()	{ return fObsP; }
+	// D, phi0, C, z0, cot(th)
+	TVectorD GetObsPar()	{ return fObsPar; }
+	// D, z0, phi0, theta, q/p, time
+	TVectorD GetObsParACTS()	{ return fObsParACTS; }
+	// Covariances
+	TMatrixDSym GetCov(){ return fCov; }
+	TMatrixDSym GetCovACTS(){ return fCovACTS; };
 };
 
