source: Sophya/trunk/SophyaLib/HiStats/histos.cc@ 1053

//
// $Id: histos.cc,v 1.6 2000-06-30 13:21:52 ansari Exp $
//

#include "machdefs.h"
#include <string.h>
#include <stdio.h>
#include <math.h>
#include "histos.h"
#include "perrors.h"
#include "poly.h"
#include "strutil.h"
#include "generalfit.h"

/*!
  \class SOPHYA::Histo
  \ingroup HiStats
  Classe d'histogrammes 1D
*/

/********* Methode *********/
/*! Constructeur par defaut */
Histo::Histo()
: data(NULL), err2(NULL),
  under(0), over(0), nHist(0), nEntries(0),
  bins(0), min(0), max(0),
  binWidth(0)
{
 END_CONSTRUCTOR
}

/********* Methode *********/
/*! Constructeur d'un histo de nBin bins allant de xMin a xMax */
Histo::Histo(float xMin, float xMax, int nBin)
: data(new float[nBin]), err2(NULL),
  under(0), over(0), nHist(0), nEntries(0),
  bins(nBin), min(xMin), max(xMax),
  binWidth((max - min)/nBin)
{
  Zero();
  END_CONSTRUCTOR
}

/********* Methode *********/
/*! Constructeur par copie */
Histo::Histo(const Histo& H)
: data(new float[H.bins]), err2(NULL),
  under(H.under), over(H.over), nHist(H.nHist), nEntries(H.nEntries),
  bins(H.bins), min(H.min), max(H.max),
  binWidth(H.binWidth)
{
  memcpy(data, H.data, bins*sizeof(float));
  if( H.err2 != NULL ) {
    err2 = new double[bins];
    memcpy(err2, H.err2, bins*sizeof(double));
  }
  END_CONSTRUCTOR
}

/********* Methode *********/
/*! Gestion de la des-allocation */
void Histo::Delete()
{
  if( data != NULL ) { delete[] data; data = NULL;}
  if( err2 != NULL ) { delete[] err2; err2 = NULL;}
  under = over = min = max = binWidth= 0.;
  nHist = 0.;
  nEntries = bins = 0;
}

/********* Methode *********/
/*! Destructeur */
Histo::~Histo()
{
Delete();
}

/********* Methode *********/
/*!
  Remise a zero
*/
void Histo::Zero()
{
  memset(data, 0, bins*sizeof(float));
  under = over = 0;
  nHist = 0;
  nEntries = 0;
  if( err2 != NULL ) memset(err2, 0, bins*sizeof(double));
}

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

/********* Methode *********/
/*!
  Operateur egal
*/
Histo& Histo::operator = (const Histo& h)
{
  if(this == &h) return *this;
  if( h.bins > bins ) Delete();
  if(!data) data = new float[h.bins];
  if( !h.err2 && err2 ) { delete [] err2; err2=NULL;}
  if( h.err2 && !err2 ) err2 = new double[h.bins];

  under = h.under;
  over = h.over;
  nHist = h.nHist;
  nEntries = h.nEntries;
  bins = h.bins;
  min = h.min;
  max = h.max;
  binWidth = h.binWidth;
  
  memcpy(data, h.data, bins*sizeof(float));
  if(err2) memcpy(err2, h.err2, bins*sizeof(double));

  return *this;
}

/********* Methode *********/
/*!
  Operateur de multiplication par une constante
*/
Histo& Histo::operator *= (double b)
{
double b2 = b*b;
for(int i=0;i<bins;i++) {
  data[i] *= b;
  if(err2) err2[i] *= b2;
}
under *= b;
over *= b;
nHist *= b;

return *this;
}

/*!
  Operateur de division par une constante
*/
Histo& Histo::operator /= (double b)
{
if (b==0.) THROW(inconsistentErr);
double b2 = b*b;
for(int i=0;i<bins;i++) {
  data[i] /= b;
  if(err2) err2[i] /= b2;
}
under /= b;
over /= b;
nHist /= b;

return *this;
}

/*!
  Operateur d'addition d'une constante
*/
Histo& Histo::operator += (double b)
{
for(int i=0;i<bins;i++) data[i] += b;
under += b;
over += b;
nHist += bins*b;

return *this;
}

/*!
  Operateur de soustraction d'une constante
*/
Histo& Histo::operator -= (double b)
{
for(int i=0;i<bins;i++) data[i] -= b;
under -= b;
over -= b;
nHist -= bins*b;

return *this;
}

/********* Methode *********/
/*!
  Operateur H += H1
*/
Histo& Histo::operator += (const Histo& a)
{
if(bins!=a.bins) THROW(sizeMismatchErr);
for(int i=0;i<bins;i++) {
  data[i] += a(i);
  if(err2 && a.err2) err2[i] += a.Error2(i);
}
under += a.under;
over += a.over;
nHist += a.nHist;
nEntries += a.nEntries;

return *this;
}

/*!
  Operateur H -= H1
*/
Histo& Histo::operator -= (const Histo& a)
{
if(bins!=a.bins) THROW(sizeMismatchErr);
for(int i=0;i<bins;i++) {
  data[i] -= a(i);
  if(err2 && a.err2) err2[i] += a.Error2(i);
}
under -= a.under;
over -= a.over;
nHist -= a.nHist;
nEntries += a.nEntries;

return *this;
}

/*!
  Operateur H *= H1
*/
Histo& Histo::operator *= (const Histo& a)
{
if(bins!=a.bins) THROW(sizeMismatchErr);
nHist = 0.;
for(int i=0;i<bins;i++) {
  if(err2 && a.err2)
      err2[i] = a(i)*a(i)*err2[i] + data[i]*data[i]*a.Error2(i);
  data[i] *= a(i);
  nHist += data[i];
}
under *= a.under;
over *= a.over;
nEntries += a.nEntries;

return *this;
}

