source: Sophya/trunk/SophyaLib/HiStats/histos2.cc@ 1074

//
//    cmv 05/08/96
//

#include "machdefs.h"

#include <string.h>
#include <new.h>
#include <stdlib.h>
#include <stdio.h>
#ifndef __mac__
#include <values.h>
#endif

#include "histos2.h"
#include "generalfit.h"

/*!
  \class SOPHYA::Histo2D
  \ingroup HiStats
  Classe d'histogrammes 2D
  \verbatim
    Remarque sur les indices:
     H(i,j)   -> i  = coord x (0<i<nx),     j  = coord y (0<j<ny)
     v(ii,jj) -> ii = ligne (0<i<NRows()),  jj = colonne (0<i<NCol())
     On fait une correspondance directe i<->ii et j<->jj
     ce qui, en representation classique des histos2D et des matrices
     entraine une inversion x<->y cad une symetrie / diagonale principale
     H(0,...)  represente ^         mais v(0,...) represente
                          |x.......               |xxxxxxxx|
                          |x.......               |........|
                          |x.......               |........|
                          |x.......               |........|
                          |x.......               |........|
                          --------->
                          colonne no 1            ligne no 1
  \endverbatim
*/

///////////////////////////////////////////////////////////////////
/*!
  Createur d'un histogramme 2D ayant nxBin,nyBin bins
  entre xMin,xMax et yMin,yMax.
*/
Histo2D::Histo2D(float xMin, float xMax, int nxBin
                ,float yMin, float yMax, int nyBin)
      : data(new float[nxBin*nyBin]), err2(NULL)
      , nHist(0), nEntries(0)
      , nx(nxBin), ny(nyBin), nxy(nxBin*nyBin)
      , xmin(xMin), xmax(xMax), ymin(yMin), ymax(yMax)
      , wbinx((xMax - xMin)/nxBin), wbiny((yMax - yMin)/nyBin)
      , hprojx(NULL), hprojy(NULL)
{
ASSERT(nxBin>0 && nyBin>0 && xMin<xMax && yMin<yMax);
for(int i=0;i<3;i++) for(int j=0;j<3;j++) over[i][j]=0.;
Zero();
b_s.H = NULL;
END_CONSTRUCTOR
}

/*!
  Constructeur par copie.
*/
Histo2D::Histo2D(const Histo2D& h)
{
int i,j;
data = new float[h.nxy];
memcpy(data, h.data, h.nxy*sizeof(float));

err2 = NULL;
if(h.err2) {
  err2 = new double[h.nxy];
  memcpy(err2, h.err2, h.nxy*sizeof(double));
}

nHist = h.nHist; nEntries = h.nEntries;
for(i=0;i<3;i++) for(j=0;j<3;j++) over[i][j]=h.over[i][j];
nx = h.nx;  ny = h.ny;  nxy = h.nxy;
xmin = h.xmin; xmax = h.xmax; ymin = h.ymin; ymax = h.ymax;
wbinx = h.wbinx; wbiny = h.wbiny;
b_s.H = NULL;

hprojx = hprojy = NULL;
if(h.hprojx) {
  SetProjX();
  *hprojx = *(h.hprojx);
}
if(h.hprojy) {
  SetProjY();
  *hprojy = *(h.hprojy); 
}

int nb;
float min,max;
nb = h.NSliX();
if(nb>0) {
  SetSliX(nb);
  for(i=0; i<NSliX();i++) *HSliX(i) = *(h.HSliX(i));
}
nb = h.NSliY();
if(nb>0) {
  SetSliY(nb);
  for(i=0; i<NSliY();i++) *HSliY(i) = *(h.HSliY(i));
}

nb = h.NBandX();
if(nb>0) {
  for(i=0; i<nb;i++) {
    h.GetBandX(i,min,max);
    SetBandX(min,max);
    *HBandX(i) = *(h.HBandX(i));
  }
  for(i=0; i<NBandX();i++) *HBandX(i) = *(h.HBandX(i));
}
nb = h.NBandY();
if(nb>0) {
  for(i=0; i<nb;i++) {
    h.GetBandY(i,min,max);
    SetBandY(min,max);
    *HBandY(i) = *(h.HBandY(i));
  }
  for(i=0; i<NBandY();i++) *HBandY(i) = *(h.HBandY(i));
}

END_CONSTRUCTOR
}

/*!
  Constructeur par defaut.
*/
Histo2D::Histo2D()
      : data(NULL), err2(NULL)
      , nHist(0), nEntries(0)
      , nx(0), ny(0), nxy(0)
      , xmin(0), xmax(0), ymin(0), ymax(0)
      , wbinx(0), wbiny(0)
      , hprojx(NULL), hprojy(NULL)
{
for(int i=0;i<3;i++) for(int j=0;j<3;j++) over[i][j]=0.;
b_s.H = NULL;
END_CONSTRUCTOR
}

///////////////////////////////////////////////////////////////////
/*!
  Desallocation de la place de l'histogramme (fct privee).
*/
void Histo2D::Delete()
{
  if( data != NULL ) { delete[] data; data = NULL;}

  if( err2 != NULL ) { delete[] err2; err2 = NULL;}

  DelProj();

  DelBandX();
  DelBandY();

  DelSliX();
  DelSliY();

  nHist = 0;
  nEntries = 0;
  nx = 0; ny = 0; nxy = 0;
  xmin = 0; xmax = 0; ymin = 0; ymax = 0;
  wbinx = 0; wbiny = 0;
  for(int i=0;i<3;i++) for(int j=0;j<3;j++) over[i][j]=0.;
  b_s.H = NULL;
}

/*!
  Destructeur.
*/
Histo2D::~Histo2D()
{
Delete();
}

///////////////////////////////////////////////////////////////////
/*!
  Remise a zero du contenu, des erreurs et des valeurs.
*/
void Histo2D::Zero()
{
  nHist = nEntries = 0;
  for(int i=0;i<3;i++) for(int j=0;j<3;j++) over[i][j]=0.;
  memset(data, 0, nxy*sizeof(float));
  memset(over, 0, 9*sizeof(float));

  if( err2 != NULL ) memset(err2, 0, nxy*sizeof(double));

  ZeroProj();

  ZeroBandX();
  ZeroBandY();

  ZeroSliX();
  ZeroSliY();
}

///////////////////////////////////////////////////////////////////
/*!
  Pour avoir le calcul des erreurs.
*/
void Histo2D::Errors()
{
 if( nxy > 0 ) {
   if(err2==NULL) err2 = new double[nxy];
   memset(err2, 0, nxy*sizeof(double));
 }
}

