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

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

#include "machdefs.h"

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

#include "histos2.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(r_8 xMin,r_8 xMax,int_4 nxBin,r_8 yMin,r_8 yMax,int_4 nyBin)
      : mData(new r_8[nxBin*nyBin]), mErr2(NULL)
      , nHist(0), nEntries(0)
      , mNx(nxBin), mNy(nyBin), mNxy(nxBin*nyBin)
      , mXmin(xMin), mXmax(xMax), mYmin(yMin), mYmax(yMax)
      , mWBinx((xMax - xMin)/nxBin), mWBiny((yMax - yMin)/nyBin)
      , mHprojx(NULL), mHprojy(NULL)
{
ASSERT(nxBin>0 && nyBin>0 && xMin<xMax && yMin<yMax);
for(int_4 i=0;i<3;i++) for(int_4 j=0;j<3;j++) mOver[i][j]=0.;
Zero();
mB_s.H = NULL;
}

/*!
  Createur d'un histogramme 2D ayant nxBin,nyBin bins
  entre xMin,xMax et yMin,yMax.
*/
Histo2D::Histo2D(r_4 xMin,r_4 xMax,int_4 nxBin,r_4 yMin,r_4 yMax,int_4 nyBin)
      : mData(new r_8[nxBin*nyBin]), mErr2(NULL)
      , nHist(0), nEntries(0)
      , mNx(nxBin), mNy(nyBin), mNxy(nxBin*nyBin)
      , mXmin((r_8)xMin), mXmax((r_8)xMax), mYmin((r_8)yMin), mYmax((r_8)yMax)
      , mWBinx((xMax - xMin)/nxBin), mWBiny((yMax - yMin)/nyBin)
      , mHprojx(NULL), mHprojy(NULL)
{
ASSERT(nxBin>0 && nyBin>0 && xMin<xMax && yMin<yMax);
for(int_4 i=0;i<3;i++) for(int_4 j=0;j<3;j++) mOver[i][j]=0.;
Zero();
mB_s.H = NULL;
}

/*!
  Constructeur par copie.
*/
Histo2D::Histo2D(const Histo2D& h)
{
int_4 i,j;
mData = new r_8[h.mNxy];
memcpy(mData, h.mData, h.mNxy*sizeof(r_8));

mErr2 = NULL;
if(h.mErr2) {
  mErr2 = new r_8[h.mNxy];
  memcpy(mErr2, h.mErr2, h.mNxy*sizeof(r_8));
}

nHist = h.nHist; nEntries = h.nEntries;
for(i=0;i<3;i++) for(j=0;j<3;j++) mOver[i][j]=h.mOver[i][j];
mNx = h.mNx;  mNy = h.mNy;  mNxy = h.mNxy;
mXmin = h.mXmin; mXmax = h.mXmax; mYmin = h.mYmin; mYmax = h.mYmax;
mWBinx = h.mWBinx; mWBiny = h.mWBiny;
mB_s.H = NULL;

mHprojx = mHprojy = NULL;
if(h.mHprojx) {
  SetProjX();
  *mHprojx = *(h.mHprojx);
}
if(h.mHprojy) {
  SetProjY();
  *mHprojy = *(h.mHprojy); 
}

int_4 nb;
r_8 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));
}

}

