source: Sophya/trunk/SophyaLib/NTools/histos2.cc@ 243

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

#include "defs.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        Histo2D
// Lib  Outils++ 
// include      histos2.h
//
//      Classe d'histogrammes 2D
//--

///////////////////////////////////////////////////////////////////
//++
// Titre        Constructeurs
//--

//++
Histo2D::Histo2D(float xMin, float xMax, int nxBin
                ,float yMin, float yMax, int nyBin)
//
//      Createur d'un histogramme 2D ayant nxBin,nyBin bins
//      entre xMin,xMax et yMin,yMax.
//--
      : 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)
{
DBASSERT(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
}

//++
Histo2D::Histo2D(const Histo2D& h)
//
//      Constructeur par copie.
//--
{
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
}

//++
Histo2D::Histo2D()
//
//      Constructeur par defaut.
//--
      : 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
}

//++
Histo2D::Histo2D(char *flnm)
//
//      Constructeur par lecture d'un fichier ppersist.
//--
      : 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)
{
b_s.H = NULL;
PInPersist s(flnm);
Read(s);
END_CONSTRUCTOR
}

///////////////////////////////////////////////////////////////////
//++
void Histo2D::Delete()
//
//      Desallocation de la place de l'histogramme (fct privee).
//--
{
  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;
}

//++
Histo2D::~Histo2D()
//
//      Destructeur.
//--
{
Delete();
}

///////////////////////////////////////////////////////////////////
//++
// Titre        Methodes generales
//--

//++
void Histo2D::Zero()
//
//      Remise a zero du contenu, des erreurs et des valeurs.
//--
{
  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();
}

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

///////////////////////////////////////////////////////////////////
//++
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;
}

///////////////////////////////////////////////////////////////////
//++
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;
}

//++
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;
}

//++
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;
}

//++
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;
}

///////////////////////////////////////////////////////////////////
//++
Histo2D operator * (const Histo2D& a, double b)
//
//--
{
  Histo2D result(a);
  return (result *= b);
}

//++
Histo2D operator * (double b, const Histo2D& a)
//
//--
{
  Histo2D result(a);
  return (result *= b);
}

//++
Histo2D operator / (const Histo2D& a, double b)
//
//--
{
  Histo2D result(a);
  return (result /= b);
}

//++
Histo2D operator + (const Histo2D& a, double b)
//
//--
{
  Histo2D result(a);
  return (result += b);
}

//++
Histo2D operator + (double b, const Histo2D& a)
//
//--
{
  Histo2D result(a);
  return (result += b);
}

//++
Histo2D operator - (const Histo2D& a, double b)
//
//--
{
  Histo2D result(a);
  return (result -= b);
}

//++
Histo2D operator - (double b, const Histo2D& a)
//
//--
{
  Histo2D result(a);
  result *= -1.;
  return (result += b);
}

///////////////////////////////////////////////////////////////////
//++
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;
}

//++
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;
}

//++
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;
}

//++
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;
}

///////////////////////////////////////////////////////////////////
//++
Histo2D operator + (const Histo2D& a, const Histo2D& b)
//
//--
{
if (b.nx!=a.nx || b.ny!=a.ny) THROW(sizeMismatchErr);
Histo2D c(a);
return (c += b);
}

//++
Histo2D operator - (const Histo2D& a, const Histo2D& b)
//
//--
{
if (b.nx!=a.nx || b.ny!=a.ny) THROW(sizeMismatchErr);
Histo2D c(a);
return (c -= b);
}

//++
Histo2D operator * (const Histo2D& a, const Histo2D& b)
//
//--
{
if (b.nx!=a.nx || b.ny!=a.ny) THROW(sizeMismatchErr);
Histo2D c(a);
return (c *= b);
}

//++
Histo2D operator / (const Histo2D& a, const Histo2D& b)
//
//--
{
if (b.nx!=a.nx || b.ny!=a.ny) THROW(sizeMismatchErr);
Histo2D c(a);
return (c /= b);
}

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

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

//++
//| 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
//--

//++
void Histo2D::GetValue(Matrix &v)
//
//      Remplissage d'un tableau avec les valeurs du contenu.
//--
{
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;
}

//++
void Histo2D::GetError2(Matrix &v)
//
//      Remplissage d'un tableau avec les valeurs du carre des erreurs.
//--
{
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;
}

//++
void Histo2D::GetError(Matrix &v)
//
//      Remplissage d'un tableau avec les valeurs des erreurs.
//--
{
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;
}

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

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

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

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

