| 1 | // This may look like C code, but it is really -*- C++ -*-
 | 
|---|
| 2 | //
 | 
|---|
| 3 | // $Id: histos.h,v 1.20 2004-10-26 16:26:39 cmv Exp $
 | 
|---|
| 4 | //
 | 
|---|
| 5 | 
 | 
|---|
| 6 | #ifndef HISTOS_SEEN
 | 
|---|
| 7 | #define HISTOS_SEEN
 | 
|---|
| 8 | 
 | 
|---|
| 9 | #include "objfio.h"
 | 
|---|
| 10 | #include <iostream>
 | 
|---|
| 11 | #include <stdio.h>
 | 
|---|
| 12 | #include "peida.h"
 | 
|---|
| 13 | #include "tvector.h"
 | 
|---|
| 14 | #include "ppersist.h"
 | 
|---|
| 15 | #include "anydataobj.h"
 | 
|---|
| 16 | 
 | 
|---|
| 17 | namespace SOPHYA {
 | 
|---|
| 18 | 
 | 
|---|
| 19 | //! 1 dimension histograms
 | 
|---|
| 20 | class Histo : public AnyDataObj {
 | 
|---|
| 21 |   friend class ObjFileIO<Histo>;
 | 
|---|
| 22 | public:
 | 
|---|
| 23 | 
 | 
|---|
| 24 |   // CREATOR / DESTRUCTOR
 | 
|---|
| 25 |   Histo();
 | 
|---|
| 26 |   Histo(r_8 xMin, r_8 xMax, int_4 nBin=100);
 | 
|---|
| 27 |   Histo(r_4 xMin, r_4 xMax, int_4 nBin=100);
 | 
|---|
| 28 |   Histo(const Histo& H);
 | 
|---|
| 29 |   virtual ~Histo();
 | 
|---|
| 30 | 
 | 
|---|
| 31 |   // OPTIONS
 | 
|---|
| 32 |   void Errors();
 | 
|---|
| 33 | 
 | 
|---|
| 34 |   // UPDATING or SETTING
 | 
|---|
| 35 |   void Zero();
 | 
|---|
| 36 |   void Add(r_8 x, r_8 w);
 | 
|---|
| 37 |   inline void Add(r_8 x) {Add(x,1.);}
 | 
|---|
| 38 |   void AddBin(int_4 numBin, r_8 w);
 | 
|---|
| 39 |   inline void AddBin(int_4 numBin) {AddBin(numBin,1.);}
 | 
|---|
| 40 |   void SetBin(r_8 x, r_8 w);
 | 
|---|
| 41 |   inline void SetBin(r_8 x) {SetBin(x,1.);}
 | 
|---|
| 42 |   void SetBin(int_4 numBin, r_8 w);
 | 
|---|
| 43 |   inline void SetBin(int_4 numBin) {SetBin(numBin,1.);}
 | 
|---|
| 44 |   void SetErr2(r_8 x, r_8 e2);
 | 
|---|
| 45 |   void SetErr2(int_4 numBin, r_8 e2);
 | 
|---|
| 46 |   void SetErr(r_8 x, r_8 e);
 | 
|---|
| 47 |   void SetErr(int_4 numBin, r_8 e);
 | 
|---|
| 48 |   virtual void UpdateHisto(bool force=false) const;
 | 
|---|
| 49 | 
 | 
|---|
| 50 |   // Operators
 | 
|---|
| 51 |   Histo& operator = (const Histo& h);
 | 
|---|
| 52 |   Histo& operator *= (r_8 b);
 | 
|---|
| 53 |   Histo& operator /= (r_8 b);
 | 
|---|
| 54 |   Histo& operator += (r_8 b);
 | 
|---|
| 55 |   Histo& operator -= (r_8 b);
 | 
|---|
| 56 |   Histo& operator += (const Histo& a);
 | 
|---|
| 57 |   Histo& operator -= (const Histo& a);
 | 
|---|
| 58 |   Histo& operator *= (const Histo& a);
 | 
|---|
| 59 |   Histo& operator /= (const Histo& a);
 | 
|---|
| 60 | 
 | 
|---|
| 61 |   // get/put dans/depuis un vector
 | 
|---|
| 62 |   void GetAbsc(TVector<r_8>& v) const;
 | 
|---|
| 63 |   void GetValue(TVector<r_8>& v) const;
 | 
|---|
| 64 |   void GetError2(TVector<r_8>& v) const;
 | 
|---|
| 65 |   void GetError(TVector<r_8>& v) const;
 | 
|---|
| 66 |   void PutValue(TVector<r_8>& v, int_4 ierr=0);
 | 
|---|
| 67 |   void PutValueAdd(TVector<r_8> &v, int_4 ierr=0);
 | 
|---|
| 68 |   void PutError2(TVector<r_8>& v);
 | 
|---|
| 69 |   void PutError2Add(TVector<r_8>& v);
 | 
|---|
| 70 |   void PutError(TVector<r_8>& v);
 | 
|---|
| 71 | 
 | 
|---|
| 72 |   // INLINES
 | 
|---|
| 73 |   //! Retourne l'abscisse minimum
 | 
|---|
| 74 |   inline r_8 XMin() const {return mMin;}
 | 
|---|
| 75 |   //! Retourne l'abscisse maximum
 | 
|---|
| 76 |   inline r_8 XMax() const {return mMax;}
 | 
|---|
| 77 |   //! Retourne le nombre de bins
 | 
|---|
| 78 |   inline int_4 NBins() const {return mBins;}
 | 
|---|
| 79 |   //! Retourne la largeur du bin
 | 
|---|
| 80 |   inline r_8 BinWidth() const {return binWidth;}
 | 
|---|
| 81 |   //! Retourne le pointeur sur le tableaux des contenus
 | 
|---|
| 82 |   inline r_8* Bins() const {return mData;}
 | 
|---|
| 83 |   //! Retourne le contenu du bin i
 | 
|---|
| 84 |   inline r_8 operator()(int_4 i) const  {return mData[i];}
 | 
|---|
| 85 |   //! Remplit le contenu du bin i
 | 
|---|
| 86 |   inline r_8& operator()(int_4 i) {return mData[i];}
 | 
|---|
| 87 |   //! retourne "true" si il y a des erreurs stoquees
 | 
|---|
| 88 |   inline bool HasErrors() const
 | 
|---|
| 89 |          {if(mErr2) return true; else return false;}
 | 
|---|
| 90 |   //! Retourne l'erreur du bin i
 | 
|---|
| 91 |   inline r_8 Error(int_4 i) const
 | 
|---|
| 92 |          {if(mErr2) {if(mErr2[i]>0.) return sqrt(mErr2[i]); else return 0.;}
 | 
|---|
| 93 |                     else return 0.;}
 | 
|---|
| 94 |   //! Retourne l'erreur au carre du bin i
 | 
|---|
| 95 |   inline r_8 Error2(int_4 i) const
 | 
|---|
| 96 |          {if(mErr2) return mErr2[i]; else return 0.;}
 | 
|---|
| 97 |   //! Remplit l'erreur au carre du bin i
 | 
|---|
| 98 |   inline r_8& Error2(int_4 i) {return mErr2[i];}
 | 
|---|
| 99 |   //! Retourne la somme ponderee
 | 
|---|
| 100 |   inline r_8 NData() const            {return nHist;}
 | 
|---|
| 101 |   //! Retourne le nombre d'entrees
 | 
|---|
| 102 |   inline uint_8 NEntries() const         {return nEntries;}
 | 
|---|
| 103 |   //! Retourne le nombre d'overflow
 | 
|---|
| 104 |   inline r_8 NOver() const            {return mOver;}
 | 
|---|
| 105 |   //! Retourne le nombre d'underflow
 | 
|---|
| 106 |   inline r_8 NUnder() const           {return mUnder;}
 | 
|---|
| 107 | 
 | 
|---|
| 108 |   //! Retourne l'abscisse du bord inferieur du bin i
 | 
|---|
| 109 |   inline r_8 BinLowEdge(int_4 i)  const {return mMin + i*binWidth;}
 | 
|---|
| 110 |   //! Retourne l'abscisse du centre du bin i
 | 
|---|
| 111 |   inline r_8 BinCenter(int_4 i)   const {return mMin + (i+0.5)*binWidth;}
 | 
|---|
| 112 |   //! Retourne l'abscisse du bord superieur du bin i
 | 
|---|
| 113 |   inline r_8 BinHighEdge(int_4 i) const {return mMin + (i+1)*binWidth;} 
 | 
|---|
| 114 |   //! Retourne le numero du bin contenant l'abscisse x
 | 
|---|
| 115 |   inline int_4 FindBin(r_8 x) const   
 | 
|---|
| 116 |          {return (int_4) floor((x - mMin) / binWidth);}
 | 
|---|
| 117 | 
 | 
|---|
| 118 |   // Info, statistique et calculs sur les histogrammes
 | 
|---|
| 119 |   int_4     BinNonNul() const;
 | 
|---|
| 120 |   int_4     ErrNonNul() const;
 | 
|---|
| 121 |   int_4     IMax() const;
 | 
|---|
| 122 |   int_4     IMin() const;
 | 
|---|
| 123 |   r_8       VMax() const;
 | 
|---|
| 124 |   r_8       VMin() const;
 | 
|---|
| 125 |   r_8       Mean() const;
 | 
|---|
| 126 |   r_8       Sigma() const;
 | 
|---|
| 127 |   r_8       MeanLH(int_4 il,int_4 ih) const;
 | 
|---|
| 128 |   r_8       SigmaLH(int_4 il,int_4 ih) const;
 | 
|---|
| 129 |   r_8       Mean(r_8 x0, r_8 dx) const;
 | 
|---|
| 130 |   r_8       Sigma(r_8 x0, r_8 dx) const;
 | 
|---|
| 131 |   r_8       CleanedMean() const;
 | 
|---|
| 132 |   r_8       CleanedMean(r_8& sigma) const;
 | 
|---|
| 133 |   int_4     BinPercent(r_8 per) const;
 | 
|---|
| 134 |   int_4     BinPercent(r_8 x,r_8 per,int_4& imin,int_4& imax) const;
 | 
|---|
| 135 |   int_4     BinPercent(r_8 x,r_8 per,r_8& xmin,r_8& xmax) const;
 | 
|---|
| 136 |   void      HInteg(r_8 norm = 0.);
 | 
|---|
| 137 |   void      HDeriv();
 | 
|---|
| 138 |   virtual void HRebin(int_4 nbinew);
 | 
|---|
| 139 | 
 | 
|---|
| 140 |   int_4     MaxiLocal(r_8& maxi,int_4& imax,r_8& maxn,int_4& imaxn) const;
 | 
|---|
| 141 |   r_8       FitMax(int_4 degree=2, r_8 frac=0.5f, int_4 debug=0) const;
 | 
|---|
| 142 |   r_8       FindWidth(r_8 xmax,r_8 frac=0.5f, int_4 debug=0) const;
 | 
|---|
| 143 |   r_8       FindWidth(r_8 frac=0.5f, int_4 debug=0) const;
 | 
|---|
| 144 |   int_4     EstimeMax(r_8& xm,int_4 SzPav = 3) const;
 | 
|---|
| 145 |   int_4     EstimeMax(int_4& im,r_8& xm,int_4 SzPav = 3) const;
 | 
|---|
| 146 |   void      EstimeWidthS(r_8 frac,r_8& widthG,r_8& widthD) const;
 | 
|---|
| 147 | 
 | 
|---|
| 148 |   // Print et Display ASCII
 | 
|---|
| 149 |   virtual void Show(ostream& os) const;
 | 
|---|
| 150 |   inline  void Show() const { Show(cout); }
 | 
|---|
| 151 |   void Print(int_4 dyn = 100, r_8 hmin = 1., r_8 hmax = -1.,
 | 
|---|
| 152 |              int_4 pflag = 0, int_4 il = 1, int_4 ih = -1) const;
 | 
|---|
| 153 |   
 | 
|---|
| 154 | protected:
 | 
|---|
| 155 |   void Delete();
 | 
|---|
| 156 | 
 | 
|---|
| 157 |   r_8*   mData;    //!< donnees
 | 
|---|
| 158 |   r_8*   mErr2;    //!< erreurs carrees
 | 
|---|
| 159 |   r_8    mUnder;   //!< underflow
 | 
|---|
| 160 |   r_8    mOver;    //!< overflow
 | 
|---|
| 161 |   r_8    nHist;    //!< somme ponderee des entrees
 | 
|---|
| 162 |   uint_8 nEntries; //!< nombre d'entrees
 | 
|---|
| 163 |   int_4  mBins;    //!< nombre de bins
 | 
|---|
| 164 |   r_8    mMin;     //!< abscisse minimum
 | 
|---|
| 165 |   r_8    mMax;     //!< abscisse maximum
 | 
|---|
| 166 |   r_8    binWidth; //!< largeur du bin
 | 
|---|
| 167 | };
 | 
|---|
| 168 | 
 | 
|---|
| 169 | /*! Prints histogram information on stream \b s (h.Show(s)) */
 | 
|---|
| 170 | inline ostream& operator << (ostream& s, Histo const & h)
 | 
|---|
| 171 |   {  h.Show(s);  return(s);  }
 | 
|---|
| 172 | 
 | 
|---|
| 173 | /*! \ingroup HiStats \fn operator<<(POuttPersist&,Histo)
 | 
|---|
| 174 |   \brief Persistance management */
 | 
|---|
| 175 | inline POutPersist& operator << (POutPersist& os, Histo & obj)
 | 
|---|
| 176 | { ObjFileIO<Histo> fio(&obj);  fio.Write(os);  return(os); }
 | 
|---|
| 177 | /*! \ingroup HiStats \fn operator<<(POuttPersist&,Histo)
 | 
|---|
| 178 |   \brief Persistance management */
 | 
|---|
| 179 | inline PInPersist& operator >> (PInPersist& is, Histo & obj)
 | 
|---|
| 180 | { ObjFileIO<Histo> fio(&obj); is.SkipToNextObject(); fio.Read(is); return(is); }
 | 
|---|
| 181 | 
 | 
|---|
| 182 | // Classe pour la gestion de persistance
 | 
|---|
| 183 | // ObjFileIO<Histo>
 | 
|---|
| 184 | 
 | 
|---|
| 185 | /*! \ingroup HiStats \fn operator*(const Histo&,r_8)
 | 
|---|
| 186 |   \brief Operateur H2 = H1 * b */
 | 
|---|
| 187 | inline Histo operator * (const Histo& a, r_8 b)
 | 
|---|
| 188 | {
 | 
|---|
| 189 |   Histo result(a);
 | 
|---|
| 190 |   return (result *= b);
 | 
|---|
| 191 | }
 | 
|---|
| 192 | 
 | 
|---|
| 193 | /*! \ingroup HiStats \fn operator*(r_8,const Histo&)
 | 
|---|
| 194 |   \brief Operateur H2 = b * H1 */
 | 
|---|
| 195 | inline Histo operator * (r_8 b, const Histo& a)
 | 
|---|
| 196 | {
 | 
|---|
| 197 |   Histo result(a);
 | 
|---|
| 198 |   return (result *= b);
 | 
|---|
| 199 | }
 | 
|---|
| 200 | 
 | 
|---|
| 201 | /*! \ingroup HiStats \fn operator/(const Histo&,r_8)
 | 
|---|
| 202 |   \brief Operateur H2 = H1 / b */
 | 
|---|
| 203 | inline Histo operator / (const Histo& a, r_8 b)
 | 
|---|
| 204 | {
 | 
|---|
| 205 |   Histo result(a);
 | 
|---|
| 206 |   return (result /= b);
 | 
|---|
| 207 | }
 | 
|---|
| 208 | 
 | 
|---|
| 209 | /*! \ingroup HiStats \fn operator+(const Histo&,r_8)
 | 
|---|
| 210 |   \brief Operateur H2 = H1 + b */
 | 
|---|
| 211 | inline Histo operator + (const Histo& a, r_8 b)
 | 
|---|
| 212 | {
 | 
|---|
| 213 |   Histo result(a);
 | 
|---|
| 214 |   return (result += b);
 | 
|---|
| 215 | }
 | 
|---|
| 216 | 
 | 
|---|
| 217 | /*! \ingroup HiStats \fn operator+(r_8,const Histo&)
 | 
|---|
| 218 |   \brief Operateur H2 = b + H1 */
 | 
|---|
| 219 | inline Histo operator + (r_8 b, const Histo& a)
 | 
|---|
| 220 | {
 | 
|---|
| 221 |   Histo result(a);
 | 
|---|
| 222 |   return (result += b);
 | 
|---|
| 223 | }
 | 
|---|
| 224 | 
 | 
|---|
| 225 | /*! \ingroup HiStats \fn operator-(const Histo&,r_8)
 | 
|---|
| 226 |   \brief Operateur H2 = H1 - b */
 | 
|---|
| 227 | inline Histo operator - (const Histo& a, r_8 b)
 | 
|---|
| 228 | {
 | 
|---|
| 229 |   Histo result(a);
 | 
|---|
| 230 |   return (result -= b);
 | 
|---|
| 231 | }
 | 
|---|
| 232 | 
 | 
|---|
| 233 | /*! \ingroup HiStats \fn operator-(r_8,const Histo&)
 | 
|---|
| 234 |   \brief Operateur H2 = b - H1 */
 | 
|---|
| 235 | inline Histo operator - (r_8 b, const Histo& a)
 | 
|---|
| 236 | {
 | 
|---|
| 237 |   Histo result(a);
 | 
|---|
| 238 |   result *= -1.;
 | 
|---|
| 239 |   return (result += b);
 | 
|---|
| 240 | }
 | 
|---|
| 241 | 
 | 
|---|
| 242 | /*! \ingroup HiStats \fn operator+(const Histo&,const Histo&)
 | 
|---|
| 243 |   \brief Operateur H = H1 + H2 */
 | 
|---|
| 244 | inline Histo operator + (const Histo& a, const Histo& b)
 | 
|---|
| 245 | {
 | 
|---|
| 246 | if (b.NBins()!=a.NBins()) throw SzMismatchError(PExcLongMessage(""));
 | 
|---|
| 247 | Histo c(a);
 | 
|---|
| 248 | return (c += b);
 | 
|---|
| 249 | }
 | 
|---|
| 250 | 
 | 
|---|
| 251 | /*! \ingroup HiStats \fn operator-(const Histo&,const Histo&)
 | 
|---|
| 252 |   \brief Operateur H = H1 - H2 */
 | 
|---|
| 253 | inline Histo operator - (const Histo& a, const Histo& b)
 | 
|---|
| 254 | {
 | 
|---|
| 255 | if (b.NBins()!=a.NBins()) throw SzMismatchError(PExcLongMessage(""));
 | 
|---|
| 256 | Histo c(a);
 | 
|---|
| 257 | return (c -= b);
 | 
|---|
| 258 | }
 | 
|---|
| 259 | 
 | 
|---|
| 260 | /*! \ingroup HiStats \fn operator*(const Histo&,const Histo&)
 | 
|---|
| 261 |   \brief Operateur H = H1 * H2 */
 | 
|---|
| 262 | inline Histo operator * (const Histo& a, const Histo& b)
 | 
|---|
| 263 | {
 | 
|---|
| 264 | if (b.NBins()!=a.NBins()) throw SzMismatchError(PExcLongMessage(""));
 | 
|---|
| 265 | Histo c(a);
 | 
|---|
| 266 | return (c *= b);
 | 
|---|
| 267 | }
 | 
|---|
| 268 | 
 | 
|---|
| 269 | /*! \ingroup HiStats \fn operator/(const Histo&,const Histo&)
 | 
|---|
| 270 |   \brief Operateur H = H1 / H2 */
 | 
|---|
| 271 | inline Histo operator / (const Histo& a, const Histo& b)
 | 
|---|
| 272 | {
 | 
|---|
| 273 | if (b.NBins()!=a.NBins()) throw SzMismatchError(PExcLongMessage(""));
 | 
|---|
| 274 | Histo c(a);
 | 
|---|
| 275 | return (c /= b);
 | 
|---|
| 276 | }
 | 
|---|
| 277 | 
 | 
|---|
| 278 | } // Fin du namespace
 | 
|---|
| 279 | 
 | 
|---|
| 280 | #endif // HISTOS_SEEN
 | 
|---|