rotationApplication.cxx

Go to the documentation of this file.

/**********************************************************************
 *                                                                    *
 * Copyright (c) 2005 , FNAL LCG ROOT MathLib Team                    *
 *                                                                    *
 *                                                                    *
 **********************************************************************/

// RotationApplication.cpp
//
// Created by:  M. Fischler, Aug 10, 2005
//
// Tests that each Rotation produces correct results when applied to
// each form of vector in each oordinate system, and incidently that products
// of rotations work properly.
//
// The strategy is to build up sequences of rotations about the X, Y, and Z
// axes, such that we can easily determine the correct vector answer.
//
// =================================================================

#include "Math/GenVector/DisplacementVector3D.h"
#include "Math/GenVector/PositionVector3D.h"
#include "Math/GenVector/Cartesian3D.h"
#include "Math/GenVector/Polar3D.h"
#include "Math/GenVector/CylindricalEta3D.h"
#include "Math/GenVector/LorentzVector.h"
#include "Math/GenVector/Rotation3D.h"
#include "Math/GenVector/AxisAngle.h"
#include "Math/GenVector/EulerAngles.h"
#include "Math/GenVector/Quaternion.h"
#include "Math/GenVector/RotationZYX.h"
#include "Math/GenVector/RotationX.h"
#include "Math/GenVector/RotationY.h"
#include "Math/GenVector/RotationZ.h"

#include "Math/GenVector/VectorUtil.h"

#include "Math/Vector3Dfwd.h"

#include "CoordinateTraits.h"
#include "RotationTraits.h"

#include <iostream>
#include <limits>
#include <cmath>
#include <vector>

using std::sin;
using std::cos;

//#define TRACE2
//#define TRACE1
//#define DEBUG

using namespace ROOT::Math;

#ifndef __CINT__

template <typename T1, typename T2 >
struct Precision {
   enum { result = std::numeric_limits<T1>::digits <= std::numeric_limits<T2>::digits   };
};

template <typename T1, typename T2, bool>
struct LessPreciseType {
   typedef T1 type;
};
template <typename T1, typename T2>
struct LessPreciseType<T1, T2, false> {
   typedef T2 type;
};


const double machEpsilon = pow (2.0, -52.0);

template <typename Scalar1, typename Scalar2>
int
closeEnough (Scalar1 s1, Scalar2 s2, std::string const & coord, double ticks) {
   int ret = 0;
   int pr = std::cout.precision(18);
   Scalar1 eps1 = std::numeric_limits<Scalar1>::epsilon();
   Scalar2 eps2 = std::numeric_limits<Scalar2>::epsilon();
   typedef typename LessPreciseType<Scalar1, Scalar2,Precision<Scalar1,Scalar2>::result>::type Scalar;
   Scalar epsilon = (eps1 >= eps2) ? eps1 : eps2;
   Scalar ss1 (s1);
   Scalar ss2 (s2);
   Scalar diff = ss1 - ss2;
   if (diff < 0) diff = -diff;
   if ( diff > ticks*epsilon ) {
      ret=3;
      std::cout << "\n\n????????\n\nAbsolute discrepancy in " << coord << "(): "
      << ss1 << " != " << ss2 << "\n"
      << "   (Allowed discrepancy is " << ticks
      << " ticks = " << ticks*epsilon
      << ")\nDifference is " << diff/epsilon << " ticks\n";
   }
   std::cout.precision (pr);
   return ret;
}

template <class V1, class V2>
int compare3D (const V1 & v1, const V2 & v2, double ticks) {
   int ret =0;
   typedef typename V1::CoordinateType CoordType1;
   typedef typename V2::CoordinateType CoordType2;

   ret |= closeEnough ( v1.x(),     v2.x(),     "x"     ,ticks);
   ret |= closeEnough ( v1.y(),     v2.y(),     "y"     ,ticks);
   ret |= closeEnough ( v1.z(),     v2.z(),     "z"     ,ticks);

   if (ret != 0) {
      std::cout << "Discrepancy detected (see above) is between:\n  "
      << CoordinateTraits<CoordType1>::name() << " and\n  "
      << CoordinateTraits<CoordType2>::name() << "\n"
      << "with v = (" << v1.x() << ", " << v1.y() << ", " << v1.z()
      << ")\nv1 is " << v1
      << "\nv2 is " << v2 << "\n\n";
   }

   return ret;
}