//++
void Histo2D::PutError(Matrix &v)
//
//      Remplissage des erreurs de l'histo avec les valeurs d'un tableau.
//--
{
int i,j;
if(v.NRows()!=nx || v.NCol()!=ny) THROW(sizeMismatchErr);
if(!err2) Errors();
for(i=0;i<nx;i++) for(j=0;j<ny;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 *********/
//++
void Histo2D::Add(float x, float y, float w)
//
//      Addition du contenu de l'histo pour x,y poids 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++;
}

///////////////////////////////////////////////////////////////////
//++
void Histo2D::IJMax(int& imax,int& jmax,int il,int ih,int jl,int jh)
//
//      Recherche du bin du maximum dans le pave [il,ih][jl,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);}
}

//++
void Histo2D::IJMin(int& imax,int& jmax,int il,int ih,int jl,int jh)
//
//      Recherche du bin du minimum dans le pave [il,ih][jl,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);}
}


//++
float Histo2D::VMax(int il,int ih,int jl,int jh) const
//
//      Recherche du maximum dans le pave [il,ih][jl,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;

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;
}

//++
float Histo2D::VMin(int il,int ih,int jl,int jh) const
//
//      Recherche du minimum dans le pave [il,ih][jl,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;

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;
}

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


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

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

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

//++
int Histo2D::EstimeMax(int im,int jm,float& xm,float& ym,int SzPav)
//
//      Determine les abscisses et ordonnees du maximum donne par im,jm
//      en moyennant dans un pave SzPav x SzPav autour du maximum.
//|  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
//--
{
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;
}

}

//++
int Histo2D::FindMax(int& im,int& jm,int SzPav,float Dz
                    ,int il,int ih,int jl,int jh)
//
//      Pour trouver le maximum de l'histogramme en tenant compte
//      des fluctuations.
//| 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
//--
{
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;
}

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

//////////////////////////////////////////////////////////
//++
int  Histo2D::Fit(GeneralFit& gfit,unsigned short typ_err)
//
//      Fit de l'histogramme par ``gfit''.
//| 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.
//--
{
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;
  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();
}

//++
Histo2D* Histo2D::FitResidus(GeneralFit& gfit)
//
//      Retourne une classe contenant les residus du fit ``gfit''.
//--
{
if(NBinX()<=0 || NBinY()<=0) return NULL;
GeneralFunction* f = gfit.GetFunction();
if(f==NULL) return NULL;
Vector par = gfit.GetParm();
Histo2D* h2 = new Histo2D(*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;
}

//++
Histo2D* Histo2D::FitFunction(GeneralFit& gfit)
//
//      Retourne une classe contenant la fonction du fit ``gfit''.
//--
{
if(NBinX()<=0 || NBinY()<=0) return NULL;
GeneralFunction* f = gfit.GetFunction();
if(f==NULL) return NULL;
Vector par = gfit.GetParm();
Histo2D* h2 = new Histo2D(*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;
}

///////////////////////////////////////////////////////////////////
//++
void Histo2D::PrintStatus()
//
//      Impression des informations sur l'histogramme.
//--
{
printf("~Histo::Print    nHist=%g nEntries=%d\n",nHist,nEntries);
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);
}

///////////////////////////////////////////////////////////////////
//++
void Histo2D::Print(float min,float max
                   ,int il,int ih,int jl,int jh)
//
//      Impression de l'histogramme sur stdout entre [il,ih] et [jl,jh].
//--
{
int ns = 35;
// numero d'index:  00000000001111111111222222222233333
//                  01234567890123456789012345678901234
// valeur entiere:  00000000001111111111222222222233333
//                  12345678901234567890123456789012345
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");

}

// Rappel des inline functions pour commentaires
//++
// inline float   XMin()
//      Retourne l'abscisse minimum.
//--
//++
// inline float   XMax()
//      Retourne l'abscisse maximum.
//--
//++
// inline float   YMin()
//      Retourne l'ordonnee minimum.
//--
//++
// inline float   YMax()
//      Retourne l'ordonnee maximum.
//--
//++
// inline int     NBinX()
//      Retourne le nombre de bins selon X.
//--
//++
// inline int     NBinY()
//      Retourne le nombre de bins selon Y.
//--
//++
// inline float   WBinX()
//      Retourne la largeur du bin selon X.
//--
//++
// inline float   WBinY()
//      Retourne la largeur du bin selon Y.
//--
//++
// inline float*  Bins() const
//      Retourne le pointeur sur le tableaux des contenus.
//--
//++
// inline float   operator()(int i,int j) const
//      Retourne le contenu du bin i,j.
//--
//++
// inline float&  operator()(int i,int j)
//      Remplit le contenu du bin i,j.
//--
//++
// inline float   Error(int i,int j)  const
//      Retourne l'erreur du bin i,j.
//--
//++
// inline double  Error2(int i,int j) const
//      Retourne l'erreur au carre du bin i,j.
//--
//++
// inline double& Error2(int i,int j)
//      Remplit l'erreur au carre du bin i,j.
//--
//++
// inline float   NData() const
//      Retourne la somme ponderee.
//--
//++
// inline int     NEntries() const
//      Retourne le nombre d'entrees.
//--
//++
// inline void BinLowEdge(int i,int j,float& x,float& y)
//      Retourne l'abscisse et l'ordonnee du coin inferieur du bin i,j.
//--
//++
// inline void BinCenter(int i,int j,float& x,float& y)
//      Retourne l'abscisse et l'ordonnee du centre du bin i,j.
//--
//++
// inline void BinHighEdge(int i,int j,float& x,float& y)
//      Retourne l'abscisse et l'ordonnee du coin superieur du bin i,j.
//--
//++
// inline void FindBin(float x,float y,int& i,int& j)
//      Retourne les numeros du bin contenant l'abscisse et l'ordonnee x,y.
//--

///////////////////////////////////////////////////////////////////
//++
// Titre        Methodes pour gerer les projections
//--

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

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

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

//++
void Histo2D::ShowProj()
//
//      Informations sur les projections.
//--
{
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;
}

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

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

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

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

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

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

// Rappel des inline functions pour commentaires
//++
// inline Histo*  HProjX()
//      Retourne le pointeur sur l'histo 1D de la projection selon X.
//--
//++
// inline Histo*  HProjY()
//      Retourne le pointeur sur l'histo 1D de la projection selon Y.
//--

///////////////////////////////////////////////////////////////////
//++
// Titre        Methodes pour gerer les bandes
//--

//++
int Histo2D::SetBandX(float ybmin,float ybmax)
//
//      Pour creer une bande en X entre ybmin et 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;
}