/*!
  Operateur H /= H1
*/
Histo& Histo::operator /= (const Histo& a)
{
if(bins!=a.bins) THROW(sizeMismatchErr);
nHist = 0.;
for(int i=0;i<bins;i++) {
  if(a(i)==0.) {
    data[i]=0.;
    if(err2) err2[i]=0.;
    continue;
  }
  if(err2 && a.err2)
      err2[i] = (err2[i] + data[i]/a(i)*data[i]/a(i)*a.Error2(i))
                /(a(i)*a(i));
  data[i] /= a(i);
  nHist += data[i];
}
if(a.under!=0.) under /= a.under; else under = 0.;
if(a.over!=0.)  over  /= a.over;  else over  = 0.;
nEntries += a.nEntries;

return *this;
}

/********* Methode *********/
/*!
  Remplissage d'un tableau avec la valeur des abscisses
*/
void Histo::GetAbsc(TVector<r_8> &v)
{
v.Realloc(bins);
for(int i=0;i<bins;i++) v(i) = BinLowEdge(i);
return;
}

/*!
  Remplissage d'un tableau avec la valeur du contenu
*/
void Histo::GetValue(TVector<r_8> &v)
{
v.Realloc(bins);
for(int i=0;i<bins;i++) v(i) = data[i];
return;
}

/*!
  Remplissage d'un tableau avec la valeur des erreurs au carre
*/
void Histo::GetError2(TVector<r_8> &v)
{
v.Realloc(bins);
if(!err2) {for(int i=0;i<bins;i++) v(i) = 0.; return;}
for(int i=0;i<bins;i++) v(i) = err2[i];
return;
}

/*!
  Remplissage d'un tableau avec la valeur des erreurs
*/
void Histo::GetError(TVector<r_8> &v)
{
v.Realloc(bins);
if(!err2) {for(int i=0;i<bins;i++) v(i) = 0.; return;}
for(int i=0;i<bins;i++) v(i) = Error(i);
return;
}

/********* Methode *********/
/*!
  Remplissage du contenu de l'histo avec les valeurs d'un vecteur
*/
void Histo::PutValue(TVector<r_8> &v, int ierr)
{
if(v.NElts()<(uint_4) bins) THROW(sizeMismatchErr);
for(int i=0;i<bins;i++) {
  data[i] = v(i);
  if(err2&&ierr) err2[i] = fabs(v(i));
}
return;
}

/*!
  Addition du contenu de l'histo avec les valeurs d'un vecteur
*/
void Histo::PutValueAdd(TVector<r_8> &v, int ierr)
{
if(v.NElts()<(uint_4) bins) THROW(sizeMismatchErr);
for(int i=0;i<bins;i++) {
  data[i] += v(i);
  if(err2 && ierr) err2[i] += fabs(v(i));
}
return;
}

/*!
  Remplissage des erreurs au carre de l'histo avec les valeurs d'un vecteur
*/
void Histo::PutError2(TVector<r_8> &v)
{
if(v.NElts()<(uint_4) bins) THROW(sizeMismatchErr);
if(!err2) Errors();
for(int i=0;i<bins;i++) err2[i] = v(i);
return;
}

/*!
  Addition des erreurs au carre de l'histo avec les valeurs d'un vecteur
*/
void Histo::PutError2Add(TVector<r_8> &v)
{
if(v.NElts()<(uint_4) bins) THROW(sizeMismatchErr);
if(!err2) Errors();
for(int i=0;i<bins;i++) if(v(i)>0.) err2[i] += v(i);
return;
}

/*!
  Remplissage des erreurs de l'histo avec les valeurs d'un vecteur
*/
void Histo::PutError(TVector<r_8> &v)
{
if(v.NElts()<(uint_4) bins) THROW(sizeMismatchErr);
if(!err2) Errors();
for(int i=0;i<bins;i++)
  if(v(i)>0.) err2[i]=v(i)*v(i); else err2[i]=-v(i)*v(i);
return;
}

/********* Methode *********/
/*!
  Addition du contenu de l'histo pour abscisse x poids w
*/
void Histo::Add(float x, float w)
{
  int numBin = FindBin(x);
  if (numBin<0) under += w;
  else if (numBin>=bins) over += w;
  else {
    data[numBin] += w;
    if(err2!=NULL) err2[numBin] += w*w;
    nHist += w;
    nEntries++;
  }
}

/********* Methode *********/
/*!
  Addition du contenu de l'histo bin numBin poids w
*/
void Histo::AddBin(int numBin, float w)
{
  if (numBin<0) under += w;
  else if (numBin>=bins) over += w;
  else {
    data[numBin] += w;
    if(err2!=NULL) err2[numBin] += w*w;
    nHist += w;
    nEntries++;
  }
}

/********* Methode *********/
/*!
  Remplissage du contenu de l'histo pour abscisse x poids w
*/
void Histo::SetBin(float x, float w)
{
  int numBin = FindBin(x);
  SetBin(numBin,w);
}

/********* Methode *********/
/*!
  Remplissage du contenu de l'histo pour numBin poids w
*/
void Histo::SetBin(int numBin, float w)
{
  if (numBin<0) under = w;
  else if (numBin>=bins) over = w;
  else {
    nHist -= data[numBin];
    data[numBin] = w;
    nHist += w;
  }
}

/********* Methode *********/
/*!
  Remplissage des erreurs au carre pour abscisse x
*/
void Histo::SetErr2(float x, double e2)
{
  int numBin = FindBin(x);
  SetErr2(numBin,e2);
}

