diff --git a/math/genvector/inc/Math/GenVector/BitReproducible.h b/math/genvector/inc/Math/GenVector/BitReproducible.h
index f677ecd47b836..027206bdd642d 100644
--- a/math/genvector/inc/Math/GenVector/BitReproducible.h
+++ b/math/genvector/inc/Math/GenVector/BitReproducible.h
@@ -32,8 +32,8 @@ class BitReproducibleException  : public std::exception
 public:
   BitReproducibleException(const std::string & w) throw() : fMsg(w) {}
   ~BitReproducibleException() throw() {}
-  const char* what() const throw() { return fMsg.c_str(); }
-private:
+  const char *what() const throw() override { return fMsg.c_str(); }
+  private:
   std::string fMsg;
 };  // DoubConvException
 
diff --git a/math/genvector/inc/Math/GenVector/Cartesian2D.h b/math/genvector/inc/Math/GenVector/Cartesian2D.h
index ccccb9e54db53..adc6e981fb424 100644
--- a/math/genvector/inc/Math/GenVector/Cartesian2D.h
+++ b/math/genvector/inc/Math/GenVector/Cartesian2D.h
@@ -89,8 +89,8 @@ public :
    Scalar X()     const { return fX;}
    Scalar Y()     const { return fY;}
    Scalar Mag2()  const { return fX*fX + fY*fY; }
-   Scalar R()     const { return std::sqrt( Mag2());}
-   Scalar Phi()   const { return (fX==0 && fY==0) ? 0.0 : atan2(fY,fX);}
+   Scalar R() const { return std::sqrt(Mag2()); }
+   Scalar Phi() const { return (fX == Scalar(0) && fY == Scalar(0)) ? Scalar(0) : std::atan2(fY, fX); }
 
    /**
        set the x coordinate value keeping y constant
@@ -124,9 +124,9 @@ public :
        rotate by an angle
     */
    void Rotate(Scalar angle) {
-      Scalar s = std::sin(angle);
-      Scalar c = std::cos(angle);
-      SetCoordinates( c*fX - s*fY, s*fX + c * fY );
+      const Scalar s = std::sin(angle);
+      const Scalar c = std::cos(angle);
+      SetCoordinates(c * fX - s * fY, s * fX + c * fY);
    }
 
    /**
@@ -161,7 +161,7 @@ public :
    template <class T2>
    explicit Cartesian2D( const Polar2D<T2> & v )
    {
-      Scalar r = v.R();    // re-using this instead of calling v.X() and v.Y()
+      const Scalar r = v.R(); // re-using this instead of calling v.X() and v.Y()
       // is the speed improvement
       fX = r * std::cos(v.Phi());
       fY = r * std::sin(v.Phi());
@@ -174,9 +174,9 @@ public :
    template <class T2>
    Cartesian2D & operator = (const Polar2D<T2> & v)
    {
-      Scalar r = v.R();
-      fX = r * std::cos(v.Phi());
-      fY = r * std::sin(v.Phi());
+      const Scalar r = v.R();
+      fX             = r * std::cos(v.Phi());
+      fY             = r * std::sin(v.Phi());
       return *this;
    }
 
diff --git a/math/genvector/inc/Math/GenVector/Cartesian3D.h b/math/genvector/inc/Math/GenVector/Cartesian3D.h
index dfd87ba911af7..0e06d8e3b4484 100644
--- a/math/genvector/inc/Math/GenVector/Cartesian3D.h
+++ b/math/genvector/inc/Math/GenVector/Cartesian3D.h
@@ -24,11 +24,10 @@
 #include "Math/Math.h"
 
 #include <limits>
-
+#include <cmath>
 
 #include "Math/GenVector/eta.h"
 
-
 namespace ROOT {
 
 namespace Math {
@@ -84,7 +83,6 @@ public :
       return *this;
    }
 
-
    /**
       Set internal data based on an array of 3 Scalar numbers
    */
@@ -111,15 +109,17 @@ public :
    Scalar Z()     const { return fZ;}
    Scalar Mag2()  const { return fX*fX + fY*fY + fZ*fZ;}
    Scalar Perp2() const { return fX*fX + fY*fY ;}
-   Scalar Rho()   const { return std::sqrt( Perp2());}
-   Scalar R()     const { return std::sqrt( Mag2());}
-   Scalar Theta() const { return (fX==0 && fY==0 && fZ==0) ?
-                             0 : atan2(Rho(),Z());}
-   Scalar Phi()   const { return (fX==0 && fY==0) ? 0 : atan2(fY,fX);}
+   Scalar Rho() const { return std::sqrt(Perp2()); }
+   Scalar R() const { return std::sqrt(Mag2()); }
+   Scalar Theta() const
+   {
+      return (fX == Scalar(0) && fY == Scalar(0) && fZ == Scalar(0)) ? Scalar(0) : std::atan2(Rho(), Z());
+   }
+   Scalar Phi() const { return (fX == Scalar(0) && fY == Scalar(0)) ? Scalar(0) : std::atan2(fY, fX); }
 
    // pseudorapidity
    Scalar Eta() const {
-      return Impl::Eta_FromRhoZ ( Rho(),fZ);
+      return Impl::Eta_FromRhoZ( Rho(), fZ );
    }
 
    /**
@@ -149,7 +149,12 @@ public :
    /**
       scale the vector by a scalar quantity a
    */
-   void Scale(Scalar a) { fX *= a; fY *= a;  fZ *= a; }
+   void Scale(Scalar a)
+   {
+      fX *= a;
+      fY *= a;
+      fZ *= a;
+   }
 
    /**
       negate the vector
@@ -190,7 +195,8 @@ public :
    template <class T2>
    explicit Cartesian3D( const Polar3D<T2> & v ) : fZ (v.Z())
    {
-      T rho = v.Rho(); // re-using this instead of calling v.X() and v.Y()
+      const T rho = v.Rho();
+      // re-using this instead of calling v.X() and v.Y()
       // is the speed improvement
       fX = rho * std::cos(v.Phi());
       fY = rho * std::sin(v.Phi());
@@ -203,9 +209,9 @@ public :
    template <class T2>
    Cartesian3D & operator = (const Polar3D<T2> & v)
    {
-      T rho = v.Rho();
-      fX = rho * std::cos(v.Phi());
-      fY = rho * std::sin(v.Phi());
+      const T rho = v.Rho();
+      fX          = rho * std::cos(v.Phi());
+      fY          = rho * std::sin(v.Phi());
       fZ = v.Z();
       return *this;
    }
diff --git a/math/genvector/inc/Math/GenVector/Cylindrical3D.h b/math/genvector/inc/Math/GenVector/Cylindrical3D.h
index 1630c02782d4b..91702b3fc7d4a 100644
--- a/math/genvector/inc/Math/GenVector/Cylindrical3D.h
+++ b/math/genvector/inc/Math/GenVector/Cylindrical3D.h
@@ -21,7 +21,7 @@
 #include "Math/GenVector/eta.h"
 
 #include <limits>
-
+#include <cmath>
 
 namespace ROOT {
 
@@ -106,11 +106,10 @@ public :
    {rho=fRho; zz=fZ; phi=fPhi;}
 
 private:
-   inline static Scalar pi() { return M_PI; }
-   inline void Restrict() {
-      if ( fPhi <= -pi() || fPhi > pi() )
-         fPhi = fPhi - std::floor( fPhi/(2*pi()) +.5 ) * 2*pi();
-      return;
+   inline static Scalar pi() { return Scalar(M_PI); }
+   inline void          Restrict()
+   {
+      if (fPhi <= -pi() || fPhi > pi()) fPhi = fPhi - std::floor(fPhi / (2 * pi()) + .5) * 2 * pi();
    }
 public:
 
@@ -120,13 +119,13 @@ public :
    Scalar Z()     const { return fZ;   }
    Scalar Phi()   const { return fPhi; }
 
-   Scalar X()     const { return fRho*std::cos(fPhi); }
-   Scalar Y()     const { return fRho*std::sin(fPhi); }
+   Scalar X() const { return fRho * std::cos(fPhi); }
+   Scalar Y() const { return fRho * std::sin(fPhi); }
 
    Scalar Mag2()  const { return fRho*fRho + fZ*fZ;   }
-   Scalar R()     const { return std::sqrt( Mag2());  }
+   Scalar R() const { return std::sqrt(Mag2()); }
    Scalar Perp2() const { return fRho*fRho;           }
-   Scalar Theta() const { return (fRho==0 && fZ==0 ) ? 0 : atan2(fRho,fZ); }
+   Scalar Theta() const { return (fRho == Scalar(0) && fZ == Scalar(0)) ? Scalar(0) : std::atan2(fRho, fZ); }
 
    // pseudorapidity - use same implementation as in Cartesian3D
    Scalar Eta() const {
diff --git a/math/genvector/inc/Math/GenVector/CylindricalEta3D.h b/math/genvector/inc/Math/GenVector/CylindricalEta3D.h
index 0b4b0e20739ab..dfa958031f0ca 100644
--- a/math/genvector/inc/Math/GenVector/CylindricalEta3D.h
+++ b/math/genvector/inc/Math/GenVector/CylindricalEta3D.h
@@ -26,6 +26,7 @@
 
 
 #include <limits>
+#include <cmath>
 
 #include "Math/Math.h"
 
@@ -69,12 +70,12 @@ public :
   explicit CylindricalEta3D( const CoordSystem & v ) :
      fRho(v.Rho() ),  fEta(v.Eta() ),  fPhi(v.Phi() )
   {
-    static Scalar bigEta =
-         -.3f *std::log(std::numeric_limits<Scalar>::epsilon());
-    if ( std::fabs(fEta) > bigEta ) {
-       fRho *= v.Z()/Z(); // This gives a small absolute adjustment in rho,
-       // which, for large eta, results in a significant
-       // improvement in the faithfullness of reproducing z.
+     static Scalar bigEta = Scalar(-0.3) * log(std::numeric_limits<Scalar>::epsilon());
+     if (std::fabs(fEta) > bigEta) {
+        // This gives a small absolute adjustment in rho,
+        // which, for large eta, results in a significant
+        // improvement in the faithfullness of reproducing z.
+        fRho *= v.Z() / Z();
     }
   }
 
@@ -124,8 +125,7 @@ public :
 private:
    inline static Scalar pi() { return M_PI; }
    inline void Restrict() {
-      if ( fPhi <= -pi() || fPhi > pi() )
-         fPhi = fPhi - std::floor( fPhi/(2*pi()) +.5 ) * 2*pi();
+      if (fPhi <= -pi() || fPhi > pi()) fPhi = fPhi - std::floor(fPhi / (2 * pi()) + .5) * 2 * pi();
       return;
    }
 public:
@@ -135,20 +135,24 @@ public :
    T Rho()   const { return fRho; }
    T Eta()   const { return fEta; }
    T Phi()   const { return fPhi; }
-   T X()     const { return fRho*std::cos(fPhi); }
-   T Y()     const { return fRho*std::sin(fPhi); }
-   T Z()     const { return fRho >  0 ? fRho*std::sinh(fEta) :
-      fEta == 0 ? 0                    :
-      fEta >  0 ? fEta - etaMax<T>()   :
-      fEta + etaMax<T>()   ; }
-   T R()     const { return fRho > 0          ? fRho*std::cosh(fEta) :
-      fEta >  etaMax<T>() ?  fEta - etaMax<T>()   :
-      fEta < -etaMax<T>() ? -fEta - etaMax<T>()   :
-      0     ; }
-   T Mag2()  const { return R()*R();              }
+   T X() const { return fRho * std::cos(fPhi); }
+   T Y() const { return fRho * std::sin(fPhi); }
+   T Z() const
+   {
+      return fRho > 0 ? fRho * std::sinh(fEta) : fEta == 0 ? 0 : fEta > 0 ? fEta - etaMax<T>() : fEta + etaMax<T>();
+   }
+   T R() const
+   {
+      return fRho > 0 ? fRho * std::cosh(fEta)
+                      : fEta > etaMax<T>() ? fEta - etaMax<T>() : fEta < -etaMax<T>() ? -fEta - etaMax<T>() : 0;
+   }
+   T Mag2() const
+   {
+      const Scalar r = R();
+      return r * r;
+   }
    T Perp2() const { return fRho*fRho;            }
-   T Theta() const { return  fRho >  0 ? 2* std::atan( std::exp( - fEta ) ) :
-         (fEta >= 0 ? 0 : pi() );  }
+   T Theta() const { return fRho > 0 ? 2 * std::atan(exp(-fEta)) : (fEta >= 0 ? 0 : pi()); }
 
    // setters (only for data members)
 
diff --git a/math/genvector/inc/Math/GenVector/DisplacementVector3D.h b/math/genvector/inc/Math/GenVector/DisplacementVector3D.h
index 3fceba90daeb1..d197c757dbb1c 100644
--- a/math/genvector/inc/Math/GenVector/DisplacementVector3D.h
+++ b/math/genvector/inc/Math/GenVector/DisplacementVector3D.h
@@ -249,11 +249,13 @@ namespace ROOT {
        */
       template <class IT>
       void GetCoordinates( IT begin) const {
-         Scalar a,b,c = 0;
-         GetCoordinates (a,b,c);
+         Scalar a = Scalar(0);
+         Scalar b = Scalar(0);
+         Scalar c = Scalar(0);
+         GetCoordinates(a, b, c);
          *begin++ = a;
          *begin++ = b;
-         *begin = c;
+         *begin   = c;
       }
 
       /**
@@ -262,8 +264,8 @@ namespace ROOT {
          then (x, y, z) are converted to that form)
        */
       DisplacementVector3D<CoordSystem, Tag>& SetXYZ (Scalar a, Scalar b, Scalar c) {
-            fCoordinates.SetXYZ(a,b,c);
-            return *this;
+         fCoordinates.SetXYZ(a, b, c);
+         return *this;
       }
 
       // ------------------- Equality -----------------
