Index: Makefile
===================================================================
--- Makefile	(revision 3051358a56b3d611ea993535575ffebc4f158441)
+++ Makefile	(revision ff9fb2d9367d146a8325d4b77b9a9a4c0bd2ada8)
@@ -408,4 +408,5 @@
 	modules/MomentumSmearing.h \
 	modules/TrackSmearing.h \
+	modules/TrackCovariance.h \
 	modules/ImpactParameterSmearing.h \
 	modules/TimeSmearing.h \
@@ -645,4 +646,12 @@
 tmp/external/Hector/H_VerticalQuadrupole.$(ObjSuf): \
 	external/Hector/H_VerticalQuadrupole.$(SrcSuf)
+tmp/external/TrackCovariance/ObsTrk.$(ObjSuf): \
+	external/TrackCovariance/ObsTrk.$(SrcSuf)
+tmp/external/TrackCovariance/SolGeom.$(ObjSuf): \
+	external/TrackCovariance/SolGeom.$(SrcSuf)
+tmp/external/TrackCovariance/SolGridCov.$(ObjSuf): \
+	external/TrackCovariance/SolGridCov.$(SrcSuf)
+tmp/external/TrackCovariance/SolTrack.$(ObjSuf): \
+	external/TrackCovariance/SolTrack.$(SrcSuf)
 tmp/modules/AngularSmearing.$(ObjSuf): \
 	modules/AngularSmearing.$(SrcSuf) \
@@ -1004,4 +1013,9 @@
 	external/ExRootAnalysis/ExRootFilter.h \
 	external/ExRootAnalysis/ExRootResult.h
+tmp/modules/TrackCovariance.$(ObjSuf): \
+	modules/TrackCovariance.$(SrcSuf) \
+	modules/TrackCovariance.h \
+	classes/DelphesClasses.h \
+	external/TrackCovariance/SolGeom.h
 tmp/modules/TrackPileUpSubtractor.$(ObjSuf): \
 	modules/TrackPileUpSubtractor.$(SrcSuf) \
@@ -1131,4 +1145,8 @@
 	tmp/external/Hector/H_VerticalKicker.$(ObjSuf) \
 	tmp/external/Hector/H_VerticalQuadrupole.$(ObjSuf) \
+	tmp/external/TrackCovariance/ObsTrk.$(ObjSuf) \
+	tmp/external/TrackCovariance/SolGeom.$(ObjSuf) \
+	tmp/external/TrackCovariance/SolGridCov.$(ObjSuf) \
+	tmp/external/TrackCovariance/SolTrack.$(ObjSuf) \
 	tmp/modules/AngularSmearing.$(ObjSuf) \
 	tmp/modules/BTagging.$(ObjSuf) \
@@ -1170,4 +1188,5 @@
 	tmp/modules/TrackCountingBTagging.$(ObjSuf) \
 	tmp/modules/TrackCountingTauTagging.$(ObjSuf) \
+	tmp/modules/TrackCovariance.$(ObjSuf) \
 	tmp/modules/TrackPileUpSubtractor.$(ObjSuf) \
 	tmp/modules/TrackSmearing.$(ObjSuf) \
@@ -1884,4 +1903,8 @@
 	@touch $@
 
+modules/TrackCovariance.h: \
+	classes/DelphesModule.h
+	@touch $@
+
 modules/ExampleModule.h: \
 	classes/DelphesModule.h
Index: doc/genMakefile.tcl
===================================================================
--- doc/genMakefile.tcl	(revision 3051358a56b3d611ea993535575ffebc4f158441)
+++ doc/genMakefile.tcl	(revision ff9fb2d9367d146a8325d4b77b9a9a4c0bd2ada8)
@@ -321,5 +321,5 @@
 dictDeps {DISPLAY_DICT} {display/DisplayLinkDef.h}
 