/********* Methode *********/
/*!
  Remplissage des erreurs au carre pour numBin poids
*/
void Histo::SetErr2(int numBin, double e2)
{
  if( err2==NULL) return;
  if ( numBin<0 || numBin>=bins ) return;
  err2[numBin] = e2;
}

/********* Methode *********/
/*!
  Remplissage des erreurs pour abscisse x
*/
void Histo::SetErr(float x, float e)
{
SetErr2(x, (double) e*e);
}

/********* Methode *********/
/*!
  Remplissage des erreurs pour numBin
*/
void Histo::SetErr(int numBin, float e)
{
SetErr2(numBin, (double) e*e);
}

/********* Methode *********/
/*!
  Numero du bin ayant le contenu maximum
*/
int Histo::IMax() const
{
  int imx=0;
  if(bins==1) return imx;
  float mx=data[0];
  for (int i=1; i<bins; i++)
    if (data[i]>mx) {imx = i; mx=data[i];}
  return imx;
}

/********* Methode *********/
/*!
  Numero du bin ayant le contenu minimum
*/
int Histo::IMin() const
{
  int imx=0;
  if(bins==1) return imx;
  float mx=data[0];
  for (int i=1; i<bins; i++)
    if (data[i]<mx) {imx = i; mx=data[i];}
  return imx;
}

/********* Methode *********/
/*!
  Valeur le contenu maximum
*/
float Histo::VMax() const
{
  float mx=data[0];
  if(bins==1) return mx;
  for (int i=1;i<bins;i++) if(data[i]>mx) mx=data[i];
  return mx;
}

/********* Methode *********/
/*!
  Valeur le contenu minimum
*/
float Histo::VMin() const
{
  float mx=data[0];
  if(bins==1) return mx;
  for (int i=1;i<bins;i++) if(data[i]<mx) mx=data[i];
  return mx;
}

/********* Methode *********/
/*!
  Valeur moyenne
*/
float Histo::Mean() const
{
  double n = 0;
  double sx = 0;
  for (int i=0; i<bins; i++) {
    double v = (data[i]>=0.) ? data[i] : -data[i];
    n += v;
    sx += BinCenter(i)*v;
  }
  if(n>0.) return sx/n; else return min;
}

/********* Methode *********/
/*!
  Valeur du sigma
*/
float Histo::Sigma() const
{
  double n = 0;
  double sx = 0;
  double sx2= 0;
  for (int i=0; i<bins; i++) {
    double v = (data[i]>=0.) ? data[i] : -data[i];
    n += v;
    sx += BinCenter(i)*v;
    sx2+= BinCenter(i)*BinCenter(i)*v;
  }
  // attention a l'erreur d'arrondi si un seul bin rempli!!
  if( n == 0 ) return 0.;
  sx2 = sx2/n - (sx/n)*(sx/n);
  if( sx2>0. ) return sqrt( sx2 );
      else return 0.;
}

/********* Methode *********/
/*!
  Valeur de la moyenne calculee entre les bin il et ih
*/
float Histo::MeanLH(int il,int ih) const
{
  if( ih < il ) { il = 0; ih = bins-1; }
  if( il < 0 ) il = 0;
  if( ih >= bins ) ih = bins-1;
  double n = 0;
  double sx = 0;
  for (int i=il; i<=ih; i++) {
    double v = (data[i]>=0.) ? data[i] : -data[i];
    n += v;
    sx += BinCenter(i)*v;
  }
  if(n>0.) return sx/n; else return BinLowEdge(il);
}

/********* Methode *********/
/*!
  Valeur de la moyenne calculee entre les bin il et ih
*/
float Histo::SigmaLH(int il,int ih) const
{
  if( ih < il ) { il = 0; ih = bins - 1; }
  if( il < 0 ) il = 0;
  if( ih >= bins ) ih = bins - 1;
  double n = 0;
  double sx = 0;
  double sx2= 0;
  for (int i=il; i<=ih; i++) {
    double v = (data[i]>=0.) ? data[i] : -data[i];
    n += v;
    sx += BinCenter(i)*v;
    sx2+= BinCenter(i)*BinCenter(i)*v;
  }
  if( n == 0 ) return 0.;
  sx2 = sx2/n - (sx/n)*(sx/n);
  if( sx2>0. ) return sqrt( sx2 );
      else return 0.;
}

/********* Methode *********/
/*!
  Valeur de la moyenne calculee entre x0-dx et x0+dx
*/
float Histo::Mean(float x0, float dx) const
{
  double sdata = 0;
  double sx = 0;
  int iMin = FindBin(x0-dx);
  if (iMin<0) iMin = 0;
  int iMax = FindBin(x0+dx);
  if (iMax>bins) iMax=bins;
  for (int i=iMin; i<iMax; i++) {
    double v = (data[i]>=0.) ? data[i] : -data[i];
    sx += BinCenter(i)*v;
    sdata += v;
  }
  if(sdata>0.) return sx/sdata; else return min;
}

/********* Methode *********/
/*!
  Valeur du sigma calcule entre x0-dx et x0+dx
*/
float Histo::Sigma(float x0, float dx) const
{
  double sx = 0;
  double sx2= 0;
  double sdata = 0;
  int iMin = FindBin(x0-dx);
  if (iMin<0) iMin = 0;
  int iMax = FindBin(x0+dx);
  if (iMax>bins) iMax=bins;
  for (int i=iMin; i<iMax; i++) {
    double v = (data[i]>=0.) ? data[i] : -data[i];
    sx += BinCenter(i)*v;
    sx2+= BinCenter(i)*BinCenter(i)*v;
    sdata += v;
  }
  if(sdata>0.) return sqrt( sx2/sdata - (sx/sdata)*(sx/sdata) );
          else return 0.;
}

