| Geant4 Cross Reference |
1 // Copyright (C) 2010, Guy Barrand. All rights reserved.
2 // See the file tools.license for terms.
3
4 #ifndef tools_histo_c1d
5 #define tools_histo_c1d
6
7 #include "base_cloud"
8
9 #include "../mnmx"
10
11 #include "h1d"
12
13 namespace tools {
14 namespace histo {
15
16 class c1d : public base_cloud {
17 public:
18 static const std::string& s_class() {
19 static const std::string s_v("tools::histo::c1d");
20 return s_v;
21 }
22 public:
23 bool set_title(const std::string& a_title){
24 m_title = a_title;
25 if(m_histo) m_histo->set_title(a_title);
26 return true;
27 }
28
29 unsigned int dimension() const {return 1;}
30 bool reset() {
31 clear();
32 delete m_histo;
33 m_histo = 0;
34 return true;
35 }
36 unsigned int entries() const {
37 return m_histo ? m_histo->all_entries() : (unsigned int)m_ws.size();
38 }
39 public:
40 double sum_of_weights() const {
41 return (m_histo ? m_histo->sum_bin_heights() : m_Sw);
42 }
43
44 bool convert_to_histogram(){
45 if( (m_cnv_x_num<=0) || (m_cnv_x_max<=m_cnv_x_min) ) {
46 // Cloud min, max should be included in the histo.
47 double dx = 0.01 * (upper_edge() - lower_edge())/BINS();
48 return convert(BINS(),lower_edge(),upper_edge() + dx);
49 } else {
50 return convert(m_cnv_x_num,m_cnv_x_min,m_cnv_x_max);
51 }
52 }
53
54 bool is_converted() const {return m_histo ? true : false;}
55 bool scale(double a_scale) {
56 if(m_histo) {
57 return m_histo->scale(a_scale);
58 } else {
59 size_t number = m_ws.size();
60 for(size_t index=0;index<number;index++) m_ws[index] *= a_scale;
61 m_Sw *= a_scale;
62 m_Sxw *= a_scale;
63 m_Sx2w *= a_scale;
64 return true;
65 }
66 }
67
68 bool set_histogram(h1d* a_histo){ //we take ownership of a_histo.
69 reset();
70 m_histo = a_histo;
71 return true;
72 }
73
74 public:
75 bool fill(double aX,double aW = 1){
76 if(!m_histo && (m_limit!=UNLIMITED()) &&
77 ((int)m_xs.size()>=m_limit)){
78 convert_to_histogram();
79 }
80
81 if(m_histo) {
82 return m_histo->fill(aX,aW);
83 } else {
84 if(m_xs.size()) {
85 m_lower_x = mn<double>(aX,m_lower_x);
86 m_upper_x = mx<double>(aX,m_upper_x);
87 } else {
88 m_lower_x = aX;
89 m_upper_x = aX;
90 }
91 m_xs.push_back(aX);
92 m_ws.push_back(aW);
93 m_Sw += aW;
94 double xw = aX * aW;
95 m_Sxw += xw;
96 m_Sx2w += aX * xw;
97 return true;
98 }
99 }
100
101 double lower_edge() const {
102 return (m_histo ? m_histo->axis().lower_edge() : m_lower_x);
103 }
104 double upper_edge() const {
105 return (m_histo ? m_histo->axis().upper_edge() : m_upper_x);
106 }
107 double value(unsigned int a_index) const {return (m_histo ?0:m_xs[a_index]);}
108 double weight(unsigned int a_index) const {return (m_histo ?0:m_ws[a_index]);}
109 double mean() const {
110 return (m_histo ? m_histo->mean() : (m_Sw?m_Sxw/m_Sw:0));
111 }
112 double rms() const {
113 double _rms = 0; //FIXME nan.
114 if(m_histo) {
115 _rms = m_histo->rms();
116 } else {
117 if(m_Sw==0) {
118 } else {
119 double _mean = m_Sxw / m_Sw;
120 _rms = ::sqrt(::fabs( (m_Sx2w / m_Sw) - _mean * _mean));
121 }
122 }
123 return _rms;
124 }
125
126 bool convert(unsigned int a_bins,double a_lower_edge,double a_upper_edge){
127 if(m_histo) return true;
128 m_histo = new histo::h1d(base_cloud::title(),a_bins,a_lower_edge,a_upper_edge);
129 if(!m_histo) return false;
130 bool status = fill_histogram(*m_histo);
131 clear();
132 return status;
133 }
134
135 bool convert(const std::vector<double>& a_edges) {
136 if(m_histo) return true;
137 m_histo = new histo::h1d(base_cloud::title(),a_edges);
138 if(!m_histo) return false;
139 bool status = fill_histogram(*m_histo);
140 clear();
141 return status;
142 }
143
144 const histo::h1d& histogram() const {
145 if(!m_histo) const_cast<c1d&>(*this).convert_to_histogram();
146 return *m_histo;
147 }
148 template <class HISTO>
149 bool fill_histogram(HISTO& a_histo) const {
150 size_t number = m_xs.size();
151 for(size_t index=0;index<number;index++) {
152 if(!a_histo.fill(m_xs[index],m_ws[index])) return false;
153 }
154 return true;
155 }
156 bool set_conversion_parameters(unsigned int aCnvXnumber,
157 double aCnvXmin,double aCnvXmax){
158 m_cnv_x_num = aCnvXnumber;
159 m_cnv_x_min = aCnvXmin;
160 m_cnv_x_max = aCnvXmax;
161 return true;
162 }
163
164 public:
165 c1d()
166 :base_cloud(UNLIMITED())
167 ,m_lower_x(0),m_upper_x(0)
168 ,m_Sxw(0),m_Sx2w(0)
169 ,m_cnv_x_num(0),m_cnv_x_min(0),m_cnv_x_max(0),m_histo(0)
170 {}
171
172 c1d(const std::string& a_title,int aLimit = base_cloud::UNLIMITED())
173 :base_cloud(aLimit)
174 ,m_lower_x(0),m_upper_x(0)
175 ,m_Sxw(0),m_Sx2w(0)
176 ,m_cnv_x_num(0),m_cnv_x_min(0),m_cnv_x_max(0),m_histo(0)
177 {
178 set_title(a_title);
179 }
180
181 virtual ~c1d(){delete m_histo;}
182 public:
183 c1d(const c1d& a_from)
184 :base_cloud(a_from)
185 ,m_xs(a_from.m_xs)
186 ,m_lower_x(a_from.m_lower_x)
187 ,m_upper_x(a_from.m_upper_x)
188 ,m_Sxw(a_from.m_Sxw)
189 ,m_Sx2w(a_from.m_Sx2w)
190 ,m_cnv_x_num(a_from.m_cnv_x_num)
191 ,m_cnv_x_min(a_from.m_cnv_x_min)
192 ,m_cnv_x_max(a_from.m_cnv_x_max)
193 ,m_histo(0)
194 {
195 if(a_from.m_histo) {
196 m_histo = new histo::h1d(*a_from.m_histo);
197 }
198 }
199
200 c1d& operator=(const c1d& a_from){
201 base_cloud::operator=(a_from);
202 if(&a_from==this) return *this;
203 m_xs = a_from.m_xs;
204 m_lower_x = a_from.m_lower_x;
205 m_upper_x = a_from.m_upper_x;
206 m_Sxw = a_from.m_Sxw;
207 m_Sx2w = a_from.m_Sx2w;
208 m_cnv_x_num = a_from.m_cnv_x_num;
209 m_cnv_x_min = a_from.m_cnv_x_min;
210 m_cnv_x_max = a_from.m_cnv_x_max;
211 delete m_histo;
212 m_histo = 0;
213 if(a_from.m_histo) {
214 m_histo = new histo::h1d(*a_from.m_histo);
215 }
216 return *this;
217 }
218
219 public: //AIDA API
220 double lowerEdge() const {return lower_edge();}
221 double upperEdge() const {return upper_edge();}
222 template <class HISTO>
223 bool fillHistogram(HISTO& a_histo) const {return fill_histogram<HISTO>(a_histo);}
224 protected:
225 void clear(){
226 m_lower_x = 0;
227 m_upper_x = 0;
228 m_Sw = 0;
229 m_Sxw = 0;
230 m_Sx2w = 0;
231 m_xs.clear();
232 m_ws.clear();
233 }
234
235 protected:
236 std::vector<double> m_xs;
237 double m_lower_x;
238 double m_upper_x;
239 double m_Sxw;
240 double m_Sx2w;
241 //
242 unsigned int m_cnv_x_num;
243 double m_cnv_x_min;
244 double m_cnv_x_max;
245 histo::h1d* m_histo;
246 };
247
248 }}
249
250 #endif