template <class V> struct correctedTicks {
   double operator()(double ticks, const V& /*v */ , const XYZVector & /* ans */)
   {
      return ticks;
   }
};

double correctTicks (double ticks,double z,double r) {
   double e = ticks*std::fabs( z*z / (r*r-z*z) );
   if (e < ticks) return ticks;
   return e;
}

template<> struct
correctedTicks< DisplacementVector3D < CylindricalEta3D<double> > > {
   double operator()(double ticks,
                     const DisplacementVector3D< CylindricalEta3D<double> >& v,
                     const XYZVector & ans) {
      double t1 = correctTicks (ticks,   v.z(),   v.r());
      double t2 = correctTicks (ticks, ans.z(), ans.r());
      return t1 > t2 ? t1 : t2;
   }
};

template<> struct
correctedTicks< PositionVector3D < CylindricalEta3D<double> > > {
   double operator()(double ticks,
                     const PositionVector3D< CylindricalEta3D<double> >& v,
                     const XYZVector & ans) {
      double t1 = correctTicks (ticks,   v.z(),   v.r());
      double t2 = correctTicks (ticks, ans.z(), ans.r());
      return t1 > t2 ? t1 : t2;
   }
};

template <class R, class V>
int testApplication(const R& r,const V& v,const XYZVector &answer,double t) {

   typedef typename V::CoordinateType CoordType;

   int ret = 0;
   correctedTicks<V> ct;
   double ticks = ct(t, v, answer);
#ifdef DEBUG
   std::cout <<">>>>> Testing application of "
   << RotationTraits<R>::name() << " " << r << "\non "
   << CoordinateTraits<typename V::CoordinateType>::name() <<
   v << " ticks = " << ticks;
#endif
#ifdef TRACE2
   std::cout << "  about to do V rv = r(v) - \n";
#endif
   V rv = r(v);
#ifdef TRACE2
   std::cout << "ok ";
#endif
   // comparison here should be == but need to use 10 ticks to be sure for 32 bits machines
   // when results are flushed in memory and approximated
   if ( compare3D(rv, r*v, 10 ) != 0) {
      std::cout << "Inconsistency between R(v) and R*v for R = "
      << RotationTraits<R>::name() << " " << r
      << "\nacting on "  << CoordinateTraits<CoordType>::name() << v << "\n";
      ret |= 9;
   }
#ifdef TRACE2
   std::cout << "+ also did rv != r*v ";
#endif
   if ( closeEnough(v.r(), rv.r(), "r", ticks) != 0 ) {
      std::cout << "Radius change  between R(v) and R*v for R = "
      << RotationTraits<R>::name() << " " << r
      << "\nacting on "
      << CoordinateTraits<CoordType>::name()
      << v << "\n";
      ret |= 17;
   }
#ifdef TRACE2
   std::cout << "\n---- about to do compare3D ----";
#endif
   ret |= compare3D (rv, answer, ticks);
#ifdef TRACE2
   std::cout << " done \n";
#endif
#ifdef DEBUG
   if (ret == 0) std::cout << " OK\n";
#endif
   return ret;
}


XYZVector rxv ( double phi, XYZVector v ) {
   double c = cos(phi);
   double s = sin(phi);
   return XYZVector ( v.x(), c*v.y()-s*v.z(), s*v.y()+c*v.z() );
}

XYZVector ryv ( double phi, XYZVector v ) {
   double c = cos(phi);
   double s = sin(phi);
   return XYZVector ( c*v.x()+s*v.z(), v.y(), -s*v.x()+c*v.z() );
}

XYZVector rzv ( double phi, XYZVector v ) {
   double c = cos(phi);
   double s = sin(phi);
   return XYZVector ( c*v.x()-s*v.y(), s*v.x()+c*v.y(), v.z() );
}

enum XYZ { X, Y, Z } ;

struct TestRotation {
   std::vector<XYZ> xyz;
   std::vector<double> phi;
   TestRotation (std::vector<XYZ> const & xyz_, std::vector<double> const & phi_)
   : xyz(xyz_), phi(phi_) {}
};