//++
int Histo2D::SetBandY(float xbmin,float xbmax)
//
//      Pour creer une bande en Y entre xbmin et 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;
}

//++
void Histo2D::DelBandX()
//
//      Destruction des histogrammes des bandes selon X.
//--
{
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());
}

//++
void Histo2D::DelBandY()
//
//      Destruction des histogrammes des bandes selon Y.
//--
{
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());
}

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

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

//++
Histo* Histo2D::HBandX(int n) const
//
//      Retourne un pointeur sur la bande numero `n' selon X.
//--
{
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;
}

//++
Histo* Histo2D::HBandY(int n) const
//
//      Retourne un pointeur sur la bande numero `n' selon Y.
//--
{
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;
}

//++
void Histo2D::GetBandX(int n,float& ybmin,float& ybmax) const
//
//      Retourne les limites de la bande numero `n' selon X.
//--
{
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;
}

//++
void Histo2D::GetBandY(int n,float& xbmin,float& xbmax) const
//
//      Retourne les limites de la bande numero `n' selon Y.
//--
{
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;
}

//++
void Histo2D::ShowBand(int lp)
//
//      Informations sur les bandes.
//--
{
  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;
  }
}
}

// Rappel des inline functions pour commentaires
//++
// inline int     NBandX()
//      Retourne le nombre de bandes selon X
//--
//++
// inline int     NBandY()
//      Retourne le nombre de bandes selon Y
//--

///////////////////////////////////////////////////////////////////
//++
// Titre        Methodes pour gerer les bandes equidistantes
//--

//++
int Histo2D::SetSliX(int nsli)
//
//      Pour creer `nsli' bandes equidistantes selon X.
//--
{
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();
}

//++
int Histo2D::SetSliY(int nsli)
//
//      Pour creer `nsli' bandes equidistantes selon Y.
//--
{
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();
}

//++
void Histo2D::DelSliX()
//
//      Destruction des bandes equidistantes selon X.
//--
{
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());
}

//++
void Histo2D::DelSliY()
//
//      Destruction des bandes equidistantes selon Y.
//--
{
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());
}

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

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

//++
Histo* Histo2D::HSliX(int n) const
//
//      Retourne un pointeur sur la bande equidistante numero `n'
//      selon X.
//--
{
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;
}

