source: Sophya/trunk/SophyaLib/HiStats/histos2.h@ 1058

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

//! 2 dimensions histograms
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);
  Histo2D& operator += (const Histo2D& a);
  Histo2D& operator -= (const Histo2D& a);
  Histo2D& operator *= (const Histo2D& a);
  Histo2D& operator /= (const Histo2D& a);

  // get/put dans/depuis une matrice / vector
  void GetXCoor(TVector<r_8>& v);
  void GetValue(TMatrix<r_8> &v);
  void GetYCoor(TVector<r_8>& v);
  void GetError2(TMatrix<r_8>& v);
  void GetError(TMatrix<r_8>& v);
  void PutValue(TMatrix<r_8>& v, int ierr=0);
  void PutValueAdd(TMatrix<r_8>& v, int ierr=0);
  void PutError2(TMatrix<r_8>& v);
  void PutError2Add(TMatrix<r_8>& v);
  void PutError(TMatrix<r_8>& v);

  // INLINES
  //! Retourne l'abscisse minimum.
  inline float   XMin() const {return xmin;}
  //! Retourne l'abscisse maximum.
  inline float   XMax() const {return xmax;}
  //! Retourne l'ordonnee minimum.
  inline float   YMin() const {return ymin;}
  //! Retourne l'ordonnee maximum.
  inline float   YMax() const {return ymax;}
  //! Retourne le nombre de bins selon X.
  inline int_4   NBinX() const {return nx;}
  //! Retourne le nombre de bins selon Y.
  inline int_4   NBinY() const {return ny;}
  //! Retourne la largeur du bin selon X.
  inline float   WBinX() const {return wbinx;}
  //! Retourne la largeur du bin selon Y.
  inline float   WBinY() const {return wbiny;}
  //! Retourne le pointeur sur le tableaux des contenus.
  inline float*  Bins() const  {return data;}
  //! Retourne le contenu du bin i,j.
  inline float   operator()(int i,int j) const {return data[j*nx+i];}
  //! Remplit le contenu du bin i,j.
  inline float&  operator()(int i,int j)       {return data[j*nx+i];}
  //! retourne "true" si il y a des erreurs stoquees
  inline bool    HasErrors() { if(err2) return true; else return false;}
  //! Retourne l'erreur du bin i,j.
  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.;}
  //! Remplit l'erreur au carre du bin i,j.
  inline double  Error2(int i,int j) const
                   {if(err2) return err2[j*nx+i]; else return 0.;}
  //! Remplit l'erreur au carre du bin i,j.
  inline double& Error2(int i,int j) {return err2[j*nx+i];}
  //! Retourne la somme ponderee.
  inline float   NData() const    {return (float) nHist;}
  //! Retourne le nombre d'entrees.
  inline int_4   NEntries() const {return nEntries;}
  //! Retourne l'abscisse et l'ordonnee du coin inferieur du bin i,j.
  inline void BinLowEdge(int i,int j,float& x,float& y)
                 {x = xmin + i*wbinx; y = ymin + j*wbiny;}
  //! Retourne l'abscisse et l'ordonnee du centre du bin i,j.
  inline void BinCenter(int i,int j,float& x,float& y)
                 {x = xmin + (i+0.5)*wbinx; y = ymin + (j+0.5)*wbiny;}
  //! Retourne l'abscisse et l'ordonnee du coin superieur du bin i,j.
  inline void BinHighEdge(int i,int j,float& x,float& y)
                 {x = xmin + (i+1)*wbinx; y = ymin + (j+1)*wbiny;}
  //! Retourne les numeros du bin contenant l'abscisse et l'ordonnee x,y.
  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();
  //! Retourne le pointeur sur l'histo 1D de la projection selon X.
  inline Histo*  HProjX() const {return hprojx;}
  //! Retourne le pointeur sur l'histo 1D de la projection selon Y.
  inline Histo*  HProjY() const {return hprojy;}
  void           ShowProj();

  // BANDES
  //! Retourne le nombre de bandes selon X
  inline int     NBandX() const {return lbandx.size();}
  //! Retourne le nombre de bandes selon Y
  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
  //! Retourne le nombre de slices selon X
  inline int     NSliX() const {return lslix.size();}
  //! Retourne le nombre de slices selon Y
  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
  //! structure de definition des bandes
  struct bande_slice {
    int num;   //!< nombre de bandes
    float min; //!< limite minimum pour remplir la bande
    float max; //!< limite maximum pour remplir la bande
    Histo* H;  //!< pointer sur l Histo 1D de la bande
    STRUCTCOMP(bande_slice)
  };
