diff --git a/test/Makefile b/test/Makefile
index 0c4ab25694cf487f302430744aa787aa668e4664..f74342317d9372e9f0e331554ddc09daa07df198 100644
--- a/test/Makefile
+++ b/test/Makefile
@@ -73,11 +73,32 @@ STRESSVECO   = stressVector.$(ObjSuf)
 STRESSVECS   = stressVector.$(SrcSuf)
 STRESSVEC    = stressVector$(ExeSuf)
 
-STRESSMATHO   = stressMathCore.$(ObjSuf)
-STRESSMATHS   = stressMathCore.$(SrcSuf)
+ifneq ($(USE_REFLEX),)
+CXXFLAGS += -DUSE_REFLEX
+STRESSMATHO   = stressMathCore.$(ObjSuf) TrackMathCoreCint.$(ObjSuf) TrackMathCoreRflx.$(ObjSuf) 
+STRESSMATHS   = stressMathCore.$(SrcSuf) TrackMathCoreCint.$(SrcSuf) TrackMathCoreRflx.$(SrcSuf)
+ifeq ($(PLATFORM),win32)
+STRESSMATHLIBS = '$(ROOTSYS)/lib/libMathCore.lib' '$(ROOTSYS)/lib/libReflex.lib' '$(ROOTSYS)/lib/libCintex.lib' 
+else
+STRESSMATHLIBS = -lMathCore -lReflex -lCintex
+endif
+
+else 
+STRESSMATHO   = stressMathCore.$(ObjSuf) TrackMathCoreCint.$(ObjSuf) 
+STRESSMATHS   = stressMathCore.$(SrcSuf) TrackMathCoreCint.$(SrcSuf)
+
+ifeq ($(PLATFORM),win32)
+STRESSMATHLIBS = '$(ROOTSYS)/lib/libMathCore.lib'
+else
+STRESSMATHLIBS = -lMathCore
+endif
+ 
+endif
+
 STRESSMATH    = stressMathCore$(ExeSuf)
 
 
+
 endif
 
 VLAZYO        = vlazy.$(ObjSuf)
@@ -306,12 +327,84 @@ else
 endif
 		@echo "$@ done"
 
+
+#ifneq ($(useReflex),)
+
+
+
+libTrackMathCoreRflx.$(DllSuf): 	TrackMathCoreRflx.$(ObjSuf) 
+ifeq ($(ARCH),aix)
+		/usr/ibmcxx/bin/makeC++SharedLib $(OutPutOpt) $@ $(LIBS) -p 0 $^
+else
+ifeq ($(ARCH),aix5)
+		/usr/vacpp/bin/makeC++SharedLib $(OutPutOpt) $@ $(LIBS) -p 0 $^
+else
+ifeq ($(PLATFORM),macosx)
+# We need to make both the .dylib and the .so
+		$(LD) $(SOFLAGS) $^ $(OutPutOpt) $@
+ifeq ($(MACOSX_MINOR),4)
+		ln -sf $@ $(subst .$(DllSuf),.so,$@)
+else
+		$(LD) -bundle -undefined $(UNDEFOPT) $(LDFLAGS) $^ \
+		   $(OutPutOpt) $(subst .$(DllSuf),.so,$@)
+endif
+else
+ifeq ($(PLATFORM),win32)
+		bindexplib $* $^ > $*.def
+		lib -nologo -MACHINE:IX86 $^ -def:$*.def \
+		   $(OutPutOpt)$(EVENTLIB)
+		$(LD) $(SOFLAGS) $(LDFLAGS) $^ $*.exp $(LIBS) \
+		   $(OutPutOpt)$@
+		$(MT_DLL)
+else
+		$(LD) $(SOFLAGS) $(LDFLAGS) $^ $(OutPutOpt) $@ $(EXPLLINKLIBS)
+endif
+endif
+endif
+endif
+		@echo "$@ done"
+
+#endif
+
+libTrackMathCoreCint.$(DllSuf): 	TrackMathCoreCint.$(ObjSuf) 
+ifeq ($(ARCH),aix)
+		/usr/ibmcxx/bin/makeC++SharedLib $(OutPutOpt) $@ $(LIBS) -p 0 $^
+else
+ifeq ($(ARCH),aix5)
+		/usr/vacpp/bin/makeC++SharedLib $(OutPutOpt) $@ $(LIBS) -p 0 $^
+else
+ifeq ($(PLATFORM),macosx)
+# We need to make both the .dylib and the .so
+		$(LD) $(SOFLAGS) $^ $(OutPutOpt) $@
+ifeq ($(MACOSX_MINOR),4)
+		ln -sf $@ $(subst .$(DllSuf),.so,$@)
+else
+		$(LD) -bundle -undefined $(UNDEFOPT) $(LDFLAGS) $^ \
+		   $(OutPutOpt) $(subst .$(DllSuf),.so,$@)
+endif
+else
+ifeq ($(PLATFORM),win32)
+		bindexplib $* $^ > $*.def
+		lib -nologo -MACHINE:IX86 $^ -def:$*.def \
+		   $(OutPutOpt)$(EVENTLIB)
+		$(LD) $(SOFLAGS) $(LDFLAGS) $^ $*.exp $(LIBS) \
+		   $(OutPutOpt)$@
+		$(MT_DLL)
+else
+		$(LD) $(SOFLAGS) $(LDFLAGS) $^ $(OutPutOpt) $@ $(EXPLLINKLIBS)
+endif
+endif
+endif
+endif
+		@echo "$@ done"
+
+
 $(STRESSMATH):   $(STRESSMATHO)
 ifeq ($(PLATFORM),win32)