/********* Methode *********/
/*!
  Valeur de la moyenne et du sigma nettoyes
*/
float Histo::CleanedMean(float& sigma) const
{
  if (!nHist) return 0;
  // on fait ca de facon bete, on pourra raffiner apres...
  float x0 = Mean();
  float s  = Sigma()+binWidth;
  
  for (int i=0; i<3; i++) {
    float xx0 = Mean(x0,5*s);
    s  = Sigma(x0,5*s)+binWidth;
    x0 = xx0;
  }
  sigma = s;
  return x0;
}

/********* Methode *********/
/*!
  Valeur de la moyenne nettoyee
*/
float Histo::CleanedMean() const
{
  if (!nHist) return 0;
  float s=0;
  return CleanedMean(s);
}

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

/********* Methode *********/
/*!
  Retourne le nombre de bins ayant une erreur non-nulle
*/
int Histo::ErrNonNul() const
{
if(err2==NULL) return -1;
int non=0;
for (int i=0;i<bins;i++) if( err2[i] != 0. ) non++;
return non;
}

/********* Methode *********/
/*!
  Renvoie le numero de bin tel que il y ait "per" pourcent entrees
  entre le bin 0 et ce bin (y compris le contenu de ce bin).
*/
int Histo::BinPercent(float per) const
{
double n = nHist*per;
double s = 0.;
int i;
for(i=0; i<bins; i++ ) {
  s += data[i];
  if( s >= n ) break;
}
if( i == bins ) i = bins-1;
return i;
}

/********* Methode *********/
/*!
  Renvoie les numeros de bins imin,imax (0=<i<bins) tels que:
  \verbatim
       notons I = bin contenant l'abscisse x
              N1 = nombre d'entrees entre bin 0 et I compris
              N2 = nombre d'entrees entre bin I et bins-1 compris
       imin = bin tel que nombre d'entrees entre imin et I = N1 * per
       imax = bin tel que nombre d'entrees entre I et imax = N2 * per
     Return: <0 si echec
             min(I-imin,imax-I) si ok
  \endverbatim
*/
int Histo::BinPercent(float x,float per,int& imin,int& imax)
{
imin = imax = -1;
if( per <= 0. ) return -1;

int I = FindBin(x);
if( I<0 || I>=bins ) return -2;

double N1 = 0.; for(int i=0; i<=I;   i++) N1 += data[i]; N1 *= per;
double N2 = 0.; {for(int i=I; i<bins; i++) N2 += data[i]; N2 *= per;}

double n = 0.;
for(imin=I; imin>=0; imin--)   { n += data[imin]; if(n>=N1) break; }
if( imin<0 ) imin = 0;
// cout<<"I="<<I<<" imin="<<imin<<" n="<<n<<" N1="<<N1<<endl;

n = 0.;
for(imax=I; imax<bins; imax++) { n += data[imax]; if(n>=N2) break; }
if( imax>=bins ) imax = bins-1;
// cout<<"I="<<I<<" imax="<<imax<<" n="<<n<<" N2="<<N2<<endl;

return ( imax-I > I-imin ) ? I-imin : imax-I ;
}

/********* Methode *********/
/*!
  Idem precedent mais renvoie xmin et xmax
*/
int Histo::BinPercent(float x,float per,float& xmin,float& xmax)
{
xmin = xmax = 0.;
int imin,imax;
int rc = BinPercent(x,per,imin,imax);
if( rc >= 0 ) { xmin = BinLowEdge(imin); xmax = BinHighEdge(imax); }
return rc;
}

/********* Methode *********/
/*!
  Remplace l'histogramme par son integrale normalise a norm:
  si norm <= 0 : pas de normalisation, integration seule
*/
void Histo::HInteg(float norm)
{
if(bins>1)
  for(int i=1; i<bins; i++) {
    data[i] += data[i-1];
    if(err2!=NULL) err2[i] += err2[i-1];
  }
if( norm <= 0. ) return;
norm /= data[bins-1];
for(int i=0; i<bins; i++) {
  data[i] *= norm;
  if(err2!=NULL) err2[i] *= norm*norm;
}
}

/********* Methode *********/
/*!
  Remplace l'histogramme par sa derivee
*/
void Histo::HDeriv()
{
if( bins <= 1 ) return;
float *temp = new float[bins];
memcpy(temp, data, bins*sizeof(float));
if(bins>=3) for(int i=1; i<bins-1; i++)
              temp[i] = (data[i+1]-data[i-1])/(2.*binWidth);
temp[0] = (data[1]-data[0])/binWidth;
temp[bins-1] = (data[bins-1]-data[bins-2])/binWidth;
memcpy(data, temp, bins*sizeof(float));
delete [] temp;
}

/********* Methode *********/
/*!
  Pour rebinner l'histogramme sur nbinew bins
*/
void Histo::HRebin(int nbinew)
{
  if( nbinew <= 0 ) return;

  // memoraisation de l'histogramme original
  Histo H(*this);

  // Le nombre de bins est il plus grand ??
  if( nbinew > bins ) {
    delete [] data; data = NULL;
    data = new float[nbinew];
    if( err2 != NULL ) {
      delete [] err2; err2 = NULL;
      err2 = new double[nbinew];
    }
  }

  // remise en forme de this
  bins = nbinew;
  this->Zero();
  binWidth = (max-min)/bins;
  // On recopie les parametres qui n'ont pas changes
  under = H.under;
  over  = H.over;


  // Remplissage
  for(int i=0;i<bins;i++) {
    float xmi = BinLowEdge(i);
    float xma = BinHighEdge(i);
    int imi =  H.FindBin(xmi);
    if( imi < 0 ) imi = 0;
    int ima =  H.FindBin(xma);
    if( ima >= H.bins ) ima = H.bins-1;
    double w = 0.;
    if( ima == imi ) {
      w = H(imi) * binWidth/H.binWidth;
    } else {
      w += H(imi) * (H.BinHighEdge(imi)-xmi)/H.binWidth;
      w += H(ima) * (xma -H.BinLowEdge(ima))/H.binWidth;
      if( ima > imi+1 )
          for(int ii=imi+1;ii<ima;ii++) w += H(ii);
    }
    AddBin(i,(float) w);
  }

}

