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