#include <stdio.h>
#include <stdlib.h>
#include "sopnamsp.h"
#include "math.h"
#include <iostream>
#include <string>

// Pour le test 
#include "histinit.h"
#include "histos.h"
#include "histos2.h"
#include "ntuple.h"


#include "fmath.h"
#include "srandgen.h"
#include "fioarr.h"


int main(int narg, char* arg[])
{

  SophyaInit();

  cout << "Hash check, OMatrix : " << hex << PIOPersist::Hash("OMatrix") << dec << endl;

  cout << "Test Creation / PPersist Save NTuple, Histo, ..." << endl;

  POutPersist so("tobjio.ppf");

  NTuple *nt;
  float xnt[4];
  char *ntn[4] = {"x","y","ex","ey"};
  int nent,i;
  nent = 1000;
  nt = new NTuple(2,ntn);  // Creation NTuple (AVEC new )
  for(i=0; i<nent; i++) {
    xnt[0] = i*50./nent; xnt[1] = 0.2*xnt[0]*xnt[0]+(NorRand()*xnt[0]/3.);
    nt->Fill(xnt);
  }

  string nom = "nt21";
  ObjFileIO<NTuple> oiont(nt);
  oiont.Write(so, nom);
  cout << "Writing " << nom << endl;

  string nx, ny, nz;
  string ex,ey,ez;
  nx = "x";  ny = "y";
  nent = 20;
  nt = new NTuple(4,ntn);  // Creation NTuple (AVEC new )
  for(i=0; i<nent; i++) {
    xnt[0] = i*50./nent; xnt[1] = 0.2*xnt[0]*xnt[0]+(NorRand()*xnt[0]/3.);
    xnt[2] = sqrt((double) xnt[0]); xnt[3] = xnt[0]*2;
    nt->Fill(xnt);
  }
  nom = "nt22";
  ObjFileIO<NTuple> oiont2(nt);
  oiont2.Write(so, nom);
  cout << "Writing " << nom << endl;

  {  //  Objet Histo-2D
  r_8 x,y,w;
  int i,j;
  Histo2D* h = new Histo2D(-10.,10.,25,-10.,10.,25);
  for(i=0; i<25; i++) for(j=0; j<25; j++) {
    h->BinCenter(i,j,x,y);
    w = 100.*exp(-0.5*(x*x/9. + y*y/9.));
    h->Add(x,y,w);
    w = 50.*exp(-0.5*((x+10.)*(x+10.)/9. + (y+10.)*(y+10.)/9.));
    h->Add(x,y,w);
    }
  nom = "h2d";
  ObjFileIO<Histo2D> oioh(h);
  oioh.Write(so, nom);
  cout << "Writing " << nom << endl;
  }

  {
  Histo *h;
  float x;


  h = new Histo(0., 200., 100);  // Creation histo (AVEC new )
  for(i=0; i<100; i++)  // Remplissage d'histo
    { x = (2*i+1.);  h->Add(x, x*(200.-x)); }
  nom = "h1d";
  ObjFileIO<Histo> oioh(h);
  oioh.Write(so, nom);
  cout << "Writing " << nom << endl;
  
  }

  {   // Vecteurs 
  Vector* vpx, *vpy;
  vpx = new Vector(20);  vpy = new Vector(20);
  int i;   double x;
  for(i=0; i<20; i++) { x = 8.*i/20;  (*vpx)(i) = x;  (*vpy)(i) = 3*cos(x)+2*sin(x); }
  nom = "vx1";
  FIO_TArray<double>  oiov1(vpx);
  oiov1.Write(so, nom);
  cout << "Writing " << nom << endl;
  nom = "vy1";
  FIO_TArray<double>  oiov2(vpy);
  oiov2.Write(so, nom);
  cout << "Writing " << nom << endl;
  vpx = new Vector(30);  
  for(i=0; i<30; i++) { x = 8.*i/30+0.1;  (*vpx)(i) = sin(2.*x)/(2.*x); }
  nom = "vx2";
  FIO_TArray<double>  oiov3(vpx);
  oiov3.Write(so, nom);
  cout << "Writing " << nom << endl;
  }

      
  {       // NTuple 3D
  float xnt[6];
  double ang;
  char *ntn[6] = {"x","y","z","ex","ey","ez"};
  int nent,i;
      
  NTuple* nt = new NTuple(6,ntn);  // Creation NTuple (AVEC new )
  nent = 3000;
  for(i=0; i<nent; i++) {
  ang = 5.*6.28*(double)i/(double)nent;
  xnt[0] = cos(ang)*1.5+NorRand()*0.1; xnt[1] = sin(ang)*1.5+NorRand()*0.1;
  xnt[2] = ang/8.+NorRand()*0.05;   xnt[3] = xnt[4] = xnt[5] = 0.;
  nt->Fill(xnt);
  }

  nom = "nt31";
  ObjFileIO<NTuple> oiont(nt);
  oiont.Write(so, nom);
  cout << "Writing " << nom << endl;
  
  nt = new NTuple(6,ntn);  // Creation NTuple (AVEC new )
  nent = 100;
  for(i=0; i<nent; i++) {
    ang = 5.*6.28*(double)i/(double)nent;
    xnt[0] = cos(-ang)*1+NorRand()*0.0; xnt[1] = sin(-ang)*1.+NorRand()*0.0;
    xnt[2] = ang/8.+NorRand()*0.0;   xnt[3] = 0.1; xnt[4] = 0.15;  xnt[5] = 0.07;
    nt->Fill(xnt);
    }
  nom = "nt32";
  ObjFileIO<NTuple> oiont2(nt);
  oiont2.Write(so, nom);
  cout << "Writing " << nom << endl;
  }

  {    // Matrices 
  Matrix* mtx = new Matrix(50, 50);
  int i,j;
  double xp,yp,rp;
  for(i=0; i<50; i++)
     for(j=0; j<50; j++) { 
     xp = i-15;   yp = j-25;
     rp = (xp*xp)/16.+(yp*yp)/16.;
     (*mtx)(j,i) = 11.*exp(-rp);
     xp = i-35;   yp = j-15;
     rp = (xp*xp)/9.+(yp*yp)/9.;
     (*mtx)(j,i) += 3.*exp(-rp);
     xp = i-45;   yp = j-40;
     rp = (xp*xp)/10.+(yp*yp)/20.;
     (*mtx)(j,i) += 6.*exp(-rp);
     }
  nom = "mtx1";
  FIO_TArray<double> oiom(mtx);
  oiom.Write(so, nom);
  cout << "Writing " << nom << endl;
  }

  cout << " End --- exit .... " << endl;
}