From 08b3a9bbf93e6e868a1db6c80854d60863f0f604 Mon Sep 17 00:00:00 2001
From: Lorenzo Moneta <Lorenzo.Moneta@cern.ch>
Date: Thu, 9 Nov 2006 21:25:33 +0000
Subject: [PATCH] add new tests for checking bug 21311 and test on performances
 of point and vector additions

git-svn-id: http://root.cern.ch/svn/root/trunk@16726 27541ba8-7e3a-0410-8455-c3a389f83636
---
 mathcore/test/Makefile          |  20 +++-
 mathcore/test/stress3D.cxx      | 206 ++++++++++++++++++++++++++++++++
 mathcore/test/testBoost.cxx     |  11 +-
 mathcore/test/testGenVector.cxx |  38 +++++-
 mathcore/test/testIterator.cxx  |  62 ++++++++++
 5 files changed, 329 insertions(+), 8 deletions(-)
 create mode 100644 mathcore/test/stress3D.cxx
 create mode 100644 mathcore/test/testIterator.cxx

diff --git a/mathcore/test/Makefile b/mathcore/test/Makefile
index 6dc9defb116..99aaae16461 100644
--- a/mathcore/test/Makefile
+++ b/mathcore/test/Makefile
@@ -40,6 +40,10 @@ BOOSTOBJ     = testBoost.$(ObjSuf)
 BOOSTSRC     = testBoost.$(SrcSuf)
 BOOST        = testBoost$(ExeSuf)
 
+ITERATOROBJ     = testIterator.$(ObjSuf)
+ITERATORSRC     = testIterator.$(SrcSuf)
+ITERATOR        = testIterator$(ExeSuf)
+
 GENVECTOROBJ     = testGenVector.$(ObjSuf)
 GENVECTORSRC     = testGenVector.$(SrcSuf)
 GENVECTOR        = testGenVector$(ExeSuf)
@@ -48,12 +52,17 @@ VECTORIOOBJ     = testVectorIO.$(ObjSuf)
 VECTORIOSRC     = testVectorIO.$(SrcSuf)
 VECTORIO        = testVectorIO$(ExeSuf)
 
+STRESS3DOBJ     = stress3D.$(ObjSuf)
+STRESS3DSRC     = stress3D.$(SrcSuf)
+STRESS3D        = stress3D$(ExeSuf)
+
 
 
-OBJS          = $(COORDINATES3DOBJ) $(COORDINATES4DOBJ) $(ROTATIONOBJ) $(BOOSTOBJ) $(GENVECTOROBJ) $(VECTORIOOBJ)
+OBJS          = $(COORDINATES3DOBJ) $(COORDINATES4DOBJ) $(ROTATIONOBJ) $(BOOSTOBJ) $(GENVECTOROBJ) $(VECTORIOOBJ) $(STRESS3DOBJ) $(ITERATOROBJ)
 
 
-PROGRAMS      = $(COORDINATES3D)  $(COORDINATES4D) $(ROTATION) $(BOOST) $(GENVECTOR) $(VECTORIO)
+PROGRAMS      = $(COORDINATES3D)  $(COORDINATES4D) $(ROTATION) $(BOOST) $(GENVECTOR) $(VECTORIO) $(STRESS3D) $(ITERATOR)
+
 
 		  
 .SUFFIXES: .$(SrcSuf) .$(ObjSuf) $(ExeSuf)
@@ -86,6 +95,13 @@ $(VECTORIO):   	  $(VECTORIOOBJ)
 		    $(LD) $(LDFLAGS) $^ $(LIBS) $(EXTRALIBS) $(EXTRAIOLIBS) $(OutPutOpt)$@
 		    @echo "$@ done"
 
+$(STRESS3D):      $(STRESS3DOBJ)
+		    $(LD) $(LDFLAGS) $^ $(LIBS) $(EXTRALIBS) $(OutPutOpt)$@
+		    @echo "$@ done"
+
+$(ITERATOR):           $(ITERATOROBJ)
+		    $(LD) $(LDFLAGS) $^ $(LIBS) $(EXTRALIBS) $(OutPutOpt)$@
+		    @echo "$@ done"
 
 
 clean:
diff --git a/mathcore/test/stress3D.cxx b/mathcore/test/stress3D.cxx
new file mode 100644
index 00000000000..3992f1542a8
--- /dev/null
+++ b/mathcore/test/stress3D.cxx
@@ -0,0 +1,206 @@
+#include "Math/Vector3D.h"
+#include "Math/Point3D.h"
+#include "TVector3.h"
+
+#include "TStopwatch.h"
+#include "TRandom3.h"
+
+
+class VectorTest { 
+
+private: 
+
+  size_t n2Loop ;
+  size_t nGen;
+
+// global data variables 
+  std::vector<double> dataX; 
+  std::vector<double> dataY;  
+  std::vector<double> dataZ;  
+
+
+
+public: 
+  
+  VectorTest(int n1, int n2) : 
+    n2Loop(n1),
+    nGen(n2),
+    dataX(std::vector<double>(n2)), 
+    dataY(std::vector<double>(n2)), 
+    dataZ(std::vector<double>(n2)) 
+  {}
+    
+
+  
+
+  void print(TStopwatch & time, std::string s) { 
+    int pr = std::cout.precision(8);
+    std::cout << s << "\t" << " time = " << time.RealTime() << "\t(sec)\t" 
+      //    << time.CpuTime() 
+	      << std::endl;
+    std::cout.precision(pr);
+  }
+
+    
+  int check(std::string name, double s1, double s2, double s3, double scale=1) {
+    double eps = 10*scale*std::numeric_limits<double>::epsilon();
+    if (  std::fabs(s1-s2) < eps*std::fabs(s1) && std::fabs(s1-s3)  < eps*std::fabs(s1) ) return 0; 
+    std::cout.precision(16);
+    std::cout << s1 << "\t" << s2 <<"\t" << s3 << "\n";
+    std::cout << "Test " << name << " failed !!\n\n"; 
+    return -1; 
+  }
+
+
+  void genData() { 
+    int n = nGen;
+
+  // generate n -4 momentum quantities 
+    TRandom3 r;
+    for (int i = 0; i < n ; ++ i) { 
+ 
+      double phi = r.Rndm()*3.1415926535897931; 
+      double eta = r.Uniform(-5.,5.); 
+      double pt = r.Exp(10.);
+    
+    // fill vectors 
+    
+      ROOT::Math::RhoEtaPhiVector q( pt, eta, phi); 
+      dataX[i] =  q.x(); 
+      dataY[i] =  q.y(); 
+      dataZ[i] =  q.z(); 
+
+    }
+  }
+
+
+
+  template <class V> 
+  void testCreate( std::vector<V *> & dataV, TStopwatch & tim, double& t,  std::string s) { 
+    
+    int n = dataX.size(); 
+    dataV.resize(n); 
+    tim.Start();
+    for (int i = 0; i < n; ++i) { 
+      dataV[i] = new V( dataX[i], dataY[i], dataZ[i] ); 
+    }
+    tim.Stop();
+    t += tim.RealTime();
+    print(tim,s);
+  }
+
+
+template <class V>
+double testVectorAddition( const std::vector<V *> & dataV, TStopwatch & tim, double& t,  std::string s) { 
+  unsigned int n = std::min(n2Loop, dataV.size() );
+  tim.Start(); 
+  V  vSum = *(dataV[0]);
+  for (unsigned int i = 1; i < n; ++i) { 
+    vSum += *(dataV[i]); 
+  }
+  tim.Stop();
+  print(tim,s);
+  t += tim.RealTime();
+  return vSum.Mag2(); 
+}  
+
+template <class P>
+double testPointAddition( const std::vector<P *> & dataP, TStopwatch & tim, double& t,  std::string s) { 
+  unsigned int n = std::min(n2Loop, dataP.size() );
+  tim.Start(); 
+  P  pSum = *(dataP[0]);
+  for (unsigned int i = 1; i < n; ++i) { 
+     P & p2 = *(dataP[i]); 
+#ifndef HEAP_CREATION
+     pSum += ROOT::Math::XYZVector(p2);
+#else
+     ROOT::Math::XYZVector * v2 = new ROOT::Math::XYZVector(p2);
+     pSum += *v2;
+#endif
+  }
+  tim.Stop();
+  print(tim,s);
+  t += tim.RealTime();
+  return pSum.Mag2(); 
+}
+
+};
+
+
+int main(int argc,const char *argv[]) { 
+
+  int ngen = 1000000;
+  if (argc > 1)  ngen = atoi(argv[1]);
+  int nloop2 = ngen;
+  if (argc > 2)  nloop2 = atoi(argv[1]);
+
+
+  TStopwatch t;
+
+  VectorTest a(ngen,nloop2);
+
+  a.genData(); 
+
+
+  double t1 = 0;
+  double t2 = 0;
+  double t3 = 0;
+
+  std::vector<TVector3 *> v1;
+  std::vector<ROOT::Math::XYZVector *> v2;
+  std::vector<ROOT::Math::XYZPoint *> v3;
+
+  double s1,s2,s3;
+
+  a.genData(); 
+  a.testCreate     (v2, t, t2,    "creation XYZVector     " );  
+
+  a.genData(); 
+  a.testCreate     (v3, t, t3,    "creation XYZPoint      " ); 
+
+  a.genData(); 
+  a.testCreate     (v1, t, t1,    "creation TVector3      " ); 
+
+  std::cout << "\n";
+
+#ifdef MORETEST
+  t1 = 0; 
+  t2 = 0; 
+  t3 = 0; 
+#endif
+
+  s1=a.testVectorAddition   (v1, t, t1, "Addition TVector3      " );  
+  s2=a.testVectorAddition   (v2, t, t2, "Addition XYZVector     "  ); 
+  s3=a.testPointAddition    (v3, t, t3, "Addition XYZPoint      " );       
+
+  a.check("Addition",s1,s2,s3);
+
+#ifdef MORETEST
+
+  s2=a.testVectorAddition   (v2, t, t2, "Addition XYZVector     "  ); 
+  s1=a.testVectorAddition   (v1, t, t1, "Addition TVector3      " );  
+  s3=a.testPointAddition    (v3, t, t3, "Addition XYZPoint      " );       
+
+  s2=a.testVectorAddition   (v2, t, t2, "Addition XYZVector     "  ); 
+  s3=a.testPointAddition    (v3, t, t3, "Addition XYZPoint      " );       
+  s1=a.testVectorAddition   (v1, t, t1, "Addition TVector3      " );  
+
+  s1=a.testVectorAddition   (v1, t, t1, "Addition TVector3      " );  
+  s3=a.testPointAddition    (v3, t, t3, "Addition XYZPoint      " );       
+  s2=a.testVectorAddition   (v2, t, t2, "Addition XYZVector     "  ); 
+
+  s3=a.testPointAddition    (v3, t, t3, "Addition XYZPoint      " );       
+  s2=a.testVectorAddition   (v2, t, t2, "Addition XYZVector     "  ); 
+  s1=a.testVectorAddition   (v1, t, t1, "Addition TVector3      " );  
+
+  s3=a.testPointAddition    (v3, t, t3, "Addition XYZPoint      " );       
+  s1=a.testVectorAddition   (v1, t, t1, "Addition TVector3      " );  
+  s2=a.testVectorAddition   (v2, t, t2, "Addition XYZVector     "  ); 
+
+#endif
+
+  std::cout << "Total Time for  TVector3        = " << t1 << "\t(sec)" << std::endl;
+  std::cout << "Total Time for  XYZVector       = " << t2 << "\t(sec)" << std::endl;
+  std::cout << "Total Time for  XYZPoint        = " << t3 << "\t(sec)" << std::endl;
+
+}
diff --git a/mathcore/test/testBoost.cxx b/mathcore/test/testBoost.cxx
index d6e032466c9..b6d3ddded53 100644
--- a/mathcore/test/testBoost.cxx
+++ b/mathcore/test/testBoost.cxx
@@ -4,6 +4,8 @@
 #include "Math/Vector3D.h"
 #include "Math/Vector4D.h"
 