-sourceDeps {DELPHES} {classes/*.cc} {modules/*.cc} {external/ExRootAnalysis/*.cc} {external/Hector/*.cc}
+sourceDeps {DELPHES} {classes/*.cc} {modules/*.cc} {external/ExRootAnalysis/*.cc} {external/Hector/*.cc} {external/TrackCovariance/*.cc}
 
 sourceDeps {FASTJET} {modules/FastJet*.cc} {modules/RunPUPPI.cc} {external/PUPPI/*.cc} {external/fastjet/*.cc} {external/fastjet/tools/*.cc} {external/fastjet/plugins/*/*.cc} {external/fastjet/contribs/*/*.cc}
Index: external/TrackCovariance/AUTHORS
===================================================================
--- external/TrackCovariance/AUTHORS	(revision ff9fb2d9367d146a8325d4b77b9a9a4c0bd2ada8)
+++ external/TrackCovariance/AUTHORS	(revision ff9fb2d9367d146a8325d4b77b9a9a4c0bd2ada8)
@@ -0,0 +1,1 @@
+Franco Bedeschi
Index: external/TrackCovariance/ObsTrk.cc
===================================================================
--- external/TrackCovariance/ObsTrk.cc	(revision ff9fb2d9367d146a8325d4b77b9a9a4c0bd2ada8)
+++ external/TrackCovariance/ObsTrk.cc	(revision ff9fb2d9367d146a8325d4b77b9a9a4c0bd2ada8)
@@ -0,0 +1,141 @@
+#include <iostream>
+
+#include <TMath.h>
+#include <TVector3.h>
+#include <TMatrixD.h>
+#include <TMatrixDSym.h>
+#include <TDecompChol.h>
+#include <TRandom.h>
+
+#include "SolGridCov.h"
+#include "ObsTrk.h"
+
+using namespace std;
+
+// 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);
+}
+
+ObsTrk::~ObsTrk()
+{
+}
+
+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;
+}
+
+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;
+}
+
+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 ObsTrk::ParToQ(TVectorD Par)
+{
+  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;
+}
Index: external/TrackCovariance/ObsTrk.h
===================================================================
--- external/TrackCovariance/ObsTrk.h	(revision ff9fb2d9367d146a8325d4b77b9a9a4c0bd2ada8)
+++ external/TrackCovariance/ObsTrk.h	(revision ff9fb2d9367d146a8325d4b77b9a9a4c0bd2ada8)
@@ -0,0 +1,57 @@
+#ifndef G__OBSTRK_H
+#define G__OBSTRK_H
+
+#include <TVector3.h>
+#include <TVectorD.h>
+#include <TMatrixDSym.h>
+
+class SolGridCov;
+
+// 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
+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; }
+};
+
+#endif
Index: external/TrackCovariance/SolGeom.cc
===================================================================
--- external/TrackCovariance/SolGeom.cc	(revision ff9fb2d9367d146a8325d4b77b9a9a4c0bd2ada8)
+++ external/TrackCovariance/SolGeom.cc	(revision ff9fb2d9367d146a8325d4b77b9a9a4c0bd2ada8)
@@ -0,0 +1,110 @@
+﻿#include <iostream>
+
+#include <TString.h>
+
+#include "SolGeom.h"
+#include "SolTrack.h"
+
+using namespace std;
+
+SolGeom::SolGeom()
+{
+  //
+  // Magnetic field
+  //
+  fB = 2.0;
+  //
+  // Create arrays
+  //
+  ftyLay = new Int_t[fNlMax];     // Layer type 1 = R (barrel) or 2 = z (forward/backward)
+  fLyLabl = new TString[fNlMax];  // Layer label
+  fxMin = new Double_t[fNlMax];   // Minimum dimension z for barrel  or R for forward
+  fxMax = new Double_t[fNlMax];   // Maximum dimension z for barrel  or R for forward
+  frPos = new Double_t[fNlMax];   // R/z location of layer
+  fthLay = new Double_t[fNlMax];  // Thickness (meters)
+  frlLay = new Double_t[fNlMax];  // Radiation length (meters)
+  fnmLay = new Int_t[fNlMax];     // Number of measurements in layers (1D or 2D)
+  fstLayU = new Double_t[fNlMax]; // Stereo angle (rad) - 0(pi/2) = axial(z) layer - Upper side
+  fstLayL = new Double_t[fNlMax]; // Stereo angle (rad) - 0(pi/2) = axial(z) layer - Lower side
+  fsgLayU = new Double_t[fNlMax]; // Resolution Upper side (meters) - 0 = no measurement
+  fsgLayL = new Double_t[fNlMax]; // Resolution Lower side (meters) - 0 = no measurement
+  fflLay = new Bool_t[fNlMax];    // measurement flag = T, scattering only = F
+  //
+  // Load geometry info in SolGeom.h
+  //
+  fNlay = 0; // Actual number of layers
+  fBlay = 0; // Nr. of barrel layers
+  fFlay = 0; // Nr. of forward/backward layers
+  fNm = 0;   // Nr. of measuring layers
+}
+
+void SolGeom::Read(const char *data)
+{
+  char strng[200];
+  int nbytes = 200;
+  FILE *fdata = fopen(data, "r");
+  if (!fdata)
+  {
+    cout << "SolGeom::GeoRead - can't open input file" << endl;
+    return;
+  }
+  Int_t tyLay;
+  char LyLabl[20];
+  float xMin;
+  float xMax;
+  float rPos;
+  float thLay;
+  float rlLay;
+  Int_t nmLay;
+  float stLayU;
+  float stLayL;
+  float sgLayU;
+  float sgLayL;
+  Int_t flLay;
+
+  fNlay = 0;
+  while (fgets(strng, nbytes, fdata) != NULL)
+  {
+    cout << strng;
+    int status = sscanf(strng, "%d %s %g %g %g %g %g %d %g %g %g %g %d",
+      &tyLay, LyLabl, &xMin, &xMax, &rPos, &thLay,
+      &rlLay, &nmLay, &stLayU, &stLayL, &sgLayU, &sgLayL, &flLay);
+    ftyLay[fNlay] = tyLay;
+    fLyLabl[fNlay] = LyLabl;
+    fxMin[fNlay] = xMin;
+    fxMax[fNlay] = xMax;
+    frPos[fNlay] = rPos;
+    fthLay[fNlay] = thLay;
+    frlLay[fNlay] = rlLay;
+    fnmLay[fNlay] = nmLay;
+    fstLayU[fNlay] = stLayU;
+    fstLayL[fNlay] = stLayL;
+    fsgLayU[fNlay] = sgLayU;
+    fsgLayL[fNlay] = sgLayL;
+    fflLay[fNlay] = flLay;
+
+    fNlay++;
+    if (tyLay == 1) fBlay++;
+    if (tyLay == 2) fFlay++;
+    if (flLay == 1) fNm++;
+
+  }
+
+  cout << "SolGeom::GeoRead completed with " << fNlay << " layers input" << endl;
+}
+
+SolGeom::~SolGeom()
+{
+  delete[] ftyLay;
+  delete[] fxMin;
+  delete[] fxMax;
+  delete[] frPos;
+  delete[] fthLay;
+  delete[] frlLay;
+  delete[] fnmLay;
+  delete[] fstLayU;
+  delete[] fstLayL;
+  delete[] fsgLayU;
+  delete[] fsgLayL;
+  delete[] fflLay;
+}
Index: external/TrackCovariance/SolGeom.h
===================================================================
--- external/TrackCovariance/SolGeom.h	(revision ff9fb2d9367d146a8325d4b77b9a9a4c0bd2ada8)
+++ external/TrackCovariance/SolGeom.h	(revision ff9fb2d9367d146a8325d4b77b9a9a4c0bd2ada8)
@@ -0,0 +1,65 @@
+#ifndef G__SOLGEOM_H
+#define G__SOLGEOM_H
+
+#include <RtypesCore.h>
+#include <TString.h>
+
+// Class to create geometry for solenoid geometry
+
+class SolGeom{
+  //
+  // Units are m
+  //
+private:
+  const Int_t fNlMax = 200; // Maximum number of layers
+
+  // B field
+  Double_t fB; // B field in Tesla
+
+  // Barrel layer properties
+  Int_t fNlay;       // Total number of layers
+  Int_t fBlay;       // Number of barrel layers
+  Int_t fFlay;       // Number of forward/backward layers
+  Int_t fNm;         // Nr. measurement layers
+  Int_t *ftyLay;     // Layer type 1 = R (barrel) or 2 = z (forward/backward)
+  TString *fLyLabl;  // Layer label
+  // Barrel: PIPE, VTXLOW, VTXHIGH, DCHCANI, DCH, DCHCANO, BSILWRP, MAG, BPRESH
+  // Fw/Bw: VTXDSK, DCHWALL, FSILWRP, FRAD, FPRESH
+  Double_t *fxMin;   // Minimum dimension z for barrel  or R for forward
+  Double_t *fxMax;   // Maximum dimension z for barrel  or R for forward
+  Double_t *frPos;   // R/z location of layer
+  Double_t *fthLay;  // Thickness (meters)
+  Double_t *frlLay;  // Radiation length (meters)
+  Int_t    *fnmLay;  // Number of measurements in layers (1D or 2D)
+  Double_t *fstLayU; // Stereo angle (rad) - 0(pi/2) = axial(z) layer - Upper side
+  Double_t *fstLayL; // Stereo angle (rad) - 0(pi/2) = axial(z) layer - Lower side
+  Double_t *fsgLayU; // Resolution Upper side (meters) - 0 = no measurement
+  Double_t *fsgLayL; // Resolution Lower side (meters) - 0 = no measurement
+  Bool_t   *fflLay;  // measurement flag = T, scattering only = F
+
+public:
+  SolGeom();
+  ~SolGeom();
+
+  void Read(const char *data);
+
+  Double_t B()                     { return fB; }
+  Int_t    Nl()                    { return fNlay; }
+  Int_t    Nm()                    { return fNm; }
+  Int_t    NBl()                   { return fBlay; }
+  TString  lLabl(Int_t nlayer)     { return fLyLabl[nlayer]; }
+  Int_t    lTyp(Int_t nlayer)      { return ftyLay[nlayer]; }
+  Double_t lxMin(Int_t nlayer)     { return fxMin[nlayer]; }
+  Double_t lxMax(Int_t nlayer)     { return fxMax[nlayer]; }
+  Double_t lPos(Int_t nlayer)      { return frPos[nlayer]; }
+  Double_t lTh(Int_t nlayer)       { return fthLay[nlayer]; }
+  Double_t lX0(Int_t nlayer)       { return frlLay[nlayer]; }
+  Int_t    lND(Int_t nlayer)       { return fnmLay[nlayer]; }
+  Double_t lStU(Int_t nlayer)      { return fstLayU[nlayer]; }
+  Double_t lStL(Int_t nlayer)      { return fstLayL[nlayer]; }
+  Double_t lSgU(Int_t nlayer)      { return fsgLayU[nlayer]; }
+  Double_t lSgL(Int_t nlayer)      { return fsgLayL[nlayer]; }
+  Bool_t   isMeasure(Int_t nlayer) { return fflLay[nlayer]; }
+};
+
+#endif
Index: external/TrackCovariance/SolGridCov.cc
===================================================================
--- external/TrackCovariance/SolGridCov.cc	(revision ff9fb2d9367d146a8325d4b77b9a9a4c0bd2ada8)
+++ external/TrackCovariance/SolGridCov.cc	(revision ff9fb2d9367d146a8325d4b77b9a9a4c0bd2ada8)
@@ -0,0 +1,157 @@
+#include <iostream>
+
+#include <TMath.h>
+#include <TVectorD.h>
+#include <TMatrixDSym.h>
+#include <TDecompChol.h>
+#include <TMatrixDSymEigen.h>
+
+#include "SolGridCov.h"
+#include "SolGeom.h"
+#include "SolTrack.h"
+
+using namespace std;
+
+SolGridCov::SolGridCov()
+{
+  // Define pt-polar angle grid
+  fNpt = 22;
+  fPta.ResizeTo(fNpt);
+  Double_t p[] = { 0.1, 0.2, 0.5, 0.7, 1., 2., 3., 4., 6., 8., 10., 15.,
+                   20., 25., 30., 40., 50., 60., 80., 100., 150., 200. };
+  for (Int_t ip = 0; ip < fNpt; ip++) fPta(ip) = p[ip];
+
+  fNang = 13;
+  fAnga.ResizeTo(fNang);
+  Double_t a[] = { 10., 15., 20., 25., 30., 35., 40., 45., 50., 60., 70., 80., 90. };
+  for (Int_t ia = 0; ia < fNang; ia++) fAnga(ia) = a[ia];
+  fCov = new TMatrixDSym[fNpt * fNang];
+  for (Int_t ip = 0; ip < fNpt; ip++)
+  {
+    for (Int_t ia = 0; ia < fNang; ia++) fCov[ip * fNang + ia].ResizeTo(5, 5);
+  }
+}
+
+SolGridCov::~SolGridCov()
+{
+  delete[] fCov;
+}
+
+void SolGridCov::Calc(SolGeom *G)
+{
+  TVectorD pta = fPta;
+  TVectorD anga = fAnga;
+  Bool_t Res = kTRUE; Bool_t MS = kTRUE; // Resolution and multiple scattering flags
+  for (Int_t ip = 0; ip < fNpt; ip++) // Loop on pt grid
+  {
+    Int_t ipt = TMath::Nint(10 * pta(ip));
+    for (Int_t ia = 0; ia < fNang; ia++) // Loop on angle grid
+    {
+      Double_t th = TMath::Pi() * (anga(ia)) / 180.;
+      Double_t x[3], p[3];
+      x[0] = 0; x[1] = 0; x[2] = 0; // Set origin
+      p[0] = pta(ip); p[1] = 0; p[2] = pta(ip) / TMath::Tan(th);
+      //
+      SolTrack *tr = new SolTrack(x, p, G); // Initialize track
+      tr->CovCalc(Res, MS); // Calculate covariance
+      fCov[ip * fNang + ia] = tr->Cov(); // Get covariance
+    }
+  }
+}
+// Find bin in grid
+Int_t SolGridCov::GetMinIndex(Double_t xval, Int_t N, TVectorD x)
+{
+  Int_t min = -1; // default for xval below the lower limit
+  if (xval < x(0))return min;
+  if (xval>x(N - 1)){ min = N; return min; }
+  for (Int_t i = 0; i < N; i++) if (xval>x(i))min = i;
+  return min;
+}
+// Force positive definitness in normalized matrix
+TMatrixDSym SolGridCov::MakePosDef(TMatrixDSym NormMat)
+{
+  // Input: symmetric matrix with 1's on diagonal
+  // Output: positive definite matrix with 1's on diagonal
+
+  // Default return value
+  TMatrixDSym rMatN = NormMat;
+  // Check the diagonal
+  Bool_t Check = kFALSE;
+  Int_t Size = NormMat.GetNcols();
+  for (Int_t i = 0; i < Size; i++)if (TMath::Abs(NormMat(i, i) - 1.0)>1.0E-15)Check = kTRUE;
+  if (Check)
+  {
+    cout << "SolGridCov::MakePosDef: input matrix doesn ot have 1 on diagonal. Abort." << endl;
+    return rMatN;
+  }
+  // Diagonalize matrix
+  TMatrixDSymEigen Eign(NormMat);
+  TMatrixD U = Eign.GetEigenVectors();
+  TVectorD lambda = Eign.GetEigenValues();
+  // Reset negative eigenvalues to small positive value
+  TMatrixDSym D(Size); D.Zero(); Double_t eps = 1.0e-13;
+  for (Int_t i = 0; i < Size; i++)
+  {
+    D(i, i) = lambda(i);
+    if (lambda(i) <= 0) D(i, i) = eps;
+  }
+  // Rebuild matrix
+  TMatrixD Ut(TMatrixD::kTransposed, U);
+  TMatrixD rMat = (U*D)*Ut;       // Now it is positive defite
+  // Restore all ones on diagonal
+  for (Int_t i1 = 0; i1 < Size; i1++)
+  {
+    Double_t rn1 = TMath::Sqrt(rMat(i1, i1));
+    for (Int_t i2 = 0; i2 <= i1; i2++)
+    {
+      Double_t rn2 = TMath::Sqrt(rMat(i2, i2));
+      rMatN(i1, i2) = 0.5*(rMat(i1, i2) + rMat(i2, i1)) / (rn1*rn2);
+      rMatN(i2, i1) = rMatN(i1, i2);
+    }
+  }
+  return rMatN;
+}
+// Interpolate covariance matrix: Bi-linear interpolation
+TMatrixDSym SolGridCov::GetCov(Double_t pt, Double_t ang)
+{
+  // pt in GeV and ang in degrees
+  Int_t minPt = GetMinIndex(pt, fNpt, fPta);
+  if (minPt == -1)minPt = 0;
+  if (minPt == fNpt - 1)minPt = fNpt - 2;
+  Double_t dpt = fPta(minPt + 1) - fPta(minPt);
+  // Put ang in 0-90 range
+  ang = TMath::Abs(ang);
+  while (ang > 90.)ang -= 90.;  // Needs to be fixed
+  Int_t minAng = GetMinIndex(ang, fNang, fAnga);
+  if (minAng == -1)minAng = 0;
+  if (minAng == fNang - 1)minAng = fNang - 2;
+  Double_t dang = fAnga(minAng + 1) - fAnga(minAng);
+  //
+  Double_t tpt = (pt - fPta(minPt)) / dpt;
+  Double_t tang = (ang - fAnga(minAng)) / dang;
+  //
+  TMatrixDSym C11 = fCov[minPt * fNang + minAng];
+  TMatrixDSym C12 = fCov[minPt * fNang + minAng + 1];
+  TMatrixDSym C21 = fCov[(minPt + 1) * fNang + minAng];
+  TMatrixDSym C22 = fCov[(minPt + 1) * fNang + minAng + 1];
+  TMatrixDSym Cv = ((1-tpt) * (1-tang)) * C11 +
+                   ((1-tpt) *    tang ) * C12 +
+                   (   tpt  * (1-tang)) * C21 +
+                   (   tpt  *    tang ) * C22;
+  // Check for positive definiteness
+  TMatrixDSym CvN = Cv;
+  TMatrixDSym DCvInv(5); DCvInv.Zero();
+  for (Int_t id = 0; id < 5; id++) DCvInv(id, id) = 1.0 / TMath::Sqrt(Cv(id, id));
+  CvN.Similarity(DCvInv); // Normalize diagonal to 1
+  TDecompChol Chl(CvN);
+  if (!Chl.Decompose())
+  {
+    cout << "SolGridCov::GetCov: Interpolated matrix not positive definite. Recovering ...." << endl;
+    TMatrixDSym rCv = MakePosDef(CvN); CvN = rCv;
+    TMatrixDSym DCv(5); DCv.Zero();
+    for (Int_t id = 0; id < 5; id++) DCv(id, id) = TMath::Sqrt(Cv(id, id));
+    Cv = CvN.Similarity(DCv); // Normalize diagonal to 1
+  }
+
+  return Cv;
+}
Index: external/TrackCovariance/SolGridCov.h
===================================================================
--- external/TrackCovariance/SolGridCov.h	(revision ff9fb2d9367d146a8325d4b77b9a9a4c0bd2ada8)
+++ external/TrackCovariance/SolGridCov.h	(revision ff9fb2d9367d146a8325d4b77b9a9a4c0bd2ada8)
@@ -0,0 +1,36 @@
+#ifndef G__SOLGRIDCOV_H
+#define G__SOLGRIDCOV_H
+
+#include <TVectorD.h>
+#include <TMatrixDSym.h>
+
+class SolGeom;
+
+// Class to create geometry for solenoid geometry
+
+class SolGridCov{
+  // Class to handle storing and retrieving/interpolation of covariance matrices
+private:
+  Int_t fNpt;        // Number of pt points in grid
+  TVectorD fPta;     // Array of pt points in GeV
+  Int_t fNang;       // Number of angle points in grid
+  TVectorD fAnga;    // Array of angle points in degrees
+  TMatrixDSym *fCov; // Pointers to grid of covariance matrices
+  // Service routines
+  Int_t GetMinIndex(Double_t xval, Int_t N, TVectorD x); // Find bin
+  TMatrixDSym MakePosDef(TMatrixDSym NormMat); // Force positive definitness
+public:
+  SolGridCov();
+  ~SolGridCov();
+
+  void Calc(SolGeom *G);
+
+  // Covariance interpolation
+  Double_t GetMinPt()  { return fPta(0); }
+  Double_t GetMaxPt()  { return fPta(fNpt - 1); }
+  Double_t GetMinAng() { return fAnga(0); }
+  Double_t GetMaxAng() { return fAnga(fNang - 1); }
+  TMatrixDSym GetCov(Double_t pt, Double_t ang);
+};
+
+#endif
Index: external/TrackCovariance/SolTrack.cc
===================================================================
--- external/TrackCovariance/SolTrack.cc	(revision ff9fb2d9367d146a8325d4b77b9a9a4c0bd2ada8)
+++ external/TrackCovariance/SolTrack.cc	(revision ff9fb2d9367d146a8325d4b77b9a9a4c0bd2ada8)
@@ -0,0 +1,476 @@
+#include <iostream>
+
+#include <TString.h>
+#include <TMath.h>
+#include <TMatrixD.h>
+#include <TMatrixDSym.h>
+#include <TDecompChol.h>
+#include <TMatrixDSymEigen.h>
+
+#include "SolGeom.h"
+#include "SolTrack.h"
+
+using namespace std;
+
+SolTrack::SolTrack(Double_t *x, Double_t *p, SolGeom *G)
+{
+  fG = G;
+  // Store momentum
+  fp[0] = p[0]; fp[1] = p[1]; fp[2] = p[2];
+  Double_t px = p[0]; Double_t py = p[1]; Double_t pz = p[2];
+  fx[0] = x[0]; fx[1] = x[1]; fx[2] = x[2];
+  Double_t xx = x[0]; Double_t yy = x[1]; Double_t zz = x[2];
+  // Store parameters
+  Double_t pt = TMath::Sqrt(px*px + py*py);
+  Double_t Charge = 1.0;            // Don't worry about charge for now
+  Double_t a = -Charge*G->B()*0.2998;     // Normalized speed of light
+  Double_t C = a / (2 * pt);          // pt in GeV, B in Tesla, C in meters
+  Double_t r2 = xx*xx + yy*yy;
+  Double_t cross = xx*py - yy*px;
+  Double_t T = TMath::Sqrt(pt*pt - 2 * a*cross + a*a*r2);
+  Double_t phi0 = TMath::ATan2((py - a*xx) / T, (px + a*yy) / T);
+  Double_t D;
+  if (pt < 10.0) D = (T - pt) / a;
+  else D = (-2 * cross + a*r2) / (T + pt);
+  Double_t B = C*TMath::Sqrt((r2 - D*D) / (1 + 2 * C*D));
+  Double_t st = TMath::ASin(B) / C;
+  Double_t ct = pz / pt;
+  Double_t z0 = zz - ct*st;
+  fpar[0] = D;
+  fpar[1] = phi0;
+  fpar[2] = C;
+  fpar[3] = z0;
+  fpar[4] = ct;
+  // Init covariances
+  fCov.ResizeTo(5, 5);
+}
+
+SolTrack::SolTrack(Double_t D, Double_t phi0, Double_t C, Double_t z0, Double_t ct, SolGeom *G)
+{
+  fG = G;
+  // Store parameters
+  fpar[0] = D;
+  fpar[1] = phi0;
+  fpar[2] = C;
+  fpar[3] = z0;
+  fpar[4] = ct;
+  // Store momentum
+  Double_t pt = G->B()*0.2998 / TMath::Abs(2 * C);
+  Double_t px = pt*TMath::Cos(phi0);
+  Double_t py = pt*TMath::Sin(phi0);
+  Double_t pz = pt*ct;
+
+  fp[0] = px; fp[1] = py; fp[2] = pz;
+  fx[0] = -D*TMath::Sin(phi0); fx[1] = D*TMath::Cos(phi0);  fx[2] = z0;
+  // Init covariances
+  fCov.ResizeTo(5, 5);
+}
+// Destructor
+SolTrack::~SolTrack()
+{
+  fCov.Clear();
+}
+// Calculate intersection with given layer
+Bool_t SolTrack::HitLayer(Int_t il, Double_t &R, Double_t &phi, Double_t &zz)
+{
+  Double_t Di = D();
+  Double_t phi0i = phi0();
+  Double_t Ci = C();
+  Double_t z0i = z0();
+  Double_t cti = ct();
+
+  R = 0; phi = 0; zz = 0;
+
+  Bool_t val = kFALSE;
+  if (fG->lTyp(il) == 1)      // Cylinder: layer at constant R
+  {
+    R = fG->lPos(il);
+    Double_t argph = (Ci*R + (1 + Ci*Di)*Di / R) / (1. + 2.*Ci*Di);
+    if (TMath::Abs(argph) < 1.0)
+    {
+      Double_t argz = Ci*TMath::Sqrt((R*R - Di*Di) / (1 + 2 * Ci*Di));
+      if (TMath::Abs(argz) < 1.0)
+      {
+        zz = z0i + cti*TMath::ASin(argz) / Ci;
+        if (zz > fG->lxMin(il) && zz < fG->lxMax(il))
+        {
+          phi = phi0i + TMath::ASin(argph);
+          val = kTRUE;
+        }
+      }
+    }
+  }
+  else if (fG->lTyp(il) == 2)   // disk: layer at constant z
+  {
+    zz = fG->lPos(il);
+    Double_t arg = Ci*(zz - z0i) / cti;
+    if (TMath::Abs(arg) < 1.0 && (zz - z0i) / cti > 0)
+    {
+      R = TMath::Sqrt(Di*Di + (1. + 2.*Ci*Di)*pow(TMath::Sin(arg), 2) / (Ci*Ci));
+      if (R > fG->lxMin(il) && R < fG->lxMax(il))
+      {
+        Double_t arg1 = (Ci*R + (1 + Ci*Di)*Di / R) / (1. + 2.*Ci*Di);
+        if (TMath::Abs(arg1) < 1.0)
+        {
+          phi = phi0i + TMath::ASin(arg1);
+          val = kTRUE;
+        }
+      }
+    }
+  }
+  //
+  return val;
+}
+// # of layers hit
+Int_t SolTrack::nHit()
+{
+  Int_t kh = 0;
+  Double_t R; Double_t phi; Double_t zz;
+  for (Int_t i = 0; i < fG->Nl(); i++)
+  if (HitLayer(i, R, phi, zz))kh++;
+
+  return kh;
+}
+// List of layers hit with intersections
+Int_t SolTrack::HitList(Int_t *&ihh, Double_t *&rhh, Double_t *&zhh)
+{
+  // Return lists of hits associated to a track including all scattering layers.
+  // Return value is the total number of measurement hits
+  // kmh = total number of measurement layers hit for given type
+  // ihh = pointer to layer number
+  // rhh = radius of hit
+  // zhh = z of hit
+
+  // ***** NB: double layers with stereo on lower layer not included
+
+  Int_t kh = 0; // Number of layers hit
+  Int_t kmh = 0; // Number of measurement layers of given type
+  for (Int_t i = 0; i < fG->Nl(); i++)
+  {
+    Double_t R; Double_t phi; Double_t zz;
+    if (HitLayer(i, R, phi, zz)) // Only barrel type layers
+    {
+      zhh[kh] = zz;
+      rhh[kh] = R;
+      ihh[kh] = i;
+      if (fG->isMeasure(i))kmh++;
+      kh++;
+    }
+  }
+
+  return kmh;
+}
+// Covariance matrix estimation
+void SolTrack::CovCalc(Bool_t Res, Bool_t MS)
+{
+  // Input flags:
+  //        Res = .TRUE. turn on resolution effects/Use standard resolutions
+  //            .FALSE. set all resolutions to 0
+  //        MS  = .TRUE. include Multiple Scattering
+  // Assumptions:
+  // 1. Measurement layers can do one or two measurements
+  // 2. On disks at constant z:
+  //    - Upper side measurement is phi
+  //    - Lower side measurement is R
+
+  // Fill list of layers hit
+  Int_t ntry = 0;
+  Int_t ntrymax = 0;
+  Int_t Nhit = nHit(); // Total number of layers hit
+  Double_t *zhh = new Double_t[Nhit]; // z of hit
+  Double_t *rhh = new Double_t[Nhit]; // r of hit
+  Double_t *dhh = new Double_t[Nhit]; // distance of hit from origin
+  Int_t    *ihh = new Int_t[Nhit];    // true index of layer
+  Int_t kmh; // Number of measurement layers hit
+
+  kmh = HitList(ihh, rhh, zhh); // hit layer list
+  Int_t mTot = 0; // Total number of measurements
+  for (Int_t i = 0; i < Nhit; i++)
+  {
+    dhh[i] = TMath::Sqrt(rhh[i] * rhh[i] + zhh[i] * zhh[i]);
+    if (fG->isMeasure(ihh[i])) mTot += fG->lND(ihh[i]); // Count number of measurements
+  }
+  // Order hit list by increasing distance from origin
+  Int_t *hord = new Int_t[Nhit]; // hit order by increasing distance from origin
+  TMath::Sort(Nhit, dhh, hord, kFALSE); // Order by increasing distance from origin
+  Double_t *zh = new Double_t[Nhit]; // d-ordered z of hit
+  Double_t *rh = new Double_t[Nhit]; // d-ordered r of hit
+  Int_t    *ih = new Int_t[Nhit];    // d-ordered true index of layer
+  for (Int_t i = 0; i < Nhit; i++)
+  {
+    Int_t il = hord[i]; // Hit layer numbering
+    zh[i] = zhh[il];
+    rh[i] = rhh[il];
+    ih[i] = ihh[il];
+  }
+  // Store interdistances and multiple scattering angles
+  Double_t sn2t = 1.0 / (1 + ct()*ct()); //sin^2 theta of track
+  Double_t cs2t = 1.0 - sn2t;            //cos^2 theta
+  Double_t snt = TMath::Sqrt(sn2t);      // sin theta
+  Double_t cst = TMath::Sqrt(cs2t);      // cos theta
+  Double_t px0 = pt() * TMath::Cos(phi0()); // Momentum at minimum approach
+  Double_t py0 = pt() * TMath::Sin(phi0());
+  Double_t pz0 = pt() * ct();
+  TMatrixDSym dik(Nhit);  dik.Zero(); // Distances between layers
+  Double_t *thms = new Double_t[Nhit]; // Scattering angles/plane
+  Double_t *cs = new Double_t[Nhit]; // Cosine of angle with layer normal
+  for (Int_t ii = 0; ii < Nhit; ii++) // Hit layer loop
+  {
+    Int_t i = ih[ii]; // Get true layer number
+    Double_t B = C()*TMath::Sqrt((rh[ii] * rh[ii] - D()*D()) / (1 + 2 * C()*D()));
+    Double_t pxi = px0*(1-2*B*B)-2*py0*B*TMath::Sqrt(1-B*B); // Momentum at scattering layer
+    Double_t pyi = py0*(1-2*B*B)+2*px0*B*TMath::Sqrt(1-B*B);
+    Double_t pzi = pz0;
+    Double_t ArgRp = (rh[ii]*C() + (1 + C() * D())*D() / rh[ii]) / (1 + 2 * C()*D());
+    Double_t phi = phi0() + TMath::ASin(ArgRp);
+    Double_t nx = TMath::Cos(phi); // Barrel layer normal
+    Double_t ny = TMath::Sin(phi);
+    Double_t nz = 0.0;
+    if (fG->lTyp(i) == 2) // this is Z layer
+    {
+      nx = 0.0;
+      ny = 0.0;
+      nz = 1.0;
+    }
+    Double_t corr = (pxi*nx + pyi * ny + pzi * nz) / p();
+    cs[ii] = corr;
+    Double_t Rlf = fG->lTh(i) / (corr*fG->lX0(i)); // Rad. length fraction
+    thms[ii] = 0.0136*TMath::Sqrt(Rlf)*(1.0 + 0.038*TMath::Log(Rlf)) / p(); // MS angle
+    if (!MS)thms[ii] = 0;
+    //
+    for (Int_t kk = 0; kk < ii; kk++) // Fill distances between layers
+    {
+      Double_t Ci = C();
+      dik(ii, kk) = (TMath::ASin(Ci*rh[ii])-TMath::ASin(Ci*rh[kk]))/(Ci*snt);
+      dik(kk, ii) = dik(ii, kk);
+    }
+    // Store momentum components for resolution correction cosines
+    Double_t *pRi = new Double_t[Nhit];
+    pRi[ii] = TMath::Abs(pxi * TMath::Cos(phi) + pyi * TMath::Sin(phi)); // Radial component
+    Double_t *pPhi = new Double_t[Nhit];
+    pPhi[ii] = TMath::Abs(pxi * TMath::Sin(phi) - pyi * TMath::Cos(phi)); // Phi component
+  }
+  // Fill measurement covariance
+  Int_t *mTl = new Int_t[mTot]; // Pointer from measurement number to true layer number
+  TMatrixDSym Sm(mTot); Sm.Zero(); // Measurement covariance
+  TMatrixD Rm(mTot, 5); // Derivative matrix
+  Int_t im = 0;
+  // Fill derivatives and error matrix with MS
+  Double_t AngMax = 90.; Double_t AngMx = AngMax * TMath::Pi() / 180.;
+  Double_t csMin = TMath::Cos(AngMx); // Set maximum angle wrt normal
+  //
+  for (Int_t ii = 0; ii < Nhit; ii++)
+  {
+    Int_t i = ih[ii]; // True layer number
+    Int_t ityp  = fG->lTyp(i); // Layer type Barrel or Z
+    Int_t nmeai = fG->lND(i); // # measurements in layer
+    if (fG->isMeasure(i) && nmeai >0 && cs[ii] > csMin)
+    {
+      Double_t Di    = D(); // Get true track parameters
+      Double_t phi0i = phi0();
+      Double_t Ci    = C();
+      Double_t z0i   = z0();
+      Double_t cti   = ct();
+      //
+      Double_t Ri    = rh[ii];
+      Double_t ArgRp = (Ri*Ci + (1 + Ci * Di)*Di / Ri) / (1 + 2 * Ci*Di);
+      Double_t ArgRz = Ci * TMath::Sqrt((Ri*Ri - Di * Di) / (1 + 2 * Ci*Di));
+      TVectorD dRphi(5); dRphi.Zero(); // R-phi derivatives @ const. R
+      TVectorD dRz(5); dRz.Zero(); // z derivatives @ const. R
+      // Derivative overflow protection
+      Double_t dMin = 0.8;
+      dRphi(0) = (1 - 2 * Ci*Ci*Ri*Ri) /
+        TMath::Max(dMin,TMath::Sqrt(1 - ArgRp * ArgRp)); // D derivative
+      dRphi(1) = Ri; // phi0 derivative
+      dRphi(2) = Ri * Ri /
+        TMath::Max(dMin,TMath::Sqrt(1 - ArgRp * ArgRp)); // C derivative
+      dRphi(3) = 0.0; // z0 derivative
+      dRphi(4) = 0.0; // cot(theta) derivative
+
+      dRz(0) = -cti * Di /
+        (Ri*TMath::Max(dMin,TMath::Sqrt(1 - Ci * Ci*Ri*Ri))); // D
+      dRz(1) = 0.0; // Phi0
+      dRz(2) = cti * (Ri*Ci / TMath::Sqrt(1-Ri*Ri*Ci*Ci) -
+        TMath::ASin(Ri*Ci)) / (Ci*Ci); // C
+      dRz(3) = 1.0; // Z0
+      dRz(4) = TMath::ASin(ArgRz) / Ci; // Cot(theta)
+
+      for (Int_t nmi = 0; nmi < nmeai; nmi++)
+      {
+        mTl[im] = i;
+        Double_t stri = 0;
+        Double_t sig = 0;
+        if (nmi + 1 == 1) // Upper layer measurements
+        {
+          stri = fG->lStU(i); // Stereo angle
+          Double_t csa = TMath::Cos(stri);
+          Double_t ssa = TMath::Sin(stri);
+          sig = fG->lSgU(i); // Resolution
+          if (ityp == 1) // Barrel type layer (Measure R-phi, stereo or z at const. R)
+          {
+            // Almost exact solution (CD<<1, D<<R)
+            Rm(im, 0) = csa * dRphi(0) - ssa * dRz(0); // D derivative
+            Rm(im, 1) = csa * dRphi(1) - ssa * dRz(1); // phi0 derivative
+            Rm(im, 2) = csa * dRphi(2) - ssa * dRz(2); // C derivative
+            Rm(im, 3) = csa * dRphi(3) - ssa * dRz(3); // z0 derivative
+            Rm(im, 4) = csa * dRphi(4) - ssa * dRz(4); // cot(theta) derivative
+          }
+          if (ityp == 2) // Z type layer (Measure phi at const. Z)
+          {
+            Rm(im, 0) = 1.0; // D derivative
+            Rm(im, 1) = rh[ii]; // phi0 derivative
+            Rm(im, 2) = rh[ii] * rh[ii]; // C derivative
+            Rm(im, 3) = 0; // z0 derivative
+            Rm(im, 4) = 0; // cot(theta) derivative
+          }
+        }
+        if (nmi + 1 == 2) // Lower layer measurements
+        {
+          stri = fG->lStL(i); // Stereo angle
+          Double_t csa = TMath::Cos(stri);
+          Double_t ssa = TMath::Sin(stri);
+          sig = fG->lSgL(i); // Resolution
+          if (ityp == 1) // Barrel type layer (measure R-phi, stereo or z at const. R)
+          {
+            // Almost exact solution (CD<<1, D<<R)
+            Rm(im, 0) = csa * dRphi(0) - ssa * dRz(0); // D derivative
+            Rm(im, 1) = csa * dRphi(1) - ssa * dRz(1); // phi0 derivative
+            Rm(im, 2) = csa * dRphi(2) - ssa * dRz(2); // C derivative
+            Rm(im, 3) = csa * dRphi(3) - ssa * dRz(3); // z0 derivative
+            Rm(im, 4) = csa * dRphi(4) - ssa * dRz(4); // cot(theta) derivative
+          }
+          if (ityp == 2) // Z type layer (Measure R at const. z)
+          {
+            Rm(im, 0) = 0; // D derivative
+            Rm(im, 1) = 0; // phi0 derivative
+            Rm(im, 2) = 0; // C derivative
+            Rm(im, 3) = -1.0 / ct(); // z0 derivative
+            Rm(im, 4) = -rh[ii] / ct(); // cot(theta) derivative
+          }
+        }
+        // Derivative calculation completed
+        // Now calculate measurement error matrix
+        Int_t km = 0;
+        for (Int_t kk = 0; kk <= ii; kk++)
+        {
+          Int_t k = ih[kk]; // True layer number
+          Int_t ktyp = fG->lTyp(k); // Layer type Barrel or
+          Int_t nmeak = fG->lND(k); // # measurements in layer
+          if (fG->isMeasure(k) && nmeak > 0 &&cs[kk] > csMin)
+          {
+            for (Int_t nmk = 0; nmk < nmeak; nmk++)
+            {
+              Double_t strk = 0;
+              if (nmk + 1 == 1) strk = fG->lStU(k); // Stereo angle
+              if (nmk + 1 == 2) strk = fG->lStL(k); // Stereo angle
+              if (im == km && Res) Sm(im, km) += sig*sig; // Detector resolution on diagonal
+              //
+              // Loop on all layers below for MS contributions
+              for (Int_t jj = 0; jj < kk; jj++)
+              {
+                Double_t di = dik(ii, jj);
+                Double_t dk = dik(kk, jj);
+                Double_t ms = thms[jj];
+                Double_t msk = ms; Double_t msi = ms;
+                if (ityp == 1) msi = ms / snt; // Barrel
+                else if (ityp == 2) msi = ms / cst; // Disk
+                if (ktyp == 1) msk = ms / snt; // Barrel
+                else if (ktyp == 2) msk = ms / cst; // Disk
+                Double_t ci = TMath::Cos(stri); Double_t si = TMath::Sin(stri);
+                Double_t ck = TMath::Cos(strk); Double_t sk = TMath::Sin(strk);
+                Sm(im, km) += di*dk*(ci*ck*ms*ms + si*sk*msi*msk); // Ms contribution
+              }
+              Sm(km, im) = Sm(im, km);
+              km++;
+            }
+          }
+        }
+        im++; mTot = im;
+      }
+    }
+  }
+  Sm.ResizeTo(mTot, mTot);
+  Rm.ResizeTo(mTot, 5);
+
+  // Calculate covariance from derivatives and measurement error matrix
+  TMatrixDSym DSmInv(mTot); DSmInv.Zero();
+  for (Int_t id = 0; id < mTot; id++) DSmInv(id, id) = 1.0 / TMath::Sqrt(Sm(id, id));
+  TMatrixDSym SmN = Sm.Similarity(DSmInv);  // Normalize diagonal to 1
+  // Protected matrix inversions
+  TDecompChol Chl(SmN);
+  TMatrixDSym SmNinv = SmN;
+  if (Chl.Decompose())
+  {
+    Bool_t OK;
+    SmNinv = Chl.Invert(OK);
+  }
+  else
+  {
+    cout << "SolTrack::CovCalc: Error matrix not positive definite. Recovering ...." << endl;
+    if (ntry < ntrymax)
+    {
+      SmNinv.Print();
+      ntry++;
+    }
+    TMatrixDSym rSmN = MakePosDef(SmN); SmN = rSmN;
+    TDecompChol rChl(SmN);
+    SmNinv = SmN;
+    Bool_t OK = rChl.Decompose();
+    SmNinv    = rChl.Invert(OK);
+  }
+  Sm = SmNinv.Similarity(DSmInv); // Error matrix inverted
+  TMatrixDSym H = Sm.SimilarityT(Rm); // Calculate half Hessian
+  // Normalize before inversion
+  const Int_t Npar = 5;
+  TMatrixDSym DHinv(Npar); DHinv.Zero();
+  for (Int_t i = 0; i < Npar; i++)DHinv(i, i) = 1.0 / TMath::Sqrt(H(i, i));
+  TMatrixDSym Hnrm = H.Similarity(DHinv);
+  // Invert and restore
+  Hnrm.Invert();
+  fCov = Hnrm.Similarity(DHinv);
+}
+
+// Force positive definitness in normalized matrix
+TMatrixDSym SolTrack::MakePosDef(TMatrixDSym NormMat)
+{
+  // Input: symmetric matrix with 1's on diagonal
+  // Output: positive definite matrix with 1's on diagonal
+
+  // Default return value
+  TMatrixDSym rMatN = NormMat;
+  // Check the diagonal
+  Bool_t Check = kFALSE;
+  Int_t Size = NormMat.GetNcols();
+  for (Int_t i = 0; i < Size; i++)if (TMath::Abs(NormMat(i, i) - 1.0)>1.0E-15)Check = kTRUE;
+  if (Check)
+  {
+    cout << "SolTrack::MakePosDef: input matrix doesn't have 1 on diagonal. Abort." << endl;
+    return rMatN;
+  }
+  // Diagonalize matrix
+  TMatrixDSymEigen Eign(NormMat);
+  TMatrixD U = Eign.GetEigenVectors();
+  TVectorD lambda = Eign.GetEigenValues();
+  // Reset negative eigenvalues to small positive value
+  TMatrixDSym D(Size); D.Zero(); Double_t eps = 1.0e-13;
+  for (Int_t i = 0; i < Size; i++)
+  {
+    D(i, i) = lambda(i);
+    if (lambda(i) <= 0) D(i, i) = eps;
+  }
+  // Rebuild matrix
+  TMatrixD Ut(TMatrixD::kTransposed, U);
+  TMatrixD rMat = (U*D)*Ut;       // Now it is positive defite
+  // Restore all ones on diagonal
+  for (Int_t i1 = 0; i1 < Size; i1++)
+  {
+    Double_t rn1 = TMath::Sqrt(rMat(i1, i1));
+    for (Int_t i2 = 0; i2 <= i1; i2++)
+    {
+      Double_t rn2 = TMath::Sqrt(rMat(i2, i2));
+      rMatN(i1, i2) = 0.5*(rMat(i1, i2) + rMat(i2, i1)) / (rn1*rn2);
+      rMatN(i2, i1) = rMatN(i1, i2);
+    }
+  }
+  return rMatN;
+}
Index: external/TrackCovariance/SolTrack.h
===================================================================
--- external/TrackCovariance/SolTrack.h	(revision ff9fb2d9367d146a8325d4b77b9a9a4c0bd2ada8)
+++ external/TrackCovariance/SolTrack.h	(revision ff9fb2d9367d146a8325d4b77b9a9a4c0bd2ada8)
@@ -0,0 +1,64 @@
+#ifndef G__SOLTRK_H
+#define G__SOLTRK_H
+
+#include <TMath.h>
+#include <TMatrixDSym.h>
+
+class SolGeom;
+
+// Class to store track information
+// Assumes that the geometry has been initialized
+class SolTrack{
+  // Track handling class
+  // Assume tracks originate from (0,0) for the time being
+private:
+  Int_t fNl;        // Actual number of layers
+  SolGeom *fG;      // Geometry
+  Double_t fp[3];   // px, py, pz momentum
+  Double_t fx[3];   //  x,  y,  z track origin
+  Double_t fpar[5]; // D, phi0, C, z0, cot(theta)
+
+  TMatrixDSym fCov; // Full covariance matrix
+
+public:
+  // Constructors
+  SolTrack(Double_t *x, Double_t *p, SolGeom *G);
+  SolTrack(Double_t D, Double_t phi0, Double_t C, Double_t z0, Double_t ct, SolGeom *G);
+  // Destructor
+  ~SolTrack();
+  // Accessors
+  // Position (at minimum approach)
+  Double_t x() { return fx[0]; }
+  Double_t y() { return fx[1]; }
+  Double_t z() { return fx[2]; }
+  // Momentum  (at minimum approach)
+  Double_t px() { return fp[0]; }
+  Double_t py() { return fp[1]; }
+  Double_t pz() { return fp[2]; }
+  Double_t pt() { return TMath::Sqrt(fp[0] * fp[0] + fp[1] * fp[1]); }
+  Double_t p()  { return TMath::Sqrt(fp[0] * fp[0] + fp[1] * fp[1] + fp[2] * fp[2]); }
+  // Track parameters
+  Double_t D()    { return fpar[0]; }
+  Double_t phi0() { return fpar[1]; }
+  Double_t C()    { return fpar[2]; }
+  Double_t z0()   { return fpar[3]; }
+  Double_t ct()   { return fpar[4]; }
+  // Covariance
+  TMatrixDSym Cov()   { return fCov; }
+  // Track parameter covariance calculation
+  void CovCalc(Bool_t Res, Bool_t MS);
+  // Parameter errors
+  Double_t s_D()    { return TMath::Sqrt(fCov(0, 0)); }
+  Double_t s_phi0() { return TMath::Sqrt(fCov(1, 1)); }
+  Double_t s_C()    { return TMath::Sqrt(fCov(2, 2)); }
+  Double_t s_pt()   { return 2 * s_C()*pt() / (0.2998*fG->B()); } // Dpt/pt
+  Double_t s_z0()   { return TMath::Sqrt(fCov(3, 3)); }
+  Double_t s_ct()   { return TMath::Sqrt(fCov(4, 4)); }
+  // Track hit management
+  Int_t nHit();
+  Bool_t HitLayer(Int_t Layer, Double_t &R, Double_t &phi, Double_t &zz);
+  Int_t HitList(Int_t *&ihh, Double_t *&rhh, Double_t *&zhh);
+  // Make normalized matrix positive definite
+  TMatrixDSym MakePosDef(TMatrixDSym NormMat);
+};
+#endif
Index: modules/ModulesLinkDef.h
===================================================================
--- modules/ModulesLinkDef.h	(revision 3051358a56b3d611ea993535575ffebc4f158441)
+++ modules/ModulesLinkDef.h	(revision ff9fb2d9367d146a8325d4b77b9a9a4c0bd2ada8)
@@ -36,4 +36,5 @@
 #include "modules/MomentumSmearing.h"
 #include "modules/TrackSmearing.h"
+#include "modules/TrackCovariance.h"
 #include "modules/ImpactParameterSmearing.h"
 #include "modules/TimeSmearing.h"
@@ -91,4 +92,5 @@
 #pragma link C++ class MomentumSmearing+;
 #pragma link C++ class TrackSmearing+;
+#pragma link C++ class TrackCovariance+;
 #pragma link C++ class ImpactParameterSmearing+;
 #pragma link C++ class TimeSmearing+;
Index: modules/TrackCovariance.cc
===================================================================
--- modules/TrackCovariance.cc	(revision ff9fb2d9367d146a8325d4b77b9a9a4c0bd2ada8)
+++ modules/TrackCovariance.cc	(revision ff9fb2d9367d146a8325d4b77b9a9a4c0bd2ada8)
@@ -0,0 +1,110 @@
+/*
+ *  Delphes: a framework for fast simulation of a generic collider experiment
+ *  Copyright (C) 2020  Universite catholique de Louvain (UCLouvain), Belgium
+ *
+ *  This program is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/** \class TrackCovariance
+ *
+ *  Smears track parameters according to appropriate covariance matrix.
+ *
+ *  \author P. Demin - UCLouvain, Louvain-la-Neuve
+ *
+ */
+
+#include "modules/TrackCovariance.h"
+
+#include "classes/DelphesClasses.h"
+
+#include "TrackCovariance/SolGeom.h"
+#include "TrackCovariance/SolGridCov.h"
+#include "TrackCovariance/ObsTrk.h"
+
+#include "TLorentzVector.h"
+#include "TMath.h"
+#include "TObjArray.h"
+
+//------------------------------------------------------------------------------
+
+TrackCovariance::TrackCovariance() :
+  fGeometry(0), fCovariance(0), fItInputArray(0)
+{
+  fGeometry = new SolGeom();
+  fCovariance = new SolGridCov();
+}
+
+//------------------------------------------------------------------------------
+
+TrackCovariance::~TrackCovariance()
+{
+  if(fGeometry) delete fGeometry;
+  if(fCovariance) delete fCovariance;
+}
+
+//------------------------------------------------------------------------------
+
+void TrackCovariance::Init()
+{
+  fBz = GetDouble("Bz", 0.0);
+  fGeometry->Read(GetString("DetectorGeometry", ""));
+
+  fCovariance->Calc(fGeometry);
+
+  // import input array
+
+  fInputArray = ImportArray(GetString("InputArray", "TrackMerger/tracks"));
+  fItInputArray = fInputArray->MakeIterator();
+
+  // create output array
+
+  fOutputArray = ExportArray(GetString("OutputArray", "tracks"));
+}
+
+//------------------------------------------------------------------------------
+
+void TrackCovariance::Finish()
+{
+  if(fItInputArray) delete fItInputArray;
+}
+
+//------------------------------------------------------------------------------
+
+void TrackCovariance::Process()
+{
+  Candidate *candidate, *mother;
+  Double_t energy;
+
+  fItInputArray->Reset();
+  while((candidate = static_cast<Candidate *>(fItInputArray->Next())))
+  {
+    const TLorentzVector &candidatePosition = candidate->Position;
+    const TLorentzVector &candidateMomentum = candidate->Momentum;
+
+    energy = candidateMomentum.E();
+
+    ObsTrk track(candidatePosition.Vect(), candidateMomentum.Vect(), candidate->Charge, fBz, fCovariance);
+
+    mother = candidate;
+    candidate = static_cast<Candidate *>(candidate->Clone());
+    candidate->Momentum.SetVect(track.GetObsP());
+    candidate->Momentum.SetE(energy);
+
+    candidate->AddCandidate(mother);
+
+    fOutputArray->Add(candidate);
+  }
+}
+
+//------------------------------------------------------------------------------
Index: modules/TrackCovariance.h
===================================================================
--- modules/TrackCovariance.h	(revision ff9fb2d9367d146a8325d4b77b9a9a4c0bd2ada8)
+++ modules/TrackCovariance.h	(revision ff9fb2d9367d146a8325d4b77b9a9a4c0bd2ada8)
@@ -0,0 +1,63 @@
+/*
+ *  Delphes: a framework for fast simulation of a generic collider experiment
+ *  Copyright (C) 2020  Universite catholique de Louvain (UCLouvain), Belgium
+ *
+ *  This program is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef TrackCovariance_h
+#define TrackCovariance_h
+
+/** \class TrackCovariance
+ *
+ *  Smears track parameters according to appropriate covariance matrix.
+ *
+ *  \author P. Demin - UCLouvain, Louvain-la-Neuve
+ *
+ */
+
+#include "classes/DelphesModule.h"
+
+class TIterator;
+class TObjArray;
+
+class SolGeom;
+class SolGridCov;
+
+class TrackCovariance: public DelphesModule
+{
+public:
+  TrackCovariance();
+  ~TrackCovariance();
+
+  void Init();
+  void Process();
+  void Finish();
+
+private:
+  Double_t fBz;
+
+  SolGeom *fGeometry;
+  SolGridCov *fCovariance;
+
+  TIterator *fItInputArray; //!
+
+  const TObjArray *fInputArray; //!
+
+  TObjArray *fOutputArray; //!
+
+  ClassDef(TrackCovariance, 1)
+};
+
+#endif
