| Geant4 Cross Reference |
1 // Copyright (C) 2010, Guy Barrand. All rights
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
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,
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_bi
31 }
32 public:
33 bool equals(const p1& a_from,const TW& a_pre
34 if(!parent::equals(a_from,a_prec,a_fabs))
35 if(m_cut_v!=a_from.m_cut_v) return false;
36 if(!numbers_are_equal(m_min_v,a_from.m_min
37 if(!numbers_are_equal(m_max_v,a_from.m_max
38 if(!vectors_are_equal(m_bin_Svw,a_from.m_b
39 if(!vectors_are_equal(m_bin_Sv2w,a_from.m_
40 return true;
41 }
42 bool equals_TH(const p1& a_from,const TW& a_
43 if(!parent::equals_TH(a_from,a_prec,a_fabs
44 if(m_cut_v!=a_from.m_cut_v) return false;
45 if(!numbers_are_equal(m_min_v,a_from.m_min
46 if(!numbers_are_equal(m_max_v,a_from.m_max
47 if(!vectors_are_equal(m_bin_Svw,a_from.m_b
48 if(!vectors_are_equal(m_bin_Sv2w,a_from.m_
49 return true;
50 }
51
52 virtual TH bin_error(int aI) const { //TH sh
53 TO offset;
54 if(!parent::_find_offset(aI,offset)) retur
55
56 //FIXME Is it correct ?
57 // TProfile::GetBinError with kERRORMEAN m
58 // Stat_t cont = fArray[bin];
59 // Stat_t sum = parent::m_bin_entries.fA
60 // Stat_t err2 = fSumw2.fArray[bin];
61 // if (sum == 0) return 0;
62 // Stat_t eprim;
63 // Stat_t contsum = cont/sum;
64 // Stat_t eprim2 = TMath::Abs(err2/sum -
65 // eprim = TMath::Sqrt(eprim2);
66 // ... ???
67 // if (fErrorMode == kERRORMEAN) return e
68
69 TW sw = parent::m_bin_Sw[offset]; //R
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 conts
74 TV _rms = ::sqrt(::fabs((sv2w/sw) - _mean
75 // rms = get_bin_rms_value.
76 return _rms/::sqrt(sw); //ROOT kERRORMEAN
77 }
78
79 public:
80 bool multiply(TW a_factor){
81 if(!parent::base_multiply(a_factor)) retur
82 for(bn_t ibin=0;ibin<parent::m_bin_number;
83 m_bin_Svw[ibin] *= a_factor;
84 }
85 return true;
86 }
87 bool scale(TW a_factor) {return multiply(a_f
88
89 TV bin_Svw(int aI) const {
90 TO offset;
91 if(!parent::_find_offset(aI,offset)) retur
92 return m_bin_Svw[offset];
93 }
94 TV bin_Sv2w(int aI) const {
95 TO offset;
96 if(!parent::_find_offset(aI,offset)) retur
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;
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_in
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 * aWe
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_
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*aWeigh
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_entrie
174 if(parent::m_dimension!=1) return false;
175 if(a_ibin>(parent::m_axes[0].m_number_of_b
176
177 bool inRange = true;
178 if(a_ibin==0) inRange = false;
179 else if(a_ibin==(parent::m_axes[0].m_numbe
180
181 TO offset = a_ibin;
182
183 parent::m_all_entries -= parent::m_bin_ent
184 if(inRange) {
185 parent::m_in_range_entries -= parent::m_
186 parent::m_in_range_Sw -= parent::m_bin_S
187 parent::m_in_range_Sw2 -= parent::m_bin_
188 parent::m_in_range_Sxw[0] -= parent::m_b
189 parent::m_in_range_Sx2w[0] -= parent::m_
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_entri
217 if(parent::m_dimension!=1) {
218 a_entries = 0;a_Sw = 0;a_Sw2 = 0;a_Sxw =
219 return false;
220 }
221 if(a_ibin>(parent::m_axes[0].m_number_of_b
222 a_entries = 0;a_Sw = 0;a_Sw2 = 0;a_Sxw =
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)) retur
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 *
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;
256 m_bin_Svw[ibin] += a_histo.m_bin_Svw[ibi
257 m_bin_Sv2w[ibin] += a_histo.m_bin_Sv2w[i
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;
264 m_bin_Svw[ibin] -= a_histo.m_bin_Svw[ibi
265 m_bin_Sv2w[ibin] -= a_histo.m_bin_Sv2w[i
266 }
267 return true;
268 }
269
270 bool gather_bins(unsigned int a_factor) { //
271 if(!a_factor) return false;
272
273 // actual bin number must be a multiple of
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,_ax
286 } else {
287 const std::vector<TC>& _edges = _axis.ed
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 ed
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;if
306 offac = offset+ifac;
307 new_h->m_bin_entries[new_offset] += pa
308 new_h->m_bin_Sw[new_offset] += parent:
309 new_h->m_bin_Sw2[new_offset] += parent
310 new_h->m_bin_Sxw[new_offset][0] += par
311 new_h->m_bin_Sx2w[new_offset][0] += pa
312
313 new_h->m_bin_Svw[new_offset] += m_bin_
314 new_h->m_bin_Sv2w[new_offset] += m_bin
315 }
316 }
317
318 //underflow :
319 new_offset = 0;
320 offac = 0;
321 new_h->m_bin_entries[new_offset] = parent:
322 new_h->m_bin_Sw[new_offset] = parent::m_bi
323 new_h->m_bin_Sw2[new_offset] = parent::m_b
324 new_h->m_bin_Sxw[new_offset][0] = parent::
325 new_h->m_bin_Sx2w[new_offset][0] = parent:
326 new_h->m_bin_Svw[new_offset] = m_bin_Svw[o
327 new_h->m_bin_Sv2w[new_offset] = m_bin_Sv2w
328
329 //overflow :
330 new_offset = new_n+1;
331 offac = n+1;
332 new_h->m_bin_entries[new_offset] = parent:
333 new_h->m_bin_Sw[new_offset] = parent::m_bi
334 new_h->m_bin_Sw2[new_offset] = parent::m_b
335 new_h->m_bin_Sxw[new_offset][0] = parent::
336 new_h->m_bin_Sx2w[new_offset][0] = parent:
337 new_h->m_bin_Svw[new_offset] = m_bin_Svw[o
338 new_h->m_bin_Sv2w[new_offset] = m_bin_Sv2w
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,
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,
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::vec
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::vec
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 aXm
410 if(!parent::configure(aXnumber,aXmin,aXmax
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_edge
421 if(!parent::configure(a_edges)) return fal
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 aXm
432 if(!parent::configure(aXnumber,aXmin,aXmax
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_edge
443 if(!parent::configure(a_edges)) return fal
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_bi
455 const vs_t& bins_sum_v2w() const {return m_b
456
457 TW get_Svw() const {
458 TW sw = 0;
459 for(TO ibin=0;ibin<parent::m_bin_number;ib
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;ib
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