//
// Cint macro to test I/O of mathcore Lorentz Vectors in a Tree and compare with a 
// TLorentzVector. A ROOT tree is written and read in both using either a XYZTVector or /// a TLorentzVector. 
// 
//  To execute the macro type in: 
//
// root[0]: .x  mathcoreVectorIO.C


#ifdef USE_REFLEX
#include "Cintex/Cintex.h"
#include "Reflex/Reflex.h"
#endif

#include "TRandom3.h"
#include "TStopwatch.h"
#include "TSystem.h"
#include "TFile.h"
#include "TTree.h"
#include "TH1D.h"
#include "TCanvas.h"

#include <iostream>

#include "TLorentzVector.h"

#include "Math/Vector4D.h"
#include "Math/Vector3D.h"
#include "Math/Point3D.h"

#include "Track.h"


#define DEBUG


//#define READONLY

//#define USE_REFLEX

#define DEFVECTOR4D(TYPE) \
typedef TYPE AVector4D; \
const std::string vector4d_type = #TYPE ;

#define DEFVECTOR3D(TYPE) \
typedef TYPE AVector3D; \
const std::string vector3d_type = #TYPE ;

#define DEFPOINT3D(TYPE) \
typedef TYPE APoint3D; \
const std::string point3d_type = #TYPE ;



//const double tol = 1.0E-16;
const double tol = 1.0E-6; // or doublr 32 or float

DEFVECTOR4D(ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<Double32_t> >);
//DEFVECTOR4D(ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<Double32_t> >);

DEFVECTOR3D(ROOT::Math::DisplacementVector3D<ROOT::Math::Cartesian3D<Double32_t> >); 

DEFPOINT3D(ROOT::Math::PositionVector3D<ROOT::Math::Cartesian3D<Double32_t> >); 


// DEFVECTOR4D(ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<Double32_t> >) 

//using namespace ROOT::Math;

template<class Vector> 
inline double getMag2(const Vector & v) { 
  return v.mag2();
}

// inline double getMag2(const VecTrackD & v) { 
//    // print the read points
//    std::cout << "VecTRackD " << std::endl;
//    for (VecTrackD::It itr = v.begin() ; itr != v.end(); ++itr) 
//       std::cout << (*itr).Pos() << std::endl; 

//   return v.mag2();
// }


inline double getMag2(const TVector3 & v) { 
  return v.Mag2();
}

inline double getMag2(const TLorentzVector & v) { 
  return v.Mag2();
}

template<class U> 
inline void setValues(ROOT::Math::DisplacementVector3D<U> & v, const double x[] ) { 
    v.SetXYZ(x[0],x[1],x[2]); 
}

template<class U> 
inline void setValues(ROOT::Math::PositionVector3D<U> & v, const double x[]) { 
    v.SetXYZ(x[0],x[1],x[2]); 
}

template<class U> 
inline void setValues(ROOT::Math::LorentzVector<U> & v, const double x[]) { 
    v.SetXYZT(x[0],x[1],x[2],x[3]); 
}
// specialization for T -classes 
inline void setValues(TVector3 & v, const double x[]) { 
    v.SetXYZ(x[0],x[1],x[2]); 
}
inline void setValues(TLorentzVector & v, const double x[]) { 
    v.SetXYZT(x[0],x[1],x[2],x[3]); 
}


template <class Vector> 
double testDummy(int n) { 

  TRandom3 R(111);

  TStopwatch timer;

  Vector v1;
  double s = 0; 
  double p[4];

  timer.Start();
  for (int i = 0; i < n; ++i) { 
        p[0] = R.Gaus(0,10);
	p[1] = R.Gaus(0,10);
	p[2] = R.Gaus(0,10);
	p[3] = R.Gaus(100,10);
	setValues(v1,p);
	s += getMag2(v1); 
  }

  timer.Stop();

  double sav = std::sqrt(s/double(n));

  std::cout << " Time for Random gen " << timer.RealTime() << "  " << timer.CpuTime() << std::endl; 
  int pr = std::cout.precision(18);  std::cout << "Average : " << sav << std::endl;   std::cout.precision(pr);

  return sav; 
}

//----------------------------------------------------------------
/// writing
//----------------------------------------------------------------

