| 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_p2
5 #define tools_histo_p2
6
7 #include "b2"
8 #include "profile_data"
9
10 namespace tools {
11 namespace histo {
12
13 template <class TC,class TO,class TN,class TW,class TH,class TV>
14 class p2 : public b2<TC,TO,TN,TW,TH> {
15 typedef b2<TC,TO,TN,TW,TH> parent;
16 public:
17 typedef profile_data<TC,TO,TN,TW,TV> pd_t;
18 typedef typename parent::bn_t bn_t;
19 typedef std::vector<TV> vs_t;
20 protected:
21 virtual TH get_bin_height(TO a_offset) const {
22 return (parent::m_bin_Sw[a_offset] ? (m_bin_Svw[a_offset]/parent::m_bin_Sw[a_offset]):0);
23 }
24 public:
25 bool equals(const p2& a_from,const TW& a_prec,TW(*a_fabs)(TW)) const {
26 if(!parent::equals(a_from,a_prec,a_fabs)) return false;
27 if(m_cut_v!=a_from.m_cut_v) return false;
28 if(!numbers_are_equal(m_min_v,a_from.m_min_v,a_prec,a_fabs)) return false;
29 if(!numbers_are_equal(m_max_v,a_from.m_max_v,a_prec,a_fabs)) return false;
30 if(!vectors_are_equal(m_bin_Svw,a_from.m_bin_Svw,a_prec,a_fabs)) return false;
31 if(!vectors_are_equal(m_bin_Sv2w,a_from.m_bin_Sv2w,a_prec,a_fabs)) return false;
32 return true;
33 }
34 bool equals_TH(const p2& a_from,const TW& a_prec,TW(*a_fabs)(TW)) const {
35 if(!parent::equals_TH(a_from,a_prec,a_fabs,false)) return false;
36 if(m_cut_v!=a_from.m_cut_v) return false;
37 if(!numbers_are_equal(m_min_v,a_from.m_min_v,a_prec,a_fabs)) return false;
38 if(!numbers_are_equal(m_max_v,a_from.m_max_v,a_prec,a_fabs)) return false;
39 if(!vectors_are_equal(m_bin_Svw,a_from.m_bin_Svw,a_prec,a_fabs)) return false;
40 if(!vectors_are_equal(m_bin_Sv2w,a_from.m_bin_Sv2w,a_prec,a_fabs)) return false;
41 return true;
42 }
43
44 virtual TH bin_error(int aI,int aJ) const { //TH should be the same as TV
45 TO offset;
46 if(!parent::_find_offset(aI,aJ,offset)) return 0;
47
48 //FIXME Is it correct ?
49 // TProfile::GetBinError with kERRORMEAN mode does :
50 // Stat_t cont = fArray[bin]; //Svw (see TProfile::Fill)
51 // Stat_t sum = parent::m_bin_entries.fArray[bin]; //Sw
52 // Stat_t err2 = fSumw2.fArray[bin]; //Sv2w
53 // if (sum == 0) return 0;
54 // Stat_t eprim;
55 // Stat_t contsum = cont/sum;
56 // Stat_t eprim2 = TMath::Abs(err2/sum - contsum*contsum);
57 // eprim = TMath::Sqrt(eprim2);
58 // ... ???
59 // if (fErrorMode == kERRORMEAN) return eprim/TMath::Sqrt(sum);
60
61 TW sw = parent::m_bin_Sw[offset]; //ROOT sum
62 if(sw==0) return 0;
63 TV svw = m_bin_Svw[offset]; //ROOT cont
64 TV sv2w = m_bin_Sv2w[offset]; //ROOT err2
65 TV mean = (svw / sw); //ROOT contsum
66 TV rms = ::sqrt(::fabs((sv2w/sw) - mean * mean)); //ROOT eprim
67 // rms = get_bin_rms_value.
68 return rms/::sqrt(sw); //ROOT kERRORMEAN mode returned value
69 }
70
71 public:
72 bool multiply(TW a_factor){
73 if(!parent::base_multiply(a_factor)) return false;
74 for(bn_t ibin=0;ibin<parent::m_bin_number;ibin++) {
75 m_bin_Svw[ibin] *= a_factor;
76 }
77 return true;
78 }
79 bool scale(TW a_factor) {return multiply(a_factor);}
80
81 TV bin_Svw(int aI,int aJ) const {
82 TO offset;
83 if(!parent::_find_offset(aI,aJ,offset)) return 0;
84 return m_bin_Svw[offset];
85 }
86 TV bin_Sv2w(int aI,int aJ) const {
87 TO offset;
88 if(!parent::_find_offset(aI,aJ,offset)) return 0;
89 return m_bin_Sv2w[offset];
90 }
91
92 bool reset() {
93 parent::base_reset();
94 for(bn_t ibin=0;ibin<parent::m_bin_number;ibin++) {
95 m_bin_Svw[ibin] = 0;
96 m_bin_Sv2w[ibin] = 0;
97 }
98 return true;
99 }
100
101 void copy_from_data(const pd_t& a_from) {
102 parent::base_from_data(a_from);
103 m_bin_Svw = a_from.m_bin_Svw;
104 m_bin_Sv2w = a_from.m_bin_Sv2w;
105 m_cut_v = a_from.m_cut_v;
106 m_min_v = a_from.m_min_v;
107 m_max_v = a_from.m_max_v;
108 }
109 pd_t get_histo_data() const {
110 pd_t hd(parent::dac());
111 hd.m_is_profile = true;
112 hd.m_bin_Svw = m_bin_Svw;
113 hd.m_bin_Sv2w = m_bin_Sv2w;
114 hd.m_cut_v = m_cut_v;
115 hd.m_min_v = m_min_v;
116 hd.m_max_v = m_max_v;
117 return hd;
118 }
119
120 bool fill(TC aX,TC aY,TV aV,TW aWeight = 1) {
121 if(parent::m_dimension!=2) return false;
122
123 if(m_cut_v) {
124 if( (aV<m_min_v) || (aV>=m_max_v) ) {
125 return true;
126 }
127 }
128
129 bn_t ibin,jbin;
130 if(!parent::m_axes[0].coord_to_absolute_index(aX,ibin)) return false;
131 if(!parent::m_axes[1].coord_to_absolute_index(aY,jbin)) return false;
132 bn_t offset = ibin + jbin * parent::m_axes[1].m_offset;
133
134 parent::m_bin_entries[offset]++;
135 parent::m_bin_Sw[offset] += aWeight;
136 parent::m_bin_Sw2[offset] += aWeight * aWeight;
137
138 TC xw = aX * aWeight;
139 TC x2w = aX * xw;
140 parent::m_bin_Sxw[offset][0] += xw;
141 parent::m_bin_Sx2w[offset][0] += x2w;
142
143 TC yw = aY * aWeight;
144 TC y2w = aY * yw;
145 parent::m_bin_Sxw[offset][1] += yw;
146 parent::m_bin_Sx2w[offset][1] += y2w;
147
148 bool inRange = true;
149 if(ibin==0) inRange = false;
150 else if(ibin==(parent::m_axes[0].m_number_of_bins+1)) inRange = false;
151
152 if(jbin==0) inRange = false;
153 else if(jbin==(parent::m_axes[1].m_number_of_bins+1)) inRange = false;
154
155 parent::m_all_entries++;
156 if(inRange) {
157 parent::m_in_range_plane_Sxyw[0] += aX * aY * aWeight;
158
159 // fast getters :
160 parent::m_in_range_entries++;
161 parent::m_in_range_Sw += aWeight;
162 parent::m_in_range_Sw2 += aWeight*aWeight;
163
164 parent::m_in_range_Sxw[0] += xw;
165 parent::m_in_range_Sx2w[0] += x2w;
166
167 parent::m_in_range_Sxw[1] += yw;
168 parent::m_in_range_Sx2w[1] += y2w;
169 }
170
171 // Profile part :
172 TV vw = aV * aWeight;
173 m_bin_Svw[offset] += vw;
174 m_bin_Sv2w[offset] += aV * vw;
175
176 return true;
177 }
178
179 TV bin_rms_value(int aI,int aJ) const {
180 TO offset;
181 if(!parent::_find_offset(aI,aJ,offset)) return 0;
182
183 TW sw = parent::m_bin_Sw[offset];
184 if(sw==0) return 0;
185 TV svw = m_bin_Svw[offset];
186 TV sv2w = m_bin_Sv2w[offset];
187 TV mean = (svw / sw);
188 return ::sqrt(::fabs((sv2w / sw) - mean * mean));
189 }
190
191 bool add(const p2& a_histo){
192 parent::base_add(a_histo);
193 for(bn_t ibin=0;ibin<parent::m_bin_number;ibin++) {
194 m_bin_Svw[ibin] += a_histo.m_bin_Svw[ibin];
195 m_bin_Sv2w[ibin] += a_histo.m_bin_Sv2w[ibin];
196 }
197 return true;
198 }
199 bool subtract(const p2& a_histo){
200 parent::base_subtract(a_histo);
201 for(bn_t ibin=0;ibin<parent::m_bin_number;ibin++) {
202 m_bin_Svw[ibin] -= a_histo.m_bin_Svw[ibin];
203 m_bin_Sv2w[ibin] -= a_histo.m_bin_Sv2w[ibin];
204 }
205 return true;
206 }
207
208 bool cut_v() const {return m_cut_v;}
209 TV min_v() const {return m_min_v;}
210 TV max_v() const {return m_max_v;}
211 public:
212 p2(const std::string& a_title,
213 bn_t aXnumber,TC aXmin,TC aXmax,
214 bn_t aYnumber,TC aYmin,TC aYmax)
215 :parent(a_title,aXnumber,aXmin,aXmax,aYnumber,aYmin,aYmax)
216 ,m_cut_v(false)
217 ,m_min_v(0)
218 ,m_max_v(0)
219 {
220 m_bin_Svw.resize(parent::m_bin_number,0);
221 m_bin_Sv2w.resize(parent::m_bin_number,0);
222 }
223
224 p2(const std::string& a_title,
225 bn_t aXnumber,TC aXmin,TC aXmax,
226 bn_t aYnumber,TC aYmin,TC aYmax,
227 TV aVmin,TV aVmax)
228 :parent(a_title,aXnumber,aXmin,aXmax,aYnumber,aYmin,aYmax)
229 ,m_cut_v(true)
230 ,m_min_v(aVmin)
231 ,m_max_v(aVmax)
232 {
233 m_bin_Svw.resize(parent::m_bin_number,0);
234 m_bin_Sv2w.resize(parent::m_bin_number,0);
235 }
236
237 p2(const std::string& a_title,
238 const std::vector<TC>& a_edges_x,
239 const std::vector<TC>& a_edges_y)
240 :parent(a_title,a_edges_x,a_edges_y)
241 ,m_cut_v(false)
242 ,m_min_v(0)
243 ,m_max_v(0)
244 {
245 m_bin_Svw.resize(parent::m_bin_number,0);
246 m_bin_Sv2w.resize(parent::m_bin_number,0);
247 }
248
249 p2(const std::string& a_title,
250 const std::vector<TC>& a_edges_x,
251 const std::vector<TC>& a_edges_y,
252 TV aVmin,TV aVmax)
253 :parent(a_title,a_edges_x,a_edges_y)
254 ,m_cut_v(true)
255 ,m_min_v(aVmin)
256 ,m_max_v(aVmax)
257 {
258 m_bin_Svw.resize(parent::m_bin_number,0);
259 m_bin_Sv2w.resize(parent::m_bin_number,0);
260 }
261
262 virtual ~p2(){}
263 public:
264 p2(const p2& a_from)
265 :parent(a_from)
266 ,m_cut_v(a_from.m_cut_v)
267 ,m_min_v(a_from.m_min_v)
268 ,m_max_v(a_from.m_max_v)
269 ,m_bin_Svw(a_from.m_bin_Svw)
270 ,m_bin_Sv2w(a_from.m_bin_Sv2w)
271 {}
272 p2& operator=(const p2& a_from){
273 parent::operator=(a_from);
274 m_cut_v = a_from.m_cut_v;
275 m_min_v = a_from.m_min_v;
276 m_max_v = a_from.m_max_v;
277 m_bin_Svw = a_from.m_bin_Svw;
278 m_bin_Sv2w = a_from.m_bin_Sv2w;
279 return *this;
280 }
281 public:
282 bool configure(bn_t aXnumber,TC aXmin,TC aXmax,bn_t aYnumber,TC aYmin,TC aYmax){
283 if(!parent::configure(aXnumber,aXmin,aXmax,aYnumber,aYmin,aYmax)) return false;
284 m_bin_Svw.clear();
285 m_bin_Sv2w.clear();
286 m_bin_Svw.resize(parent::m_bin_number,0);
287 m_bin_Sv2w.resize(parent::m_bin_number,0);
288 m_cut_v = false;
289 m_min_v = 0;
290 m_max_v = 0;
291 return true;
292 }
293 bool configure(const std::vector<TC>& a_edges_x,const std::vector<TC>& a_edges_y) {
294 if(!parent::configure(a_edges_x,a_edges_y)) return false;
295 m_bin_Svw.clear();
296 m_bin_Sv2w.clear();
297 m_bin_Svw.resize(parent::m_bin_number,0);
298 m_bin_Sv2w.resize(parent::m_bin_number,0);
299 m_cut_v = false;
300 m_min_v = 0;
301 m_max_v = 0;
302 return true;
303 }
304 bool configure(bn_t aXnumber,TC aXmin,TC aXmax,bn_t aYnumber,TC aYmin,TC aYmax,TV aVmin,TV aVmax){
305 if(!parent::configure(aXnumber,aXmin,aXmax,aYnumber,aYmin,aYmax)) return false;
306 m_bin_Svw.clear();
307 m_bin_Sv2w.clear();
308 m_bin_Svw.resize(parent::m_bin_number,0);
309 m_bin_Sv2w.resize(parent::m_bin_number,0);
310 m_cut_v = true;
311 m_min_v = aVmin;
312 m_max_v = aVmax;
313 return true;
314 }
315 bool configure(const std::vector<TC>& a_edges_x,const std::vector<TC>& a_edges_y,TV aVmin,TV aVmax) {
316 if(!parent::configure(a_edges_x,a_edges_y)) return false;
317 m_bin_Svw.clear();
318 m_bin_Sv2w.clear();
319 m_bin_Svw.resize(parent::m_bin_number,0);
320 m_bin_Sv2w.resize(parent::m_bin_number,0);
321 m_cut_v = true;
322 m_min_v = aVmin;
323 m_max_v = aVmax;
324 return true;
325 }
326 public:
327 const vs_t& bins_sum_vw() const {return m_bin_Svw;}
328 const vs_t& bins_sum_v2w() const {return m_bin_Sv2w;}
329
330 TW get_Svw() const {
331 TW sw = 0;
332 for(TO ibin=0;ibin<parent::m_bin_number;ibin++) {
333 if(!histo::is_out(parent::m_axes,ibin)) {
334 sw += m_bin_Svw[ibin];
335 }
336 }
337 return sw;
338 }
339 TW get_Sv2w() const {
340 TW sw = 0;
341 for(TO ibin=0;ibin<parent::m_bin_number;ibin++) {
342 if(!histo::is_out(parent::m_axes,ibin)) {
343 sw += m_bin_Sv2w[ibin];
344 }
345 }
346 return sw;
347 }
348 protected:
349 bool m_cut_v;
350 TV m_min_v;
351 TV m_max_v;
352 vs_t m_bin_Svw;
353 vs_t m_bin_Sv2w;
354 };
355
356 }}
357
358 #endif
359