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

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


#include "fmath.h"
#include "nbrandom.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;

  try {
  POutPersist so("tobjio2.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";
  cout << "Writing " << nom << endl;
  so.PutObject(*nt, nom);

  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";
  cout << "Writing " << nom << endl;
  so.PutObject(*nt, nom);

  {  //  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";
  cout << "Writing " << nom << endl;
  so.PutObject(*h, nom);
  }


  {    // Matrices 
  TArray<double> * mtx = new TArray<double>(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,0) = 11.*exp(-rp);
     xp = i-35;   yp = j-15;
     rp = (xp*xp)/9.+(yp*yp)/9.;
     (*mtx)(j,i,0) += 3.*exp(-rp);
     xp = i-45;   yp = j-40;
     rp = (xp*xp)/10.+(yp*yp)/20.;
     (*mtx)(j,i,0) += 6.*exp(-rp);
     }
  nom = "mtx1";
  cout << "Writing " << nom << endl;
  so.PutObject(*mtx, nom);

  }
  }
  catch (PThrowable exc) {
    cerr << " Catched Exception PThrowable - Msg= " << exc.Msg() << endl;
  }
  cout << " End --- exit .... " << endl;
}