template<class Vector> 
double write(int n, const std::string & file_name, const std::string & vector_type, int compress = 0) { 

  TStopwatch timer;

  TRandom3 R(111);

  std::cout << "writing a tree with " << vector_type << std::endl; 

  std::string fname = file_name + ".root";
  TFile f1(fname.c_str(),"RECREATE","",compress);

  // create tree
  std::string tree_name="Tree with" + vector_type; 
  TTree t1("t1",tree_name.c_str());

  Vector *v1 = new Vector();
  std::cout << "typeID written : " << typeid(*v1).name() << std::endl;

  t1.Branch("Vector branch",vector_type.c_str(),&v1);

  timer.Start();
  double p[4]; 
  double s = 0; 
  for (int i = 0; i < n; ++i) { 
        p[0] = R.Gaus(0,10);
	p[1] = R.Gaus(0,10);
	p[2] = R.Gaus(0,10);
	p[3] = R.Gaus(100,10);
	//CylindricalEta4D<double> & c = v1->Coordinates();
	//c.SetValues(Px,pY,pZ,E);
	setValues(*v1,p);
	t1.Fill();
	s += getMag2(*v1); 
  }

  f1.Write();
  timer.Stop();

  double sav = std::sqrt(s/double(n));


#ifdef DEBUG
  t1.Print();
  std::cout << " Time for Writing " << file_name << "\t: " << timer.RealTime() << "  " << timer.CpuTime() << std::endl; 
  int pr = std::cout.precision(18);  std::cout << "Average : " << sav << std::endl;   std::cout.precision(pr);
#endif
  
  return sav; 
}


//----------------------------------------------------------------
/// reading
//----------------------------------------------------------------

template<class Vector> 
double read(const std::string & file_name) { 

  TStopwatch timer;

  std::string fname = file_name + ".root";

  TFile f1(fname.c_str());
  if (f1.IsZombie() ) { 
    std::cout << " Error opening file " << file_name << std::endl; 
    return -1; 
  }

  //TFile f1("mathcoreVectorIO_D32.root");

  // create tree
  TTree *t1 = (TTree*)f1.Get("t1");

  Vector *v1 = 0;

  std::cout << "reading typeID  : " << typeid(*v1).name() << std::endl;

  t1->SetBranchAddress("Vector branch",&v1);

  timer.Start();
  int n = (int) t1->GetEntries();
  std::cout << " Tree Entries " << n << std::endl; 
  double s=0;
  for (int i = 0; i < n; ++i) { 
    t1->GetEntry(i);
    s += getMag2(*v1);
  }


  timer.Stop();

  double sav = std::sqrt(s/double(n));

#ifdef DEBUG
  std::cout << " Time for Reading " << file_name << "\t: " << timer.RealTime() << "  " << timer.CpuTime() << std::endl; 
  int pr = std::cout.precision(18);  std::cout << "Average : " << sav << std::endl;   std::cout.precision(pr);
#endif

  return sav; 
}

template<class TrackType>
double writeTrack(int n, const std::string & file_name, int compress = 0) { 

  std::cout << "\n";
  std::cout << "  Test writing .." << file_name << " .....\n";
  std::cout << "**************************************************\n";

  TRandom3 R(111);

  TStopwatch timer;

  //std::cout << "writing a tree with " << vector_type << std::endl; 

  std::string fname = file_name + ".root";
  TFile f1(fname.c_str(),"RECREATE","",compress);

  // create tree
  std::string tree_name="Tree with TrackD"; 
  TTree t1("t1",tree_name.c_str());

  TrackType *track = new TrackType();
  std::cout << "typeID written : " << typeid(*track).name() << std::endl;

  
  //t1.Branch("Vector branch",&track,16000,0);
  // default split model
  t1.Branch("Vector branch",&track);

  timer.Start();
  double s = 0; 
  ROOT::Math::XYZTVector q; 
  ROOT::Math::XYZPoint p; 
  typedef typename TrackType::VectorType V;
  typedef typename TrackType::PointType  P;

  for (int i = 0; i < n; ++i) { 
    q.SetXYZT( R.Gaus(0,10),
	       R.Gaus(0,10),
	       R.Gaus(0,10),
	       R.Gaus(100,10) ); 
    p.SetXYZ( q.X(), q.Y(), q.Z() ); 
	
    track->Set( V(q), P(p) ); 
    
    t1.Fill();
    s += getMag2( *track ); 
  }

  f1.Write();
  timer.Stop();

  double sav = std::sqrt(s/double(n));


#ifdef DEBUG
  t1.Print();
  std::cout << " Time for Writing " << file_name << "\t: " << timer.RealTime() << "  " << timer.CpuTime() << std::endl; 
  int pr = std::cout.precision(18);  std::cout << "Average : " << sav << std::endl;   std::cout.precision(pr);
#endif
  
  return sav; 

}