Rotation3D rrr (TestRotation const & t) {
#ifdef TRACE1
   std::cout << "---- rrr ----";
#endif
   Rotation3D r;
   for (unsigned int i=0; i<t.xyz.size(); ++i) {
      switch ( t.xyz[i] ) {
         case X:
         r = r*RotationX(t.phi[i]);
         break;
         case Y:
         r = r*RotationY(t.phi[i]);
         break;
         case Z:
         r = r*RotationZ(t.phi[i]);
         break;
      }
   }
#ifdef TRACE1
   std::cout << " done\n";
#endif
   return r;
}

XYZVector ans ( TestRotation const & t, XYZVector const & v_in) {
   XYZVector v(v_in);
   for (int i=t.xyz.size()-1; i >= 0; --i) {
      switch ( t.xyz[i] ) {
         case X:
         v = rxv ( t.phi[i],v );
         break;
         case Y:
         v = ryv ( t.phi[i],v );
         break;
         case Z:
         v = rzv ( t.phi[i],v );
         break;
      }
   }
   return v;
}

const double pi = 3.1415926535897932385;

std::vector<TestRotation>
makeTestRotations () {
#ifdef TRACE1
   std::cout << "---- makeTestRotations ----";
#endif
   std::vector<TestRotation> t;
   std::vector<XYZ>    xyz;
   std::vector<double> phi;



   xyz.clear();      phi.clear();
   xyz.push_back(X); phi.push_back( pi/2 );
   t.push_back(TestRotation(xyz,phi));

   xyz.clear();      phi.clear();
   xyz.push_back(Y); phi.push_back( pi/2 );
   t.push_back(TestRotation(xyz,phi));

   xyz.clear();      phi.clear();
   xyz.push_back(Z); phi.push_back( pi/2 );
   t.push_back(TestRotation(xyz,phi));

   xyz.clear();      phi.clear();
   xyz.push_back(X); phi.push_back( -pi/6 );
   t.push_back(TestRotation(xyz,phi));

   xyz.clear();      phi.clear();
   xyz.push_back(Y); phi.push_back( pi/6 );
   t.push_back(TestRotation(xyz,phi));

   xyz.clear();      phi.clear();
   xyz.push_back(Z); phi.push_back( pi/3 );
   t.push_back(TestRotation(xyz,phi));

   xyz.clear();      phi.clear();
   xyz.push_back(X); phi.push_back( -pi/6 );
   xyz.push_back(Y); phi.push_back(  pi/3 );
   t.push_back(TestRotation(xyz,phi));

   xyz.clear();      phi.clear();
   xyz.push_back(X); phi.push_back( -pi/6 );
   xyz.push_back(Y); phi.push_back(  pi/4 );
   xyz.push_back(Z); phi.push_back( -pi/5 );
   t.push_back(TestRotation(xyz,phi));

   xyz.clear();      phi.clear();
   xyz.push_back(Y); phi.push_back( pi );
   xyz.push_back(X); phi.push_back( -pi/2 );
   xyz.push_back(Z); phi.push_back( -pi/1.5 );
   xyz.push_back(Y); phi.push_back( -pi/3 );
   t.push_back(TestRotation(xyz,phi));

   xyz.clear();      phi.clear();
   xyz.push_back(Z); phi.push_back(  1.3 );
   xyz.push_back(Y); phi.push_back( -1.1 );
   xyz.push_back(X); phi.push_back(  0.4 );
   xyz.push_back(Y); phi.push_back(  0.7 );
   t.push_back(TestRotation(xyz,phi));

   xyz.clear();      phi.clear();
   xyz.push_back(X); phi.push_back(  1.3 );
   xyz.push_back(Z); phi.push_back( -1.1 );
   xyz.push_back(Y); phi.push_back(  0.4 );
   xyz.push_back(Z); phi.push_back(  0.7 );
   t.push_back(TestRotation(xyz,phi));

   xyz.clear();      phi.clear();
   xyz.push_back(Y); phi.push_back(  1.3 );
   xyz.push_back(X); phi.push_back( -1.1 );
   xyz.push_back(Z); phi.push_back(  0.4 );
   xyz.push_back(X); phi.push_back(  0.7 );
   t.push_back(TestRotation(xyz,phi));

   xyz.clear();      phi.clear();
   xyz.push_back(Z); phi.push_back(  .03 );
   xyz.push_back(Y); phi.push_back( -.05 );
   xyz.push_back(X); phi.push_back(  0.04 );
   xyz.push_back(Y); phi.push_back(  0.07 );
   xyz.push_back(Z); phi.push_back( -0.02 );
   t.push_back(TestRotation(xyz,phi));

#ifdef TRACE1
   std::cout << " done\n";
#endif
   return t;
}