///////////////////////////////////////////////////////////////////
/*!
  Operateur H2 = H1
*/
Histo2D& Histo2D::operator = (const Histo2D& h)
{
  int i,j,nb;
  float min,max;

  if(this == &h) return *this;
  if( h.nxy > nxy ) Delete();
  if(!data) data = new float[h.nxy];
  if( !h.err2 && err2 ) { delete [] err2; err2=NULL;}
  if( h.err2 && !err2 ) err2 = new double[h.nxy];

  for(i=0;i<3;i++) for(j=0;j<3;j++) over[i][j] = h.over[i][j];
  nHist = h.nHist;
  nEntries = h.nEntries;
  nx = h.nx;  ny = h.ny;  nxy = h.nxy;
  xmin = h.xmin; xmax = h.xmax; wbinx = h.wbinx;
  ymin = h.ymin; ymax = h.ymax; wbiny = h.wbiny;
  
  memcpy(data, h.data, nxy*sizeof(float));
  if(err2) memcpy(err2, h.err2, nxy*sizeof(double));

  DelProjX();
  if(h.hprojx) {
    SetProjX();
    *hprojx = *(h.hprojx);
  }
  DelProjY();
  if(h.hprojy) {
    SetProjY();
    *hprojy = *(h.hprojy);
  }

  DelSliX();
  nb = h.NSliX();
  if(nb>0) {
    SetSliX(nb);
    for(i=0; i<NSliX();i++) *HSliX(i) = *(h.HSliX(i));
  }
  DelSliY();
  nb = h.NSliY();
  if(nb>0) {
    SetSliY(nb);
    for(i=0; i<NSliY();i++) *HSliY(i) = *(h.HSliY(i));
  }

  DelBandX();
  nb = h.NBandX();
  if(nb>0) {
    for(i=0; i<nb;i++) {
      h.GetBandX(i,min,max);
      SetBandX(min,max);
      *HBandX(i) = *(h.HBandX(i));
    }
    for(i=0; i<NBandX();i++) *HBandX(i) = *(h.HBandX(i));
  }
  DelBandY();
  nb = h.NBandY();
  if(nb>0) {
    for(i=0; i<nb;i++) {
      h.GetBandY(i,min,max);
      SetBandY(min,max);
      *HBandY(i) = *(h.HBandY(i));
    }
    for(i=0; i<NBandY();i++) *HBandY(i) = *(h.HBandY(i));
  }

  return *this;
}

///////////////////////////////////////////////////////////////////
/*!
  Operateur H *= b
*/
Histo2D& Histo2D::operator *= (double b)
{
int i,j;
double b2 = b*b;
for(i=0;i<nxy;i++) {
  data[i] *= b;
  if(err2) err2[i] *= b2;
}
for(i=0;i<3;i++) for(j=0;j<3;j++) over[i][j] *= b;
nHist *= b;

if(hprojx) *hprojx *= b;
if(hprojy) *hprojy *= b;
if(NSliX()>0) for(i=0; i<NSliX();i++) *HSliX(i) *= b;
if(NSliY()>0) for(i=0; i<NSliY();i++) *HSliY(i) *= b;
if(NBandX()>0) for(i=0; i<NBandX();i++) *HBandX(i) *= b;
if(NBandY()>0) for(i=0; i<NBandY();i++) *HBandY(i) *= b;

return *this;
}

/*!
  Operateur H /= b
*/
Histo2D& Histo2D::operator /= (double b)
{
int i,j;
if (b==0.) THROW(inconsistentErr);
double b2 = b*b;
for(i=0;i<nxy;i++) {
  data[i] /= b;
  if(err2) err2[i] /= b2;
}
for(i=0;i<3;i++) for(j=0;j<3;j++) over[i][j] /= b;
nHist /= b;

if(hprojx) *hprojx /= b;
if(hprojy) *hprojy /= b;
if(NSliX()>0) for(i=0; i<NSliX();i++) *HSliX(i) /= b;
if(NSliY()>0) for(i=0; i<NSliY();i++) *HSliY(i) /= b;
if(NBandX()>0) for(i=0; i<NBandX();i++) *HBandX(i) /= b;
if(NBandY()>0) for(i=0; i<NBandY();i++) *HBandY(i) /= b;

return *this;
}

/*!
  Operateur H += b
*/
Histo2D& Histo2D::operator += (double b)
{
int i,j;
float min,max;
for(i=0;i<nxy;i++) data[i] += b;
for(i=0;i<3;i++) for(j=0;j<3;j++) over[i][j] += b;
nHist += nxy*b;

if(hprojx) *hprojx += b*ny;
if(hprojy) *hprojy += b*nx;
if(NSliX()>0) for(i=0; i<NSliX();i++) *HSliX(i) += b*ny/NSliX();
if(NSliY()>0) for(i=0; i<NSliY();i++) *HSliY(i) += b*nx/NSliY();
if(NBandX()>0) for(i=0; i<NBandX();i++) {
  GetBandX(i,min,max);
  *HBandX(i) += b*(max-min)/(ymax-ymin)*ny;
}
if(NBandY()>0) for(i=0; i<NBandY();i++) {
  GetBandY(i,min,max);
  *HBandY(i) += b*(max-min)/(xmax-xmin)*nx;
}

return *this;
}

/*!
  Operateur H -= b
*/
Histo2D& Histo2D::operator -= (double b)
{
int i,j;
float min,max;
for(i=0;i<nxy;i++) data[i] -= b;
for(i=0;i<3;i++) for(j=0;j<3;j++) over[i][j] -= b;
nHist -= nxy*b;

if(hprojx) *hprojx -= b*ny;
if(hprojy) *hprojy -= b*nx;
if(NSliX()>0) for(i=0; i<NSliX();i++) *HSliX(i) -= b*ny/NSliX();
if(NSliY()>0) for(i=0; i<NSliY();i++) *HSliY(i) -= b*nx/NSliY();
if(NBandX()>0) for(i=0; i<NBandX();i++) {
  GetBandX(i,min,max);
  *HBandX(i) -= b*(max-min)/(ymax-ymin)*ny;
}
if(NBandY()>0) for(i=0; i<NBandY();i++) {
  GetBandY(i,min,max);
  *HBandY(i) -= b*(max-min)/(xmax-xmin)*nx;
}

return *this;
}

///////////////////////////////////////////////////////////////////
/*!
  Operateur H += H1
*/
Histo2D& Histo2D::operator += (const Histo2D& a)
{
int i,j;
if(nx!=a.nx || ny!=a.ny) THROW(sizeMismatchErr);
for(i=0;i<nxy;i++) {
  data[i] += a.data[i];
  if(err2 && a.err2) err2[i] += a.err2[i];
}
for(i=0;i<3;i++) for(j=0;j<3;j++) over[i][j] += a.over[i][j];
nHist += a.nHist;
nEntries += a.nEntries;

if(hprojx && a.hprojx) *hprojx += *(a.hprojx);
if(hprojy && a.hprojy) *hprojy += *(a.hprojy);
ZeroSliX();  ZeroSliY();
ZeroBandX(); ZeroBandY();

return *this;
}

/*!
  Operateur H -= H1
*/
Histo2D& Histo2D::operator -= (const Histo2D& a)
{
int i,j;
if(nx!=a.nx || ny!=a.ny) THROW(sizeMismatchErr);
for(i=0;i<nxy;i++) {
  data[i] -= a.data[i];
  if(err2 && a.err2) err2[i] += a.err2[i];
}
for(i=0;i<3;i++) for(j=0;j<3;j++) over[i][j] += a.over[i][j];
nHist -= a.nHist;
nEntries += a.nEntries;

if(hprojx && a.hprojx) *hprojx -= *(a.hprojx);
if(hprojy && a.hprojy) *hprojy -= *(a.hprojy);
ZeroSliX();  ZeroSliY();
ZeroBandX(); ZeroBandY();

return *this;
}

/*!
  Operateur H *= H1
*/
Histo2D& Histo2D::operator *= (const Histo2D& a)
{
int i,j;
if(nx!=a.nx || ny!=a.ny) THROW(sizeMismatchErr);
nHist = 0.;
for(i=0;i<nxy;i++) {
  if(err2 && a.err2)
      err2[i] = a.data[i]*a.data[i]*err2[i] + data[i]*data[i]*a.err2[i];
  data[i] *= a.data[i];
  nHist += data[i];
}
for(i=0;i<3;i++) for(j=0;j<3;j++) over[i][j] *= a.over[i][j];
nEntries += a.nEntries;

if(hprojx && a.hprojx) *hprojx *= *(a.hprojx);
if(hprojy && a.hprojy) *hprojy *= *(a.hprojy);
ZeroSliX();  ZeroSliY();
ZeroBandX(); ZeroBandY();

return *this;
}

