source: Sophya/trunk/SophyaLib/HiStats/hist2err.h@ 3865

#ifndef HIST2ERR_SEEN
#define HIST2ERR_SEEN

#include "objfio.h"
#include <iostream>
#include <stdio.h>
#include "tmatrix.h"

namespace SOPHYA {

//  Forward class declaration for Fits handler
template <class T>  class FitsHandler;

//! 2 dimensions histograms with errors given by user
class Histo2DErr : public AnyDataObj {
  friend class ObjFileIO<Histo2DErr>;
  friend class FitsHandler<Histo2DErr>;
public:

  //! Create or destroy
  Histo2DErr(void);
  Histo2DErr(r_8 xmin,r_8 xmax,int_4 nx,r_8 ymin,r_8 ymax,int_4 ny);
  Histo2DErr(const Histo2DErr& H);
  virtual ~Histo2DErr(void);

  //! Updating or Setting
  void ReCenterBinX(void);
  void ReCenterBinY(void);
  void ReCenterBin(void);
  void ReCenterBinW(void);
  void Zero(void);

  //! Getting Info
  int_4 NBinX(void) {return nx_;}
  int_4 NBinY(void) {return ny_;}
  r_8 XMin(void) {return xmin_;}
  r_8 YMin(void) {return ymin_;}
  r_8 XMax(void) {return xmax_;}
  r_8 YMax(void) {return ymax_;}
  r_8 WBinX(void) {return dx_;}
  r_8 WBinY(void) {return dy_;}

  //! Retourne le contenu du bin
  inline r_8 operator()(int_4 i,int_4 j) const
    {
    if(i<0 || i>=nx_ || j<0 || j>=ny_) return 0.;
    return data_(i,j);
    }
  inline r_8& operator()(int_4 i,int_4 j)
    {
    return data_(i,j);
    }

  //! Retourne le nombre l'erreur au carre dans le bin
  inline r_8 Error2(int_4 i,int_4 j) const
    {
    if(i<0 || i>=nx_ || j<0 || j>=ny_) return 0.;
    return err2_(i,j);
    }
  inline r_8 Error(int_4 i,int_4 j) const
    {
    if(i<0 || i>=nx_ || j<0 || j>=ny_) return 0.;
    if(err2_(i,j)>0.) return sqrt(err2_(i,j)); else return 0.;
    }

  //! Retourne le nombre d'entree dans le bin
  inline r_8 NEntBin(int_4 i,int_4 j) const
    {
    if(i<0 || i>=nx_ || j<0 || j>=ny_) return 0.;
    return ndata_(i,j);
    }

  //! Retourne l'abscisse et l'ordonnee du coin inferieur du bin i,j.
  inline void BinLowEdge(int_4 i,int_4 j,r_8& x,r_8& y) const
              {x = xmin_ + i*dx_; y = ymin_ + j*dy_;}
  //! Retourne l'abscisse et l'ordonnee du centre du bin i,j.
  inline void BinCenter(int_4 i,int_4 j,r_8& x,r_8& y) const
              {x = xmin_ + (i+0.5)*dx_; y = ymin_ + (j+0.5)*dy_;}
  //! Retourne l'abscisse et l'ordonnee du coin superieur du bin i,j.
  inline void BinHighEdge(int_4 i,int_4 j,r_8& x,r_8& y) const
              {x = xmin_ + (i+1)*dx_; y = ymin_ + (j+1)*dy_;}
  //! Retourne les numeros du bin contenant l'abscisse et l'ordonnee x,y.
  inline void FindBin(r_8 x,r_8 y,int_4& i,int_4& j) const
    {i=(int_4) floor((x-xmin_)/dx_); j=(int_4) floor((y-ymin_)/dy_);}

  //! Addition du contenu de l'histo pour abscisse x poids w et l'erreur e2
  inline void Add(r_8 x, r_8 y, r_8 w, r_8 e)
    {
    int_4 i,j; FindBin(x,y,i,j);
    if(i<0 || i>=nx_ || j<0 || j>=ny_) return;
    data_(i,j) += w; ndata_(i,j) += 1.; err2_(i,j) += e*e;
    }
  inline void Add(r_8 x,r_8 y, r_8 w) {Add(x,y,w,w);}

  //! Addition du contenu de l'histo pour le bin numBin poids w et l'erreur e
  inline void AddBin(int_4 i,int_4 j, r_8 w, r_8 e)
    {
    if(i<0 || i>=nx_ || j<0 || j>=ny_) return;
    data_(i,j) += w; ndata_(i,j) += 1.; err2_(i,j) += e*e;
    }
  inline void AddBin(int_4 i,int_4 j, r_8 w) {AddBin(i,j,w,w);}

  //! remplissage contenu de l'histo pour le bin i poids w et l'erreur e
  inline void SetBin(int_4 i,int_4 j, r_8 w, r_8 e, r_8 nb)
    {
    if(i<0 || i>=nx_ || j<0 || j>=ny_) return;
    data_(i,j)  = w;
    err2_(i,j)  = e*e;
    ndata_(i,j) = nb;
    }
  //! remplissage de la valeur pour le bin i
  inline void SetBin(int_4 i,int_4 j, r_8 w)
    {
    if(i<0 || i>=nx_ || j<0 || j>=ny_) return;
    data_(i,j) = w;
    }
  //! remplissage de l'erreur carree pour le bin i
  void SetErr2(int_4 i,int_4 j, r_8 e2)
    {
    if(i<0 || i>=nx_) return;
    err2_(i,j) = e2;
    }
  //! remplissage nombre d'entrees pour le bin i
  void SetNentB(int_4 i,int_4 j, r_8 nb)
    {
    if(i<0 || i>=nx_ || j<0 || j>=ny_) return;
    ndata_(i,j) = nb;
    }

  //! Compute the Mean histogram
  void ToMean(void);
  void FromMean(void);
  int_4 NMean(void) {return mMean;}
  void SetMean(int_4 nmean) {mMean = nmean;}

  //! Replace the errors by the variance
  void ToVariance(void);
  void FromVariance(void);

  //! Fill an histogram with an histogram
  void FillFrHErr(Histo2DErr& hfrom);

  //! Recuperation des matrices elementaires
  void GetData(TMatrix<r_8>& data) {data = data_;}
  void GetError2(TMatrix<r_8>& err2) {err2 = err2_;}
  void GetNData(TMatrix<r_8>& ndata) {ndata = ndata_;}

  // Operators
  Histo2DErr& operator = (const Histo2DErr& h);
  Histo2DErr& operator *= (r_8 b);

  // Print
  void Show(ostream& os) const;
  void Show() const { Show(cout); }

  // Write ASCII
  int WriteASCII(string fname);
  int ReadASCII(string fname);

protected:
  void CreateOrResize(r_8 xmin,r_8 xmax,int_4 nx,r_8 ymin,r_8 ymax,int_4 ny);

  r_4 xmin_,xmax_,dx_,ymin_,ymax_,dy_;
  int_4 nx_,ny_;
  TMatrix<r_8> data_;
  TMatrix<r_8> err2_;
  TMatrix<r_8> ndata_;
  int_4 mMean;  //!< Nombre d'appels a ToMean/Variance(+1) ou FromMean/Variance(-1)
};

/*! Prints histogram information on stream \b s (h.Show(s)) */
inline ostream& operator << (ostream& s, Histo2DErr const & h)
  {  h.Show(s);  return(s);  }

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

} // Fin du namespace

#endif // HIST2ERR_SEEN