std::vector<XYZVector> makeTestVectors () {
#ifdef TRACE1
   std::cout << "---- makeTestVectors ----";
#endif
   std::vector<XYZVector> vs;
   vs.push_back(XYZVector ( 1, 0, 0 ));
   vs.push_back(XYZVector ( 0, 1, 0 ));
   vs.push_back(XYZVector ( 0, 0, 1 ));
   vs.push_back(XYZVector ( -1, 0, 0 ));
   vs.push_back(XYZVector ( 0, -1, 0 ));
   vs.push_back(XYZVector ( 0, 0, -1 ));
   vs.push_back(XYZVector ( 1, 2, 3 ));
   vs.push_back(XYZVector ( 2, -1, 3 ));
   vs.push_back(XYZVector ( -3, 1, -2 ));
   vs.push_back(XYZVector ( 0, .00001, -2 ));

#ifdef TRACE1
   std::cout << " done\n";
#endif
   return vs;
}

template <class R, class V>
int doTest (TestRotation const & testRotation, XYZVector const & testVector,
            double ticks) {
#ifdef TRACE1
   std::cout << "---- doTest ----";
#endif
   int ret = 0;
   R r ( rrr(testRotation) );
   V v(testVector);
   XYZVector rv = ans (testRotation, testVector);
   ret |= testApplication (r, v, rv, ticks);
#ifdef TRACE1
   std::cout << " done\n";
#endif

   if (ret == 0) std::cout << ".";

   return ret;
}

template <class R, class C>
int doTestL (TestRotation const & testRotation, XYZVector const & testVector,
             double ticks) {
#ifdef TRACE1
   std::cout << "---- doTestL ----";
#endif
   int ret = 0;
   R r ( rrr(testRotation) );
   LorentzVector<C> v;
   double x = testVector.X();
   double y = testVector.Y();
   double z = testVector.Z();
   double t = std::sqrt (x*x + y*y + z*z + 1);
   v.SetXYZT ( x, y, z, t );
   XYZVector rv = ans (testRotation, testVector);
   ret |= testApplication (r, v, rv, ticks);
   LorentzVector<C> rv2 = r(v);
   ret |= closeEnough (t, rv2.E(), "t", ticks);
#ifdef TRACE1
   std::cout << " done\n";
#endif
   return ret;
}

struct ForeignVector {
   typedef Cartesian3D<> CoordinateType;
   XYZVector v;
   template <class V>
   explicit ForeignVector (V const & v_) : v(v_) {}
   ForeignVector (double xx, double yy, double zz) : v(xx,yy,zz) {}
   double x() const { return v.x(); }
   double y() const { return v.y(); }
   double z() const { return v.z(); }
   double r() const { return v.r(); }
   bool operator==(ForeignVector const & rhs) {return v == rhs.v;}
   bool operator!=(ForeignVector const & rhs) {return v != rhs.v;}
};
std::ostream & operator<< (std::ostream& os, const ForeignVector& v) {
   return os << v.v;
}