/*!
  Operateur H /= H1
*/
Histo2D& Histo2D::operator /= (const Histo2D& a)
{
int i,j;
if(nx!=a.nx || ny!=a.ny) THROW(sizeMismatchErr);
nHist = 0.;
for(i=0;i<nxy;i++) {
  if(a.data[i]==0.) {
    data[i]=0.;
    if(err2) err2[i]=0.;
    continue;
  }
  if(err2 && a.err2)
      err2[i] = (err2[i] + data[i]/a.data[i]*data[i]/a.data[i]*a.err2[i])
                /(a.data[i]*a.data[i]);
  data[i] /= a.data[i];
  nHist += data[i];
}
for(i=0;i<3;i++) for(j=0;j<3;j++) 
  if(a.over[i][j]!=0.) over[i][j] *= a.over[i][j]; else over[i][j] = 0.;
nEntries += a.nEntries;

if(hprojx && a.hprojx) *hprojx /= *(a.hprojx);
if(hprojy && a.hprojy) *hprojy /= *(a.hprojy);
ZeroSliX();  ZeroSliY();
ZeroBandX(); ZeroBandY();

return *this;
}

///////////////////////////////////////////////////////////////////
/*!
  Remplissage d'un tableau avec les valeurs des abscisses.
*/
void Histo2D::GetXCoor(TVector<r_8> &v)
{
float x,y;
v.Realloc(nx);
for(int i=0;i<nx;i++) {BinLowEdge(i,0,x,y); v(i) = x;}
return;
}

/*!
  Remplissage d'un tableau avec les valeurs des ordonnees.
*/
void Histo2D::GetYCoor(TVector<r_8> &v)
{
float x,y;
v.Realloc(ny);
for(int i=0;i<ny;i++) {BinLowEdge(0,i,x,y); v(i) = y;}
return;
}

/*!
  Remplissage d'un tableau avec les valeurs du contenu.
*/
void Histo2D::GetValue(TMatrix<r_8> &v)
{
v.Realloc(nx,ny);
for(int i=0;i<nx;i++)
  for(int j=0;j<ny;j++) v(i,j) = (*this)(i,j);
return;
}

/*!
  Remplissage d'un tableau avec les valeurs du carre des erreurs.
*/
void Histo2D::GetError2(TMatrix<r_8> &v)
{
int i,j;
v.Realloc(nx,ny);
if(!err2)
  {for(i=0;i<nx;i++) for(j=0;j<ny;j++) v(i,j) = 0.; return;}
for(i=0;i<nx;i++) for(j=0;j<ny;j++) v(i,j) = Error2(i,j);
return;
}

/*!
  Remplissage d'un tableau avec les valeurs des erreurs.
*/
void Histo2D::GetError(TMatrix<r_8> &v)
{
int i,j;
v.Realloc(nx,ny);
if(!err2)
  {for(i=0;i<nx;i++) for(j=0;j<ny;j++) v(i,j) = 0.; return;}
for(i=0;i<nx;i++) for(j=0;j<ny;j++) v(i,j) = Error(i,j);
return;
}

///////////////////////////////////////////////////////////////////
/*!
  Remplissage du contenu de l'histo avec les valeurs d'un tableau.
*/
void Histo2D::PutValue(TMatrix<r_8> &v, int ierr)
{
int i,j;
//if(v.NRows()!=(uint_4)nx || v.NCol()!=(uint_4)ny) THROW(sizeMismatchErr);
uint_4 nnx = (v.NRows()<(uint_4)nx)? v.NRows(): (uint_4)nx;
uint_4 nny = (v.NCol() <(uint_4)ny)? v.NCol() : (uint_4)ny;
if(nnx>0 && nny>0) for(i=0;i<nnx;i++) for(j=0;j<nny;j++) {
  (*this)(i,j) = v(i,j);
  if(err2 && ierr) Error2(i,j) = fabs(v(i,j));
}
return;
}

/*!
  Addition du contenu de l'histo avec les valeurs d'un tableau.
*/
void Histo2D::PutValueAdd(TMatrix<r_8> &v, int ierr)
{
int i,j;
//if(v.NRows()!=(uint_4)nx || v.NCol()!=(uint_4)ny) THROW(sizeMismatchErr);
uint_4 nnx = (v.NRows()<(uint_4)nx)? v.NRows(): (uint_4)nx;
uint_4 nny = (v.NCol() <(uint_4)ny)? v.NCol() : (uint_4)ny;
if(nnx>0 && nny>0) for(i=0;i<nnx;i++) for(j=0;j<nny;j++) {
  (*this)(i,j) += v(i,j);
  if(err2 && ierr) Error2(i,j) += fabs(v(i,j));
}
return;
}

/*!
  Remplissage des erreurs au carre de l'histo
  avec les valeurs d'un tableau.
*/
void Histo2D::PutError2(TMatrix<r_8> &v)
{
int i,j;
//if(v.NRows()!=(uint_4)nx || v.NCol()!=(uint_4)ny) THROW(sizeMismatchErr);
uint_4 nnx = (v.NRows()<(uint_4)nx)? v.NRows(): (uint_4)nx;
uint_4 nny = (v.NCol() <(uint_4)ny)? v.NCol() : (uint_4)ny;
if(nnx>0 && nny>0) {
  if(!err2) Errors();
  for(i=0;i<nnx;i++) for(j=0;j<nny;j++) Error2(i,j) = v(i,j);
}
return;
}

/*!
  Addition des erreurs au carre de l'histo
  avec les valeurs d'un tableau.
*/
void Histo2D::PutError2Add(TMatrix<r_8> &v)
{
int i,j;
//if(v.NRows()!=(uint_4)nx || v.NCol()!=(uint_4)ny) THROW(sizeMismatchErr);
uint_4 nnx = (v.NRows()<(uint_4)nx)? v.NRows(): (uint_4)nx;
uint_4 nny = (v.NCol() <(uint_4)ny)? v.NCol() : (uint_4)ny;
if(nnx>0 && nny>0) {
  if(!err2) Errors();
  for(i=0;i<nnx;i++) for(j=0;j<nny;j++)
          if(v(i,j)>0.) Error2(i,j) += v(i,j);
}
return;
}

/*!
  Remplissage des erreurs de l'histo avec les valeurs d'un tableau.
*/
void Histo2D::PutError(TMatrix<r_8> &v)
{
int i,j;
//if(v.NRows()!=(uint_4)nx || v.NCol()!=(uint_4)ny) THROW(sizeMismatchErr);
uint_4 nnx = (v.NRows()<(uint_4)nx)? v.NRows(): (uint_4)nx;
uint_4 nny = (v.NCol() <(uint_4)ny)? v.NCol() : (uint_4)ny;
if(nnx>0 && nny>0) {
  if(!err2) Errors();
  for(i=0;i<nnx;i++) for(j=0;j<nny;j++)
    if(v(i,j)>0.) Error2(i,j)=v(i,j)*v(i,j); else Error2(i,j)= -v(i,j)*v(i,j);
}
return;
}

