used Vector::Scalar instead of double in utility functions for compiling with Vc types

git-svn-id: http://root.cern.ch/svn/root/branches/dev/rootvc@48449 27541ba8-7e3a-0410-8455-c3a389f83636

math/genvector/inc/Math/GenVector/VectorUtil.h

-// @(#)root/mathcore:$Id$
+// @(#)root/mathcore:$Id: 9ef2a4a7bd1b62c1293920c2af2f64791c75bdd8 $
 // Authors: W. Brown, M. Fischler, L. Moneta    2005  
 
 
@@ -58,8 +58,8 @@ namespace ROOT {
 	 \f[ \Delta \phi = \phi_2 - \phi_1 \f]
       */
       template <class Vector1, class Vector2> 
-      inline double DeltaPhi( const Vector1 & v1, const Vector2 & v2) { 
-	double dphi = v2.Phi() - v1.Phi(); 
+      inline typename Vector1::Scalar DeltaPhi( const Vector1 & v1, const Vector2 & v2) { 
+	typename Vector1::Scalar dphi = v2.Phi() - v1.Phi(); 
 	if ( dphi > M_PI ) {
 	  dphi -= 2.0*M_PI;
 	} else if ( dphi <= -M_PI ) {
@@ -79,9 +79,9 @@ namespace ROOT {
        \f[ \Delta R2 = ( \Delta \phi )^2 + ( \Delta \eta )^2  \f]
     */ 
       template <class Vector1, class Vector2> 
-      inline double DeltaR2( const Vector1 & v1, const Vector2 & v2) { 
-	double dphi = DeltaPhi(v1,v2); 
-	double deta = v2.Eta() - v1.Eta(); 
+      inline typename Vector1::Scalar DeltaR2( const Vector1 & v1, const Vector2 & v2) { 
+	typename Vector1::Scalar dphi = DeltaPhi(v1,v2); 
+	typename Vector1::Scalar deta = v2.Eta() - v1.Eta(); 
 	return dphi*dphi + deta*deta; 
       }
 
@@ -92,9 +92,9 @@ namespace ROOT {
        \param v2  Vector 2
        \return   Angle between the two vectors
        \f[ \Delta R = \sqrt{  ( \Delta \phi )^2 + ( \Delta \eta )^2 } \f]
-    */ 
+    */
       template <class Vector1, class Vector2> 
-      inline double DeltaR( const Vector1 & v1, const Vector2 & v2) { 
+      inline typename Vector1::Scalar DeltaR( const Vector1 & v1, const Vector2 & v2) { 
 	return std::sqrt( DeltaR2(v1,v2) ); 
       }
 
@@ -214,18 +214,33 @@ namespace ROOT {
        \f[ M_{12} = \sqrt{ (\vec{v1} + \vec{v2} ) \cdot (\vec{v1} + \vec{v2} ) } \f]
     */ 
       template <class Vector1, class Vector2> 
-      inline double InvariantMass( const Vector1 & v1, const Vector2 & v2) { 
-	double ee = (v1.E() + v2.E() );
-	double xx = (v1.X() + v2.X() );
-	double yy = (v1.Y() + v2.Y() );
-	double zz = (v1.Z() + v2.Z() );
-	double mm2 = ee*ee - xx*xx - yy*yy - zz*zz; 
+      inline typename Vector1::Scalar InvariantMass( const Vector1 & v1, const Vector2 & v2) { 
+         typedef typename  Vector1::Scalar Scalar;
+         Scalar ee = (v1.E() + v2.E() );
+         Scalar xx = (v1.X() + v2.X() );
+         Scalar yy = (v1.Y() + v2.Y() );
+         Scalar zz = (v1.Z() + v2.Z() );
+         Scalar mm2 = ee*ee - xx*xx - yy*yy - zz*zz; 
 	return mm2 < 0.0 ? -std::sqrt(-mm2) : std::sqrt(mm2);
 // 	PxPyPzE4D<double> q(xx,yy,zz,ee); 
 // 	return q.M();
 	//return ( v1 + v2).mag(); 
       }
 
+      template <class Vector1, class Vector2> 
+      inline typename Vector1::Scalar InvariantMass2( const Vector1 & v1, const Vector2 & v2) { 
+         typedef typename  Vector1::Scalar Scalar;
+         Scalar ee = (v1.E() + v2.E() );
+         Scalar xx = (v1.X() + v2.X() );
+         Scalar yy = (v1.Y() + v2.Y() );
+         Scalar zz = (v1.Z() + v2.Z() );
+         Scalar mm2 = ee*ee - xx*xx - yy*yy - zz*zz; 
+         return mm2 ; // < 0.0 ? -std::sqrt(-mm2) : std::sqrt(mm2);
+// 	PxPyPzE4D<double> q(xx,yy,zz,ee); 
+// 	return q.M();
+	//return ( v1 + v2).mag(); 
+      }
+
       // rotation and transformations
       
       
@@ -339,15 +354,16 @@ namespace ROOT {
 	  X(), Y(), Z(), T()  and SetXYZT methods.
 	  The beta of the boost must be <= 1 or a nul Lorentz Vector will be returned	  	  
       */
-      template <class LVector> 
-      LVector boostX(const LVector & v, double beta) { 
+       template <class LVector, class T> 
+      LVector boostX(const LVector & v, T beta) { 
          if (beta >= 1) {
             GenVector::Throw ("Beta Vector supplied to set Boost represents speed >= c");
             return LVector();
          }    
-	double gamma = 1.0/ std::sqrt(1.0 - beta*beta); 
-	double x2 = gamma * v.X() + gamma * beta * v.T();
-	double t2 = gamma * beta * v.X() + gamma * v.T(); 
+	T gamma = 1.0/ std::sqrt(1.0 - beta*beta); 
+        typename LVector::Scalar x2 = gamma * v.X() + gamma * beta * v.T();
+        typename LVector::Scalar t2 = gamma * beta * v.X() + gamma * v.T(); 
+
 	LVector lv; 
 	lv.SetXYZT(x2,v.Y(),v.Z(),t2);
 	return lv;