template <class R>
int doTestOfR (TestRotation const & testRotation, XYZVector const & testVector){
#ifdef TRACE1
   std::cout << "---- doTestofR ----\n";
#endif
   int ret = 0;
   const double ticks = 100;  // move from 32 to 100
#ifdef DEBUG
   std::cout << ">>>>> DisplacementVector3D< Cartesian3D<double> \n";
#endif
   ret |= doTest <R, DisplacementVector3D< Cartesian3D<double> > >
   (testRotation,testVector,ticks);
#ifdef DEBUG
   std::cout << ">>>>> DisplacementVector3D< Polar3D<double> \n";
#endif
   ret |= doTest <R, DisplacementVector3D< Polar3D<double> > >
   (testRotation,testVector,ticks);
#ifdef DEBUG
   std::cout << ">>>>> DisplacementVector3D< CylindricalEta3D<double> \n";
#endif
   ret |= doTest <R, DisplacementVector3D< CylindricalEta3D<double> > >
   (testRotation,testVector,ticks);
#ifdef DEBUG
   std::cout << ">>>>> PositionVector3D< Cartesian3D<double> \n";
#endif
   ret |= doTest <R, PositionVector3D< Cartesian3D<double> > >
   (testRotation,testVector,ticks);
#ifdef DEBUG
   std::cout << ">>>>> PositionVector3D< Polar3D<double> \n";
#endif
   ret |= doTest <R, PositionVector3D< Polar3D<double> > >
   (testRotation,testVector,ticks);
#ifdef DEBUG
   std::cout << ">>>>> PositionVector3D< CylindricalEta3D<double> \n";
#endif
   ret |= doTest <R, PositionVector3D< CylindricalEta3D<double> > >
   (testRotation,testVector,ticks);
#ifdef DEBUG
   std::cout << ">>>>> ForeignVector\n";
#endif
   ret |= doTest <R, ForeignVector >
   (testRotation,testVector,ticks);
#ifdef DEBUG
   std::cout << ">>>>> LorentzVector<PxPyPzE4D<double> >\n";
#endif
   ret |= doTestL <R, PxPyPzE4D<double> >
   (testRotation,testVector,ticks);
#ifdef TRACE1
   std::cout << " ---- doTestofR ---- done\n";
#endif

   if (ret == 0) std::cout << ".";

   // TODO - other 4D coordinates

   return ret;
}


int exerciseTestCase (TestRotation const & testRotation,
                      XYZVector const & testVector)
{

   std::cout << ">>>>> Rotation Tests of " << testVector << "\t\t: " ;

#ifdef TRACE1
   std::cout << "---- exerciseTestCase ----";
#endif
   int ret = 0;
   ret |= doTestOfR <Rotation3D>  (testRotation,testVector);
   ret |= doTestOfR <AxisAngle>   (testRotation,testVector);
   ret |= doTestOfR <EulerAngles> (testRotation,testVector);
   ret |= doTestOfR <Quaternion>  (testRotation,testVector);
   ret |= doTestOfR <RotationZYX> (testRotation,testVector);
#ifdef TRACE1
   std::cout << " done\n";
#endif

   if (ret == 0)
   std::cout << "\t OK\n";
   else {
      std::cout << "\t Failed!\n ";
      std::cerr << "\n>>>>> Rotation Tests of " << testVector << "\t\t:\t FAILED \n";
   }

   return ret;
}

// ===== Axial test section =============

template <class R, class V>
int doTestA (XYZVector const & testVector, double ticks) {
#ifdef TRACE1
   std::cout << "---- doTestA ----";
#endif
   int ret = 0;
   V v(testVector);
   XYZVector rv;
   for (double angle = -4.0;  angle < 4.0; angle += .15) {
      RotationX rx (angle);
      rv = rxv (angle, testVector);
      ret |= testApplication (rx, v, rv, ticks);
      RotationY ry (angle);
      rv = ryv (angle, testVector);
      ret |= testApplication (ry, v, rv, ticks);
      RotationZ rz (angle);
      rv = rzv (angle, testVector);
      ret |= testApplication (rz, v, rv, ticks);
   }
#ifdef TRACE1
   std::cout << " done\n";
#endif
   if (ret == 0) std::cout << ".";
   return ret;
}

template <class R, class C>
int doTestLA (XYZVector const & testVector, double ticks) {
#ifdef TRACE1
   std::cout << "---- doTestLA ----";
#endif
   int ret = 0;
   LorentzVector<C> v;
   double x = testVector.X();
   double y = testVector.Y();
   double z = testVector.Z();
   double t = std::sqrt (x*x + y*y + z*z + 1);
   v.SetXYZT ( x, y, z, t );
   XYZVector rv;
   for (double angle = -4.0;  angle < 4.0; angle += .15) {
      //std::cout << "\n============ angle is " << angle << "\n";
      RotationX rx (angle);
      rv = rxv (angle, testVector);
      ret |= testApplication (rx, v, rv, ticks);
      RotationY ry (angle);
      rv = ryv (angle, testVector);
      ret |= testApplication (ry, v, rv, ticks);
      RotationZ rz (angle);
      rv = rzv (angle, testVector);
      ret |= testApplication (rz, v, rv, ticks);
   }
#ifdef TRACE1
   std::cout << " done\n";
#endif

   if (ret == 0) std::cout << ".";
   return ret;
}