/********* Methode *********/
/*!
  Retourne le nombre de maxima locaux, la valeur du maximum (maxi)
  et sa position (imax), ainsi que la valeur du second maximum
  local (maxn) et sa position (imaxn).
  Attention: si un maximum a la meme valeur sur plusieurs bins
  consecutifs, le bin le plus a droite est pris.
*/
int Histo::MaxiLocal(float& maxi,int& imax,float& maxn,int& imaxn)
{
int nml = 0;
imax = imaxn = -1;
maxi = maxn = -1.f;

bool up = true;
bool down = false;
for(int i=0;i<bins;i++) {
  if( !up ) if( data[i] > data[i-1] ) up = true;
  if( up && !down ) {
    if( i == bins-1 ) down=true;
      else if( data[i] > data[i+1] ) down=true;
  }

  if( up && down ) {
    nml++;
    up = down = false;
    if( imax < 0 ) {
      imax = i; maxi = data[i];
    } else if( data[i] >= maxi ) {
      imaxn = imax; maxn = maxi;
      imax = i; maxi = data[i];
    } else {
      if( imaxn < 0 || data[i] >= maxn ) { imaxn = i; maxn = data[i]; }
    }
  }

}
return nml;
}

/********* Methode *********/
/*!
  Fit de la position du maximum de l'histo par un polynome
  de degre `degree' a `frac' pourcent du maximum.
  L'histo est suppose etre remplit de valeurs positives
*/
float Histo::FitMax(int degree, float frac, int debug) const
{
  if (degree < 2 || degree > 3) degree = 3;
  if (frac <= 0. || frac >= 1.) frac = 0.5;

  if (debug > 1)
    cout<<"Histo::FitMax : Nb Entrees histo ="<<NEntries()<<endl;

  if (NEntries() < 1) THROW(inconsistentErr);

  int iMax = IMax();
  float hmax = (*this)(iMax);
  float xCenter = BinCenter(iMax);

  if(debug>1)
    cout<<"Max histo i="<<iMax<<" x="<<xCenter<<" v="<<hmax<<endl;

  /* find limits at frac*hmax */

  float limit = frac*hmax;

  volatile int iLow = iMax;
  while (iLow>0 && (*this)(iLow)>limit) iLow--;

  volatile int iHigh = iMax;
  while (iHigh<NBins()-1 && (*this)(iHigh)>limit) iHigh++;

  int nLowHigh;
  for(;;) {
    nLowHigh = 0;
    for (int i=iLow; i<=iHigh; i++) if((*this)(i)>0) {
      if(!err2) nLowHigh++;
        else if(Error2(i)>0.) nLowHigh++;
    }
    if (debug > 1) cout <<"Limites histo "<<iLow<<" - "<<iHigh
                        <<" ("<<nLowHigh<<" non nuls)"<<endl;
    if( nLowHigh >= degree+1 ) break;
    iLow--;  iHigh++;
    if(iLow<0 && iHigh>=NBins()) {
      if (debug>1)
        cout<<"Mode histogramme = "<<xCenter
            <<"   BinCenter("<<iMax<<")"<<endl;
      return xCenter;
    }
    if(iLow  < 0       ) iLow  = 0;
    if(iHigh >= NBins()) iHigh = NBins()-1;
  }

  TVector<r_8> xFit(nLowHigh);
  TVector<r_8> yFit(nLowHigh);
  TVector<r_8> e2Fit(nLowHigh);
  TVector<r_8> errcoef(1);
  int ii = 0;
  for (int j=iLow; j<=iHigh; j++) {
    if ((*this)(j)>0) {
      if(!err2) {
        xFit(ii) = BinCenter(j)-xCenter;
        yFit(ii) = (*this)(j);
        e2Fit(ii) = yFit(ii);
        ii++;
      } else if(Error2(j)>0.) {
        xFit(ii) = BinCenter(j)-xCenter;
        yFit(ii) = (*this)(j);
        e2Fit(ii) = Error2(j);
        ii++;
      }
    }
  }
  if(debug>4) {
    int k;
    for(k=0;k<nLowHigh;k++) cout<<" "<<xFit(k);  cout<<endl;
    for(k=0;k<nLowHigh;k++) cout<<" "<<yFit(k);  cout<<endl;
    for(k=0;k<nLowHigh;k++) cout<<" "<<e2Fit(k); cout<<endl;
  }
  if( ii != nLowHigh ) THROW(inconsistentErr);
  Poly pol(degree);
  TRY {
    pol.Fit(xFit, yFit, e2Fit, degree, errcoef);
    if (debug>1) cout << "resultat fit = " << pol << endl;
    pol.Derivate();
  } CATCHALL {
    THROW_SAME;
  } ENDTRY

  int DPolDeg = pol.Degre();
  float fd = 0.;

  if (DPolDeg == 0) {
    // on est dans le cas d'un fit de droite
    if( pol[0] > 0. ) {
      // on a fitte une droite de pente >0.
      fd = xFit(nLowHigh-1) + binWidth/2. + xCenter;
    } else if( pol[0] < 0. ) {
      // on a fitte une droite de pente <0.
      fd = xFit(0) - binWidth/2. + xCenter;
    } else {
      // on a fitte une droite de pente =0. (creneau)
      fd = (xFit(0)+xFit(nLowHigh-1))/2. + xCenter;
    }
  } else if (DPolDeg == 1) {
    // on est dans le cas d'un fit de parabole
    double r=0;
    if (pol.Root1(r)==0) THROW(inconsistentErr);
    fd = r + xCenter;
  } else if (DPolDeg == 2) {
    // on est dans le cas d'un fit de cubique
    double r1=0;
    double r2=0;
    if (pol.Root2(r1,r2) == 0) THROW(inconsistentErr);
    pol.Derivate();
    fd = (pol(r1)<0) ? r1 + xCenter : r2 + xCenter;
  } else {
    // on est dans un cas non prevu
    THROW(inconsistentErr);
  }

  if(fd>max) fd = max;
  if(fd<min) fd = min;

  if (debug>1)
    cout << "Mode histogramme = " << fd
         << "    (DerPol.degre =" << DPolDeg
         << " pol.coeff[0] =" << pol[0]
         << ")" << endl;

  return fd;
}


