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