| 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_p1
5 #define tools_histo_p1
6
7 #include "b1"
8 #include "profile_data"
9
10 namespace tools {
11 namespace histo {
12
13 //TC is for a coordinate.
14 //TO is for an offset/index to identify/find a bind.
15 //TN is for a number of entries.
16 //TW is for a weight.
17 //TH is for a height. Should be the same as TV.
18 //TV is for a value (in general same as TC).
19
20 template <class TC,class TO,class TN,class TW,class TH,class TV>
21 class p1 : public b1<TC,TO,TN,TW,TH> {
22 typedef b1<TC,TO,TN,TW,TH> parent;
23 public:
24 typedef profile_data<TC,TO,TN,TW,TV> pd_t;
25 typedef typename parent::bn_t bn_t;
26 typedef typename parent::axis_t axis_t;
27 typedef std::vector<TV> vs_t;
28 protected:
29 virtual TH get_bin_height(TO a_offset) const {
30 return (parent::m_bin_Sw[a_offset] ? (m_bin_Svw[a_offset]/parent::m_bin_Sw[a_offset]):0);
31 }
32 public:
33 bool equals(const p1& a_from,const TW& a_prec,TW(*a_fabs)(TW)) const {
34 if(!parent::equals(a_from,a_prec,a_fabs)) return false;
35 if(m_cut_v!=a_from.m_cut_v) return false;
36 if(!numbers_are_equal(m_min_v,a_from.m_min_v,a_prec,a_fabs)) return false;
37 if(!numbers_are_equal(m_max_v,a_from.m_max_v,a_prec,a_fabs)) return false;
38 if(!vectors_are_equal(m_bin_Svw,a_from.m_bin_Svw,a_prec,a_fabs)) return false;
39 if(!vectors_are_equal(m_bin_Sv2w,a_from.m_bin_Sv2w,a_prec,a_fabs)) return false;
40 return true;
41 }
42 bool equals_TH(const p1& a_from,const TW& a_prec,TW(*a_fabs)(TW)) const {
43 if(!parent::equals_TH(a_from,a_prec,a_fabs,false)) return false;
44 if(m_cut_v!=a_from.m_cut_v) return false;
45 if(!numbers_are_equal(m_min_v,a_from.m_min_v,a_prec,a_fabs)) return false;
46 if(!numbers_are_equal(m_max_v,a_from.m_max_v,a_prec,a_fabs)) return false;
47 if(!vectors_are_equal(m_bin_Svw,a_from.m_bin_Svw,a_prec,a_fabs)) return false;
48 if(!vectors_are_equal(m_bin_Sv2w,a_from.m_bin_Sv2w,a_prec,a_fabs)) return false;
49 return true;
50 }
51
52 virtual TH bin_error(int aI) const { //TH should be the same as TV
53 TO offset;
54 if(!parent::_find_offset(aI,offset)) return 0;
55
56 //FIXME Is it correct ?
57 // TProfile::GetBinError with kERRORMEAN mode does :
58 // Stat_t cont = fArray[bin]; //Svw (see TProfile::Fill)
59 // Stat_t sum = parent::m_bin_entries.fArray[bin]; //Sw
60 // Stat_t err2 = fSumw2.fArray[bin]; //Sv2w
61 // if (sum == 0) return 0;
62 // Stat_t eprim;
63 // Stat_t contsum = cont/sum;
64 // Stat_t eprim2 = TMath::Abs(err2/sum - contsum*contsum);
65 // eprim = TMath::Sqrt(eprim2);
66 // ... ???
67 // if (fErrorMode == kERRORMEAN) return eprim/TMath::Sqrt(sum);
68
69 TW sw = parent::m_bin_Sw[offset]; //ROOT sum
70 if(sw==0) return 0;
71 TV svw = m_bin_Svw[offset]; //ROOT cont
72 TV sv2w = m_bin_Sv2w[offset]; //ROOT err2
73 TV _mean = (svw / sw); //ROOT contsum
74 TV _rms = ::sqrt(::fabs((sv2w/sw) - _mean * _mean)); //ROOT eprim
75 // rms = get_bin_rms_value.
76 return _rms/::sqrt(sw); //ROOT kERRORMEAN mode returned value
77 }
78
79 public:
80 bool multiply(TW a_factor){
81 if(!parent::base_multiply(a_factor)) return false;
82 for(bn_t ibin=0;ibin<parent::m_bin_number;ibin++) {
83 m_bin_Svw[ibin] *= a_factor;
84 }
85 return true;
86 }
87 bool scale(TW a_factor) {return multiply(a_factor);}
88
89 TV bin_Svw(int aI) const {
90 TO offset;
91 if(!parent::_find_offset(aI,offset)) return 0;
92 return m_bin_Svw[offset];
93 }
94 TV bin_Sv2w(int aI) const {
95 TO offset;
96 if(!parent::_find_offset(aI,offset)) return 0;
97 return m_bin_Sv2w[offset];
98 }
99
100 bool reset() {
101 parent::base_reset();
102 for(bn_t ibin=0;ibin<parent::m_bin_number;ibin++) {
103 m_bin_Svw[ibin] = 0;
104 m_bin_Sv2w[ibin] = 0;
105 }
106 return true;
107 }
108
109 void copy_from_data(const pd_t& a_from) {
110 parent::base_from_data(a_from);
111 m_bin_Svw = a_from.m_bin_Svw;
112 m_bin_Sv2w = a_from.m_bin_Sv2w;
113 m_cut_v = a_from.m_cut_v;
114 m_min_v = a_from.m_min_v;
115 m_max_v = a_from.m_max_v;
116 }
117 pd_t get_histo_data() const {
118 pd_t hd(parent::dac());
119 hd.m_is_profile = true;
120 hd.m_bin_Svw = m_bin_Svw;
121 hd.m_bin_Sv2w = m_bin_Sv2w;
122 hd.m_cut_v = m_cut_v;
123 hd.m_min_v = m_min_v;
124 hd.m_max_v = m_max_v;
125 return hd;
126 }
127
128 bool fill(TC aX,TV aV,TW aWeight = 1) {
129 if(parent::m_dimension!=1) return false;
130
131 if(m_cut_v) {
132 if( (aV<m_min_v) || (aV>=m_max_v) ) {
133 return true;
134 }
135 }
136
137 bn_t ibin;
138 if(!parent::m_axes[0].coord_to_absolute_index(aX,ibin)) return false;
139 TO offset = ibin;
140
141 parent::m_bin_entries[offset]++;
142 parent::m_bin_Sw[offset] += aWeight;
143 parent::m_bin_Sw2[offset] += aWeight * aWeight;
144
145 TC xw = aX * aWeight;
146 TC x2w = aX * xw;
147 parent::m_bin_Sxw[offset][0] += xw;
148 parent::m_bin_Sx2w[offset][0] += x2w;
149
150 bool inRange = true;
151 if(ibin==0) inRange = false;
152 else if(ibin==(parent::m_axes[0].m_number_of_bins+1)) inRange = false;
153
154 parent::m_all_entries++;
155 if(inRange) {
156 // fast getters :
157 parent::m_in_range_entries++;
158 parent::m_in_range_Sw += aWeight;
159 parent::m_in_range_Sw2 += aWeight*aWeight;
160
161 parent::m_in_range_Sxw[0] += xw;
162 parent::m_in_range_Sx2w[0] += x2w;
163 }
164
165 // Profile part :
166 TV vw = aV * aWeight;
167 m_bin_Svw[offset] += vw;
168 m_bin_Sv2w[offset] += aV * vw;
169
170 return true;
171 }
172
173 bool set_bin_content(bn_t a_ibin,TN a_entries,TW a_Sw,TW a_Sw2,TC a_Sxw,TC a_Sx2w,TC a_Svw,TC a_Sv2w) {
174 if(parent::m_dimension!=1) return false;
175 if(a_ibin>(parent::m_axes[0].m_number_of_bins+1)) return false;
176
177 bool inRange = true;
178 if(a_ibin==0) inRange = false;
179 else if(a_ibin==(parent::m_axes[0].m_number_of_bins+1)) inRange = false;
180
181 TO offset = a_ibin;
182
183 parent::m_all_entries -= parent::m_bin_entries[offset];
184 if(inRange) {
185 parent::m_in_range_entries -= parent::m_bin_entries[offset];
186 parent::m_in_range_Sw -= parent::m_bin_Sw[offset];
187 parent::m_in_range_Sw2 -= parent::m_bin_Sw2[offset];
188 parent::m_in_range_Sxw[0] -= parent::m_bin_Sxw[offset][0];
189 parent::m_in_range_Sx2w[0] -= parent::m_bin_Sx2w[offset][0];
190 }
191
192 parent::m_bin_entries[offset] = a_entries;
193 parent::m_bin_Sw[offset] = a_Sw;
194 parent::m_bin_Sw2[offset] = a_Sw2;
195
196 parent::m_bin_Sxw[offset][0] = a_Sxw;
197 parent::m_bin_Sx2w[offset][0] = a_Sx2w;
198
199 parent::m_all_entries += a_entries;
200 if(inRange) {
201 parent::m_in_range_entries += a_entries;
202 parent::m_in_range_Sw += a_Sw;
203 parent::m_in_range_Sw2 += a_Sw2;
204
205 parent::m_in_range_Sxw[0] += a_Sxw;
206 parent::m_in_range_Sx2w[0] += a_Sx2w;
207 }
208
209 // Profile part :
210 m_bin_Svw[offset] = a_Svw;
211 m_bin_Sv2w[offset] = a_Sv2w;
212
213 return true;
214 }
215
216 bool get_bin_content(bn_t a_ibin,TN& a_entries,TW& a_Sw,TW& a_Sw2,TC& a_Sxw,TC& a_Sx2w,TC& a_Svw,TC& a_Sv2w) {
217 if(parent::m_dimension!=1) {
218 a_entries = 0;a_Sw = 0;a_Sw2 = 0;a_Sxw = 0;a_Sx2w = 0;
219 return false;
220 }
221 if(a_ibin>(parent::m_axes[0].m_number_of_bins+1)) {
222 a_entries = 0;a_Sw = 0;a_Sw2 = 0;a_Sxw = 0;a_Sx2w = 0;
223 return false;
224 }
225
226 TO offset = a_ibin;
227
228 a_entries = parent::m_bin_entries[offset];
229 a_Sw = parent::m_bin_Sw[offset];
230 a_Sw2 = parent::m_bin_Sw2[offset];
231
232 a_Sxw = parent::m_bin_Sxw[offset][0];
233 a_Sx2w = parent::m_bin_Sx2w[offset][0];
234
235 // Profile part :
236 a_Svw = m_bin_Svw[offset];
237 a_Sv2w = m_bin_Sv2w[offset];
238
239 return true;
240 }
241
242 TV bin_rms_value(int aI) const {
243 TO offset;
244 if(!parent::_find_offset(aI,offset)) return 0;
245 TW sw = parent::m_bin_Sw[offset];
246 if(sw==0) return 0;
247 TV svw = m_bin_Svw[offset];
248 TV sv2w = m_bin_Sv2w[offset];
249 TV _mean = (svw / sw);
250 return ::sqrt(::fabs((sv2w / sw) - _mean * _mean));
251 }
252
253 bool add(const p1& a_histo){
254 parent::base_add(a_histo);
255 for(bn_t ibin=0;ibin<parent::m_bin_number;ibin++) {
256 m_bin_Svw[ibin] += a_histo.m_bin_Svw[ibin];
257 m_bin_Sv2w[ibin] += a_histo.m_bin_Sv2w[ibin];
258 }
259 return true;
260 }
261 bool subtract(const p1& a_histo){
262 parent::base_subtract(a_histo);
263 for(bn_t ibin=0;ibin<parent::m_bin_number;ibin++) {
264 m_bin_Svw[ibin] -= a_histo.m_bin_Svw[ibin];
265 m_bin_Sv2w[ibin] -= a_histo.m_bin_Sv2w[ibin];
266 }
267 return true;
268 }
269
270 bool gather_bins(unsigned int a_factor) { //for exa 2,3.
271 if(!a_factor) return false;
272
273 // actual bin number must be a multiple of a_factor.
274
275 const axis_t& _axis = parent::axis();
276
277 bn_t n = _axis.bins();
278 if(!n) return false;
279
280 bn_t new_n = n/a_factor;
281 if(a_factor*new_n!=n) return false;
282
283 p1* new_h = 0;
284 if(_axis.is_fixed_binning()) {
285 new_h = new p1(parent::m_title,new_n,_axis.lower_edge(),_axis.upper_edge());
286 } else {
287 const std::vector<TC>& _edges = _axis.edges();
288 std::vector<TC> new_edges(new_n+1);
289 for(bn_t ibin=0;ibin<new_n;ibin++) {
290 new_edges[ibin] = _edges[ibin*a_factor];
291 }
292 new_edges[new_n] = _edges[n]; //upper edge.
293 new_h = new p1(parent::m_title,new_edges);
294 }
295 if(!new_h) return false;
296
297 new_h->m_cut_v = m_cut_v;
298 new_h->m_min_v = m_min_v;
299 new_h->m_max_v = m_max_v;
300
301 bn_t offset,new_offset,offac;
302 for(bn_t ibin=0;ibin<new_n;ibin++) {
303 new_offset = ibin+1;
304 offset = a_factor*ibin+1;
305 for(unsigned int ifac=0;ifac<a_factor;ifac++) {
306 offac = offset+ifac;
307 new_h->m_bin_entries[new_offset] += parent::m_bin_entries[offac];
308 new_h->m_bin_Sw[new_offset] += parent::m_bin_Sw[offac];
309 new_h->m_bin_Sw2[new_offset] += parent::m_bin_Sw2[offac];
310 new_h->m_bin_Sxw[new_offset][0] += parent::m_bin_Sxw[offac][0];
311 new_h->m_bin_Sx2w[new_offset][0] += parent::m_bin_Sx2w[offac][0];
312
313 new_h->m_bin_Svw[new_offset] += m_bin_Svw[offac];
314 new_h->m_bin_Sv2w[new_offset] += m_bin_Sv2w[offac];
315 }
316 }
317
318 //underflow :
319 new_offset = 0;
320 offac = 0;
321 new_h->m_bin_entries[new_offset] = parent::m_bin_entries[offac];
322 new_h->m_bin_Sw[new_offset] = parent::m_bin_Sw[offac];
323 new_h->m_bin_Sw2[new_offset] = parent::m_bin_Sw2[offac];
324 new_h->m_bin_Sxw[new_offset][0] = parent::m_bin_Sxw[offac][0];
325 new_h->m_bin_Sx2w[new_offset][0] = parent::m_bin_Sx2w[offac][0];
326 new_h->m_bin_Svw[new_offset] = m_bin_Svw[offac];
327 new_h->m_bin_Sv2w[new_offset] = m_bin_Sv2w[offac];
328
329 //overflow :
330 new_offset = new_n+1;
331 offac = n+1;
332 new_h->m_bin_entries[new_offset] = parent::m_bin_entries[offac];
333 new_h->m_bin_Sw[new_offset] = parent::m_bin_Sw[offac];
334 new_h->m_bin_Sw2[new_offset] = parent::m_bin_Sw2[offac];
335 new_h->m_bin_Sxw[new_offset][0] = parent::m_bin_Sxw[offac][0];
336 new_h->m_bin_Sx2w[new_offset][0] = parent::m_bin_Sx2w[offac][0];
337 new_h->m_bin_Svw[new_offset] = m_bin_Svw[offac];
338 new_h->m_bin_Sv2w[new_offset] = m_bin_Sv2w[offac];
339
340 *this = *new_h;
341 return true;
342 }
343
344 bool cut_v() const {return m_cut_v;}
345 TV min_v() const {return m_min_v;}
346 TV max_v() const {return m_max_v;}
347
348 public:
349 p1(const std::string& a_title,bn_t aXnumber,TC aXmin,TC aXmax)
350 :parent(a_title,aXnumber,aXmin,aXmax)
351 ,m_cut_v(false)
352 ,m_min_v(0)
353 ,m_max_v(0)
354 {
355 m_bin_Svw.resize(parent::m_bin_number,0);
356 m_bin_Sv2w.resize(parent::m_bin_number,0);
357 }
358
359 p1(const std::string& a_title,bn_t aXnumber,TC aXmin,TC aXmax,TV aVmin,TV aVmax)
360 :parent(a_title,aXnumber,aXmin,aXmax)
361 ,m_cut_v(true)
362 ,m_min_v(aVmin)
363 ,m_max_v(aVmax)
364 {
365 m_bin_Svw.resize(parent::m_bin_number,0);
366 m_bin_Sv2w.resize(parent::m_bin_number,0);
367 }
368
369 p1(const std::string& a_title,const std::vector<TC>& a_edges)
370 :parent(a_title,a_edges)
371 ,m_cut_v(false)
372 ,m_min_v(0)
373 ,m_max_v(0)
374 {
375 m_bin_Svw.resize(parent::m_bin_number,0);
376 m_bin_Sv2w.resize(parent::m_bin_number,0);
377 }
378
379 p1(const std::string& a_title,const std::vector<TC>& a_edges,TV aVmin,TV aVmax)
380 :parent(a_title,a_edges)
381 ,m_cut_v(true)
382 ,m_min_v(aVmin)
383 ,m_max_v(aVmax)
384 {
385 m_bin_Svw.resize(parent::m_bin_number,0);
386 m_bin_Sv2w.resize(parent::m_bin_number,0);
387 }
388
389 virtual ~p1(){}
390 public:
391 p1(const p1& a_from)
392 :parent(a_from)
393 ,m_cut_v(a_from.m_cut_v)
394 ,m_min_v(a_from.m_min_v)
395 ,m_max_v(a_from.m_max_v)
396 ,m_bin_Svw(a_from.m_bin_Svw)
397 ,m_bin_Sv2w(a_from.m_bin_Sv2w)
398 {}
399 p1& operator=(const p1& a_from){
400 parent::operator=(a_from);
401 m_cut_v = a_from.m_cut_v;
402 m_min_v = a_from.m_min_v;
403 m_max_v = a_from.m_max_v;
404 m_bin_Svw = a_from.m_bin_Svw;
405 m_bin_Sv2w = a_from.m_bin_Sv2w;
406 return *this;
407 }
408 public:
409 bool configure(bn_t aXnumber,TC aXmin,TC aXmax){
410 if(!parent::configure(aXnumber,aXmin,aXmax)) return false;
411 m_bin_Svw.clear();
412 m_bin_Sv2w.clear();
413 m_bin_Svw.resize(parent::m_bin_number,0);
414 m_bin_Sv2w.resize(parent::m_bin_number,0);
415 m_cut_v = false;
416 m_min_v = 0;
417 m_max_v = 0;
418 return true;
419 }
420 bool configure(const std::vector<TC>& a_edges) {
421 if(!parent::configure(a_edges)) return false;
422 m_bin_Svw.clear();
423 m_bin_Sv2w.clear();
424 m_bin_Svw.resize(parent::m_bin_number,0);
425 m_bin_Sv2w.resize(parent::m_bin_number,0);
426 m_cut_v = false;
427 m_min_v = 0;
428 m_max_v = 0;
429 return true;
430 }
431 bool configure(bn_t aXnumber,TC aXmin,TC aXmax,TV aVmin,TV aVmax){
432 if(!parent::configure(aXnumber,aXmin,aXmax)) return false;
433 m_bin_Svw.clear();
434 m_bin_Sv2w.clear();
435 m_bin_Svw.resize(parent::m_bin_number,0);
436 m_bin_Sv2w.resize(parent::m_bin_number,0);
437 m_cut_v = true;
438 m_min_v = aVmin;
439 m_max_v = aVmax;
440 return true;
441 }
442 bool configure(const std::vector<TC>& a_edges,TV aVmin,TV aVmax) {
443 if(!parent::configure(a_edges)) return false;
444 m_bin_Svw.clear();
445 m_bin_Sv2w.clear();
446 m_bin_Svw.resize(parent::m_bin_number,0);
447 m_bin_Sv2w.resize(parent::m_bin_number,0);
448 m_cut_v = true;
449 m_min_v = aVmin;
450 m_max_v = aVmax;
451 return true;
452 }
453 public:
454 const vs_t& bins_sum_vw() const {return m_bin_Svw;}
455 const vs_t& bins_sum_v2w() const {return m_bin_Sv2w;}
456
457 TW get_Svw() const {
458 TW sw = 0;
459 for(TO ibin=0;ibin<parent::m_bin_number;ibin++) {
460 if(!histo::is_out(parent::m_axes,ibin)) {
461 sw += m_bin_Svw[ibin];
462 }
463 }
464 return sw;
465 }
466 TW get_Sv2w() const {
467 TW sw = 0;
468 for(TO ibin=0;ibin<parent::m_bin_number;ibin++) {
469 if(!histo::is_out(parent::m_axes,ibin)) {
470 sw += m_bin_Sv2w[ibin];
471 }
472 }
473 return sw;
474 }
475 protected:
476 bool m_cut_v;
477 TV m_min_v;
478 TV m_max_v;
479 vs_t m_bin_Svw;
480 vs_t m_bin_Sv2w;
481 };
482
483 }}
484
485 #endif
486