@@ -333,11 +335,24 @@ namespace ROOT {
       Scalar Perp2() const { return fCoordinates.Perp2();}
 
       /**
-         return unit vector parallel to this
+         return unit vector parallel to this (scalar)
+      */
+      template <typename SCALAR = Scalar, typename std::enable_if<std::is_arithmetic<SCALAR>::value>::type * = nullptr>
+      DisplacementVector3D Unit() const
+      {
+         const auto tot = R();
+         return tot == 0 ? *this : DisplacementVector3D(*this) / tot;
+      }
+
+      /**
+         return unit vector parallel to this (vector)
       */
-      DisplacementVector3D  Unit() const {
-        Scalar tot = R();
-        return tot == 0 ? *this : DisplacementVector3D(*this) / tot;
+      template <typename SCALAR = Scalar, typename std::enable_if<!std::is_arithmetic<SCALAR>::value>::type * = nullptr>
+      DisplacementVector3D Unit() const
+      {
+         SCALAR tot            = R();
+         tot(tot == SCALAR(0)) = SCALAR(1);
+         return DisplacementVector3D(*this) / tot;
       }
 
       // ------ Setting of individual elements present in coordinate system ------
@@ -619,36 +634,46 @@ namespace ROOT {
 
     // ------------- I/O to/from streams -------------
 
-    template< class char_t, class traits_t, class T, class U >
-      inline
-      std::basic_ostream<char_t,traits_t> &
-      operator << ( std::basic_ostream<char_t,traits_t> & os
-                  , DisplacementVector3D<T,U> const & v
-                  )
+    template <class char_t, class traits_t, class T, class U,
+              typename std::enable_if<std::is_arithmetic<typename DisplacementVector3D<T, U>::Scalar>::value>::type * =
+                 nullptr>
+    std::basic_ostream<char_t, traits_t> &operator<<(std::basic_ostream<char_t, traits_t> &os,
+                                                     DisplacementVector3D<T, U> const &v)
     {
-      if( !os )  return os;
-
-      typename T::Scalar a, b, c;
-      v.GetCoordinates(a, b, c);
-
-      if( detail::get_manip( os, detail::bitforbit ) )  {
-        detail::set_manip( os, detail::bitforbit, '\00' );
-        typedef GenVector_detail::BitReproducible BR;
-        BR::Output(os, a);
-        BR::Output(os, b);
-        BR::Output(os, c);
-      }
-      else  {
-        os << detail::get_manip( os, detail::open  ) << a
-           << detail::get_manip( os, detail::sep   ) << b
-           << detail::get_manip( os, detail::sep   ) << c
-           << detail::get_manip( os, detail::close );
+       if (os) {
+
+          typename T::Scalar a, b, c;
+          v.GetCoordinates(a, b, c);
+
+          if (detail::get_manip(os, detail::bitforbit)) {
+             detail::set_manip(os, detail::bitforbit, '\00');
+             typedef GenVector_detail::BitReproducible BR;
+             BR::Output(os, a);
+             BR::Output(os, b);
+             BR::Output(os, c);
+          } else {
+             os << detail::get_manip(os, detail::open) << a << detail::get_manip(os, detail::sep) << b
+                << detail::get_manip(os, detail::sep) << c << detail::get_manip(os, detail::close);
+          }
       }
-
       return os;
-
     }  // op<< <>()
 
+    template <class char_t, class traits_t, class T, class U,
+              typename std::enable_if<!std::is_arithmetic<typename DisplacementVector3D<T, U>::Scalar>::value>::type * =
+                 nullptr>
+    std::basic_ostream<char_t, traits_t> &operator<<(std::basic_ostream<char_t, traits_t> &os,
+                                                     DisplacementVector3D<T, U> const &v)
+    {
+       if (os) {
+          os << "{ ";
+          for (std::size_t i = 0; i < PositionVector3D<T, U>::Scalar::Size; ++i) {
+             os << "(" << v.x()[i] << "," << v.y()[i] << "," << v.z()[i] << ") ";
+          }
+          os << "}";
+       }
+       return os;
+    } // op<< <>()
 
     template< class char_t, class traits_t, class T, class U >
       inline
diff --git a/math/genvector/inc/Math/GenVector/LorentzVector.h b/math/genvector/inc/Math/GenVector/LorentzVector.h
index 026a3593c9133..932100493118e 100644
--- a/math/genvector/inc/Math/GenVector/LorentzVector.h
+++ b/math/genvector/inc/Math/GenVector/LorentzVector.h
@@ -24,7 +24,7 @@
 
 #include "Math/GenVector/GenVectorIO.h"
 
-
+#include <cmath>
 
 namespace ROOT {
 
@@ -488,9 +488,9 @@ namespace ROOT {
           // TODO - It would be good to check that E > Pz and use the Throw()
           //        mechanism or at least load a NAN if not.
           //        We should then move the code to a .cpp file.
-          Scalar ee = E();
-          Scalar ppz = Pz();
-          return .5f* std::log( (ee+ppz)/(ee-ppz) );
+          const Scalar ee  = E();
+          const Scalar ppz = Pz();
+          return Scalar(0.5) * std::log((ee + ppz) / (ee - ppz));
        }
 
        /**
@@ -499,9 +499,9 @@ namespace ROOT {
        Scalar ColinearRapidity() const {
           // TODO - It would be good to check that E > P and use the Throw()
           //        mechanism or at least load a NAN if not.
-          Scalar ee = E();
-          Scalar pp = P();
-          return .5f* std::log( (ee+pp)/(ee-pp) );
+          const Scalar ee = E();
+          const Scalar pp = P();
+          return Scalar(0.5) * std::log((ee + pp) / (ee - pp));
        }
 
        /**
@@ -597,8 +597,8 @@ namespace ROOT {
            Return Gamma scalar value
        */
        Scalar Gamma() const {
-          Scalar v2 = P2();
-          Scalar t2 = E()*E();
+          const Scalar v2 = P2();
+          const Scalar t2 = std::pow(E(), 2);
           if (E() == 0) {
              if ( P2() == 0) {
                 return 1;
@@ -614,7 +614,7 @@ namespace ROOT {
           else if ( t2 == v2 ) {
              GenVector::Throw ("LorentzVector::Gamma() - gamma computed for a lightlike LorentzVector. Infinite result");
           }
-          return 1./std::sqrt(1. - v2/t2 );
+          return Scalar(1) / std::sqrt(Scalar(1) - v2 / t2);
        } /* gamma */
 
 
diff --git a/math/genvector/inc/Math/GenVector/Plane3D.h b/math/genvector/inc/Math/GenVector/Plane3D.h
index a0eda70fd02ce..5871d7598b182 100644
--- a/math/genvector/inc/Math/GenVector/Plane3D.h
+++ b/math/genvector/inc/Math/GenVector/Plane3D.h
@@ -17,6 +17,8 @@
 #ifndef ROOT_Math_GenVector_Plane3D
 #define ROOT_Math_GenVector_Plane3D  1
 
+#include <type_traits>
+
 #include "Math/GenVector/DisplacementVector3D.h"
 #include "Math/GenVector/PositionVector3D.h"
 
@@ -26,228 +28,265 @@ namespace ROOT {
 
 namespace Math {
 
-
+namespace Impl {
 
 //_______________________________________________________________________________
-   /**
-      Class describing a geometrical plane in 3 dimensions.
-      A Plane3D is a 2 dimensional surface spanned by two linearly independent vectors.
-      The plane is described by the equation
-      \f$ a*x + b*y + c*z + d = 0 \f$ where (a,b,c) are the components of the
-      normal vector to the plane \f$ n = (a,b,c)  \f$ and \f$ d = - n \dot x \f$, where x is any point
-      belonging to plane.
-      More information on the mathematics describing a plane in 3D is available on
-      <A HREF=http://mathworld.wolfram.com/Plane.html>MathWord</A>.
-      The Plane3D class contains the 4 scalar values in double which represent the
-      four coefficients, fA, fB, fC, fD. fA, fB, fC are the normal components normalized to 1,
-      i.e. fA**2 + fB**2 + fC**2 = 1
-
-      @ingroup GenVector
-   */
-   class Plane3D {
+/**
+   Class describing a geometrical plane in 3 dimensions.
+   A Plane3D is a 2 dimensional surface spanned by two linearly independent vectors.
+   The plane is described by the equation
+   \f$ a*x + b*y + c*z + d = 0 \f$ where (a,b,c) are the components of the
+   normal vector to the plane \f$ n = (a,b,c)  \f$ and \f$ d = - n \dot x \f$, where x is any point
+   belonging to plane.
+   More information on the mathematics describing a plane in 3D is available on
+   <A HREF=http://mathworld.wolfram.com/Plane.html>MathWord</A>.
+   The Plane3D class contains the 4 scalar values in T which represent the
+   four coefficients, fA, fB, fC, fD. fA, fB, fC are the normal components normalized to 1,
+   i.e. fA**2 + fB**2 + fC**2 = 1
 
-   public:
+   @ingroup GenVector
+*/
 
-      // ------ ctors ------
+template <typename T = double>
+class Plane3D {
 
-      typedef double Scalar;
+public:
+   // ------ ctors ------
 
-      typedef  DisplacementVector3D<Cartesian3D<double>, DefaultCoordinateSystemTag > Vector;
-      typedef  PositionVector3D<Cartesian3D<double>, DefaultCoordinateSystemTag >     Point;
+   typedef T Scalar;
 
+   typedef DisplacementVector3D<Cartesian3D<T>, DefaultCoordinateSystemTag> Vector;
+   typedef PositionVector3D<Cartesian3D<T>, DefaultCoordinateSystemTag>     Point;
 
+   /**
+      default constructor create plane z = 0
+   */
+   Plane3D() : fA(0), fB(0), fC(1), fD(0) {}
 
-      /**
-         default constructor create plane z = 0
-      */
-      Plane3D ( ) : fA(0), fB(0), fC(1.), fD(0) { }
-
-      /**
-       generic constructors from the four scalar values describing the plane
-       according to the equation ax + by + cz + d = 0
-         \param a scalar value
-         \param b scalar value
-         \param c scalar value
-         \param d sxcalar value
-      */
-      Plane3D(const Scalar & a, const Scalar & b, const Scalar & c, const Scalar & d);
+   /**
+    generic constructors from the four scalar values describing the plane
+    according to the equation ax + by + cz + d = 0
+      \param a scalar value
+      \param b scalar value
+      \param c scalar value
+      \param d sxcalar value
+   */
+   Plane3D(const Scalar &a, const Scalar &b, const Scalar &c, const Scalar &d) : fA(a), fB(b), fC(c), fD(d)
+   {
+      // renormalize a,b,c to unit
+      Normalize();
+   }
 
-      /**
-       constructor a Plane3D from a normal vector and a point coplanar to the plane
-       \param n normal expressed as a ROOT::Math::DisplacementVector3D<Cartesian3D<double> >
-       \param p point  expressed as a  ROOT::Math::PositionVector3D<Cartesian3D<double> >
-      */
-      Plane3D(const Vector & n, const Point & p )
-      {
-         BuildFromVecAndPoint( n, p );
-      }
+   /**
+    constructor a Plane3D from a normal vector and a point coplanar to the plane
+    \param n normal expressed as a ROOT::Math::DisplacementVector3D<Cartesian3D<T> >
+    \param p point  expressed as a  ROOT::Math::PositionVector3D<Cartesian3D<T> >
+   */
+   Plane3D(const Vector &n, const Point &p) { BuildFromVecAndPoint(n, p); }
 
+   /**
+    Construct from a generic DisplacementVector3D (normal vector) and PositionVector3D (point coplanar to
+    the plane)
+    \param n normal expressed as a generic ROOT::Math::DisplacementVector3D
+    \param p point  expressed as a generic ROOT::Math::PositionVector3D
+   */
+   template <class T1, class T2, class U>
+   Plane3D(const DisplacementVector3D<T1, U> &n, const PositionVector3D<T2, U> &p)
+   {
+      BuildFromVecAndPoint(Vector(n), Point(p));
+   }
 
-      /**
-       Construct from a generic DisplacementVector3D (normal vector) and PositionVector3D (point coplanar to
-       the plane)
-       \param n normal expressed as a generic ROOT::Math::DisplacementVector3D
-       \param p point  expressed as a generic ROOT::Math::PositionVector3D
-      */
-      template<class T1, class T2, class U>
-      Plane3D( const  DisplacementVector3D<T1,U> & n, const  PositionVector3D<T2,U> & p)
-      {
-         BuildFromVecAndPoint( Vector(n), Point(p) );
-      }
+   /**
+    constructor from three Cartesian point belonging to the plane
+    \param p1 point1  expressed as a generic ROOT::Math::PositionVector3D
+    \param p2 point2  expressed as a generic ROOT::Math::PositionVector3D
+    \param p3 point3  expressed as a generic ROOT::Math::PositionVector3D
+   */
+   Plane3D(const Point &p1, const Point &p2, const Point &p3) { BuildFrom3Points(p1, p2, p3); }
 
-      /**
-       constructor from three Cartesian point belonging to the plane
-       \param p1 point1  expressed as a generic ROOT::Math::PositionVector3D
-       \param p2 point2  expressed as a generic ROOT::Math::PositionVector3D
-       \param p3 point3  expressed as a generic ROOT::Math::PositionVector3D
-      */
-      Plane3D(const Point & p1, const Point & p2, const Point & p3  ) {
-         BuildFrom3Points(p1,p2,p3);
-      }
+   /**
+    constructor from three generic point belonging to the plane
+    \param p1 point1 expressed as  ROOT::Math::DisplacementVector3D<Cartesian3D<T> >
+    \param p2 point2 expressed as  ROOT::Math::DisplacementVector3D<Cartesian3D<T> >
+    \param p3 point3 expressed as  ROOT::Math::DisplacementVector3D<Cartesian3D<T> >
+   */
+   template <class T1, class T2, class T3, class U>
+   Plane3D(const PositionVector3D<T1, U> &p1, const PositionVector3D<T2, U> &p2, const PositionVector3D<T3, U> &p3)
+   {
+      BuildFrom3Points(Point(p1.X(), p1.Y(), p1.Z()), Point(p2.X(), p2.Y(), p2.Z()), Point(p3.X(), p3.Y(), p3.Z()));
+   }
 
-      /**
-       constructor from three generic point belonging to the plane
-       \param p1 point1 expressed as  ROOT::Math::DisplacementVector3D<Cartesian3D<double> >
-       \param p2 point2 expressed as  ROOT::Math::DisplacementVector3D<Cartesian3D<double> >
-       \param p3 point3 expressed as  ROOT::Math::DisplacementVector3D<Cartesian3D<double> >
-      */
-      template <class T1, class T2, class T3, class U>
-      Plane3D(const  PositionVector3D<T1,U> & p1, const  PositionVector3D<T2,U> & p2, const  PositionVector3D<T3,U> & p3  )
-      {
-         BuildFrom3Points( Point(p1.X(), p1.Y(), p1.Z()),
-                           Point(p2.X(), p2.Y(), p2.Z()),
-                           Point(p3.X(), p3.Y(), p3.Z()) );
-      }
+   // compiler-generated copy ctor and dtor are fine.
 
+   // ------ assignment ------
 
+   /**
+      Assignment operator from other Plane3D class
+   */
+   Plane3D &operator=(const Plane3D &plane)
+   {
+      fA = plane.fA;
+      fB = plane.fB;
+      fC = plane.fC;
+      fD = plane.fD;
+      return *this;
+   }
 
-      // compiler-generated copy ctor and dtor are fine.
+   /**
+      Return the a coefficient of the plane equation \f$ a*x + b*y + c*z + d = 0 \f$. It is also the
+      x-component of the vector perpendicular to the plane.
+   */
+   Scalar A() const { return fA; }
 
-      // ------ assignment ------
+   /**
+      Return the b coefficient of the plane equation \f$ a*x + b*y + c*z + d = 0 \f$. It is also the
+      y-component of the vector perpendicular to the plane
+   */
+   Scalar B() const { return fB; }
 
-      /**
-         Assignment operator from other Plane3D class
-      */
-      Plane3D & operator= ( const Plane3D & plane) {
-         fA = plane.fA;
-         fB = plane.fB;
-         fC = plane.fC;
-         fD = plane.fD;
-         return *this;
-      }
+   /**
+      Return the c coefficient of the plane equation \f$ a*x + b*y + c*z + d = 0 \f$. It is also the
+      z-component of the vector perpendicular to the plane
+   */
+   Scalar C() const { return fC; }
 
-      /**
-         Return the a coefficient of the plane equation \f$ a*x + b*y + c*z + d = 0 \f$. It is also the
-         x-component of the vector perpendicular to the plane.
-      */
-      Scalar A() const { return fA; }
-
-      /**
-         Return the b coefficient of the plane equation \f$ a*x + b*y + c*z + d = 0 \f$. It is also the
-         y-component of the vector perpendicular to the plane
-      */
-      Scalar B() const { return fB; }
-
-      /**
-         Return the c coefficient of the plane equation \f$ a*x + b*y + c*z + d = 0 \f$. It is also the
-         z-component of the vector perpendicular to the plane
-      */
-      Scalar C() const { return fC; }
-
-      /**
-         Return the d coefficient of the plane equation \f$ a*x + b*y + c*z + d = 0 \f$. It is also
-         the distance from the origin (HesseDistance)
-      */
-      Scalar D() const { return fD; }
-
-      /**
-         Return normal vector to the plane as Cartesian DisplacementVector
-      */
-      Vector Normal() const {
-         return Vector(fA, fB, fC);
-      }
+   /**
+      Return the d coefficient of the plane equation \f$ a*x + b*y + c*z + d = 0 \f$. It is also
+      the distance from the origin (HesseDistance)
+   */
+   Scalar D() const { return fD; }
 
-      /**
-       Return the Hesse Distance (distance from the origin) of the plane or
-       the d coefficient expressed in normalize form
-      */
-      Scalar HesseDistance() const {
-         return fD;
-      }
+   /**
+      Return normal vector to the plane as Cartesian DisplacementVector
+   */
+   Vector Normal() const { return Vector(fA, fB, fC); }
 
+   /**
+    Return the Hesse Distance (distance from the origin) of the plane or
+    the d coefficient expressed in normalize form
+   */
+   Scalar HesseDistance() const { return fD; }
 
-      /**
-       Return the signed distance to a Point.
-       The distance is signed positive if the Point is in the same side of the
-       normal vector to the plane.
-       \param p Point expressed in Cartesian Coordinates
-       */
-      Scalar Distance(const Point & p) const;
-
-      /**
-       Return the distance to a Point described with generic coordinates
-       \param p Point expressed as generic ROOT::Math::PositionVector3D
-       */
-      template <class T, class U>
-      Scalar Distance(const PositionVector3D<T,U> & p) const {
-         return Distance( Point(p.X(), p.Y(), p.Z() ) );
-      }
+   /**
+    Return the signed distance to a Point.
+    The distance is signed positive if the Point is in the same side of the
+    normal vector to the plane.
+    \param p Point expressed in Cartesian Coordinates
+    */
+   Scalar Distance(const Point &p) const { return fA * p.X() + fB * p.Y() + fC * p.Z() + fD; }
 
-      /**
-       Return the projection of a Cartesian point to a plane
-       \param p Point expressed as PositionVector3D<Cartesian3D<double> >
-       */
-      Point ProjectOntoPlane(const Point & p) const;
-
-      /**
-       Return the projection of a point to a plane
-       \param p Point expressed as generic ROOT::Math::PositionVector3D
-       */
-      template <class T, class U>
-      PositionVector3D<T,U> ProjectOntoPlane(const PositionVector3D<T,U> & p) const {
-         Point pxyz = ProjectOntoPlane(Point(p.X(), p.Y(), p.Z() ) );
-         PositionVector3D<T,U> p2;
-         p2.SetXYZ( pxyz.X(), pxyz.Y(), pxyz.Z() );
-         return p2;
-      }
+   /**
+    Return the distance to a Point described with generic coordinates
+    \param p Point expressed as generic ROOT::Math::PositionVector3D
+    */
+   template <class T1, class U>
+   Scalar Distance(const PositionVector3D<T1, U> &p) const
+   {
+      return Distance(Point(p.X(), p.Y(), p.Z()));
+   }
 
+   /**
+    Return the projection of a Cartesian point to a plane
+    \param p Point expressed as PositionVector3D<Cartesian3D<T> >
+    */
+   Point ProjectOntoPlane(const Point &p) const
+   {
+      const Scalar d = Distance(p);
+      return XYZPoint(p.X() - fA * d, p.Y() - fB * d, p.Z() - fC * d);
+   }
 
+   /**
+    Return the projection of a point to a plane
+    \param p Point expressed as generic ROOT::Math::PositionVector3D
+    */
+   template <class T1, class U>
+   PositionVector3D<T1, U> ProjectOntoPlane(const PositionVector3D<T1, U> &p) const
+   {
+      const Point pxyz = ProjectOntoPlane(Point(p.X(), p.Y(), p.Z()));
+      return PositionVector3D<T, U>(pxyz.X(), pxyz.Y(), pxyz.Z());
+   }
+
+   // ------------------- Equality -----------------
 
-      // ------------------- Equality -----------------
+   /**
+      Exact equality
+   */
+   bool operator==(const Plane3D &rhs) const { return (fA == rhs.fA && fB == rhs.fB && fC == rhs.fC && fD == rhs.fD); }
+   bool operator!=(const Plane3D &rhs) const { return !(operator==(rhs)); }
 
-      /**
-         Exact equality
-      */
-      bool operator==(const Plane3D & rhs) const {
-         return  fA  == rhs.fA &&  fB == rhs.fB  &&  fC == rhs.fC && fD == rhs.fD;
-      }
-      bool operator!= (const Plane3D & rhs) const {
-         return !(operator==(rhs));
+protected:
+   /**
+      Normalize the normal (a,b,c) plane components
+   */
+   template <typename SCALAR = T, typename std::enable_if<std::is_arithmetic<SCALAR>::value>::type * = nullptr>
+   void Normalize()
+   {
+      // normalize the plane
+      const SCALAR s = std::sqrt(fA * fA + fB * fB + fC * fC);
+      // what to do if s = 0 ?
+      if (s == SCALAR(0)) {
+         fD = SCALAR(0);
+      } else {
+         const SCALAR w = Scalar(1) / s;
+         fA *= w;
+         fB *= w;
+         fC *= w;
+         fD *= w;
       }
+   }
 
-   protected:
-
-      /**
-         Normalize the normal (a,b,c) plane components
-      */
-      void Normalize();
-
-
-   private:
-
-      // internal method to construct class from a vector and a point
-      void BuildFromVecAndPoint(const Vector & n, const Point & p);
-      // internal method to construct class from 3 points
-      void BuildFrom3Points(const Point & p1, const Point & p2, const Point & p3);
-
-      // plane data members the four scalar which  satisfies fA*x + fB*y + fC*z + fD = 0
-      // for every point (x,y,z) belonging to the plane.
-      // fA**2 + fB**2 + fC** =1 plane is stored in normalized form
-      Scalar  fA;
-      Scalar  fB;
-      Scalar  fC;
-      Scalar  fD;
-
+   /**
+     Normalize the normal (a,b,c) plane components
+   */
+   template <typename SCALAR = T, typename std::enable_if<!std::is_arithmetic<SCALAR>::value>::type * = nullptr>
+   void Normalize()
+   {
+      // normalize the plane
+      SCALAR s = sqrt(fA * fA + fB * fB + fC * fC);
+      // what to do if s = 0 ?
+      const auto m = (s == SCALAR(0));
+      // set zero entries to 1 in the vector to avoid /0 later on
+      s(m)           = SCALAR(1);
+      fD(m)          = SCALAR(0);
+      const SCALAR w = SCALAR(1) / s;
+      fA *= w;
+      fB *= w;
+      fC *= w;
+      fD *= w;
+   }
+
+private:
+   // internal method to construct class from a vector and a point
+   void BuildFromVecAndPoint(const Vector &n, const Point &p)
+   {
+      // build from a normal vector and a point
+      fA = n.X();
+      fB = n.Y();
+      fC = n.Z();
+      fD = -n.Dot(p);
+      Normalize();
+   }
+
+   // internal method to construct class from 3 points
+   void BuildFrom3Points(const Point &p1, const Point &p2, const Point &p3)
+   {
+      // plane from thre points
+      // normal is (x3-x1) cross (x2 -x1)
+      const Vector n = (p2 - p1).Cross(p3 - p1);
+      fA             = n.X();
+      fB             = n.Y();
+      fC             = n.Z();
+      fD             = -n.Dot(p1);
+      Normalize();
+   }
+
+   // plane data members the four scalar which  satisfies fA*x + fB*y + fC*z + fD = 0
+   // for every point (x,y,z) belonging to the plane.
+   // fA**2 + fB**2 + fC** =1 plane is stored in normalized form
+   Scalar fA;
+   Scalar fB;
+   Scalar fC;
+   Scalar fD;
 
    };  // Plane3D<>
 
@@ -255,9 +294,19 @@ namespace Math {
       Stream Output and Input
    */
    // TODO - I/O should be put in the manipulator form
-
-   std::ostream & operator<< (std::ostream & os, const Plane3D & p);
-
+   template <typename T>
+   std::ostream &operator<<(std::ostream &os, const Plane3D<T> &p)
+   {
+      os << "\n"
+         << p.Normal().X() << "  " << p.Normal().Y() << "  " << p.Normal().Z() << "  " << p.HesseDistance() << "\n";
+      return os;
+   }
+
+   } // end namespace Impl
+
+   // typedefs for double and float versions
+   typedef Impl::Plane3D<double> Plane3D;
+   typedef Impl::Plane3D<float>  Plane3DF;
 
 } // end namespace Math
 
diff --git a/math/genvector/inc/Math/GenVector/Polar2D.h b/math/genvector/inc/Math/GenVector/Polar2D.h
index 895778d3a55ab..64eb7e71206f2 100644
--- a/math/genvector/inc/Math/GenVector/Polar2D.h
+++ b/math/genvector/inc/Math/GenVector/Polar2D.h
@@ -97,8 +97,8 @@ public :
 
    Scalar R()     const { return fR;}
    Scalar Phi()   const { return fPhi; }
-   Scalar X()     const { return fR*std::cos(fPhi);}
-   Scalar Y()     const { return fR*std::sin(fPhi);}
+   Scalar X() const { return fR * std::cos(fPhi); }
+   Scalar Y() const { return fR * std::sin(fPhi); }
    Scalar Mag2()  const { return fR*fR;}
 
 
@@ -134,9 +134,7 @@ public :
       restrict abgle hi to be between -PI and PI
     */
    inline void Restrict() {
-      if ( fPhi <= -pi() || fPhi > pi() )
-         fPhi = fPhi - std::floor( fPhi/(2*pi()) +.5 ) * 2*pi();
-      return;
+      if (fPhi <= -pi() || fPhi > pi()) fPhi = fPhi - std::floor(fPhi / (2 * pi()) + .5) * 2 * pi();
    }
 
 public:
diff --git a/math/genvector/inc/Math/GenVector/Polar3D.h b/math/genvector/inc/Math/GenVector/Polar3D.h
index 3fa1e66e15922..a92c761ba9d1c 100644
--- a/math/genvector/inc/Math/GenVector/Polar3D.h
+++ b/math/genvector/inc/Math/GenVector/Polar3D.h
@@ -23,6 +23,7 @@
 
 #include "Math/GenVector/eta.h"
 
+#include <cmath>
 
 namespace ROOT {
 
@@ -109,12 +110,12 @@ public :
    Scalar R()     const { return fR;}
    Scalar Phi()   const { return fPhi; }
    Scalar Theta() const { return fTheta; }
-   Scalar Rho()   const { return fR*std::sin(fTheta); }
-   Scalar X()     const { return Rho()*std::cos(fPhi);}
-   Scalar Y()     const { return Rho()*std::sin(fPhi);}
-   Scalar Z()     const { return fR*std::cos(fTheta); }
+   Scalar Rho() const { return fR * std::sin(fTheta); }
+   Scalar X() const { return Rho() * std::cos(fPhi); }
+   Scalar Y() const { return Rho() * std::sin(fPhi); }
+   Scalar Z() const { return fR * std::cos(fTheta); }
    Scalar Mag2()  const { return fR*fR;}
-   Scalar Perp2() const { return Rho()*Rho(); }
+   Scalar Perp2() const { return std::pow(Rho(), 2); }
 
    // pseudorapidity
    Scalar Eta() const
@@ -156,9 +157,7 @@ public :
 private:
    inline static Scalar pi()  { return M_PI; }
    inline void Restrict() {
-      if ( fPhi <= -pi() || fPhi > pi() )
-         fPhi = fPhi - std::floor( fPhi/(2*pi()) +.5 ) * 2*pi();
-      return;
+      if (fPhi <= -pi() || fPhi > pi()) fPhi = fPhi - std::floor(fPhi / (2 * pi()) + .5) * 2 * pi();
    }
 
 public:
@@ -208,8 +207,8 @@ public :
 
    // The following make this coordinate system look enough like a CLHEP
    // vector that an assignment member template can work with either
-   T x() const { return X();}
-   T y() const { return Y();}
+   T x() const { return X(); }
+   T y() const { return Y(); }
    T z() const { return Z(); }
 
    // ============= Specializations for improved speed ==================
diff --git a/math/genvector/inc/Math/GenVector/PositionVector3D.h b/math/genvector/inc/Math/GenVector/PositionVector3D.h
index e27d9dce48887..9a1aa5cb08adc 100644
--- a/math/genvector/inc/Math/GenVector/PositionVector3D.h
+++ b/math/genvector/inc/Math/GenVector/PositionVector3D.h
@@ -227,8 +227,10 @@ namespace ROOT {
        */
       template <class IT>
       void GetCoordinates( IT begin ) const {
-         Scalar a,b,c = 0;
-         GetCoordinates (a,b,c);
+         Scalar a = Scalar(0);
+         Scalar b = Scalar(0);
+         Scalar c = Scalar(0);
+         GetCoordinates(a, b, c);
          *begin++ = a;
          *begin++ = b;
          *begin   = c;
@@ -575,36 +577,48 @@ namespace ROOT {
 
     // ------------- I/O to/from streams -------------
 
-    template< class char_t, class traits_t, class T, class U >
-      inline
-      std::basic_ostream<char_t,traits_t> &
-      operator << ( std::basic_ostream<char_t,traits_t> & os
-                  , PositionVector3D<T,U> const & v
-                  )
+    template <
+       class char_t, class traits_t, class T, class U,
+       typename std::enable_if<std::is_arithmetic<typename PositionVector3D<T, U>::Scalar>::value>::type * = nullptr>
+    std::basic_ostream<char_t, traits_t> &operator<<(std::basic_ostream<char_t, traits_t> &os,
+                                                     PositionVector3D<T, U> const &v)
     {
-      if( !os )  return os;
-
-      typename T::Scalar a, b, c;
-      v.GetCoordinates(a, b, c);
-
-      if( detail::get_manip( os, detail::bitforbit ) )  {
-        detail::set_manip( os, detail::bitforbit, '\00' );
-        typedef GenVector_detail::BitReproducible BR;
-        BR::Output(os, a);
-        BR::Output(os, b);
-        BR::Output(os, c);
-      }
-      else  {
-        os << detail::get_manip( os, detail::open  ) << a
-           << detail::get_manip( os, detail::sep   ) << b
-           << detail::get_manip( os, detail::sep   ) << c
-           << detail::get_manip( os, detail::close );
+       if (os) {
+
+          typename T::Scalar a = 0;
+          typename T::Scalar b = 0;
+          typename T::Scalar c = 0;
+          v.GetCoordinates(a, b, c);
+
+          if (detail::get_manip(os, detail::bitforbit)) {
+             detail::set_manip(os, detail::bitforbit, '\00');
+             typedef GenVector_detail::BitReproducible BR;
+             BR::Output(os, a);
+             BR::Output(os, b);
+             BR::Output(os, c);
+          } else {
+             os << detail::get_manip(os, detail::open) << a << detail::get_manip(os, detail::sep) << b
+                << detail::get_manip(os, detail::sep) << c << detail::get_manip(os, detail::close);
+          }
       }
-
       return os;
-
     }  // op<< <>()
 
+    template <
+       class char_t, class traits_t, class T, class U,
+       typename std::enable_if<!std::is_arithmetic<typename PositionVector3D<T, U>::Scalar>::value>::type * = nullptr>
+    std::basic_ostream<char_t, traits_t> &operator<<(std::basic_ostream<char_t, traits_t> &os,
+                                                     PositionVector3D<T, U> const &v)
+    {
+       if (os) {
+          os << "{ ";
+          for (std::size_t i = 0; i < PositionVector3D<T, U>::Scalar::Size; ++i) {
+             os << "(" << v.x()[i] << "," << v.y()[i] << "," << v.z()[i] << ") ";
+          }
+          os << "}";
+       }
+       return os;
+    } // op<< <>()
 
     template< class char_t, class traits_t, class T, class U >
       inline
diff --git a/math/genvector/inc/Math/GenVector/PtEtaPhiE4D.h b/math/genvector/inc/Math/GenVector/PtEtaPhiE4D.h
index 539e6ecca22fe..bc9a01d479f62 100644
--- a/math/genvector/inc/Math/GenVector/PtEtaPhiE4D.h
+++ b/math/genvector/inc/Math/GenVector/PtEtaPhiE4D.h
@@ -32,7 +32,7 @@
 #include <iostream>
 #endif
 
-
+#include <cmath>
 
 namespace ROOT {
 
@@ -137,15 +137,13 @@ public :
 
    // other coordinate representation
 
-   Scalar Px() const { return fPt*cos(fPhi);}
+   Scalar Px() const { return fPt * std::cos(fPhi); }
    Scalar X () const { return Px();         }
-   Scalar Py() const { return fPt*sin(fPhi);}
+   Scalar Py() const { return fPt * std::sin(fPhi); }
    Scalar Y () const { return Py();         }
    Scalar Pz() const {
-      return fPt >   0 ? fPt*std::sinh(fEta)     :
-         fEta == 0 ? 0                       :
-         fEta >  0 ? fEta - etaMax<Scalar>() :
-         fEta + etaMax<Scalar>() ;
+      return fPt > 0 ? fPt * std::sinh(fEta)
+                     : fEta == 0 ? 0 : fEta > 0 ? fEta - etaMax<Scalar>() : fEta + etaMax<Scalar>();
    }
    Scalar Z () const { return Pz(); }
 
@@ -153,29 +151,36 @@ public :
        magnitude of momentum
    */
    Scalar P() const {
-      return  fPt  > 0                 ?  fPt*std::cosh(fEta)       :
-         fEta >  etaMax<Scalar>() ?  fEta - etaMax<Scalar>()   :
-         fEta < -etaMax<Scalar>() ? -fEta - etaMax<Scalar>()   :
-         0                         ;
+      return fPt > 0 ? fPt * std::cosh(fEta)
+                     : fEta > etaMax<Scalar>() ? fEta - etaMax<Scalar>()
+                                               : fEta < -etaMax<Scalar>() ? -fEta - etaMax<Scalar>() : 0;
    }
    Scalar R() const { return P(); }
 
    /**
        squared magnitude of spatial components (momentum squared)
    */
-   Scalar P2() const { Scalar p = P(); return p*p; }
+   Scalar P2() const
+   {
+      const Scalar p = P();
+      return p * p;
+   }
 
    /**
       vector magnitude squared (or mass squared)
    */
-   Scalar M2() const { Scalar p = P(); return fE*fE - p*p; }
+   Scalar M2() const
+   {
+      const Scalar p = P();
+      return fE * fE - p * p;
+   }
    Scalar Mag2() const { return M2(); }
 
    /**
       invariant mass
    */
    Scalar M() const    {
-      Scalar mm = M2();
+      const Scalar mm = M2();
       if (mm >= 0) {
          return std::sqrt(mm);
       } else {
@@ -201,7 +206,7 @@ public :
       transverse mass
    */
    Scalar Mt() const {
-      Scalar mm = Mt2();
+      const Scalar mm = Mt2();
       if (mm >= 0) {
          return std::sqrt(mm);
       } else {
@@ -224,15 +229,16 @@ public :
    /**
        transverse energy squared
    */
-   Scalar Et2() const { Scalar et = Et(); return et*et; }
-
+   Scalar Et2() const
+   {
+      const Scalar et = Et();
+      return et * et;
+   }
 
 private:
    inline static Scalar pi() { return M_PI; }
    inline void Restrict() {
-      if ( fPhi <= -pi() || fPhi > pi() )
-         fPhi = fPhi - std::floor( fPhi/(2*pi()) +.5 ) * 2*pi();
-      return;
+      if (fPhi <= -pi() || fPhi > pi()) fPhi = fPhi - std::floor(fPhi / (2 * pi()) + .5) * 2 * pi();
    }
 public:
 
@@ -240,9 +246,8 @@ public :
       polar angle
    */
    Scalar Theta() const {
-      if (fPt  >  0) return 2* std::atan( exp( - fEta ) );
-      if (fEta >= 0) return 0;
-      return pi();
+     return ( fPt  > 0  ? Scalar(2) * std::atan(std::exp(-fEta)) :
+              fEta >= 0 ? 0 : pi() );
    }
 
    // --------- Set Coordinates of this system  ---------------
diff --git a/math/genvector/inc/Math/GenVector/PtEtaPhiM4D.h b/math/genvector/inc/Math/GenVector/PtEtaPhiM4D.h
index 5055c33e0eb2d..277713424ae14 100644
--- a/math/genvector/inc/Math/GenVector/PtEtaPhiM4D.h
+++ b/math/genvector/inc/Math/GenVector/PtEtaPhiM4D.h
@@ -30,6 +30,8 @@
 #include <iostream>
 #endif
 
+#include <cmath>
+
 namespace ROOT {
 
 namespace Math {
@@ -141,22 +143,20 @@ public :
        in this coordinate system it can be negagative if set that way.
    */
    Scalar M()   const { return fM;   }
-   Scalar Mag() const    { return M(); }
+   Scalar Mag() const { return M(); }
 
    Scalar Perp()const { return Pt(); }
    Scalar Rho() const { return Pt(); }
 
    // other coordinate representation
 
-   Scalar Px() const { return fPt*cos(fPhi);}
+   Scalar Px() const { return fPt * std::cos(fPhi); }
    Scalar X () const { return Px();         }
-   Scalar Py() const { return fPt*sin(fPhi);}
+   Scalar Py() const { return fPt * std::sin(fPhi); }
    Scalar Y () const { return Py();         }
    Scalar Pz() const {
-      return fPt >   0 ? fPt*std::sinh(fEta)     :
-         fEta == 0 ? 0                       :
-         fEta >  0 ? fEta - etaMax<Scalar>() :
-         fEta + etaMax<Scalar>() ;
+      return fPt > 0 ? fPt * std::sinh(fEta)
+                     : fEta == 0 ? 0 : fEta > 0 ? fEta - etaMax<Scalar>() : fEta + etaMax<Scalar>();
    }
    Scalar Z () const { return Pz(); }
 
@@ -164,17 +164,20 @@ public :
        magnitude of momentum
    */
    Scalar P() const {
-      return  fPt  > 0                 ?  fPt*std::cosh(fEta)       :
-         fEta >  etaMax<Scalar>() ?  fEta - etaMax<Scalar>()   :
-         fEta < -etaMax<Scalar>() ? -fEta - etaMax<Scalar>()   :
-         0                         ;
+      return fPt > 0 ? fPt * std::cosh(fEta)
+                     : fEta > etaMax<Scalar>() ? fEta - etaMax<Scalar>()
+                                               : fEta < -etaMax<Scalar>() ? -fEta - etaMax<Scalar>() : 0;
    }
    Scalar R() const { return P(); }
 
    /**
        squared magnitude of spatial components (momentum squared)
    */
-   Scalar P2() const { Scalar p = P(); return p*p; }
+   Scalar P2() const
+   {
+      const Scalar p = P();
+      return p * p;
+   }
 
    /**
        energy squared
@@ -188,7 +191,7 @@ public :
    /**
        Energy (timelike component of momentum-energy 4-vector)
    */
-   Scalar E()   const { return std::sqrt(E2() ); }
+   Scalar E() const { return std::sqrt(E2()); }
 
    Scalar T()   const { return E();  }
 
@@ -216,7 +219,7 @@ public :
       transverse mass - will be negative if Mt2() is negative
    */
    Scalar Mt() const {
-      Scalar mm = Mt2();
+      const Scalar mm = Mt2();
       if (mm >= 0) {
          return std::sqrt(mm);
       } else {
@@ -231,32 +234,30 @@ public :
    */
    Scalar Et2() const {
       // a bit faster than et * et
-      return 2. * E2()/ ( std::cosh(2 * fEta) + 1 );
+      return 2. * E2() / (std::cosh(2 * fEta) + 1);
    }
 
    /**
       transverse energy
    */
    Scalar Et() const {
-      return E() / std::cosh(fEta);
+     return E() / std::cosh(fEta);
    }
 
 private:
    inline static Scalar pi() { return M_PI; }
    inline void RestrictPhi() {
-      if ( fPhi <= -pi() || fPhi > pi() )
-         fPhi = fPhi - std::floor( fPhi/(2*pi()) +.5 ) * 2*pi();
-      return;
+      if (fPhi <= -pi() || fPhi > pi()) fPhi = fPhi - std::floor(fPhi / (2 * pi()) + .5) * 2 * pi();
    }
    // restrict the value of negative mass to avoid unphysical negative E2 values
    // M2 must be less than P2 for the tachionic particles - otherwise use positive values
    inline void RestrictNegMass() {
-      if ( fM >=0 ) return;
-      if ( P2() - fM*fM  < 0 ) {
-         GenVector::Throw ("PtEtaPhiM4D::unphysical value of mass, set to closest physical value");
-         fM = - P();
+      if (fM < 0) {
+         if (P2() - fM * fM < 0) {
+            GenVector::Throw("PtEtaPhiM4D::unphysical value of mass, set to closest physical value");
+            fM = -P();
+         }
       }
-      return;
    }
 
 public:
@@ -265,9 +266,8 @@ public :
       polar angle
    */
    Scalar Theta() const {
-      if (fPt  >  0) return 2* std::atan( exp( - fEta ) );
-      if (fEta >= 0) return 0;
-      return pi();
+     return ( fPt  > 0  ? Scalar(2) * std::atan(std::exp(-fEta)) :
+              fEta >= 0 ? 0 : pi() );
    }
 
    // --------- Set Coordinates of this system  ---------------
diff --git a/math/genvector/inc/Math/GenVector/PxPyPzE4D.h b/math/genvector/inc/Math/GenVector/PxPyPzE4D.h
index ec2a48d983f3d..6361482864951 100644
--- a/math/genvector/inc/Math/GenVector/PxPyPzE4D.h
+++ b/math/genvector/inc/Math/GenVector/PxPyPzE4D.h
@@ -147,8 +147,9 @@ public :
    /**
       invariant mass
    */
-   Scalar M() const    {
-      Scalar mm = M2();
+   Scalar M() const
+   {
+      const Scalar mm = M2();
       if (mm >= 0) {
          return std::sqrt(mm);
       } else {
@@ -168,7 +169,7 @@ public :
    /**
       Transverse spatial component (P_perp or rho)
    */
-   Scalar Pt()   const { return std::sqrt(Perp2());}
+   Scalar Pt() const { return std::sqrt(Perp2()); }
    Scalar Perp() const { return Pt();}
    Scalar Rho()  const { return Pt();}
 
@@ -181,7 +182,7 @@ public :
       transverse mass
    */
    Scalar Mt() const {
-      Scalar mm = Mt2();
+      const Scalar mm = Mt2();
       if (mm >= 0) {
          return std::sqrt(mm);
       } else {
@@ -204,7 +205,7 @@ public :
       transverse energy
    */
    Scalar Et() const {
-      Scalar etet = Et2();
+      const Scalar etet = Et2();
       return fT < 0.0 ? -std::sqrt(etet) : std::sqrt(etet);
    }
 
@@ -212,14 +213,14 @@ public :
       azimuthal angle
    */
    Scalar Phi() const  {
-      return (fX == 0.0 && fY == 0.0) ? 0 : std::atan2(fY,fX);
+      return (fX == 0.0 && fY == 0.0) ? 0 : std::atan2(fY, fX);
    }
 
    /**
       polar angle
    */
    Scalar Theta() const {
-      return (fX == 0.0 && fY == 0.0 && fZ == 0.0) ? 0 : std::atan2(Pt(),fZ);
+      return (fX == 0.0 && fY == 0.0 && fZ == 0.0) ? 0 : std::atan2(Pt(), fZ);
    }
 
    /**
diff --git a/math/genvector/inc/Math/GenVector/PxPyPzM4D.h b/math/genvector/inc/Math/GenVector/PxPyPzM4D.h
index 711c14990ba22..157edffd51821 100644
--- a/math/genvector/inc/Math/GenVector/PxPyPzM4D.h
+++ b/math/genvector/inc/Math/GenVector/PxPyPzM4D.h
@@ -154,7 +154,7 @@ public :
    /**
       Energy
     */
-   Scalar E()  const { return std::sqrt(E2() ); }
+   Scalar E() const { return std::sqrt(E2()); }
 
    Scalar T() const { return E();}
 
@@ -198,7 +198,7 @@ public :
    /**
       Transverse spatial component (P_perp or rho)
    */
-   Scalar Pt()   const { return std::sqrt(Perp2());}
+   Scalar Pt() const { return std::sqrt(Perp2()); }
    Scalar Perp() const { return Pt();}
    Scalar Rho()  const { return Pt();}
 
@@ -211,7 +211,7 @@ public :
       transverse mass
    */
    Scalar Mt() const {
-      Scalar mm = Mt2();
+      const Scalar mm = Mt2();
       if (mm >= 0) {
          return std::sqrt(mm);
       } else {
@@ -234,7 +234,7 @@ public :
       transverse energy
    */
    Scalar Et() const {
-      Scalar etet = Et2();
+      const Scalar etet = Et2();
       return std::sqrt(etet);
    }
 
@@ -242,14 +242,14 @@ public :
       azimuthal angle
    */
    Scalar Phi() const  {
-      return (fX == 0.0 && fY == 0.0) ? 0.0 : std::atan2(fY,fX);
+      return (fX == 0.0 && fY == 0.0) ? 0.0 : std::atan2(fY, fX);
    }
 
    /**
       polar angle
    */
    Scalar Theta() const {
-      return (fX == 0.0 && fY == 0.0 && fZ == 0.0) ? 0 : std::atan2(Pt(),fZ);
+      return (fX == 0.0 && fY == 0.0 && fZ == 0.0) ? 0 : std::atan2(Pt(), fZ);
    }
 
    /**
diff --git a/math/genvector/inc/Math/GenVector/RotationX.h b/math/genvector/inc/Math/GenVector/RotationX.h
index e4768eef07f3c..46724188cd275 100644
--- a/math/genvector/inc/Math/GenVector/RotationX.h
+++ b/math/genvector/inc/Math/GenVector/RotationX.h
@@ -94,13 +94,13 @@ class RotationX {
    /**
       Get the angle
    */
-   void GetAngle ( Scalar & angle ) const { angle = atan2 (fSin,fCos); }
+   void GetAngle(Scalar &angle) const { angle = std::atan2(fSin, fCos); }
    void GetComponents ( Scalar & angle ) const { GetAngle(angle); }
 
    /**
       Angle of rotation
    */
-   Scalar Angle () const { return atan2 (fSin,fCos); }
+   Scalar Angle() const { return std::atan2(fSin, fCos); }
 
    /**
       Sine or Cosine of the rotation angle
diff --git a/math/genvector/inc/Math/GenVector/RotationY.h b/math/genvector/inc/Math/GenVector/RotationY.h
index af148e4245859..0181a472f52a5 100644
--- a/math/genvector/inc/Math/GenVector/RotationY.h
+++ b/math/genvector/inc/Math/GenVector/RotationY.h
@@ -94,13 +94,13 @@ class RotationY {
    /**
       Get the angle
    */
-   void GetAngle ( Scalar & angle ) const { angle = atan2 (fSin,fCos); }
+   void GetAngle(Scalar &angle) const { angle = std::atan2(fSin, fCos); }
    void GetComponents ( Scalar & angle ) const { GetAngle(angle); }
 
    /**
       Angle of rotation
    */
-   Scalar Angle () const { return atan2 (fSin,fCos); }
+   Scalar Angle() const { return std::atan2(fSin, fCos); }
 
    /**
       Sine or Cosine of the rotation angle
diff --git a/math/genvector/inc/Math/GenVector/RotationZ.h b/math/genvector/inc/Math/GenVector/RotationZ.h
index 1752f1f2cde9f..819f819cac62d 100644
--- a/math/genvector/inc/Math/GenVector/RotationZ.h
+++ b/math/genvector/inc/Math/GenVector/RotationZ.h
@@ -94,13 +94,13 @@ class RotationZ {
    /**
       Get the angle
    */
-   void GetAngle ( Scalar & angle ) const { angle = atan2 (fSin,fCos); }
+   void GetAngle(Scalar &angle) const { angle = std::atan2(fSin, fCos); }
    void GetComponents ( Scalar & angle ) const { GetAngle(angle); }
 
    /**
       Angle of rotation
    */
-   Scalar Angle () const { return atan2 (fSin,fCos); }
+   Scalar Angle() const { return std::atan2(fSin, fCos); }
 
    /**
       Sine or Cosine of the rotation angle
diff --git a/math/genvector/inc/Math/GenVector/Transform3D.h b/math/genvector/inc/Math/GenVector/Transform3D.h
index 1071594a8f1d1..54a2dee8580d6 100644
--- a/math/genvector/inc/Math/GenVector/Transform3D.h
+++ b/math/genvector/inc/Math/GenVector/Transform3D.h
@@ -36,6 +36,8 @@
 #include "Math/GenVector/RotationZfwd.h"
 
 #include <iostream>
+#include <type_traits>
+#include <cmath>
 
 //#include "Math/Vector3Dfwd.h"
 
@@ -45,9 +47,7 @@ namespace ROOT {
 
 namespace Math {
 
-
-   class Plane3D;
-
+namespace Impl {
 
 //_________________________________________________________________________________________
 /**
@@ -74,14 +74,14 @@ namespace Math {
 
 */
 
+template <typename T = double>
 class Transform3D {
 
-
 public:
+   typedef T Scalar;
 
-   typedef  DisplacementVector3D<Cartesian3D<double>, DefaultCoordinateSystemTag >  Vector;
-   typedef  PositionVector3D<Cartesian3D<double>, DefaultCoordinateSystemTag >      Point;
-
+   typedef DisplacementVector3D<Cartesian3D<T>, DefaultCoordinateSystemTag> Vector;
+   typedef PositionVector3D<Cartesian3D<T>, DefaultCoordinateSystemTag>     Point;
 
    enum ETransform3DMatrixIndex {
       kXX = 0, kXY = 1, kXZ = 2, kDX = 3,
@@ -119,10 +119,7 @@ class Transform3D {
    /**
       Construct from a rotation and then a translation described by a Translation3D class
    */
-   Transform3D( const Rotation3D & r, const Translation3D & t)
-   {
-      AssignFrom( r, t.Vect() );
-   }
+   Transform3D(const Rotation3D &r, const Translation3D<T> &t) { AssignFrom(r, t.Vect()); }
 
    /**
       Construct from a rotation (any rotation object)  and then a translation
@@ -143,7 +140,7 @@ class Transform3D {
       Rotation3D class
    */
    template <class ARotation>
-   Transform3D( const ARotation & r, const Translation3D & t)
+   Transform3D(const ARotation &r, const Translation3D<T> &t)
    {
       AssignFrom( Rotation3D(r), t.Vect() );
    }
@@ -214,11 +211,7 @@ class Transform3D {
       Construct from a translation only, represented by a Translation3D class
       and with an identity rotation
    */
-   explicit Transform3D( const Translation3D & t) {
-      AssignFrom(t.Vect());
-   }
-
-
+   explicit Transform3D(const Translation3D<T> &t) { AssignFrom(t.Vect()); }
 
    //#if !defined(__MAKECINT__) && !defined(G__DICTIONARY)  // this is ambigous with double * , double *
 
@@ -239,23 +232,197 @@ class Transform3D {
    }
 #endif
 
-
+public:
    /**
       Construct transformation from one coordinate system defined by three
       points (origin + two axis) to
       a new coordinate system defined by other three points (origin + axis)
+      Scalar version.
       @param fr0  point defining origin of original reference system
       @param fr1  point defining first axis of original reference system
       @param fr2  point defining second axis of original reference system
       @param to0  point defining origin of transformed reference system
       @param to1  point defining first axis transformed reference system
       @param to2  point defining second axis transformed reference system
+   */
+   template <typename SCALAR = T, typename std::enable_if<std::is_arithmetic<SCALAR>::value>::type * = nullptr>
+   Transform3D(const Point &fr0, const Point &fr1, const Point &fr2, const Point &to0, const Point &to1,
+               const Point &to2)
+   {
+      // takes impl. from CLHEP ( E.Chernyaev). To be checked
+
+      Vector x1 = (fr1 - fr0).Unit();
+      Vector y1 = (fr2 - fr0).Unit();
+      Vector x2 = (to1 - to0).Unit();
+      Vector y2 = (to2 - to0).Unit();
+
+      //   C H E C K   A N G L E S
+
+      const T cos1 = x1.Dot(y1);
+      const T cos2 = x2.Dot(y2);
+
+      if (std::fabs(T(1) - cos1) <= T(0.000001) || std::fabs(T(1) - cos2) <= T(0.000001)) {
+         std::cerr << "Transform3D: Error : zero angle between axes" << std::endl;
+         SetIdentity();
+      } else {
+         if (std::fabs(cos1 - cos2) > T(0.000001)) {
+            std::cerr << "Transform3D: Warning: angles between axes are not equal" << std::endl;
+         }
+
+         //   F I N D   R O T A T I O N   M A T R I X
+
+         Vector z1 = (x1.Cross(y1)).Unit();
+         y1        = z1.Cross(x1);
+
+         Vector z2 = (x2.Cross(y2)).Unit();
+         y2        = z2.Cross(x2);
+
+         T x1x = x1.x();
+         T x1y = x1.y();
+         T x1z = x1.z();
+         T y1x = y1.x();
+         T y1y = y1.y();
+         T y1z = y1.z();
+         T z1x = z1.x();
+         T z1y = z1.y();
+         T z1z = z1.z();
+
+         T x2x = x2.x();
+         T x2y = x2.y();
+         T x2z = x2.z();
+         T y2x = y2.x();
+         T y2y = y2.y();
+         T y2z = y2.z();
+         T z2x = z2.x();
+         T z2y = z2.y();
+         T z2z = z2.z();
+
+         T detxx = (y1y * z1z - z1y * y1z);
+         T detxy = -(y1x * z1z - z1x * y1z);
+         T detxz = (y1x * z1y - z1x * y1y);
+         T detyx = -(x1y * z1z - z1y * x1z);
+         T detyy = (x1x * z1z - z1x * x1z);
+         T detyz = -(x1x * z1y - z1x * x1y);
+         T detzx = (x1y * y1z - y1y * x1z);
+         T detzy = -(x1x * y1z - y1x * x1z);
+         T detzz = (x1x * y1y - y1x * x1y);
+
+         T txx = x2x * detxx + y2x * detyx + z2x * detzx;
+         T txy = x2x * detxy + y2x * detyy + z2x * detzy;
+         T txz = x2x * detxz + y2x * detyz + z2x * detzz;
+         T tyx = x2y * detxx + y2y * detyx + z2y * detzx;
+         T tyy = x2y * detxy + y2y * detyy + z2y * detzy;
+         T tyz = x2y * detxz + y2y * detyz + z2y * detzz;
+         T tzx = x2z * detxx + y2z * detyx + z2z * detzx;
+         T tzy = x2z * detxy + y2z * detyy + z2z * detzy;
+         T tzz = x2z * detxz + y2z * detyz + z2z * detzz;
+
+         //   S E T    T R A N S F O R M A T I O N
+
+         T dx1 = fr0.x(), dy1 = fr0.y(), dz1 = fr0.z();
+         T dx2 = to0.x(), dy2 = to0.y(), dz2 = to0.z();
+
+         SetComponents(txx, txy, txz, dx2 - txx * dx1 - txy * dy1 - txz * dz1, tyx, tyy, tyz,
+                       dy2 - tyx * dx1 - tyy * dy1 - tyz * dz1, tzx, tzy, tzz, dz2 - tzx * dx1 - tzy * dy1 - tzz * dz1);
+      }
+   }
 
+   /**
+      Construct transformation from one coordinate system defined by three
+      points (origin + two axis) to
+      a new coordinate system defined by other three points (origin + axis)
+      Vectorised version.
+      @param fr0  point defining origin of original reference system
+      @param fr1  point defining first axis of original reference system
+      @param fr2  point defining second axis of original reference system
+      @param to0  point defining origin of transformed reference system
+      @param to1  point defining first axis transformed reference system
+      @param to2  point defining second axis transformed reference system
    */
-   Transform3D
-   (const Point & fr0, const Point & fr1, const Point & fr2,
-    const Point & to0, const Point & to1, const Point & to2 );
+   template <typename SCALAR = T, typename std::enable_if<!std::is_arithmetic<SCALAR>::value>::type * = nullptr>
+   Transform3D(const Point &fr0, const Point &fr1, const Point &fr2, const Point &to0, const Point &to1,
+               const Point &to2)
+   {
+      // takes impl. from CLHEP ( E.Chernyaev). To be checked
+
+      Vector x1 = (fr1 - fr0).Unit();
+      Vector y1 = (fr2 - fr0).Unit();
+      Vector x2 = (to1 - to0).Unit();
+      Vector y2 = (to2 - to0).Unit();
 
+      //   C H E C K   A N G L E S
+
+      const T cos1 = x1.Dot(y1);
+      const T cos2 = x2.Dot(y2);
+
+      const auto m1 = (abs(T(1) - cos1) <= T(0.000001) || abs(T(1) - cos2) <= T(0.000001));
+
+      const auto m2 = (abs(cos1 - cos2) > T(0.000001));
+      if (any_of(m2)) {
+         std::cerr << "Transform3D: Warning: angles between axes are not equal" << std::endl;
+      }
+
+      //   F I N D   R O T A T I O N   M A T R I X
+
+      Vector z1 = (x1.Cross(y1)).Unit();
+      y1        = z1.Cross(x1);
+
+      Vector z2 = (x2.Cross(y2)).Unit();
+      y2        = z2.Cross(x2);
+
+      T x1x = x1.x();
+      T x1y = x1.y();
+      T x1z = x1.z();
+      T y1x = y1.x();
+      T y1y = y1.y();
+      T y1z = y1.z();
+      T z1x = z1.x();
+      T z1y = z1.y();
+      T z1z = z1.z();
+
+      T x2x = x2.x();
+      T x2y = x2.y();
+      T x2z = x2.z();
+      T y2x = y2.x();
+      T y2y = y2.y();
+      T y2z = y2.z();
+      T z2x = z2.x();
+      T z2y = z2.y();
+      T z2z = z2.z();
+
+      T detxx = (y1y * z1z - z1y * y1z);
+      T detxy = -(y1x * z1z - z1x * y1z);
+      T detxz = (y1x * z1y - z1x * y1y);
+      T detyx = -(x1y * z1z - z1y * x1z);
+      T detyy = (x1x * z1z - z1x * x1z);
+      T detyz = -(x1x * z1y - z1x * x1y);
+      T detzx = (x1y * y1z - y1y * x1z);
+      T detzy = -(x1x * y1z - y1x * x1z);
+      T detzz = (x1x * y1y - y1x * x1y);
+
+      T txx = x2x * detxx + y2x * detyx + z2x * detzx;
+      T txy = x2x * detxy + y2x * detyy + z2x * detzy;
+      T txz = x2x * detxz + y2x * detyz + z2x * detzz;
+      T tyx = x2y * detxx + y2y * detyx + z2y * detzx;
+      T tyy = x2y * detxy + y2y * detyy + z2y * detzy;
+      T tyz = x2y * detxz + y2y * detyz + z2y * detzz;
+      T tzx = x2z * detxx + y2z * detyx + z2z * detzx;
+      T tzy = x2z * detxy + y2z * detyy + z2z * detzy;
+      T tzz = x2z * detxz + y2z * detyz + z2z * detzz;
+
+      //   S E T    T R A N S F O R M A T I O N
+
+      T dx1 = fr0.x(), dy1 = fr0.y(), dz1 = fr0.z();
+      T dx2 = to0.x(), dy2 = to0.y(), dz2 = to0.z();
+
+      SetComponents(txx, txy, txz, dx2 - txx * dx1 - txy * dy1 - txz * dz1, tyx, tyy, tyz,
+                    dy2 - tyx * dx1 - tyy * dy1 - tyz * dz1, tzx, tzy, tzz, dz2 - tzx * dx1 - tzy * dy1 - tzz * dz1);
+
+      if (any_of(m1)) {
+         std::cerr << "Transform3D: Error : zero angle between axes" << std::endl;
+         SetIdentity(m1);
+      }
+   }
 
    // use compiler generated copy ctor, copy assignmet and dtor
 
@@ -273,9 +440,7 @@ class Transform3D {
    /**
       Raw constructor from 12 Scalar components
    */
-   Transform3D(double  xx, double  xy, double  xz, double dx,
-               double  yx, double  yy, double  yz, double dy,
-               double  zx, double  zy, double  zz, double dz)
+   Transform3D(T xx, T xy, T xz, T dx, T yx, T yy, T yz, T dy, T zx, T zy, T zz, T dz)
    {
       SetComponents (xx, xy, xz, dx, yx, yy, yz, dy, zx, zy, zz, dz);
    }
@@ -287,8 +452,9 @@ class Transform3D {
       The 3x3 sub-block is assumed to be the rotation part and the translations vector
       are described by the 4-th column
    */
-   template<class ForeignMatrix>
-   Transform3D & operator= (const ForeignMatrix & m) {
+   template <class ForeignMatrix>
+   Transform3D<T> &operator=(const ForeignMatrix &m)
+   {
       SetComponents(m);
       return *this;
    }
@@ -336,7 +502,7 @@ class Transform3D {
    */
    template<class IT>
    void GetComponents(IT begin) const {
-      std::copy ( fM, fM+12, begin );
+      std::copy(fM, fM + 12, begin);
    }
 
    /**
@@ -370,10 +536,8 @@ class Transform3D {
    /**
       Set the components from 12 scalars
    */
-   void
-   SetComponents (double  xx, double  xy, double  xz, double dx,
-                  double  yx, double  yy, double  yz, double dy,
-                  double  zx, double  zy, double  zz, double dz) {
+   void SetComponents(T xx, T xy, T xz, T dx, T yx, T yy, T yz, T dy, T zx, T zy, T zz, T dz)
+   {
       fM[kXX]=xx;  fM[kXY]=xy;  fM[kXZ]=xz;  fM[kDX]=dx;
       fM[kYX]=yx;  fM[kYY]=yy;  fM[kYZ]=yz;  fM[kDY]=dy;
       fM[kZX]=zx;  fM[kZY]=zy;  fM[kZZ]=zz;  fM[kDZ]=dz;
@@ -382,10 +546,8 @@ class Transform3D {
    /**
       Get the components into 12 scalars
    */
-   void
-   GetComponents (double &xx, double &xy, double &xz, double &dx,
-                  double &yx, double &yy, double &yz, double &dy,
-                  double &zx, double &zy, double &zz, double &dz) const {
+   void GetComponents(T &xx, T &xy, T &xz, T &dx, T &yx, T &yy, T &yz, T &dy, T &zx, T &zy, T &zz, T &dz) const
+   {
       xx=fM[kXX];  xy=fM[kXY];  xz=fM[kXZ];  dx=fM[kDX];
       yx=fM[kYX];  yy=fM[kYY];  yz=fM[kYZ];  dy=fM[kDY];
       zx=fM[kZX];  zy=fM[kZY];  zz=fM[kZZ];  dz=fM[kDZ];
@@ -425,9 +587,7 @@ class Transform3D {
    */
    template <class AnyRotation>
    AnyRotation Rotation() const {
-      return AnyRotation(Rotation3D(fM[kXX], fM[kXY], fM[kXZ],
-                                    fM[kYX], fM[kYY], fM[kYZ],
-                                    fM[kZX], fM[kZY], fM[kZZ] ) );
+      return AnyRotation(Rotation3D(fM[kXX], fM[kXY], fM[kXZ], fM[kYX], fM[kYY], fM[kYZ], fM[kZX], fM[kZY], fM[kZZ]));
    }
 
    /**
@@ -441,9 +601,7 @@ class Transform3D {
    /**
       Get the translation representing the 3D transformation in a Cartesian vector
    */
-   Translation3D Translation() const {
-      return Translation3D( fM[kDX], fM[kDY], fM[kDZ] );
-   }
+   Translation3D<T> Translation() const { return Translation3D<T>(fM[kDX], fM[kDY], fM[kDZ]); }
 
    /**
       Get the translation representing the 3D transformation in any vector
@@ -483,10 +641,18 @@ class Transform3D {
    /**
       Transformation operation for Position Vector in any coordinate system
    */
-   template<class CoordSystem >
-   PositionVector3D<CoordSystem> operator() (const PositionVector3D <CoordSystem> & p) const {
-      Point xyzNew = operator() ( Point(p) );
-      return  PositionVector3D<CoordSystem> (xyzNew);
+   template <class CoordSystem>
+   PositionVector3D<CoordSystem> operator()(const PositionVector3D<CoordSystem> &p) const
+   {
+      return PositionVector3D<CoordSystem>(operator()(Point(p)));
+   }
+   /**
+      Transformation operation for Position Vector in any coordinate system
+   */
+   template <class CoordSystem>
+   PositionVector3D<CoordSystem> operator*(const PositionVector3D<CoordSystem> &v) const
+   {
+      return operator()(v);
    }
 
    /**
@@ -494,16 +660,24 @@ class Transform3D {
    */
    template<class CoordSystem >
    DisplacementVector3D<CoordSystem> operator() (const DisplacementVector3D <CoordSystem> & v) const {
-      Vector xyzNew = operator() ( Vector(v) );
-      return  DisplacementVector3D<CoordSystem> (xyzNew);
+      return DisplacementVector3D<CoordSystem>(operator()(Vector(v)));
+   }
+   /**
+      Transformation operation for Displacement Vector in any coordinate system
+   */
+   template <class CoordSystem>
+   DisplacementVector3D<CoordSystem> operator*(const DisplacementVector3D<CoordSystem> &v) const
+   {
+      return operator()(v);
    }
 
    /**
       Transformation operation for points between different coordinate system tags
    */
-   template<class CoordSystem, class Tag1, class Tag2 >
-   void Transform (const PositionVector3D <CoordSystem,Tag1> & p1, PositionVector3D <CoordSystem,Tag2> & p2  ) const {
-      Point xyzNew = operator() ( Point(p1.X(), p1.Y(), p1.Z()) );
+   template <class CoordSystem, class Tag1, class Tag2>
+   void Transform(const PositionVector3D<CoordSystem, Tag1> &p1, PositionVector3D<CoordSystem, Tag2> &p2) const
+   {
+      const Point xyzNew = operator()(Point(p1.X(), p1.Y(), p1.Z()));
       p2.SetXYZ( xyzNew.X(), xyzNew.Y(), xyzNew.Z() );
    }
 
@@ -511,9 +685,10 @@ class Transform3D {
    /**
       Transformation operation for Displacement Vector of different coordinate systems
    */
-   template<class CoordSystem,  class Tag1, class Tag2 >
-   void Transform (const DisplacementVector3D <CoordSystem,Tag1> & v1, DisplacementVector3D <CoordSystem,Tag2> & v2  ) const {
-      Vector xyzNew = operator() ( Vector(v1.X(), v1.Y(), v1.Z() ) );
+   template <class CoordSystem, class Tag1, class Tag2>
+   void Transform(const DisplacementVector3D<CoordSystem, Tag1> &v1, DisplacementVector3D<CoordSystem, Tag2> &v2) const
+   {
+      const Vector xyzNew = operator()(Vector(v1.X(), v1.Y(), v1.Z()));
       v2.SetXYZ( xyzNew.X(), xyzNew.Y(), xyzNew.Z() );
    }
 
@@ -522,112 +697,260 @@ class Transform3D {
    */
    template <class CoordSystem >
    LorentzVector<CoordSystem> operator() (const LorentzVector<CoordSystem> & q) const {
-      Vector xyzNew = operator() ( Vector(q.Vect() ) );
-      return  LorentzVector<CoordSystem> (xyzNew.X(), xyzNew.Y(), xyzNew.Z(), q.E() );
+      const Vector xyzNew = operator()(Vector(q.Vect()));
+      return LorentzVector<CoordSystem>(xyzNew.X(), xyzNew.Y(), xyzNew.Z(), q.E());
    }
-
    /**
-      Transformation on a 3D plane
+      Transformation operation for a Lorentz Vector in any  coordinate system
    */
-   Plane3D operator() (const Plane3D & plane) const;
-
-
-   // skip transformation for arbitrary vectors - not really defined if point or displacement vectors
+   template <class CoordSystem>
+   LorentzVector<CoordSystem> operator*(const LorentzVector<CoordSystem> &q) const
+   {
+      return operator()(q);
+   }
 
-   // same but with operator *
    /**
-      Transformation operation for Vectors. Apply same rules as operator()
-      depending on type of vector.
-      Will work only for DisplacementVector3D, PositionVector3D and LorentzVector
+      Transformation on a 3D plane
    */
-   template<class AVector >
-   AVector operator * (const AVector & v) const {
-      return operator() (v);
+   Plane3D<T> operator()(const Plane3D<T> &plane) const
+   {
+      // transformations on a 3D plane
+      const Vector n = plane.Normal();
+      // take a point on the plane. Use origin projection on the plane
+      // ( -ad, -bd, -cd) if (a**2 + b**2 + c**2 ) = 1
+      const T d = plane.HesseDistance();
+      Point   p(-d * n.X(), -d * n.Y(), -d * n.Z());
+      return Plane3D<T>(operator()(n), operator()(p));
    }
 
-
+   // skip transformation for arbitrary vectors - not really defined if point or displacement vectors
 
    /**
       multiply (combine) with another transformation in place
    */
-   inline Transform3D & operator *= (const Transform3D  & t);
+   inline Transform3D<T> &operator*=(const Transform3D<T> &t);
 
    /**
       multiply (combine) two transformations
    */
-   inline Transform3D operator * (const Transform3D  & t) const;
+   inline Transform3D<T> operator*(const Transform3D<T> &t) const;
+
+   /**
+       Invert the transformation in place (scalar)
+   */
+   template <typename SCALAR = T, typename std::enable_if<std::is_arithmetic<SCALAR>::value>::type * = nullptr>
+   void Invert()
+   {
+      //
+      // Name: Transform3D::inverse                     Date:    24.09.96
+      // Author: E.Chernyaev (IHEP/Protvino)            Revised:
+      //
+      // Function: Find inverse affine transformation.
+
+      T detxx = fM[kYY] * fM[kZZ] - fM[kYZ] * fM[kZY];
+      T detxy = fM[kYX] * fM[kZZ] - fM[kYZ] * fM[kZX];
+      T detxz = fM[kYX] * fM[kZY] - fM[kYY] * fM[kZX];
+      T det   = fM[kXX] * detxx - fM[kXY] * detxy + fM[kXZ] * detxz;
+      if (det == T(0)) {
+         std::cerr << "Transform3D::inverse error: zero determinant" << std::endl;
+         return;
+      }
+      det = T(1) / det;
+      detxx *= det;
+      detxy *= det;
+      detxz *= det;
+      T detyx = (fM[kXY] * fM[kZZ] - fM[kXZ] * fM[kZY]) * det;
+      T detyy = (fM[kXX] * fM[kZZ] - fM[kXZ] * fM[kZX]) * det;
+      T detyz = (fM[kXX] * fM[kZY] - fM[kXY] * fM[kZX]) * det;
+      T detzx = (fM[kXY] * fM[kYZ] - fM[kXZ] * fM[kYY]) * det;
+      T detzy = (fM[kXX] * fM[kYZ] - fM[kXZ] * fM[kYX]) * det;
+      T detzz = (fM[kXX] * fM[kYY] - fM[kXY] * fM[kYX]) * det;
+      SetComponents(detxx, -detyx, detzx, -detxx * fM[kDX] + detyx * fM[kDY] - detzx * fM[kDZ], -detxy, detyy, -detzy,
+                    detxy * fM[kDX] - detyy * fM[kDY] + detzy * fM[kDZ], detxz, -detyz, detzz,
+                    -detxz * fM[kDX] + detyz * fM[kDY] - detzz * fM[kDZ]);
+   }
 
    /**
-       Invert the transformation in place
+       Invert the transformation in place (vectorised)
    */
-   void Invert();
+   template <typename SCALAR = T, typename std::enable_if<!std::is_arithmetic<SCALAR>::value>::type * = nullptr>
+   void Invert()
+   {
+      //
+      // Name: Transform3D::inverse                     Date:    24.09.96
+      // Author: E.Chernyaev (IHEP/Protvino)            Revised:
+      //
+      // Function: Find inverse affine transformation.
+
+      T          detxx    = fM[kYY] * fM[kZZ] - fM[kYZ] * fM[kZY];
+      T          detxy    = fM[kYX] * fM[kZZ] - fM[kYZ] * fM[kZX];
+      T          detxz    = fM[kYX] * fM[kZY] - fM[kYY] * fM[kZX];
+      T          det      = fM[kXX] * detxx - fM[kXY] * detxy + fM[kXZ] * detxz;
+      const auto detZmask = (det == T(0));
+      if (any_of(detZmask)) {
+         std::cerr << "Transform3D::inverse error: zero determinant" << std::endl;
+         det(detZmask) = T(1);
+      }
+      det = T(1) / det;
+      detxx *= det;
+      detxy *= det;
+      detxz *= det;
+      T detyx = (fM[kXY] * fM[kZZ] - fM[kXZ] * fM[kZY]) * det;
+      T detyy = (fM[kXX] * fM[kZZ] - fM[kXZ] * fM[kZX]) * det;
+      T detyz = (fM[kXX] * fM[kZY] - fM[kXY] * fM[kZX]) * det;
+      T detzx = (fM[kXY] * fM[kYZ] - fM[kXZ] * fM[kYY]) * det;
+      T detzy = (fM[kXX] * fM[kYZ] - fM[kXZ] * fM[kYX]) * det;
+      T detzz = (fM[kXX] * fM[kYY] - fM[kXY] * fM[kYX]) * det;
+      // Set det=0 cases to 0
+      if (any_of(detZmask)) {
+         detxx(detZmask) = T(0);
+         detxy(detZmask) = T(0);
+         detxz(detZmask) = T(0);
+         detyx(detZmask) = T(0);
+         detyy(detZmask) = T(0);
+         detyz(detZmask) = T(0);
+         detzx(detZmask) = T(0);
+         detzy(detZmask) = T(0);
+         detzz(detZmask) = T(0);
+      }
+      // set final components
+      SetComponents(detxx, -detyx, detzx, -detxx * fM[kDX] + detyx * fM[kDY] - detzx * fM[kDZ], -detxy, detyy, -detzy,
+                    detxy * fM[kDX] - detyy * fM[kDY] + detzy * fM[kDZ], detxz, -detyz, detzz,
+                    -detxz * fM[kDX] + detyz * fM[kDY] - detzz * fM[kDZ]);
+   }
 
    /**
       Return the inverse of the transformation.
    */
-   Transform3D Inverse() const {
-      Transform3D t(*this);
+   Transform3D<T> Inverse() const
+   {
+      Transform3D<T> t(*this);
       t.Invert();
       return t;
    }
 
-
    /**
       Equality operator. Check equality for each element
-      To do: use double tolerance
+      To do: use T tolerance
    */
-   bool operator == (const Transform3D & rhs) const {
-      if( fM[0] != rhs.fM[0] )  return false;
-      if( fM[1] != rhs.fM[1] )  return false;
-      if( fM[2] != rhs.fM[2] )  return false;
-      if( fM[3] != rhs.fM[3] )  return false;
-      if( fM[4] != rhs.fM[4] )  return false;
-      if( fM[5] != rhs.fM[5] )  return false;
-      if( fM[6] != rhs.fM[6] )  return false;
-      if( fM[7] != rhs.fM[7] )  return false;
-      if( fM[8] != rhs.fM[8] )  return false;
-      if( fM[9] != rhs.fM[9] )  return false;
-      if( fM[10]!= rhs.fM[10] ) return false;
-      if( fM[11]!= rhs.fM[11] ) return false;
-      return true;
+   bool operator==(const Transform3D<T> &rhs) const
+   {
+      return (fM[0] == rhs.fM[0] && fM[1] == rhs.fM[1] && fM[2] == rhs.fM[2] && fM[3] == rhs.fM[3] &&
+              fM[4] == rhs.fM[4] && fM[5] == rhs.fM[5] && fM[6] == rhs.fM[6] && fM[7] == rhs.fM[7] &&
+              fM[8] == rhs.fM[8] && fM[9] == rhs.fM[9] && fM[10] == rhs.fM[10] && fM[11] == rhs.fM[11]);
    }
 
    /**
       Inequality operator. Check equality for each element
-      To do: use double tolerance
+      To do: use T tolerance
    */
-   bool operator != (const Transform3D & rhs) const {
-      return ! operator==(rhs);
-   }
-
+   bool operator!=(const Transform3D<T> &rhs) const { return !operator==(rhs); }
 
 protected:
 
    /**
       make transformation from first a rotation then a translation
    */
-   void  AssignFrom( const Rotation3D & r, const Vector & v);
+   void AssignFrom(const Rotation3D &r, const Vector &v)
+   {
+      // assignment  from rotation + translation
+
+      T rotData[9];
+      r.GetComponents(rotData, rotData + 9);
+      // first raw
+      for (int i = 0; i < 3; ++i) fM[i] = rotData[i];
+      // second raw
+      for (int i = 0; i < 3; ++i) fM[kYX + i] = rotData[3 + i];
+      // third raw
+      for (int i = 0; i < 3; ++i) fM[kZX + i] = rotData[6 + i];
+
+      // translation data
+      T vecData[3];
+      v.GetCoordinates(vecData, vecData + 3);
+      fM[kDX] = vecData[0];
+      fM[kDY] = vecData[1];
+      fM[kDZ] = vecData[2];
+   }
 
    /**
       make transformation from only rotations (zero translation)
    */
-   void  AssignFrom( const Rotation3D & r);
+   void AssignFrom(const Rotation3D &r)
+   {
+      // assign from only a rotation  (null translation)
+      T rotData[9];
+      r.GetComponents(rotData, rotData + 9);
+      for (int i = 0; i < 3; ++i) {
+         for (int j = 0; j < 3; ++j) fM[4 * i + j] = rotData[3 * i + j];
+         // empty vector data
+         fM[4 * i + 3] = T(0);
+      }
+   }
 
    /**
       make transformation from only translation (identity rotations)
    */
-   void  AssignFrom( const Vector & v);
+   void AssignFrom(const Vector &v)
+   {
+      // assign from a translation only (identity rotations)
+      fM[kXX] = T(1);
+      fM[kXY] = T(0);
+      fM[kXZ] = T(0);
+      fM[kDX] = v.X();
+      fM[kYX] = T(0);
+      fM[kYY] = T(1);
+      fM[kYZ] = T(0);
+      fM[kDY] = v.Y();
+      fM[kZX] = T(0);
+      fM[kZY] = T(0);
+      fM[kZZ] = T(1);
+      fM[kDZ] = v.Z();
+   }
 
    /**
       Set identity transformation (identity rotation , zero translation)
    */
-   void SetIdentity() ;
-
-private:
-
+   void SetIdentity()
+   {
+      // set identity ( identity rotation and zero translation)
+      fM[kXX] = T(1);
+      fM[kXY] = T(0);
+      fM[kXZ] = T(0);
+      fM[kDX] = T(0);
+      fM[kYX] = T(0);
+      fM[kYY] = T(1);
+      fM[kYZ] = T(0);
+      fM[kDY] = T(0);
+      fM[kZX] = T(0);
+      fM[kZY] = T(0);
+      fM[kZZ] = T(1);
+      fM[kDZ] = T(0);
+   }
 
-   double fM[12];    // transformation elements (3x4 matrix)
+   /**
+      Set identity transformation (identity rotation , zero translation)
+      vectorised version that sets using a mask
+   */
+   template <typename SCALAR = T, typename std::enable_if<!std::is_arithmetic<SCALAR>::value>::type * = nullptr>
+   void SetIdentity(const typename SCALAR::mask_type m)
+   {
+      // set identity ( identity rotation and zero translation)
+      fM[kXX](m) = T(1);
+      fM[kXY](m) = T(0);
+      fM[kXZ](m) = T(0);
+      fM[kDX](m) = T(0);
+      fM[kYX](m) = T(0);
+      fM[kYY](m) = T(1);
+      fM[kYZ](m) = T(0);
+      fM[kDY](m) = T(0);
+      fM[kZX](m) = T(0);
+      fM[kZY](m) = T(0);
+      fM[kZZ](m) = T(1);
+      fM[kDZ](m) = T(0);
+   }
 
+private:
+   T fM[12]; // transformation elements (3x4 matrix)
 };
 
 
@@ -635,7 +958,8 @@ class Transform3D {
 
 // inline functions (combination of transformations)
 
-inline Transform3D & Transform3D::operator *= (const Transform3D  & t)
+template <class T>
+inline Transform3D<T> &Transform3D<T>::operator*=(const Transform3D<T> &t)
 {
    // combination of transformations
 
@@ -657,27 +981,25 @@ inline Transform3D & Transform3D::operator *= (const Transform3D  & t)
    return *this;
 }
 
-
-
-inline Transform3D Transform3D::operator * (const Transform3D  & t) const
+template <class T>
+inline Transform3D<T> Transform3D<T>::operator*(const Transform3D<T> &t) const
 {
    // combination of transformations
 
-   return Transform3D(fM[kXX]*t.fM[kXX]+fM[kXY]*t.fM[kYX]+fM[kXZ]*t.fM[kZX],
-                      fM[kXX]*t.fM[kXY]+fM[kXY]*t.fM[kYY]+fM[kXZ]*t.fM[kZY],
-                      fM[kXX]*t.fM[kXZ]+fM[kXY]*t.fM[kYZ]+fM[kXZ]*t.fM[kZZ],
-                      fM[kXX]*t.fM[kDX]+fM[kXY]*t.fM[kDY]+fM[kXZ]*t.fM[kDZ]+fM[kDX],
-
-                      fM[kYX]*t.fM[kXX]+fM[kYY]*t.fM[kYX]+fM[kYZ]*t.fM[kZX],
-                      fM[kYX]*t.fM[kXY]+fM[kYY]*t.fM[kYY]+fM[kYZ]*t.fM[kZY],
-                      fM[kYX]*t.fM[kXZ]+fM[kYY]*t.fM[kYZ]+fM[kYZ]*t.fM[kZZ],
-                      fM[kYX]*t.fM[kDX]+fM[kYY]*t.fM[kDY]+fM[kYZ]*t.fM[kDZ]+fM[kDY],
+   return Transform3D<T>(fM[kXX] * t.fM[kXX] + fM[kXY] * t.fM[kYX] + fM[kXZ] * t.fM[kZX],
+                         fM[kXX] * t.fM[kXY] + fM[kXY] * t.fM[kYY] + fM[kXZ] * t.fM[kZY],
+                         fM[kXX] * t.fM[kXZ] + fM[kXY] * t.fM[kYZ] + fM[kXZ] * t.fM[kZZ],
+                         fM[kXX] * t.fM[kDX] + fM[kXY] * t.fM[kDY] + fM[kXZ] * t.fM[kDZ] + fM[kDX],
 
-                      fM[kZX]*t.fM[kXX]+fM[kZY]*t.fM[kYX]+fM[kZZ]*t.fM[kZX],
-                      fM[kZX]*t.fM[kXY]+fM[kZY]*t.fM[kYY]+fM[kZZ]*t.fM[kZY],
-                      fM[kZX]*t.fM[kXZ]+fM[kZY]*t.fM[kYZ]+fM[kZZ]*t.fM[kZZ],
-                      fM[kZX]*t.fM[kDX]+fM[kZY]*t.fM[kDY]+fM[kZZ]*t.fM[kDZ]+fM[kDZ]  );
+                         fM[kYX] * t.fM[kXX] + fM[kYY] * t.fM[kYX] + fM[kYZ] * t.fM[kZX],
+                         fM[kYX] * t.fM[kXY] + fM[kYY] * t.fM[kYY] + fM[kYZ] * t.fM[kZY],
+                         fM[kYX] * t.fM[kXZ] + fM[kYY] * t.fM[kYZ] + fM[kYZ] * t.fM[kZZ],
+                         fM[kYX] * t.fM[kDX] + fM[kYY] * t.fM[kDY] + fM[kYZ] * t.fM[kDZ] + fM[kDY],
 
+                         fM[kZX] * t.fM[kXX] + fM[kZY] * t.fM[kYX] + fM[kZZ] * t.fM[kZX],
+                         fM[kZX] * t.fM[kXY] + fM[kZY] * t.fM[kYY] + fM[kZZ] * t.fM[kZY],
+                         fM[kZX] * t.fM[kXZ] + fM[kZY] * t.fM[kYZ] + fM[kZZ] * t.fM[kZZ],
+                         fM[kZX] * t.fM[kDX] + fM[kZY] * t.fM[kDY] + fM[kZZ] * t.fM[kDZ] + fM[kDZ]);
 }
 
 
@@ -693,36 +1015,52 @@ inline Transform3D Transform3D::operator * (const Transform3D  & t) const
    combine a translation and a rotation to give a transform3d
    First the translation then the rotation
  */
-inline Transform3D operator * (const Rotation3D & r, const Translation3D & t) {
-   return Transform3D( r, r(t.Vect()) );
+template <class T>
+inline Transform3D<T> operator*(const Rotation3D &r, const Translation3D<T> &t)
+{
+   return Transform3D<T>(r, r(t.Vect()));
 }
-inline Transform3D operator * (const RotationX & r, const Translation3D & t) {
+template <class T>
+inline Transform3D<T> operator*(const RotationX &r, const Translation3D<T> &t)
+{
    Rotation3D r3(r);
-   return Transform3D( r3, r3(t.Vect()) );
+   return Transform3D<T>(r3, r3(t.Vect()));
 }
-inline Transform3D operator * (const RotationY & r, const Translation3D & t) {
+template <class T>
+inline Transform3D<T> operator*(const RotationY &r, const Translation3D<T> &t)
+{
    Rotation3D r3(r);
-   return Transform3D( r3, r3(t.Vect()) );
+   return Transform3D<T>(r3, r3(t.Vect()));
 }
-inline Transform3D operator * (const RotationZ & r, const Translation3D & t) {
+template <class T>
+inline Transform3D<T> operator*(const RotationZ &r, const Translation3D<T> &t)
+{
    Rotation3D r3(r);
-   return Transform3D( r3, r3(t.Vect()) );
+   return Transform3D<T>(r3, r3(t.Vect()));
 }
-inline Transform3D operator * (const RotationZYX & r, const Translation3D & t) {
+template <class T>
+inline Transform3D<T> operator*(const RotationZYX &r, const Translation3D<T> &t)
+{
    Rotation3D r3(r);
-   return Transform3D( r3, r3(t.Vect()) );
+   return Transform3D<T>(r3, r3(t.Vect()));
 }
-inline Transform3D operator * (const AxisAngle & r, const Translation3D & t) {
+template <class T>
+inline Transform3D<T> operator*(const AxisAngle &r, const Translation3D<T> &t)
+{
    Rotation3D r3(r);
-   return Transform3D( r3, r3(t.Vect()) );
+   return Transform3D<T>(r3, r3(t.Vect()));
 }
-inline Transform3D operator * (const EulerAngles & r, const Translation3D & t) {
+template <class T>
+inline Transform3D<T> operator*(const EulerAngles &r, const Translation3D<T> &t)
+{
    Rotation3D r3(r);
-   return Transform3D( r3, r3(t.Vect()) );
+   return Transform3D<T>(r3, r3(t.Vect()));
 }
-inline Transform3D operator * (const Quaternion & r, const Translation3D & t) {
+template <class T>
+inline Transform3D<T> operator*(const Quaternion &r, const Translation3D<T> &t)
+{
    Rotation3D r3(r);
-   return Transform3D( r3, r3(t.Vect()) );
+   return Transform3D<T>(r3, r3(t.Vect()));
 }
 
 // ------ combination of a  rotation (first)  and then a translation ------
@@ -731,29 +1069,45 @@ inline Transform3D operator * (const Quaternion & r, const Translation3D & t) {
    combine a rotation and a translation to give a transform3d
    First a rotation then the translation
  */
-inline Transform3D operator * (const Translation3D & t, const Rotation3D & r) {
-   return Transform3D( r, t.Vect());
+template <class T>
+inline Transform3D<T> operator*(const Translation3D<T> &t, const Rotation3D &r)
+{
+   return Transform3D<T>(r, t.Vect());
 }
-inline Transform3D operator * (const Translation3D & t, const RotationX & r) {
-   return Transform3D( Rotation3D(r) , t.Vect());
+template <class T>
+inline Transform3D<T> operator*(const Translation3D<T> &t, const RotationX &r)
+{
+   return Transform3D<T>(Rotation3D(r), t.Vect());
 }
-inline Transform3D operator * (const Translation3D & t, const RotationY & r) {
-   return Transform3D( Rotation3D(r) , t.Vect());
+template <class T>
+inline Transform3D<T> operator*(const Translation3D<T> &t, const RotationY &r)
+{
+   return Transform3D<T>(Rotation3D(r), t.Vect());
 }
-inline Transform3D operator * (const Translation3D & t, const RotationZ & r) {
-   return Transform3D( Rotation3D(r) , t.Vect());
+template <class T>
+inline Transform3D<T> operator*(const Translation3D<T> &t, const RotationZ &r)
+{
+   return Transform3D<T>(Rotation3D(r), t.Vect());
 }
-inline Transform3D operator * (const Translation3D & t, const RotationZYX & r) {
-   return Transform3D( Rotation3D(r) , t.Vect());
+template <class T>
+inline Transform3D<T> operator*(const Translation3D<T> &t, const RotationZYX &r)
+{
+   return Transform3D<T>(Rotation3D(r), t.Vect());
 }
-inline Transform3D operator * (const Translation3D & t, const EulerAngles & r) {
-   return Transform3D( Rotation3D(r) , t.Vect());
+template <class T>
+inline Transform3D<T> operator*(const Translation3D<T> &t, const EulerAngles &r)
+{
+   return Transform3D<T>(Rotation3D(r), t.Vect());
 }
-inline Transform3D operator * (const Translation3D & t, const Quaternion & r) {
-   return Transform3D( Rotation3D(r) , t.Vect());
+template <class T>
+inline Transform3D<T> operator*(const Translation3D<T> &t, const Quaternion &r)
+{
+   return Transform3D<T>(Rotation3D(r), t.Vect());
 }
-inline Transform3D operator * (const Translation3D & t, const AxisAngle & r) {
-   return Transform3D( Rotation3D(r) , t.Vect());
+template <class T>
+inline Transform3D<T> operator*(const Translation3D<T> &t, const AxisAngle &r)
+{
+   return Transform3D<T>(Rotation3D(r), t.Vect());
 }
 
 // ------ combination of a Transform3D and a pure translation------
@@ -762,17 +1116,21 @@ inline Transform3D operator * (const Translation3D & t, const AxisAngle & r) {
    combine a transformation and a translation to give a transform3d
    First the translation then the transform3D
  */
-inline Transform3D operator * (const Transform3D & t, const Translation3D & d) {
+template <class T>
+inline Transform3D<T> operator*(const Transform3D<T> &t, const Translation3D<T> &d)
+{
    Rotation3D r = t.Rotation();
-   return Transform3D( r, r( d.Vect() ) + t.Translation().Vect()  );
+   return Transform3D<T>(r, r(d.Vect()) + t.Translation().Vect());
 }
 
 /**
    combine a translation and a transformation to give a transform3d
    First the transformation then the translation
  */
-inline Transform3D operator * (const Translation3D & d, const Transform3D & t) {
-   return Transform3D( t.Rotation(), t.Translation().Vect() + d.Vect());
+template <class T>
+inline Transform3D<T> operator*(const Translation3D<T> &d, const Transform3D<T> &t)
+{
+   return Transform3D<T>(t.Rotation(), t.Translation().Vect() + d.Vect());
 }
 
 // ------ combination of a Transform3D and any rotation------
@@ -782,29 +1140,45 @@ inline Transform3D operator * (const Translation3D & d, const Transform3D & t) {
    combine a transformation and a rotation to give a transform3d
    First the rotation then the transform3D
  */
-inline Transform3D operator * (const Transform3D & t, const Rotation3D & r) {
-   return Transform3D( t.Rotation()*r ,  t.Translation()  );
+template <class T>
+inline Transform3D<T> operator*(const Transform3D<T> &t, const Rotation3D &r)
+{
+   return Transform3D<T>(t.Rotation() * r, t.Translation());
 }
-inline Transform3D operator * (const Transform3D & t, const RotationX & r) {
-   return Transform3D( t.Rotation()*r ,  t.Translation()  );
+template <class T>
+inline Transform3D<T> operator*(const Transform3D<T> &t, const RotationX &r)
+{
+   return Transform3D<T>(t.Rotation() * r, t.Translation());
 }
-inline Transform3D operator * (const Transform3D & t, const RotationY & r) {
-   return Transform3D( t.Rotation()*r ,  t.Translation()  );
+template <class T>
+inline Transform3D<T> operator*(const Transform3D<T> &t, const RotationY &r)
+{
+   return Transform3D<T>(t.Rotation() * r, t.Translation());
 }
-inline Transform3D operator * (const Transform3D & t, const RotationZ & r) {
-   return Transform3D( t.Rotation()*r ,  t.Translation()  );
+template <class T>
+inline Transform3D<T> operator*(const Transform3D<T> &t, const RotationZ &r)
+{
+   return Transform3D<T>(t.Rotation() * r, t.Translation());
 }
-inline Transform3D operator * (const Transform3D & t, const RotationZYX & r) {
-   return Transform3D( t.Rotation()*r ,  t.Translation()  );
+template <class T>
+inline Transform3D<T> operator*(const Transform3D<T> &t, const RotationZYX &r)
+{
+   return Transform3D<T>(t.Rotation() * r, t.Translation());
 }
-inline Transform3D operator * (const Transform3D & t, const EulerAngles & r) {
-   return Transform3D( t.Rotation()*r ,  t.Translation()  );
+template <class T>
+inline Transform3D<T> operator*(const Transform3D<T> &t, const EulerAngles &r)
+{
+   return Transform3D<T>(t.Rotation() * r, t.Translation());
 }
-inline Transform3D operator * (const Transform3D & t, const AxisAngle & r) {
-   return Transform3D( t.Rotation()*r ,  t.Translation()  );
+template <class T>
+inline Transform3D<T> operator*(const Transform3D<T> &t, const AxisAngle &r)
+{
+   return Transform3D<T>(t.Rotation() * r, t.Translation());
 }
-inline Transform3D operator * (const Transform3D & t, const Quaternion & r) {
-   return Transform3D( t.Rotation()*r ,  t.Translation()  );
+template <class T>
+inline Transform3D<T> operator*(const Transform3D<T> &t, const Quaternion &r)
+{
+   return Transform3D<T>(t.Rotation() * r, t.Translation());
 }
 
 
@@ -813,36 +1187,52 @@ inline Transform3D operator * (const Transform3D & t, const Quaternion & r) {
    combine a rotation and a transformation to give a transform3d
    First the transformation then the rotation
  */
-inline Transform3D operator * (const Rotation3D & r, const Transform3D & t) {
-   return Transform3D( r * t.Rotation(), r * t.Translation().Vect() );
+template <class T>
+inline Transform3D<T> operator*(const Rotation3D &r, const Transform3D<T> &t)
+{
+   return Transform3D<T>(r * t.Rotation(), r * t.Translation().Vect());
 }
-inline Transform3D operator * (const RotationX & r, const Transform3D & t) {
+template <class T>
+inline Transform3D<T> operator*(const RotationX &r, const Transform3D<T> &t)
+{
    Rotation3D r3d(r);
-   return Transform3D( r3d * t.Rotation(), r3d * t.Translation().Vect() );
+   return Transform3D<T>(r3d * t.Rotation(), r3d * t.Translation().Vect());
 }
-inline Transform3D operator * (const RotationY & r, const Transform3D & t) {
+template <class T>
+inline Transform3D<T> operator*(const RotationY &r, const Transform3D<T> &t)
+{
    Rotation3D r3d(r);
-   return Transform3D( r3d * t.Rotation(), r3d * t.Translation().Vect() );
+   return Transform3D<T>(r3d * t.Rotation(), r3d * t.Translation().Vect());
 }
-inline Transform3D operator * (const RotationZ & r, const Transform3D & t) {
+template <class T>
+inline Transform3D<T> operator*(const RotationZ &r, const Transform3D<T> &t)
+{
    Rotation3D r3d(r);
-   return Transform3D( r3d * t.Rotation(), r3d * t.Translation().Vect() );
+   return Transform3D<T>(r3d * t.Rotation(), r3d * t.Translation().Vect());
 }
-inline Transform3D operator * (const RotationZYX & r, const Transform3D & t) {
+template <class T>
+inline Transform3D<T> operator*(const RotationZYX &r, const Transform3D<T> &t)
+{
    Rotation3D r3d(r);
-   return Transform3D( r3d * t.Rotation(), r3d * t.Translation().Vect() );
+   return Transform3D<T>(r3d * t.Rotation(), r3d * t.Translation().Vect());
 }
-inline Transform3D operator * (const EulerAngles & r, const Transform3D & t) {
+template <class T>
+inline Transform3D<T> operator*(const EulerAngles &r, const Transform3D<T> &t)
+{
    Rotation3D r3d(r);
-   return Transform3D( r3d * t.Rotation(), r3d * t.Translation().Vect() );
+   return Transform3D<T>(r3d * t.Rotation(), r3d * t.Translation().Vect());
 }
-inline Transform3D operator * (const AxisAngle & r, const Transform3D & t) {
+template <class T>
+inline Transform3D<T> operator*(const AxisAngle &r, const Transform3D<T> &t)
+{
    Rotation3D r3d(r);
-   return Transform3D( r3d * t.Rotation(), r3d * t.Translation().Vect() );
+   return Transform3D<T>(r3d * t.Rotation(), r3d * t.Translation().Vect());
 }
-inline Transform3D operator * (const Quaternion & r, const Transform3D & t) {
+template <class T>
+inline Transform3D<T> operator*(const Quaternion &r, const Transform3D<T> &t)
+{
    Rotation3D r3d(r);
-   return Transform3D( r3d * t.Rotation(), r3d * t.Translation().Vect() );
+   return Transform3D<T>(r3d * t.Rotation(), r3d * t.Translation().Vect());
 }
 
 
@@ -852,10 +1242,27 @@ inline Transform3D operator * (const Quaternion & r, const Transform3D & t) {
 /**
    print the 12 components of the Transform3D
  */
-std::ostream & operator<< (std::ostream & os, const Transform3D & t);
+template <class T>
+std::ostream &operator<<(std::ostream &os, const Transform3D<T> &t)
+{
+   // TODO - this will need changing for machine-readable issues
+   //        and even the human readable form needs formatting improvements
+
+   T m[12];
+   t.GetComponents(m, m + 12);
+   os << "\n" << m[0] << "  " << m[1] << "  " << m[2] << "  " << m[3];
+   os << "\n" << m[4] << "  " << m[5] << "  " << m[6] << "  " << m[7];
+   os << "\n" << m[8] << "  " << m[9] << "  " << m[10] << "  " << m[11] << "\n";
+   return os;
+}
+
+} // end namespace Impl
 
+// typedefs for double and float versions
+typedef Impl::Transform3D<double> Transform3D;
+typedef Impl::Transform3D<float>  Transform3DF;
 
-   } // end namespace Math
+} // end namespace Math
 
 } // end namespace ROOT
 
diff --git a/math/genvector/inc/Math/GenVector/Translation3D.h b/math/genvector/inc/Math/GenVector/Translation3D.h
index b7d0d8feff7aa..c0c6a4a2a5447 100644
--- a/math/genvector/inc/Math/GenVector/Translation3D.h
+++ b/math/genvector/inc/Math/GenVector/Translation3D.h
@@ -19,21 +19,20 @@
 
 #include "Math/GenVector/DisplacementVector3D.h"
 
+#include "Math/GenVector/Plane3D.h"
+
 #include "Math/GenVector/PositionVector3Dfwd.h"
 
 #include "Math/GenVector/LorentzVectorfwd.h"
 
 #include <iostream>
-
-
+#include <type_traits>
 
 namespace ROOT {
 
 namespace Math {
 
-
-   class Plane3D;
-
+namespace Impl {
 
 //____________________________________________________________________________________________________
 /**
@@ -48,15 +47,13 @@ namespace Math {
 
 */
 
+template <typename T = double>
 class Translation3D {
 
+public:
+   typedef T Scalar;
 
-  public:
-
-    typedef  DisplacementVector3D<Cartesian3D<double>, DefaultCoordinateSystemTag >  Vector;
-
-
-
+   typedef DisplacementVector3D<Cartesian3D<T>, DefaultCoordinateSystemTag> Vector;
 
    /**
        Default constructor ( zero translation )
@@ -76,10 +73,7 @@ class Translation3D {
    /**
       Construct from x,y,z values representing the translation
    */
-   Translation3D(double dx, double dy, double dz) :
-      fVect( Vector(dx, dy, dz) )
-   {  }
-
+   Translation3D(T dx, T dy, T dz) : fVect(Vector(dx, dy, dz)) {}
 
    /**
       Construct from any Displacement vector in ant tag and coordinate system
@@ -97,9 +91,9 @@ class Translation3D {
       @param p2  point defining origin of transformed reference system
 
    */
-   template<class CoordSystem, class Tag>
-   Translation3D (const  PositionVector3D<CoordSystem,Tag> & p1, const PositionVector3D<CoordSystem,Tag> & p2 ) :
-      fVect(p2-p1)
+   template <class CoordSystem, class Tag>
+   Translation3D(const PositionVector3D<CoordSystem, Tag> &p1, const PositionVector3D<CoordSystem, Tag> &p2)
+      : fVect(p2 - p1)
    { }
 
 
@@ -143,28 +137,17 @@ class Translation3D {
    /**
       Set the components from 3 scalars
    */
-   void
-   SetComponents (double  dx, double  dy, double  dz ) {
-      fVect.SetCoordinates(dx,dy,dz);
-   }
+   void SetComponents(T dx, T dy, T dz) { fVect.SetCoordinates(dx, dy, dz); }
 
    /**
       Get the components into 3 scalars
    */
-   void
-   GetComponents (double &dx, double &dy, double &dz) const {
-      fVect.GetCoordinates(dx,dy,dz);
-   }
-
+   void GetComponents(T &dx, T &dy, T &dz) const { fVect.GetCoordinates(dx, dy, dz); }
 
    /**
       Set the XYZ vector components from 3 scalars
    */
-   void
-   SetXYZ (double  dx, double  dy, double  dz ) {
-      fVect.SetXYZ(dx,dy,dz);
-   }
-
+   void SetXYZ(T dx, T dy, T dz) { fVect.SetXYZ(dx, dy, dz); }
 
    // operations on points and vectors
 
@@ -174,11 +157,15 @@ class Translation3D {
    */
    template<class CoordSystem, class Tag >
    PositionVector3D<CoordSystem,Tag> operator() (const PositionVector3D <CoordSystem,Tag> & p) const {
-      PositionVector3D<CoordSystem,Tag>  tmp;
-      tmp.SetXYZ (p.X() + fVect.X(),
-                  p.Y() + fVect.Y(),
-                  p.Z() + fVect.Z() ) ;
-      return tmp;
+      return PositionVector3D<CoordSystem, Tag>(p.X() + fVect.X(), p.Y() + fVect.Y(), p.Z() + fVect.Z());
+   }
+   /**
+     Transformation operation
+   */
+   template <class CoordSystem, class Tag>
+   PositionVector3D<CoordSystem, Tag> operator*(const PositionVector3D<CoordSystem, Tag> &v) const
+   {
+      return operator()(v);
    }
 
    /**
@@ -189,6 +176,14 @@ class Translation3D {
    DisplacementVector3D<CoordSystem,Tag> operator() (const DisplacementVector3D <CoordSystem,Tag> & v) const {
       return  v;
    }
+   /**
+     Transformation operation
+   */
+   template <class CoordSystem, class Tag>
+   DisplacementVector3D<CoordSystem, Tag> operator*(const DisplacementVector3D<CoordSystem, Tag> &v) const
+   {
+      return operator()(v);
+   }
 
    /**
       Transformation operation for points between different coordinate system tags
@@ -200,47 +195,53 @@ class Translation3D {
       p2 =  operator()(tmp);
     }
 
-
    /**
       Transformation operation for Displacement Vector of different coordinate systems
    */
-   template<class CoordSystem,  class Tag1, class Tag2 >
-   void Transform (const DisplacementVector3D <CoordSystem,Tag1> & v1, DisplacementVector3D <CoordSystem,Tag2> & v2  ) const {
-      // just copy v1 in v2
-      v2.SetXYZ(v1.X(), v1.Y(), v1.Z() );
+    template <class CoordSystem, class Tag1, class Tag2>
+    void Transform(const DisplacementVector3D<CoordSystem, Tag1> &v1, DisplacementVector3D<CoordSystem, Tag2> &v2) const
+    {
+       // just copy v1 in v2
+       v2.SetXYZ(v1.X(), v1.Y(), v1.Z());
    }
 
    /**
       Transformation operation for a Lorentz Vector in any  coordinate system
       A LorentzVector contains a displacement vector so no translation applies as well
    */
-   template <class CoordSystem >
-   LorentzVector<CoordSystem> operator() (const LorentzVector<CoordSystem> & q) const {
-      return  q;
+   template <class CoordSystem>
+   LorentzVector<CoordSystem> operator()(const LorentzVector<CoordSystem> &q) const
+   {
+      return q;
    }
-
    /**
-      Transformation on a 3D plane
+     Transformation operation
    */
-   Plane3D operator() (const Plane3D & plane) const;
-
+   template <class CoordSystem>
+   LorentzVector<CoordSystem> operator*(const LorentzVector<CoordSystem> &q) const
+   {
+      return operator()(q);
+   }
 
    /**
-      Transformation operation for Vectors. Apply same rules as operator()
-      depending on type of vector.
-      Will work only for DisplacementVector3D, PositionVector3D and LorentzVector
+      Transformation on a 3D plane
    */
-   template<class AVector >
-   AVector operator * (const AVector & v) const {
-      return operator() (v);
+   Plane3D<T> operator()(const Plane3D<T> &plane) const
+   {
+      // transformations on a 3D plane
+      const Vector n = plane.Normal();
+      // take a point on the plane. Use origin projection on the plane
+      // ( -ad, -bd, -cd) if (a**2 + b**2 + c**2 ) = 1
+      const T                          d = plane.HesseDistance();
+      PositionVector3D<Cartesian3D<T>> p(-d * n.X(), -d * n.Y(), -d * n.Z());
+      return PLANE(operator()(n), operator()(p));
    }
 
-
-
    /**
       multiply (combine) with another transformation in place
    */
-   Translation3D & operator *= (const Translation3D  & t) {
+   Translation3D<T> &operator*=(const Translation3D<T> &t)
+   {
       fVect+= t.Vect();
       return *this;
    }
@@ -248,9 +249,7 @@ class Translation3D {
    /**
       multiply (combine) two transformations
    */
-   Translation3D operator * (const Translation3D  & t) const {
-      return Translation3D( fVect + t.Vect() );
-   }
+   Translation3D<T> operator*(const Translation3D<T> &t) const { return Translation3D<T>(fVect + t.Vect()); }
 
    /**
        Invert the transformation in place
@@ -262,24 +261,18 @@ class Translation3D {
    /**
       Return the inverse of the transformation.
    */
-   Translation3D Inverse() const {
-      return Translation3D( -fVect.X(), -fVect.Y(),-fVect.Z() );
-   }
-
+   Translation3D<T> Inverse() const { return Translation3D<T>(-fVect.X(), -fVect.Y(), -fVect.Z()); }
 
    /**
       Equality/inequality operators
    */
-   bool operator == (const Translation3D & rhs) const {
+   bool operator==(const Translation3D<T> &rhs) const
+   {
       if( fVect != rhs.fVect )  return false;
       return true;
    }
 
-   bool operator != (const Translation3D & rhs) const {
-      return ! operator==(rhs);
-   }
-
-
+   bool operator!=(const Translation3D<T> &rhs) const { return !operator==(rhs); }
 
 private:
 
@@ -295,11 +288,26 @@ class Translation3D {
 
 // TODO - I/O should be put in the manipulator form
 
-std::ostream & operator<< (std::ostream & os, const Translation3D & t);
+template <class T>
+std::ostream &operator<<(std::ostream &os, const Translation3D<T> &t)
+{
+   // TODO - this will need changing for machine-readable issues
+   //        and even the human readable form needs formatiing improvements
+
+   T m[3];
+   t.GetComponents(m, m + 3);
+   return os << "\n" << m[0] << "  " << m[1] << "  " << m[2] << "\n";
+}
 
 // need a function Transform = Translation * Rotation ???
 
-   } // end namespace Math
+} // end namespace Impl
+
+// typedefs for double and float versions
+typedef Impl::Translation3D<double> Translation3D;
+typedef Impl::Translation3D<float>  Translation3DF;
+
+} // end namespace Math
 
 } // end namespace ROOT
 
diff --git a/math/genvector/inc/Math/GenVector/VectorUtil.h b/math/genvector/inc/Math/GenVector/VectorUtil.h
index 07a2f6c99d516..4f3ac65427737 100644
--- a/math/genvector/inc/Math/GenVector/VectorUtil.h
+++ b/math/genvector/inc/Math/GenVector/VectorUtil.h
@@ -251,8 +251,8 @@ namespace ROOT {
           */
          template <class Vector>
          Vector RotateX(const Vector & v, double alpha) {
-            double sina = sin(alpha);
-            double cosa = cos(alpha);
+            double sina = std::sin(alpha);
+            double cosa = std::cos(alpha);
             double y2 = v.Y() * cosa - v.Z()*sina;
             double z2 = v.Z() * cosa + v.Y() * sina;
             Vector vrot;
@@ -268,8 +268,8 @@ namespace ROOT {
           */
          template <class Vector>
          Vector RotateY(const Vector & v, double alpha) {
-            double sina = sin(alpha);
-            double cosa = cos(alpha);
+            double sina = std::sin(alpha);
+            double cosa = std::cos(alpha);
             double x2 = v.X() * cosa + v.Z() * sina;
             double z2 = v.Z() * cosa - v.X() * sina;
             Vector vrot;
@@ -285,8 +285,8 @@ namespace ROOT {
           */
          template <class Vector>
          Vector RotateZ(const Vector & v, double alpha) {
-            double sina = sin(alpha);
-            double cosa = cos(alpha);
+            double sina = std::sin(alpha);
+            double cosa = std::cos(alpha);
             double x2 = v.X() * cosa - v.Y() * sina;
             double y2 = v.Y() * cosa + v.X() * sina;
             Vector vrot;
diff --git a/math/genvector/inc/Math/GenVector/eta.h b/math/genvector/inc/Math/GenVector/eta.h
index 434471771e734..24ee558efb9e9 100644
--- a/math/genvector/inc/Math/GenVector/eta.h
+++ b/math/genvector/inc/Math/GenVector/eta.h
@@ -50,14 +50,14 @@ namespace ROOT {
 
               // value to control Taylor expansion of sqrt
               static const Scalar big_z_scaled =
-                 std::pow(std::numeric_limits<Scalar>::epsilon(),static_cast<Scalar>(-.25));
+                 std::pow(std::numeric_limits<Scalar>::epsilon(), static_cast<Scalar>(-.25));
 
               Scalar z_scaled = z/rho;
               if (std::fabs(z_scaled) < big_z_scaled) {
-                 return std::log(z_scaled+std::sqrt(z_scaled*z_scaled+1.0));
+                 return std::log(z_scaled + std::sqrt(z_scaled * z_scaled + 1.0));
               } else {
                  // apply correction using first order Taylor expansion of sqrt
-                 return  z>0 ? std::log(2.0*z_scaled + 0.5/z_scaled) : -std::log(-2.0*z_scaled);
+                 return z > 0 ? std::log(2.0 * z_scaled + 0.5 / z_scaled) : -std::log(-2.0 * z_scaled);
               }
            }
            // case vector has rho = 0
@@ -80,8 +80,7 @@ namespace ROOT {
         */
         template<typename Scalar>
         inline Scalar Eta_FromTheta(Scalar theta, Scalar r) {
-
-           Scalar tanThetaOver2 = std::tan( theta/2.);
+           Scalar tanThetaOver2 = std::tan(theta / 2.);
            if (tanThetaOver2 == 0) {
               return r + etaMax<Scalar>();
            }
diff --git a/math/genvector/src/Plane3D.cxx b/math/genvector/src/Plane3D.cxx
deleted file mode 100644
index 6b8d25c02e65d..0000000000000
--- a/math/genvector/src/Plane3D.cxx
+++ /dev/null
@@ -1,109 +0,0 @@
-// @(#)root/mathcore:$Id$
-// Authors: W. Brown, M. Fischler, L. Moneta    2005
-
-/**********************************************************************
- *                                                                    *
- * Copyright (c) 2005 , LCG ROOT MathLib Team                         *
- *                                                                    *
- *                                                                    *
- **********************************************************************/
-
-// implementation file for class Plane3D
-//
-// Created by: Lorenzo Moneta  December 2 2005
-//
-//
-
-#include "Math/GenVector/Plane3D.h"
-
-#include <cmath>
-
-
-
-
-namespace ROOT {
-
-namespace Math {
-
-
-typedef Plane3D::Scalar Scalar;
-typedef Plane3D::Point  XYZPoint;
-typedef Plane3D::Vector XYZVector;
-
-// ========== Constructors and Assignment =====================
-
-
-// constructor from 4 scalars numbers (a,b,c,d)
-Plane3D::Plane3D(const Scalar & a, const Scalar & b, const Scalar & c, const Scalar & d) :
-   fA(a), fB(b), fC(c), fD(d)
-{
-   //renormalize a,b,c to unit
-   Normalize();
-}
-
-// internal method to construct from a normal vector and a point
-void Plane3D::BuildFromVecAndPoint(const XYZVector & n, const XYZPoint & p )
-{
-   // build from a normal vector and a point
-   fA =  n.X();
-   fB =  n.Y();
-   fC =  n.Z();
-   fD = - n.Dot(p);
-   Normalize();
-}
-
-// internl method to construct from three points
-void Plane3D::BuildFrom3Points( const XYZPoint & p1, const XYZPoint & p2, const XYZPoint & p3 ) {
-
-   // plane from thre points
-   // normal is (x3-x1) cross (x2 -x1)
-   XYZVector n = (p2-p1).Cross(p3-p1);
-   fA = n.X();
-   fB = n.Y();
-   fC = n.Z();
-   fD = - n.Dot(p1);
-   Normalize();
-}
-
-// distance plane- point
-Scalar Plane3D::Distance(const XYZPoint & p) const {
-   return fA*p.X() + fB*p.Y() + fC*p.Z() + fD;
-}
-
-void Plane3D::Normalize() {
-   // normalize the plane
-   Scalar s = std::sqrt( fA*fA + fB*fB + fC*fC );
-   // what to do if s = 0 ??
-   if ( s == 0) { fD = 0; return; }
-   Scalar w = 1./s;
-   fA *= w;
-   fB *= w;
-   fC *= w;
-   fD *= w;
-}
-
-
-// projection of a point onto the plane
-XYZPoint Plane3D::ProjectOntoPlane(const XYZPoint & p) const {
-   Scalar d = Distance(p);
-   return XYZPoint( p.X() - fA*d, p.Y() - fB*d, p.Z() - fC*d);
-}
-
-
-// output
-std::ostream & operator<< (std::ostream & os, const Plane3D & p) {
-   os << "\n" << p.Normal().X()
-   << "  " << p.Normal().Y()
-   << "  " << p.Normal().Z()
-   << "  " << p.HesseDistance()
-   << "\n";
-   return os;
-}
-
-
-
-
-
-}  // end namespace Math
-}  // end namespace ROOT
-
diff --git a/math/genvector/src/Transform3D.cxx b/math/genvector/src/Transform3D.cxx
deleted file mode 100644
index 0bc7c83e8e84b..0000000000000
--- a/math/genvector/src/Transform3D.cxx
+++ /dev/null
@@ -1,228 +0,0 @@
-// @(#)root/mathcore:$Id$
-// Authors: W. Brown, M. Fischler, L. Moneta    2005
-
-/**********************************************************************
- *                                                                    *
- * Copyright (c) 2005 , LCG ROOT MathLib Team                         *
- *                                                                    *
- *                                                                    *
- **********************************************************************/
-
-// implementation file for class Transform3D
-//
-// Created by: Lorenzo Moneta  October 27 2005
-//
-//
-#include "Math/GenVector/GenVectorIO.h"
-
-#include "Math/GenVector/Transform3D.h"
-#include "Math/GenVector/Plane3D.h"
-
-#include <cmath>
-#include <algorithm>
-
-
-
-
-namespace ROOT {
-
-namespace Math {
-
-
-typedef Transform3D::Point  XYZPoint;
-typedef Transform3D::Vector XYZVector;
-
-
-// ========== Constructors and Assignment =====================
-
-
-
-// construct from two ref frames
-Transform3D::Transform3D(const XYZPoint & fr0, const XYZPoint & fr1, const XYZPoint & fr2,
-                         const XYZPoint & to0, const XYZPoint & to1, const XYZPoint & to2 )
-{
-   // takes impl. from CLHEP ( E.Chernyaev). To be checked
-
-   XYZVector x1,y1,z1, x2,y2,z2;
-   x1 = (fr1 - fr0).Unit();
-   y1 = (fr2 - fr0).Unit();
-   x2 = (to1 - to0).Unit();
-   y2 = (to2 - to0).Unit();
-
-   //   C H E C K   A N G L E S
-
-   double cos1, cos2;
-   cos1 = x1.Dot(y1);
-   cos2 = x2.Dot(y2);
-
-   if (std::fabs(1.0-cos1) <= 0.000001 || std::fabs(1.0-cos2) <= 0.000001) {
-      std::cerr << "Transform3D: Error : zero angle between axes" << std::endl;
-      SetIdentity();
-   } else {
-      if (std::fabs(cos1-cos2) > 0.000001) {
-         std::cerr << "Transform3D: Warning: angles between axes are not equal"
-         << std::endl;
-      }
-
-      //   F I N D   R O T A T I O N   M A T R I X
-
-      z1 = (x1.Cross(y1)).Unit();
-      y1  = z1.Cross(x1);
-
-      z2 = (x2.Cross(y2)).Unit();
-      y2  = z2.Cross(x2);
-
-      double x1x = x1.x();      double x1y = x1.y();      double x1z = x1.z();
-      double y1x = y1.x();      double y1y = y1.y();      double y1z = y1.z();
-      double z1x = z1.x();      double z1y = z1.y();      double z1z = z1.z();
-
-      double x2x = x2.x();      double x2y = x2.y();      double x2z = x2.z();
-      double y2x = y2.x();      double y2y = y2.y();      double y2z = y2.z();
-      double z2x = z2.x();      double z2y = z2.y();      double z2z = z2.z();
-
-      double detxx =  (y1y *z1z  - z1y *y1z );
-      double detxy = -(y1x *z1z  - z1x *y1z );
-      double detxz =  (y1x *z1y  - z1x *y1y );
-      double detyx = -(x1y *z1z  - z1y *x1z );
-      double detyy =  (x1x *z1z  - z1x *x1z );
-      double detyz = -(x1x *z1y  - z1x *x1y );
-      double detzx =  (x1y *y1z  - y1y *x1z );
-      double detzy = -(x1x *y1z  - y1x *x1z );
-      double detzz =  (x1x *y1y  - y1x *x1y );
-
-      double txx = x2x *detxx + y2x *detyx + z2x *detzx;
-      double txy = x2x *detxy + y2x *detyy + z2x *detzy;
-      double txz = x2x *detxz + y2x *detyz + z2x *detzz;
-      double tyx = x2y *detxx + y2y *detyx + z2y *detzx;
-      double tyy = x2y *detxy + y2y *detyy + z2y *detzy;
-      double tyz = x2y *detxz + y2y *detyz + z2y *detzz;
-      double tzx = x2z *detxx + y2z *detyx + z2z *detzx;
-      double tzy = x2z *detxy + y2z *detyy + z2z *detzy;
-      double tzz = x2z *detxz + y2z *detyz + z2z *detzz;
-
-      //   S E T    T R A N S F O R M A T I O N
-
-      double dx1 = fr0.x(), dy1 = fr0.y(), dz1 = fr0.z();
-      double dx2 = to0.x(), dy2 = to0.y(), dz2 = to0.z();
-
-      SetComponents(txx, txy, txz, dx2-txx*dx1-txy*dy1-txz*dz1,
-                    tyx, tyy, tyz, dy2-tyx*dx1-tyy*dy1-tyz*dz1,
-                    tzx, tzy, tzz, dz2-tzx*dx1-tzy*dy1-tzz*dz1);
-   }
-}
-
-
-// inversion (from CLHEP)
-void Transform3D::Invert()
-{
-   //
-   // Name: Transform3D::inverse                     Date:    24.09.96
-   // Author: E.Chernyaev (IHEP/Protvino)            Revised:
-   //
-   // Function: Find inverse affine transformation.
-
-   double detxx = fM[kYY]*fM[kZZ] - fM[kYZ]*fM[kZY];
-   double detxy = fM[kYX]*fM[kZZ] - fM[kYZ]*fM[kZX];
-   double detxz = fM[kYX]*fM[kZY] - fM[kYY]*fM[kZX];
-   double det   = fM[kXX]*detxx - fM[kXY]*detxy + fM[kXZ]*detxz;
-   if (det == 0) {
-      std::cerr << "Transform3D::inverse error: zero determinant" << std::endl;
-      return;
-   }
-   det = 1./det; detxx *= det; detxy *= det; detxz *= det;
-   double detyx = (fM[kXY]*fM[kZZ] - fM[kXZ]*fM[kZY] )*det;
-   double detyy = (fM[kXX]*fM[kZZ] - fM[kXZ]*fM[kZX] )*det;
-   double detyz = (fM[kXX]*fM[kZY] - fM[kXY]*fM[kZX] )*det;
-   double detzx = (fM[kXY]*fM[kYZ] - fM[kXZ]*fM[kYY] )*det;
-   double detzy = (fM[kXX]*fM[kYZ] - fM[kXZ]*fM[kYX] )*det;
-   double detzz = (fM[kXX]*fM[kYY] - fM[kXY]*fM[kYX] )*det;
-   SetComponents
-      (detxx, -detyx,  detzx, -detxx*fM[kDX]+detyx*fM[kDY]-detzx*fM[kDZ],
-       -detxy,  detyy, -detzy,  detxy*fM[kDX]-detyy*fM[kDY]+detzy*fM[kDZ],
-       detxz, -detyz,  detzz, -detxz*fM[kDX]+detyz*fM[kDY]-detzz*fM[kDZ]);
-}
-
-
-
-
-void Transform3D::SetIdentity()
-{
-   //set identity ( identity rotation and zero translation)
-   fM[kXX] = 1.0;  fM[kXY] = 0.0; fM[kXZ] = 0.0; fM[kDX] = 0.0;
-   fM[kYX] = 0.0;  fM[kYY] = 1.0; fM[kYZ] = 0.0; fM[kDY] = 0.0;
-   fM[kZX] = 0.0;  fM[kZY] = 0.0; fM[kZZ] = 1.0; fM[kDZ] = 0.0;
-}
-
-
-void Transform3D::AssignFrom (const Rotation3D  & r,  const XYZVector & v)
-{
-   // assignment  from rotation + translation
-
-   double rotData[9];
-   r.GetComponents(rotData, rotData +9);
-   // first raw
-   for (int i = 0; i < 3; ++i)
-      fM[i] = rotData[i];
-   // second raw
-   for (int i = 0; i < 3; ++i)
-      fM[kYX+i] = rotData[3+i];
-   // third raw
-   for (int i = 0; i < 3; ++i)
-      fM[kZX+i] = rotData[6+i];
-
-   // translation data
-   double vecData[3];
-   v.GetCoordinates(vecData, vecData+3);
-   fM[kDX] = vecData[0];
-   fM[kDY] = vecData[1];
-   fM[kDZ] = vecData[2];
-}
-
-
-void Transform3D::AssignFrom(const Rotation3D & r)
-{
-   // assign from only a rotation  (null translation)
-   double rotData[9];
-   r.GetComponents(rotData, rotData +9);
-   for (int i = 0; i < 3; ++i) {
-      for (int j = 0; j < 3; ++j)
-         fM[4*i + j] = rotData[3*i+j];
-      // empty vector data
-      fM[4*i + 3] = 0;
-   }
-}
-
-void Transform3D::AssignFrom(const XYZVector & v)
-{
-   // assign from a translation only (identity rotations)
-   fM[kXX] = 1.0;  fM[kXY] = 0.0; fM[kXZ] = 0.0; fM[kDX] = v.X();
-   fM[kYX] = 0.0;  fM[kYY] = 1.0; fM[kYZ] = 0.0; fM[kDY] = v.Y();
-   fM[kZX] = 0.0;  fM[kZY] = 0.0; fM[kZZ] = 1.0; fM[kDZ] = v.Z();
-}
-
-Plane3D Transform3D::operator() (const Plane3D & plane) const
-{
-   // transformations on a 3D plane
-   XYZVector n = plane.Normal();
-   // take a point on the plane. Use origin projection on the plane
-   // ( -ad, -bd, -cd) if (a**2 + b**2 + c**2 ) = 1
-   double d = plane.HesseDistance();
-   XYZPoint p( - d * n.X() , - d *n.Y(), -d *n.Z() );
-   return Plane3D ( operator() (n), operator() (p) );
-}
-
-std::ostream & operator<< (std::ostream & os, const Transform3D & t)
-{
-   // TODO - this will need changing for machine-readable issues
-   //        and even the human readable form needs formatting improvements
-
-   double m[12];
-   t.GetComponents(m, m+12);
-   os << "\n" << m[0] << "  " << m[1] << "  " << m[2] << "  " << m[3] ;
-   os << "\n" << m[4] << "  " << m[5] << "  " << m[6] << "  " << m[7] ;
-   os << "\n" << m[8] << "  " << m[9] << "  " << m[10]<< "  " << m[11] << "\n";
-   return os;
-}
-
-}  // end namespace Math
-}  // end namespace ROOT
diff --git a/math/genvector/src/Translation3D.cxx b/math/genvector/src/Translation3D.cxx
deleted file mode 100644
index 2a17f3855cefb..0000000000000
--- a/math/genvector/src/Translation3D.cxx
+++ /dev/null
@@ -1,62 +0,0 @@
-// @(#)root/mathcore:$Id$
-// Authors: W. Brown, M. Fischler, L. Moneta    2005
-
-/**********************************************************************
- *                                                                    *
- * Copyright (c) 2005 , LCG ROOT MathLib Team                         *
- *                                                                    *
- *                                                                    *
- **********************************************************************/
-
-// implementation file for class Translation3D
-//
-// Created by: Lorenzo Moneta  October 27 2005
-//
-//
-
-#include "Math/GenVector/Translation3D.h"
-#include "Math/GenVector/Plane3D.h"
-#include "Math/GenVector/PositionVector3D.h"
-
-#include <cmath>
-#include <algorithm>
-
-
-
-
-namespace ROOT {
-
-namespace Math {
-
-
-typedef Translation3D::Vector XYZVector;
-typedef PositionVector3D<Cartesian3D<double> > XYZPoint;
-
-
-// ========== Constructors and Assignment =====================
-
-
-Plane3D Translation3D::operator() (const Plane3D & plane) const
-{
-   // transformations on a 3D plane
-   XYZVector n = plane.Normal();
-   // take a point on the plane. Use origin projection on the plane
-   // ( -ad, -bd, -cd) if (a**2 + b**2 + c**2 ) = 1
-   double d = plane.HesseDistance();
-   XYZPoint p( - d * n.X() , - d *n.Y(), -d *n.Z() );
-   return Plane3D ( operator() (n), operator() (p) );
-}
-
-std::ostream & operator<< (std::ostream & os, const Translation3D & t)
-{
-   // TODO - this will need changing for machine-readable issues
-   //        and even the human readable form needs formatiing improvements
-
-   double m[3];
-   t.GetComponents(m, m+3);
-   os << "\n" << m[0] << "  " << m[1] << "  " << m[2] << "\n";
-   return os;
-}
-
-}  // end namespace Math
-}  // end namespace ROOT
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 214c10ddc0d7e..ce569ac3889a3 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -310,6 +310,20 @@ if(ROOT_vc_FOUND)
   ROOT_ADD_TEST(test-Vc COMMAND testVc FAILREGEX "FAILED|Error in")
 endif()
 
+#--Vc GenVector test-----------------------------------------------------------------------------------
+if(ROOT_vc_FOUND)
+  include_directories(${Vc_INCLUDE_DIRS})
+  
+  ROOT_ADD_CXX_FLAG(vc_flags -Wno-unused-parameter)
+  ROOT_ADD_CXX_FLAG(vc_flags -Wno-missing-braces)
+  ROOT_ADD_CXX_FLAG(vc_flags -Wno-conditional-uninitialized)
+  ROOT_ADD_CXX_FLAG(vc_flags -Wno-mismatched-tags)
+  ROOT_ADD_CXX_FLAG(vc_flags -Wno-undefined-var-template)
+  set_source_files_properties(testGenVectorVc.cxx PROPERTIES COMPILE_FLAGS "${vc_flags}")
+  ROOT_EXECUTABLE(testGenVectorVc testGenVectorVc.cxx LIBRARIES Physics GenVector ${Vc_LIBRARIES} BUILTINS Vc)
+  ROOT_ADD_TEST(test-GenVector-Vc COMMAND testGenVectorVc FAILREGEX "FAILED|Error in")
+endif()
+
 #--g2root------------------------------------------------------------------------------------------
 if(TARGET g2root)
   ROOT_ADD_TEST(test-g2root COMMAND g2root)
diff --git a/test/testGenVectorVc.cxx b/test/testGenVectorVc.cxx
new file mode 100644
index 0000000000000..6eed12b70e773
--- /dev/null
+++ b/test/testGenVectorVc.cxx
@@ -0,0 +1,365 @@
+
+// Vc. Must be before the ROOT includes for std:: math functions to work...
+#include <Vc/Vc>
+
+// ROOT
+#include "Math/GenVector/PositionVector3D.h"
+#include "Math/GenVector/DisplacementVector3D.h"
+#include "Math/GenVector/Plane3D.h"
+#include "Math/GenVector/Transform3D.h"
+#include "TStopwatch.h"
+
+// STL
+#include <random>
+#include <vector>
+#include <iostream>
+#include <string>
+#include <typeinfo>
+#include <cmath>
+
+// note scale here is > 1 as SIMD and scalar floating point calculations not
+// expected to be bit wise identical
+int compare(double v1, double v2, const std::string &name = "", double scale = 1000.0)
+{
+   //  ntest = ntest + 1;
+
+   // numerical double limit for epsilon
+   const double eps     = scale * std::numeric_limits<double>::epsilon();
+   int          iret    = 0;
+   double       delta   = v2 - v1;
+   double       d       = 0;
+   if (delta < 0) delta = -delta;
+   if (v1 == 0 || v2 == 0) {
+      if (delta > eps) {
+         iret = 1;
+      }
+   }
+   // skip case v1 or v2 is infinity
+   else {
+      d = v1;
+
+      if (v1 < 0) d = -d;
+      // add also case when delta is small by default
+      if (delta / d > eps && delta > eps) iret = 1;
+   }
+
+   if (iret != 0) {
+      int pr = std::cout.precision(18);
+      std::cout << "\nDiscrepancy in " << name << "() : " << v1 << " != " << v2 << " discr = " << int(delta / d / eps)
+                << "   (Allowed discrepancy is " << eps << ")\n";
+      std::cout.precision(pr);
+   }
+   return iret;
+}
+
+// randomn generator
+static std::default_random_engine gen;
+// Distributions for each member
+static std::uniform_real_distribution<double> p_x(-800, 800), p_y(-600, 600), p_z(10000, 10500);
+static std::uniform_real_distribution<double> d_x(-0.2, 0.2), d_y(-0.1, 0.1), d_z(0.95, 0.99);
+static std::uniform_real_distribution<double> c_x(3100, 3200), c_y(10, 15), c_z(3200, 3300);
+static std::uniform_real_distribution<double> r_rad(8500, 8600);
+static std::uniform_real_distribution<double> p0(-0.002, 0.002), p1(-0.2, 0.2), p2(0.97, 0.99), p3(-1300, 1300);
+
+template <typename POINT, typename VECTOR, typename PLANE, typename FTYPE>
+class Data {
+public:
+   typedef std::vector<Data, Vc::Allocator<Data>> Vector;
+
+public:
+   POINT  position;
+   VECTOR direction;
+   POINT  CoC;
+   PLANE  plane;
+   FTYPE  radius{0};
+
+public:
+   template <typename INDATA>
+   Data(const INDATA &ind)
+      : position(ind.position.x(), ind.position.y(), ind.position.z()),
+        direction(ind.direction.x(), ind.direction.y(), ind.direction.z()), CoC(ind.CoC.x(), ind.CoC.y(), ind.CoC.z()),
+        plane(ind.plane.A(), ind.plane.B(), ind.plane.C(), ind.plane.D()), radius(ind.radius)
+   {
+   }
+   Data()
+      : position(p_x(gen), p_y(gen), p_z(gen)), direction(d_x(gen), d_y(gen), d_z(gen)),
+        CoC(c_x(gen), c_y(gen), c_z(gen)), plane(p0(gen), p1(gen), p2(gen), p3(gen)), radius(r_rad(gen))
+   {
+   }
+};
+
+template <typename INDATA, typename OUTDATA>
+void fill(const INDATA &in, OUTDATA &out)
+{
+   out.clear();
+   out.reserve(in.size());
+   for (const auto &i : in) {
+      out.emplace_back(i);
+   }
+}
+
+template <typename POINT, typename VECTOR, typename FTYPE>
+inline
+   typename std::enable_if<std::is_arithmetic<typename POINT::Scalar>::value &&
+                              std::is_arithmetic<typename VECTOR::Scalar>::value && std::is_arithmetic<FTYPE>::value,
+                           bool>::type
+   reflectSpherical(POINT &position, VECTOR &direction, const POINT &CoC, const FTYPE radius)
+{
+   constexpr FTYPE zero(0), two(2.0), four(4.0), half(0.5);
+   const FTYPE     a     = direction.Mag2();
+   const VECTOR    delta = position - CoC;
+   const FTYPE     b     = two * direction.Dot(delta);
+   const FTYPE     c     = delta.Mag2() - radius * radius;
+   const FTYPE     discr = b * b - four * a * c;
+   const bool      OK    = discr > zero;
+   if (OK) {
+      const FTYPE dist = half * (std::sqrt(discr) - b) / a;
+      // change position to the intersection point
+      position += dist * direction;
+      // reflect the vector
+      // r = u - 2(u.n)n, r=reflection, u=incident, n=normal
+      const VECTOR normal = position - CoC;
+      direction -= (two * normal.Dot(direction) / normal.Mag2()) * normal;
+   }
+   return OK;
+}
+
+template <typename POINT, typename VECTOR, typename FTYPE>
+inline
+   typename std::enable_if<!std::is_arithmetic<typename POINT::Scalar>::value &&
+                              !std::is_arithmetic<typename VECTOR::Scalar>::value && !std::is_arithmetic<FTYPE>::value,
+                           typename FTYPE::mask_type>::type
+   reflectSpherical(POINT &position, VECTOR &direction, const POINT &CoC, const FTYPE radius)
+{
+   const FTYPE               two(2.0), four(4.0), half(0.5);
+   const FTYPE               a     = direction.Mag2();
+   const VECTOR              delta = position - CoC;
+   const FTYPE               b     = two * direction.Dot(delta);
+   const FTYPE               c     = delta.Mag2() - radius * radius;
+   FTYPE                     discr = b * b - four * a * c;
+   typename FTYPE::mask_type OK    = discr > FTYPE::Zero();
+   if (any_of(OK)) {
+      // Zero out the negative values in discr, to prevent sqrt(-ve)
+      discr(!OK) = FTYPE::Zero();
+      // compute the distance
+      const FTYPE dist = half * (sqrt(discr) - b) / a;
+      // change position to the intersection point
+      position += dist * direction;
+      // reflect the vector
+      // r = u - 2(u.n)n, r=reflection, u=incident, n=normal
+      const VECTOR normal = position - CoC;
+      direction -= (two * normal.Dot(direction) / normal.Mag2()) * normal;
+   }
+   // return the mask indicating which results should be used
+   return OK;
+}
+
+template <typename POINT, typename VECTOR, typename PLANE>
+inline typename std::enable_if<std::is_arithmetic<typename POINT::Scalar>::value &&
+                                  std::is_arithmetic<typename VECTOR::Scalar>::value &&
+                                  std::is_arithmetic<typename PLANE::Scalar>::value,
+                               bool>::type
+reflectPlane(POINT &position, VECTOR &direction, const PLANE &plane)
+{
+   constexpr typename POINT::Scalar two(2.0);
+   const bool                       OK = true;
+   // Plane normal
+   const auto &normal = plane.Normal();
+   // compute distance to the plane
+   const auto scalar   = direction.Dot(normal);
+   const auto distance = -(plane.Distance(position)) / scalar;
+   // change position to reflection point and update direction
+   position += distance * direction;
+   direction -= two * scalar * normal;
+   return OK;
+}
+
+template <typename POINT, typename VECTOR, typename PLANE, typename FTYPE = typename POINT::Scalar>
+inline typename std::enable_if<!std::is_arithmetic<typename POINT::Scalar>::value &&
+                                  !std::is_arithmetic<typename VECTOR::Scalar>::value &&
+                                  !std::is_arithmetic<typename PLANE::Scalar>::value,
+                               typename FTYPE::mask_type>::type
+reflectPlane(POINT &position, VECTOR &direction, const PLANE &plane)
+{
+   const typename POINT::Scalar    two(2.0);
+   const typename FTYPE::mask_type OK(true);
+   // Plane normal
+   const VECTOR normal = plane.Normal();
+   // compute distance to the plane
+   const FTYPE scalar   = direction.Dot(normal);
+   const FTYPE distance = -(plane.Distance(position)) / scalar;
+   // change position to reflection point and update direction
+   position += distance * direction;
+   direction -= two * scalar * normal;
+   return OK;
+}
+
+template <typename T>
+using Point = ROOT::Math::PositionVector3D<ROOT::Math::Cartesian3D<T>, ROOT::Math::DefaultCoordinateSystemTag>;
+template <typename T>
+using Vector = ROOT::Math::DisplacementVector3D<ROOT::Math::Cartesian3D<T>, ROOT::Math::DefaultCoordinateSystemTag>;
+template <typename T>
+using Plane = ROOT::Math::Impl::Plane3D<T>;
+
+int main(int /*argc*/, char ** /*argv*/)
+{
+   int ret = 0;
+
+   {
+
+      const unsigned int nPhotons = 100;
+      std::cout << "Creating " << nPhotons << " random photons ..." << std::endl;
+
+      // Scalar Types
+      Data<Point<double>, Vector<double>, Plane<double>, double>::Vector scalar_data(nPhotons);
+
+      // Vc Types
+      Data<Point<Vc::double_v>, Vector<Vc::double_v>, Plane<Vc::double_v>, Vc::double_v>::Vector vc_data;
+      // Clone the exact random values from the Scalar vector
+      // Note we are making the same number of entries in the container, but each entry is a vector entry
+      // with Vc::double_t::Size entries.
+      fill(scalar_data, vc_data);
+
+      // Loop over the two containers and compare
+      std::cout << "Ray Tracing :-" << std::endl;
+      for (size_t i = 0; i < nPhotons; ++i) {
+         auto &sc = scalar_data[i];
+         auto &vc = vc_data[i];
+
+         // ray tracing
+         reflectSpherical(sc.position, sc.direction, sc.CoC, sc.radius);
+         reflectPlane(sc.position, sc.direction, sc.plane);
+         reflectSpherical(vc.position, vc.direction, vc.CoC, vc.radius);
+         reflectPlane(vc.position, vc.direction, vc.plane);
+
+         std::cout << "Position  " << sc.position << " " << vc.position << std::endl;
+         std::cout << "Direction " << sc.direction << " " << vc.direction << std::endl;
+
+         for (std::size_t j = 0; j < Vc::double_v::Size; ++j) {
+            ret |= compare(sc.position.x(), vc.position.x()[j]);
+            ret |= compare(sc.position.y(), vc.position.y()[j]);
+            ret |= compare(sc.position.z(), vc.position.z()[j]);
+            ret |= compare(sc.direction.x(), vc.direction.x()[j]);
+            ret |= compare(sc.direction.y(), vc.direction.y()[j]);
+            ret |= compare(sc.direction.z(), vc.direction.z()[j]);
+         }
+      }
+
+      // Now test Transformation3D
+      std::cout << "Transforms :-" << std::endl;
+      for (size_t i = 0; i < nPhotons; ++i) {
+         auto &sc = scalar_data[i];
+         auto &vc = vc_data[i];
+
+         // make 6 random scalar Points
+         Point<double> sp1(p_x(gen), p_y(gen), p_z(gen));
+         Point<double> sp2(p_x(gen), p_y(gen), p_z(gen));
+         Point<double> sp3(p_x(gen), p_y(gen), p_z(gen));
+         Point<double> sp4(p_x(gen), p_y(gen), p_z(gen));
+         Point<double> sp5(p_x(gen), p_y(gen), p_z(gen));
+         Point<double> sp6(p_x(gen), p_y(gen), p_z(gen));
+         // clone to Vc versions
+         Point<Vc::double_v> vp1(sp1.x(), sp1.y(), sp1.z());
+         Point<Vc::double_v> vp2(sp2.x(), sp2.y(), sp2.z());
+         Point<Vc::double_v> vp3(sp3.x(), sp3.y(), sp3.z());
+         Point<Vc::double_v> vp4(sp4.x(), sp4.y(), sp4.z());
+         Point<Vc::double_v> vp5(sp5.x(), sp5.y(), sp5.z());
+         Point<Vc::double_v> vp6(sp6.x(), sp6.y(), sp6.z());
+
+         // Make transformations from points
+         // note warnings about axis not having the same angles expected here...
+         // point is to check scalar and vector versions do the same thing
+         ROOT::Math::Impl::Transform3D<double>       st(sp1, sp2, sp3, sp4, sp5, sp6);
+         ROOT::Math::Impl::Transform3D<Vc::double_v> vt(vp1, vp2, vp3, vp4, vp5, vp6);
+
+         // transform the vectors
+         const auto sv = st * sc.direction;
+         const auto vv = vt * vc.direction;
+         std::cout << "Transformed Direction " << sv << " " << vv << std::endl;
+
+         // invert the transformations
+         st.Invert();
+         vt.Invert();
+
+         // Move the points back
+         const auto sv_i = st * sv;
+         const auto vv_i = vt * vv;
+         std::cout << "Transformed Back Direction " << sc.direction << " " << sv_i << " " << vv_i << std::endl;
+
+         for (std::size_t j = 0; j < Vc::double_v::Size; ++j) {
+            ret |= compare(sv.x(), vv.x()[j]);
+            ret |= compare(sv.y(), vv.y()[j]);
+            ret |= compare(sv.z(), vv.z()[j]);
+            ret |= compare(sc.direction.x(), vv_i.x()[j]);
+            ret |= compare(sc.direction.y(), vv_i.y()[j]);
+            ret |= compare(sc.direction.z(), vv_i.z()[j]);
+         }
+
+         ret |= compare(sc.direction.x(), sv_i.x());
+         ret |= compare(sc.direction.y(), sv_i.y());
+         ret |= compare(sc.direction.z(), sv_i.z());
+      }
+   }
+
+   // now run some timing tests
+   {
+      const unsigned int nPhotons = 96000; // Must be multiple of 16 to avoid padding issues below...
+
+      const unsigned int nTests = 1000; // number of tests to run
+
+      // scalar data
+      Data<Point<double>, Vector<double>, Plane<double>, double>::Vector scalar_data(nPhotons);
+      // vector data with total equal number of photons (including vectorised size)
+      Data<Point<Vc::double_v>, Vector<Vc::double_v>, Plane<Vc::double_v>, Vc::double_v>::Vector vc_data(
+         nPhotons / Vc::double_v::Size);
+
+      TStopwatch t;
+
+      double best_time_scalar{9e30}, best_time_vector{9e30};
+
+      // time the scalar implementation
+      for (unsigned int i = 0; i < nTests; ++i) {
+         t.Start();
+         for (auto &sc : scalar_data) {
+            reflectSpherical(sc.position, sc.direction, sc.CoC, sc.radius);
+            reflectPlane(sc.position, sc.direction, sc.plane);
+         }
+         t.Stop();
+         const auto time = t.RealTime();
+         if (time < best_time_scalar) {
+            best_time_scalar = time;
+         }
+      }
+
+      // time the Vc implementation
+      for (unsigned int i = 0; i < nTests; ++i) {
+         t.Start();
+         for (auto &vc : vc_data) {
+            reflectSpherical(vc.position, vc.direction, vc.CoC, vc.radius);
+            reflectPlane(vc.position, vc.direction, vc.plane);
+         }
+         t.Stop();
+         const auto time = t.RealTime();
+         if (time < best_time_vector) {
+            best_time_vector = time;
+         }
+      }
+
+      std::cout << "Scalar best time        = " << best_time_scalar << std::endl;
+      std::cout << "Vectorised Vc best time = " << best_time_vector << std::endl;
+      std::cout << "Vectorised Vc SIMD size = " << Vc::double_v::Size << std::endl;
+      std::cout << "Vectorised Vc speedup   = " << best_time_scalar / best_time_vector << std::endl;
+
+      // assert that the vector time is roughly Vc::double_v::Size times smaller than the scalar time
+      // allow 25% for 'safety'
+      // if (std::fabs((best_time_vector * Vc::double_v::Size) - best_time_scalar) > 0.25 * best_time_scalar) {
+      //   ++ret;
+      // }
+   }
+
+   if (ret)
+      std::cerr << "test FAILED !!! " << std::endl;
+   else
+      std::cout << "test OK " << std::endl;
+   return ret;
+}