template <class R>
int doTestOfAxial (XYZVector const & testVector){
#ifdef TRACE1
   std::cout << "---- doTestOfAxial ----\n";
#endif
   int ret = 0;
   const double ticks = 32;
#ifdef DEBUG
   std::cout << ">>>>> DisplacementVector3D< Cartesian3D<double> \n";
#endif
   ret |= doTestA <R, DisplacementVector3D< Cartesian3D<double> > >
   (testVector,ticks);
#ifdef DEBUG
   std::cout << ">>>>> DisplacementVector3D< Polar3D<double> \n";
#endif
   ret |= doTestA <R, DisplacementVector3D< Polar3D<double> > >
   (testVector,ticks);
#ifdef DEBUG
   std::cout << ">>>>> DisplacementVector3D< CylindricalEta3D<double> \n";
#endif
   ret |= doTestA <R, DisplacementVector3D< CylindricalEta3D<double> > >
   (testVector,ticks);
#ifdef DEBUG
   std::cout << ">>>>> PositionVector3D< Cartesian3D<double> \n";
#endif
   ret |= doTestA <R, PositionVector3D< Cartesian3D<double> > >
   (testVector,ticks);
#ifdef DEBUG
   std::cout << ">>>>> PositionVector3D< Polar3D<double> \n";
#endif
   ret |= doTestA <R, PositionVector3D< Polar3D<double> > > (testVector,ticks);
#ifdef DEBUG
   std::cout << ">>>>> PositionVector3D< CylindricalEta3D<double> \n";
#endif
   ret |= doTestA <R, PositionVector3D< CylindricalEta3D<double> > >
   (testVector,ticks);
#ifdef DEBUG
   std::cout << ">>>>> ForeignVector\n";
#endif
   ret |= doTestA <R, ForeignVector > (testVector,ticks);
#ifdef DEBUG
   std::cout << ">>>>> LorentzVector<PxPyPzE4D<double> >\n";
#endif
   ret |= doTestLA <R, PxPyPzE4D<double> > (testVector,ticks);
#ifdef TRACE1
   std::cout << " ---- doTestofR ---- done\n";
#endif
   // TODO - other 4D coordinates

   if (ret == 0) std::cout << ".";

   return ret;
}

int exerciseAxialTest (XYZVector const & testVector)
{

#ifdef TRACE1
   std::cout << "---- exerciseAxialTest ----";
#endif

   std::cout << ">>>>> Axial Rotation Tests of " << testVector << "\t\t: ";

   int ret = 0;
   ret |= doTestOfAxial <RotationX> (testVector);
   ret |= doTestOfAxial <RotationY> (testVector);
   ret |= doTestOfAxial <RotationZ> (testVector);
#ifdef TRACE1
   std::cout << " done\n";
#endif

   if (ret == 0)
   std::cout << "\t OK\n";
   else {
      std::cout << "\t Failed!\n ";
      std::cerr << "\n>>>>> Axial Rotation Tests of " << testVector << "\t\t:\t FAILED \n";
   }

   return ret;
}


#endif // endif on __CINT__

// ======================================


int rotationApplication (bool forceRun = false) {
   int ret = 0;

   bool skipTests = false;
#if defined(__i386__)
   // do not run by default tests on 32 bit architecture
   // since extended precision will make it difficult
   skipTests = true;
#endif

   if (skipTests && !forceRun) {
      std::cout << "Skip the tests - it is probably a 32 bit arch - return 0" << std::endl;
      return 0;
   }

   std::vector<TestRotation> testRotations = makeTestRotations();
   std::vector<XYZVector>    testVectors   = makeTestVectors();
   for ( std::vector<TestRotation>::const_iterator n =  testRotations.begin();
        n !=  testRotations.end(); ++n ) {
      for ( std::vector<XYZVector>::const_iterator m =  testVectors.begin();
           m !=  testVectors.end(); ++m ) {
         ret |= exerciseTestCase (*n, *m);
      }
   }
   for ( std::vector<XYZVector>::const_iterator vp =  testVectors.begin();
        vp !=  testVectors.end(); ++vp ) {
      ret |= exerciseAxialTest (*vp);
   }

   return ret;
}

int main() {
   int ret =  rotationApplication(true);
   if (ret)  std::cerr << "test FAILED !!! " << std::endl;
   else   std::cout << "test OK " << std::endl;
   return ret;
}