/*!
  Constructeur par defaut.
*/
Histo2D::Histo2D()
      : mData(NULL), mErr2(NULL)
      , nHist(0), nEntries(0)
      , mNx(0), mNy(0), mNxy(0)
      , mXmin(0), mXmax(0), mYmin(0), mYmax(0)
      , mWBinx(0), mWBiny(0)
      , mHprojx(NULL), mHprojy(NULL)
{
for(int_4 i=0;i<3;i++) for(int_4 j=0;j<3;j++) mOver[i][j]=0.;
mB_s.H = NULL;
}

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

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

  DelProj();

  DelBandX();
  DelBandY();

  DelSliX();
  DelSliY();

  nHist = 0;
  nEntries = 0;
  mNx = 0; mNy = 0; mNxy = 0;
  mXmin = 0; mXmax = 0; mYmin = 0; mYmax = 0;
  mWBinx = 0; mWBiny = 0;
  for(int_4 i=0;i<3;i++) for(int_4 j=0;j<3;j++) mOver[i][j]=0.;
  mB_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_4 i=0;i<3;i++) for(int_4 j=0;j<3;j++) mOver[i][j]=0.;
  memset(mData, 0, mNxy*sizeof(r_8));
  memset(mOver, 0, 9*sizeof(r_8));

  if( mErr2 != NULL ) memset(mErr2, 0, mNxy*sizeof(r_8));

  ZeroProj();

  ZeroBandX();
  ZeroBandY();

  ZeroSliX();
  ZeroSliY();
}

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

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

  if(this == &h) return *this;
  if( h.mNxy > mNxy ) Delete();
  if(!mData) mData = new r_8[h.mNxy];
  if( !h.mErr2 && mErr2 ) { delete [] mErr2; mErr2=NULL;}
  if( h.mErr2 && !mErr2 ) mErr2 = new r_8[h.mNxy];

  for(i=0;i<3;i++) for(j=0;j<3;j++) mOver[i][j] = h.mOver[i][j];
  nHist = h.nHist;
  nEntries = h.nEntries;
  mNx = h.mNx;  mNy = h.mNy;  mNxy = h.mNxy;
  mXmin = h.mXmin; mXmax = h.mXmax; mWBinx = h.mWBinx;
  mYmin = h.mYmin; mYmax = h.mYmax; mWBiny = h.mWBiny;
  
  memcpy(mData, h.mData, mNxy*sizeof(r_8));
  if(mErr2) memcpy(mErr2, h.mErr2, mNxy*sizeof(r_8));

  DelProjX();
  if(h.mHprojx) {
    SetProjX();
    *mHprojx = *(h.mHprojx);
  }
  DelProjY();
  if(h.mHprojy) {
    SetProjY();
    *mHprojy = *(h.mHprojy);
  }

  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 *= (r_8 b)
{
int_4 i,j;
r_8 b2 = b*b;
for(i=0;i<mNxy;i++) {
  mData[i] *= b;
  if(mErr2) mErr2[i] *= b2;
}
for(i=0;i<3;i++) for(j=0;j<3;j++) mOver[i][j] *= b;
nHist *= b;

if(mHprojx) *mHprojx *= b;
if(mHprojy) *mHprojy *= 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 /= (r_8 b)
{
int_4 i,j;
if (b==0.) throw ParmError("Histo2D::operator / (0) ");
r_8 b2 = b*b;
for(i=0;i<mNxy;i++) {
  mData[i] /= b;
  if(mErr2) mErr2[i] /= b2;
}
for(i=0;i<3;i++) for(j=0;j<3;j++) mOver[i][j] /= b;
nHist /= b;

if(mHprojx) *mHprojx /= b;
if(mHprojy) *mHprojy /= 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 += (r_8 b)
{
int_4 i,j;
r_8 min,max;
for(i=0;i<mNxy;i++) mData[i] += b;
for(i=0;i<3;i++) for(j=0;j<3;j++) mOver[i][j] += b;
nHist += mNxy*b;

if(mHprojx) *mHprojx += b*mNy;
if(mHprojy) *mHprojy += b*mNx;
if(NSliX()>0) for(i=0; i<NSliX();i++) *HSliX(i) += b*mNy/NSliX();
if(NSliY()>0) for(i=0; i<NSliY();i++) *HSliY(i) += b*mNx/NSliY();
if(NBandX()>0) for(i=0; i<NBandX();i++) {
  GetBandX(i,min,max);
  *HBandX(i) += b*(max-min)/(mYmax-mYmin)*mNy;
}
if(NBandY()>0) for(i=0; i<NBandY();i++) {
  GetBandY(i,min,max);
  *HBandY(i) += b*(max-min)/(mXmax-mXmin)*mNx;
}

return *this;
}

/*!
  Operateur H -= b
*/
Histo2D& Histo2D::operator -= (r_8 b)
{
int_4 i,j;
r_8 min,max;
for(i=0;i<mNxy;i++) mData[i] -= b;
for(i=0;i<3;i++) for(j=0;j<3;j++) mOver[i][j] -= b;
nHist -= mNxy*b;

if(mHprojx) *mHprojx -= b*mNy;
if(mHprojy) *mHprojy -= b*mNx;
if(NSliX()>0) for(i=0; i<NSliX();i++) *HSliX(i) -= b*mNy/NSliX();
if(NSliY()>0) for(i=0; i<NSliY();i++) *HSliY(i) -= b*mNx/NSliY();
if(NBandX()>0) for(i=0; i<NBandX();i++) {
  GetBandX(i,min,max);
  *HBandX(i) -= b*(max-min)/(mYmax-mYmin)*mNy;
}
if(NBandY()>0) for(i=0; i<NBandY();i++) {
  GetBandY(i,min,max);
  *HBandY(i) -= b*(max-min)/(mXmax-mXmin)*mNx;
}

return *this;
}

///////////////////////////////////////////////////////////////////
/*!
  Operateur H += H1
*/
Histo2D& Histo2D::operator += (const Histo2D& a)
{
int_4 i,j;
if(mNx!=a.mNx || mNy!=a.mNy) throw SzMismatchError("Histo2D::operator += ");
for(i=0;i<mNxy;i++) {
  mData[i] += a.mData[i];
  if(mErr2 && a.mErr2) mErr2[i] += a.mErr2[i];
}
for(i=0;i<3;i++) for(j=0;j<3;j++) mOver[i][j] += a.mOver[i][j];
nHist += a.nHist;
nEntries += a.nEntries;

if(mHprojx && a.mHprojx) *mHprojx += *(a.mHprojx);
if(mHprojy && a.mHprojy) *mHprojy += *(a.mHprojy);
ZeroSliX();  ZeroSliY();
ZeroBandX(); ZeroBandY();

return *this;
}

/*!
  Operateur H -= H1
*/
Histo2D& Histo2D::operator -= (const Histo2D& a)
{
int_4 i,j;
if(mNx!=a.mNx || mNy!=a.mNy) throw SzMismatchError("Histo2D::operator -= ");
for(i=0;i<mNxy;i++) {
  mData[i] -= a.mData[i];
  if(mErr2 && a.mErr2) mErr2[i] += a.mErr2[i];
}
for(i=0;i<3;i++) for(j=0;j<3;j++) mOver[i][j] += a.mOver[i][j];
nHist -= a.nHist;
nEntries += a.nEntries;

if(mHprojx && a.mHprojx) *mHprojx -= *(a.mHprojx);
if(mHprojy && a.mHprojy) *mHprojy -= *(a.mHprojy);
ZeroSliX();  ZeroSliY();
ZeroBandX(); ZeroBandY();

return *this;
}

/*!
  Operateur H *= H1
*/
Histo2D& Histo2D::operator *= (const Histo2D& a)
{
int_4 i,j;
if(mNx!=a.mNx || mNy!=a.mNy) throw SzMismatchError("Histo2D::operator *= ");
nHist = 0.;
for(i=0;i<mNxy;i++) {
  if(mErr2 && a.mErr2)
      mErr2[i] = a.mData[i]*a.mData[i]*mErr2[i] + mData[i]*mData[i]*a.mErr2[i];
  mData[i] *= a.mData[i];
  nHist += mData[i];
}
for(i=0;i<3;i++) for(j=0;j<3;j++) mOver[i][j] *= a.mOver[i][j];
nEntries += a.nEntries;

if(mHprojx && a.mHprojx) *mHprojx *= *(a.mHprojx);
if(mHprojy && a.mHprojy) *mHprojy *= *(a.mHprojy);
ZeroSliX();  ZeroSliY();
ZeroBandX(); ZeroBandY();

return *this;
}

/*!
  Operateur H /= H1
*/
Histo2D& Histo2D::operator /= (const Histo2D& a)
{
int_4 i,j;
if(mNx!=a.mNx || mNy!=a.mNy) throw SzMismatchError("Histo2D::operator /= ");
nHist = 0.;
for(i=0;i<mNxy;i++) {
  if(a.mData[i]==0.) {
    mData[i]=0.;
    if(mErr2) mErr2[i]=0.;
    continue;
  }
  if(mErr2 && a.mErr2)
      mErr2[i] = (mErr2[i] + mData[i]/a.mData[i]*mData[i]/a.mData[i]*a.mErr2[i])
                /(a.mData[i]*a.mData[i]);
  mData[i] /= a.mData[i];
  nHist += mData[i];
}
for(i=0;i<3;i++) for(j=0;j<3;j++) 
  if(a.mOver[i][j]!=0.) mOver[i][j] *= a.mOver[i][j]; else mOver[i][j] = 0.;
nEntries += a.nEntries;

if(mHprojx && a.mHprojx) *mHprojx /= *(a.mHprojx);
if(mHprojy && a.mHprojy) *mHprojy /= *(a.mHprojy);
ZeroSliX();  ZeroSliY();
ZeroBandX(); ZeroBandY();

return *this;
}

///////////////////////////////////////////////////////////////////
/*!
  Remplissage d'un tableau avec les valeurs des abscisses.
*/
void Histo2D::GetXCoor(TVector<r_8> &v) const
{
r_8 x,y;
v.Realloc(mNx);
for(int_4 i=0;i<mNx;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) const
{
r_8 x,y;
v.Realloc(mNy);
for(int_4 i=0;i<mNy;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) const
{
v.Realloc(mNx,mNy);
for(int_4 i=0;i<mNx;i++)
  for(int_4 j=0;j<mNy;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) const
{
int_4 i,j;
v.Realloc(mNx,mNy);
if(!mErr2)
  {for(i=0;i<mNx;i++) for(j=0;j<mNy;j++) v(i,j) = 0.; return;}
for(i=0;i<mNx;i++) for(j=0;j<mNy;j++) v(i,j) = Error2(i,j);
return;
}

/*!
  Remplissage d'un tableau avec les valeurs des erreurs.
*/
void Histo2D::GetError(TMatrix<r_8> &v) const
{
int_4 i,j;
v.Realloc(mNx,mNy);
if(!mErr2)
  {for(i=0;i<mNx;i++) for(j=0;j<mNy;j++) v(i,j) = 0.; return;}
for(i=0;i<mNx;i++) for(j=0;j<mNy;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_4 ierr)
{
//if(v.NRows()!=(uint_4)mNx || v.NCol()!=(uint_4)mNy) throw SzMismatchError("Histo2D::PutValue()");
uint_4 nnx = (v.NRows()<(uint_4)mNx)? v.NRows(): (uint_4)mNx;
uint_4 nny = (v.NCol() <(uint_4)mNy)? v.NCol() : (uint_4)mNy;
if(nnx>0 && nny>0) for(uint_4 i=0;i<nnx;i++) for(uint_4 j=0;j<nny;j++) {
  (*this)(i,j) = v(i,j);
  if(mErr2 && 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_4 ierr)
{
//if(v.NRows()!=(uint_4)mNx || v.NCol()!=(uint_4)mNy) throw SzMismatchError("Histo2D::PutValueAdd ");
uint_4 nnx = (v.NRows()<(uint_4)mNx)? v.NRows(): (uint_4)mNx;
uint_4 nny = (v.NCol() <(uint_4)mNy)? v.NCol() : (uint_4)mNy;
if(nnx>0 && nny>0) for(uint_4 i=0;i<nnx;i++) for(uint_4 j=0;j<nny;j++) {
  (*this)(i,j) += v(i,j);
  if(mErr2 && 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)
{
//if(v.NRows()!=(uint_4)mNx || v.NCol()!=(uint_4)mNy) throw SzMismatchError("Histo2D::PutError2 ");
uint_4 nnx = (v.NRows()<(uint_4)mNx)? v.NRows(): (uint_4)mNx;
uint_4 nny = (v.NCol() <(uint_4)mNy)? v.NCol() : (uint_4)mNy;
if(nnx>0 && nny>0) {
  if(!mErr2) Errors();
  for(uint_4 i=0;i<nnx;i++) for(uint_4 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)
{
//if(v.NRows()!=(uint_4)mNx || v.NCol()!=(uint_4)mNy) throw SzMismatchError("Histo2D::PutError2Add ");
uint_4 nnx = (v.NRows()<(uint_4)mNx)? v.NRows(): (uint_4)mNx;
uint_4 nny = (v.NCol() <(uint_4)mNy)? v.NCol() : (uint_4)mNy;
if(nnx>0 && nny>0) {
  if(!mErr2) Errors();
  for(uint_4 i=0;i<nnx;i++) for(uint_4 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)
{
//if(v.NRows()!=(uint_4)mNx || v.NCol()!=(uint_4)mNy) throw SzMismatchError("Histo2D::PutError ");
uint_4 nnx = (v.NRows()<(uint_4)mNx)? v.NRows(): (uint_4)mNx;
uint_4 nny = (v.NCol() <(uint_4)mNy)? v.NCol() : (uint_4)mNy;
if(nnx>0 && nny>0) {
  if(!mErr2) Errors();
  for(uint_4 i=0;i<nnx;i++) for(uint_4 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(r_8 x, r_8 y, r_8 w)
{
list<bande_slice>::iterator it;
int_4 i,j;
FindBin(x,y,i,j);

if( mHprojx != NULL ) mHprojx->Add(x,w);
if( mHprojy != NULL ) mHprojy->Add(y,w);

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

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

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

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

if( i<0 || i>=mNx || j<0 || j>=mNy ) {
  if(i<0) i=0; else if(i>=mNx) i=2; else i=1;
  if(j<0) j=0; else if(j>=mNy) j=2; else j=1;
  mOver[i][j] += w;
  mOver[1][1] += w;
  return;
}

mData[j*mNx+i] += w;
if(mErr2!=NULL) mErr2[j*mNx+i] += w*w;
nHist += w;
nEntries++;
}

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

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

r_8 mx=(*this)(il,jl);
for (int_4 i=il; i<=ih; i++)
   for (int_4 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_4& imax,int_4& jmax,int_4 il,int_4 ih,int_4 jl,int_4 jh) const
{
if( il > ih ) { il = 0; ih = mNx-1; }
if( jl > jh ) { jl = 0; jh = mNy-1; }
if( il < 0 ) il = 0;
if( jl < 0 ) jl = 0;
if( ih >= mNx ) ih = mNx-1;
if( jh >= mNy ) jh = mNy-1;

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

r_8 mx=(*this)(il,jl);
for (int_4 i=il; i<=ih; i++)
   for (int_4 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].
*/
r_8 Histo2D::VMax(int_4 il,int_4 ih,int_4 jl,int_4 jh) const
{
if( il > ih ) { il = 0; ih = mNx-1; }
if( jl > jh ) { jl = 0; jh = mNy-1; }
if( il < 0 ) il = 0;
if( jl < 0 ) jl = 0;
if( ih >= mNx ) ih = mNx-1;
if( jh >= mNy ) jh = mNy-1;

r_8 mx=(*this)(il,jl);
if(mNxy==1) return mx;
for (int_4 i=il; i<=ih; i++)
   for (int_4 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].
*/
r_8 Histo2D::VMin(int_4 il,int_4 ih,int_4 jl,int_4 jh) const
{
if( il > ih ) { il = 0; ih = mNx-1; }
if( jl > jh ) { jl = 0; jh = mNy-1; }
if( il < 0 ) il = 0;
if( jl < 0 ) jl = 0;
if( ih >= mNx ) ih = mNx-1;
if( jh >= mNy ) jh = mNy-1;

r_8 mx=(*this)(il,jl);
if(mNxy==1) return mx;
for (int_4 i=il; i<=ih; i++)
   for (int_4 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
     mOver[3][3]:       20 | 21 | 22
                           |    |
                       --------------
                           |    |
                        10 | 11 | 12        11 = all overflow+underflow
                           |    |
                       --------------
                           |    |
                        00 | 01 | 02
  \endverbatim
*/
r_8 Histo2D::NOver(int_4 i,int_4 j) const
{
if( i < 0 || i>=3 || j < 0 || j>=3 ) return mOver[1][1];
return mOver[i][j];
}


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

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

///////////////////////////////////////////////////////////////////
/*!
  Idem EstimeMax(int...) mais retourne x,y.
*/
int_4 Histo2D::EstimeMax(r_8& xm,r_8& ym,int_4 SzPav
                      ,int_4 il,int_4 ih,int_4 jl,int_4 jh) const
{
int_4 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_4 Histo2D::EstimeMax(int_4 im,int_4 jm,r_8& xm,r_8& ym,int_4 SzPav) const
{
xm = ym = 0;
if( SzPav <= 0 ) return -1;
if( im < 0 || im >= mNx ) return -1;
if( jm < 0 || jm >= mNy ) return -1;

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

int_4 rc = 0;
r_8 dxm = 0, dym = 0, wx = 0;
for(int_4 i=im-SzPav;i<=im+SzPav;i++) {
  if( i<0 || i>= mNx ) {rc=1; continue;}
  for(int_4 j=jm-SzPav;j<=jm+SzPav;j++) {
    if( j<0 || j>= mNy ) {rc=1; continue;}
    r_8 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_4 Histo2D::FindMax(int_4& im,int_4& jm,int_4 SzPav,r_8 Dz
                    ,int_4 il,int_4 ih,int_4 jl,int_4 jh) const
{
if( il > ih ) { il = 0; ih = mNx-1; }
if( jl > jh ) { jl = 0; jh = mNy-1; }
if( il < 0 ) il = 0;
if( jl < 0 ) jl = 0;
if( ih >= mNx ) ih = mNx-1;
if( jh >= mNy ) jh = mNy-1;
if( SzPav < 0 ) SzPav = 0;
  else { if( SzPav%2 == 0 ) SzPav++; SzPav = (SzPav-1)/2;}
if( Dz < 0 ) Dz = 0.;
r_8 max = VMax(il,ih,jl,jh) - Dz;
int_4 nmax = 0;
r_8 sumx = -1.e20;
for(int_4 i=il;i<=ih;i++) for(int_4 j=jl;j<=jh;j++) {
  if( (*this)(i,j) < max) continue;
  nmax++;
  r_8 sum = 0.;
  for(int_4 ii=i-SzPav;ii<=i+SzPav;ii++) {
    if( ii<0 || ii >= mNx ) continue;
    for(int_4 jj=j-SzPav;jj<=j+SzPav;jj++) {
      if( jj<0 || jj >= mNy ) continue;
      sum += (*this)(ii,jj);
    }
  }
  if(nmax==1 || sum>sumx) {im=i; jm=j; sumx=sum;}
}
if( nmax <= 0 ) { IJMax(im,jm,il,ih,jl,jh); return 1;}
return nmax;
}

///////////////////////////////////////////////////////////////////
/*!
  Impression des informations sur l'histogramme.
*/
void Histo2D::Show(ostream & os) const
{
  os << "Histo2D::Show    nHist=" << nHist << " nEntries=" << nEntries;
  if(HasErrors()) os << "  Errors=Y \n"; else os << "  Errors=N\n";
  os << "mOver: [ " ;
  for(int j=2; j>=0; j--) {
    for(int i=0; i<3; i++) 
      os << mOver[j][i] << " " ;
    if (j != 0) os << "// ";
  }
  os << "]\n";
  os << "  nx=" << mNx << " xmin=" << mXmin << " xmax=" << mXmax 
     << " binx=" << mWBinx ;
  os << "  ny=" << mNy << " ymin=" << mYmin << " ymax=" << mYmax 
     << " biny=" << mWBiny << endl;
   
  //printf("mOver: [ %g %g %g // %g %g %g  // %g %g %g ]\n"
  //      ,mOver[2][0],mOver[2][1],mOver[2][2]
  //      ,mOver[1][0],mOver[1][1],mOver[1][2]
  //      ,mOver[0][0],mOver[0][1],mOver[0][2]);
  //printf("  nx=%d xmin=%g xmax=%g binx=%g  ",mNx,mXmin,mXmax,mWBinx);
  //printf("  ny=%d ymin=%g ymax=%g biny=%g\n",mNy,mYmin,mYmax,mWBiny);
}

///////////////////////////////////////////////////////////////////
/*!
  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(r_8 min,r_8 max
                   ,int_4 il,int_4 ih,int_4 jl,int_4 jh) const
{
int_4 ns = 35;
const char *s =    "+23456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";

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

PrintStatus();

if( il != 0 || ih != mNx-1 || jl != 0 || jh != mNy-1 ) {
  r_8 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( mNx-1 >= 100 ) {
  printf("     ");
  for(int_4 i=il;i<=ih;i++) printf("%1d",(int_4) (i%1000)/100);
  printf("\n");
}
if( mNx-1 >= 10 ) {
  printf("     ");
  for(int_4 i=il;i<=ih;i++) printf("%1d",(int_4) (i%100)/10);
  printf("\n");
}
printf("     ");
for(int_4 i=il;i<=ih;i++) printf("%1d",i%10);
printf("\n");
printf("     "); {for(int_4 i=il;i<=ih;i++) printf("-"); printf("\n");}

// imprime histogramme
for(int_4 j=jh;j>=jl;j--) {
  printf("%3d: ",j);
  for(int_4 i=il;i<=ih;i++) {
    int_4 h;
    if( 1<=max-min && max-min<=35 ) h = (int_4)( (*this)(i,j) - min ) - 1;
      else h = (int_4)( ((*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_4 i=il;i<=ih;i++) printf("-"); printf("\n");}
if( mNx-1 >= 100 ) {
  printf("     ");
  for(int_4 i=il;i<=ih;i++) printf("%1d",(int_4) (i%1000)/100);
  printf("\n");
}
if( mNx-1 >= 10 ) {
  printf("     ");
  for(int_4 i=il;i<=ih;i++) printf("%1d",(int_4) (i%100)/10);
  printf("\n");
}
printf("     ");
{for(int_4 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( mHprojx != NULL ) DelProjX();
mHprojx = new Histo(mXmin,mXmax,mNx);
if( mErr2 != NULL && mHprojx != NULL ) mHprojx->Errors();
}

/*!
  Pour creer la projection Y.
*/
void Histo2D::SetProjY()
{
if( mHprojy != NULL ) DelProjY();
mHprojy = new Histo(mYmin,mYmax,mNy);
if( mErr2 != NULL && mHprojy != NULL  ) mHprojy->Errors();
}

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

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

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

/*!
  Destruction de l'histogramme de la projection selon X.
*/
void Histo2D::DelProjY()
{
if( mHprojy == NULL ) return;
delete mHprojy;
mHprojy = 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( mHprojx == NULL ) return;
mHprojx->Zero();
}

/*!
  Remise a zero de la projection selon Y.
*/
void Histo2D::ZeroProjY()
{
if( mHprojy == NULL ) return;
mHprojy->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_4 Histo2D::SetBandX(r_8 ybmin,r_8 ybmax)
{
mB_s.num = mLBandx.size();
mB_s.min = ybmin;
mB_s.max = ybmax;
mB_s.H = new Histo(mXmin,mXmax,mNx);
mLBandx.push_back(mB_s);
mB_s.H = NULL;
return mLBandx.size()-1;
}

/*!
  Pour creer une bande en Y entre xbmin et xbmax.
*/
int_4 Histo2D::SetBandY(r_8 xbmin,r_8 xbmax)
{
mB_s.num = mLBandy.size();
mB_s.min = xbmin;
mB_s.max = xbmax;
mB_s.H = new Histo(mYmin,mYmax,mNy);
mLBandy.push_back(mB_s);
mB_s.H = NULL;
return mLBandy.size()-1;
}

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

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

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

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

/*!
  Retourne un pointeur sur la bande numero `n' selon X.
*/
Histo* Histo2D::HBandX(int_4 n) const
{
if( mLBandx.size() <= 0 || n < 0 || n >= (int_4) mLBandx.size() ) return NULL;
for(list<bande_slice>::const_iterator i = mLBandx.begin(); i != mLBandx.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_4 n) const
{
if( mLBandy.size() <= 0 || n < 0 || n >= (int_4) mLBandy.size() ) return NULL;
for(list<bande_slice>::const_iterator i = mLBandy.begin(); i != mLBandy.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_4 n,r_8& ybmin,r_8& ybmax) const
{
ybmin = 0.; ybmax = 0.;
if( mLBandx.size() <= 0 || n < 0 || n >= (int_4) mLBandx.size() ) return;
for(list<bande_slice>::const_iterator i = mLBandx.begin(); i != mLBandx.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_4 n,r_8& xbmin,r_8& xbmax) const
{
xbmin = 0.; xbmax = 0.;
if( mLBandy.size() <= 0 || n < 0 || n >= (int_4) mLBandy.size() ) return;
for(list<bande_slice>::const_iterator i = mLBandy.begin(); i != mLBandy.end(); i++)
  if( (*i).num == n ) { xbmin = (*i).min; xbmax = (*i).max; return;}
return;
}

/*!
  Informations sur les bandes.
*/
void Histo2D::ShowBand(int_4 lp) const
{
  cout << ">>>> Nombre de bande X : " << mLBandx.size() << endl;
if( lp>0 && mLBandx.size()>0 ) {
  list<bande_slice>::const_iterator i;
  for(i = mLBandx.begin(); i != mLBandx.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 : " << mLBandy.size() << endl;
if( lp>0 && mLBandy.size()>0 ) {
  list<bande_slice>::const_iterator i;
  for(i = mLBandy.begin(); i != mLBandy.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_4 Histo2D::SetSliX(int_4 nsli)
{
if( nsli <= 0 ) return -1;
if( nsli >  mNy ) nsli = mNy;
if( mLSlix.size() > 0 ) DelSliX();
r_8 w = (mYmax-mYmin)/nsli;

for(int_4 i=0; i<nsli; i++ ) {
  mB_s.num = i;
  mB_s.min = mYmin + i*w;
  mB_s.max = mB_s.min + w;
  mB_s.H = new Histo(mXmin,mXmax,mNx);
  mLSlix.push_back(mB_s);
  mB_s.H = NULL;
}
return (int_4) mLSlix.size();
}

/*!
  Pour creer `nsli' bandes equidistantes selon Y.
*/
int_4 Histo2D::SetSliY(int_4 nsli)
{
if( nsli <= 0 ) return -1;
if( nsli >  mNx ) nsli = mNx;
if( mLSliy.size() > 0 ) DelSliY();
r_8 w = (mXmax-mXmin)/nsli;

for(int_4 i=0; i<nsli; i++ ) {
  mB_s.num = i;
  mB_s.min = mXmin + i*w;
  mB_s.max = mB_s.min + w;
  mB_s.H = new Histo(mYmin,mYmax,mNy);
  mLSliy.push_back(mB_s);
  mB_s.H = NULL;
}
return (int_4) mLSliy.size();
}

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

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

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

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

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

/*!
  Informations sur les bandes equidistantes.
*/
void Histo2D::ShowSli(int_4 lp) const
{
list<bande_slice>::const_iterator i;
cout << ">>>> Nombre de slice X : " << mLSlix.size() << endl;
if( lp>0 && mLSlix.size() > 0 )
  for(i = mLSlix.begin(); i != mLSlix.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 : " << mLSliy.size() << endl;
if( lp>0 && mLSliy.size()>0 )
  for(i = mLSliy.begin(); i != mLSliy.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 
// --------------------------------------------------------
///////////////////////////////////////////////////////////


DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
void  ObjFileIO<Histo2D>::ReadSelf(PInPersist& is)
{
char strg[256];

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

r_8 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->mNx);
is.Get(dobj->mNy);
is.Get(dobj->mNxy);
is.Get(errok);
is.Get(dobj->nEntries);
is.Get(dobj->nHist);

is.Get(dobj->mXmin);
is.Get(dobj->mXmax);
is.Get(dobj->mYmin);
is.Get(dobj->mYmax);
is.Get(dobj->mWBinx);
is.Get(dobj->mWBiny);

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

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

// Lecture histo2D
dobj->mData = new r_8[dobj->mNxy];
is.GetLine(strg, 255);
{for(int_4 j=0;j<dobj->mNy;j++) is.Get(dobj->mData+j*dobj->mNx,dobj->mNx);}

// Lecture erreurs
if(errok) {
  is.GetLine(strg, 255);
  dobj->mErr2 = new r_8[dobj->mNxy];
  for(int_4 j=0;j<dobj->mNy;j++) is.Get(dobj->mErr2+j*dobj->mNx,dobj->mNx);
}

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

// Lecture des slices
if(nslix>0) {
  is.GetLine(strg, 255);
  dobj->SetSliX(nslix);
  ASSERT (nslix==dobj->NSliX());
  for(int_4 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_4 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_4 j=0; j<nbanx; j++) {
    is.Get(min); is.Get(max);
    dobj->SetBandX(min,max);
  }}
  ASSERT (nbanx==dobj->NBandX());
  {for(int_4 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_4 j=0; j<nbany; j++) {
    is.Get(min); is.Get(max);
    dobj->SetBandY(min,max);
  }}
  ASSERT (nbany==dobj->NBandY());
  {for(int_4 j=0; j<dobj->NBandY(); j++) {
    ObjFileIO<Histo> fio_h(dobj->HBandY(j));
    fio_h.Read(is);
  }}
}

return;
}

DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
void ObjFileIO<Histo2D>::WriteSelf(POutPersist& os) const
{
if (dobj == NULL)   return;
char strg[256];

// Que faut-il ecrire?
int_4 errok = (dobj->mErr2) ? 1 : 0;
int_4 projx = (dobj->mHprojx) ? 1 : 0;
int_4 projy = (dobj->mHprojy) ? 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->mNx,dobj->mNy,dobj->mNxy,errok);
os.PutLine(strg);
sprintf(strg,"nHist=%g nEntries=%d",dobj->nHist,dobj->nEntries);
os.PutLine(strg);
sprintf(strg,"wbinx=%g wbiny=%g",dobj->mWBinx,dobj->mWBiny);
os.PutLine(strg);
sprintf(strg,"xmin=%g xmax=%g ymin=%g ymax=%g"
       ,dobj->mXmin,dobj->mXmax,dobj->mYmin,dobj->mYmax);
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,"mOver %g %g %g %g %g %g %g %g %g"
       ,dobj->mOver[0][0],dobj->mOver[0][1],dobj->mOver[0][2]
       ,dobj->mOver[1][0],dobj->mOver[1][1],dobj->mOver[1][2]
       ,dobj->mOver[2][0],dobj->mOver[2][1],dobj->mOver[2][2]);
os.PutLine(strg);

// Ecriture variables de definitions
os.Put(dobj->mNx);
os.Put(dobj->mNy);
os.Put(dobj->mNxy);
os.Put(errok);
os.Put(dobj->nEntries);
os.Put(dobj->nHist);

os.Put(dobj->mXmin);
os.Put(dobj->mXmax);
os.Put(dobj->mYmin);
os.Put(dobj->mYmax);
os.Put(dobj->mWBinx);
os.Put(dobj->mWBiny);

os.Put(&(dobj->mOver[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->mNxy,dobj->mNx,dobj->mNy);
os.PutLine(strg);
{for(int_4 j=0;j<dobj->mNy;j++) os.Put(dobj->mData+j*dobj->mNx,dobj->mNx);}

// Ecriture erreurs
if(errok) {
  sprintf(strg,"Histo2D: Tableau des erreurs %d = %d * %d"
         ,dobj->mNxy,dobj->mNx,dobj->mNy);
  os.PutLine(strg);
  for(int_4 j=0;j<dobj->mNy;j++) os.Put(dobj->mErr2+j*dobj->mNx,dobj->mNx);
}

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

// Ecriture des slices
if(nslix>0) {
  sprintf(strg,"Histo2D: Slices X %d",nslix);
  os.PutLine(strg);
  for(int_4 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_4 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->mLBandx.begin(); it != dobj->mLBandx.end(); it++) {
    r_8 min = (*it).min; r_8 max = (*it).max;
    os.Put(min); os.Put(max);
  }
  for(it = dobj->mLBandx.begin(); it != dobj->mLBandx.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->mLBandy.begin(); it != dobj->mLBandy.end(); it++) {
    r_8 min = (*it).min; r_8 max = (*it).max;
    os.Put(min); os.Put(max);
  }
  for(it = dobj->mLBandy.begin(); it != dobj->mLBandy.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