| Geant4 Cross Reference |
1 // Copyright (C) 2010, Guy Barrand. All rights
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,
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_bi
23 }
24 public:
25 bool equals(const p2& a_from,const TW& a_pre
26 if(!parent::equals(a_from,a_prec,a_fabs))
27 if(m_cut_v!=a_from.m_cut_v) return false;
28 if(!numbers_are_equal(m_min_v,a_from.m_min
29 if(!numbers_are_equal(m_max_v,a_from.m_max
30 if(!vectors_are_equal(m_bin_Svw,a_from.m_b
31 if(!vectors_are_equal(m_bin_Sv2w,a_from.m_
32 return true;
33 }
34 bool equals_TH(const p2& a_from,const TW& a_
35 if(!parent::equals_TH(a_from,a_prec,a_fabs
36 if(m_cut_v!=a_from.m_cut_v) return false;
37 if(!numbers_are_equal(m_min_v,a_from.m_min
38 if(!numbers_are_equal(m_max_v,a_from.m_max
39 if(!vectors_are_equal(m_bin_Svw,a_from.m_b
40 if(!vectors_are_equal(m_bin_Sv2w,a_from.m_
41 return true;
42 }
43
44 virtual TH bin_error(int aI,int aJ) const {
45 TO offset;
46 if(!parent::_find_offset(aI,aJ,offset)) re
47
48 //FIXME Is it correct ?
49 // TProfile::GetBinError with kERRORMEAN m
50 // Stat_t cont = fArray[bin];
51 // Stat_t sum = parent::m_bin_entries.fA
52 // Stat_t err2 = fSumw2.fArray[bin];
53 // if (sum == 0) return 0;
54 // Stat_t eprim;
55 // Stat_t contsum = cont/sum;
56 // Stat_t eprim2 = TMath::Abs(err2/sum -
57 // eprim = TMath::Sqrt(eprim2);
58 // ... ???
59 // if (fErrorMode == kERRORMEAN) return e
60
61 TW sw = parent::m_bin_Sw[offset]; //R
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 contsu
66 TV rms = ::sqrt(::fabs((sv2w/sw) - mean *
67 // rms = get_bin_rms_value.
68 return rms/::sqrt(sw); //ROOT kERRORMEAN m
69 }
70
71 public:
72 bool multiply(TW a_factor){
73 if(!parent::base_multiply(a_factor)) retur
74 for(bn_t ibin=0;ibin<parent::m_bin_number;
75 m_bin_Svw[ibin] *= a_factor;
76 }
77 return true;
78 }
79 bool scale(TW a_factor) {return multiply(a_f
80
81 TV bin_Svw(int aI,int aJ) const {
82 TO offset;
83 if(!parent::_find_offset(aI,aJ,offset)) re
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)) re
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;
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_in
131 if(!parent::m_axes[1].coord_to_absolute_in
132 bn_t offset = ibin + jbin * parent::m_axes
133
134 parent::m_bin_entries[offset]++;
135 parent::m_bin_Sw[offset] += aWeight;
136 parent::m_bin_Sw2[offset] += aWeight * aWe
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_
151
152 if(jbin==0) inRange = false;
153 else if(jbin==(parent::m_axes[1].m_number_
154
155 parent::m_all_entries++;
156 if(inRange) {
157 parent::m_in_range_plane_Sxyw[0] += aX *
158
159 // fast getters :
160 parent::m_in_range_entries++;
161 parent::m_in_range_Sw += aWeight;
162 parent::m_in_range_Sw2 += aWeight*aWeigh
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)) re
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 *
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;
194 m_bin_Svw[ibin] += a_histo.m_bin_Svw[ibi
195 m_bin_Sv2w[ibin] += a_histo.m_bin_Sv2w[i
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;
202 m_bin_Svw[ibin] -= a_histo.m_bin_Svw[ibi
203 m_bin_Sv2w[ibin] -= a_histo.m_bin_Sv2w[i
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,aYnumbe
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,aYnumbe
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 aXm
283 if(!parent::configure(aXnumber,aXmin,aXmax
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_edge
294 if(!parent::configure(a_edges_x,a_edges_y)
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 aXm
305 if(!parent::configure(aXnumber,aXmin,aXmax
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_edge
316 if(!parent::configure(a_edges_x,a_edges_y)
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_bi
328 const vs_t& bins_sum_v2w() const {return m_b
329
330 TW get_Svw() const {
331 TW sw = 0;
332 for(TO ibin=0;ibin<parent::m_bin_number;ib
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;ib
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