| Geant4 Cross Reference |
1 // Copyright (C) 2010, Guy Barrand. All rights
2 // See the file tools.license for terms.
3
4 #ifndef tools_histo_base_histo
5 #define tools_histo_base_histo
6
7 #ifdef TOOLS_MEM
8 #include "../mem"
9 #endif
10
11 #include "histo_data"
12
13 #include <cmath>
14 #include <map> //for annotations
15 #include <ostream>
16
17 namespace tools {
18 namespace histo {
19
20 //TC is for a coordinate.
21 //TO is for an offset used to identify a bin.
22 //TN is for a number of entries.
23 //TW is for a weight.
24 //TH is for a height.
25
26 template <class TC,class TO,class TN,class TW,
27 class base_histo : protected histo_data<TC,TO,
28 typedef histo_data<TC,TO,TN,TW> parent;
29 private:
30 static const std::string& s_class() {
31 static const std::string s_v("tools::histo
32 return s_v;
33 }
34 public:
35 typedef histo_data<TC,TO,TN,TW> hd_t;
36 typedef axis<TC,TO> axis_t;
37 typedef typename axis_t::bn_t bn_t;
38 typedef unsigned int dim_t;
39 typedef TC coordinate_t;
40 typedef TO offset_t;
41 typedef TN num_entries_t;
42 typedef TW weight_t;
43 typedef TH height_t;
44 protected:
45 virtual TH get_bin_height(TO) const = 0; //
46 protected:
47 void base_from_data(const hd_t& a_from) {par
48 #ifdef tools_histo_base_histo //for backward c
49 hd_t base_get_data() const {
50 hd_t hd;
51 hd = *this;
52 return hd;
53 }
54 #endif
55 public:
56 const hd_t& dac() const {return *this;} //da
57 protected:
58 base_histo():parent() {
59 #ifdef TOOLS_MEM
60 mem::increment(s_class().c_str());
61 #endif
62 }
63 protected:
64 virtual ~base_histo() {
65 #ifdef TOOLS_MEM
66 mem::decrement(s_class().c_str());
67 #endif
68 }
69 protected:
70 base_histo(const base_histo& a_from):parent(
71 #ifdef TOOLS_MEM
72 mem::increment(s_class().c_str());
73 #endif
74 }
75
76 base_histo& operator=(const base_histo& a_fr
77 if(&a_from==this) return *this;
78 parent::operator=(a_from);
79 return *this;
80 }
81
82 public:
83 bool equals(const base_histo& a_from,const T
84 return parent::equals(a_from,a_prec,a_fabs
85 }
86
87 const std::string& title() const {return par
88 bool set_title(const std::string& a_title){p
89 dim_t dimension() const {return parent::m_di
90 dim_t number_of_planes() const {return dim_p
91
92 TN entries() const {
93 return parent::m_in_range_entries; //not s
94 }
95 TN all_entries() const {
96 return parent::m_all_entries; //works also
97 }
98
99 TN extra_entries() const {
100 return parent::m_all_entries-parent::m_in_
101 }
102 TW equivalent_bin_entries() const {
103 TW sw = 0;
104 TW sw2 = 0;
105 for(TO ibin=0;ibin<parent::m_bin_number;ib
106 if(!histo::is_out(parent::m_axes,ibin))
107 sw += parent::m_bin_Sw[ibin];
108 sw2 += parent::m_bin_Sw2[ibin];
109 }
110 }
111 if(sw2==0) return 0;
112 return (sw * sw)/sw2;
113 }
114 TH sum_bin_heights() const {
115 TH sh = 0;
116 for(TO ibin=0;ibin<parent::m_bin_number;ib
117 if(!histo::is_out(parent::m_axes,ibin))
118 sh += get_bin_height(ibin);
119 }
120 }
121 return sh;
122 }
123 TH sum_all_bin_heights() const {
124 TH sh = 0;
125 for(TO ibin=0;ibin<parent::m_bin_number;ib
126 sh += get_bin_height(ibin);
127 }
128 return sh;
129 }
130
131 TH sum_extra_bin_heights() const {
132 TH sh = 0;
133 for(TO ibin=0;ibin<parent::m_bin_number;ib
134 if(histo::is_out(parent::m_axes,ibin)) {
135 sh += get_bin_height(ibin);
136 }
137 }
138 return sh;
139 }
140
141 TH min_bin_height() const {
142 TH value = 0;
143 bool first = true;
144 for(TO ibin=0;ibin<parent::m_bin_number;ib
145 if(!histo::is_out(parent::m_axes,ibin))
146 TH vbin = get_bin_height(ibin);
147 if(first) {
148 first = false;
149 value = vbin;
150 } else {
151 if(vbin<=value) value = vbin;
152 }
153 }
154 }
155 return value;
156 }
157
158 TH max_bin_height() const {
159 TH value = 0;
160 bool first = true;
161 for(TO ibin=0;ibin<parent::m_bin_number;ib
162 if(!histo::is_out(parent::m_axes,ibin))
163 TH vbin = get_bin_height(ibin);
164 if(first) {
165 first = false;
166 value = vbin;
167 } else {
168 if(vbin>=value) value = vbin;
169 }
170 }
171 }
172 return value;
173 }
174
175 bool min_bin_height_with_entries(TH& a_value
176 TH value = 0;
177 bool first = true;
178 for(TO ibin=0;ibin<parent::m_bin_number;ib
179 if(!histo::is_out(parent::m_axes,ibin) &
180 TH vbin = get_bin_height(ibin);
181 if(first) {
182 first = false;
183 value = vbin;
184 } else {
185 if(vbin<=value) value = vbin;
186 }
187 }
188 }
189 a_value = value;
190 return first?false:true; //return true if
191 }
192
193 bool max_bin_height_with_entries(TH& a_value
194 TH value = 0;
195 bool first = true;
196 for(TO ibin=0;ibin<parent::m_bin_number;ib
197 if(!histo::is_out(parent::m_axes,ibin) &
198 TH vbin = get_bin_height(ibin);
199 if(first) {
200 first = false;
201 value = vbin;
202 } else {
203 if(vbin>=value) value = vbin;
204 }
205 }
206 }
207 a_value = value;
208 return first?false:true; //return true if
209 }
210
211 bool has_entries_per_bin() const { //to dete
212 // it assumes that update_fast_getters() h
213 if(parent::m_in_range_entries) return true
214 // may be a from-root histo :
215 if(parent::m_in_range_Sw) return false;
216 // no in range entries and weight :
217 return true; //for exa not filled = ok.
218 }
219
220 public: //histo_data
221 bool get_ith_axis_Sxw(dim_t a_axis,TC& a_val
222 a_value = 0;
223 if(a_axis>=parent::m_dimension) return fal
224 for(TO ibin=0;ibin<parent::m_bin_number;ib
225 if(!histo::is_out(parent::m_axes,ibin))
226 a_value += parent::m_bin_Sxw[ibin][a_a
227 }
228 }
229 return true;
230 }
231
232 bool get_ith_axis_Sx2w(dim_t a_axis,TC& a_va
233 a_value = 0;
234 if(a_axis>=parent::m_dimension) return fal
235 for(TO ibin=0;ibin<parent::m_bin_number;ib
236 if(!histo::is_out(parent::m_axes,ibin))
237 a_value += parent::m_bin_Sx2w[ibin][a_
238 }
239 }
240 return true;
241 }
242
243 TW get_in_range_Sw() const {return parent::m
244 TW get_in_range_Sw2() const {return parent::
245
246 void get_Sw_Sw2(TW& a_sw,TW& a_sw2) const {
247 a_sw = 0;
248 a_sw2 = 0;
249 for(TO ibin=0;ibin<parent::m_bin_number;ib
250 if(!histo::is_out(parent::m_axes,ibin))
251 a_sw += parent::m_bin_Sw[ibin];
252 a_sw2 += parent::m_bin_Sw2[ibin];
253 }
254 }
255 }
256 void get_all_Sw_Sw2(TW& a_sw,TW& a_sw2) cons
257 a_sw = 0;
258 a_sw2 = 0;
259 for(TO ibin=0;ibin<parent::m_bin_number;ib
260 a_sw += parent::m_bin_Sw[ibin];
261 a_sw2 += parent::m_bin_Sw2[ibin];
262 }
263 }
264
265 /*
266 TW get_all_Sw() const {
267 TW sw = 0;
268 for(TO ibin=0;ibin<m_bin_number;ibin++) sw
269 return sw;
270 }
271 TN get_all_entries() const {
272 TN number = 0;
273 for(TO ibin=0;ibin<m_bin_number;ibin++) {
274 number += m_bin_entries[ibin];
275 }
276 return number;
277 }
278
279 // for tools/wroot/streamers :
280 TN get_entries() const {
281 TN number = 0;
282 for(TO ibin=0;ibin<m_bin_number;ibin++) {
283 if(!histo::is_out(m_axes,ibin)) {
284 number += m_bin_entries[ibin];
285 }
286 }
287 return number;
288 }
289 */
290 protected:
291 enum {AxisX=0,AxisY=1,AxisZ=2};
292
293 bool configure(dim_t a_dim,
294 const std::vector<bn_t>& aNum
295 const std::vector<TC>& aMins,
296 const std::vector<TC>& aMaxs)
297 // Clear :
298 parent::m_bin_entries.clear();
299 parent::m_bin_Sw.clear();
300 parent::m_bin_Sw2.clear();
301 parent::m_bin_Sxw.clear();
302 parent::m_bin_Sx2w.clear();
303 parent::m_in_range_Sxw.clear();
304 parent::m_in_range_Sx2w.clear();
305 parent::m_axes.clear();
306 parent::m_in_range_plane_Sxyw.clear();
307 parent::m_annotations.clear();
308
309 parent::m_bin_number = 0;
310 parent::m_dimension = 0;
311 parent::m_all_entries = 0;
312 parent::m_in_range_entries = 0;
313 parent::m_in_range_Sw = 0;
314 parent::m_in_range_Sw2 = 0;
315 parent::m_in_range_Sxw.resize(a_dim,0);
316 parent::m_in_range_Sx2w.resize(a_dim,0);
317
318 // Some checks :
319 if(!a_dim) return false;
320 parent::m_axes.resize(a_dim);
321 // Setup axes :
322 for(dim_t iaxis=0;iaxis<a_dim;iaxis++) {
323 if(!parent::m_axes[iaxis].configure(aNum
324 // do not do :
325 // m_axes.clear()
326 // so that :
327 // b1::axis(),b2::axis_[x,y]()
328 // do not crash in case of a bad booki
329 //m_axes.clear();
330 return false;
331 }
332 }
333
334 parent::m_dimension = a_dim;
335
336 base_allocate(); //set m_bin_number.
337
338 return true;
339 }
340
341 bool configure(dim_t a_dim,const std::vector
342 // Clear :
343 parent::m_bin_entries.clear();
344 parent::m_bin_Sw.clear();
345 parent::m_bin_Sw2.clear();
346 parent::m_bin_Sxw.clear();
347 parent::m_bin_Sx2w.clear();
348 parent::m_in_range_Sxw.clear();
349 parent::m_in_range_Sx2w.clear();
350 parent::m_axes.clear();
351 parent::m_in_range_plane_Sxyw.clear();
352 parent::m_annotations.clear();
353
354 parent::m_bin_number = 0;
355 parent::m_dimension = 0;
356 parent::m_all_entries = 0;
357 parent::m_in_range_entries = 0;
358 parent::m_in_range_Sw = 0;
359 parent::m_in_range_Sw2 = 0;
360 parent::m_in_range_Sxw.resize(a_dim,0);
361 parent::m_in_range_Sx2w.resize(a_dim,0);
362
363 // Some checks :
364 if(!a_dim) return false;
365 parent::m_axes.resize(a_dim);
366 // Setup axes :
367 for(dim_t iaxis=0;iaxis<a_dim;iaxis++) {
368 if(!parent::m_axes[iaxis].configure(a_ed
369 //m_axes.clear();
370 return false;
371 }
372 }
373
374 parent::m_dimension = a_dim;
375
376 base_allocate(); //set m_bin_number.
377
378 return true;
379 }
380
381 void base_reset() {
382 // Reset content (different of clear that
383 for(TO ibin=0;ibin<parent::m_bin_number;ib
384 parent::m_bin_entries[ibin] = 0;
385 parent::m_bin_Sw[ibin] = 0;
386 parent::m_bin_Sw2[ibin] = 0;
387 for(dim_t iaxis=0;iaxis<parent::m_dimens
388 parent::m_bin_Sxw[ibin][iaxis] = 0;
389 parent::m_bin_Sx2w[ibin][iaxis] = 0;
390 }
391 }
392 parent::m_in_range_plane_Sxyw.assign(dim_p
393 //profile not done here.
394 parent::reset_fast_getters();
395 }
396
397 protected:
398 void base_allocate() {
399 dim_t iaxis;
400 // Add two bins for the [under,out]flow da
401 TO n_bin = 1;
402 for(iaxis=0;iaxis<parent::m_dimension;iaxi
403 n_bin *= (parent::m_axes[iaxis].bins() +
404 }
405
406 parent::m_bin_entries.resize(n_bin,0);
407 parent::m_bin_Sw.resize(n_bin,0);
408 parent::m_bin_Sw2.resize(n_bin,0);
409
410 std::vector<TC> empty;
411 empty.resize(parent::m_dimension,0);
412 parent::m_bin_Sxw.resize(n_bin,empty);
413 parent::m_bin_Sx2w.resize(n_bin,empty);
414
415 parent::m_bin_number = n_bin; // All bins
416
417 parent::m_axes[0].m_offset = 1;
418 for(iaxis=1;iaxis<parent::m_dimension;iaxi
419 parent::m_axes[iaxis].m_offset = parent:
420 }
421
422 parent::m_in_range_plane_Sxyw.resize(dim_p
423 }
424
425 public:
426 // to access data from methods :
427 const std::vector<TN>& bins_entries() const
428 const std::vector<TW>& bins_sum_w() const {r
429 const std::vector<TW>& bins_sum_w2() const {
430 const std::vector< std::vector<TC> >& bins_s
431 const std::vector< std::vector<TC> >& bins_s
432 const std::vector<TC>& in_range_planes_xyw()
433
434 public:
435 const axis_t& get_axis(int a_index) const {r
436 offset_t get_bins() const {return parent::m_
437 const std::string& get_title() const {return
438 dim_t get_dimension() const {return parent::
439 bool is_valid() const {return (parent::m_dim
440
441 public: //annotations :
442 typedef std::map<std::string,std::string> an
443 const annotations_t& annotations() const {re
444 annotations_t annotations() {return parent::
445
446 void add_annotation(const std::string& a_key
447 parent::m_annotations[a_key] = a_value; //
448 }
449 bool annotation(const std::string& a_key,std
450 annotations_t::const_iterator it = parent:
451 if(it==parent::m_annotations.end()) {a_val
452 a_value = (*it).second;
453 return true;
454 }
455
456 void set_annotations(const annotations_t& a_
457
458 void hprint_annotations(std::ostream& a_out)
459 a_out << " * ANNOTATIONS :" << std::endl;
460 annotations_t::const_iterator it;
461 for(it=parent::m_annotations.begin();it!=p
462 //out << " * (" << index << ") "
463 a_out << " * " << (*it).first << " = "
464 }
465 }
466
467 protected:
468 bool is_compatible(const base_histo& a_histo
469 if(parent::m_dimension!=a_histo.m_dimensio
470 for(dim_t iaxis=0;iaxis<parent::m_dimensio
471 if(!parent::m_axes[iaxis].is_compatible(
472 }
473 return true;
474 }
475
476 void base_add(const base_histo& a_histo){
477 // The only histogram operation that makes
478 for(TO ibin=0;ibin<parent::m_bin_number;ib
479 parent::m_bin_entries[ibin] += a_histo.m
480 parent::m_bin_Sw[ibin] += a_histo.m_bin_
481 parent::m_bin_Sw2[ibin] += a_histo.m_bin
482 for(dim_t iaxis=0;iaxis<parent::m_dimens
483 parent::m_bin_Sxw[ibin][iaxis] += a_hi
484 parent::m_bin_Sx2w[ibin][iaxis] += a_h
485 }
486 }
487 {size_t nplane = parent::m_in_range_plane_S
488 for(size_t iplane=0;iplane<nplane;iplane++
489 parent::m_in_range_plane_Sxyw[iplane] +=
490 parent::update_fast_getters();
491 }
492
493 void base_subtract(const base_histo& a_histo
494 //ill-defined operation. We keep that beca
495 // We build a new histo with one entry in
496 for(TO ibin=0;ibin<parent::m_bin_number;ib
497 parent::m_bin_entries[ibin] = 1;
498
499 parent::m_bin_Sw[ibin] -= a_histo.m_bin_
500 // Yes, it is a += in the below.
501 parent::m_bin_Sw2[ibin] += a_histo.m_bin
502 for(dim_t iaxis=0;iaxis<parent::m_dimens
503 parent::m_bin_Sxw[ibin][iaxis] -= a_hi
504 parent::m_bin_Sx2w[ibin][iaxis] -= a_h
505 }
506 }
507
508 //{for(dim_t iplane=0;iplane<nplane;iplane++)
509
510 parent::update_fast_getters();
511 }
512
513 bool base_multiply(const base_histo& a_histo
514 //ill-defined operation. We keep that beca
515
516 // We build a new histo with one entry in
517 // this.w * a_histo.w
518 // The current histo is overriden with thi
519 // The m_bin_Sw2 computation is consistent
520
521 if(!is_compatible(a_histo)) return false;
522
523 std::vector<int> is(parent::m_dimension);
524 for(TO ibin=0;ibin<parent::m_bin_number;ib
525 TW swa = parent::m_bin_Sw[ibin];
526 TW sw2a = parent::m_bin_Sw2[ibin];
527 TW swb = a_histo.m_bin_Sw[ibin];
528 TW sw2b = a_histo.m_bin_Sw2[ibin];
529 TW sw = swa * swb;
530 parent::m_bin_entries[ibin] = 1;
531 parent::m_bin_Sw[ibin] = sw;
532 parent::m_bin_Sw2[ibin] = sw2a * swb * s
533 histo::get_indices(parent::m_axes,ibin,i
534 for(dim_t iaxis=0;iaxis<parent::m_dimens
535 TC x = parent::m_axes[iaxis].bin_cente
536 parent::m_bin_Sxw[ibin][iaxis] = x * s
537 parent::m_bin_Sx2w[ibin][iaxis] = x *
538 }
539 }
540
541 //{for(dim_t iplane=0;iplane<nplane;iplane++)
542
543 parent::update_fast_getters();
544 return true;
545 }
546
547 bool base_divide(const base_histo& a_histo)
548 //ill-defined operation. We keep that beca
549
550 // We build a new histo with one entry in
551 // this.w / a_histo.w
552 // The current histo is overriden with thi
553 // The m_bin_Sw2 computation is consistent
554
555 if(!is_compatible(a_histo)) return false;
556
557 std::vector<int> is(parent::m_dimension);
558 for(TO ibin=0;ibin<parent::m_bin_number;ib
559 histo::get_indices(parent::m_axes,ibin,i
560 TW swa = parent::m_bin_Sw[ibin];
561 TW swb = a_histo.m_bin_Sw[ibin];
562 TW sw2a = parent::m_bin_Sw2[ibin];
563 TW sw2b = a_histo.m_bin_Sw2[ibin];
564 if(swb!=0) {
565 parent::m_bin_entries[ibin] = 1;
566 TW sw = swa / swb;
567 parent::m_bin_Sw[ibin] = sw;
568 TW swb2 = swb * swb;
569 parent::m_bin_Sw2[ibin] = sw2a / swb2
570 for(dim_t iaxis=0;iaxis<parent::m_dime
571 TC x = parent::m_axes[iaxis].bin_cen
572 parent::m_bin_Sxw[ibin][iaxis] = x *
573 parent::m_bin_Sx2w[ibin][iaxis] = x
574 }
575 } else {
576 parent::m_bin_entries[ibin] = 0;
577 parent::m_bin_Sw[ibin] = 0;
578 parent::m_bin_Sw2[ibin] = 0;
579 for(dim_t iaxis=0;iaxis<parent::m_dime
580 parent::m_bin_Sxw[ibin][iaxis] = 0;
581 parent::m_bin_Sx2w[ibin][iaxis] = 0;
582 }
583 }
584 }
585
586 //{for(dim_t iplane=0;iplane<nplane;iplane++)
587
588 parent::update_fast_getters();
589 return true;
590 }
591
592 bool base_multiply(TW a_factor) {
593 if(a_factor<0) return false;
594 TW factor2 = a_factor * a_factor;
595 for(TO ibin=0;ibin<parent::m_bin_number;ib
596 parent::m_bin_Sw[ibin] *= a_factor;
597 parent::m_bin_Sw2[ibin] *= factor2;
598 for(dim_t iaxis=0;iaxis<parent::m_dimens
599 parent::m_bin_Sxw[ibin][iaxis] *= a_fa
600 parent::m_bin_Sx2w[ibin][iaxis] *= a_f
601 }
602 }
603 {size_t nplane = parent::m_in_range_plane_S
604 for(size_t iplane=0;iplane<nplane;iplane++
605 parent::update_fast_getters();
606 return true;
607 }
608
609 bool get_ith_axis_mean(dim_t a_axis,TC& a_va
610 a_value = 0;
611 if(a_axis>=parent::m_dimension) return fal
612 TW sw = 0;
613 TC sxw = 0;
614 for(TO ibin=0;ibin<parent::m_bin_number;ib
615 if(!histo::is_out(parent::m_axes,ibin))
616 sw += parent::m_bin_Sw[ibin];
617 sxw += parent::m_bin_Sxw[ibin][a_axis]
618 }
619 }
620 if(sw==0) return false;
621 a_value = sxw/sw;
622 return true;
623 }
624
625 bool get_ith_axis_rms(dim_t a_axis,TC& a_val
626 a_value = 0;
627 if(a_axis>=parent::m_dimension) return fal
628 TW sw = 0;
629 TC sxw = 0;
630 TC sx2w = 0;
631 for(TO ibin=0;ibin<parent::m_bin_number;ib
632 if(!histo::is_out(parent::m_axes,ibin))
633 sw += parent::m_bin_Sw[ibin];
634 sxw += parent::m_bin_Sxw[ibin][a_axis]
635 sx2w += parent::m_bin_Sx2w[ibin][a_axi
636 }
637 }
638 if(sw==0) return false;
639 TC mean = sxw/sw;
640 a_value = ::sqrt(::fabs((sx2w / sw) - mean
641 return true;
642 }
643
644 TN get_bin_entries(const std::vector<int>& a
645 if(parent::m_bin_number==0) return 0;
646 TO offset;
647 if(!histo::get_offset(parent::m_axes,aIs,o
648 return parent::m_bin_entries[offset];
649 }
650 };
651
652 // predefined annotation keys :
653 inline const std::string& key_axis_x_title() {
654 static const std::string s_v("axis_x.title")
655 return s_v;
656 }
657 inline const std::string& key_axis_y_title() {
658 static const std::string s_v("axis_y.title")
659 return s_v;
660 }
661 inline const std::string& key_axis_z_title() {
662 static const std::string s_v("axis_z.title")
663 return s_v;
664 }
665
666 }}
667
668 #endif