-		$(LD) $(LDFLAGS) $^ $(LIBS) ../lib/libMathCore.lib  $(OutPutOpt)$@
+		$(LD) $(LDFLAGS) $^ $(LIBS) $(STRESSMATHLIBS)  $(OutPutOpt)$@
 		$(MT_EXE)
 else
-		$(LD) $(LDFLAGS) $^ $(LIBS) -lMathCore  $(OutPutOpt)$@
+		$(LD) $(LDFLAGS) $^ $(LIBS) $(STRESSMATHLIBS)  $(OutPutOpt)$@
 endif
 		@echo "$@ done"
 
@@ -558,5 +651,16 @@ guiviewerDict.$(SrcSuf): guiviewer.h guiviewerLinkDef.h
 	@echo "Generating dictionary $@..."
 	@rootcint -f $@ -c $^
 
+stressMathCore.$(ObjSuf): 	TrackMathCore.h
+TrackMathCoreCint.$(SrcSuf): 	TrackMathCore.h TrackMathCoreLinkDef.h
+	@echo "Generating dictionary $@ using rootcint ..."
+	@rootcint -f $@ -c $^	
+
+TrackMathCoreRflx.$(SrcSuf): 	TrackMathCore.h TrackMathCoreRflx.xml
+			@echo "Generating dictionary $@ using gccxml ..."
+			genreflex TrackMathCore.h --selection_file=TrackMathCoreRflx.xml -o TrackMathCoreRflx.cxx  -I$(ROOTSYS)/include
+
+
+
 .$(SrcSuf).$(ObjSuf):
 	$(CXX)  $(CXXFLAGS) -c $<