#ifdef __DECCXX
protected:
#endif

  void Delete();

  float*         data; //!< donnees
  double*        err2; //!< erreurs carrees

  float          over[3][3]; //!< overflow table
  double         nHist;      //!< somme ponderee des entrees
  int_4          nEntries;   //!< nombre d'entrees

  int_4          nx;    //!< nombre de bins en X
  int_4          ny;    //!< nombre de bins en Y
  int_4          nxy;   //!< nombre de bins total
  float          xmin;  //!< abscisse minimum
  float          xmax;  //!< abscisse maximum
  float          ymin;  //!< ordonnee minimum
  float          ymax;  //!< ordonnee maximum
  float          wbinx; //!< largeur du bin en X
  float          wbiny; //!< largeur du bin en Y

  bande_slice    b_s;

  Histo*         hprojx; //!< pointer sur Histo des proj X
  Histo*         hprojy; //!< pointer sur Histo des proj Y

  list<bande_slice>  lbandx; //!< liste des bandes selon X
  list<bande_slice>  lbandy; //!< liste des bandes selon Y
  
  list<bande_slice>  lslix; //!< liste des slices selon X
  list<bande_slice>  lsliy; //!< liste des slices selon Y

};

/////////////////////////////////////////////////////////////////////////
// Classe pour la gestion de persistance

/*! \ingroup HiStats \fn operator<<(POuttPersist&,Histo2D)
  \brief Persistance management */
inline POutPersist& operator << (POutPersist& os, Histo2D & obj)
{ ObjFileIO<Histo2D> fio(&obj);  fio.Write(os);  return(os); }
/*! \ingroup HiStats \fn operator<<(POuttPersist&,Histo2D)
  \brief Persistance management */
inline PInPersist& operator >> (PInPersist& is, Histo2D & obj)
{ ObjFileIO<Histo2D> fio(&obj);  fio.Read(is);  return(is); }

// Classe pour la gestion de persistance
// ObjFileIO<Histo2D>

/*! \ingroup HiStats \fn operator*(const Histo2D&,double)
  \brief Operateur H2 = H1 * b */
inline Histo2D operator * (const Histo2D& a, double b)
{
  Histo2D result(a);
  return (result *= b);
}

/*! \ingroup HiStats \fn operator*(double,const Histo2D&)
  \brief Operateur H2 = b * H1 */
inline Histo2D operator * (double b, const Histo2D& a)
{
  Histo2D result(a);
  return (result *= b);
}

/*! \ingroup HiStats \fn operator/(const Histo2D&,double)
  \brief Operateur H2 = H1 / b */
inline Histo2D operator / (const Histo2D& a, double b)
{
  Histo2D result(a);
  return (result /= b);
}

/*! \ingroup HiStats \fn operator+(const Histo2D&,double)
  \brief Operateur H2 = H1 + b */
inline Histo2D operator + (const Histo2D& a, double b)
{
  Histo2D result(a);
  return (result += b);
}

/*! \ingroup HiStats \fn operator+(double,const Histo2D&)
  \brief Operateur H2 = b + H1 */
inline Histo2D operator + (double b, const Histo2D& a)
{
  Histo2D result(a);
  return (result += b);
}

/*! \ingroup HiStats \fn operator-(const Histo2D&,double)
  \brief Operateur H2 = H1 - b */
inline Histo2D operator - (const Histo2D& a, double b)
{
  Histo2D result(a);
  return (result -= b);
}

/*! \ingroup HiStats \fn operator-(double,const Histo2D&)
  \brief Operateur H2 = b - H1 */
inline Histo2D operator - (double b, const Histo2D& a)
{
  Histo2D result(a);
  result *= -1.;
  return (result += b);
}

/*! \ingroup HiStats \fn operator+(const Histo2D&,const Histo2D&)
  \brief Operateur H = H1 + H2 */

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

/*! \ingroup HiStats \fn operator-(const Histo2D&,const Histo2D&)
  \brief Operateur H = H1 - H2 */
inline 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);
}

/*! \ingroup HiStats \fn operator*(const Histo2D&,const Histo2D&)
  \brief Operateur H = H1 * H2 */
inline 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);
}

/*! \ingroup HiStats \fn operator/(const Histo2D&,const Histo2D&)
  \brief Operateur H = H1 / H2 */
inline 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);
}

} // Fin du namespace

#endif // HISTOS2_SEEN