///////////////////////////////////////////////////////////////////
/********* Methode *********/
/*!
  Addition du contenu de l'histo pour x,y poids w.
*/
void Histo2D::Add(float x, float y, float w)
{
list<bande_slice>::iterator it;
int i,j;
FindBin(x,y,i,j);

if( hprojx != NULL ) hprojx->Add(x,w);
if( hprojy != NULL ) hprojy->Add(y,w);

if(lbandx.size()>0)
  for( it = lbandx.begin(); it != lbandx.end(); it++)
    if( (*it).min <= y && y < (*it).max ) (*it).H->Add(x,w);

if(lbandy.size()>0)
  for( it = lbandy.begin(); it != lbandy.end(); it++)
    if( (*it).min <= x && x < (*it).max ) (*it).H->Add(y,w);

if(lslix.size()>0)
  for( it = lslix.begin(); it != lslix.end(); it++)
    if( (*it).min <= y && y < (*it).max ) (*it).H->Add(x,w);

if(lsliy.size()>0)
  for( it = lsliy.begin(); it != lsliy.end(); it++)
    if( (*it).min <= x && x < (*it).max ) (*it).H->Add(y,w);

if( i<0 || i>=nx || j<0 || j>=ny ) {
  if(i<0) i=0; else if(i>=nx) i=2; else i=1;
  if(j<0) j=0; else if(j>=ny) j=2; else j=1;
  over[i][j] += w;
  over[1][1] += w;
  return;
}

data[j*nx+i] += w;
if(err2!=NULL) err2[j*nx+i] += w*w;
nHist += w;
nEntries++;
}

///////////////////////////////////////////////////////////////////
/*!
  Recherche du bin du maximum dans le pave [il,ih][jl,jh].
*/
void Histo2D::IJMax(int& imax,int& jmax,int il,int ih,int jl,int jh)
{
if( il > ih ) { il = 0; ih = nx-1; }
if( jl > jh ) { jl = 0; jh = ny-1; }
if( il < 0 ) il = 0;
if( jl < 0 ) jl = 0;
if( ih >= nx ) ih = nx-1;
if( jh >= ny ) jh = ny-1;

imax = jmax = 0;
if(nxy==1) return;

float mx=(*this)(il,jl);
for (int i=il; i<=ih; i++)
   for (int j=jl; j<=jh; j++)
      if ((*this)(i,j)>mx) {imax = i; jmax = j; mx=(*this)(i,j);}
}

/*!
  Recherche du bin du minimum dans le pave [il,ih][jl,jh].
*/
void Histo2D::IJMin(int& imax,int& jmax,int il,int ih,int jl,int jh)
{
if( il > ih ) { il = 0; ih = nx-1; }
if( jl > jh ) { jl = 0; jh = ny-1; }
if( il < 0 ) il = 0;
if( jl < 0 ) jl = 0;
if( ih >= nx ) ih = nx-1;
if( jh >= ny ) jh = ny-1;

imax = jmax = 0;
if(nxy==1) return;

float mx=(*this)(il,jl);
for (int i=il; i<=ih; i++)
   for (int j=jl; j<=jh; j++)
      if ((*this)(i,j)<mx) {imax = i; jmax = j; mx=(*this)(i,j);}
}


/*!
  Recherche du maximum dans le pave [il,ih][jl,jh].
*/
float Histo2D::VMax(int il,int ih,int jl,int jh) const
{
if( il > ih ) { il = 0; ih = nx-1; }
if( jl > jh ) { jl = 0; jh = ny-1; }
if( il < 0 ) il = 0;
if( jl < 0 ) jl = 0;
if( ih >= nx ) ih = nx-1;
if( jh >= ny ) jh = ny-1;

float mx=(*this)(il,jl);
if(nxy==1) return mx;
for (int i=il; i<=ih; i++)
   for (int j=jl; j<=jh; j++)
      if ((*this)(i,j)>mx) mx=(*this)(i,j);
return mx;
}

/*!
  Recherche du minimum dans le pave [il,ih][jl,jh].
*/
float Histo2D::VMin(int il,int ih,int jl,int jh) const
{
if( il > ih ) { il = 0; ih = nx-1; }
if( jl > jh ) { jl = 0; jh = ny-1; }
if( il < 0 ) il = 0;
if( jl < 0 ) jl = 0;
if( ih >= nx ) ih = nx-1;
if( jh >= ny ) jh = ny-1;

float mx=(*this)(il,jl);
if(nxy==1) return mx;
for (int i=il; i<=ih; i++)
   for (int j=jl; j<=jh; j++)
      if ((*this)(i,j)<mx) mx=(*this)(i,j);
return mx;
}

///////////////////////////////////////////////////////////////////
/*!
  Renvoie les under.overflow dans les 8 quadrants.
  \verbatim
     over[3][3]:        20 | 21 | 22
                           |    |
                       --------------
                           |    |
                        10 | 11 | 12        11 = all overflow+underflow
                           |    |
                       --------------
                           |    |
                        00 | 01 | 02
  \endverbatim
*/
float Histo2D::NOver(int i,int j) const
{
if( i < 0 || i>=3 || j < 0 || j>=3 ) return over[1][1];
return over[i][j];
}


///////////////////////////////////////////////////////////////////
/*!
  Retourne le nombre de bins non-nuls.
*/
int Histo2D::BinNonNul() const
{
int non=0;
for (int i=0;i<nxy;i++) if( data[i] != 0. ) non++;
return non;
}

/*!
  Retourne le nombre de bins avec erreurs non-nulles.
*/
int Histo2D::ErrNonNul() const
{
if(err2==NULL) return -1;
int non=0;
for (int i=0;i<nxy;i++) if( err2[i] != 0. ) non++;
return non;
}

///////////////////////////////////////////////////////////////////
/*!
  Idem EstimeMax(int...) mais retourne x,y.
*/
int Histo2D::EstimeMax(float& xm,float& ym,int SzPav
                      ,int il,int ih,int jl,int jh)
{
int im,jm;
IJMax(im,jm,il,ih,jl,jh);
return EstimeMax(im,jm,xm,ym,SzPav);
}

/*!
  Determine les abscisses et ordonnees du maximum donne par im,jm
  en moyennant dans un pave SzPav x SzPav autour du maximum.
  \verbatim
     Return:
       0 = si fit maximum reussi avec SzPav pixels
       1 = si fit maximum reussi avec moins que SzPav pixels
                                   dans au moins 1 direction
       2 = si fit maximum echoue et renvoit BinCenter()
      -1 = si echec: SzPav <= 0 ou im,jm hors limites
  \endverbatim
*/
int Histo2D::EstimeMax(int im,int jm,float& xm,float& ym,int SzPav)
{
xm = ym = 0;
if( SzPav <= 0 ) return -1;
if( im < 0 || im >= nx ) return -1;
if( jm < 0 || jm >= ny ) return -1;

if( SzPav%2 == 0 ) SzPav++;
SzPav = (SzPav-1)/2;

int rc = 0;
double dxm = 0, dym = 0, wx = 0;
for(int i=im-SzPav;i<=im+SzPav;i++) {
  if( i<0 || i>= nx ) {rc=1; continue;}
  for(int j=jm-SzPav;j<=jm+SzPav;j++) {
    if( j<0 || j>= ny ) {rc=1; continue;}
    float x,y;
    BinCenter(i,j,x,y);
    dxm += x * (*this)(i,j);
    dym += y * (*this)(i,j);
    wx  += (*this)(i,j);
  }
}

if( wx > 0. ) {
  xm = dxm/wx;
  ym = dym/wx;
  return rc;
} else {
  BinCenter(im,jm,xm,ym);
  return 2;
}

}

