// dummy track class for testing I/o of matric

//#ifdef USE_REFLEX
//#define G__DICTIONARY
//#endif

#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 containing a vector and a point
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;
      assert(itr == end);
   }

   enum {  kSize =        3 + 4    };

   static std::string Type() { return "TrackD"; } 

   static bool IsD32() { return false; }

   TrackD & operator += (const TrackD & t) { 
      fPos += Vector3D(t.fPos);
      fVec += t.fVec;
      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];
      return s;
  }

   void Print() const { 
     std::cout << "Point  " << fPos << std::endl;
     std::cout << "Vec    " << fVec << std::endl;
   }

private:

   Point3D    fPos;
   Vector4D   fVec; 
      
}; 

// track class based on  of Double32


class  TrackD32 { 

public:

   typedef Matrix4D::const_iterator const_iterator; 

   TrackD32() {}

   enum {  kSize =        3 + 4    };

   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;
      assert(itr == end);
      
   }


   TrackD32 & operator += (const TrackD32 & t) { 
      fPos += Vector3D(t.fPos);
      fVec += t.fVec;
      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];
      return s;
  }

   void Print() const { 
     std::cout << "Point  " << fPos << std::endl;
     std::cout << "Vec    " << fVec << std::endl;
   }

private:

   Point3D32     fPos;
   Vector4D32    fVec; 

      
}; 


// track class  (containing a vector and a point) and matrices

class  TrackErrD { 

public:

   typedef Matrix4D::const_iterator const_iterator; 

   TrackErrD() {}
   
   TrackErrD(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 "TrackErrD"; } 

   static bool IsD32() { return false; }

   TrackErrD & operator += (const TrackErrD & 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  TrackErrD32 { 

public:

   typedef Matrix4D::const_iterator const_iterator; 

   TrackErrD32() {}

   enum {  kSize =  
      3 + 4 + Matrix4D32::kSize + SymMatrix6D32::rep_type::kSize
   };

   static std::string Type() { return "TrackErrD32"; } 

   static bool IsD32() { return true; }
   
   TrackErrD32(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);
      
   }


   TrackErrD32 & operator += (const TrackErrD32 & 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

template<class T> 
class VecTrack {

public: 

   typedef typename std::vector<T>::const_iterator It;
   
 
   VecTrack() {
     // create klen empty trackD
       fTrks.reserve(kLen);
       for (int i = 0; i < kLen; ++i) { 
         fTrks.push_back(T() ); 
       }
   }

   VecTrack(double * ibegin, double * iend) 
     {
       fTrks.reserve(kLen);
       double * itr = ibegin;
       for (int i = 0; i < kLen; ++i) { 
         fTrks.push_back(T(itr, itr + T::kSize) ); 
         itr += T::kSize;
       }
       assert( itr == iend); 
     }

   enum {  kLen = 3, kSize =  kLen*T::kSize };

   static std::string Type() { return "VecTrack<"+T::Type()+">"; } 

   static bool IsD32() { return false; }


   VecTrack & operator += (const VecTrack & 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<T>  fTrks;

};

// for instantiating the template VecTrackD class for reflex
struct Dummy { 

  VecTrack<TrackD>     v1; 
  VecTrack<TrackErrD> v2; 
};