/********* Methode *********/
/*!
  Calcul de la largeur a frac pourcent du maximum
  autour du bin du maximum.
  L'histo est suppose etre remplit de valeurs positives
*/
float Histo::FindWidth(float frac, int debug) const
{
  float xmax = BinCenter( IMax() );
  return FindWidth(xmax,frac,debug);
}

/********* Methode *********/
/*!
  Calcul de la largeur a frac pourcent de la valeur du bin xmax.
  L'histo est suppose etre remplit de valeurs positives
*/
float Histo::FindWidth(float xmax,float frac, int debug) const
{
  if (frac <= 0 || frac >= 1.) frac = 0.5;

  if (debug > 1)
    cout << "Histo::FindWidth a " << frac
         << " de xmax= " << xmax
         << " , ndata= " << NData()
         << " , nent= " << NEntries()
         << " , nbin= " << NBins() << endl;

  if (NEntries() < 1) THROW(inconsistentErr);
  if (NBins() < 3) THROW(inconsistentErr);

  int iMax = FindBin(xmax);
  if (iMax<0 || iMax>=NBins()) THROW(inconsistentErr);
  float hmax = data[iMax];
  float limit = frac*hmax;
  if (debug > 1)
    cout << "  Max histo[" << iMax << "] = " << hmax
         << ", limite " << limit << endl;

  int iLow = iMax;
  while (iLow>=0 && data[iLow]>limit) iLow--;
  if( iLow < 0 ) iLow = 0;

  int iHigh = iMax;
  while (iHigh<NBins() && data[iHigh]>limit) iHigh++;
  if( iHigh >=NBins() ) iHigh = NBins()-1;

  float xlow   = BinCenter(iLow);
  float ylow   = data[iLow];

  float xlow1=xlow, ylow1=ylow;
  if(iLow+1<NBins()) {
    xlow1  = BinCenter(iLow+1);
    ylow1  = data[iLow+1];
  }

  float xhigh  = BinCenter(iHigh);
  float yhigh  = data[iHigh];

  float xhigh1=xhigh, yhigh1=yhigh;
  if(iHigh-1>=0) {
    xhigh1 = BinCenter(iHigh-1);
    yhigh1 = data[iHigh-1];
  }

  float xlpred,xhpred,wd;

  if(ylow1>ylow) {
    xlpred = xlow + (xlow1-xlow)/(ylow1-ylow)*(limit-ylow);
    if(xlpred < xlow) xlpred = xlow;
  } else xlpred = xlow;

  if(yhigh1>yhigh) {
    xhpred = xhigh + (xhigh1-xhigh)/(yhigh1-yhigh)*(limit-yhigh);
    if(xhpred > xhigh) xhpred = xhigh;
  } else xhpred = xhigh;

  wd = xhpred - xlpred;

  if (debug > 1) {
    cout << "Limites histo: " << "  Width " << wd << endl;
    cout << "  Low: [" << iLow  << "]=" << ylow << " pred-> " << xlpred
         << " - High: [" << iHigh  << "]=" << yhigh << " pred-> " << xhpred
         << endl;
  }

  return wd;
}


/********* Methode *********/
/*!
  Cf suivant mais im est le bin du maximum de l'histo
*/
int Histo::EstimeMax(float& xm,int SzPav)
{
int im = IMax();
return EstimeMax(im,xm,SzPav);
}

/********* Methode *********/
/*!
  Determine l'abscisses du maximum donne par im
  en moyennant dans un pave SzPav  autour du maximum
  \verbatim
     Return:
        0 = si fit maximum reussi avec SzPav pixels
        1 = si fit maximum reussi avec moins que SzPav pixels
        2 = si fit maximum echoue et renvoit BinCenter()
       -1 = si echec: SzPav <= 0 ou im hors limites
  \endverbatim
*/
int Histo::EstimeMax(int& im,float& xm,int SzPav)
{
xm = 0;
if( SzPav <= 0 ) return -1;
if( im < 0 || im >= bins ) return -1;

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

int rc = 0;
double dxm = 0, wx = 0;
for(int i=im-SzPav;i<=im+SzPav;i++) {
  if( i<0 || i>= bins ) {rc=1; continue;}
  dxm += BinCenter(i) * (*this)(i);
  wx  += (*this)(i);
}

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

}