+#include <iterator> 
+
 using namespace ROOT::Math;
 
 int main() { 
@@ -28,9 +30,12 @@ int main() {
   Polar3DVector bv(0.99999,1.,2);
   std::cout << "BoostVector " << XYZVector(bv) << " beta boost = " << XYZVector(bv).R() << std::endl; 
   Boost b(bv);
-  double d[3];
-  b.GetComponents(d,d+3);
-  std::cout << "Boost Components" << d[0] << " " << d[1] << "  " << d[2] << std::endl;
+  std::cout << "Boost Components : "; 
+  std::ostream_iterator<double> oi(std::cout,"\t");
+  b.GetComponents(oi);
+  std::cout << std::endl;
+  
+
 
 
   vb1 = b(v);
diff --git a/mathcore/test/testGenVector.cxx b/mathcore/test/testGenVector.cxx
index e435409ad58..288fec85708 100644
--- a/mathcore/test/testGenVector.cxx
+++ b/mathcore/test/testGenVector.cxx
@@ -1,3 +1,4 @@
+
 #include "Math/Vector3D.h"
 #include "Math/Point3D.h"
 #include "Math/EulerAngles.h"
@@ -14,8 +15,6 @@
 
 #include "Math/VectorUtil.h"
 
-#include <iostream>
-
 using namespace ROOT::Math;
 using namespace ROOT::Math::VectorUtil;
 
@@ -176,6 +175,13 @@ int testPoint3D() {
   vg4 = pg - pl;
 #endif
 
+  // operator - 
+  XYZPoint q1(1.,2.,3.);
+  XYZPoint q2 = -1.* q1; 
+  XYZVector v2 = -XYZVector(q1);
+  iret |= compare(XYZVector(q2) == v2,true,"reflection");
+
+
   if (iret == 0) std::cout << "\t\t\t\t\tOK\n"; 
   else std::cout << "\t\t\t\tFAILED\n"; 
 
@@ -264,7 +270,14 @@ int testRotations3D() {
   bool comp = (rotInv == rot ); 
   iret |= compare(comp,true,"inversion");
 
-    
+  // rotation and scaling of points
+  XYZPoint q1(1.,2,3); double a = 3;
+  XYZPoint qr1 =  rot( a * q1);
+  XYZPoint qr2 =  a * rot( q1);
+  iret |= compare(qr1.X(), qr2.X(),"x diff",10 );
+  iret |= compare(qr1.Y(), qr2.Y(),"y diff",10 );
+  iret |= compare(qr1.Z(), qr2.Z(),"z diff",10 );
+
 
   if (iret == 0) std::cout << "\tOK\n"; 
   else std::cout << "\t FAILED\n";
@@ -337,6 +350,25 @@ int testTransform3D() {
   iret |= compare(v3.Y(), v2.Y(),"y diff",10 );
   iret |= compare(v3.Z(), v2.Z(),"z diff",10 );
 
+  XYZPoint q1(1,2,3);
+  XYZPoint q2(-1,-2,-3);
+  XYZPoint q3 = q1 +  XYZVector(q2);
+  //std::cout << q3 << std::endl; 
+  XYZPoint qt3 = t3(q3); 
+  //std::cout << qt3 << std::endl; 
+  XYZPoint qt1 = t3(q1);
+  XYZVector vt2 = t3( XYZVector(q2) );
+  XYZPoint qt4 = qt1 + vt2; 
+  iret |= compare(qt3.X(), qt4.X(),"x diff",10 );
+  iret |= compare(qt3.Y(), qt4.Y(),"y diff",10 );
+  iret |= compare(qt3.Z(),  qt4.Z(),"z diff",10 );
+     //std::cout << qt4 << std::endl; 
+
+  // this fails
+//  double a = 3;
+  //XYZPoint q4 = a*q1;
+//   std::cout << t3( a * q1) << std::endl;
+//   std::cout << a * t3(q1) << std::endl;
 
 
   if (iret == 0) std::cout << "\t\t\t\tOK\n"; 
diff --git a/mathcore/test/testIterator.cxx b/mathcore/test/testIterator.cxx
new file mode 100644
index 00000000000..21fdf52fc16
--- /dev/null
+++ b/mathcore/test/testIterator.cxx
@@ -0,0 +1,62 @@
+
+#include "Math/Vector3D.h"
+#include "Math/Point3D.h"
+#include "Math/Vector4D.h"
+#include "Math/EulerAngles.h"
+
+#include "Math/Transform3D.h"
+#include "Math/LorentzRotation.h"
+#include "Math/Boost.h"
+
+#include "Math/Rotation3D.h"
+#include "Math/RotationX.h"
+#include "Math/RotationY.h"
+#include "Math/RotationZ.h"
+#include "Math/Quaternion.h"
+#include "Math/AxisAngle.h"
+#include "Math/EulerAngles.h"
+
+#include "Math/VectorUtil.h"
+
+#include <iostream>
+#include <iterator> 
+#include <vector>
+
+using namespace ROOT::Math;
+using namespace ROOT::Math::VectorUtil;
+
+template <class Rot> 
+void printRot(const Rot & rot) {
+   std::cout << "rot:  ( " ;
+   std::ostream_iterator<double> oi(std::cout,"  ");
+   rot.GetComponents(oi); 
+   std::cout << ") " << std::endl;
+}   
+template<class V> 
+void printVec(const V & v ) {
+   std::cout << "vec : ( " ;
+   std::ostream_iterator<double> oi(std::cout,"  ");
+   v.GetCoordinates(oi); 
+   std::cout << ") " << std::endl;
+}   
+int main() { 
+
+   std::vector<double> data(16); 
+   XYZVector v(1.,2.,3);
+   printVec(v); 
+   XYZPoint p(v); printVec(p);
+   XYZTVector q(1.,2,3,4); printVec(q); 
+
+   AxisAngle ar(v,3.); printRot(ar); 
+   EulerAngles er(ar); printRot(er); 
+   Quaternion qr(er); printRot(qr); 
+   Rotation3D rr(qr) ; printRot(rr); 
+ 
+   Transform3D t(rr,v); printRot(t); 
+
+   
+   Boost b(0.9*v/(v.R())); printRot(b); 
+   LorentzRotation lr(rr); printRot(lr);
+   LorentzRotation lr2(b); printRot(lr2);
+
+}
-- 
GitLab