/*!
  Pour trouver le maximum de l'histogramme en tenant compte
  des fluctuations.
  \verbatim
    Methode:
    1-/ On recherche le bin maximum MAX de l'histogramme
    2-/ On considere que tous les pixels compris entre [MAX-Dz,MAX]
        peuvent etre des pixels maxima.
    3-/ On identifie le bin maximum en choissisant le pixel du 2-/
        tel que la somme des pixels dans un pave SzPav x SzPav soit maximale.
    INPUT:
     SzPav = taille du pave pour departager
     Dz    = tolerance pour identifier tous les pixels "maximum"
    OUTPUT:
     im,jm = pixel maximum trouve
    RETURN:
      <0 =  Echec
      >0 =  nombre de pixels possibles pour le maximum
  \endverbatim
*/
int Histo2D::FindMax(int& im,int& jm,int SzPav,float Dz
                    ,int il,int ih,int jl,int jh)
{
if( il > ih ) { il = 0; ih = nx-1; }
if( jl > jh ) { jl = 0; jh = ny-1; }
if( il < 0 ) il = 0;
if( jl < 0 ) jl = 0;
if( ih >= nx ) ih = nx-1;
if( jh >= ny ) jh = ny-1;
if( SzPav < 0 ) SzPav = 0;
  else { if( SzPav%2 == 0 ) SzPav++; SzPav = (SzPav-1)/2;}
if( Dz < 0 ) Dz = 0.;
float max = VMax(il,ih,jl,jh) - Dz;
int nmax = 0;
float sumx = -MAXFLOAT;
for(int i=il;i<=ih;i++) for(int j=jl;j<=jh;j++) {
  if( (*this)(i,j) < max) continue;
  nmax++;
  float sum = 0.;
  for(int ii=i-SzPav;ii<=i+SzPav;ii++) {
    if( ii<0 || ii >= nx ) continue;
    for(int jj=j-SzPav;jj<=j+SzPav;jj++) {
      if( jj<0 || jj >= ny ) continue;
      sum += (*this)(ii,jj);
    }
  }
  if( sum > sumx ) { im = i; jm = j; sumx = sum;}
}
if( nmax <= 0 ) { IJMax(im,jm,il,ih,jl,jh); return 1;}
return nmax;
}

//////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////
/*!
  Fit de l'histogramme par ``gfit''.
  \verbatim
    typ_err = 0 :
       - erreur attachee au bin si elle existe
       - sinon 1
    typ_err = 1 :
       - erreur attachee au bin si elle existe
       - sinon max( sqrt(abs(bin) ,1 )
    typ_err = 2 :
       - erreur forcee a 1
    typ_err = 3 :
       - erreur forcee a max( sqrt(abs(bin) ,1 )
    typ_err = 4 :
       - erreur forcee a 1, nulle si bin a zero.
    typ_err = 5 :
       - erreur forcee a max( sqrt(abs(bin) ,1 ),
         nulle si bin a zero.
  \endverbatim
*/
int  Histo2D::Fit(GeneralFit& gfit,unsigned short typ_err)
{
if(NBinX()*NBinY()<=0) return -1000;
if(typ_err>5) typ_err=0;

GeneralFitData mydata(2,NBinX()*NBinY());

for(int i=0;i<NBinX();i++) for(int j=0;j<NBinY();j++) {
  float x,y;
  BinCenter(i,j,x,y);
  double f = (double) (*this)(i,j);
  double saf = sqrt(fabs(f)); if(saf<1.) saf=1.;
  double e=0.;
  if(typ_err==0)      {if(HasErrors()) e=Error(i,j); else e=1.;}
  else if(typ_err==1) {if(HasErrors()) e=Error(i,j); else e=saf;}
  else if(typ_err==2) e=1.;
  else if(typ_err==3) e=saf;
  else if(typ_err==4) e=(f==0.)?0.:1.;
  else if(typ_err==5) e=(f==0.)?0.:saf;
  mydata.AddData2((double) x,(double) y,f,e);
}

gfit.SetData(&mydata);

return gfit.Fit();
}

/*!
  Retourne une classe contenant les residus du fit ``gfit''.
*/
Histo2D Histo2D::FitResidus(GeneralFit& gfit)
{
if(NBinX()<=0 || NBinY()<=0)
  throw(SzMismatchError("Histo2D::FitResidus: size mismatch\n"));
GeneralFunction* f = gfit.GetFunction();
if(f==NULL)
  throw(NullPtrError("Histo2D::FitResidus: NULL pointer\n"));
TVector<r_8> par = gfit.GetParm();
Histo2D h2(*this);
for(int i=0;i<NBinX();i++) for(int j=0;j<NBinY();j++) {
  float xc,yc;
  BinCenter(i,j,xc,yc);
  double x[2] = {(double)xc,(double)yc};
  h2(i,j) -= (float) f->Value(x,par.Data());
}
return h2;
}

/*!
  Retourne une classe contenant la fonction du fit ``gfit''.
*/
Histo2D Histo2D::FitFunction(GeneralFit& gfit)
{
if(NBinX()<=0 || NBinY()<=0)
  throw(SzMismatchError("Histo2D::FitFunction: size mismatch\n"));
GeneralFunction* f = gfit.GetFunction();
if(f==NULL)
  throw(NullPtrError("Histo2D::FitFunction: NULL pointer\n"));
TVector<r_8> par = gfit.GetParm();
Histo2D h2(*this);
for(int i=0;i<NBinX();i++) for(int j=0;j<NBinY();j++) {
  float xc,yc;
  BinCenter(i,j,xc,yc);
  double x[2] = {(double)xc,(double)yc};
  h2(i,j) = (float) f->Value(x,par.Data());
}
return h2;
}

///////////////////////////////////////////////////////////////////
/*!
  Impression des informations sur l'histogramme.
*/
void Histo2D::PrintStatus()
{
printf("~Histo::Print    nHist=%g nEntries=%d",nHist,nEntries);
if(HasErrors()) printf("  Errors=1\n"); else printf("  Errors=0\n");
printf("over: [ %g %g %g // %g %g %g  // %g %g %g ]\n"
      ,over[2][0],over[2][1],over[2][2]
      ,over[1][0],over[1][1],over[1][2]
      ,over[0][0],over[0][1],over[0][2]);
printf("  nx=%d xmin=%g xmax=%g binx=%g  ",nx,xmin,xmax,wbinx);
printf("  ny=%d ymin=%g ymax=%g biny=%g\n",ny,ymin,ymax,wbiny);
}

///////////////////////////////////////////////////////////////////
/*!
  Impression de l'histogramme sur stdout entre [il,ih] et [jl,jh].
  \verbatim
    numero d'index:  00000000001111111111222222222233333
                     01234567890123456789012345678901234
    valeur entiere:  00000000001111111111222222222233333
                     12345678901234567890123456789012345
  \endverbatim
*/
void Histo2D::Print(float min,float max
                   ,int il,int ih,int jl,int jh)
{
int ns = 35;
const char *s =    "+23456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";

if( il > ih ) { il = 0; ih = nx-1; }
if( jl > jh ) { jl = 0; jh = ny-1; }
if( il < 0 ) il = 0;
if( jl < 0 ) jl = 0;
if( ih >= nx ) ih = nx-1;
if( jh >= ny ) jh = ny-1;

PrintStatus();

if( il != 0 || ih != nx-1 || jl != 0 || jh != ny-1 ) {
  float xl,xh,yl,yh;
  BinLowEdge(il,jl,xl,yl);
  BinHighEdge(ih,jh,xh,yh);
  printf("  impression");
  printf("  en X: %d=[%d,%d] xmin=%g xmax=%g  "
        ,ih-il+1,il,ih,xl,xh);
  printf("  en Y: %d=[%d,%d] ymin=%g ymax=%g\n"
        ,jh-jl+1,jl,jh,yl,yh);
}

if(min >= max) { if(min != 0.) min = VMin(il,ih,jl,jh); else min=0.;
                 max = VMax(il,ih,jl,jh); }
if(min>max) return;
if(min==max) {min -= 1.; max += 1.;}
printf("   min=%g max=%g\n",min,max);

// imprime numero de bin en colonne
printf("\n");
if( nx-1 >= 100 ) {
  printf("     ");
  for(int i=il;i<=ih;i++) printf("%1d",(int) (i%1000)/100);
  printf("\n");
}
if( nx-1 >= 10 ) {
  printf("     ");
  for(int i=il;i<=ih;i++) printf("%1d",(int) (i%100)/10);
  printf("\n");
}
printf("     ");
for(int i=il;i<=ih;i++) printf("%1d",i%10);
printf("\n");
printf("     "); {for(int i=il;i<=ih;i++) printf("-"); printf("\n");}

// imprime histogramme
for(int j=jh;j>=jl;j--) {
  printf("%3d: ",j);
  for(int i=il;i<=ih;i++) {
    int h;
    if( 1<=max-min && max-min<=35 ) h = (int)( (*this)(i,j) - min ) - 1;
      else h = (int)( ((*this)(i,j)-min)/(max-min) * ns ) - 1;
    char c;
    if(h<0 && (*this)(i,j)>min) c = '.';
      else if(h<0) c = ' ';
        else if(h>=ns) c = '*';
          else c = s[h];
    printf("%c",c);
  }
  printf("\n");
}

// imprime numero de bin en colonne
printf("     "); {for(int i=il;i<=ih;i++) printf("-"); printf("\n");}
if( nx-1 >= 100 ) {
  printf("     ");
  for(int i=il;i<=ih;i++) printf("%1d",(int) (i%1000)/100);
  printf("\n");
}
if( nx-1 >= 10 ) {
  printf("     ");
  for(int i=il;i<=ih;i++) printf("%1d",(int) (i%100)/10);
  printf("\n");
}
printf("     ");
{for(int i=il;i<=ih;i++) printf("%1d",i%10);}
printf("\n");

}