/********* Methode *********/
/*!
  Determine la largeur a frac% du maximum a gauche (widthG)
  et a droite (widthD)
*/
void Histo::EstimeWidthS(float frac,float& widthG,float& widthD)
{
int i;
widthG = widthD = -1.;
if( bins<=1 || frac<=0. || frac>=1. ) return;

int imax = 0;
float maxi = data[0];
for(i=1;i<bins;i++) if(data[i]>maxi) {imax=i; maxi=data[i];}
float xmax = BinCenter(imax);
float maxf = maxi * frac;

// recherche du sigma a gauche
widthG = 0.;
for(i=imax;i>=0;i--) if( data[i] <= maxf ) break;
if(i<0) i=0;
if(i<imax) {
  if( data[i+1] != data[i] ) {
    widthG = BinCenter(i) + binWidth
           * (maxf-data[i])/(data[i+1]-data[i]);
    widthG = xmax - widthG;
  } else widthG = xmax - BinCenter(i);
}

// recherche du sigma a droite
widthD = 0.;
for(i=imax;i<bins;i++) if( data[i] <= maxf ) break;
if(i>=bins) i=bins-1;
if(i>imax) {
  if( data[i] != data[i-1] ) {
    widthD = BinCenter(i) - binWidth
           * (maxf-data[i])/(data[i-1]-data[i]);
    widthD -= xmax;
  } else widthD = BinCenter(i) - xmax;
}

}

//////////////////////////////////////////////////////////
/*!
  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  Histo::Fit(GeneralFit& gfit,unsigned short typ_err)
{
if(NBins()<=0) return -1000;
if(typ_err>5) typ_err=0;

GeneralFitData mydata(1,NBins());

for(int i=0;i<NBins();i++) {
  double x = (double) BinCenter(i);
  double f = (double) (*this)(i);
  double saf = sqrt(fabs((double) f)); if(saf<1.) saf=1.;
  double e=0.;
  if(typ_err==0)      {if(HasErrors()) e=Error(i); else e=1.;}
  else if(typ_err==1) {if(HasErrors()) e=Error(i); 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.AddData1(x,f,e);
}

gfit.SetData(&mydata);

return gfit.Fit();
}

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

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

/********* Methode *********/
/*!
  Impression de l'histogramme dans le fichier fp
  \verbatim
      hdyn = nombre de colonnes pour imprimer les etoiles
          si =0 alors defaut(100),
          si <0 pas de print histo seulement infos
      hmin = minimum de la dynamique
      hmax = maximum de la dynamique
          si hmax<=hmin           : hmin=VMin() hmax=VMax()
          si hmax<=hmin et hmin=0 : hmin=0 hmax=VMax()
          sinon : hmin hmax
      pflag < 0 : on imprime les informations (nbin,min,...) sans l'histogramme
            = 0 : on imprime  "BinCenter(i) data[i]"  (note "... ...")
            bit 0 on : (v=1): numero_bin "... ..."
            bit 1 on : (v=2): "... ..." erreur
  \endverbatim
*/
void Histo::PrintF(FILE * fp, int hdyn,float hmin, float hmax,int pflag,
              int il, int ih)
{

 
 if( il > ih ) { il = 0; ih = bins-1; }
 if( il < 0 ) il = 0;
 if( ih >= bins ) ih = bins-1;

 double dhmin = (double) hmin;
 double dhmax = (double) hmax;
 double hb,hbmn,hbmx;

 if(hdyn==0) hdyn = 100;

 cout << "~Histo::Print  "
      << " nHist=" << nHist << "  nEntries=" << nEntries
      << "  under=" << under << "  over=" << over << endl;
 cout << "       bins=" << bins
      << "  min=" << min << "  max=" << max
      << "  binWidth=" << binWidth << endl;
 cout << "       mean=" << Mean() << " r.m.s=" << Sigma() << endl;

 if(hdyn<0 || pflag<0 ) return;

 if(dhmax<=dhmin) { if(hmin != 0.) dhmin = (double) VMin(); else dhmin=0.;
                  dhmax = (double) VMax(); }
 if(dhmin>dhmax) return;
 if(dhmin==dhmax) {dhmin -= 1.; dhmax += 1.;}
 double wb = (dhmax-dhmin) / (double) hdyn;

 // determination de la position de la valeur zero
 int i0 = (int)(-dhmin/wb);

 // si le zero est dans la dynamique,
 //  il doit imperativement etre une limite de bin
 if( 0 <= i0 && i0 < hdyn ) {
   hbmn = dhmin + i0*wb;
   hbmx = hbmn + wb;
   if( hbmn != 0. ) {
     hbmn *= -1.;
     if( hbmn < hbmx ) {
       // le zero est + pres du bord negatif du bin
       dhmin += hbmn;
       dhmax += hbmn;
     } else {
       // le zero est + pres du bord positif du bin
       dhmin -= hbmx;
       dhmax -= hbmx;
     }
     wb = (dhmax-dhmin) / (double) hdyn;
     i0 = (int)(-dhmin/wb);
   }
 }

 cout <<"       bin minimum="<<dhmin<<"   bin maximum="<<dhmax
      <<"   binw="<<wb<< endl;

 char* s = new char[hdyn+1];
 s[hdyn] = '\0';

 // premiere ligne
 {for(int i=0;i<hdyn;i++) s[i] = '=';}
 if( 0 <= i0 && i0 < hdyn ) s[i0] = '0';
 if(pflag&1) fprintf( fp, "====");
 fprintf( fp, "======================");
 if(pflag&2 && err2!=NULL) fprintf( fp, "===========");
 fprintf( fp, "==%s\n",s);

 // histogramme
 {for(int i=il;i<=ih;i++) {
   for(int k=0;k<hdyn;k++) s[k] = ' ';
   hb = (*this)(i);

   //choix du bin (hdyn bin entre [dhmin,dhmax[
   int ii = (int)( (hb-dhmin)/wb );
   if(ii<0) ii = 0; else if(ii>=hdyn) ii = hdyn-1;

   // limite du bin
   hbmn = dhmin + ii*wb;
   hbmx = hbmn + wb;

   // remplissage de s[] en tenant compte des courbes +/-
   if(i0<0) {
     // courbe entierement positive
     for(int k=0;k<=ii;k++) s[k] = 'X';
   } else if(i0>=hdyn) {
     // courbe entierement negative
     for(int k=hdyn-1;k>=ii;k--) s[k] = 'X';
   } else {
     // courbe positive et negative
     s[i0] = '|';
     if(ii>i0)        for(int k=i0+1;k<=ii;k++) s[k] = 'X';
       else if(ii<i0) for(int k=ii;k<i0;k++)    s[k] = 'X';
          else s[ii] = 'X';
   }

   // valeur a mettre dans le bin le plus haut/bas
   int ib;
   if(hb>0.) ib = (int)( 10.*(hb-hbmn)/(hbmx-hbmn) );
     else if(hb<0.) ib = (int)( 10.*(hbmx-hb)/(hbmx-hbmn) );
       else ib = -1;
   if(ib==-1) s[ii] = '|';
     else if(ib==0) s[ii] = '.';
       else if(ib>0 && ib<10) s[ii] = (char)((int) '0' + ib);

   // traitement des saturations +/-
   if( hb < dhmin ) s[0] = '*';
    else if( hb > dhmax ) s[hdyn-1] = '*';

   if(pflag&1) fprintf( fp, "%3d ",i);
   fprintf( fp, "%10.4g %10.4g ",BinCenter(i),hb);
   if(pflag&2 && err2!=NULL) fprintf( fp, "%10.4g ",Error(i));
   fprintf( fp, "= %s\n",s);
 }}

 // derniere ligne
 for(int i=0;i<hdyn;i++) s[i] = '=';
 if( 0 <= i0 && i0 < hdyn ) s[i0] = '0';
 if(pflag&1) fprintf( fp, "====");
 fprintf( fp, "======================");
 if(pflag&2 && err2!=NULL) fprintf( fp, "===========");
 fprintf( fp, "==%s\n",s);

 // valeur basse des bins (sur ["ndig-1" digits + signe] =  ndig char (>=3))
 const int ndig = 7;
 char sn[2*ndig+10];
 hb = ( fabs(dhmin) > fabs(dhmax) ) ? fabs(dhmin) : fabs(dhmax);
 int n;
 if( hb > 0. ) n = (int)(log10(hb*1.00000001)); else n = 1;
 double scaledig = pow(10.,(double) ndig-2);
 double expo = scaledig/pow(10.,(double) n);
 // cout <<"n="<<n<<" hb="<<hb<<" scaledig="<<scaledig<<" expo="<<expo<<endl;
 for(int k=0;k<ndig;k++) {
   for(int i=0;i<hdyn;i++) s[i] = ' ';
   {for(int i=0;i<hdyn;i++) {
     n = (int)( (dhmin + i*wb)*expo );
     for(int j=0;j<2*ndig+10;j++) sn[j] = ' ';
     sprintf(sn,"%d%c",n,'\0');
     strip(sn,'B',' ');
     // cout <<"n="<<n<<" sn=("<<sn<<")  l="<<strlen(sn)<<" k="<<k;
     if( (int) strlen(sn) > k ) s[i] = sn[k];
     // cout <<"  s=("<<s<<")"<<endl;
   }}
   if(pflag&1) fprintf( fp, "    ");
   fprintf( fp, "                      ");
   if(pflag&2 && err2!=NULL) fprintf( fp, "           ");
   fprintf( fp, "  %s\n",s);
 }
 fprintf( fp, "  (valeurs a multiplier par %.0e)\n",1./expo);

 delete[] s;
}

