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

// This may look like C code, but it is really -*- C++ -*-
//
// $Id: histos.h,v 1.7 2000-06-30 13:21:53 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;

//! 1 dimension histograms
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);
  Histo& operator += (const Histo& a);
  Histo& operator -= (const Histo& a);
  Histo& operator *= (const Histo& a);
  Histo& operator /= (const Histo& a);

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

  // INLINES
  //! Retourne l'abscisse minimum
  inline float          XMin() const             {return min;}
  //! Retourne l'abscisse maximum
  inline float          XMax() const             {return max;}
  //! Retourne le nombre de bins
  inline int_4          NBins() const            {return bins;}
  //! Retourne la largeur du bin
  inline float          BinWidth() const         {return binWidth;}
  //! Retourne le pointeur sur le tableaux des contenus
  inline float*         Bins() const             {return data;}
  //! Retourne le contenu du bin i
  inline float          operator()(int i) const  {return data[i];}
  //! Remplit le contenu du bin i
  inline float&         operator()(int i)        {return data[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
  inline float          Error(int i) const
                          {if(err2) {if(err2[i]>0.) return sqrt(err2[i]); else return 0.;}
                           else return 0.;}
  //! Retourne l'erreur au carre du bin i
  inline double         Error2(int i) const
                          {if(err2) return err2[i]; else return 0.;}
  //! Remplit l'erreur au carre du bin i
  inline double&        Error2(int i) {return err2[i];}
  //! Retourne la somme ponderee
  inline float          NData() const            {return (float) nHist;}
  //! Retourne le nombre d'entrees
  inline float          NEntries() const         {return nEntries;}
  //! Retourne le nombre d'overflow
  inline float          NOver() const            {return over;}
  //! Retourne le nombre d'underflow
  inline float          NUnder() const           {return under;}

  //! Retourne l'abscisse du bord inferieur du bin i
  inline float          BinLowEdge(int i)  const {return min + i*binWidth;}
  //! Retourne l'abscisse du centre du bin i
  inline float          BinCenter(int i)   const {return min + (i+0.5)*binWidth;}
  //! Retourne l'abscisse du bord superieur du bin i
  inline float          BinHighEdge(int i) const {return min + (i+1)*binWidth;} 
  //! Retourne le numero du bin contenant l'abscisse x
  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;     //!< donnees
  double*        err2;     //!< erreurs carrees
  float          under;    //!< underflow
  float          over;     //!< overflow
  double         nHist;    //!< somme ponderee des entrees
  int_4          nEntries; //!< nombre d'entrees
  int_4          bins;     //!< nombre de bins
  float          min;      //!< abscisse minimum
  float          max;      //!< abscisse maximum
  float          binWidth; //!< largeur du bin
};


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

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

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

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

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

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

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

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

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

/*! \ingroup HiStats \fn operator+(const Histo&,const Histo&)
  \brief Operateur H = H1 + H2 */
inline Histo operator + (const Histo& a, const Histo& b)
{
if (b.NBins()!=a.NBins()) THROW(sizeMismatchErr);
Histo c(a);
return (c += b);
}

/*! \ingroup HiStats \fn operator-(const Histo&,const Histo&)
  \brief Operateur H = H1 - H2 */
inline Histo operator - (const Histo& a, const Histo& b)
{
if (b.NBins()!=a.NBins()) THROW(sizeMismatchErr);
Histo c(a);
return (c -= b);
}

/*! \ingroup HiStats \fn operator*(const Histo&,const Histo&)
  \brief Operateur H = H1 * H2 */
inline Histo operator * (const Histo& a, const Histo& b)
{
if (b.NBins()!=a.NBins()) THROW(sizeMismatchErr);
Histo c(a);
return (c *= b);
}

/*! \ingroup HiStats \fn operator/(const Histo&,const Histo&)
  \brief Operateur H = H1 / H2 */
inline Histo operator / (const Histo& a, const Histo& b)
{
if (b.NBins()!=a.NBins()) THROW(sizeMismatchErr);
Histo c(a);
return (c /= b);
}

} // Fin du namespace

#endif // HISTOS_SEEN