///////////////////////////////////////////////////////////////////
// Titre        Methodes pour gerer les projections

/*!
  Pour creer la projection X.
*/
void Histo2D::SetProjX()
{
if( hprojx != NULL ) DelProjX();
hprojx = new Histo(xmin,xmax,nx);
if( err2 != NULL && hprojx != NULL ) hprojx->Errors();
}

/*!
  Pour creer la projection Y.
*/
void Histo2D::SetProjY()
{
if( hprojy != NULL ) DelProjY();
hprojy = new Histo(ymin,ymax,ny);
if( err2 != NULL && hprojy != NULL  ) hprojy->Errors();
}

/*!
  Pour creer les projections X et Y.
*/
void Histo2D::SetProj()
{
SetProjX();
SetProjY();
}

/*!
  Informations sur les projections.
*/
void Histo2D::ShowProj()
{
if( hprojx != NULL ) cout << ">>>> Projection X set : "<< hprojx <<endl;
  else cout << ">>>> NO Projection X set"<<endl;
if( hprojy != NULL ) cout << ">>>> Projection Y set : "<< hprojy <<endl;
  else cout << ">>>> NO Projection Y set"<<endl;
}

/*!
  Destruction de l'histogramme de la projection selon X.
*/
void Histo2D::DelProjX()
{
if( hprojx == NULL ) return;
delete hprojx;
hprojx = NULL;
}

/*!
  Destruction de l'histogramme de la projection selon X.
*/
void Histo2D::DelProjY()
{
if( hprojy == NULL ) return;
delete hprojy;
hprojy = NULL;
}

/*!
  Destruction des histogrammes des projections selon X et Y.
*/
void Histo2D::DelProj()
{
DelProjX();
DelProjY();
}

/*!
  Remise a zero de la projection selon X.
*/
void Histo2D::ZeroProjX()
{
if( hprojx == NULL ) return;
hprojx->Zero();
}

/*!
  Remise a zero de la projection selon Y.
*/
void Histo2D::ZeroProjY()
{
if( hprojy == NULL ) return;
hprojy->Zero();
}

/*!
  Remise a zero des projections selon X et Y.
*/
void Histo2D::ZeroProj()
{
ZeroProjX();
ZeroProjY();
}

///////////////////////////////////////////////////////////////////
// Titre        Methodes pour gerer les bandes

/*!
  Pour creer une bande en X entre ybmin et ybmax.
*/
int Histo2D::SetBandX(float ybmin,float ybmax)
{
b_s.num = lbandx.size();
b_s.min = ybmin;
b_s.max = ybmax;
b_s.H = new Histo(xmin,xmax,nx);
lbandx.push_back(b_s);
b_s.H = NULL;
return lbandx.size()-1;
}

/*!
  Pour creer une bande en Y entre xbmin et xbmax.
*/
int Histo2D::SetBandY(float xbmin,float xbmax)
{
b_s.num = lbandy.size();
b_s.min = xbmin;
b_s.max = xbmax;
b_s.H = new Histo(ymin,ymax,ny);
lbandy.push_back(b_s);
b_s.H = NULL;
return lbandy.size()-1;
}

/*!
  Destruction des histogrammes des bandes selon X.
*/
void Histo2D::DelBandX()
{
if( lbandx.size() <= 0 ) return;
for(list<bande_slice>::iterator i = lbandx.begin(); i != lbandx.end(); i++)
       if( (*i).H != NULL ) {delete (*i).H; (*i).H=NULL;}
lbandx.erase(lbandx.begin(),lbandx.end());
}

/*!
  Destruction des histogrammes des bandes selon Y.
*/
void Histo2D::DelBandY()
{
if( lbandy.size() <= 0 ) return;
for(list<bande_slice>::iterator i = lbandy.begin(); i != lbandy.end(); i++)
       if( (*i).H != NULL ) {delete (*i).H; (*i).H=NULL;}
lbandy.erase(lbandy.begin(),lbandy.end());
}

/*!
  Remise a zero des bandes selon X.
*/
void Histo2D::ZeroBandX()
{
if( lbandx.size() <= 0 ) return;
list<bande_slice>::iterator i;
for(i = lbandx.begin(); i != lbandx.end(); i++)
              (*i).H->Zero();
}

/*!
  Remise a zero des bandes selon Y.
*/
void Histo2D::ZeroBandY()
{
if( lbandy.size() <= 0 ) return;
list<bande_slice>::iterator i;
for(i = lbandy.begin(); i != lbandy.end(); i++)
              (*i).H->Zero();
}

/*!
  Retourne un pointeur sur la bande numero `n' selon X.
*/
Histo* Histo2D::HBandX(int n) const
{
if( lbandx.size() <= 0 || n < 0 || n >= (int) lbandx.size() ) return NULL;
for(list<bande_slice>::const_iterator i = lbandx.begin(); i != lbandx.end(); i++)
              if( (*i).num == n ) return (*i).H;
return NULL;
}

/*!
  Retourne un pointeur sur la bande numero `n' selon Y.
*/
Histo* Histo2D::HBandY(int n) const
{
if( lbandy.size() <= 0 || n < 0 || n >= (int) lbandy.size() ) return NULL;
for(list<bande_slice>::const_iterator i = lbandy.begin(); i != lbandy.end(); i++)
              if( (*i).num == n ) return (*i).H;
return NULL;
}

/*!
  Retourne les limites de la bande numero `n' selon X.
*/
void Histo2D::GetBandX(int n,float& ybmin,float& ybmax) const
{
ybmin = 0.; ybmax = 0.;
if( lbandx.size() <= 0 || n < 0 || n >= (int) lbandx.size() ) return;
for(list<bande_slice>::const_iterator i = lbandx.begin(); i != lbandx.end(); i++)
  if( (*i).num == n ) { ybmin = (*i).min; ybmax = (*i).max; return;}
return;
}

