Index: /trunk/SophyaLib/HiStats/Makefile
===================================================================
--- /trunk/SophyaLib/HiStats/Makefile	(revision 763)
+++ /trunk/SophyaLib/HiStats/Makefile	(revision 763)
@@ -0,0 +1,103 @@
+include ../Mgr/Makefile.h
+ 
+all: $(LIB)libHiStats.a
+clean:
+	rm -f $(OBJ)hisprof.o $(OBJ)histinit.o $(OBJ)histos.o $(OBJ)histos2.o $(OBJ)ntuple.o $(OBJ)xntuple.o
+	rm -f $(LIB)libHiStats.a
+$(LIB)libHiStats.a : $(OBJ)hisprof.o $(OBJ)histinit.o $(OBJ)histos.o $(OBJ)histos2.o $(OBJ)ntuple.o $(OBJ)xntuple.o
+	$(ARCXX) $(ARCXXFLAGS) $@ $?
+ifeq ($(CXX),cxx)
+	$(CXX) $? $(CPPFLAGS) $(CXXFLAGS) -L$(LIB) -lHiStats -L$(EXTLIBPATH) -o $(OBJ)xx.x -Hf
+endif 
+ 
+$(OBJ)hisprof.o: hisprof.cc $(INC)/machdefs.h \
+ $(INC)/pclassids.h hisprof.h \
+ $(INC)/peida.h \
+ $(INC)/utils.h \
+ $(INC)/perrors.h \
+ $(INC)/pexceptions.h \
+ $(INC)/fmath.h \
+ $(INC)/tvector.h \
+ $(INC)/tmatrix.h \
+ $(INC)/ppersist.h \
+ $(INC)/gnumd5.h \
+ $(INC)/anydataobj.h \
+ $(INC)/ndatablock.h histos.h \
+ $(INC)/objfio.h
+$(OBJ)histinit.o: histinit.cc $(INC)/machdefs.h \
+ $(INC)/pclassids.h histinit.h \
+ $(INC)/ntoolsinit.h \
+ $(INC)/tarrinit.h \
+ $(INC)/sophyainit.h histos.h \
+ $(INC)/objfio.h \
+ $(INC)/anydataobj.h \
+ $(INC)/ppersist.h \
+ $(INC)/pexceptions.h \
+ $(INC)/gnumd5.h \
+ $(INC)/peida.h \
+ $(INC)/utils.h \
+ $(INC)/perrors.h \
+ $(INC)/fmath.h \
+ $(INC)/tvector.h \
+ $(INC)/tmatrix.h \
+ $(INC)/ndatablock.h histos2.h hisprof.h \
+ ntuple.h $(INC)/ntupintf.h \
+ $(INC)/dvlist.h xntuple.h
+$(OBJ)histos.o: histos.cc $(INC)/machdefs.h \
+ $(INC)/pclassids.h histos.h \
+ $(INC)/objfio.h \
+ $(INC)/anydataobj.h \
+ $(INC)/ppersist.h \
+ $(INC)/pexceptions.h \
+ $(INC)/gnumd5.h \
+ $(INC)/peida.h \
+ $(INC)/utils.h \
+ $(INC)/perrors.h \
+ $(INC)/fmath.h \
+ $(INC)/tvector.h \
+ $(INC)/tmatrix.h \
+ $(INC)/ndatablock.h \
+ $(INC)/poly.h \
+ $(INC)/strutil.h \
+ $(INC)/generalfit.h \
+ $(INC)/generaldata.h \
+ $(INC)/ntupintf.h
+$(OBJ)histos2.o: histos2.cc $(INC)/machdefs.h \
+ $(INC)/pclassids.h histos2.h \
+ $(INC)/peida.h \
+ $(INC)/utils.h \
+ $(INC)/perrors.h \
+ $(INC)/pexceptions.h \
+ $(INC)/fmath.h histos.h \
+ $(INC)/objfio.h \
+ $(INC)/anydataobj.h \
+ $(INC)/ppersist.h \
+ $(INC)/gnumd5.h \
+ $(INC)/tvector.h \
+ $(INC)/tmatrix.h \
+ $(INC)/ndatablock.h \
+ $(INC)/generalfit.h \
+ $(INC)/generaldata.h \
+ $(INC)/ntupintf.h \
+ $(INC)/poly.h
+$(OBJ)ntuple.o: ntuple.cc $(INC)/strutil.h \
+ $(INC)/perrors.h \
+ $(INC)/pexceptions.h \
+ $(INC)/machdefs.h \
+ $(INC)/pclassids.h ntuple.h \
+ $(INC)/objfio.h \
+ $(INC)/anydataobj.h \
+ $(INC)/ppersist.h \
+ $(INC)/gnumd5.h \
+ $(INC)/ntupintf.h \
+ $(INC)/dvlist.h
+$(OBJ)xntuple.o: xntuple.cc $(INC)/ppersist.h \
+ $(INC)/machdefs.h \
+ $(INC)/pclassids.h \
+ $(INC)/pexceptions.h \
+ $(INC)/gnumd5.h \
+ $(INC)/perrors.h xntuple.h \
+ $(INC)/objfio.h \
+ $(INC)/anydataobj.h \
+ $(INC)/ntupintf.h \
+ $(INC)/dvlist.h
Index: /trunk/SophyaLib/HiStats/hbook.h
===================================================================
--- /trunk/SophyaLib/HiStats/hbook.h	(revision 763)
+++ /trunk/SophyaLib/HiStats/hbook.h	(revision 763)
@@ -0,0 +1,83 @@
+/* link to hbook     */
+/*      CMV  94      */
+
+#ifndef HBOOK_H_SEEN
+#define HBOOK_H_SEEN
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+void hbandx_(int_4 *,float *,float *,float *);
+void hbandy_(int_4 *,float *,float *,float *);
+
+void hbarx_(int_4 *);
+void hbary_(int_4 *);
+
+void hbfun1_(int_4 *,char *,int_4 *,float *,float *,float fun(float *),int_4);
+void hbfun2_(int_4 *,char *,int_4 *,float *,float *,int_4 *,float *,float *
+            ,float fun(float *,float *),int_4);
+void hbigbi_(int_4 *,int_4 *);
+void hbook1_(int_4 *,char *,int_4 *,float *,float *,float *,int_4);
+void hbook2_(int_4 *,char *,int_4 *,float *,float *
+            ,int_4 *,float *,float *,float *,int_4);
+void hbpro_(int_4 *,float *);
+void hbprof_(int_4 *,char *,int_4 *,float *,float *
+            ,float *,float *,char *,int_4,int_4);
+void hbprox_(int_4 *,float *);
+void hbproy_(int_4 *,float *);
+void hbslix_(int_4 *,int_4 *,float *);
+void hbsliy_(int_4 *,int_4 *,float *);
+void hcdir_(char *,char *,int_4,int_4);
+void hdelet_(int_4 *);
+void hfill_(int_4 *,float *,float *,float *);
+void hfithn_(int_4 *,char *,char *,int_4 *,float *
+            ,float *,float *,float *,float *,float *,int_4,int_4);
+void hfn_(int_4 *,float *);
+void hfunc_(int_4 *,float fun(float *));
+void hgive_(int_4 *,char *,int_4 *,float *,float *
+                          ,int_4 *,float *,float *,int_4 *,int_4 *);
+void hgn_(int_4 *,int_4 *,int_4 *,float *,int_4 *);
+void hgnf_(int_4 *,int_4 *,float *,int_4 *);
+void hgnpar_(int_4 *,char *,int_4);
+float hi_(int_4 *,int_4 *);
+void hidopt_(int_4 *,char *,int_4);
+float hie_(int_4 *,int_4 *);
+float hif_(int_4 *,int_4 *);
+float hij_(int_4 *,int_4 *,int_4 *);
+void hijxy_(int_4 *,int_4 *,int_4 *,float *,float *);
+void hix_(int_4 *,int_4 *,float *);
+void hldir_(char *,char *,int_4,int_4);
+void hmdir_(char *,char *,int_4,int_4);
+float hmax_(int_4 *);
+void hmaxim_(int_4 *,float *);
+float hmin_(int_4 *);
+void hminim_(int_4 *,float *);
+void hnorma_(int_4 *,float *);
+void hnoent_(int_4 *,int_4 *);
+void hpak_(int_4 *,float *);
+void hpake_(int_4 *,float *);
+void hprint_(int_4 *);
+void hrend_(char *,int_4);
+void hrin_(int_4 *,int_4 *,int_4 *);
+float hrndm1_(int_4 *);
+void hrndm2_(int_4 *,float *,float *);
+void hropen_(int_4 *,char *,char *,char *,int_4 *,int_4 *,int_4,int_4,int_4);
+void hrout_(int_4 *,int_4 *,char *,int_4);
+void hscale_(int_4 *,float *);
+float hstati_(int_4 *,int_4 *,char *,int_4 *,int_4);
+float hsum_(int_4 *);
+void hunpak_(int_4 *,float *,char *,int_4 *,int_4);
+void hunpke_(int_4 *,float *,char *,int_4 *,int_4);
+float hx_(int_4 *,float *);
+float hxe_(int_4 *,float *);
+float hxi_(int_4 *,float *,int_4 *);
+float hxyij_(int_4 *,float *,float *,int_4 *,int_4 *);
+float hxy_(int_4 *,float *,float *);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif
+
Index: /trunk/SophyaLib/HiStats/hisprof.cc
===================================================================
--- /trunk/SophyaLib/HiStats/hisprof.cc	(revision 763)
+++ /trunk/SophyaLib/HiStats/hisprof.cc	(revision 763)
@@ -0,0 +1,372 @@
+#include "machdefs.h"
+#include <string.h>
+#include <stdio.h>
+#include <math.h>
+#include "hisprof.h"
+#include "perrors.h"
+
+//++
+// Class	HProf
+// Lib	Outils++ 
+// include	hisprof.h
+//
+//	Classe de profil d'histogrammes.
+//--
+
+//++
+// Titre	Constructeurs
+//--
+
+/********* Methode *********/
+//++
+HProf::HProf()
+//
+//	Constructeur par defaut.
+//--
+: Histo()
+, SumY(NULL), SumY2(NULL), SumW(NULL), Ok(false), YMin(1.), YMax(-1.), Opt(0)
+{
+ END_CONSTRUCTOR
+}
+
+/********* Methode *********/
+//++
+HProf::HProf(float xMin, float xMax, int nBin, float yMin, float yMax)
+//
+//	Constructeur. Histogramme de profil de ``nBin'' bins entre ``xMin''
+//	et ``xMax'' avec coupure d'acceptance sur y entre ``yMin'' et ``yMax''.
+//	Si yMin>=yMax alors pas de coupure d'acceptance sur y.
+//	Par defaut l'erreur du profil represente la dispersion dans le bin,
+//	SetErrOpt(1) permet de demander de calculer l'erreur sur la moyenne.
+//--
+: Histo(xMin,xMax,nBin)
+, SumY(new double[nBin]), SumY2(new double[nBin]), SumW(new double[nBin])
+, Ok(false), YMin(yMin), YMax(yMax), Opt(0)
+{
+  Histo::Errors();
+  Zero();
+  END_CONSTRUCTOR
+}
+
+/********* Methode *********/
+//++
+HProf::HProf(const HProf& H)
+//
+//	Constructeur par copie.
+//--
+: Histo(H)
+, SumY(new double[H.bins]), SumY2(new double[H.bins]), SumW(new double[H.bins])
+, Ok(H.Ok), YMin(H.YMin), YMax(H.YMax), Opt(H.Opt)
+{
+  memcpy(SumY,  H.SumY,  bins*sizeof(double));
+  memcpy(SumY2, H.SumY2, bins*sizeof(double));
+  memcpy(SumW,  H.SumW,  bins*sizeof(double));
+  END_CONSTRUCTOR
+}
+
+/********* Methode *********/
+void HProf::Delete()
+{
+  if( SumY  != NULL ) { delete[] SumY;  SumY = NULL;}
+  if( SumY2 != NULL ) { delete[] SumY2; SumY2 = NULL;}
+  if( SumW  != NULL ) { delete[] SumW;  SumW = NULL;}
+  Ok = false;
+}
+
+/********* Methode *********/
+HProf::~HProf()
+{
+Delete();
+}
+
+//++
+// Titre	Methodes
+//--
+
+/********* Methode *********/
+//++
+void HProf::Zero()
+//
+//	Remise a zero
+//--
+{
+  memset(SumY,  0, bins*sizeof(double));
+  memset(SumY2, 0, bins*sizeof(double));
+  memset(SumW,  0, bins*sizeof(double));
+  Ok = false;
+  Histo::Zero();
+}
+
+/********* Methode *********/
+//++
+void HProf::SetErrOpt(bool spread)
+//
+//	Pour changer l'option de calcul de l'erreur du profile.
+//	Si ``spread'' = true alors l'erreur represente la dispersion
+//	des donnees dans le bin, si = false elle represente l'erreur
+//	sur la moyenne du bin.
+//| - Pour le bin (j):
+//| H(j) = sum(y), E(j) = sum(y^2), L(j) = sum(w)
+//| ->  s(j) = sqrt(E(j)/L(j) - (H(j)/L(j))^2)  dispersion
+//| ->  e(j) = s(j)/sqrt(L(j)) erreur sur la moyenne
+//| spread=true:  opt=0 : dispersion des donnees dans le bin = s(j)
+//| spread=false: opt=1 : erreur sur la moyenne du bin = e(j)
+//--
+{
+int opt = (spread) ? 0 : 1;
+if(opt!=Opt) {Opt=opt; Ok=false;}
+}
+
+/********* Methode *********/
+//++
+void HProf::UpdateHisto() const
+//
+//	Pour mettre a jour l'histogramme de profil.
+//	Il est important d'appeler cette methode si on veut
+//	utiliser les methodes de la classe Histo qui ne sont
+//	pas redefinies dans la classe HProf.
+//	En effet, pour des raisons de precision la classe
+//	HProf travaille avec des tableaux en double precision
+//	et seulement au moment ou l'on a besoin de l'histogramme
+//	ce dernier est rempli avec les valeurs demandees (moyenne
+//	et dispersion/erreur sur la moyenne).
+//--
+{
+
+float m,e2;
+for(int i=0;i<bins;i++) {
+  if(SumW[i]<=0.) {
+    m = e2 = 0.;
+  } else {
+    m  = SumY[i]/SumW[i];
+    e2 = SumY2[i]/SumW[i] - m*m;
+    if(Opt) e2 /= SumW[i];
+  }
+  data[i] = m;
+  err2[i] = e2;
+}
+// Attention, a cause de "WriteSelf const" UpdateHisto doit etre "const".
+// Comme on veut modifier Ok, on est oblige de faire cette entourloupe:
+HProf *buff = (HProf *) this;
+buff->Ok = true;
+}
+
+/********* Methode *********/
+//++
+void HProf::Add(float x, float y, float w)
+//
+//	Addition au contenu de l'histo pour abscisse x de la valeur y et poids w
+//--
+{
+  if(YMax>YMin && (y<YMin || YMax<y)) return;
+  int numBin = FindBin(x);
+  if (numBin<0) under += w;
+  else if (numBin>=bins) over += w;
+  else {
+    Ok = false;
+    SumY[numBin]  += y;
+    SumY2[numBin] += y*y;
+    SumW[numBin]  += w;
+    nHist += w;
+    nEntries++;
+  }
+}
+
+/********* Methode *********/
+//++
+void HProf::AddBin(int numBin, float y, float w)
+//
+//	Addition au contenu de l'histo bin numBin de la valeur y poids w
+//--
+{
+  if(YMax>YMin && (y<YMin || YMax<y)) return;
+  if (numBin<0) under += w;
+  else if (numBin>=bins) over += w;
+  else {
+    Ok = false;
+    SumY[numBin]  += y;
+    SumY2[numBin] += y*y;
+    SumW[numBin]  += w;
+    nHist += w;
+    nEntries++;
+  }
+}
+
+/********* Methode *********/
+//++
+HProf& HProf::operator = (const HProf& h)
+//
+//--
+{
+if(this == &h) return *this;
+if( h.bins > bins ) Delete();
+Histo *hthis = (Histo *) this;
+*hthis = (Histo) h;
+if(!SumY)  SumY  = new double[bins];
+if(!SumY2) SumY2 = new double[bins];
+if(!SumW)  SumW  = new double[bins];
+memcpy(SumY,  h.SumY,  bins*sizeof(double));
+memcpy(SumY2, h.SumY2, bins*sizeof(double));
+memcpy(SumW,  h.SumW,  bins*sizeof(double));
+Ok = h.Ok;
+YMin = h.YMin;
+YMax = h.YMax;
+Opt = h.Opt;
+
+return *this;
+}
+
+/********* Methode *********/
+//++
+HProf& HProf::operator += (const HProf& a)
+//
+//	Attention dans cette addition il n'y a pas de gestion
+//	des YMin et YMax. L'addition est faite meme si les limites
+//	en Y de ``a'' sont differentes de celles de ``this''.
+//--
+{
+if(bins!=a.bins) THROW(sizeMismatchErr);
+Histo *hthis = (Histo *) this;
+*hthis += (Histo) a;
+for(int i=0;i<bins;i++) {
+  SumY[i]  += a.SumY[i];
+  SumY2[i] += a.SumY2[i];
+  SumW[i]  += a.SumW[i];
+}
+Ok = false;
+
+return *this;
+}
+
+/********* Methode *********/
+//++
+HProf operator + (const HProf& a, const HProf& b)
+//
+//	Meme commentaire que pour l'operateur +=
+//--
+{
+if (b.NBins()!=a.NBins()) THROW(sizeMismatchErr);
+HProf c(a);
+return (c += b);
+}
+
+// Rappel des inlines functions pour commentaires
+//++
+//   inline Histo GetHisto()
+//	Retourne l'histogramme de profil.
+//--
+//++
+// inline void GetMean(Vector& v)
+//	Retourne le contenu de la moyenne dans le vecteur v
+//--
+//++
+// inline void GetError2(Vector& v)
+//	Retourne le contenu au carre de la dispersion/erreur dans le vecteur v
+//--
+//++
+// inline float operator()(int i) const
+//	Retourne le contenu du bin i
+//--
+//++
+// inline float Error2(int i) const
+//	Retourne le carre de la dispersion/erreur du bin i
+//--
+//++
+// inline float Error(int i) const
+//	Retourne la dispersion/erreur du bin i
+//--
+//++
+// inline int Fit(GeneralFit& gfit)
+//	Fit du profile par ``gfit''.
+//--
+//++
+// inline Histo* FitResidus(GeneralFit& gfit)
+//	Retourne l'Histogramme des residus par ``gfit''.
+//--
+//++
+// inline Histo* FitFunction(GeneralFit& gfit)
+//	Retourne l'Histogramme de la fonction fittee par ``gfit''.
+//--
+//++
+// inline void Print(int dyn,float hmin,float hmax,int pflag,int il,int ih)
+//	Print, voir detail dans Histo::Print
+//--
+
+
+///////////////////////////////////////////////////////////
+// --------------------------------------------------------
+//   Les objets delegues pour la gestion de persistance 
+// --------------------------------------------------------
+///////////////////////////////////////////////////////////
+
+
+void ObjFileIO<HProf>::ReadSelf(PInPersist& is)
+{
+char strg[256];
+
+if(dobj==NULL) dobj=new HProf;
+  else         dobj->Delete();
+
+// Lecture entete
+is.GetLine(strg,255);
+
+// Ecriture des valeurs
+is.Get(dobj->bins);
+is.Get(dobj->YMin);
+is.Get(dobj->YMax);
+is.Get(dobj->Opt);
+dobj->Ok = true;
+
+// Ecriture des donnees propres a l'histogramme de profil.
+is.GetLine(strg,255);
+dobj->SumY  = new double[dobj->bins];
+dobj->SumY2 = new double[dobj->bins];
+dobj->SumW  = new double[dobj->bins];
+is.Get(dobj->SumY,  dobj->bins);
+is.Get(dobj->SumY2, dobj->bins);
+is.Get(dobj->SumW,  dobj->bins);
+
+// Ecriture de l'histogramme
+is >> (Histo&)(*dobj);
+return;
+}
+
+void ObjFileIO<HProf>::WriteSelf(POutPersist& os) const
+{
+if (dobj == NULL)   return;
+char strg[256];
+
+dobj->UpdateHisto();
+
+// Ecriture entete pour identifier facilement
+sprintf(strg,"HProf: YMin=%f  YMax=%f  Opt=%1d",dobj->YMin,dobj->YMax,dobj->Opt);
+os.PutLine(strg);
+
+// Ecriture des valeurs
+os.Put(dobj->bins);
+os.Put(dobj->YMin);
+os.Put(dobj->YMax);
+os.Put(dobj->Opt);
+
+// Ecriture des donnees propres a l'histogramme de profil.
+sprintf(strg,"HProf: SumY SumY2 SumW");
+os.PutLine(strg);
+os.Put(dobj->SumY,  dobj->bins);
+os.Put(dobj->SumY2, dobj->bins);
+os.Put(dobj->SumW,  dobj->bins);
+
+// Ecriture de l'histogramme
+// FIO_Histo fio_h((Histo&)*dobj);
+// fio_h.Write(os);
+os << (Histo&)(*dobj);
+
+return;
+}
+
+#ifdef __CXX_PRAGMA_TEMPLATES__
+#pragma define_template ObjFileIO<HProf>
+#endif
+
+#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
+template class ObjFileIO<HProf>;
+#endif
Index: /trunk/SophyaLib/HiStats/hisprof.h
===================================================================
--- /trunk/SophyaLib/HiStats/hisprof.h	(revision 763)
+++ /trunk/SophyaLib/HiStats/hisprof.h	(revision 763)
@@ -0,0 +1,78 @@
+#ifndef HISPROF_SEEN
+#define HISPROF_SEEN
+
+#include <stdio.h>
+#include "peida.h"
+#include "tvector.h"
+#include "ppersist.h"
+#include "histos.h"
+
+namespace SOPHYA {
+
+class HProf : public Histo {
+  friend class ObjFileIO<HProf>;
+public:
+
+  // CREATOR / DESTRUCTOR
+  HProf();
+  HProf(float xMin, float xMax, int nBin=100, float yMin=1., float yMax=-1.);
+  HProf(const HProf& H);
+  virtual ~HProf();
+
+  // UPDATING or SETTING
+  void   Zero();
+  void   UpdateHisto() const;
+  void   SetErrOpt(bool spread = true);
+  void   Add(float x, float y, float w = 1.);
+  void   AddBin(int numBin, float y, float w = 1.);
+  
+  // Acces a l information
+  inline Histo GetHisto() {if(!Ok) UpdateHisto(); return (Histo) *this;}
+  inline void GetMean(Vector& v) {if(!Ok) UpdateHisto(); Histo::GetValue(v);}
+  inline void GetError2(Vector& v) {if(!Ok) UpdateHisto(); Histo::GetError2(v);}
+  inline float operator()(int i) const {if(!Ok) UpdateHisto(); return data[i];}
+  inline float Error2(int i) const {if(!Ok) UpdateHisto(); return (float) err2[i];}
+  inline float Error(int i) const
+         {if(!Ok) UpdateHisto(); return err2[i]>0. ? (float) sqrt(err2[i]) : 0.f;}
+
+  // Operators
+  HProf& operator = (const HProf& h);
+  HProf& operator += (const HProf& a);
+  friend HProf operator + (const HProf& a, const HProf& b);
+
+  // Fit
+  inline int    Fit(GeneralFit& gfit)
+         {if(!Ok) UpdateHisto(); return Histo::Fit(gfit,0);}
+  inline Histo  FitResidus(GeneralFit& gfit)
+         {if(!Ok) UpdateHisto(); return Histo::FitResidus(gfit);}
+  inline Histo  FitFunction(GeneralFit& gfit)
+         {if(!Ok) UpdateHisto(); return Histo::FitFunction(gfit);}
+
+  // Print
+  inline void Print(int dyn=100,float hmin=1.,float hmax=-1.,int pflag=0,int il=1,int ih=-1)
+                           {if(!Ok) UpdateHisto(); Histo::Print(dyn,hmin,hmax,pflag,il,ih);}
+
+protected:
+  void Delete();
+
+  double*        SumY;
+  double*        SumY2;
+  double*        SumW;
+  bool           Ok;
+  float          YMin;
+  float          YMax;
+  uint_2         Opt;
+};
+
+
+inline POutPersist& operator << (POutPersist& os, HProf & obj)
+{ ObjFileIO<HProf> fio(&obj);  fio.Write(os);  return(os); }
+inline PInPersist& operator >> (PInPersist& is, HProf & obj)
+{ ObjFileIO<HProf> fio(&obj);  fio.Read(is);  return(is); }
+// Classe pour la gestion de persistance
+// ObjFileIO<HProf>
+
+
+} // Fin du namespace
+
+#endif // HISPROF_SEEN
Index: /trunk/SophyaLib/HiStats/histinit.cc
===================================================================
--- /trunk/SophyaLib/HiStats/histinit.cc	(revision 763)
+++ /trunk/SophyaLib/HiStats/histinit.cc	(revision 763)
@@ -0,0 +1,41 @@
+// Classe d'initialisation du module HiStats 
+//                        (Histogram and Statistics)
+
+#include "machdefs.h"
+#include "histinit.h"
+
+#include "histos.h"
+#include "histos2.h"
+#include "hisprof.h"
+#include "ntuple.h"
+#include "xntuple.h"
+
+int HiStatsInitiator::FgInit = 0;
+
+HiStatsInitiator::HiStatsInitiator()
+  : NToolsInitiator()
+{
+  FgInit++;
+  if (FgInit > 1)  return;
+
+//   Enregistrement des classes PPersist du modules HiStats
+  
+  PPRegister(ObjFileIO<Histo>);
+  PPRegister(ObjFileIO<Histo2D>);
+  PPRegister(ObjFileIO<HProf>);
+
+  PPRegister(ObjFileIO<NTuple>);
+  PPRegister(ObjFileIO<XNTuple>);
+
+}
+
+HiStatsInitiator::~HiStatsInitiator()
+{
+  FgInit--;
+}
+
+
+// On met un objet initiator en statique, pour les loaders qui savent 
+// appeler le constructeur des objets statiques   Reza 08/98
+static HiStatsInitiator histatsinit;
+
Index: /trunk/SophyaLib/HiStats/histinit.h
===================================================================
--- /trunk/SophyaLib/HiStats/histinit.h	(revision 763)
+++ /trunk/SophyaLib/HiStats/histinit.h	(revision 763)
@@ -0,0 +1,27 @@
+// This may look like C code, but it is really -*- C++ -*-
+// Classe d'initialiseur du module HiStats (Histrogram/Stats)
+
+#ifndef HISTATSINIT_H_SEEN
+#define HISTATSINIT_H_SEEN
+
+
+// Classe d''initialisation pour le module de NTools
+#include "ntoolsinit.h"
+
+namespace SOPHYA {
+
+class HiStatsInitiator : public NToolsInitiator {
+private:
+  static int FgInit;
+public:
+  			HiStatsInitiator();
+  virtual		~HiStatsInitiator();
+};
+
+#undef SophyaInit
+#define SophyaInit() HiStatsInitiator histatinitiator
+
+} // Fin du namespace
+
+#endif
+
Index: /trunk/SophyaLib/HiStats/histos.cc
===================================================================
--- /trunk/SophyaLib/HiStats/histos.cc	(revision 763)
+++ /trunk/SophyaLib/HiStats/histos.cc	(revision 763)
@@ -0,0 +1,1827 @@
+//
+// $Id: histos.cc,v 1.1.1.1 2000-03-02 16:51:32 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	Histo
+// Lib	Outils++ 
+// include	histos.h
+//
+//	Classe d'histogrammes 1D
+//--
+
+//++
+// Titre	Constructeurs
+//--
+
+/********* Methode *********/
+//++
+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 *********/
+//++
+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 *********/
+//++
+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 *********/
+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 *********/
+Histo::~Histo()
+{
+Delete();
+}
+
+//++
+// Titre	Methodes
+//--
+
+/********* Methode *********/
+//++
+void Histo::Zero()
+//
+//	Remise a zero
+//--
+{
+  memset(data, 0, bins*sizeof(float));
+  under = over = 0;
+  nHist = 0;
+  nEntries = 0;
+  if( err2 != NULL ) memset(err2, 0, bins*sizeof(double));
+}
+
+/********* Methode *********/
+//++
+void Histo::Errors()
+//
+//	Pour avoir le calcul des erreurs
+//--
+{
+ if( bins > 0 ) {
+   if(err2==NULL) err2 = new double[bins];
+   memset(err2, 0, bins*sizeof(double));
+ }
+}
+
+/********* Methode *********/
+//++
+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 *********/
+//++
+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;
+}
+
+//++
+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;
+}
+
+//++
+Histo& Histo::operator += (double b)
+//
+//--
+{
+for(int i=0;i<bins;i++) data[i] += b;
+under += b;
+over += b;
+nHist += bins*b;
+
+return *this;
+}
+
+//++
+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 *********/
+//++
+Histo operator * (const Histo& a, double b)
+//
+//--
+{
+  Histo result(a);
+  return (result *= b);
+}
+
+//++
+Histo operator * (double b, const Histo& a)
+//
+//--
+{
+  Histo result(a);
+  return (result *= b);
+}
+
+//++
+Histo operator / (const Histo& a, double b)
+//
+//--
+{
+  Histo result(a);
+  return (result /= b);
+}
+
+//++
+Histo operator + (const Histo& a, double b)
+//
+//--
+{
+  Histo result(a);
+  return (result += b);
+}
+
+//++
+Histo operator + (double b, const Histo& a)
+//
+//--
+{
+  Histo result(a);
+  return (result += b);
+}
+
+//++
+Histo operator - (const Histo& a, double b)
+//
+//--
+{
+  Histo result(a);
+  return (result -= b);
+}
+
+//++
+Histo operator - (double b, const Histo& a)
+//
+//--
+{
+  Histo result(a);
+  result *= -1.;
+  return (result += b);
+}
+
+/********* Methode *********/
+//++
+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);
+    else err2[i] = 0.;
+}
+under += a.under;
+over += a.over;
+nHist += a.nHist;
+nEntries += a.nEntries;
+
+return *this;
+}
+
+//++
+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);
+    else err2[i] = 0.;
+}
+under -= a.under;
+over -= a.over;
+nHist -= a.nHist;
+nEntries += a.nEntries;
+
+return *this;
+}
+
+//++
+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);
+    else err2[i] = 0.;
+  data[i] *= a(i);
+  nHist += data[i];
+}
+under *= a.under;
+over *= a.over;
+nEntries += a.nEntries;
+
+return *this;
+}
+
+//++
+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));
+    else err2[i] = 0.;
+  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 *********/
+//++
+Histo operator + (const Histo& a, const Histo& b)
+//
+//--
+{
+if (b.NBins()!=a.NBins()) THROW(sizeMismatchErr);
+Histo c(a);
+return (c += b);
+}
+
+//++
+Histo operator - (const Histo& a, const Histo& b)
+//
+//--
+{
+if (b.NBins()!=a.NBins()) THROW(sizeMismatchErr);
+Histo c(a);
+return (c -= b);
+}
+
+//++
+Histo operator * (const Histo& a, const Histo& b)
+//
+//--
+{
+if (b.NBins()!=a.NBins()) THROW(sizeMismatchErr);
+Histo c(a);
+return (c *= b);
+}
+
+//++
+Histo operator / (const Histo& a, const Histo& b)
+//
+//--
+{
+if (b.NBins()!=a.NBins()) THROW(sizeMismatchErr);
+Histo c(a);
+return (c /= b);
+}
+
+/********* Methode *********/
+//++
+void Histo::GetAbsc(Vector &v)
+//
+//	Remplissage d'un tableau avec la valeur des abscisses
+//--
+{
+v.Realloc(bins);
+for(int i=0;i<bins;i++) v(i) = BinLowEdge(i);
+return;
+}
+
+//++
+void Histo::GetValue(Vector &v)
+//
+//	Remplissage d'un tableau avec la valeur du contenu
+//--
+{
+v.Realloc(bins);
+for(int i=0;i<bins;i++) v(i) = data[i];
+return;
+}
+
+//++
+void Histo::GetError2(Vector &v)
+//
+//	Remplissage d'un tableau avec la valeur des erreurs au carre
+//--
+{
+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;
+}
+
+//++
+void Histo::GetError(Vector &v)
+//
+//	Remplissage d'un tableau avec la valeur des erreurs
+//--
+{
+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 *********/
+//++
+void Histo::PutValue(Vector &v, int ierr)
+//
+//	Remplissage du contenu de l'histo avec les valeurs d'un vecteur
+//--
+{
+if(v.NElts()<bins) THROW(sizeMismatchErr);
+for(int i=0;i<bins;i++) {
+  data[i] = v(i);
+  if(err2&&ierr) err2[i] = fabs(v(i));
+}
+return;
+}
+
+//++
+void Histo::PutValueAdd(Vector &v, int ierr)
+//
+//	Addition du contenu de l'histo avec les valeurs d'un vecteur
+//--
+{
+if(v.NElts()<bins) THROW(sizeMismatchErr);
+for(int i=0;i<bins;i++) {
+  data[i] += v(i);
+  if(err2&&ierr) err2[i] += fabs(v(i));
+}
+return;
+}
+
+//++
+void Histo::PutError2(Vector &v)
+//
+//	Remplissage des erreurs au carre de l'histo avec les valeurs d'un vecteur
+//--
+{
+if(v.NElts()<bins) THROW(sizeMismatchErr);
+if(!err2) Errors();
+for(int i=0;i<bins;i++) err2[i] = v(i);
+return;
+}
+
+//++
+void Histo::PutError2Add(Vector &v)
+//
+//	Addition des erreurs au carre de l'histo avec les valeurs d'un vecteur
+//--
+{
+if(v.NElts()<bins) THROW(sizeMismatchErr);
+if(!err2) Errors();
+for(int i=0;i<bins;i++) if(v(i)>0.) err2[i] += v(i);
+return;
+}
+
+//++
+void Histo::PutError(Vector &v)
+//
+//	Remplissage des erreurs de l'histo avec les valeurs d'un vecteur
+//--
+{
+if(v.NElts()<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 *********/
+//++
+void Histo::Add(float x, float w)
+//
+//	Addition du contenu de l'histo pour abscisse x poids 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 *********/
+//++
+void Histo::AddBin(int numBin, float w)
+//
+//	Addition du contenu de l'histo bin numBin poids 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 *********/
+//++
+void Histo::SetBin(float x, float w)
+//
+//	Remplissage du contenu de l'histo pour abscisse x poids w
+//--
+{
+  int numBin = FindBin(x);
+  SetBin(numBin,w);
+}
+
+/********* Methode *********/
+//++
+void Histo::SetBin(int numBin, float w)
+//
+//	Remplissage du contenu de l'histo pour numBin poids w
+//--
+{
+  if (numBin<0) under = w;
+  else if (numBin>=bins) over = w;
+  else {
+    nHist -= data[numBin];
+    data[numBin] = w;
+    nHist += w;
+  }
+}
+
+/********* Methode *********/
+//++
+void Histo::SetErr2(float x, double e2)
+//
+//	Remplissage des erreurs au carre pour abscisse x
+//--
+{
+  int numBin = FindBin(x);
+  SetErr2(numBin,e2);
+}
+
+/********* Methode *********/
+//++
+void Histo::SetErr2(int numBin, double e2)
+//
+//	Remplissage des erreurs au carre pour numBin poids
+//--
+{
+  if( err2==NULL) return;
+  if ( numBin<0 || numBin>=bins ) return;
+  err2[numBin] = e2;
+}
+
+/********* Methode *********/
+//++
+void Histo::SetErr(float x, float e)
+//
+//	Remplissage des erreurs pour abscisse x
+//--
+{
+SetErr2(x, (double) e*e);
+}
+
+/********* Methode *********/
+//++
+void Histo::SetErr(int numBin, float e)
+//
+//	Remplissage des erreurs pour numBin
+//--
+{
+SetErr2(numBin, (double) e*e);
+}
+
+/********* Methode *********/
+//++
+int Histo::IMax() const
+//
+//	Numero du bin ayant le contenu maximum
+//--
+{
+  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 *********/
+//++
+int Histo::IMin() const
+//
+//	Numero du bin ayant le contenu minimum
+//--
+{
+  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 *********/
+//++
+float Histo::VMax() const
+//
+//	Valeur le contenu maximum
+//--
+{
+  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 *********/
+//++
+float Histo::VMin() const
+//
+//	Valeur le contenu minimum
+//--
+{
+  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 *********/
+//++
+float Histo::Mean() const
+//
+//	Valeur moyenne
+//--
+{
+  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 *********/
+//++
+float Histo::Sigma() const
+//
+//	Valeur du sigma
+//--
+{
+  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 *********/
+//++
+float Histo::MeanLH(int il,int ih) const
+//
+//	Valeur de la moyenne calculee entre les bin il et ih
+//--
+{
+  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 *********/
+//++
+float Histo::SigmaLH(int il,int ih) const
+//
+//	Valeur de la moyenne calculee entre les bin il et ih
+//--
+{
+  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 *********/
+//++
+float Histo::Mean(float x0, float dx) const
+//
+//	Valeur de la moyenne calculee entre x0-dx et x0+dx
+//--
+{
+  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 *********/
+//++
+float Histo::Sigma(float x0, float dx) const
+//
+//	Valeur du sigma calcule entre x0-dx et x0+dx
+//--
+{
+  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 *********/
+//++
+float Histo::CleanedMean(float& sigma) const
+//
+//	Valeur de la moyenne et du sigma nettoyes
+//--
+{
+  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 *********/
+//++
+float Histo::CleanedMean() const
+//
+//	Valeur de la moyenne nettoyee
+//--
+{
+  if (!nHist) return 0;
+  float s=0;
+  return CleanedMean(s);
+}
+
+/********* Methode *********/
+//++
+int Histo::BinNonNul() const
+//
+//	Retourne le nombre de bins non-nul
+//--
+{
+int non=0;
+for (int i=0;i<bins;i++) if( data[i] != 0. ) non++;
+return non;
+}
+
+/********* Methode *********/
+//++
+int Histo::ErrNonNul() const
+//
+//	Retourne le nombre de bins ayant une erreur non-nulle
+//--
+{
+if(err2==NULL) return -1;
+int non=0;
+for (int i=0;i<bins;i++) if( err2[i] != 0. ) non++;
+return non;
+}
+
+/********* Methode *********/
+//++
+int Histo::BinPercent(float per) const
+//
+//	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).
+//--
+{
+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 *********/
+//++
+int Histo::BinPercent(float x,float per,int& imin,int& imax)
+//
+//	Renvoie les numeros de bins imin,imax (0=<i<bins) tels que:
+//|    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
+//--
+{
+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 *********/
+//++
+int Histo::BinPercent(float x,float per,float& xmin,float& xmax)
+//
+//	Idem precedent mais renvoie xmin et 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 *********/
+//++
+void Histo::HInteg(float norm)
+//
+//	Remplace l'histogramme par son integrale normalise a norm:
+//	si norm <= 0 : pas de normalisation, integration seule
+//--
+{
+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 *********/
+//++
+void Histo::HDeriv()
+//
+//	Remplace l'histogramme par sa derivee
+//--
+{
+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 *********/
+//++
+void Histo::HRebin(int nbinew)
+//
+//	Pour rebinner l'histogramme sur nbinew bins
+//--
+{
+  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 *********/
+//++
+int Histo::MaxiLocal(float& maxi,int& imax,float& maxn,int& imaxn)
+//
+//	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 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 *********/
+//++
+float Histo::FitMax(int degree, float frac, int debug) const
+//
+//	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
+//--
+{
+  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;
+  }
+
+  Vector xFit(nLowHigh);
+  Vector yFit(nLowHigh);
+  Vector e2Fit(nLowHigh);
+  Vector 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 *********/
+//++
+float Histo::FindWidth(float frac, int debug) const
+//
+//	Calcul de la largeur a frac pourcent du maximum
+//	autour du bin du maximum.
+//	L'histo est suppose etre remplit de valeurs positives
+//--
+{
+  float xmax = BinCenter( IMax() );
+  return FindWidth(xmax,frac,debug);
+}
+
+/********* Methode *********/
+//++
+float Histo::FindWidth(float xmax,float frac, int debug) const
+//
+//	Calcul de la largeur a frac pourcent de la valeur du bin xmax.
+//	L'histo est suppose etre remplit de valeurs positives
+//--
+{
+  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 *********/
+//++
+int Histo::EstimeMax(float& xm,int SzPav)
+//
+//	Cf suivant mais im est le bin du maximum de l'histo
+//--
+{
+int im = IMax();
+return EstimeMax(im,xm,SzPav);
+}
+
+/********* Methode *********/
+//++
+int Histo::EstimeMax(int& im,float& xm,int SzPav)
+//
+//	Determine l'abscisses du maximum donne par im
+//	en moyennant dans un pave SzPav  autour du maximum
+//|  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
+//--
+{
+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 *********/
+//++
+void Histo::EstimeWidthS(float frac,float& widthG,float& widthD)
+//
+//	Determine la largeur a frac% du maximum a gauche (widthG)
+//	et a droite (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;
+}
+
+}
+
+//////////////////////////////////////////////////////////
+//++
+int  Histo::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(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;
+  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();
+}
+
+//++
+Histo Histo::FitResidus(GeneralFit& gfit)
+//
+//	Retourne une classe contenant les residus du fit ``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"));
+Vector 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;
+}
+
+//++
+Histo Histo::FitFunction(GeneralFit& gfit)
+//
+//	Retourne une classe contenant la fonction du fit ``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"));
+Vector 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 *********/
+//++
+void Histo::PrintF(FILE * fp, int hdyn,float hmin, float hmax,int pflag,
+	      int il, int ih)
+//
+//	Impression de l'histogramme dans le fichier fp
+//|   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
+//--
+{
+
+ 
+ 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 *********/
+//++
+void Histo::Print(int hdyn,float hmin, float hmax,int pflag,
+	     int il, int ih)
+//
+//	Impression de l'histogramme sur stdout
+//--
+{
+ Histo::PrintF(stdout, hdyn, hmin, hmax, pflag, il, ih);
+}
+
+
+// Rappel des inlines functions pour commentaires
+//++
+// inline float          XMin() const
+//	Retourne l'abscisse minimum
+//--
+//++
+// inline float          XMax() const
+//	Retourne l'abscisse maximum
+//--
+//++
+// inline int            NBins() const
+//	Retourne le nombre de bins
+//--
+//++
+// inline float          BinWidth() const
+//	Retourne la largeur du bin
+//--
+//++
+// inline float*         Bins() const
+//	Retourne le pointeur sur le tableaux des contenus
+//--
+//++
+// inline float          operator()(int i) const
+//	Retourne le contenu du bin i
+//--
+//++
+// inline float&         operator()(int i)
+//	Remplit le contenu du bin i
+//--
+//++
+// inline float          Error(int i) const
+//	Retourne l'erreur du bin i
+//--
+//++
+// inline double         Error2(int i) const
+//	Retourne l'erreur au carre du bin i
+//--
+//++
+// inline double&        Error2(int i)
+//	Remplit l'erreur au carre du bin i
+//--
+//++
+// inline float          NData() const
+//	Retourne la somme ponderee
+//--
+//++
+// inline float          NEntries()
+//	Retourne le nombre d'entrees
+//--
+//++
+// inline float          NOver() const
+//	Retourne le nombre d'overflow
+//--
+//++
+// inline float          NUnder() const
+//	Retourne le nombre d'underflow
+//--
+//++
+// inline float          BinLowEdge(int i)  const
+//	Retourne l'abscisse du bord inferieur du bin i
+//--
+//++
+// inline float          BinCenter(int i)   const
+//	Retourne l'abscisse du centre du bin i
+//--
+//++
+// inline float          BinHighEdge(int i) const
+//	Retourne l'abscisse du bord superieur du bin i
+//--
+//++
+// inline int            FindBin(float x) const   
+//	Retourne le numero du bin contenant l'abscisse x
+//--
+
+
+///////////////////////////////////////////////////////////
+// --------------------------------------------------------
+//   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
Index: /trunk/SophyaLib/HiStats/histos.h
===================================================================
--- /trunk/SophyaLib/HiStats/histos.h	(revision 763)
+++ /trunk/SophyaLib/HiStats/histos.h	(revision 763)
@@ -0,0 +1,172 @@
+// This may look like C code, but it is really -*- C++ -*-
+//
+// $Id: histos.h,v 1.1.1.1 2000-03-02 16:51:32 ansari Exp $
+//
+
+#ifndef HISTOS_SEEN
+#define HISTOS_SEEN
+
+#include "objfio.h"
+#include <iostream.h>
+#include <stdio.h>
+#include "peida.h"
+#include "tvector.h"
+#include "ppersist.h"
+#include "anydataobj.h"
+
+namespace SOPHYA {
+
+class GeneralFit;
+
+class Histo : public AnyDataObj {
+  friend class ObjFileIO<Histo>;
+public:
+
+  // CREATOR / DESTRUCTOR
+  Histo();
+  Histo(float xMin, float xMax, int nBin=100);
+  Histo(const Histo& H);
+  virtual        ~Histo();
+
+  // OPTIONS
+  void   Errors();
+
+  // UPDATING or SETTING
+  void   Zero();
+  void   Add(float x, float w = 1.);
+  void   AddBin(int numBin, float w = 1.);
+  void   SetBin(float x, float w = 1.);
+  void   SetBin(int numBin, float w = 1.);
+  void   SetErr2(float x, double e2);
+  void   SetErr2(int numBin, double e2);
+  void   SetErr(float x, float e);
+  void   SetErr(int numBin, float e);
+
+  // Operators
+  Histo& operator = (const Histo& h);
+  Histo& operator *= (double b);
+  Histo& operator /= (double b);
+  Histo& operator += (double b);
+  Histo& operator -= (double b);
+  friend Histo operator * (const Histo& a, double b);
+  friend Histo operator * (double b, const Histo& a);
+  friend Histo operator / (const Histo& a, double b);
+  friend Histo operator + (const Histo& a, double b);
+  friend Histo operator + (double b, const Histo& a);
+  friend Histo operator - (const Histo& a, double b);
+  friend Histo operator - (double b, const Histo& a);
+  Histo& operator += (const Histo& a);
+  Histo& operator -= (const Histo& a);
+  Histo& operator *= (const Histo& a);
+  Histo& operator /= (const Histo& a);
+  friend Histo operator + (const Histo& a, const Histo& b);
+  friend Histo operator - (const Histo& a, const Histo& b);
+  friend Histo operator * (const Histo& a, const Histo& b);
+  friend Histo operator / (const Histo& a, const Histo& b);
+
+  // get/put dans/depuis un vector
+  void GetAbsc(Vector& v);
+  void GetValue(Vector& v);
+  void GetError2(Vector& v);
+  void GetError(Vector& v);
+  void PutValue(Vector& v, int ierr=0);
+  void PutValueAdd(Vector &v, int ierr=0);
+  void PutError2(Vector& v);
+  void PutError2Add(Vector& v);
+  void PutError(Vector& v);
+
+  // INLINES
+  inline float          XMin() const             {return min;}
+  inline float          XMax() const             {return max;}
+  inline int_4          NBins() const            {return bins;}
+  inline float          BinWidth() const         {return binWidth;}
+  inline float*         Bins() const             {return data;}
+  inline float          operator()(int i) const  {return data[i];}
+  inline float&         operator()(int i)        {return data[i];}
+  inline bool           HasErrors() { if(err2) return true; else return false;}
+  inline float          Error(int i) const
+                          {if(err2) {if(err2[i]>0.) return sqrt(err2[i]); else return 0.;}
+                           else return 0.;}
+  inline double         Error2(int i) const
+                          {if(err2) return err2[i]; else return 0.;}
+  inline double&        Error2(int i) {return err2[i];}
+  inline float          NData() const            {return (float) nHist;}
+  inline float          NEntries() const         {return nEntries;}
+  inline float          NOver() const            {return over;}
+  inline float          NUnder() const           {return under;}
+
+  inline float          BinLowEdge(int i)  const {return min + i*binWidth;}
+  inline float          BinCenter(int i)   const {return min + (i+0.5)*binWidth;}
+  inline float          BinHighEdge(int i) const {return min + (i+1)*binWidth;} 
+  inline int_4          FindBin(float x) const   
+                            {return (int_4) floorf((x - min) / binWidth);}
+
+  // Info, statistique et calculs sur les histogrammes
+  int            BinNonNul() const;
+  int            ErrNonNul() const;
+  int            IMax() const;
+  int            IMin() const;
+  float          VMax() const;
+  float          VMin() const;
+  float          Mean() const;
+  float          Sigma() const;
+  float          MeanLH(int il,int ih) const;
+  float          SigmaLH(int il,int ih) const;
+  float          Mean(float x0, float dx) const;
+  float          Sigma(float x0, float dx) const;
+  float          CleanedMean() const;
+  float          CleanedMean(float& sigma) const;
+  int            BinPercent(float per) const;
+  int            BinPercent(float x,float per,int& imin,int& imax);
+  int            BinPercent(float x,float per,float& xmin,float& xmax);
+  void           HInteg(float norm = 0.);
+  void           HDeriv();
+  void           HRebin(int nbinew);
+
+  int            MaxiLocal(float& maxi,int& imax,float& maxn,int& imaxn);
+  float          FitMax(int degree=2, float frac=0.5f, int debug=0) const;
+  float          FindWidth(float xmax,float frac=0.5f, int debug=0) const;
+  float          FindWidth(float frac=0.5f, int debug=0) const;
+  int            EstimeMax(float& xm,int SzPav = 3);
+  int            EstimeMax(int& im,float& xm,int SzPav = 3);
+  void           EstimeWidthS(float frac,float& widthG,float& widthD);
+
+  // Fit
+  int    Fit(GeneralFit& gfit,unsigned short typ_err=0);
+  Histo  FitResidus(GeneralFit& gfit);
+  Histo  FitFunction(GeneralFit& gfit);
+
+  // Print et Display ASCII
+  void           PrintF(FILE * fp, int dyn = 100, float hmin = 1., float hmax = -1.,
+			int pflag = 0, int il = 1, int ih = -1);
+  void           Print(int dyn = 100, float hmin = 1., float hmax = -1.,
+			int pflag = 0, int il = 1, int ih = -1);
+
+protected:
+  void Delete();
+
+  float*         data;		 
+  double*        err2;		 
+  float          under;
+  float          over;
+  double         nHist;
+  int_4          nEntries;
+  int_4          bins;
+  float          min;
+  float          max;
+  float          binWidth;
+};
+
+
+inline POutPersist& operator << (POutPersist& os, Histo & obj)
+{ ObjFileIO<Histo> fio(&obj);  fio.Write(os);  return(os); }
+inline PInPersist& operator >> (PInPersist& is, Histo & obj)
+{ ObjFileIO<Histo> fio(&obj);  fio.Read(is);  return(is); }
+
+// Classe pour la gestion de persistance
+// ObjFileIO<Histo>
+
+
+} // Fin du namespace
+
+#endif // HISTOS_SEEN
Index: /trunk/SophyaLib/HiStats/histos2.cc
===================================================================
--- /trunk/SophyaLib/HiStats/histos2.cc	(revision 763)
+++ /trunk/SophyaLib/HiStats/histos2.cc	(revision 763)
@@ -0,0 +1,2038 @@
+//
+//    cmv 05/08/96
+//
+
+#include "machdefs.h"
+
+#include <string.h>
+#include <new.h>
+#include <stdlib.h>
+#include <stdio.h>
+#ifndef __mac__
+#include <values.h>
+#endif
+
+#include "histos2.h"
+#include "generalfit.h"
+
+//++
+// Class	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)
+{
+ASSERT(nxBin>0 && nyBin>0 && xMin<xMax && yMin<yMax);
+for(int i=0;i<3;i++) for(int j=0;j<3;j++) over[i][j]=0.;
+Zero();
+b_s.H = NULL;
+END_CONSTRUCTOR
+}
+
+//++
+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
+}
+
+///////////////////////////////////////////////////////////////////
+//++
+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.NBinX()!=a.NBinX() || b.NBinY()!=a.NBinY()) THROW(sizeMismatchErr);
+Histo2D c(a);
+return (c += b);
+}
+
+//++
+Histo2D operator - (const Histo2D& a, const Histo2D& b)
+//
+//--
+{
+if (b.NBinX()!=a.NBinX() || b.NBinY()!=a.NBinY()) THROW(sizeMismatchErr);
+Histo2D c(a);
+return (c -= b);
+}
+
+//++
+Histo2D operator * (const Histo2D& a, const Histo2D& b)
+//
+//--
+{
+if (b.NBinX()!=a.NBinX() || b.NBinY()!=a.NBinY()) THROW(sizeMismatchErr);
+Histo2D c(a);
+return (c *= b);
+}
+
+//++
+Histo2D operator / (const Histo2D& a, const Histo2D& b)
+//
+//--
+{
+if (b.NBinX()!=a.NBinX() || b.NBinY()!=a.NBinY()) 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)
+  throw(SzMismatchError("Histo2D::FitResidus: size mismatch\n"));
+GeneralFunction* f = gfit.GetFunction();
+if(f==NULL)
+  throw(NullPtrError("Histo2D::FitResidus: NULL pointer\n"));
+Vector par = gfit.GetParm();
+Histo2D h2(*this);
+for(int i=0;i<NBinX();i++) for(int j=0;j<NBinY();j++) {
+  float xc,yc;
+  BinCenter(i,j,xc,yc);
+  double x[2] = {(double)xc,(double)yc};
+  h2(i,j) -= (float) f->Value(x,par.Data());
+}
+return h2;
+}
+
+//++
+Histo2D Histo2D::FitFunction(GeneralFit& gfit)
+//
+//	Retourne une classe contenant la fonction du fit ``gfit''.
+//--
+{
+if(NBinX()<=0 || NBinY()<=0)
+  throw(SzMismatchError("Histo2D::FitFunction: size mismatch\n"));
+GeneralFunction* f = gfit.GetFunction();
+if(f==NULL)
+  throw(NullPtrError("Histo2D::FitFunction: NULL pointer\n"));
+Vector par = gfit.GetParm();
+Histo2D h2(*this);
+for(int i=0;i<NBinX();i++) for(int j=0;j<NBinY();j++) {
+  float xc,yc;
+  BinCenter(i,j,xc,yc);
+  double x[2] = {(double)xc,(double)yc};
+  h2(i,j) = (float) f->Value(x,par.Data());
+}
+return h2;
+}
+
+///////////////////////////////////////////////////////////////////
+//++
+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
+//--
+
+
+///////////////////////////////////////////////////////////
+// --------------------------------------------------------
+//   Les objets delegues pour la gestion de persistance 
+// --------------------------------------------------------
+///////////////////////////////////////////////////////////
+
+
+void  ObjFileIO<Histo2D>::ReadSelf(PInPersist& is)
+{
+char strg[256];
+
+if(dobj==NULL) dobj=new Histo2D;
+  else         dobj->Delete();
+
+float min,max;
+int_4 errok, projx, projy, nslix, nsliy, nbanx, nbany;
+
+// Lecture entete
+is.GetLine(strg, 255);
+is.GetLine(strg, 255);
+is.GetLine(strg, 255);
+is.GetLine(strg, 255);
+is.GetLine(strg, 255);
+is.GetLine(strg, 255);
+
+// Lecture variables de definitions
+is.Get(dobj->nx);
+is.Get(dobj->ny);
+is.Get(dobj->nxy);
+is.Get(errok);
+is.Get(dobj->nEntries);
+is.Get(dobj->nHist);
+
+is.Get(dobj->xmin);
+is.Get(dobj->xmax);
+is.Get(dobj->ymin);
+is.Get(dobj->ymax);
+is.Get(dobj->wbinx);
+is.Get(dobj->wbiny);
+
+is.Get(&(dobj->over[0][0]),9);
+
+is.Get(projx);
+is.Get(projy);
+is.Get(nslix);
+is.Get(nsliy);
+is.Get(nbanx);
+is.Get(nbany);
+
+// Lecture histo2D
+dobj->data = new float[dobj->nxy];
+is.GetLine(strg, 255);
+{for(int j=0;j<dobj->ny;j++) is.Get(dobj->data+j*dobj->nx,dobj->nx);}
+
+// Lecture erreurs
+if(errok) {
+  is.GetLine(strg, 255);
+  dobj->err2 = new double[dobj->nxy];
+  for(int j=0;j<dobj->ny;j++) is.Get(dobj->err2+j*dobj->nx,dobj->nx);
+}
+
+// Lecture des projections
+if(projx) {
+  is.GetLine(strg, 255);
+  dobj->SetProjX();
+  ObjFileIO<Histo> fio_h(dobj->hprojx);
+  fio_h.Read(is);
+}
+if(projy) {
+  is.GetLine(strg, 255);
+  dobj->SetProjY();
+  ObjFileIO<Histo> fio_h(dobj->hprojy);
+  fio_h.Read(is);
+}
+
+// Lecture des slices
+if(nslix>0) {
+  is.GetLine(strg, 255);
+  dobj->SetSliX(nslix);
+  ASSERT (nslix==dobj->NSliX());
+  for(int j=0;j<dobj->NSliX();j++)
+    {ObjFileIO<Histo> fio_h(dobj->HSliX(j)); fio_h.Read(is);}
+}
+if(nsliy>0) {
+  is.GetLine(strg, 255);
+  dobj->SetSliY(nsliy);
+  ASSERT (nsliy==dobj->NSliY());
+  for(int j=0;j<dobj->NSliY();j++)
+    {ObjFileIO<Histo> fio_h(dobj->HSliY(j)); fio_h.Read(is);}
+}
+
+// Lecture des bandes
+if( nbanx>0 ) {
+  is.GetLine(strg, 255);
+  {for(int j=0; j<nbanx; j++) {
+    is.Get(min); is.Get(max);
+    dobj->SetBandX(min,max);
+  }}
+  ASSERT (nbanx==dobj->NBandX());
+  {for(int j=0; j<dobj->NBandX(); j++) {
+    ObjFileIO<Histo> fio_h(dobj->HBandX(j));
+    fio_h.Read(is);
+  }}
+}
+if( nbany>0 ) {
+  is.GetLine(strg, 255);
+  {for(int j=0; j<nbany; j++) {
+    is.Get(min); is.Get(max);
+    dobj->SetBandY(min,max);
+  }}
+  ASSERT (nbany==dobj->NBandY());
+  {for(int j=0; j<dobj->NBandY(); j++) {
+    ObjFileIO<Histo> fio_h(dobj->HBandY(j));
+    fio_h.Read(is);
+  }}
+}
+
+return;
+}
+
+void ObjFileIO<Histo2D>::WriteSelf(POutPersist& os) const
+{
+if (dobj == NULL)   return;
+char strg[256];
+
+// Que faut-il ecrire?
+int_4 errok = (dobj->err2) ? 1 : 0;
+int_4 projx = (dobj->hprojx) ? 1 : 0;
+int_4 projy = (dobj->hprojy) ? 1 : 0;
+int_4 nslix = dobj->NSliX();
+int_4 nsliy = dobj->NSliY();
+int_4 nbanx = dobj->NBandX();
+int_4 nbany = dobj->NBandY();
+
+// Ecriture entete pour identifier facilement
+sprintf(strg,"nx=%d  ny=%d  nxy=%d errok=%1d"
+       ,dobj->nx,dobj->ny,dobj->nxy,errok);
+os.PutLine(strg);
+sprintf(strg,"nHist=%g nEntries=%d",dobj->nHist,dobj->nEntries);
+os.PutLine(strg);
+sprintf(strg,"wbinx=%g wbiny=%g",dobj->wbinx,dobj->wbiny);
+os.PutLine(strg);
+sprintf(strg,"xmin=%g xmax=%g ymin=%g ymax=%g"
+       ,dobj->xmin,dobj->xmax,dobj->ymin,dobj->ymax);
+os.PutLine(strg);
+sprintf(strg,"projx/y=%d %d nbandx/y=%d %d nbslix/y=%d %d"
+       ,projx,projy,nbanx,nbany,nslix,nsliy);
+os.PutLine(strg);
+sprintf(strg,"over %g %g %g %g %g %g %g %g %g"
+       ,dobj->over[0][0],dobj->over[0][1],dobj->over[0][2]
+       ,dobj->over[1][0],dobj->over[1][1],dobj->over[1][2]
+       ,dobj->over[2][0],dobj->over[2][1],dobj->over[2][2]);
+os.PutLine(strg);
+
+// Ecriture variables de definitions
+os.Put(dobj->nx);
+os.Put(dobj->ny);
+os.Put(dobj->nxy);
+os.Put(errok);
+os.Put(dobj->nEntries);
+os.Put(dobj->nHist);
+
+os.Put(dobj->xmin);
+os.Put(dobj->xmax);
+os.Put(dobj->ymin);
+os.Put(dobj->ymax);
+os.Put(dobj->wbinx);
+os.Put(dobj->wbiny);
+
+os.Put(&(dobj->over[0][0]),9);
+
+os.Put(projx);
+os.Put(projy);
+os.Put(nslix);
+os.Put(nsliy);
+os.Put(nbanx);
+os.Put(nbany);
+
+// Ecriture histo2D
+sprintf(strg,"Histo2D: Tableau des donnees %d = %d * %d"
+       ,dobj->nxy,dobj->nx,dobj->ny);
+os.PutLine(strg);
+{for(int j=0;j<dobj->ny;j++) os.Put(dobj->data+j*dobj->nx,dobj->nx);}
+
+// Ecriture erreurs
+if(errok) {
+  sprintf(strg,"Histo2D: Tableau des erreurs %d = %d * %d"
+         ,dobj->nxy,dobj->nx,dobj->ny);
+  os.PutLine(strg);
+  for(int j=0;j<dobj->ny;j++) os.Put(dobj->err2+j*dobj->nx,dobj->nx);
+}
+
+// Ecriture des projections
+if(projx) {
+  sprintf(strg,"Histo2D: Projection X");
+  os.PutLine(strg);
+  ObjFileIO<Histo> fio_h(dobj->hprojx); fio_h.Write(os);
+}
+if(projy) {
+  sprintf(strg,"Histo2D: Projection Y");
+  os.PutLine(strg);
+  ObjFileIO<Histo> fio_h(dobj->hprojy); fio_h.Write(os);
+}
+
+// Ecriture des slices
+if(nslix>0) {
+  sprintf(strg,"Histo2D: Slices X %d",nslix);
+  os.PutLine(strg);
+  for(int j=0;j<nslix;j++) {
+    Histo* h = dobj->HSliX(j);
+    ObjFileIO<Histo> fio_h(h); fio_h.Write(os);
+  }
+}
+if(nsliy>0) {
+  sprintf(strg,"Histo2D: Slices Y %d",nsliy);
+  os.PutLine(strg);
+  for(int j=0;j<nsliy;j++) {
+    Histo* h = dobj->HSliY(j);
+    ObjFileIO<Histo> fio_h(h); fio_h.Write(os);
+  }
+}
+
+// Ecriture des bandes
+if( nbanx>0 ) {
+  sprintf(strg,"Histo2D: Bandes X %d",nbanx);
+  os.PutLine(strg);
+  list<Histo2D::bande_slice>::const_iterator it;
+  for(it = dobj->lbandx.begin(); it != dobj->lbandx.end(); it++) {
+    float min = (*it).min; float max = (*it).max;
+    os.Put(min); os.Put(max);
+  }
+  for(it = dobj->lbandx.begin(); it != dobj->lbandx.end(); it++) {
+    Histo* h = (*it).H;
+    ObjFileIO<Histo> fio_h(h); fio_h.Write(os);
+  }
+}
+if( nbany>0 ) {
+  sprintf(strg,"Histo2D: Bandes Y %d",nbany);
+  os.PutLine(strg);
+  list<Histo2D::bande_slice>::const_iterator it;
+  for(it = dobj->lbandy.begin(); it != dobj->lbandy.end(); it++) {
+    float min = (*it).min; float max = (*it).max;
+    os.Put(min); os.Put(max);
+  }
+  for(it = dobj->lbandy.begin(); it != dobj->lbandy.end(); it++) {
+    Histo* h = (*it).H;
+    ObjFileIO<Histo> fio_h(h); fio_h.Write(os);
+  }
+}
+
+return;
+}
+
+
+#ifdef __CXX_PRAGMA_TEMPLATES__
+#pragma define_template ObjFileIO<Histo2D>
+#endif
+
+#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
+template class ObjFileIO<Histo2D>;
+#endif
Index: /trunk/SophyaLib/HiStats/histos2.h
===================================================================
--- /trunk/SophyaLib/HiStats/histos2.h	(revision 763)
+++ /trunk/SophyaLib/HiStats/histos2.h	(revision 763)
@@ -0,0 +1,229 @@
+//
+// histogrammes 2D     cmv 2/8/96
+//
+
+#ifndef HISTOS2_SEEN
+#define HISTOS2_SEEN
+
+#include "machdefs.h"
+#include <string>
+
+#include <list>
+#if defined(__KCC__)
+#include <list.h>
+#endif
+
+#include "peida.h"
+#include "utils.h"
+#include "histos.h"
+
+namespace SOPHYA {
+
+class GeneralFit;
+
+class Histo2D : public AnyDataObj {
+  friend class ObjFileIO<Histo2D>;
+public:
+
+  // CREATOR / DESTRUCTOR
+  Histo2D(float xMin, float xMax, int nxBin
+         ,float yMin, float yMax, int nyBin);
+  Histo2D(const Histo2D& h);
+  Histo2D();
+  virtual  ~Histo2D();
+
+  // OPTIONS
+  void   Errors();  
+  
+  // UPDATING
+  void   Zero();
+  void   Add(float x, float y, float w = 1.);
+
+  // Operators
+  Histo2D& operator = (const Histo2D& h);
+  Histo2D& operator *= (double b);
+  Histo2D& operator /= (double b);
+  Histo2D& operator += (double b);
+  Histo2D& operator -= (double b);
+  friend Histo2D operator * (const Histo2D& a, double b);
+  friend Histo2D operator * (double b, const Histo2D& a);
+  friend Histo2D operator / (const Histo2D& a, double b);
+  friend Histo2D operator + (const Histo2D& a, double b);
+  friend Histo2D operator + (double b, const Histo2D& a);
+  friend Histo2D operator - (const Histo2D& a, double b);
+  friend Histo2D operator - (double b, const Histo2D& a);
+  Histo2D& operator += (const Histo2D& a);
+  Histo2D& operator -= (const Histo2D& a);
+  Histo2D& operator *= (const Histo2D& a);
+  Histo2D& operator /= (const Histo2D& a);
+  friend Histo2D operator + (const Histo2D& a, const Histo2D& b);
+  friend Histo2D operator - (const Histo2D& a, const Histo2D& b);
+  friend Histo2D operator * (const Histo2D& a, const Histo2D& b);
+  friend Histo2D operator / (const Histo2D& a, const Histo2D& b);
+
+  // get/put dans/depuis une matrice / vector
+  void GetXCoor(Vector& v);
+  void GetValue(Matrix &v);
+  void GetYCoor(Vector& v);
+  void GetError2(Matrix& v);
+  void GetError(Matrix& v);
+  void PutValue(Matrix& v, int ierr=0);
+  void PutValueAdd(Matrix& v, int ierr=0);
+  void PutError2(Matrix& v);
+  void PutError2Add(Matrix& v);
+  void PutError(Matrix& v);
+
+  // INLINES
+  inline float   XMin() const {return xmin;}
+  inline float   XMax() const {return xmax;}
+  inline float   YMin() const {return ymin;}
+  inline float   YMax() const {return ymax;}
+  inline int_4   NBinX() const {return nx;}
+  inline int_4   NBinY() const {return ny;}
+  inline float   WBinX() const {return wbinx;}
+  inline float   WBinY() const {return wbiny;}
+  inline float*  Bins() const  {return data;}
+  inline float   operator()(int i,int j) const {return data[j*nx+i];}
+  inline float&  operator()(int i,int j)       {return data[j*nx+i];}
+  inline bool    HasErrors() { if(err2) return true; else return false;}
+  inline float   Error(int i,int j)  const
+                   {if(err2)
+                      {if(err2[j*nx+i]>0.) return sqrt(err2[j*nx+i]); else return 0.;}
+                    else return 0.;}
+  inline double  Error2(int i,int j) const
+                   {if(err2) return err2[j*nx+i]; else return 0.;}
+  inline double& Error2(int i,int j) {return err2[j*nx+i];}
+  inline float   NData() const    {return (float) nHist;}
+  inline int_4   NEntries() const {return nEntries;}
+  inline void BinLowEdge(int i,int j,float& x,float& y)
+                 {x = xmin + i*wbinx; y = ymin + j*wbiny;}
+  inline void BinCenter(int i,int j,float& x,float& y)
+                 {x = xmin + (i+0.5)*wbinx; y = ymin + (j+0.5)*wbiny;}
+  inline void BinHighEdge(int i,int j,float& x,float& y)
+                 {x = xmin + (i+1)*wbinx; y = ymin + (j+1)*wbiny;}
+  inline void FindBin(float x,float y,int& i,int& j)
+    { i = (int) floorf((x - xmin)/wbinx); j = (int) floorf((y - ymin)/wbiny);}
+
+  // Info, statistique et calculs sur les histogrammes
+  float   NOver(int i=-1,int j=-1) const;
+  int     BinNonNul() const;
+  int     ErrNonNul() const;
+  void    IJMax(int& imax,int& jmax,int il=1,int ih= -1,int jl=1,int jh= -1);
+  void    IJMin(int& imax,int& jmax,int il=1,int ih= -1,int jl=1,int jh= -1);
+  float   VMax(int il=1,int ih= -1,int jl=1,int jh= -1) const;
+  float   VMin(int il=1,int ih= -1,int jl=1,int jh= -1) const;
+  int     EstimeMax(float& xm,float& ym,int SzPav = 3
+                   ,int il=1,int ih= -1,int jl=1,int jh= -1);
+  int     EstimeMax(int im,int jm,float& xm,float& ym,int SzPav = 3);
+  int     FindMax(int& im,int& jm,int SzPav = 3,float Dz = 0.
+                 ,int il=1,int ih= -1,int jl=1,int jh= -1);
+
+  // Fit
+  int      Fit(GeneralFit& gfit,unsigned short typ_err=0);
+  Histo2D  FitResidus(GeneralFit& gfit);
+  Histo2D  FitFunction(GeneralFit& gfit);
+
+  // Print et Display ASCII
+  void    PrintStatus();
+  void    Print(float min=1.,float max=-1.
+               ,int il=1,int ih= -1,int jl=1,int jh= -1);
+
+  // PROJECTIONS
+  void           SetProjX();
+  void           SetProjY();
+  void           SetProj();
+  void           DelProjX();
+  void           DelProjY();
+  void           DelProj();
+  void           ZeroProjX();
+  void           ZeroProjY();
+  void           ZeroProj();
+  inline Histo*  HProjX() const {return hprojx;}
+  inline Histo*  HProjY() const {return hprojy;}
+  void           ShowProj();
+
+  // BANDES
+  inline int     NBandX() const {return lbandx.size();}
+  inline int     NBandY() const {return lbandy.size();}
+  int            SetBandX(float ybmin,float ybmax);
+  int            SetBandY(float xbmin,float xbmax);
+  void           DelBandX();
+  void           DelBandY();
+  void           ZeroBandX();
+  void           ZeroBandY();
+  Histo*         HBandX(int n) const;
+  Histo*         HBandY(int n) const;
+  void           GetBandX(int n,float& ybmin,float& ybmax) const;
+  void           GetBandY(int n,float& xbmin,float& xbmax) const;
+  void           ShowBand(int lp = 0);
+
+  // SLICES
+  inline int     NSliX() const {return lslix.size();}
+  inline int     NSliY() const {return lsliy.size();}
+  int            SetSliX(int nsli);
+  int            SetSliY(int nsli);
+  void           DelSliX();
+  void           DelSliY();
+  void           ZeroSliX();
+  void           ZeroSliY();
+  Histo*         HSliX(int n) const;
+  Histo*         HSliY(int n) const;
+  void           ShowSli(int lp = 0);
+
+#ifndef __DECCXX
+protected:
+#endif
+  struct bande_slice {
+    int num;
+    float min, max;
+    Histo* H;
+    STRUCTCOMP(bande_slice)
+  };
+#ifdef __DECCXX
+protected:
+#endif
+
+  void Delete();
+
+  float*         data;
+  double*        err2;
+
+  float          over[3][3];
+  double         nHist;
+  int_4          nEntries;
+
+  int_4          nx,ny,nxy;
+  float          xmin;
+  float          xmax;
+  float          ymin;
+  float          ymax;
+  float          wbinx;
+  float          wbiny;
+
+  bande_slice    b_s;
+
+  Histo*         hprojx;
+  Histo*         hprojy;
+
+  list<bande_slice>  lbandx;
+  list<bande_slice>  lbandy;
+  
+  list<bande_slice>  lslix;
+  list<bande_slice>  lsliy;
+
+};
+
+/////////////////////////////////////////////////////////////////////////
+// Classe pour la gestion de persistance
+
+inline POutPersist& operator << (POutPersist& os, Histo2D & obj)
+{ ObjFileIO<Histo2D> fio(&obj);  fio.Write(os);  return(os); }
+inline PInPersist& operator >> (PInPersist& is, Histo2D & obj)
+{ ObjFileIO<Histo2D> fio(&obj);  fio.Read(is);  return(is); }
+
+// Classe pour la gestion de persistance
+// ObjFileIO<Histo2D>
+
+} // Fin du namespace
+
+#endif // HISTOS2_SEEN
Index: /trunk/SophyaLib/HiStats/ntuple.cc
===================================================================
--- /trunk/SophyaLib/HiStats/ntuple.cc	(revision 763)
+++ /trunk/SophyaLib/HiStats/ntuple.cc	(revision 763)
@@ -0,0 +1,611 @@
+#include <stdio.h>
+#include <string.h>
+
+#include "strutil.h"
+#include "perrors.h"
+#include "ntuple.h"
+
+
+#define BADVAL -1.e19
+#define LENNAME 8
+#define LENNAME1  (LENNAME+1)
+
+//++
+// Class	NTuple
+// Lib	Outils++ 
+// include	ntuple.h
+//
+//	Classe de ntuples
+//--
+//++
+// Links        Parents
+// PPersist
+// NTupleInterface
+//--
+
+/* --Methode-- */
+//++
+NTuple::NTuple()
+//
+//	Createur par defaut
+//--
+{
+mNVar = mNEnt = mBlk = mNBlk = 0;
+mVar = NULL;
+mVarD = NULL;
+mNames = NULL;
+mInfo = NULL;
+}
+
+
+//++
+NTuple::NTuple(int nvar, char** noms, int blk)
+//
+//	Createur d'un ntuple de `nvar' variables dont les
+//	noms sont dans le tableau de cahines de caracteres `noms'
+//	avec `blk' d'evenements par blocks.
+//--
+{
+mNVar = mNEnt = mBlk = mNBlk = 0;
+mVar = NULL;
+mVarD = NULL;
+mNames = NULL;
+mInfo = NULL;
+if (nvar <= 0)  THROW(sizeMismatchErr); 
+mNVar = nvar;
+mVar = new r_4[nvar];
+mVarD = new r_8[nvar];
+if (blk < 10) blk = 10;
+mBlk = blk;
+// On prend des noms de LENNAME char pour le moment
+mNames = new char[nvar*LENNAME1];
+r_4* pt = new r_4[nvar*blk];
+mNBlk = 1;
+mPtr.push_back(pt);
+int i;
+for(i=0; i<nvar; i++)
+  { strncpy(mNames+i*LENNAME1, noms[i], LENNAME);  
+  mNames[i*LENNAME1+LENNAME] = '\0'; }
+return;
+}
+
+//                                       cmv 8/10/99
+//++
+NTuple::NTuple(const NTuple& NT)
+//
+//	Createur par copie (clone).
+//--
+: mNVar(0), mNEnt(0), mBlk(0), mNBlk(0)
+, mVar(NULL), mVarD(NULL), mNames(NULL), mInfo(NULL)
+{
+if(NT.mNVar<=0) return; // cas ou NT est cree par defaut
+mNVar = NT.mNVar;
+mBlk = NT.mBlk;
+mVar = new r_4[NT.mNVar];
+mVarD = new r_8[NT.mNVar];
+mNames = new char[NT.mNVar*LENNAME1];
+
+int i;
+r_4* pt = new r_4[mNVar*mBlk];
+mNBlk = 1; mPtr.push_back(pt);
+
+for(i=0;i<mNVar;i++) strcpy(mNames+i*LENNAME1,NT.NomIndex(i)); 
+
+if(NT.mInfo!=NULL) {mInfo = new DVList; *mInfo = *(NT.mInfo);}
+
+if(NT.mNEnt<=0) return;
+for(i=0;i<NT.mNEnt;i++) {r_4* x=NT.GetVec(i,NULL); Fill(x);}
+
+return;
+}
+
+/* --Methode-- */
+//++
+NTuple::NTuple(char *flnm)
+//
+//	Createur lecture fichier ppersist.
+//--
+{
+mNVar = mNEnt = mBlk = mNBlk = 0;
+mVar = NULL;
+mVarD = NULL;
+mNames = NULL;
+mInfo = NULL;
+PInPersist s(flnm);
+ObjFileIO<NTuple> fiont(this);
+fiont.Read(s);
+}
+
+/* --Methode-- */
+NTuple::~NTuple()
+{
+Clean();
+}
+
+/* --Methode-- */
+void NTuple::Clean()
+{
+if (mVar) delete[] mVar;
+if (mVarD) delete[] mVarD;
+if (mNames)  delete[] mNames;
+if (mInfo) delete mInfo;
+int i;
+if(mNBlk>0) for(i=0; i<mNBlk; i++)  delete[] mPtr[i];
+mPtr.erase(mPtr.begin(), mPtr.end());
+mNVar = mNEnt = mBlk = mNBlk = 0;
+mVar = NULL;
+mVarD = NULL;
+mNames = NULL;
+mInfo = NULL;
+return;
+}
+
+/* --Methode--        cmv 08/10/99 */
+//++
+NTuple& NTuple::operator = (const NTuple& NT)
+//
+//	Operateur egal (clone).
+//--
+{
+if(this == &NT) return *this;
+Clean();
+if(NT.mNVar<=0) return *this; // cas ou NT est cree par defaut
+mNVar = NT.mNVar;
+mBlk = NT.mBlk;
+mVar = new r_4[NT.mNVar];
+mVarD = new r_8[NT.mNVar];
+mNames = new char[NT.mNVar*LENNAME1];
+
+int i;
+r_4* pt = new r_4[mNVar*mBlk];
+mNBlk = 1; mPtr.push_back(pt);
+
+for(i=0;i<mNVar;i++) strcpy(mNames+i*LENNAME1,NT.NomIndex(i)); 
+
+if(NT.mInfo!=NULL) {mInfo = new DVList; *mInfo = *(NT.mInfo);}
+
+if(NT.mNEnt<=0) return *this;
+for(i=0;i<NT.mNEnt;i++) {r_4* x=NT.GetVec(i,NULL); Fill(x);}
+
+// En fait il faudrait un createur par copie qui partage les donnees
+// quand l'objet est temporaire... trop complique A FAIRE !  cmv.
+return *this;
+}
+
+/* --Methode-- */
+//++
+void  NTuple::Fill(r_4* x)
+//
+//	Remplit le ntuple avec le tableau cd reels `x'.
+//--
+{
+int numb = mNEnt/mBlk;
+if (numb >= mNBlk) {
+  r_4* pt = new r_4[mNVar*mBlk];
+  mNBlk++;
+  mPtr.push_back(pt);
+}
+int offb = mNEnt-numb*mBlk;
+memcpy((mPtr[numb]+offb*mNVar), x, mNVar*sizeof(r_4));
+mNEnt++;
+return;
+}
+
+
+/* --Methode-- */
+//++
+float NTuple::GetVal(int n, int k)  const
+//
+//	Retourne la valeur de la variable `k' de l'evenement `n'.
+//--
+{
+if (n >= mNEnt)   return(BADVAL);
+if ( (k < 0) || (k >= mNVar) )    return(BADVAL);
+int numb = n/mBlk;
+int offb = n-numb*mBlk;
+return(*(mPtr[numb]+offb*mNVar+k));
+}
+
+
+/* --Methode-- */
+//++
+int NTuple::IndexNom(const char* nom)  const
+//
+//	Retourne le numero de la variable de nom `nom'.
+//--
+{
+int i;
+for(i=0; i<mNVar; i++)  
+  if ( strcmp(nom, mNames+i*LENNAME1) == 0)  return(i);
+return(-1);
+}
+
+
+static char nomretour[2*LENNAME1];
+/* --Methode-- */
+//++
+char* NTuple::NomIndex(int k)  const
+//
+//	Retourne le nom de la variable numero 'k'
+//--
+{
+nomretour[0] = '\0';
+if ((k >= 0) && (k < mNVar))  strcpy(nomretour, mNames+k*LENNAME1);
+return(nomretour);
+}
+
+  
+/* --Methode-- */
+//++
+r_4* NTuple::GetVec(int n, r_4* ret)   const
+//
+//	Retourne l'evenement `n' dans le vecteur `ret'.
+//--
+{
+int i;
+if (ret == NULL)   ret = mVar;
+if (n >= mNEnt) {
+  for(i=0; i<mNVar; i++)   ret[i] = BADVAL;
+  return(ret);
+}
+  
+int numb = n/mBlk;
+int offb = n-numb*mBlk;
+memcpy(ret, (mPtr[numb]+offb*mNVar), mNVar*sizeof(r_4));
+return(ret);
+}
+
+/* --Methode-- */
+//++
+r_8* NTuple::GetVecD(int n, r_8* ret)   const
+//
+//	Retourne l'evenement `n' dans le vecteur `ret'.
+//--
+{
+int i;
+if (ret == NULL)   ret = mVarD;
+float *  fr = GetVec(n);
+for(i=0; i<mNVar; i++)   ret[i] = fr[i];
+return(ret);
+}
+
+
+
+/* --Methode-- */
+//++
+DVList&  NTuple::Info()
+//
+//	Renvoie une référence sur l'objet DVList Associé
+//--
+{
+if (mInfo == NULL)  mInfo = new DVList;
+return(*mInfo);
+}
+
+/* --Methode-- */
+//++
+void  NTuple::Print(int num, int nmax)  const
+//
+//	Imprime `nmax' evenements a partir du numero `num'.
+//--
+{
+int i,j;
+
+printf("Num     ");
+for(i=0; i<mNVar; i++)  printf("%8s ", mNames+i*LENNAME1);
+putchar('\n');
+
+if (nmax <= 0)  nmax = 1;
+if (num < 0)  num = 0;
+nmax += num;
+if (nmax > mNEnt) nmax = mNEnt;
+for(i=num; i<nmax; i++) {
+  GetVec(i, NULL); 
+  printf("%6d  ", i);  
+  for(j=0; j<mNVar; j++)  printf("%8g ", (float)mVar[j]);
+  putchar('\n');
+}
+return;
+}
+
+/* --Methode-- */
+//++
+void  NTuple::Show(ostream& os)  const
+//
+//	Imprime l'information generale sur le ntuple.
+//--
+{
+os << "NTuple: NVar= " << mNVar << " NEnt=" << mNEnt  
+   << " (Blk Sz,Nb= " << mBlk << " ," << mNBlk << ")\n";
+os << "            Variables       Min      Max       \n";
+int i;
+double min, max;
+char buff[128];
+for(i=0; i<mNVar; i++) {
+  GetMinMax(i, min, max);
+  sprintf(buff, "%3d  %16s  %10lg  %10lg \n", i, mNames+i*LENNAME1, min, max);
+  os << (string)buff ;
+  }
+os << endl;
+}
+
+
+/* --Methode-- */
+//++
+int  NTuple::FillFromASCIIFile(string const& fn, float defval)
+//
+//	Remplit le ntuple a partir d'un fichier ASCII.
+//	Renvoie le nombre de lignes ajoutees.
+//--
+{
+if (NbColumns() < 1)  { 
+  cout << "NTuple::FillFromASCIIFile() Ntuple has " << NbColumns() << " columns" << endl;
+  return(-1);
+  }
+FILE * fip = fopen(fn.c_str(), "r");
+if (fip == NULL) {
+  cout << "NTuple::FillFromASCIIFile() Error opening file " << fn << endl;
+  return(-2);
+  }
+
+char lineb[4096];
+char *line;
+char buff[64];
+char* ccp;
+int i,j,l,kk;
+int postab, posb;
+float* xv = new float[NbColumns()];
+
+int nlread = 0;
+int nvar = NbColumns();
+// On boucle sur toutes les lignes
+while (fgets(lineb,4096,fip) != NULL) {
+  lineb[4095] = '\0';
+  j = 0; line = lineb;
+//  On enleve les espaces et tab de debut 
+  while ( (line[j] != '\0') && ((line[j] == ' ') || (line[j] == '\t')) )  j++;
+  line = lineb+j; 
+// Il faut que le premier caractere non-espace soit un digit, ou + ou - ou .
+  if (!( isdigit(line[0]) || (line[0] == '+') || (line[0] == '-') || (line[0] == '.') ))  continue;
+  ccp = line;
+  for(kk=0; kk<nvar; kk++)  xv[kk] = defval;
+  for(kk=0; kk<nvar; kk++) {
+// Les mots sont separes par des espaces ou des tab
+    postab = posc(ccp, '\t' );
+    posb = posc(ccp, ' ' );
+    if (postab >= 0) { 
+       if (posb < 0) posb = postab;
+       else if (postab < posb)  posb = postab; 
+       } 
+    if (posb >= 0)  ccp[posb] = '\0';
+    if ( isdigit(line[0]) || (line[0] == '+') || (line[0] == '-') || (line[0] == '.') ) 
+      xv[kk] = atof(ccp);
+    if (posb < 0)  break;
+    ccp += posb+1;   j = 0;
+    while ( (ccp[j] != '\0') && ((ccp[j] == ' ') || (ccp[j] == '\t')) )  j++;
+    ccp += j;
+    }
+  Fill(xv);
+  nlread++;
+  }
+
+delete[] xv;
+cout << "NTuple::FillFromASCIIFile( " << fn << ") " << nlread << " fill from file " << endl;
+return(nlread);
+}
+
+
+// ------- Implementation de  l interface NTuple  ---------
+
+/* --Methode-- */
+uint_4 NTuple::NbLines() const
+{
+return(NEntry());
+}
+/* --Methode-- */
+uint_4 NTuple::NbColumns() const
+{
+return(NVar());
+}
+
+/* --Methode-- */
+r_8 * NTuple::GetLineD(int n) const
+{
+return(GetVecD(n));
+}
+
+/* --Methode-- */
+r_8 NTuple::GetCell(int n, int k) const
+{
+return(GetVal(n, k));
+}
+
+/* --Methode-- */
+r_8 NTuple::GetCell(int n, string const & nom) const
+{
+return(GetVal(n, nom.c_str()));
+}
+
+/* --Methode-- */
+//++
+void  NTuple::GetMinMax(int k, double& min, double& max)  const
+//
+//	Retourne le minimum et le maximum de la variable `k'.
+//--
+{
+min = 9.e19; max = -9.e19;
+if ( (k < 0) || (k >= mNVar) )    return;
+int jb,ib,i;
+double x;
+i=0;
+for(jb=0; jb< mNBlk; jb++)
+  for(ib=0; ib< mBlk; ib++) {
+    if (i >= mNEnt)  break;
+    i++;
+    x = *(mPtr[jb]+ib*mNVar+k);
+    if(i==1) {min = x; max = x;}
+    if (x < min)  min = x;
+    if (x > max)  max = x;
+  }
+return;
+}
+
+/* --Methode-- */
+void NTuple::GetMinMax(string const & nom, double& min, double& max)   const
+{
+GetMinMax(IndexNom(nom.c_str()), min, max);
+}
+
+/* --Methode-- */
+int NTuple::ColumnIndex(string const & nom)  const
+{
+return(IndexNom(nom.c_str()));
+}
+
+/* --Methode-- */
+string NTuple::ColumnName(int k) const
+{
+return(NomIndex(k));
+}
+
+/* --Methode-- */
+//++
+string NTuple::VarList_C(const char* nomx)  const
+//
+//	Retourne une chaine de caracteres avec la declaration des noms de 
+//	variables. si "nomx!=NULL" , des instructions d'affectation
+//	a partir d'un tableau "nomx[i]" sont ajoutees. 
+//--
+{
+string rets="";
+int i;
+for(i=0; i<mNVar; i++) {
+  if ( (i%5 == 0) && (i > 0) )  rets += ";";  
+  if (i%5 == 0)   rets += "\ndouble "; 
+  else rets += ",";
+  rets += mNames+i*LENNAME1;
+  }
+rets += "; \n";
+if (nomx) { 
+  char buff[256];
+  for(i=0; i<mNVar; i++) {
+    sprintf(buff,"%s=%s[%d]; ",  mNames+i*LENNAME1, nomx, i);
+    rets += buff;
+    if ( (i%3 == 0) && (i > 0) )  rets += "\n"; 
+    }
+  }
+
+return(rets);
+}
+
+
+/* --Methode-- */
+//++
+string NTuple::LineHeaderToString() const 
+//	Retourne une chaine de caracteres avec la liste des noms de  
+//	variables, utilisables pour une impression
+//--
+{
+char buff[32];
+string rets=" Num    ";
+for(int i=0; i<mNVar; i++) {
+  sprintf(buff, "%8s ", mNames+i*LENNAME1);
+  rets += buff;
+  }
+rets += '\n';
+return(rets);
+}
+
+/* --Methode-- */
+//++
+string NTuple::LineToString(int n) const
+//	Retourne une chaine de caracteres avec le contenu de la ligne "n"
+//	utilisable pour une impression
+//--
+{
+char buff[32];
+double* val;
+val = GetLineD(n);
+sprintf(buff,"%6d: ",n);  
+string rets=buff;
+int i;
+for(i=0; i<mNVar; i++) {
+  sprintf(buff, "%8.3g ", val[i]);
+  rets += buff;
+  }
+rets += '\n';
+return(rets);
+}
+
+
+/* --Methode-- */
+//++
+void   ObjFileIO<NTuple>::WriteSelf(POutPersist& s)  const
+//
+//	Ecriture ppersist du ntuple.
+//--
+{
+char strg[256];
+if (dobj->mInfo)  sprintf(strg, "NVar=%6d  NEnt=%9d  BlkSz=%6d NBlk=%6d  HasInfo", 
+                          (int)dobj->mNVar, (int)dobj->mNEnt, (int)dobj->mBlk, (int)dobj->mNBlk);
+else sprintf(strg, "NVar=%6d  NEnt=%9d  BlkSz=%6d NBlk=%6d ", 
+                   (int)dobj->mNVar, (int)dobj->mNEnt, (int)dobj->mBlk, (int)dobj->mNBlk);
+s.PutLine(strg);
+s.PutI4(dobj->mNVar);
+s.PutBytes(dobj->mNames, dobj->mNVar*LENNAME1);
+s.PutI4(dobj->mNEnt);
+s.PutI4(dobj->mBlk);
+s.PutI4(dobj->mNBlk);
+if (dobj->mInfo)  s << (*(dobj->mInfo));
+int jb;
+for(jb=0; jb<dobj->mNBlk; jb++)
+  s.PutR4s(dobj->mPtr[jb], dobj->mNVar*dobj->mBlk); 
+return;
+}
+
+/* --Methode-- */
+//++
+void  ObjFileIO<NTuple>::ReadSelf(PInPersist& s)
+//
+//	Lecture ppersist du ntuple.
+//--
+{
+
+dobj->Clean();
+
+char strg[256];
+s.GetLine(strg, 255);
+// Pour savoir s'il y avait un DVList Info associe
+bool hadinfo = false;
+if (strncmp(strg+strlen(strg)-7, "HasInfo", 7) == 0)  hadinfo = true;
+
+s.GetI4(dobj->mNVar);
+dobj->mNames = new char[dobj->mNVar*LENNAME1];
+dobj->mVar = new r_4[dobj->mNVar];
+dobj->mVarD = new r_8[dobj->mNVar];
+s.GetBytes(dobj->mNames, dobj->mNVar*LENNAME1);
+s.GetI4(dobj->mNEnt);
+s.GetI4(dobj->mBlk);
+s.GetI4(dobj->mNBlk);
+
+if (hadinfo) {    // Lecture eventuelle du DVList Info
+  if (dobj->mInfo == NULL)  dobj->mInfo = new DVList;
+  s >> (*(dobj->mInfo));
+  }
+
+int jb; 
+for(jb=0; jb<dobj->mNBlk; jb++) {
+  r_4* pt = new r_4[dobj->mNVar*dobj->mBlk];
+  dobj->mPtr.push_back(pt);
+  s.GetR4s(dobj->mPtr[jb], dobj->mNVar*dobj->mBlk); 
+}
+
+}
+
+#ifdef __CXX_PRAGMA_TEMPLATES__
+#pragma define_template ObjFileIO<NTuple>
+#endif
+
+#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
+template class ObjFileIO<NTuple>;
+#endif
Index: /trunk/SophyaLib/HiStats/ntuple.h
===================================================================
--- /trunk/SophyaLib/HiStats/ntuple.h	(revision 763)
+++ /trunk/SophyaLib/HiStats/ntuple.h	(revision 763)
@@ -0,0 +1,109 @@
+// This may look like C code, but it is really -*- C++ -*-
+// Class NTuple
+//  CMV+Reza     Juillet 97
+//  CEA-DAPNIA      LAL-IN2P3/CNRS
+
+#ifndef NTUPLE_H_SEEN
+#define NTUPLE_H_SEEN
+
+#include "objfio.h"
+
+#include <iostream.h>
+#include <string>
+#include <vector>
+
+#include "ntupintf.h"
+#include "ppersist.h"
+#include "dvlist.h"
+
+namespace SOPHYA {
+
+class NTuple : public AnyDataObj , public NTupleInterface {
+public:
+//  enum {classId = ClassId_NTuple };
+
+  		    NTuple(int nvar, char** noms, int blk=512);
+		    NTuple();
+                    NTuple(const NTuple& NT);
+		    NTuple(char* flnm);
+  virtual           ~NTuple();
+
+  NTuple&           operator = (const NTuple& NT);
+ 
+  void              Fill(r_4* x);
+
+  inline int_4      NEntry() const  { return(mNEnt); }  
+  inline int_4      NVar() const { return(mNVar); }  
+  inline int_4      BLock() const { return(mBlk); }  
+// $CHECK$ Reza 21/10/99 Pourquoi faire BLock() ?
+
+  float             GetVal(int n, int k)   const;
+  inline float	    GetVal(int n, const char* nom) const 
+                            { return(GetVal(n, IndexNom(nom)) ); }
+  int               IndexNom(const char* nom)  const ;
+  char*             NomIndex(int k) const;
+
+  r_4*	            GetVec(int n, r_4* ret=NULL)  const ;
+  r_8*	            GetVecD(int n, r_8* ret=NULL)  const ;
+
+// Impression, I/O
+  void              Print(int num, int nmax=1)  const ;
+  void              Show(ostream& os) const;
+  inline void       Show() const { Show(cout); }
+
+  DVList&           Info();
+
+// Remplissage depuis fichier ASCII 
+  int		    FillFromASCIIFile(string const& fn, float defval=0.);
+
+
+// Declaration de l interface NTuple
+  virtual uint_4        NbLines() const ;
+  virtual uint_4        NbColumns() const ;
+  virtual r_8 *         GetLineD(int n) const ;
+  virtual r_8		GetCell(int n, int k) const ;
+  virtual r_8		GetCell(int n, string const & nom) const ;
+  virtual void		GetMinMax(int k, double& min, double& max)   const ; 
+  virtual void		GetMinMax(string const & nom, double& min, double& max)   const ; 
+  virtual int           ColumnIndex(string const & nom)  const ;
+  virtual string        ColumnName(int k) const;
+  virtual string	VarList_C(const char* nomx=NULL) const ;
+  virtual string        LineHeaderToString() const;
+  virtual string        LineToString(int n) const;  
+
+//  Pour la gestion de persistance
+  friend class ObjFileIO<NTuple> ;
+
+private:
+  void  Clean();
+
+  int_4 mNVar, mNEnt, mBlk, mNBlk;
+  r_4* mVar;
+  r_8* mVarD;
+  char* mNames;
+
+  vector<r_4*> mPtr;
+
+  DVList* mInfo;        // Infos (variables) attachees au NTuple
+  
+};
+
+inline ostream& operator << (ostream& s, NTuple const & nt)
+  {  nt.Show(s);  return(s);  }
+
+inline POutPersist& operator << (POutPersist& os, NTuple & obj)
+{ ObjFileIO<NTuple> fio(&obj);  fio.Write(os);  return(os); }
+inline PInPersist& operator >> (PInPersist& is, NTuple & obj)
+{ ObjFileIO<NTuple> fio(&obj);  fio.Read(is);  return(is); }
+
+// Classe pour la gestion de persistance
+// ObjFileIO<NTuple>
+
+#ifdef __MWERKS__
+__MSL_FIX_ITERATORS__(r_4*);
+#endif
+
+} // namespace SOPHYA
+
+#endif
+  
Index: /trunk/SophyaLib/HiStats/objlist.list
===================================================================
--- /trunk/SophyaLib/HiStats/objlist.list	(revision 763)
+++ /trunk/SophyaLib/HiStats/objlist.list	(revision 763)
@@ -0,0 +1,6 @@
+hisprof.o
+histinit.o
+histos.o
+histos2.o
+ntuple.o
+xntuple.o
Index: /trunk/SophyaLib/HiStats/xntuple.cc
===================================================================
--- /trunk/SophyaLib/HiStats/xntuple.cc	(revision 763)
+++ /trunk/SophyaLib/HiStats/xntuple.cc	(revision 763)
@@ -0,0 +1,932 @@
+// -*- C++ -*- 
+// 
+// xntuple.cc 
+// N. Regnault   -      98-99 
+// Extended NTuples 
+// 
+#include <string.h>
+#include "ppersist.h"
+#include "perrors.h"
+#include "xntuple.h"
+
+#define LENNAME 31 
+#define MAXLEN 128
+#define BADVAL -1.e19
+
+//++
+// Class	XNTuple
+// Lib		Outils++ 
+// include	xntuple.h
+//
+//	Classe de ntuples
+//--
+//++
+// Links	Parents
+// PPersist
+// NTupleInterface
+//--
+//++
+// Links	Voir aussi
+// NTuple
+//--
+
+
+char* XNTuple::glob_swp = NULL ; 
+long int  XNTuple::NbNT = 0 ; 
+
+
+NTBlk::NTBlk(int ndvar, int nfvar, int nivar, int nsvar, int strsz, int sz) 
+    : sw(0), // pas swappe 
+      swoff(-1), 
+      ddata(NULL), fdata(NULL), 
+      idata(NULL), sdata(NULL)
+{
+    if(ndvar) ddata = new double[ndvar*sz] ; 
+    if(nfvar) fdata = new float[nfvar*sz] ; 
+    if(nivar) idata = new int_4[nivar*sz] ; 
+    if(nsvar) sdata = new char[nsvar*(strsz+1)*sz] ; 
+}
+
+NTBlk::~NTBlk()
+{
+    free() ; 
+}
+
+void NTBlk::free() 
+{
+    if(ddata) delete[] ddata ; ddata = NULL ; 
+    if(fdata) delete[] fdata ; fdata = NULL ; 
+    if(idata) delete[] idata ; idata = NULL ; 
+    if(sdata) delete[] sdata ; sdata = NULL ; 
+}
+
+// 
+// relit bloc ds fichier swap 
+// (seulement si bloc n'est plus en memoire) 
+// 
+void XNTuple::read_blk(NTBlk* blk) const 
+{
+    // deja en memoire ? 
+    if( !blk->sw ) return ;  
+    
+    blk->free() ; 
+    
+    // On va --> blk 
+    fseek(swf, blk->swoff, SEEK_SET) ; 
+    
+    if(mD) { 
+	blk->ddata = new double[mD*mBlkSz] ; 
+	fread(blk->ddata, sizeof(double), mD*mBlkSz, swf) ; 
+    }
+    if(mF) { 
+	blk->fdata = new float[mF*mBlkSz] ; 
+	fread(blk->fdata, sizeof(float), mF*mBlkSz, swf) ; 
+    }
+    if(mI) { 
+	blk->idata = new int_4[mI*mBlkSz] ; 
+	fread(blk->idata, sizeof(int), mI*mBlkSz, swf) ; 
+    }
+    if(mS) { 
+	blk->sdata = new char[mS*(mStrSz+1)*mBlkSz] ; 
+	fread(blk->sdata, sizeof(char), mS*mBlkSz*(mStrSz+1), swf) ; 
+    }
+}
+
+// 
+// swappe le bloc present en memoire 
+// depuis le + longtemps 
+// 
+void XNTuple::swap() const 
+{
+    // fichier de swap ouvert ? 
+    if(!swf) { 
+	swf_name = new char[swp.length()+64] ; 
+	strncpy(swf_name, swp.c_str(), swp.length()+1) ; 
+	strcat(swf_name, "NTupleXXXXXX") ; 
+	//	swf_name = strdup("NTupleXXXXXX") ; 
+	mktemp(swf_name) ; 
+	swf = fopen(swf_name, "w+") ; 
+	if(!swf) THROW(fileErr) ; 
+    }
+    
+    // bloc a swapper 
+    NTBlk* blk = sw.front() ; 
+    
+    // bloc dans fichiers de swap ? 
+    //    if( blk->swoff > 0 ) 
+    write_blk(blk) ; 
+    blk->free() ; 
+    
+    blk->sw = 1 ; // on marque le bloc swappe 
+    sw.pop_front() ; // on le vire de la liste de swap 
+    mNSwBlk++ ; 
+}
+
+
+// 
+// unswap() 
+// 
+void XNTuple::get_blk(int i) const 
+{
+    if(i<0 || i>=mNBlk) return ; 
+    read_blk(ptr[i]) ; 
+    ptr[i]->sw = 0 ; 
+    sw.push_back(ptr[i]) ; 
+    mNSwBlk-- ; 
+}
+
+
+// 
+// ecrit le bloc en bout de fichier 
+// (seulement si swoff < 0 ) 
+// 
+void XNTuple::write_blk(NTBlk* blk) const 
+{
+    // deja swappe ? 
+    if( blk->swoff >= 0 ) return ; 
+    fseek(swf, (off_t)0, SEEK_END) ; 
+    // position debut du bloc 
+    blk->swoff = ftell(swf) ; 
+    if(blk->ddata) fwrite(blk->ddata, sizeof(double), mD*mBlkSz, swf) ; 
+    if(blk->fdata) fwrite(blk->fdata, sizeof(float),  mF*mBlkSz, swf) ; 
+    if(blk->idata) fwrite(blk->idata, sizeof(int),    mI*mBlkSz, swf) ; 
+    if(blk->sdata) fwrite(blk->sdata, sizeof(char),   mS*mBlkSz*(mStrSz+1), swf) ; 
+}
+
+
+void  XNTuple::add_blk() 
+{
+    // l'ancien bloc en ecriture peut maintenant 
+    // etre swappe (s'il existe) 
+    if(mBlk>=0) sw.push_back(ptr[mBlk]) ; 
+    
+    // si pas de place, on swappe 
+    if( mNBlk+1 >= mMaxBlk ) swap() ; 
+    
+    // nouveau bloc 
+    NTBlk* tmp = new NTBlk(mD, mF, mI, mS, mStrSz, mBlkSz) ; 
+    ptr.push_back(tmp) ; 
+    
+    //    mNBlk++ ; 
+    //    mOff = 0 ; 
+    //    mBlk++ ; 
+}
+
+
+XNTuple::XNTuple()
+    : mNEnt(0), mNBlk(0), mNSwBlk(0), 
+      mBlkSz(0), mBlk(-1), mOff(0), 
+      mMaxBlk(0), mStrSz(0), 
+      mD(0), mF(0), mI(0), mS(0), 
+      mVarD(NULL), mMin(NULL), mMax(NULL), 
+      mNVars(0), mNames(NULL), 
+  swf(NULL), swf_name(NULL), mInfo(NULL)
+{
+    if(!glob_swp) SetSwapPath("/tmp/") ; 
+    swp = glob_swp ; 
+    NbNT++ ; 
+}
+
+//++
+// Titre	Constructeurs
+//--
+//++
+// XNTuple(int ndvar, int nfvar, int nivar, int nsvar, char** vnames, -
+//	   int blk=512, int maxblk=100, int strsz=30)
+//	Constructeur - Création d'un XNTuple de "ndvar" variables de type "double", -
+//                      "nfvar" de type "float", "nivar" de type "int" et "nsvar" de type -
+//                      "char*". 
+//|            * "blk"    = taille blocs de données (en nombre d'entrées) 
+//|            * "maxblk" = nombre max de blocs présents en mémoire (non swappés) 
+//|            * "strsz"  = taille des données de type char* 
+// XNTuple(XNTuple const& nt)
+//	Constructeur de copie.
+// XNTuple(string const& flnm)
+//	Constructeur - lecture à partir d'un fichier PPersist.
+//--
+
+XNTuple::XNTuple(int ndvar, int nfvar, int nivar, int nsvar, 
+		   char** vnames, 
+		   int blk, int maxblk, int strsz)
+    : mNEnt(0), mNBlk(0), mNSwBlk(0), 
+      mBlkSz(blk), mBlk(-1), mOff(0), 
+      mMaxBlk(maxblk), mStrSz(strsz), 
+      mD(ndvar), mF(nfvar), mI(nivar), mS(nsvar), 
+      mVarD(NULL), mMin(NULL), mMax(NULL), 
+      mNVars(ndvar+nfvar+nivar+nsvar), mNames(NULL), 
+      swf(NULL), swf_name(NULL), mInfo(NULL)
+{
+    if(!glob_swp) SetSwapPath("/tmp/") ; 
+    swp = glob_swp ; 
+    
+    mVarD = new r_8[mNVars] ; 
+    mMin  = new r_8[mD+mF+mI] ; 
+    int i;
+    for(i = 0 ; i < (mD+mF+mI) ; i++) mMin[i] =  9E19 ; 
+    mMax  = new r_8[mD+mF+mI] ; 
+    for(i = 0 ; i < (mD+mF+mI) ; i++) mMax[i] = -9E19 ; 
+    
+    if(mNVars) mNames = new char[mNVars*(LENNAME+1)] ; 
+    memset(mNames, 0, mNVars*(LENNAME+1));
+    for(i = 0 ; i < mNVars ; i++)
+      strncpy(mNames+i*(LENNAME+1), vnames[i], LENNAME);
+
+    NbNT++ ; 
+}
+
+
+
+XNTuple::XNTuple(string const& flnm) 
+    : mNEnt(0), mNBlk(0), mNSwBlk(0), 
+      mBlkSz(0), mBlk(-1), mOff(0), 
+      mMaxBlk(0), mStrSz(0), 
+      mD(0), mF(0), mI(0), mS(0), 
+      mVarD(NULL), mMin(NULL), mMax(NULL), 
+      mNVars(0), mNames(NULL), 
+      swf(NULL), swf_name(NULL), mInfo(NULL)
+{
+    if(!glob_swp) SetSwapPath("/tmp/") ; 
+    swp = glob_swp ; 
+    
+    PInPersist s(flnm);
+    ObjFileIO<XNTuple> fiont(this);
+    fiont.Read(s);
+    NbNT++ ; 
+}
+
+
+XNTuple::XNTuple(XNTuple const& nt) 
+    : mNEnt(0), mNBlk(0), mNSwBlk(0), 
+      mBlkSz(0), mBlk(-1), mOff(0), 
+      mMaxBlk(0), mStrSz(0), 
+      mD(0), mF(0), mI(0), mS(0), 
+      mVarD(NULL), mMin(NULL), mMax(NULL), 
+      mNVars(0), mNames(NULL), 
+      swf(NULL), swf_name(NULL), mInfo(NULL) 
+{
+    if(!glob_swp) SetSwapPath("/tmp/") ; 
+    swp = glob_swp ; 
+    Copy(nt) ; 
+    NbNT++ ; 
+}
+
+
+
+XNTuple::~XNTuple()
+{
+    clean() ; 
+    NbNT-- ; 
+    if(NbNT==0) delete[] glob_swp ; 
+}
+
+
+
+void XNTuple::clean()
+{
+    // On libere tous les blocs 
+    for(int i = 0 ; i < mNBlk ; i++) 
+	if(!ptr[i]->sw) delete ptr[i] ; 
+    ptr.erase(ptr.begin(), ptr.end()) ; 
+    sw.erase(sw.begin(), sw.end()) ; 
+    
+    // on ferme le fichier de swap 
+    if(swf) fclose(swf) ; swf = NULL ; 
+    remove(swf_name) ; 
+    delete swf_name ; swf_name = NULL ; 
+    
+    // tout le monde remis a 0 
+    mNEnt   = 0 ;     mNBlk   = 0 ; 
+    mNSwBlk = 0 ;     mBlkSz  = 0 ;  
+    mBlk    = -1 ;    mOff    = 0 ; 
+    mMaxBlk = 0 ;     mStrSz  = 0 ; 
+    mD      = 0 ;     mF      = 0 ; 
+    mI      = 0 ;     mS      = 0 ; 
+    mNVars  = 0 ; 
+    
+    
+    if(mVarD) { delete[] mVarD ; mVarD = NULL ; } 
+    if(mMin)  { delete[] mMin  ; mMin  = NULL ; } 
+    if(mMax)  { delete[] mMax  ; mMax  = NULL ; } 
+
+    if(mNames) {delete[] mNames ; mNames = NULL ; }
+    if (mInfo) delete mInfo;
+
+}
+
+
+//++
+// Titre	Méthodes
+//--
+//++
+//  void Fill(double* d_data, float* f_data, int* i_data, char** s_data) 
+//	Remplissage d'une ligne dans le NTuple
+//  void Show(ostream& os) const 
+//	Impression de la liste des variables avec min-max sur le flot "oß
+//  void Show() const 
+//	Identique à "Show(cout)"
+//  XNTuple&   operator = (XNTuple const& nt) 
+//	Opérateur égal (=) , copie "nt" dans le premier NTuple
+//--
+void XNTuple::Fill(double* d_data, float* f_data, int* i_data, char** s_data) 
+{
+    // place disponible dans bloc courant ? 
+    if( mOff==mBlkSz || mOff==0 ) {
+	add_blk() ; mOff = 0 ; mNBlk++ ; mBlk++ ; 
+    }
+    
+    if( mD && !ptr[mBlk]->ddata || 
+	mF && !ptr[mBlk]->fdata || 
+	mI && !ptr[mBlk]->idata || 
+	mS && !ptr[mBlk]->sdata )
+	THROW(parmErr) ; 
+    
+    double x ; 
+    // copie variables et update mMin, mMax 
+    if(mD) {
+	memcpy(ptr[mBlk]->ddata+mOff*mD, d_data, mD*sizeof(double)) ; 
+	for(int i = 0 ; i < mD ; i++) {
+	    x = d_data[i] ; 
+	    if(x<mMin[i]) mMin[i] = x ; 
+	    if(x>mMax[i]) mMax[i] = x ; 
+	}
+    }
+    
+    if(mF) {
+	memcpy(ptr[mBlk]->fdata+mOff*mF, f_data, mF*sizeof(float)) ; 
+	for(int i = 0 ; i < mF ; i++) {
+	    x = f_data[i] ; 
+	    if(x<mMin[i+mD]) mMin[i+mD] = x ; 
+	    if(x>mMax[i+mD]) mMax[i+mD] = x ; 
+	}
+    }
+
+    if(mI) {
+	memcpy(ptr[mBlk]->idata+mOff*mI, i_data, mI*sizeof(int)) ; 
+	for(int i = 0 ; i < mI ; i++) {
+	    x = i_data[i] ; 
+	    if(x<mMin[i+mD+mF]) mMin[i+mD+mF] = x ; 
+	    if(x>mMax[i+mD+mF]) mMax[i+mD+mF] = x ; 
+	}
+    }
+    for(int i = 0 ; i < mS ; i++) 
+	memcpy(ptr[mBlk]->sdata+(mOff*mS+i)*(mStrSz+1), 
+	       s_data[i], (strlen(s_data[i])+1)*sizeof(char)) ; 
+    mOff++ ; 
+    mNEnt++ ; 
+}
+
+// 
+// A quel index correspond mon nom ? 
+// 
+int XNTuple::IndexNom(char const* nom) const 
+{
+    int i ; 
+    for(i = 0 ; i < (mD+mF+mI+mS) ; i++)
+	if( !strncmp( mNames+i*(LENNAME+1), nom, LENNAME+1) ) 
+	    return i ; 
+    return -1 ; 
+}
+
+string XNTuple::NomIndex(int k) const 
+{
+    if( k<0 || k>=mNVars ) return "" ; 
+    return mNames + k*(LENNAME+1) ; 
+}
+
+
+// 
+// 
+// 
+double XNTuple::GetDVal(int i, int k) const 
+{
+    if( i<0 || i>=mNEnt || k<0 || k>=mD ) 
+	THROW(rangeCheckErr) ; 
+    
+    // Bloc ? 
+    int blk = (int)(i/mBlkSz) ; 
+    int off = i%mBlkSz ; 
+    
+    // bloc swappe ? 
+    if( ptr[blk]->sw ) { get_blk(blk) ; swap() ; } 
+    return ptr[blk]->ddata[off*mD+k] ; 
+}
+
+// 
+// 
+// 
+float XNTuple::GetFVal(int i, int k) const 
+{
+    if( i<0 || i>=mNEnt || k<mD || k>=(mD+mF) ) 
+	THROW(rangeCheckErr) ; 
+    k -= mD ; 
+    
+    // Bloc ? 
+    int blk = (int)(i/mBlkSz) ; 
+    int off = i%mBlkSz ; 
+    
+    // bloc swappe ? 
+    if( ptr[blk]->sw ) { get_blk(blk) ; swap() ; } 
+    return ptr[blk]->fdata[off*mF+k] ; 
+}
+
+// 
+// 
+// 
+int XNTuple::GetIVal(int i, int k) const 
+{
+    if( i<0 || i>=mNEnt || k<(mD+mF) || k>=(mD+mF+mI) )
+	THROW(rangeCheckErr) ; 
+    k -= (mD+mF) ; 
+
+    // Bloc ? 
+    int blk = (int)(i/mBlkSz) ; 
+    int off = i%mBlkSz ; 
+    
+    // bloc swappe ? 
+    if( ptr[blk]->sw ) { get_blk(blk) ; swap() ; } 
+
+    return ptr[blk]->idata[off*mI+k] ; 
+}
+
+// 
+// 
+// 
+string  XNTuple::GetSVal(int i, int k) const 
+{
+    if( i<0 || i>=mNEnt || k<(mD+mF+mI) || k>=(mD+mF+mI+mS) )
+	THROW(rangeCheckErr) ; 
+    k -= (mD+mF+mI) ; 
+    
+    // Bloc ? 
+    int blk = (int)(i/mBlkSz) ; 
+    int off = i%mBlkSz ; 
+    
+    // bloc swappe ? 
+    if( ptr[blk]->sw ) { get_blk(blk) ; swap() ; } 
+    
+    // copie de la chaine 
+    // string ret = strdup(ptr[blk]->sdata + (off*mS+k)*(mStrSz+1)) ; // $CHECK$ EA fuite de memoire
+       // attention, strdup fait un malloc et on ne fait jamais de free...
+       // a quoi sert ce strdup ?????????
+    string ret = ptr[blk]->sdata + (off*mS+k)*(mStrSz+1) ; 
+    //    return ptr[blk]->sdata[k] ; 
+    return ret ; 
+}
+
+
+// 
+// Copie bloc a bloc, avec meme parametres 
+//  
+void XNTuple::Copy(XNTuple const& nt) 
+{
+    clean() ; 
+    // Parametres 
+    mNEnt = nt.mNEnt ; 
+    mBlkSz = nt.mBlkSz ; 
+    mOff = nt.mOff ; 
+    mMaxBlk = nt.mMaxBlk ; 
+    mStrSz = nt.mStrSz ; 
+    mD = nt.mD ; 
+    mF = nt.mF ; 
+    mI = nt.mI ; 
+    mS = nt.mS ; 
+    mNVars = nt.mNVars ; 
+    
+    // noms 
+    if(mNVars) {
+	mNames = new char[mNVars*(LENNAME+1)] ; 
+        mVarD = new double[mNVars];
+	memcpy(mNames, nt.mNames, mNVars*(LENNAME+1)*sizeof(char)) ; 
+    }
+    // MinMax 
+    if(nt.mMin) { 
+	mMin = new double[(mD+mF+mI)] ; 
+	memcpy(mMin, nt.mMin, (mD+mF+mI)*sizeof(double)) ; 
+    }
+    if(nt.mMax) { 
+	mMax = new double[(mD+mF+mI)] ; 
+	memcpy(mMax, nt.mMax, (mD+mF+mI)*sizeof(double)) ; 
+    }
+    
+    //dup blocs 
+    mNBlk = 0 ; 
+    mBlk  = -1 ; 
+    for(int i = 0 ; i < nt.mNBlk ; i++) {
+	add_blk() ; mBlk++ ; mNBlk++ ; 
+	// si nt.ptr[i] swappe, on le relit 
+	if(nt.ptr[i]->sw) nt.read_blk(nt.ptr[i]) ; 
+	if(mD) 
+	    memcpy(ptr[i]->ddata, nt.ptr[i]->ddata, mD*mBlkSz*sizeof(double)) ; 
+	if(mF)
+	    memcpy(ptr[i]->fdata, nt.ptr[i]->fdata, mF*mBlkSz*sizeof(float)) ; 
+	if(mI)
+	    memcpy(ptr[i]->idata, nt.ptr[i]->idata, mF*mBlkSz*sizeof(int)) ; 
+	if(mS)
+	    memcpy(ptr[i]->sdata, nt.ptr[i]->sdata, mS*mBlkSz*sizeof(char)*(mStrSz+1)) ; 
+	if(nt.ptr[i]->sw) nt.ptr[i]->free() ; 
+    }
+
+    // DVList Info block 
+    if(nt.mInfo!=NULL) {mInfo = new DVList; *mInfo = *(nt.mInfo);}
+
+}
+
+//++
+// DVList& Info()
+//	Renvoie une référence sur l'objet DVList Associé
+//--
+
+DVList&  XNTuple::Info()
+{
+  if (mInfo == NULL)  mInfo = new DVList;
+  return(*mInfo);
+}
+
+void XNTuple::Print(int num, int nmax) const 
+{
+    printf("XNTuple::Print() : \n") ; 
+    printf("  Entrees = %d, Blocs  = %d, Bloc Size = %d\n", 
+	   mNEnt, mNBlk, mBlkSz) ; 
+    int i,j;
+    printf("  D_Vars  = %d : ", mD) ; 
+    for(i = 0 ; i < mD ; i++) 
+	printf("%s ", NomIndex(i).c_str() ) ; 
+    printf("\n") ; 
+    
+    printf("  F_Vars  = %d : ", mF) ; 
+    for(i = 0 ; i < mF ; i++) 
+	printf("%s ", NomIndex(i+mD).c_str() ) ; 
+    printf("\n") ; 
+    
+    printf("  I_Vars  = %d : ", mI) ; 
+    for(i = 0 ; i < mI ; i++) 
+	printf("%s ", NomIndex(i+mD+mF).c_str() ) ; 
+    printf("\n") ; 
+    
+    printf("  S_Vars  = %d : ", mS) ; 
+    for(i = 0 ; i < mS ; i++) 
+	printf("%s ", NomIndex(i+mD+mF+mI).c_str() ) ; 
+    printf("\n") ; 
+    
+    for(i = num ; i < num+nmax ; i++) {
+	for(j = 0 ; j < mD ; j++) printf("%f ", GetDVal(i,j)) ; 
+	printf(" -- ") ; 
+	for(j = 0 ; j < mF ; j++) printf("%f ", GetFVal(i,j+mD)) ; 
+	printf(" -- ") ; 
+	for(j = 0 ; j < mI ; j++) printf("%d ", GetIVal(i,j+mD+mF)) ; 
+	printf(" -- ") ; 
+	for(j = 0 ; j < mS ; j++) printf("%s ", GetSVal(i,j+mD+mF+mI).c_str()) ; 
+	printf("\n") ; 
+    }
+}
+
+
+void   XNTuple::Show(ostream& os) const 
+{
+    os << "XNTuple: NVar= " << mNVars << " NEnt= " << mNEnt 
+       << " (BlkSz,NBlk= " << mBlkSz << ", " << mNBlk << ")" 
+       << " (mNSwBlk= " << mNSwBlk << ")" << endl ; 
+    os << "(Sw File= " ; if(swf_name) os << swf_name ; os << ")" << endl ; 
+    
+    char* buff = new char[80] ; 
+    sprintf(buff, "Variables   :  Type         Min         Max  \n") ; 
+    os << buff ; 
+    
+    double min, max ;
+    int i; 
+    for(i = 0 ; i < mD ; i++) {
+	GetMinMax(i, min, max) ; 
+	sprintf(buff, "  %-10s:   D    %12.6g   %12.6g \n", 
+		NomIndex(i).c_str(), min, max) ; 
+	os << buff ; 
+    }
+    for(i = 0 ; i < mF ; i++) {
+	GetMinMax(i+mD, min, max) ; 
+	sprintf(buff, "  %-10s:   F    %12.6g   %12.6g \n", 
+		NomIndex(i+mD).c_str(), min, max) ; 
+	os << buff ; 
+    }
+    for(i = 0 ; i < mI ; i++) {
+	GetMinMax(i+mD+mF, min, max) ; 
+	sprintf(buff, "  %-10s:   I    %12.6g   %12.6g \n", 
+		NomIndex(i+mD+mF).c_str(), min, max) ; 
+	os << buff ; 
+    }
+    for(i = 0 ; i < mS ; i++) {
+	sprintf(buff, "  %-10s:   S      -----        ----- \n", 
+		NomIndex(i+mD+mF+mI).c_str()) ; 
+	os << buff ; 
+    }
+    delete[] buff ; 
+}
+
+int    XNTuple::FillFromASCIIFile(string const& fn, double ddval,  float dfval,
+                                        int dival, const char * dsval)
+//	Remplit le ntuple a partir d'un fichier ASCII.
+//	Renvoie le nombre de lignes ajoutees.
+{
+// a faire
+return(0);
+}
+
+
+void  XNTuple::SetSwapPath(char* p) 
+{
+    if(!glob_swp) glob_swp = new char[MAXLEN+1] ; 
+    strncpy(glob_swp, p, MAXLEN) ; 
+}
+
+// 
+// 
+// Interface NTuple 
+// 
+// 
+uint_4 XNTuple::NbLines() const
+{
+    return(NEntry());
+}
+
+uint_4 XNTuple::NbColumns() const
+{
+    return(NVar());
+}
+
+r_8* XNTuple::GetLineD(int n) const 
+{
+    // memcpy() impossible
+    // il faut faire des GetVal 
+    int i;
+    for(i = 0 ; i < mD ; i++) 
+	mVarD[i] = GetDVal(n, i) ; 
+    for(i = 0 ; i < mF ; i++) 
+	mVarD[i+mD] = (double)GetFVal(n, i+mD) ; 
+    for(i = 0 ; i < mI ; i++) 
+	mVarD[i+mD+mF] = (double)GetIVal(n, i+mD+mF) ; 
+    for(i = 0 ; i < mS ; i++) 
+	mVarD[i+mD+mF+mI] = atof(GetSVal(n, i+mD+mF+mI).c_str()) ; 
+    return mVarD ; 
+}
+
+
+r_8  XNTuple::GetCell(int n, int k) const 
+{
+    if( k<0 || k>=(mD+mF+mI+mS) ) return BADVAL ; 
+    if(k<mD) return GetDVal(n,k) ; 
+    if(k<(mF+mD)) return GetFVal(n,k) ; 
+    if(k<(mF+mD+mI)) return GetIVal(n,k) ; 
+    if(k<(mD+mF+mI+mS)) return atof(GetSVal(n,k).c_str()) ; 
+    return BADVAL ; 
+}
+
+
+r_8  XNTuple::GetCell(int n, string const& nom) const 
+{
+    int k = IndexNom(nom.c_str()) ; 
+    if(k<0) return BADVAL ; 
+    if(k<mD) return GetDVal(n,k) ; 
+    if(k<(mF+mD)) return GetFVal(n,k) ; 
+    if(k<(mF+mD+mI)) return GetIVal(n,k) ; 
+    if(k<(mD+mF+mI+mS)) return atof(GetSVal(n,k).c_str()) ; 
+    return BADVAL ; 
+}
+
+string   XNTuple::GetCelltoString(int n, int k) const
+{
+    char buff[32];
+    if(k<0) return (""); 
+    else if(k<mD) 
+      { sprintf(buff, "%g", GetDVal(n,k)); return(buff) ; } 
+    else if(k<(mF+mD)) 
+      { sprintf(buff, "%g", (double)GetFVal(n,k)); return(buff) ; } 
+    else if(k<(mF+mD+mI)) 
+      { sprintf(buff, "%d", GetIVal(n,k)); return(buff) ; } 
+    else if(k<(mD+mF+mI+mS)) 
+      return(GetSVal(n,k).c_str()) ; 
+    return("");
+}
+
+void XNTuple::GetMinMax(int k, double& min, double& max)  const
+{
+    min = 9E19 ; max = -9E19 ; 
+    // variables string non traitees 
+    if( k<0 || k>=(mD+mF+mI) ) return ; 
+    
+    min = mMin[k] ; 
+    max = mMax[k] ; 
+    return ; 
+}
+
+
+void XNTuple::GetMinMax(string const & nom, double& min, double& max)   const
+{
+    GetMinMax(IndexNom(nom.c_str()), min, max) ; 
+}
+
+
+int XNTuple::ColumnIndex(string const& nom) const 
+{
+    return IndexNom(nom.c_str()) ; 
+}
+
+
+string XNTuple::ColumnName(int k) const 
+{
+    return NomIndex(k) ; 
+}
+
+
+string XNTuple::VarList_C(const char* nomx) const 
+{
+    string rets = "" ; 
+    rets += "/* double type variables */ \n";
+    // variables double 
+    int i;
+    for(i = 0 ; i < mD ; i++) {
+	if( i>0 && (i%5)==0 ) rets += ";" ; 
+	if( (i%5)==0 ) rets+="\ndouble " ; 
+	else rets += ", " ; 
+	rets += NomIndex(i) ; 
+    }
+    if(mD) rets += ";" ; 
+    
+    // variables float 
+    rets += "/* float type variables */ \n";
+    for(i = 0 ; i < mF ; i++) {
+	if( i>0 && (i%5)==0 ) rets += ";" ; 
+	if( (i%5)==0 ) rets+="\ndouble " ; 
+	else rets += ", " ; 
+	rets += NomIndex(i+mD) ; 
+    }
+    if(mF) rets += ";" ; 
+    
+    // variables int 
+    rets += "/* int type variables */ \n";
+    for(i = 0 ; i < mI ; i++) {
+	if( i>0 && (i%5)==0 ) rets += ";" ; 
+	if( (i%5)==0 ) rets+="\ndouble " ; 
+	else rets += ", " ; 
+	rets += NomIndex(i+mD+mF) ; 
+    }
+    if(mI) rets += ";" ; 
+    
+    // variables string
+    rets += "/* string type variables */ \n"; 
+    for(i = 0 ; i < mS ; i++) {
+	if( i>0 && (i%5)==0 ) rets += ";" ; 
+	if( (i%5)==0 ) rets+="\ndouble " ; 
+	else rets += ", " ; 
+	rets += NomIndex(i+mD+mF+mI) ; 
+    }
+    rets += "; \n" ; 
+
+    if(nomx) {
+	char buff[256] ; 
+	for(i = 0 ; i < mNVars ; i++) {
+	    sprintf(buff, "%s=%s[%d]; ", NomIndex(i).c_str(), nomx, i) ; 
+	    rets+=buff ; 
+	    if( (i%3 == 0) && (i>0) ) rets += "\n" ; 
+	}
+    }
+    return rets ; 
+}
+
+
+string XNTuple::LineHeaderToString() const 
+{
+    char buff[32];
+    string rets=" Num    ";
+    for(int i=0; i<mNVars; i++) {
+	sprintf(buff, "%8s ", NomIndex(i).c_str() );
+	rets += buff;
+    }
+    rets += '\n';
+    return(rets);
+
+}
+
+
+string XNTuple::LineToString(int n) const 
+{
+    char buff[32];
+    double* val;
+    val = GetLineD(n);
+    sprintf(buff,"%6d: ",n); 
+    string rets=buff;
+    int i;
+    for(i=0; i<(mD+mF+mI); i++) {
+	sprintf(buff, "%8.3g ", val[i]);
+	rets += buff;
+    }
+    
+    for(i = mD+mF+mI ; i < mNVars ; i++) {
+	sprintf(buff, "%8s ", GetSVal(n,i).c_str() ) ; 
+	rets += buff ; 
+    }
+    
+    rets += '\n';
+    return(rets);
+}
+
+
+void   ObjFileIO<XNTuple>::WriteSelf(POutPersist& ppout) const 
+{
+  if (dobj == NULL)   return;
+//  On ecrit 3 uint_4 .... 
+//  0: Numero de version,  1 : non nul -> has info,  2 : reserve
+  uint_4 itab[3];
+  itab[0] = 1;  // Numero de version a 1
+  itab[1] = itab[2] = 0;
+  if (dobj->mInfo)  itab[1] = 1;
+  ppout.Put(itab,3);
+  if (dobj->mInfo)  if (dobj->mInfo)  ppout << (*(dobj->mInfo));
+
+  // variables internes 
+  ppout.PutI4(dobj->mNEnt) ; 
+  ppout.PutI4(dobj->mNBlk) ; 
+  ppout.PutI4(dobj->mBlkSz) ; 
+  ppout.PutI4(dobj->mBlk) ; 
+  ppout.PutI4(dobj->mOff) ; 
+  ppout.PutI4(dobj->mMaxBlk) ; 
+  ppout.PutI4(dobj->mStrSz) ; 
+  ppout.PutI4(dobj->mD) ; 
+  ppout.PutI4(dobj->mF) ; 
+  ppout.PutI4(dobj->mI) ; 
+  ppout.PutI4(dobj->mS) ; 
+  ppout.PutI4(dobj->mNVars) ; 
+    
+  // Noms 
+  ppout.PutBytes(dobj->mNames, dobj->mNVars*(LENNAME+1)) ; 
+
+  // MinMax 
+  ppout.PutR8s(dobj->mMin, (dobj->mD+dobj->mF+dobj->mI)) ; 
+  ppout.PutR8s(dobj->mMax, (dobj->mD+dobj->mF+dobj->mI)) ; 
+    
+  // Ecriture blocs N'ecrire que si datas existent 
+  for(int i = 0 ; i < dobj->mNBlk ; i++) {
+    // si bloc swappe, on le relit en douce ... 
+    if( dobj->ptr[i]->sw )  dobj->read_blk(dobj->ptr[i]) ; 
+    ppout.PutR8s(dobj->ptr[i]->ddata, dobj->mD*dobj->mBlkSz) ; 
+    ppout.PutR4s(dobj->ptr[i]->fdata, dobj->mF*dobj->mBlkSz) ; 
+    ppout.PutI4s(dobj->ptr[i]->idata, dobj->mI*dobj->mBlkSz) ; 
+    ppout.PutBytes(dobj->ptr[i]->sdata, dobj->mS*dobj->mBlkSz*(dobj->mStrSz+1)) ; 
+    if(dobj->ptr[i]->sw) dobj->ptr[i]->free() ; 
+    }
+}
+
+
+void   ObjFileIO<XNTuple>::ReadSelf(PInPersist& ppin) 
+{
+  if (dobj == NULL) dobj = new XNTuple;
+  else    dobj->clean() ; 
+//  On lit 3 uint_4 .... 
+//  0: Numero de version,  1 : non nul -> has info,  2 : reserve
+  uint_4 itab[3];
+  ppin.Get(itab,3);
+  if (itab[1] != 0) { // Lecture eventuelle du DVList Info
+    if (dobj->mInfo == NULL)  dobj->mInfo = new DVList;
+    ppin >> (*(dobj->mInfo));
+  }
+
+  ppin.GetI4(dobj->mNEnt) ; 
+  ppin.GetI4(dobj->mNBlk) ; 
+  ppin.GetI4(dobj->mBlkSz) ; 
+  ppin.GetI4(dobj->mBlk) ; 
+  ppin.GetI4(dobj->mOff) ; 
+  ppin.GetI4(dobj->mMaxBlk) ; 
+  ppin.GetI4(dobj->mStrSz) ; 
+  ppin.GetI4(dobj->mD) ; 
+  ppin.GetI4(dobj->mF) ; 
+  ppin.GetI4(dobj->mI) ; 
+  ppin.GetI4(dobj->mS) ; 
+  ppin.GetI4(dobj->mNVars) ; 
+    
+  // Noms 
+  dobj->mNames = new char[dobj->mNVars*(LENNAME+1)] ;
+  dobj->mVarD = new double[dobj->mNVars]; 
+  ppin.GetBytes(dobj->mNames, dobj->mNVars*(LENNAME+1)) ; 
+
+  // MinMax 
+  dobj->mMin = new r_8[dobj->mD+dobj->mF+dobj->mI] ; 
+  ppin.GetR8s(dobj->mMin, dobj->mD+dobj->mF+dobj->mI) ; 
+  dobj->mMax = new r_8[dobj->mD+dobj->mF+dobj->mI] ; 
+  ppin.GetR8s(dobj->mMax, dobj->mD+dobj->mF+dobj->mI) ; 
+    
+  // lecture Blocs 
+  int nblk = dobj->mNBlk ; 
+  dobj->mBlk = -1 ; dobj->mNBlk = 0 ; 
+  for(int i = 0 ; i < nblk ; i++) {
+    dobj->add_blk() ; dobj->mBlk++ ; dobj->mNBlk++ ; 
+    ppin.GetR8s(dobj->ptr[i]->ddata, dobj->mD*dobj->mBlkSz) ; 
+    ppin.GetR4s(dobj->ptr[i]->fdata, dobj->mF*dobj->mBlkSz) ; 
+    ppin.GetI4s(dobj->ptr[i]->idata, dobj->mI*dobj->mBlkSz) ; 
+    ppin.GetBytes(dobj->ptr[i]->sdata, dobj->mS*dobj->mBlkSz*(dobj->mStrSz+1)) ; 
+  } 
+}
+
+#ifdef __CXX_PRAGMA_TEMPLATES__
+#pragma define_template ObjFileIO<XNTuple>
+#endif
+
+#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
+template class ObjFileIO<XNTuple>;
+#endif
Index: /trunk/SophyaLib/HiStats/xntuple.h
===================================================================
--- /trunk/SophyaLib/HiStats/xntuple.h	(revision 763)
+++ /trunk/SophyaLib/HiStats/xntuple.h	(revision 763)
@@ -0,0 +1,165 @@
+// -*- C++ -*- 
+// 
+// XNTuple 
+// 
+// N. Regnault 03/98 
+//             07/99 (NTupleInterface)
+// 
+// 
+#ifndef XNTUPLE_H_SEEN
+#define XNTUPLE_H_SEEN
+
+#include "objfio.h"
+
+#include <string.h>
+#include <vector>
+#include <list>
+#include <map>
+
+#include "ntupintf.h"
+#include "ppersist.h"
+#include "dvlist.h"
+
+namespace SOPHYA {
+// 
+// Bloc de donnees 
+// 
+struct NTBlk
+{
+    NTBlk(int ndvar, int nfvar, int nivar, int nsvar, int strsz, int sz) ; 
+    ~NTBlk() ; 
+    
+    void free() ;     // ne libere que les zones de donnees 
+
+    int      sw ;     // swappe ? 
+    long int swoff ;  // offset dans fichier de swap 
+    float*   fdata ;  
+    double*  ddata ;  
+    int_4*   idata ;  
+    char*    sdata ;  
+} ; 
+
+
+class XNTuple : public AnyDataObj , public NTupleInterface {
+public:
+  //    enum {classId = ClassId_XNTuple };
+
+    XNTuple() ; 
+    XNTuple(int ndvar, int nfvar, int nivar, int nsvar, 
+	    char** vnames, 
+	    int blk=512, int maxblk=100, int strsz=30) ; 
+    XNTuple(XNTuple const& nt) ; 
+    XNTuple(string const & flnm); 
+    virtual           ~XNTuple();
+    
+    void              Fill(double* d_data, float* f_data, int* i_data, char** sdata);
+    
+    inline int_4      NEntry() const { return mNEnt ; } 
+    inline int_4      NDVar() const  { return mD ; } 
+    inline int_4      NFVar() const  { return mF ; } 
+    inline int_4      NIVar() const  { return mI ; } 
+    inline int_4      NSVar() const  { return mS ; } 
+    inline int_4      NVar() const   { return mNVars ; } 
+    int               IndexNom(char const* nom) const ; 
+    string            NomIndex(int k) const ; 
+    
+    double            GetDVal(int i, int k) const ; 
+    float             GetFVal(int i, int k) const ; 
+    int               GetIVal(int i, int k) const ; 
+    string            GetSVal(int i, int k) const ; 
+    //    MuTyV             GetVal(char const* nom) ; 
+    
+    void              Copy(XNTuple const& nt) ; 
+    inline XNTuple&   operator=(XNTuple const& nt) { Copy(nt) ; return *this ; } 
+    
+    void              Print(int num, int nmax=1) const ;
+    void              Show(ostream& os) const ; 
+    inline void       Show() const { Show(cout) ; } 
+
+    DVList&           Info() ; 
+
+// Remplissage depuis fichier ASCII 
+   int		      FillFromASCIIFile(string const& fn, double ddfval=0., float dfval=0.,
+                                        int dival=0, const char * dsval="");
+        
+    static  void      SetSwapPath(char* p) ; 
+    string&           SwapPath() { return swp ; } 
+    
+// Declaration de l interface NTuple
+    virtual uint_4    NbLines() const ;
+    virtual uint_4    NbColumns() const ;
+    virtual r_8 *     GetLineD(int n) const ;
+    virtual r_8	      GetCell(int n, int k) const ;
+    virtual r_8	      GetCell(int n, string const & nom) const ;
+    virtual string    GetCelltoString(int n, int k) const ;
+    virtual void      GetMinMax(int k, double& min, double& max)   const ; 
+    virtual void      GetMinMax(string const & nom, double& min, double& max)   const ; 
+    virtual int       ColumnIndex(string const & nom)  const ;
+    virtual string    ColumnName(int k) const;
+    virtual string    VarList_C(const char* nomx=NULL) const ;
+    virtual string    LineHeaderToString() const;
+    virtual string    LineToString(int n) const;  
+
+//  Pour la gestion de persistance
+    friend class ObjFileIO<XNTuple> ;
+   
+private:
+    void    clean();    // remet tout a zero 
+    void    swap() const ;    // swappe le + ancien blk en memoire 
+    void    read_blk(NTBlk*) const ;  // relit le bloc 
+    void    write_blk(NTBlk*) const ; // ecrit un bloc 
+    void    get_blk(int i) const ;  // == unswap() 
+    void    add_blk() ; // ajout d'un nouveau bloc 
+    
+    int_4   mNEnt ;   // nb total d'entrees 
+    int_4   mNBlk ;   // nb total de blocs 
+    mutable int_4   mNSwBlk ; // nb blocs swappes 
+    int_4   mBlkSz ;  // taille bloc 
+    int_4   mBlk ;    // bloc courant en ecriture 
+    int_4   mOff ;    // offset dans bloc en ecriture ; pointe sur le premiere case vide 
+    int_4   mMaxBlk ; // nb max de blocs autorises en memoire 
+    int_4   mStrSz ;  // taille des variables char** 
+    int_4   mNVars ;  // mD + mF + mI + mS 
+    int_4   mD, mF, mI, mS ; // nb variables de chaque type 
+    
+    r_8* mVarD, *mMin, *mMax ; 
+    
+    // 
+    // Noms + Prefix --> indiquer le type 
+    // 
+    //    char* mDNames, *mFNames, *mINames, *mSNames ; 
+    char* mNames ; 
+    //    typedef less<string> SCmp ; 
+    //    map<string,int,SCmp> mNames ; 
+    
+    // Donnees 
+    mutable vector<NTBlk*> ptr ; 
+    mutable list<NTBlk*> sw ; // liste --> swap() 
+
+    mutable FILE* swf ; // fichier --> swap 
+    mutable char* swf_name ; 
+    string  swp ; // swap_path 
+
+    DVList* mInfo;        // Infos (variables) attachees au NTuple
+
+    static  char*  glob_swp ; // swap_path global 
+    static  long int NbNT ; 
+}; 
+
+
+inline ostream& operator << (ostream& s, XNTuple const & nt)
+  {  nt.Show(s);  return(s);  }
+
+inline POutPersist& operator << (POutPersist& os, XNTuple & obj)
+{ ObjFileIO<XNTuple> fio(&obj);  fio.Write(os);  return(os); }
+inline PInPersist& operator >> (PInPersist& is, XNTuple & obj)
+{ ObjFileIO<XNTuple> fio(&obj);  fio.Read(is);  return(is); }
+
+// Classe pour la gestion de persistance
+// ObjFileIO<XNTuple>
+
+
+} // namespace SOPHYA
+
+#endif
+