int testResult(double w1, double r1, const std::string & type) { 

  int iret = 0; 
  std::cout << w1 << "  r " << r1 << std::endl; 
  // do like this to avoid nan's
  if (!( fabs(w1-r1)  < tol)) { 
    std::cout << "\nERROR: Differeces found  when reading " << std::endl;
    int pr = std::cout.precision(18);  std::cout << w1 << "   !=    " << r1 << std::endl; std::cout.precision(pr);
    iret = -1;
  }

  std::cout << "\n*********************************************************************************************\n"; 
  std::cout << "Test :\t " << type << "\t\t";
  if (iret ==0) 
    std::cout << "OK" << std::endl; 
  else 
    std::cout << "FAILED" << std::endl; 
  std::cout << "********************************************************************************************\n\n"; 

  return iret; 
}



int testVectorIO(bool readOnly = false) { 

  int iret = 0; 

// #ifdef __CINT__
//   gSystem->Load("libMathCore");  
//   gSystem->Load("libPhysics");  
//   using namespace ROOT::Math;
// #endif

#ifdef USE_REFLEX

// #ifdef __CINT__
//   gSystem->Load("libReflex");  
//   gSystem->Load("libCintex");  
// #endif

  std::cout << "Use Reflex dictionary " << std::endl; 

  gSystem->Load("libReflex");  
  gSystem->Load("libCintex");  

  ROOT::Cintex::Cintex::SetDebug(1);
  ROOT::Cintex::Cintex::Enable();



  //iret |= gSystem->Load("libSmatrixRflx");  
  //iret |= gSystem->Load("libMathAddRflx");  
  //iret |= gSystem->Load("libMathRflx");  
  if (iret |= 0) { 
    std::cerr <<"Failing to Load Reflex dictionaries " << std::endl;
    return iret; 
  }

#endif   // endif on using reflex

  
  int nEvents = 100000;
  //int nEvents = 100;

  double w1, r1 = 0;
  std::string fname; 
  
  testDummy<AVector4D>(nEvents);
 
  fname = "lorentzvector";
  if (readOnly) {
     w1 = 98.995527276968474;
     fname += "_prev";
  }
  else
     w1 = write<AVector4D>(nEvents,fname,vector4d_type);

  r1 = read<AVector4D>(fname);
  iret |= testResult(w1,r1,vector4d_type); 

  fname = "displacementvector";
  if (readOnly) {
     w1 = 17.3281570633214095; 
     fname += "_prev";
  }
  else 
     w1 = write<AVector3D>(nEvents,fname,vector3d_type);

  r1 = read<AVector3D>(fname);
  iret |= testResult(w1,r1,vector3d_type); 

  fname = "positionvector";
  if (readOnly)
     fname += "_prev";
  else 
     w1 = write<APoint3D>(nEvents,fname,point3d_type);
  
  r1 = read<APoint3D>(fname);
  iret |= testResult(w1,r1,point3d_type); 


  // test TrackD
  fname = "track";
  // load track dictionary 
  gSystem->Load("libTrackDict");  

  //nEvents = 10000; 

  if (readOnly) {
     fname += "_prev";
  }
  else 
    w1 = writeTrack<TrackD>( nEvents,fname); 
  
  r1 = read<TrackD>(fname);
  iret |= testResult(w1,r1,"TrackD"); 

  // test TrackD32

  fname = "trackD32";
  // load track dictionary 
  //  gSystem->Load("libTrackDict");  

  if (readOnly) {
     fname += "_prev";
  }
  else 
    w1 = writeTrack<TrackD32>( nEvents,fname); 
  
  r1 = read<TrackD32>(fname);
  iret |= testResult(w1,r1,"TrackD32"); 


  // test vector of tracks
  fname = "vectrack";

  nEvents = 10000; 

  if (readOnly) {
     fname += "_prev";
  }
  else 
    w1 = writeTrack<VecTrackD>( nEvents,fname); 
  
  r1 = read<VecTrackD>(fname);
  iret |= testResult(w1,r1,"VecTrackD"); 

  // test ClusterD (cotaining a vector of points) 
  fname = "cluster";

  nEvents = 10000; 

  if (readOnly) {
     fname += "_prev";
  }
  else 
    w1 = writeTrack<ClusterD>( nEvents,fname); 
  
  r1 = read<ClusterD>(fname);
  iret |= testResult(w1,r1,"ClusterD"); 
  
  return iret; 
}

int main(int argc, char ** ) { 

   bool readOnly = false;
   if (argc > 1) readOnly = true; 

   int iret =  testVectorIO(readOnly);
   if (iret != 0) 
      std::cerr << "testVectorIO:\t FAILED ! " << std::endl;
   else 
      std::cerr << "testVectorIO:\t OK ! " << std::endl;
   
   return iret; 

}