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