| 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_b2
5 #define tools_histo_b2
6
7 #include "base_histo"
8
9 #include <ostream>
10
11 namespace tools {
12 namespace histo {
13
14 template <class TC,class TO,class TN,class TW,class TH>
15 class b2 : public base_histo<TC,TO,TN,TW,TH> {
16 typedef base_histo<TC,TO,TN,TW,TH> parent;
17 protected:
18 enum {AxisX=0,AxisY=1};
19 public:
20 typedef base_histo<TC,TO,TN,TW,TH> base_histo_t;
21 typedef typename parent::axis_t axis_t;
22 typedef typename parent::bn_t bn_t;
23 public:
24 virtual TH bin_error(int,int) const = 0; //for print
25 public:
26 // Partition :
27 TC mean_x() const {
28 if(parent::m_in_range_Sw==0) return 0;
29 return parent::m_in_range_Sxw[0]/parent::m_in_range_Sw;
30 }
31
32 TC mean_y() const {
33 if(parent::m_in_range_Sw==0) return 0;
34 return parent::m_in_range_Sxw[1]/parent::m_in_range_Sw;
35 }
36
37 TC rms_x() const {
38 if(parent::m_in_range_Sw==0) return 0;
39 TC mean = parent::m_in_range_Sxw[0]/parent::m_in_range_Sw;
40 return ::sqrt(::fabs((parent::m_in_range_Sx2w[0] / parent::m_in_range_Sw) - mean * mean));
41 }
42
43 TC rms_y() const {
44 if(parent::m_in_range_Sw==0) return 0;
45 TC mean = parent::m_in_range_Sxw[1]/parent::m_in_range_Sw;
46 return ::sqrt(::fabs((parent::m_in_range_Sx2w[1] / parent::m_in_range_Sw) - mean * mean));
47 }
48
49 int coord_to_index_x(TC aCoord) const {
50 return axis_x().coord_to_index(aCoord);
51 }
52 int coord_to_index_y(TC aCoord) const {
53 return axis_y().coord_to_index(aCoord);
54 }
55
56 // bins :
57 TN bin_entries(int aI,int aJ) const {
58 TO offset;
59 if(!_find_offset(aI,aJ,offset)) return 0;
60 return parent::m_bin_entries[offset];
61 }
62
63 TW bin_Sw(int aI,int aJ) const {
64 TO offset;
65 if(!_find_offset(aI,aJ,offset)) return 0;
66 return parent::m_bin_Sw[offset];
67 }
68
69 TW bin_Sw2(int aI,int aJ) const {
70 TO offset;
71 if(!_find_offset(aI,aJ,offset)) return 0;
72 return parent::m_bin_Sw2[offset];
73 }
74 TC bin_Sxw(int aI,int aJ) const {
75 TO offset;
76 if(!_find_offset(aI,aJ,offset)) return 0;
77 return parent::m_bin_Sxw[offset][AxisX];
78 }
79 TC bin_Sx2w(int aI,int aJ) const {
80 TO offset;
81 if(!_find_offset(aI,aJ,offset)) return 0;
82 return parent::m_bin_Sx2w[offset][AxisX];
83 }
84 TC bin_Syw(int aI,int aJ) const {
85 TO offset;
86 if(!_find_offset(aI,aJ,offset)) return 0;
87 return parent::m_bin_Sxw[offset][AxisY];
88 }
89 TC bin_Sy2w(int aI,int aJ) const {
90 TO offset;
91 if(!_find_offset(aI,aJ,offset)) return 0;
92 return parent::m_bin_Sx2w[offset][AxisY];
93 }
94
95 TH bin_height(int aI,int aJ) const {
96 TO offset;
97 if(!_find_offset(aI,aJ,offset)) return 0;
98 return this->get_bin_height(offset);
99 }
100
101 TC bin_center_x(int aI) const {
102 return parent::m_axes[0].bin_center(aI);
103 }
104 TC bin_center_y(int aJ) const {
105 return parent::m_axes[1].bin_center(aJ);
106 }
107
108 TC bin_mean_x(int aI,int aJ) const {
109 TO offset;
110 if(!_find_offset(aI,aJ,offset)) return 0;
111 TW sw = parent::m_bin_Sw[offset];
112 if(sw==0) return 0;
113 return parent::m_bin_Sxw[offset][AxisX]/sw;
114 }
115
116 TC bin_mean_y(int aI,int aJ) const {
117 TO offset;
118 if(!_find_offset(aI,aJ,offset)) return 0;
119 TW sw = parent::m_bin_Sw[offset];
120 if(sw==0) return 0;
121 return parent::m_bin_Sxw[offset][AxisY]/sw;
122 }
123
124 TC bin_rms_x(int aI,int aJ) const {
125 TO offset;
126 if(!_find_offset(aI,aJ,offset)) return 0;
127 TW sw = parent::m_bin_Sw[offset];
128 if(sw==0) return 0;
129 TC sxw = parent::m_bin_Sxw[offset][AxisX];
130 TC sx2w = parent::m_bin_Sx2w[offset][AxisX];
131 TC mean = sxw/sw;
132 return ::sqrt(::fabs((sx2w / sw) - mean * mean));
133 }
134
135 TC bin_rms_y(int aI,int aJ) const {
136 TO offset;
137 if(!_find_offset(aI,aJ,offset)) return 0;
138 TW sw = parent::m_bin_Sw[offset];
139 if(sw==0) return 0;
140 TC sxw = parent::m_bin_Sxw[offset][AxisY];
141 TC sx2w = parent::m_bin_Sx2w[offset][AxisY];
142 TC mean = sxw/sw;
143 return ::sqrt(::fabs((sx2w / sw) - mean * mean));
144 }
145
146 // Axes :
147 const axis_t& axis_x() const {return parent::m_axes[0];}
148 const axis_t& axis_y() const {return parent::m_axes[1];}
149 axis_t& axis_x() {return parent::m_axes[0];} //touchy
150 axis_t& axis_y() {return parent::m_axes[1];} //touchy
151 // Projection :
152
153 TN bin_entries_x(int aI) const {
154 if(!parent::m_dimension) return 0;
155 bn_t ibin;
156 if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) return 0;
157 bn_t ybins = parent::m_axes[1].bins()+2;
158 TO offset;
159 TN _entries = 0;
160 for(bn_t jbin=0;jbin<ybins;jbin++) {
161 offset = ibin + jbin * parent::m_axes[1].m_offset;
162 _entries += parent::m_bin_entries[offset];
163 }
164 return _entries;
165 }
166
167 TW bin_height_x(int aI) const {
168 if(!parent::m_dimension) return 0;
169 bn_t ibin;
170 if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) return 0;
171 bn_t ybins = parent::m_axes[1].bins()+2;
172 TO offset;
173 TW sw = 0;
174 for(bn_t jbin=0;jbin<ybins;jbin++) {
175 offset = ibin + jbin * parent::m_axes[1].m_offset;
176 sw += this->get_bin_height(offset);
177 }
178 return sw;
179 }
180
181 TN bin_entries_y(int aJ) const {
182 if(!parent::m_dimension) return 0;
183 bn_t jbin;
184 if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) return 0;
185 bn_t xbins = parent::m_axes[0].bins()+2;
186 TO offset;
187 TN _entries = 0;
188 for(bn_t ibin=0;ibin<xbins;ibin++) {
189 offset = ibin + jbin * parent::m_axes[1].m_offset;
190 _entries += parent::m_bin_entries[offset];
191 }
192 return _entries;
193 }
194
195 TW bin_height_y(int aJ) const {
196 if(!parent::m_dimension) return 0;
197 bn_t jbin;
198 if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) return 0;
199 bn_t xbins = parent::m_axes[0].bins()+2;
200 TO offset;
201 TW sw = 0;
202 for(bn_t ibin=0;ibin<xbins;ibin++) {
203 offset = ibin + jbin * parent::m_axes[1].m_offset;
204 sw += this->get_bin_height(offset);
205 }
206 return sw;
207 }
208
209 TC Sxyw() const {return parent::m_in_range_plane_Sxyw[0];}
210 public:
211 //NOTE : print is a Python keyword.
212 void hprint(std::ostream& a_out) {
213 // A la HPRINT.
214 a_out << parent::dimension() << parent::title() << std::endl;
215 a_out
216 << " * ENTRIES = " << parent::all_entries() << std::endl;
217
218 // 6 | 7 | 8
219 // -----------
220 // 3 | 4 | 5
221 // -----------
222 // 0 | 1 | 2
223
224 TW height_0 = bin_height(axis_t::UNDERFLOW_BIN,
225 axis_t::UNDERFLOW_BIN);
226 TW height_2 = bin_height(axis_t::OVERFLOW_BIN,
227 axis_t::UNDERFLOW_BIN);
228 TW height_6 = bin_height(axis_t::UNDERFLOW_BIN,
229 axis_t::OVERFLOW_BIN);
230 TW height_8 = bin_height(axis_t::OVERFLOW_BIN,
231 axis_t::OVERFLOW_BIN);
232
233 bn_t i,j;
234
235 TW height_1 = 0;
236 TW height_7 = 0;
237 for(i=0;i<axis_x().bins();i++){
238 height_1 += bin_height(i,axis_t::UNDERFLOW_BIN);
239 height_7 += bin_height(i,axis_t::OVERFLOW_BIN);
240 }
241
242 TW height_3 = 0;
243 TW height_5 = 0;
244 for(j=0;j<axis_y().bins();j++){
245 height_3 += bin_height(axis_t::UNDERFLOW_BIN,j);
246 height_5 += bin_height(axis_t::OVERFLOW_BIN,j);
247 }
248
249 TW height_4 = 0;
250 for(i=0;i<axis_x().bins();i++){
251 for(j=0;j<axis_y().bins();j++){
252 height_4 += bin_height(i,j);
253 }
254 }
255
256 a_out
257 << " " << height_6 << " " << height_7 << " " << height_8 << std::endl;
258 a_out
259 << " " << height_3 << " " << height_4 << " " << height_5 << std::endl;
260 a_out
261 << " " << height_0 << " " << height_1 << " " << height_2 << std::endl;
262
263 // Some bins :
264 bn_t xbins = axis_x().bins();
265 bn_t ybins = axis_y().bins();
266 a_out
267 << " * ENTRIES[0,0] = "
268 << bin_entries(0,0)
269 << " * HEIGHT[0,0] = "
270 << bin_height(0,0)
271 << " * ERROR[0,0] = "
272 << bin_error(0,0)
273 << std::endl;
274 a_out
275 << " * ENTRIES[N/2,N/2] = "
276 << bin_entries(xbins/2,ybins/2)
277 << " * HEIGHT[N/2,N/2] = "
278 << bin_height(xbins/2,ybins/2)
279 << " * ERROR[N/2,N/2] = "
280 << bin_error(xbins/2,ybins/2)
281 << std::endl;
282 a_out
283 << " * ENTRIES[N-1,N-1] = "
284 << bin_entries(xbins-1,ybins-1)
285 << " * HEIGHT[N-1,N-1] = "
286 << bin_height(xbins-1,ybins-1)
287 << " * ERROR[N-1,N-1] = "
288 << bin_error(xbins-1,ybins-1)
289 << std::endl;
290 }
291 protected:
292 b2(const std::string& a_title,bn_t aXnumber,TC aXmin,TC aXmax,bn_t aYnumber,TC aYmin,TC aYmax) {
293 parent::m_title = a_title;
294 std::vector<bn_t> nbins;
295 nbins.push_back(aXnumber);
296 nbins.push_back(aYnumber);
297 std::vector<TC> mins;
298 mins.push_back(aXmin);
299 mins.push_back(aYmin);
300 std::vector<TC> maxs;
301 maxs.push_back(aXmax);
302 maxs.push_back(aYmax);
303 parent::configure(2,nbins,mins,maxs);
304 }
305
306 b2(const std::string& a_title,const std::vector<TC>& a_edges_x,const std::vector<TC>& a_edges_y) {
307 parent::m_title = a_title;
308 std::vector< std::vector<TC> > edges(2);
309 edges[0] = a_edges_x;
310 edges[1] = a_edges_y;
311 parent::configure(2,edges);
312 }
313
314 virtual ~b2(){}
315 protected:
316 b2(const b2& a_from):parent(a_from) {}
317 b2& operator=(const b2& a_from) {parent::operator=(a_from);return *this;}
318 public:
319 bool configure(bn_t aXnumber,TC aXmin,TC aXmax,bn_t aYnumber,TC aYmin,TC aYmax){
320 std::vector<bn_t> nbins;
321 nbins.push_back(aXnumber);
322 nbins.push_back(aYnumber);
323 std::vector<TC> mins;
324 mins.push_back(aXmin);
325 mins.push_back(aYmin);
326 std::vector<TC> maxs;
327 maxs.push_back(aXmax);
328 maxs.push_back(aYmax);
329 return parent::configure(2,nbins,mins,maxs);
330 }
331
332 bool configure(const std::vector<TC>& a_edges_x,const std::vector<TC>& a_edges_y){
333 std::vector< std::vector<TC> > edges(2);
334 edges[0] = a_edges_x;
335 edges[1] = a_edges_y;
336 return parent::configure(2,edges);
337 }
338
339 protected:
340 bool _find_offset(int aI,int aJ,TO& a_offset) const {
341 if(parent::m_dimension!=2) {a_offset=0;return false;}
342 bn_t ibin,jbin;
343 if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) {a_offset=0;return false;}
344 if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) {a_offset=0;return false;}
345 a_offset = ibin + jbin * parent::m_axes[1].m_offset;
346 return true;
347 }
348 };
349
350 }}
351
352 #endif
353
354
355
356