diff --git a/test/TrackMathCore.h b/test/TrackMathCore.h
new file mode 100644
index 0000000000000000000000000000000000000000..ff82857680b896d451c64682e2afc7bb5e73306f
--- /dev/null
+++ b/test/TrackMathCore.h
@@ -0,0 +1,233 @@
+// dummy track class for testing I/o of matric
+
+#include "Math/Point3D.h"
+#include "Math/Vector4D.h"
+#include "Math/SMatrix.h"
+#include "Rtypes.h" // for Double32_t
+
+#include <vector>
+#include <string>
+#include <iostream>
+#include <cassert>
+
+
+
+typedef ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> >        Vector4D; 
+typedef ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<Double32_t> >    Vector4D32;
+
+typedef ROOT::Math::DisplacementVector3D<ROOT::Math::Cartesian3D<double> >     Vector3D; 
+
+typedef ROOT::Math::PositionVector3D<ROOT::Math::Cartesian3D<double> >        Point3D; 
+typedef ROOT::Math::PositionVector3D<ROOT::Math::Cartesian3D<Double32_t> >    Point3D32; 
+
+typedef ROOT::Math::SMatrix<double,4,4,ROOT::Math::MatRepStd<double,4,4> >      Matrix4D;
+typedef ROOT::Math::SMatrix<Double32_t,4,4,ROOT::Math::MatRepStd<Double32_t,4,4> >      Matrix4D32;
+
+typedef ROOT::Math::SMatrix<double,6,6,ROOT::Math::MatRepSym<double,6> >      SymMatrix6D;
+typedef ROOT::Math::SMatrix<Double32_t,6,6,ROOT::Math::MatRepSym<Double32_t,6> >      SymMatrix6D32;;
+
+
+// track class 
+class  TrackD { 
+
+public:
+
+   typedef Matrix4D::const_iterator const_iterator; 
+
+   TrackD() {}
+   
+   TrackD(double * begin, double * end)  {
+      double * itr = begin; 
+      fPos.SetCoordinates(itr, itr+3); itr +=3;
+      fVec.SetCoordinates(itr,itr+4); itr+=4;
+      fMat = Matrix4D(itr,itr+16); itr += 16;
+      fSymMat = SymMatrix6D(itr,itr+21); 
+      assert(itr+21 == end);
+   }
+
+   enum {  kSize =  
+      3 + 4 + Matrix4D32::kSize + SymMatrix6D32::rep_type::kSize
+   };
+
+   static std::string Type() { return "TrackD"; } 
+
+   static bool IsD32() { return false; }
+
+   TrackD & operator += (const TrackD & t) { 
+      fPos += Vector3D(t.fPos);
+      fVec += t.fVec;
+      fMat += t.fMat;
+      fSymMat += t.fSymMat;
+      return *this;
+   }
+
+   double Sum() const { 
+      double s = 0; 
+      double d[4]; 
+      fPos.GetCoordinates(d,d+3);
+      for (int i=0; i<3;++i)
+         s+= d[i];
+      fVec.GetCoordinates(d,d+4);
+      for (int i=0; i<4;++i)
+         s+= d[i];
+      for (const_iterator itr = fMat.begin(); itr != fMat.end(); ++itr)
+         s+= *itr;
+      for (const_iterator itr = fSymMat.begin(); itr != fSymMat.end(); ++itr)
+         s+= *itr;
+      return s;
+  }
+
+   void Print() const { 
+     std::cout << "Point  " << fPos << std::endl;
+     std::cout << "Vec    " << fVec << std::endl;
+     std::cout << "Mat    " << fMat << std::endl;
+     std::cout << "SymMat " << fSymMat << std::endl;
+   }
+
+private:
+
+   Point3D    fPos;
+   Vector4D   fVec; 
+   Matrix4D fMat; 
+   SymMatrix6D fSymMat; 
+      
+}; 
+
+// track class based on  of Double32
+
+
+class  TrackD32 { 
+
+public:
+
+   typedef Matrix4D::const_iterator const_iterator; 
+
+   TrackD32() {}
+
+   enum {  kSize =  
+      3 + 4 + Matrix4D32::kSize + SymMatrix6D32::rep_type::kSize
+   };
+
+   static std::string Type() { return "TrackD32"; } 
+
+   static bool IsD32() { return true; }
+   
+   TrackD32(double * begin, double * end)  {
+      double * itr = begin; 
+      fPos.SetCoordinates(itr, itr+3); itr +=3;
+      fVec.SetCoordinates(itr,itr+4); itr+=4;
+      fMat = Matrix4D32(itr,itr+16); itr += 16;
+      fSymMat = SymMatrix6D32(itr,itr+21);
+      assert(itr+21 == end);
+      
+   }
+
+
+   TrackD32 & operator += (const TrackD32 & t) { 
+      fPos += Vector3D(t.fPos);
+      fVec += t.fVec;
+      fMat += t.fMat;
+      fSymMat += t.fSymMat;
+      return *this;
+   }
+
+   double Sum() const { 
+      double s = 0; 
+      double d[4]; 
+      fPos.GetCoordinates(d,d+3);
+      for (int i=0; i<3;++i)
+         s+= d[i];
+      fVec.GetCoordinates(d,d+4);
+      for (int i=0; i<4;++i)
+         s+= d[i];
+      for (const_iterator itr = fMat.begin(); itr != fMat.end(); ++itr)
+         s+= *itr;
+      for (const_iterator itr = fSymMat.begin(); itr != fSymMat.end(); ++itr)
+         s+= *itr;
+      return s;
+  }
+
+   void Print() const { 
+     std::cout << "Point  " << fPos << std::endl;
+     std::cout << "Vec    " << fVec << std::endl;
+     std::cout << "Mat    " << fMat << std::endl;
+     std::cout << "SymMat " << fSymMat << std::endl;
+   }
+
+private:
+
+   Point3D32    fPos;
+   Vector4D32   fVec; 
+   Matrix4D32 fMat; 
+   SymMatrix6D32 fSymMat; 
+      
+}; 
+// class containning a vector of tracks
+
+class VecTrackD {
+
+public: 
+
+   typedef std::vector<TrackD>::const_iterator It;
+   
+ 
+   VecTrackD() {
+     // create klen empty trackD
+       fTrks.reserve(kLen);
+       for (int i = 0; i < kLen; ++i) { 
+         fTrks.push_back(TrackD() ); 
+       }
+   }
+
+   VecTrackD(double * begin, double * end) 
+     {
+       fTrks.reserve(kLen);
+       double * itr = begin;
+       for (int i = 0; i < kLen; ++i) { 
+         fTrks.push_back(TrackD(itr, itr + TrackD::kSize) ); 
+         itr += TrackD::kSize;
+       }
+       assert( itr == end); 
+     }
+
+   enum {  kLen = 10, kSize =  kLen*TrackD::kSize };
+
+   static std::string Type() { return "VecTrackD"; } 
+
+   static bool IsD32() { return false; }
+
+
+   VecTrackD & operator += (const VecTrackD & v) { 
+    for (unsigned int i = 0; i < fTrks.size() ; ++i) 
+       fTrks[i] += v.fTrks[i];
+
+
+    return *this;
+   }
+
+
+  It begin() const { return fTrks.begin(); } 
+  It end() const  { return fTrks.end(); } 
+
+
+  double Sum() const { 
+    double s = 0; 
+    for (unsigned int i = 0; i < fTrks.size() ; ++i) 
+      s += fTrks[i].Sum(); 
+    
+    return s; 
+  }
+
+  
+   void Print() const { 
+     for (unsigned int i = 0; i < fTrks.size() ; ++i) {
+       std::cout << "\n======> Track ========<" << i << std::endl;
+       fTrks[i].Print();
+     }
+   }
+
+private:
+
+  std::vector<TrackD>  fTrks;
+
+};
diff --git a/test/TrackMathCoreLinkDef.h b/test/TrackMathCoreLinkDef.h
new file mode 100644
index 0000000000000000000000000000000000000000..8bf15d2a5855c8819c228c737c0679ee2945fcf1
--- /dev/null
+++ b/test/TrackMathCoreLinkDef.h
@@ -0,0 +1,21 @@
+#ifdef __CINT__
+
+
+#pragma link off all globals;
+#pragma link off all classes;
+#pragma link off all functions;
+
+
+
+#pragma link C++ class TrackD+;
+#pragma link C++ class TrackD32+;
+
+#pragma extra_include "vector";
+#include <vector>
+
+#pragma link C++ class vector<TrackD >+;
+
+#pragma link C++ class VecTrackD+;
+
+
+#endif
diff --git a/test/TrackMathCoreRflx.xml b/test/TrackMathCoreRflx.xml
new file mode 100644
index 0000000000000000000000000000000000000000..448fcb1215aefd72d619155b40c1c18e761b91dd
--- /dev/null
+++ b/test/TrackMathCoreRflx.xml
@@ -0,0 +1,16 @@
+<lcgdict>
+
+
+
+
+<class name ="TrackD">
+</class>
+ <class name ="TrackD32">
+  <field name="fVec" iotype="ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<Double32_t> >" />
+  <field name="fPos" iotype="ROOT::Math::PositionVector3D<ROOT::Math::Cartesian3D<Double32_t>,ROOT::Math::DefaultCoordinateSystemTag>" />
+</class>
+<class pattern="std::vector<TrackD*>" />
+<class name="VecTrackD" />
+
+ 
+</lcgdict>
diff --git a/test/stressMathCore.cxx b/test/stressMathCore.cxx
index 38cec5cb502a4efd82a18a04a3d0de7355f4e6e2..11f705b7937be723aee45bcaa7f831f78bc79131 100644
--- a/test/stressMathCore.cxx
+++ b/test/stressMathCore.cxx
@@ -1,4 +1,4 @@
-// @(#)root/test:$Id: stressVector.cxx 19826 2007-09-19 19:56:11Z rdm $
+// @(#)root/test:$Id$
 // Author: Lorenzo Moneta   06/2005 
 ///////////////////////////////////////////////////////////////////////////////////
 //