/********* Methode *********/
/*!
  Impression de l'histogramme sur stdout
*/
void Histo::Print(int hdyn,float hmin, float hmax,int pflag,
             int il, int ih)
{
 Histo::PrintF(stdout, hdyn, hmin, hmax, pflag, il, ih);
}


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

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

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

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

// Lecture des valeurs
is.Get(dobj->bins);
is.Get(dobj->nEntries);
int_4 errok;
is.Get(errok);

is.Get(dobj->binWidth);
is.Get(dobj->min);
is.Get(dobj->max);
is.Get(dobj->under);
is.Get(dobj->over);

is.Get(dobj->nHist);

// Lecture des donnees Histo 1D
dobj->data = new float[dobj->bins];
is.GetLine(strg, 255);
is.Get(dobj->data, dobj->bins);

// Lecture des erreurs
if(errok) {
  is.GetLine(strg, 255);
  dobj->err2 = new double[dobj->bins];
  is.Get(dobj->err2, dobj->bins);
}

return;
}

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

// Que faut-il ecrire?
int_4 errok = (dobj->err2) ? 1 : 0;

// Ecriture entete pour identifier facilement
sprintf(strg,"bins=%6d  NEnt=%15d  errok=%1d",dobj->bins,dobj->nEntries,errok);
os.PutLine(strg);
sprintf(strg,"binw=%g  min=%g max=%g",dobj->binWidth,dobj->min,dobj->max);
os.PutLine(strg);
sprintf(strg, "under=%g over=%g nHist=%g",dobj->under,dobj->over,dobj->nHist);
os.PutLine(strg);

// Ecriture des valeurs
os.Put(dobj->bins);
os.Put(dobj->nEntries);
os.Put(errok);

os.Put(dobj->binWidth);
os.Put(dobj->min);
os.Put(dobj->max);
os.Put(dobj->under);
os.Put(dobj->over);

os.Put(dobj->nHist);

// Ecriture des donnees Histo 1D
sprintf(strg,"Histo: Tableau des donnees %d",dobj->bins);
os.PutLine(strg);
os.Put(dobj->data, dobj->bins);

// Ecriture des erreurs
if(errok) {
  sprintf(strg,"Histo: Tableau des erreurs %d",dobj->bins);
  os.PutLine(strg);
  os.Put(dobj->err2, dobj->bins);
}

return;
}

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

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