/*!
  Retourne les limites de la bande numero `n' selon Y.
*/
void Histo2D::GetBandY(int n,float& xbmin,float& xbmax) const
{
xbmin = 0.; xbmax = 0.;
if( lbandy.size() <= 0 || n < 0 || n >= (int) lbandy.size() ) return;
for(list<bande_slice>::const_iterator i = lbandy.begin(); i != lbandy.end(); i++)
  if( (*i).num == n ) { xbmin = (*i).min; xbmax = (*i).max; return;}
return;
}

/*!
  Informations sur les bandes.
*/
void Histo2D::ShowBand(int lp)
{
  cout << ">>>> Nombre de bande X : " << lbandx.size() << endl;
if( lp>0 && lbandx.size()>0 ) {
  list<bande_slice>::iterator i;
  for(i = lbandx.begin(); i != lbandx.end(); i++) {
    cout<<"  "<<(*i).num<<" de ymin="<<(*i).min<<" a ymax="<<(*i).max;
    if(lp>1) cout << "   H=" << (*i).H;
    cout << endl;
  }
}

cout << ">>>> Nombre de bande Y : " << lbandy.size() << endl;
if( lp>0 && lbandy.size()>0 ) {
  list<bande_slice>::iterator i;
  for(i = lbandy.begin(); i != lbandy.end(); i++) {
    cout<<"  "<<(*i).num<<" de xmin="<<(*i).min<<" a xmax="<<(*i).max;
    if(lp>1) cout << "   H=" << (*i).H;
    cout << endl;
  }
}
}

///////////////////////////////////////////////////////////////////
// Titre  Methodes pour gerer les bandes equidistantes ou slices

/*!
  Pour creer `nsli' bandes equidistantes selon X.
*/
int Histo2D::SetSliX(int nsli)
{
if( nsli <= 0 ) return -1;
if( nsli >  ny ) nsli = ny;
if( lslix.size() > 0 ) DelSliX();
float w = (ymax-ymin)/nsli;

for(int i=0; i<nsli; i++ ) {
  b_s.num = i;
  b_s.min = ymin + i*w;
  b_s.max = b_s.min + w;
  b_s.H = new Histo(xmin,xmax,nx);
  lslix.push_back(b_s);
  b_s.H = NULL;
}
return (int) lslix.size();
}

/*!
  Pour creer `nsli' bandes equidistantes selon Y.
*/
int Histo2D::SetSliY(int nsli)
{
if( nsli <= 0 ) return -1;
if( nsli >  nx ) nsli = nx;
if( lsliy.size() > 0 ) DelSliY();
float w = (xmax-xmin)/nsli;

for(int i=0; i<nsli; i++ ) {
  b_s.num = i;
  b_s.min = xmin + i*w;
  b_s.max = b_s.min + w;
  b_s.H = new Histo(ymin,ymax,ny);
  lsliy.push_back(b_s);
  b_s.H = NULL;
}
return (int) lsliy.size();
}

/*!
  Destruction des bandes equidistantes selon X.
*/
void Histo2D::DelSliX()
{
if( lslix.size() <= 0 ) return;
for(list<bande_slice>::iterator i = lslix.begin(); i != lslix.end(); i++)
       if( (*i).H != NULL ) {delete (*i).H; (*i).H=NULL;}
lslix.erase(lslix.begin(),lslix.end());
}

/*!
  Destruction des bandes equidistantes selon Y.
*/
void Histo2D::DelSliY()
{
if( lsliy.size() <= 0 ) return;
for(list<bande_slice>::iterator i = lsliy.begin(); i != lsliy.end(); i++)
       if( (*i).H != NULL ) {delete (*i).H; (*i).H=NULL;}
lsliy.erase(lsliy.begin(),lsliy.end());
}

/*!
  Remise a zero des bandes equidistantes selon X.
*/
void Histo2D::ZeroSliX()
{
if( lslix.size() <= 0 ) return;
list<bande_slice>::iterator i;
for(i = lslix.begin(); i != lslix.end(); i++)
              (*i).H->Zero();
}

/*!
  Remise a zero des bandes equidistantes selon Y.
*/
void Histo2D::ZeroSliY()
{
if( lsliy.size() <= 0 ) return;
list<bande_slice>::iterator i;
for(i = lsliy.begin(); i != lsliy.end(); i++)
              (*i).H->Zero();
}

/*!
  Retourne un pointeur sur la bande equidistante numero `n'
  selon X.
*/
Histo* Histo2D::HSliX(int n) const
{
if( lslix.size() <= 0 || n < 0 || n >= (int) lslix.size() ) return NULL;
for(list<bande_slice>::const_iterator i = lslix.begin(); i != lslix.end(); i++)
              if( (*i).num == n ) return (*i).H;
return NULL;
}

/*!
  Retourne un pointeur sur la bande equidistante numero `n'
  selon Y.
*/
Histo* Histo2D::HSliY(int n) const
{
if( lsliy.size() <= 0 || n < 0 || n >= (int) lsliy.size() ) return NULL;
for(list<bande_slice>::const_iterator i = lsliy.begin(); i != lsliy.end(); i++)
              if( (*i).num == n ) return (*i).H;
return NULL;
}

/*!
  Informations sur les bandes equidistantes.
*/
void Histo2D::ShowSli(int lp)
{
list<bande_slice>::iterator i;
cout << ">>>> Nombre de slice X : " << lslix.size() << endl;
if( lp>0 && lslix.size() > 0 )
  for(i = lslix.begin(); i != lslix.end(); i++) {
    cout<<"  "<<(*i).num<<" de ymin="<<(*i).min<<" a ymax="<<(*i).max;
    if(lp>1) cout << "   H=" << (*i).H;
    cout << endl;
  }

cout << ">>>> Nombre de slice Y : " << lsliy.size() << endl;
if( lp>0 && lsliy.size()>0 )
  for(i = lsliy.begin(); i != lsliy.end(); i++) {
    cout<<"  "<<(*i).num<<" de xmin="<<(*i).min<<" a xmax="<<(*i).max;
    if(lp>1) cout << "   H=" << (*i).H;
    cout << endl;
  }
}

///////////////////////////////////////////////////////////
// --------------------------------------------------------
//   Les objets delegues pour la gestion de persistance 
// --------------------------------------------------------
///////////////////////////////////////////////////////////