//++
Histo* Histo2D::HSliY(int n) const
//
//      Retourne un pointeur sur la bande equidistante numero `n'
//      selon Y.
//--
{
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;
}

//++
void Histo2D::ShowSli(int lp)
//
//      Informations sur les bandes equidistantes.
//--
{
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;
  }
}

// Rappel des inline functions pour commentaires
//++
// inline int     NSliX()
//      Retourne le nombre de slices selon X
//--
//++
// inline int     NSliY()
//      Retourne le nombre de slices selon Y
//--


///////////////////////////////////////////////////////////////////
//++
// Titre        Methodes pour ecriture ppersist
//--

//++
void  Histo2D::WriteSelf(POutPersist& s)  const
//
//      Ecriture fichier de type ppersist.
//--
{
char strg[256];
int j;
float min,max;
Histo* h;

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

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

// Ecriture variables de definitions
s.PutI4(nx);
s.PutI4(ny);
s.PutI4(nxy);
s.PutI4(errok);
s.PutI4(nEntries);
s.PutR8(nHist);

s.PutR4(xmin);
s.PutR4(xmax);
s.PutR4(ymin);
s.PutR4(ymax);
s.PutR4(wbinx);
s.PutR4(wbiny);

s.PutR4s(&over[0][0],9);

s.PutI4(projx);
s.PutI4(projy);
s.PutI4(nslix);
s.PutI4(nsliy);
s.PutI4(nbanx);
s.PutI4(nbany);

// Ecriture histo2D
sprintf(strg,"Histo2D: Tableau des donnees %d = %d * %d",nxy,nx,ny);
s.PutLine(strg);
for(j=0;j<ny;j++) s.PutR4s(data+j*nx,nx);


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

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

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

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

return;
}

//++
void  Histo2D::ReadSelf(PInPersist& s)
//
//      Lecture fichier de type ppersist.
//--
{
Delete();

int j;
float min,max;
Histo* h;
char strg[256];
int_4 errok, projx, projy, nslix, nsliy, nbanx, nbany;

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


// Lecture variables de definitions
s.GetI4(nx);
s.GetI4(ny);
s.GetI4(nxy);
s.GetI4(errok);
s.GetI4(nEntries);
s.GetR8(nHist);

s.GetR4(xmin);
s.GetR4(xmax);
s.GetR4(ymin);
s.GetR4(ymax);
s.GetR4(wbinx);
s.GetR4(wbiny);

s.GetR4s(&over[0][0],9);

s.GetI4(projx);
s.GetI4(projy);
s.GetI4(nslix);
s.GetI4(nsliy);
s.GetI4(nbanx);
s.GetI4(nbany);

// Lecture histo2D
data = new float[nxy];
s.GetLine(strg, 255);
for(j=0;j<ny;j++) s.GetR4s(data+j*nx,nx);

// Lecture erreurs
if(errok) {
  s.GetLine(strg, 255);
  err2 = new double[nxy];
  for(j=0;j<ny;j++) s.GetR8s(err2+j*nx,nx);
}

// Lecture des projections
if(projx) {
  s.GetLine(strg, 255);
  SetProjX();
  hprojx->Read(s);
}
if(projy) {
  s.GetLine(strg, 255);
  SetProjY();
  hprojy->Read(s);
}

// Lecture des slices
if(nslix>0) {
  s.GetLine(strg, 255);
  SetSliX(nslix);
  DBASSERT (nslix==NSliX());
  for(j=0;j<NSliX();j++) HSliX(j)->Read(s);
}
if(nsliy>0) {
  s.GetLine(strg, 255);
  SetSliY(nsliy);
  DBASSERT (nsliy==NSliY());
  for(j=0;j<NSliY();j++) HSliY(j)->Read(s);
}

// Lecture des bandes
if( nbanx>0 ) {
  s.GetLine(strg, 255);
  for( j=0; j<nbanx; j++) {
    s.GetR4(min);
    s.GetR4(max);
    SetBandX(min,max);
  }
  DBASSERT (nbanx==NBandX());
  for( j=0; j<NBandX(); j++) {
    h = HBandX(j);
    h->Read(s);
  }
}
if( nbany>0 ) {
  s.GetLine(strg, 255);
  for( j=0; j<nbany; j++) {
    s.GetR4(min);
    s.GetR4(max);
    SetBandY(min,max);
  }
  DBASSERT (nbany==NBandY());
  for( j=0; j<NBandY(); j++) {
    h = HBandY(j);
    h->Read(s);
  }
}

return;
}