@@ -36,13 +36,24 @@
 #include "Math/SVector.h"
 #include "Math/SMatrix.h"
 
+#include "TrackMathCore.h"
+
+
+//#define USE_REFLEX
+#ifdef USE_REFLEX
+#include "Cintex/Cintex.h"
+#include "Reflex/Reflex.h"
+#endif
+
 using namespace ROOT::Math; 
 
 #endif
 
 
+
 bool debug = true;  // print out reason of test failures
 bool debugTime = false; // print out separate timings for vectors 
+bool removeFiles = true; // remove Output root files 
 
 void PrintTest(std::string name) { 
    std::cout << std::left << std::setw(40) << name; 
@@ -428,11 +439,12 @@ int testStatFunctions(int /* nfunc */) {
 
 template<class V> 
 struct VecType { 
-   static std::string name() { return "Vector";}
+   static std::string name() { return "MathCoreVector";}
 }; 
 template<>
 struct VecType<XYVector> {
    static std::string name() { return "XYVector";}
+   static std::string name32() { return "ROOT::Math::DisplacementVector2D<ROOT::Math::Cartesian2D<Double32_t> >";}
 }; 
 template<>
 struct VecType<Polar2DVector> {
@@ -441,6 +453,7 @@ struct VecType<Polar2DVector> {
 template<>
 struct VecType<XYZVector> {
    static std::string name() { return "XYZVector";}
+   static std::string name32() { return "ROOT::Math::DisplacementVector3D<ROOT::Math::Cartesian3D<Double32_t> >";}
 }; 
 template<>
 struct VecType<Polar3DVector> {
@@ -457,6 +470,7 @@ struct VecType<RhoZPhiVector> {
 template<>
 struct VecType<PxPyPzEVector> {
    static std::string name() { return "PxPyPzEVector";}
+   static std::string name32() { return "ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<Double32_t> >";}
 }; 
 template<>
 struct VecType<PtEtaPhiEVector> {
@@ -470,6 +484,29 @@ template<>
 struct VecType<PxPyPzMVector> {
    static std::string name() { return "PxPyPzMVector";}
 }; 
+template<>
+struct VecType<SVector<double,3> > {
+   static std::string name() { return "SVector3";}
+   static std::string name32() { return "ROOT::Math::SVector<Double32_t,3>";}
+}; 
+template<>
+struct VecType<SVector<double,4> > {
+   static std::string name() { return "SVector4";}
+   static std::string name32() { return "ROOT::Math::SVector<Double32_t,4>";}
+}; 
+template<>
+struct VecType<TrackD> {
+   static std::string name() { return "TrackD";}
+}; 
+template<>
+struct VecType<TrackD32> {
+   static std::string name() { return "TrackD32";}
+}; 
+template<>
+struct VecType<VecTrackD> {
+   static std::string name() { return "VecTrackD";}
+}; 
+
 
 // generic (2 dim)
 template<class V, int Dim>  
@@ -810,6 +847,7 @@ public:
       if (typeName == "") {
          typeName = "ROOT::Math::" + VecType<V>::name();
       }
+      //std::cout << typeName << std::endl;
 
       TBranch * br = tree.Branch("Vector branch",typeName.c_str(),&v1);
       if (br == 0) { 
@@ -825,7 +863,11 @@ public:
          tree.Fill();
       }
 
-      return file.Write();
+      //tree.Print(); // debug
+
+      file.Write();
+      file.Close();
+      return file.GetSize();
    }
 
    template<class V> 
@@ -868,7 +910,9 @@ public:
          tree->GetEntry(i);
          dataV.push_back(*v1); 
       }
-      gSystem->Unlink(fname.c_str());
+
+      if (removeFiles) gSystem->Unlink(fname.c_str());
+
       
       return 0; 
    }
@@ -906,9 +950,23 @@ public:
       return tot;
    }  
 
-};
+   template <class V>
+   double testAdditionTR( const std::vector<V > & dataV) { 
+      V v0;
+      Timer t;
+      for (int i = 0; i < nGen; ++i) { 
+         v0 += dataV[i]; 
+      }
+      //v0.Print();
+      return v0.Sum();
+   }  
 
 
+};
+
+//--------------------------------------------------------------------------------------
+// test of all physics vector (GenVector's)
+//--------------------------------------------------------------------------------------
 
 template<class V1, class V2, int Dim> 
 int testVector(int ngen, bool testio=false) { 
@@ -930,7 +988,7 @@ int testVector(int ngen, bool testio=false) {
    double sref1, sref2 = 0; 
 
    a.testCreate(v1);             iret |= a.check(VecType<V1>::name()+" creation",v1.size(),ngen);
-   s1 = a.testAddition(v1);    iret |= a.check(VecType<V1>::name()+" addition",s1,a.Sum(),Dim*4);
+   s1 = a.testAddition(v1);    iret |= a.check(VecType<V1>::name()+" addition",s1,a.Sum(),Dim*8);
    sref1 = s1; 
    v1.clear();
    assert(v1.size() == 0);
@@ -959,21 +1017,38 @@ int testVector(int ngen, bool testio=false) {
    int ir = 0;
    if (!testio) return iret; 
 
-   fsize = a.testWrite(v1);  iret |= a.check(VecType<V1>::name()+" write",fsize>100,1);
+   double estSize = ngen*8 * Dim + 10 * ngen;  //add extra bytes
+   scale = 0.1 / std::numeric_limits<double>::epsilon();
+   fsize = a.testWrite(v1);  iret |= a.check(VecType<V1>::name()+" write",fsize,estSize,scale);
    ir = a.testRead(v1);   iret |= a.check(VecType<V1>::name()+" read",ir,0);
    s1 = a.testAddition(v1);       iret |= a.check(VecType<V1>::name()+" after read",s1,sref1);
 
    // test io vector 2
-   fsize = a.testWrite(v2);  iret |= a.check(VecType<V2>::name()+" write",fsize>100,1);
+   fsize = a.testWrite(v2);  iret |= a.check(VecType<V2>::name()+" write",fsize,estSize,scale);
    ir = a.testRead(v2);      iret |= a.check(VecType<V2>::name()+" read",ir,0);
    s2 = a.testAddition(v2);       iret |= a.check(VecType<V2>::name()+" after read",s2,sref2);
 
+
+   // test io of double 32 for vector 1
+   if (Dim==2) return iret; 
+
+   estSize = ngen*4 * Dim + 10 * ngen;  //add extra bytes
+   scale = 0.1 / std::numeric_limits<double>::epsilon();
+   
+   std::string typeName = VecType<V1>::name32();
+   fsize = a.testWrite(v1,typeName);  iret |= a.check(VecType<V1>::name() +"_D32 write",fsize,estSize,scale);
+   ir = a.testRead(v1);   iret |= a.check(VecType<V1>::name()+"_D32 read",ir,0);
+   s1 = a.testAddition(v1);       iret |= a.check(VecType<V1>::name()+"_D32 after read",s1,sref1,1.E9);
+
+
    return iret; 
 
 }
 
-
+//--------------------------------------------------------------------------------------
 // test of Svector of dim 3 or 4
+//--------------------------------------------------------------------------------------
+
 template<int Dim> 
 int testVector34(int ngen, bool testio=false) { 
    
@@ -994,7 +1069,7 @@ int testVector34(int ngen, bool testio=false) {
 
    std::string name = "SVector<double," + Util::ToString(Dim) + ">"; 
    a.testCreate(v1);             iret |= a.check(name+" creation",v1.size(),ngen);
-   s1 = a.testAddition(v1);    iret |= a.check(name+" addition",s1,a.Sum(),Dim*4);
+   s1 = a.testAddition(v1);    iret |= a.check(name+" addition",s1,a.Sum(),Dim*8);
    sref1 = s1; 
 
    // test the io
@@ -1003,14 +1078,32 @@ int testVector34(int ngen, bool testio=false) {
    if (!testio) return iret; 
 
    std::string typeName = "ROOT::Math::"+name;
-   fsize = a.testWrite(v1,typeName);  iret |= a.check(name+" write",fsize>100,1);
+
+   double estSize = ngen*8 * Dim + 10 * ngen;
+   double scale = 0.1 / std::numeric_limits<double>::epsilon();
+   fsize = a.testWrite(v1,typeName);  iret |= a.check(name+" write",fsize,estSize,scale);
    ir = a.testRead(v1);   iret |= a.check(name+" read",ir,0);
    s1 = a.testAddition(v1);       iret |= a.check(name+" after read",s1,sref1);
 
+   //std::cout << "File size = " << fsize << " estimated " << 8 * Dim * ngen << std::endl;
+
+   // test Double32
+   estSize = ngen*4 * Dim + 10 * ngen;  //add extra bytes
+   scale = 0.1 / std::numeric_limits<double>::epsilon();
+   
+   typeName = VecType<SV>::name32();
+   fsize = a.testWrite(v1,typeName);  iret |= a.check(VecType<SV>::name() +"_D32 write",fsize,estSize,scale);
+   ir = a.testRead(v1);   iret |= a.check(VecType<SV>::name()+"_D32 read",ir,0);
+   s1 = a.testAddition(v1);       iret |= a.check(VecType<SV>::name()+"_D32 after read",s1,sref1,1.E9);
+
+
    return iret; 
 }
 
+//--------------------------------------------------------------------------------------
 // test of generic Svector 
+//--------------------------------------------------------------------------------------
+
 template<int Dim> 
 int testSVector(int ngen, bool testio=false) { 
    
@@ -1040,11 +1133,15 @@ int testSVector(int ngen, bool testio=false) {
    int ir = 0;
    if (!testio) return iret; 
 
+
    std::string typeName = "ROOT::Math::"+name;
-   fsize = a.testWrite(v1,typeName);  iret |= a.check(name+" write",fsize>100,1);
+   double estSize = ngen*8 * Dim + 10. * ngen;
+   double scale = 0.1 / std::numeric_limits<double>::epsilon();
+   fsize = a.testWrite(v1,typeName);  iret |= a.check(name+" write",fsize,estSize,scale);
    ir = a.testRead(v1);   iret |= a.check(name+" read",ir,0);
    s1 = a.testAdditionSV(v1);       iret |= a.check(name+" after read",s1,sref1);
 
+
    return iret; 
 }
 
@@ -1054,6 +1151,9 @@ struct RepStd {
    static std::string name() {
       return "ROOT::Math::MatRepStd<double," +  Util::ToString(D1) + "," + Util::ToString(D2) + "> "; 
    }
+   static std::string name32() {
+      return "ROOT::Math::MatRepStd<Double32_t," +  Util::ToString(D1) + "," + Util::ToString(D2) + "> "; 
+   }
    static std::string sname() { return ""; } 
 };
 template<int D1>
@@ -1062,10 +1162,16 @@ struct RepSym {
    static std::string name() {
       return "ROOT::Math::MatRepSym<double," +  Util::ToString(D1) + "> "; 
    }
+   static std::string name32() {
+      return "ROOT::Math::MatRepSym<Double32_t," +  Util::ToString(D1) + "> "; 
+   }
    static std::string sname() { return "MatRepSym"; } 
 };
 
+//--------------------------------------------------------------------------------------
 // test of generic SMatrix
+//--------------------------------------------------------------------------------------
+
 template<int D1, int D2,class Rep > 
 int testSMatrix(int ngen, bool testio=false) { 
    
@@ -1097,17 +1203,99 @@ int testSMatrix(int ngen, bool testio=false) {
    sref1 = s1; 
 
    // test the io
+   if (!testio) return iret; 
    double fsize = 0;
    int ir = 0;
-   if (!testio) return iret; 
 
    // the full name is needed for sym matrices
    std::string typeName = "ROOT::Math::"+name0 + "," + Rep::name()  + ">";
-   //std::string typeName = "ROOT::Math::"+name;
-   fsize = a.testWrite(v1,typeName);  iret |= a.check(name+" write",fsize>100,1);
+
+   double estSize = ngen*8 * Dim + 10 * ngen;
+   double scale = 0.1 / std::numeric_limits<double>::epsilon();
+   fsize = a.testWrite(v1,typeName);  iret |= a.check(name+" write",fsize,estSize,scale);
    ir = a.testRead(v1);   iret |= a.check(name+" read",ir,0);
    s1 = a.testAdditionSV(v1);       iret |= a.check(name+" after read",s1,sref1);
 
+
+   //std::cout << "File size = " << fsize << " estimated " << 8 * Dim * ngen << std::endl;
+
+   // test storing as Double32_t
+   // dictionay exist only for square matrices between 3 and 6 
+   if (D1 != D2) return iret; 
+   if (D1 < 3 || D1 > 6) return iret; 
+
+   double fsize32 = 0;
+
+   name0 = "SMatrix<Double32_t," + Util::ToString(D1) + "," + Util::ToString(D2);
+   name = name0 + "," + Rep::sname() +  ">";
+   typeName = "ROOT::Math::"+name0+ "," + Rep::name32()  + ">";
+
+
+   estSize = ngen* 4 * Dim + 10 * ngen;
+   scale = 0.1 / std::numeric_limits<double>::epsilon();
+   fsize32 = a.testWrite(v1,typeName);     iret |= a.check(name+" write",fsize32,estSize,scale);
+   ir = a.testRead(v1);   iret |= a.check(name+" read",ir,0);
+   //we read back float (scale errors then)
+   s1 = a.testAdditionSV(v1);       iret |= a.check(name+" after read",s1,sref1,1.E9);
+
+
+   return iret; 
+}
+
+
+//--------------------------------------------------------------------------------------
+// test of a track an object containing vector and matrices
+//--------------------------------------------------------------------------------------
+template<class T>
+int testTrack(int ngen) { 
+   
+   int iret = 0;
+
+   
+   const int Dim = T::kSize;
+   std::string name = T::Type();
+   
+   std::cout << "Dim is " << Dim << std::endl;
+   
+   VectorTest<Dim> a(ngen); 
+   a.genDataN();
+
+   std::vector<T> v1; 
+   v1.reserve(ngen); 
+
+   double s1 = 0; 
+   //double scale = 1; 
+   double sref1  = 0; 
+
+ 
+   a.testCreateSV(v1);             iret |= a.check(name+" creation",v1.size(),ngen);
+   s1 = a.testAdditionTR(v1);    iret |= a.check(name+" addition",s1,a.Sum(),Dim*4);
+   sref1 = s1; 
+   std::cout << " reference sum is " << sref1 << std::endl;
+
+   double fsize = 0;
+   int ir = 0;
+
+   // the full name is needed for sym matrices
+   std::string typeName = name;
+
+   int wsize = 8; 
+   if (T::IsD32() ) wsize = 4;
+   double estSize = ngen*wsize * Dim + 10 * ngen;
+   double scale = 0.1 / std::numeric_limits<double>::epsilon();
+   fsize = a.testWrite(v1,typeName);  iret |= a.check(name+" write",fsize,estSize,scale);
+   ir = a.testRead(v1);   iret |= a.check(name+" read",ir,0);
+   scale = 1; 
+   if (T::IsD32() ) scale = 1.E9; // use float numbers
+   s1 = a.testAdditionTR(v1);       iret |= a.check(name+" after read",s1,sref1,scale);
+   if (T::IsD32() ) { 
+     // check at double precision type must fail otherwise Double's are stored
+     bool pdebug = debug; debug = false; 
+     int ret = compare(" ",s1,sref1);
+     iret |= a.check(name+" Double32 test",ret,1);
+     debug = pdebug; 
+   } 
+
    return iret; 
 }
 
@@ -1155,6 +1343,7 @@ int testSMatrix(int ngen,bool io) {
 
    // asymetric matrices we have only 4x3 and 3x4
    iret |= testSMatrix<3,4,RepStd<3,4> >(ngen,io); 
+
    iret |= testSMatrix<4,3,RepStd<4,3> >(ngen,io); 
    iret |= testSMatrix<3,3,RepStd<3,3> >(ngen,io); 
    // sym matrix
@@ -1163,6 +1352,39 @@ int testSMatrix(int ngen,bool io) {
    return iret; 
 }
 
+int testCompositeObj(int ngen) { 
+
+   int iret = 0; 
+
+   std::cout <<"******************************************************************************\n";
+   std::cout << "\tTest of a Composite Object (containing Vector's and Matrices)\n";
+   std::cout <<"******************************************************************************\n";
+
+#ifndef USE_REFLEX
+   std::cout << "Test Using CINT library\n\n"; 
+
+//    iret = gSystem->Load("libTrackMathCoreCint");
+//    if (iret !=0) { 
+//       std::cerr <<"Error Loading libTrackMathCoreCint" << std::endl;
+//    }
+
+   iret |= testTrack<TrackD>(ngen); 
+   iret |= testTrack<TrackD32>(ngen); 
+   iret |= testTrack<VecTrackD>(ngen); 
+
+#else 
+   std::cout << "Test Using Reflex library\n\n"; 
+
+   ROOT::Cintex::Cintex::SetDebug(1);
+   ROOT::Cintex::Cintex::Enable();
+
+   iret |= testTrack<TrackD>(ngen); 
+   iret |= testTrack<TrackD32>(ngen); 
+   iret |= testTrack<VecTrackD>(ngen); 
+#endif
+
+   return iret;
+}
 
 int stressMathCore(double nscale = 1) { 
 
@@ -1177,7 +1399,7 @@ int stressMathCore(double nscale = 1) {
 //    if (iret !=0) return iret; 
 
 
-   TBenchmark bm; 
+   TBenchmark bm;
    bm.Start("stressMathCore");
 
    std::cout << nscale*100 << std::endl;
@@ -1189,6 +1411,8 @@ int stressMathCore(double nscale = 1) {
 
    iret |= testSMatrix(int(nscale*1000),io); 
 
+   iret |= testCompositeObj(int(nscale*1000)); 
+
    bm.Stop("stressMathCore");
    std::cout <<"******************************************************************************\n";
    bm.Print("stressMathCore");