void  ObjFileIO<Histo2D>::ReadSelf(PInPersist& is)
{
char strg[256];

if(dobj==NULL) dobj=new Histo2D;
  else         dobj->Delete();

float min,max;
int_4 errok, projx, projy, nslix, nsliy, nbanx, nbany;

// Lecture entete
is.GetLine(strg, 255);
is.GetLine(strg, 255);
is.GetLine(strg, 255);
is.GetLine(strg, 255);
is.GetLine(strg, 255);
is.GetLine(strg, 255);

// Lecture variables de definitions
is.Get(dobj->nx);
is.Get(dobj->ny);
is.Get(dobj->nxy);
is.Get(errok);
is.Get(dobj->nEntries);
is.Get(dobj->nHist);

is.Get(dobj->xmin);
is.Get(dobj->xmax);
is.Get(dobj->ymin);
is.Get(dobj->ymax);
is.Get(dobj->wbinx);
is.Get(dobj->wbiny);

is.Get(&(dobj->over[0][0]),9);

is.Get(projx);
is.Get(projy);
is.Get(nslix);
is.Get(nsliy);
is.Get(nbanx);
is.Get(nbany);

// Lecture histo2D
dobj->data = new float[dobj->nxy];
is.GetLine(strg, 255);
{for(int j=0;j<dobj->ny;j++) is.Get(dobj->data+j*dobj->nx,dobj->nx);}

// Lecture erreurs
if(errok) {
  is.GetLine(strg, 255);
  dobj->err2 = new double[dobj->nxy];
  for(int j=0;j<dobj->ny;j++) is.Get(dobj->err2+j*dobj->nx,dobj->nx);
}

// Lecture des projections
if(projx) {
  is.GetLine(strg, 255);
  dobj->SetProjX();
  ObjFileIO<Histo> fio_h(dobj->hprojx);
  fio_h.Read(is);
}
if(projy) {
  is.GetLine(strg, 255);
  dobj->SetProjY();
  ObjFileIO<Histo> fio_h(dobj->hprojy);
  fio_h.Read(is);
}

// Lecture des slices
if(nslix>0) {
  is.GetLine(strg, 255);
  dobj->SetSliX(nslix);
  ASSERT (nslix==dobj->NSliX());
  for(int j=0;j<dobj->NSliX();j++)
    {ObjFileIO<Histo> fio_h(dobj->HSliX(j)); fio_h.Read(is);}
}
if(nsliy>0) {
  is.GetLine(strg, 255);
  dobj->SetSliY(nsliy);
  ASSERT (nsliy==dobj->NSliY());
  for(int j=0;j<dobj->NSliY();j++)
    {ObjFileIO<Histo> fio_h(dobj->HSliY(j)); fio_h.Read(is);}
}

// Lecture des bandes
if( nbanx>0 ) {
  is.GetLine(strg, 255);
  {for(int j=0; j<nbanx; j++) {
    is.Get(min); is.Get(max);
    dobj->SetBandX(min,max);
  }}
  ASSERT (nbanx==dobj->NBandX());
  {for(int j=0; j<dobj->NBandX(); j++) {
    ObjFileIO<Histo> fio_h(dobj->HBandX(j));
    fio_h.Read(is);
  }}
}
if( nbany>0 ) {
  is.GetLine(strg, 255);
  {for(int j=0; j<nbany; j++) {
    is.Get(min); is.Get(max);
    dobj->SetBandY(min,max);
  }}
  ASSERT (nbany==dobj->NBandY());
  {for(int j=0; j<dobj->NBandY(); j++) {
    ObjFileIO<Histo> fio_h(dobj->HBandY(j));
    fio_h.Read(is);
  }}
}

return;
}

void ObjFileIO<Histo2D>::WriteSelf(POutPersist& os) const
{
if (dobj == NULL)   return;
char strg[256];

// Que faut-il ecrire?
int_4 errok = (dobj->err2) ? 1 : 0;
int_4 projx = (dobj->hprojx) ? 1 : 0;
int_4 projy = (dobj->hprojy) ? 1 : 0;
int_4 nslix = dobj->NSliX();
int_4 nsliy = dobj->NSliY();
int_4 nbanx = dobj->NBandX();
int_4 nbany = dobj->NBandY();

// Ecriture entete pour identifier facilement
sprintf(strg,"nx=%d  ny=%d  nxy=%d errok=%1d"
       ,dobj->nx,dobj->ny,dobj->nxy,errok);
os.PutLine(strg);
sprintf(strg,"nHist=%g nEntries=%d",dobj->nHist,dobj->nEntries);
os.PutLine(strg);
sprintf(strg,"wbinx=%g wbiny=%g",dobj->wbinx,dobj->wbiny);
os.PutLine(strg);
sprintf(strg,"xmin=%g xmax=%g ymin=%g ymax=%g"
       ,dobj->xmin,dobj->xmax,dobj->ymin,dobj->ymax);
os.PutLine(strg);
sprintf(strg,"projx/y=%d %d nbandx/y=%d %d nbslix/y=%d %d"
       ,projx,projy,nbanx,nbany,nslix,nsliy);
os.PutLine(strg);
sprintf(strg,"over %g %g %g %g %g %g %g %g %g"
       ,dobj->over[0][0],dobj->over[0][1],dobj->over[0][2]
       ,dobj->over[1][0],dobj->over[1][1],dobj->over[1][2]
       ,dobj->over[2][0],dobj->over[2][1],dobj->over[2][2]);
os.PutLine(strg);

// Ecriture variables de definitions
os.Put(dobj->nx);
os.Put(dobj->ny);
os.Put(dobj->nxy);
os.Put(errok);
os.Put(dobj->nEntries);
os.Put(dobj->nHist);

os.Put(dobj->xmin);
os.Put(dobj->xmax);
os.Put(dobj->ymin);
os.Put(dobj->ymax);
os.Put(dobj->wbinx);
os.Put(dobj->wbiny);

os.Put(&(dobj->over[0][0]),9);

os.Put(projx);
os.Put(projy);
os.Put(nslix);
os.Put(nsliy);
os.Put(nbanx);
os.Put(nbany);

// Ecriture histo2D
sprintf(strg,"Histo2D: Tableau des donnees %d = %d * %d"
       ,dobj->nxy,dobj->nx,dobj->ny);
os.PutLine(strg);
{for(int j=0;j<dobj->ny;j++) os.Put(dobj->data+j*dobj->nx,dobj->nx);}

// Ecriture erreurs
if(errok) {
  sprintf(strg,"Histo2D: Tableau des erreurs %d = %d * %d"
         ,dobj->nxy,dobj->nx,dobj->ny);
  os.PutLine(strg);
  for(int j=0;j<dobj->ny;j++) os.Put(dobj->err2+j*dobj->nx,dobj->nx);
}

// Ecriture des projections
if(projx) {
  sprintf(strg,"Histo2D: Projection X");
  os.PutLine(strg);
  ObjFileIO<Histo> fio_h(dobj->hprojx); fio_h.Write(os);
}
if(projy) {
  sprintf(strg,"Histo2D: Projection Y");
  os.PutLine(strg);
  ObjFileIO<Histo> fio_h(dobj->hprojy); fio_h.Write(os);
}

// Ecriture des slices
if(nslix>0) {
  sprintf(strg,"Histo2D: Slices X %d",nslix);
  os.PutLine(strg);
  for(int j=0;j<nslix;j++) {
    Histo* h = dobj->HSliX(j);
    ObjFileIO<Histo> fio_h(h); fio_h.Write(os);
  }
}
if(nsliy>0) {
  sprintf(strg,"Histo2D: Slices Y %d",nsliy);
  os.PutLine(strg);
  for(int j=0;j<nsliy;j++) {
    Histo* h = dobj->HSliY(j);
    ObjFileIO<Histo> fio_h(h); fio_h.Write(os);
  }
}

// Ecriture des bandes
if( nbanx>0 ) {
  sprintf(strg,"Histo2D: Bandes X %d",nbanx);
  os.PutLine(strg);
  list<Histo2D::bande_slice>::const_iterator it;
  for(it = dobj->lbandx.begin(); it != dobj->lbandx.end(); it++) {
    float min = (*it).min; float max = (*it).max;
    os.Put(min); os.Put(max);
  }
  for(it = dobj->lbandx.begin(); it != dobj->lbandx.end(); it++) {
    Histo* h = (*it).H;
    ObjFileIO<Histo> fio_h(h); fio_h.Write(os);
  }
}
if( nbany>0 ) {
  sprintf(strg,"Histo2D: Bandes Y %d",nbany);
  os.PutLine(strg);
  list<Histo2D::bande_slice>::const_iterator it;
  for(it = dobj->lbandy.begin(); it != dobj->lbandy.end(); it++) {
    float min = (*it).min; float max = (*it).max;
    os.Put(min); os.Put(max);
  }
  for(it = dobj->lbandy.begin(); it != dobj->lbandy.end(); it++) {
    Histo* h = (*it).H;
    ObjFileIO<Histo> fio_h(h); fio_h.Write(os);
  }
}

return;
}


#ifdef __CXX_PRAGMA_TEMPLATES__
#pragma define_template ObjFileIO<Histo2D>
#endif

#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
template class ObjFileIO<Histo2D>;
#endif