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

// This may look like C code, but it is really -*- C++ -*-
//
// $Id: histos.h,v 1.9 2000-07-07 09:44:16 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();
  virtual 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