| Geant4 Cross Reference |
1 // Copyright (C) 2010, Guy Barrand. All rights 1 // Copyright (C) 2010, Guy Barrand. All rights reserved.
2 // See the file tools.license for terms. 2 // See the file tools.license for terms.
3 3
4 #ifndef tools_histo_h3 4 #ifndef tools_histo_h3
5 #define tools_histo_h3 5 #define tools_histo_h3
6 6
7 #include "b3" 7 #include "b3"
8 8
9 namespace tools { 9 namespace tools {
10 namespace histo { 10 namespace histo {
11 11
12 template <class TC,class TO,class TN,class TW, 12 template <class TC,class TO,class TN,class TW,class TH>
13 class h3 : public b3<TC,TO,TN,TW,TH> { 13 class h3 : public b3<TC,TO,TN,TW,TH> {
14 typedef b3<TC,TO,TN,TW,TH> parent; 14 typedef b3<TC,TO,TN,TW,TH> parent;
15 public: 15 public:
16 typedef histo_data<TC,TO,TN,TW> hd_t; 16 typedef histo_data<TC,TO,TN,TW> hd_t;
17 typedef typename parent::bn_t bn_t; 17 typedef typename parent::bn_t bn_t;
18 protected: 18 protected:
19 virtual TH get_bin_height(TO a_offset) const 19 virtual TH get_bin_height(TO a_offset) const { //TH should be the same as TW
20 return parent::m_bin_Sw[a_offset]; 20 return parent::m_bin_Sw[a_offset];
21 } 21 }
22 public: 22 public:
23 23
24 virtual TH bin_error(int aI,int aJ,int aK) c 24 virtual TH bin_error(int aI,int aJ,int aK) const {
25 TO offset; 25 TO offset;
26 if(!parent::_find_offset(aI,aJ,aK,offset)) 26 if(!parent::_find_offset(aI,aJ,aK,offset)) return 0;
27 return ::sqrt(parent::m_bin_Sw2[offset]); 27 return ::sqrt(parent::m_bin_Sw2[offset]);
28 } 28 }
29 29
30 public: 30 public:
31 bool multiply(TW a_factor){return parent::ba 31 bool multiply(TW a_factor){return parent::base_multiply(a_factor);}
32 bool scale(TW a_factor) {return multiply(a_f 32 bool scale(TW a_factor) {return multiply(a_factor);}
33 33
34 void copy_from_data(const hd_t& a_from) {par 34 void copy_from_data(const hd_t& a_from) {parent::base_from_data(a_from);}
35 hd_t get_histo_data() const {return *this;} 35 hd_t get_histo_data() const {return *this;} //deprecated. Keep it for g4tools.
36 36
37 bool reset() { 37 bool reset() {
38 parent::base_reset(); 38 parent::base_reset();
39 return true; 39 return true;
40 } 40 }
41 41
42 bool fill(TC aX,TC aY,TC aZ,TW aWeight = 1) 42 bool fill(TC aX,TC aY,TC aZ,TW aWeight = 1) {
43 if(parent::m_dimension!=3) return false; 43 if(parent::m_dimension!=3) return false;
44 44
45 bn_t ibin,jbin,kbin; 45 bn_t ibin,jbin,kbin;
46 if(!parent::m_axes[0].coord_to_absolute_in 46 if(!parent::m_axes[0].coord_to_absolute_index(aX,ibin)) return false;
47 if(!parent::m_axes[1].coord_to_absolute_in 47 if(!parent::m_axes[1].coord_to_absolute_index(aY,jbin)) return false;
48 if(!parent::m_axes[2].coord_to_absolute_in 48 if(!parent::m_axes[2].coord_to_absolute_index(aZ,kbin)) return false;
49 TO offset = ibin + jbin * parent::m_axes[1 49 TO offset = ibin + jbin * parent::m_axes[1].m_offset + kbin * parent::m_axes[2].m_offset;
50 50
51 parent::m_bin_entries[offset]++; 51 parent::m_bin_entries[offset]++;
52 parent::m_bin_Sw[offset] += aWeight; 52 parent::m_bin_Sw[offset] += aWeight;
53 parent::m_bin_Sw2[offset] += aWeight * aWe 53 parent::m_bin_Sw2[offset] += aWeight * aWeight;
54 54
55 TC xw = aX * aWeight; 55 TC xw = aX * aWeight;
56 TC x2w = aX * xw; 56 TC x2w = aX * xw;
57 parent::m_bin_Sxw[offset][0] += xw; 57 parent::m_bin_Sxw[offset][0] += xw;
58 parent::m_bin_Sx2w[offset][0] += x2w; 58 parent::m_bin_Sx2w[offset][0] += x2w;
59 59
60 TC yw = aY * aWeight; 60 TC yw = aY * aWeight;
61 TC y2w = aY * yw; 61 TC y2w = aY * yw;
62 parent::m_bin_Sxw[offset][1] += yw; 62 parent::m_bin_Sxw[offset][1] += yw;
63 parent::m_bin_Sx2w[offset][1] += y2w; 63 parent::m_bin_Sx2w[offset][1] += y2w;
64 64
65 TC zw = aZ * aWeight; 65 TC zw = aZ * aWeight;
66 TC z2w = aZ * zw; 66 TC z2w = aZ * zw;
67 parent::m_bin_Sxw[offset][2] += zw; 67 parent::m_bin_Sxw[offset][2] += zw;
68 parent::m_bin_Sx2w[offset][2] += z2w; 68 parent::m_bin_Sx2w[offset][2] += z2w;
69 69
70 bool inRange = true; 70 bool inRange = true;
71 if(ibin==0) inRange = false; 71 if(ibin==0) inRange = false;
72 else if(ibin==(parent::m_axes[0].m_number_ 72 else if(ibin==(parent::m_axes[0].m_number_of_bins+1)) inRange = false;
73 73
74 if(jbin==0) inRange = false; 74 if(jbin==0) inRange = false;
75 else if(jbin==(parent::m_axes[1].m_number_ 75 else if(jbin==(parent::m_axes[1].m_number_of_bins+1)) inRange = false;
76 76
77 if(kbin==0) inRange = false; 77 if(kbin==0) inRange = false;
78 else if(kbin==(parent::m_axes[2].m_number_ 78 else if(kbin==(parent::m_axes[2].m_number_of_bins+1)) inRange = false;
79 79
80 parent::m_all_entries++; 80 parent::m_all_entries++;
81 if(inRange) { 81 if(inRange) {
82 parent::m_in_range_plane_Sxyw[0] += aX * 82 parent::m_in_range_plane_Sxyw[0] += aX * aY * aWeight;
83 parent::m_in_range_plane_Sxyw[1] += aY * 83 parent::m_in_range_plane_Sxyw[1] += aY * aZ * aWeight;
84 parent::m_in_range_plane_Sxyw[2] += aZ * 84 parent::m_in_range_plane_Sxyw[2] += aZ * aX * aWeight;
85 85
86 // fast getters : 86 // fast getters :
87 parent::m_in_range_entries++; 87 parent::m_in_range_entries++;
88 parent::m_in_range_Sw += aWeight; 88 parent::m_in_range_Sw += aWeight;
89 parent::m_in_range_Sw2 += aWeight*aWeigh 89 parent::m_in_range_Sw2 += aWeight*aWeight;
90 90
91 parent::m_in_range_Sxw[0] += xw; 91 parent::m_in_range_Sxw[0] += xw;
92 parent::m_in_range_Sx2w[0] += x2w; 92 parent::m_in_range_Sx2w[0] += x2w;
93 93
94 parent::m_in_range_Sxw[1] += yw; 94 parent::m_in_range_Sxw[1] += yw;
95 parent::m_in_range_Sx2w[1] += y2w; 95 parent::m_in_range_Sx2w[1] += y2w;
96 96
97 parent::m_in_range_Sxw[2] += zw; 97 parent::m_in_range_Sxw[2] += zw;
98 parent::m_in_range_Sx2w[2] += z2w; 98 parent::m_in_range_Sx2w[2] += z2w;
99 } 99 }
100 100
101 return true; 101 return true;
102 } 102 }
103 103
104 bool set_bin_content(bn_t a_ibin,bn_t a_jbin 104 bool set_bin_content(bn_t a_ibin,bn_t a_jbin,bn_t a_kbin,
105 TN a_entries,TW a_Sw,TW 105 TN a_entries,TW a_Sw,TW a_Sw2,
106 TC a_Sxw,TC a_Sx2w,TC a 106 TC a_Sxw,TC a_Sx2w,TC a_Syw,TC a_Sy2w,TC a_Szw,TC a_Sz2w) {
107 if(parent::m_dimension!=3) return false; 107 if(parent::m_dimension!=3) return false;
108 if(a_ibin>(parent::m_axes[0].m_number_of_b 108 if(a_ibin>(parent::m_axes[0].m_number_of_bins+1)) return false;
109 if(a_jbin>(parent::m_axes[1].m_number_of_b 109 if(a_jbin>(parent::m_axes[1].m_number_of_bins+1)) return false;
110 if(a_kbin>(parent::m_axes[2].m_number_of_b 110 if(a_kbin>(parent::m_axes[2].m_number_of_bins+1)) return false;
111 111
112 bool inRange = true; 112 bool inRange = true;
113 if(a_ibin==0) inRange = false; 113 if(a_ibin==0) inRange = false;
114 else if(a_ibin==(parent::m_axes[0].m_numbe 114 else if(a_ibin==(parent::m_axes[0].m_number_of_bins+1)) inRange = false;
115 if(a_jbin==0) inRange = false; 115 if(a_jbin==0) inRange = false;
116 else if(a_jbin==(parent::m_axes[1].m_numbe 116 else if(a_jbin==(parent::m_axes[1].m_number_of_bins+1)) inRange = false;
117 if(a_kbin==0) inRange = false; 117 if(a_kbin==0) inRange = false;
118 else if(a_kbin==(parent::m_axes[2].m_numbe 118 else if(a_kbin==(parent::m_axes[2].m_number_of_bins+1)) inRange = false;
119 119
120 TO offset = a_ibin + a_jbin * parent::m_ax 120 TO offset = a_ibin + a_jbin * parent::m_axes[1].m_offset + a_kbin * parent::m_axes[2].m_offset;
121 121
122 parent::m_all_entries -= parent::m_bin_ent 122 parent::m_all_entries -= parent::m_bin_entries[offset];
123 if(inRange) { 123 if(inRange) {
124 parent::m_in_range_entries -= parent::m_ 124 parent::m_in_range_entries -= parent::m_bin_entries[offset];
125 parent::m_in_range_Sw -= parent::m_bin_S 125 parent::m_in_range_Sw -= parent::m_bin_Sw[offset];
126 parent::m_in_range_Sw2 -= parent::m_bin_ 126 parent::m_in_range_Sw2 -= parent::m_bin_Sw2[offset];
127 parent::m_in_range_Sxw[0] -= parent::m_b 127 parent::m_in_range_Sxw[0] -= parent::m_bin_Sxw[offset][0];
128 parent::m_in_range_Sx2w[0] -= parent::m_ 128 parent::m_in_range_Sx2w[0] -= parent::m_bin_Sx2w[offset][0];
129 parent::m_in_range_Sxw[1] -= parent::m_b 129 parent::m_in_range_Sxw[1] -= parent::m_bin_Sxw[offset][1];
130 parent::m_in_range_Sx2w[1] -= parent::m_ 130 parent::m_in_range_Sx2w[1] -= parent::m_bin_Sx2w[offset][1];
131 parent::m_in_range_Sxw[2] -= parent::m_b 131 parent::m_in_range_Sxw[2] -= parent::m_bin_Sxw[offset][2];
132 parent::m_in_range_Sx2w[2] -= parent::m_ 132 parent::m_in_range_Sx2w[2] -= parent::m_bin_Sx2w[offset][2];
133 } 133 }
134 134
135 parent::m_bin_entries[offset] = a_entries; 135 parent::m_bin_entries[offset] = a_entries;
136 parent::m_bin_Sw[offset] = a_Sw; 136 parent::m_bin_Sw[offset] = a_Sw;
137 parent::m_bin_Sw2[offset] = a_Sw2; 137 parent::m_bin_Sw2[offset] = a_Sw2;
138 138
139 parent::m_bin_Sxw[offset][0] = a_Sxw; 139 parent::m_bin_Sxw[offset][0] = a_Sxw;
140 parent::m_bin_Sx2w[offset][0] = a_Sx2w; 140 parent::m_bin_Sx2w[offset][0] = a_Sx2w;
141 parent::m_bin_Sxw[offset][1] = a_Syw; 141 parent::m_bin_Sxw[offset][1] = a_Syw;
142 parent::m_bin_Sx2w[offset][1] = a_Sy2w; 142 parent::m_bin_Sx2w[offset][1] = a_Sy2w;
143 parent::m_bin_Sxw[offset][2] = a_Szw; 143 parent::m_bin_Sxw[offset][2] = a_Szw;
144 parent::m_bin_Sx2w[offset][2] = a_Sz2w; 144 parent::m_bin_Sx2w[offset][2] = a_Sz2w;
145 145
146 parent::m_all_entries += a_entries; 146 parent::m_all_entries += a_entries;
147 if(inRange) { 147 if(inRange) {
148 //parent::m_in_range_plane_Sxyw[0,1,2] ? 148 //parent::m_in_range_plane_Sxyw[0,1,2] ??? ill-defined.
149 149
150 parent::m_in_range_entries += a_entries; 150 parent::m_in_range_entries += a_entries;
151 parent::m_in_range_Sw += a_Sw; 151 parent::m_in_range_Sw += a_Sw;
152 parent::m_in_range_Sw2 += a_Sw2; 152 parent::m_in_range_Sw2 += a_Sw2;
153 153
154 parent::m_in_range_Sxw[0] += a_Sxw; 154 parent::m_in_range_Sxw[0] += a_Sxw;
155 parent::m_in_range_Sx2w[0] += a_Sx2w; 155 parent::m_in_range_Sx2w[0] += a_Sx2w;
156 parent::m_in_range_Sxw[1] += a_Syw; 156 parent::m_in_range_Sxw[1] += a_Syw;
157 parent::m_in_range_Sx2w[1] += a_Sy2w; 157 parent::m_in_range_Sx2w[1] += a_Sy2w;
158 parent::m_in_range_Sxw[2] += a_Szw; 158 parent::m_in_range_Sxw[2] += a_Szw;
159 parent::m_in_range_Sx2w[2] += a_Sz2w; 159 parent::m_in_range_Sx2w[2] += a_Sz2w;
160 } 160 }
161 161
162 return true; 162 return true;
163 } 163 }
164 164
165 bool get_bin_content(bn_t a_ibin,bn_t a_jbin 165 bool get_bin_content(bn_t a_ibin,bn_t a_jbin,bn_t a_kbin,
166 TN& a_entries,TW& a_Sw, 166 TN& a_entries,TW& a_Sw,TW& a_Sw2,
167 TC& a_Sxw,TC& a_Sx2w, 167 TC& a_Sxw,TC& a_Sx2w,
168 TC& a_Syw,TC& a_Sy2w, 168 TC& a_Syw,TC& a_Sy2w,
169 TC& a_Szw,TC& a_Sz2w) { 169 TC& a_Szw,TC& a_Sz2w) {
170 if(parent::m_dimension!=3) { 170 if(parent::m_dimension!=3) {
171 a_entries = 0;a_Sw = 0;a_Sw2 = 0; 171 a_entries = 0;a_Sw = 0;a_Sw2 = 0;
172 a_Sxw = 0;a_Sx2w = 0; 172 a_Sxw = 0;a_Sx2w = 0;
173 a_Syw = 0;a_Sy2w = 0; 173 a_Syw = 0;a_Sy2w = 0;
174 a_Szw = 0;a_Sz2w = 0; 174 a_Szw = 0;a_Sz2w = 0;
175 return false; 175 return false;
176 } 176 }
177 if(a_ibin>(parent::m_axes[0].m_number_of_b 177 if(a_ibin>(parent::m_axes[0].m_number_of_bins+1)) {
178 a_entries = 0;a_Sw = 0;a_Sw2 = 0; 178 a_entries = 0;a_Sw = 0;a_Sw2 = 0;
179 a_Sxw = 0;a_Sx2w = 0; 179 a_Sxw = 0;a_Sx2w = 0;
180 a_Syw = 0;a_Sy2w = 0; 180 a_Syw = 0;a_Sy2w = 0;
181 a_Szw = 0;a_Sz2w = 0; 181 a_Szw = 0;a_Sz2w = 0;
182 return false; 182 return false;
183 } 183 }
184 if(a_jbin>(parent::m_axes[1].m_number_of_b 184 if(a_jbin>(parent::m_axes[1].m_number_of_bins+1)) {
185 a_entries = 0;a_Sw = 0;a_Sw2 = 0; 185 a_entries = 0;a_Sw = 0;a_Sw2 = 0;
186 a_Sxw = 0;a_Sx2w = 0; 186 a_Sxw = 0;a_Sx2w = 0;
187 a_Syw = 0;a_Sy2w = 0; 187 a_Syw = 0;a_Sy2w = 0;
188 a_Szw = 0;a_Sz2w = 0; 188 a_Szw = 0;a_Sz2w = 0;
189 return false; 189 return false;
190 } 190 }
191 if(a_kbin>(parent::m_axes[2].m_number_of_b 191 if(a_kbin>(parent::m_axes[2].m_number_of_bins+1)) {
192 a_entries = 0;a_Sw = 0;a_Sw2 = 0; 192 a_entries = 0;a_Sw = 0;a_Sw2 = 0;
193 a_Sxw = 0;a_Sx2w = 0; 193 a_Sxw = 0;a_Sx2w = 0;
194 a_Syw = 0;a_Sy2w = 0; 194 a_Syw = 0;a_Sy2w = 0;
195 a_Szw = 0;a_Sz2w = 0; 195 a_Szw = 0;a_Sz2w = 0;
196 return false; 196 return false;
197 } 197 }
198 198
199 TO offset = a_ibin + a_jbin * parent::m_ax 199 TO offset = a_ibin + a_jbin * parent::m_axes[1].m_offset + a_kbin * parent::m_axes[2].m_offset;
200 200
201 a_entries = parent::m_bin_entries[offset]; 201 a_entries = parent::m_bin_entries[offset];
202 a_Sw = parent::m_bin_Sw[offset]; 202 a_Sw = parent::m_bin_Sw[offset];
203 a_Sw2 = parent::m_bin_Sw2[offset]; 203 a_Sw2 = parent::m_bin_Sw2[offset];
204 204
205 a_Sxw = parent::m_bin_Sxw[offset][0]; 205 a_Sxw = parent::m_bin_Sxw[offset][0];
206 a_Sx2w = parent::m_bin_Sx2w[offset][0]; 206 a_Sx2w = parent::m_bin_Sx2w[offset][0];
207 a_Syw = parent::m_bin_Sxw[offset][1]; 207 a_Syw = parent::m_bin_Sxw[offset][1];
208 a_Sy2w = parent::m_bin_Sx2w[offset][1]; 208 a_Sy2w = parent::m_bin_Sx2w[offset][1];
209 a_Szw = parent::m_bin_Sxw[offset][2]; 209 a_Szw = parent::m_bin_Sxw[offset][2];
210 a_Sz2w = parent::m_bin_Sx2w[offset][2]; 210 a_Sz2w = parent::m_bin_Sx2w[offset][2];
211 211
212 return true; 212 return true;
213 } 213 }
214 214
215 bool add(const h3& a_histo){ 215 bool add(const h3& a_histo){
216 parent::base_add(a_histo); 216 parent::base_add(a_histo);
217 return true; 217 return true;
218 } 218 }
219 bool subtract(const h3& a_histo){ 219 bool subtract(const h3& a_histo){
220 parent::base_subtract(a_histo); 220 parent::base_subtract(a_histo);
221 return true; 221 return true;
222 } 222 }
223 223
224 bool multiply(const h3& a_histo) { 224 bool multiply(const h3& a_histo) {
225 return parent::base_multiply(a_histo); 225 return parent::base_multiply(a_histo);
226 } 226 }
227 227
228 bool divide(const h3& a_histo) { 228 bool divide(const h3& a_histo) {
229 return parent::base_divide(a_histo); 229 return parent::base_divide(a_histo);
230 } 230 }
231 231
232 bool equals_TH(const h3& a_from,const TW& a_ 232 bool equals_TH(const h3& a_from,const TW& a_prec,TW(*a_fabs)(TW)) const {
233 if(!parent::equals_TH(a_from,a_prec,a_fabs 233 if(!parent::equals_TH(a_from,a_prec,a_fabs,true)) return false;
234 return true; 234 return true;
235 } 235 }
236 236
237 void not_a_profile() const {} 237 void not_a_profile() const {}
238 238
239 public: //CERN-ROOT API 239 public: //CERN-ROOT API
240 bool Fill(TC aX,TC aY,TC aZ,TW aWeight = 1) 240 bool Fill(TC aX,TC aY,TC aZ,TW aWeight = 1) {return fill(aX,aY,aZ,aWeight);}
241 241
242 public: 242 public:
243 /* 243 /*
244 // Slices : 244 // Slices :
245 h2d* slice_xy(int aKbeg,int aKend) const; 245 h2d* slice_xy(int aKbeg,int aKend) const;
246 h2d* slice_yz(int aIbeg,int aIend) const; 246 h2d* slice_yz(int aIbeg,int aIend) const;
247 h2d* slice_xz(int aJbeg,int aJend) const; 247 h2d* slice_xz(int aJbeg,int aJend) const;
248 */ 248 */
249 public: 249 public:
250 h3(const std::string& a_title, 250 h3(const std::string& a_title,
251 bn_t aXnumber,TC aXmin,TC aXmax, 251 bn_t aXnumber,TC aXmin,TC aXmax,
252 bn_t aYnumber,TC aYmin,TC aYmax, 252 bn_t aYnumber,TC aYmin,TC aYmax,
253 bn_t aZnumber,TC aZmin,TC aZmax) 253 bn_t aZnumber,TC aZmin,TC aZmax)
254 :parent(a_title,aXnumber,aXmin,aXmax, 254 :parent(a_title,aXnumber,aXmin,aXmax,
255 aYnumber,aYmin,aYmax, 255 aYnumber,aYmin,aYmax,
256 aZnumber,aZmin,aZmax) 256 aZnumber,aZmin,aZmax)
257 {} 257 {}
258 258
259 h3(const std::string& a_title, 259 h3(const std::string& a_title,
260 const std::vector<TC>& a_edges_x, 260 const std::vector<TC>& a_edges_x,
261 const std::vector<TC>& a_edges_y, 261 const std::vector<TC>& a_edges_y,
262 const std::vector<TC>& a_edges_z) 262 const std::vector<TC>& a_edges_z)
263 :parent(a_title,a_edges_x,a_edges_y,a_edges_ 263 :parent(a_title,a_edges_x,a_edges_y,a_edges_z)
264 {} 264 {}
265 265
266 virtual ~h3(){} 266 virtual ~h3(){}
267 public: 267 public:
268 h3(const h3& a_from): parent(a_from){} 268 h3(const h3& a_from): parent(a_from){}
269 h3& operator=(const h3& a_from){ 269 h3& operator=(const h3& a_from){
270 parent::operator=(a_from); 270 parent::operator=(a_from);
271 return *this; 271 return *this;
272 } 272 }
273 }; 273 };
274 274
275 }} 275 }}
276 276
277 #endif 277 #endif
278 278
279 279
280 280
281 281