Geant4 Cross Reference |
1 // Copyright (C) 2010, Guy Barrand. All rights 2 // See the file tools.license for terms. 3 4 #ifndef tools_histo_b3 5 #define tools_histo_b3 6 7 #include "base_histo" 8 9 #include <ostream> 10 11 namespace tools { 12 namespace histo { 13 14 template <class TC,class TO,class TN,class TW, 15 class b3 : public base_histo<TC,TO,TN,TW,TH> { 16 typedef base_histo<TC,TO,TN,TW,TH> parent; 17 public: 18 typedef base_histo<TC,TO,TN,TW,TH> base_hist 19 typedef typename parent::axis_t axis_t; 20 typedef typename parent::bn_t bn_t; 21 protected: 22 enum {AxisX=0,AxisY=1,AxisZ=2}; 23 public: 24 virtual TH bin_error(int,int,int) const = 0; 25 public: 26 // Partition : 27 int coord_to_index_x(TC aCoord) const { 28 return axis_x().coord_to_index(aCoord); 29 } 30 int coord_to_index_y(TC aCoord) const { 31 return axis_y().coord_to_index(aCoord); 32 } 33 int coord_to_index_z(TC aCoord) const { 34 return axis_z().coord_to_index(aCoord); 35 } 36 37 TC mean_x() const { 38 if(parent::m_in_range_Sw==0) return 0; 39 return parent::m_in_range_Sxw[0]/parent::m 40 } 41 42 TC mean_y() const { 43 if(parent::m_in_range_Sw==0) return 0; 44 return parent::m_in_range_Sxw[1]/parent::m 45 } 46 47 TC mean_z() const { 48 if(parent::m_in_range_Sw==0) return 0; 49 return parent::m_in_range_Sxw[2]/parent::m 50 } 51 52 TC rms_x() const { 53 if(parent::m_in_range_Sw==0) return 0; 54 TC mean = parent::m_in_range_Sxw[0]/parent 55 return ::sqrt(::fabs((parent::m_in_range_S 56 } 57 58 TC rms_y() const { 59 if(parent::m_in_range_Sw==0) return 0; 60 TC mean = parent::m_in_range_Sxw[1]/parent 61 return ::sqrt(::fabs((parent::m_in_range_S 62 } 63 64 TC rms_z() const { 65 if(parent::m_in_range_Sw==0) return 0; 66 TC mean = parent::m_in_range_Sxw[2]/parent 67 return ::sqrt(::fabs((parent::m_in_range_S 68 } 69 70 // bins : 71 TN bin_entries(int aI,int aJ,int aK) const { 72 TO offset; 73 if(!_find_offset(aI,aJ,aK,offset)) return 74 return parent::m_bin_entries[offset]; 75 } 76 77 TH bin_height(int aI,int aJ,int aK) const { 78 TO offset; 79 if(!_find_offset(aI,aJ,aK,offset)) return 80 return this->get_bin_height(offset); 81 } 82 83 TC bin_center_x(int aI) const {return parent 84 TC bin_center_y(int aJ) const {return parent 85 TC bin_center_z(int aK) const {return parent 86 87 TC bin_mean_x(int aI,int aJ,int aK) const { 88 TO offset; 89 if(!_find_offset(aI,aJ,aK,offset)) return 90 TW sw = parent::m_bin_Sw[offset]; 91 if(sw==0) return 0; 92 return parent::m_bin_Sxw[offset][AxisX]/sw 93 } 94 95 TC bin_mean_y(int aI,int aJ,int aK) const { 96 TO offset; 97 if(!_find_offset(aI,aJ,aK,offset)) return 98 TW sw = parent::m_bin_Sw[offset]; 99 if(sw==0) return 0; 100 return parent::m_bin_Sxw[offset][AxisY]/sw 101 } 102 103 TC bin_mean_z(int aI,int aJ,int aK) const { 104 TO offset; 105 if(!_find_offset(aI,aJ,aK,offset)) return 106 TW sw = parent::m_bin_Sw[offset]; 107 if(sw==0) return 0; 108 return parent::m_bin_Sxw[offset][AxisZ]/sw 109 } 110 111 TC bin_rms_x(int aI,int aJ,int aK) const { 112 TO offset; 113 if(!_find_offset(aI,aJ,aK,offset)) return 114 TW sw = parent::m_bin_Sw[offset]; 115 if(sw==0) return 0; 116 TC sxw = parent::m_bin_Sxw[offset][AxisX]; 117 TC sx2w = parent::m_bin_Sx2w[offset][AxisX 118 TC mean = sxw/sw; 119 return ::sqrt(::fabs((sx2w / sw) - mean * 120 } 121 122 TC bin_rms_y(int aI,int aJ,int aK) const { 123 TO offset; 124 if(!_find_offset(aI,aJ,aK,offset)) return 125 TW sw = parent::m_bin_Sw[offset]; 126 if(sw==0) return 0; 127 TC sxw = parent::m_bin_Sxw[offset][AxisY]; 128 TC sx2w = parent::m_bin_Sx2w[offset][AxisY 129 TC mean = sxw/sw; 130 return ::sqrt(::fabs((sx2w / sw) - mean * 131 } 132 133 TC bin_rms_z(int aI,int aJ,int aK) const { 134 TO offset; 135 if(!_find_offset(aI,aJ,aK,offset)) return 136 TW sw = parent::m_bin_Sw[offset]; 137 if(sw==0) return 0; 138 TC sxw = parent::m_bin_Sxw[offset][AxisZ]; 139 TC sx2w = parent::m_bin_Sx2w[offset][AxisZ 140 TC mean = sxw/sw; 141 return ::sqrt(::fabs((sx2w / sw) - mean * 142 } 143 144 // Axes : 145 const axis_t& axis_x() const {return parent: 146 const axis_t& axis_y() const {return parent: 147 const axis_t& axis_z() const {return parent: 148 axis_t& axis_x() {return parent::m_axes[0];} 149 axis_t& axis_y() {return parent::m_axes[1];} 150 axis_t& axis_z() {return parent::m_axes[2];} 151 152 // Projection : 153 TN bin_entries_x(int aI) const { 154 if(!parent::m_dimension) return 0; 155 bn_t ibin; 156 if(!parent::m_axes[0].in_range_to_absolute 157 bn_t jbin,kbin,offset; 158 bn_t ybins = parent::m_axes[1].bins()+2; 159 bn_t zbins = parent::m_axes[2].bins()+2; 160 TO yoffset = parent::m_axes[1].m_offset; 161 TO zoffset = parent::m_axes[2].m_offset; 162 TO joffset = ibin; 163 TN _entries = 0; 164 for(jbin=0;jbin<ybins;jbin++) { 165 //joffset = ibin + jbin * parent::m_axes 166 offset = joffset; 167 for(kbin=0;kbin<zbins;kbin++) { 168 //offset = joffset + kbin * parent::m_ 169 _entries += parent::m_bin_entries[offs 170 offset += zoffset; 171 } 172 joffset += yoffset; 173 } 174 return _entries; 175 } 176 177 TN bin_entries_y(int aJ) const { 178 if(!parent::m_dimension) return 0; 179 bn_t jbin; 180 if(!parent::m_axes[1].in_range_to_absolute 181 bn_t ibin,kbin; 182 TO offset; 183 bn_t xbins = parent::m_axes[0].bins()+2; 184 bn_t zbins = parent::m_axes[2].bins()+2; 185 TO yoffset = parent::m_axes[1].m_offset; 186 TO zoffset = parent::m_axes[2].m_offset; 187 TO joffset = jbin * yoffset; 188 TN _entries = 0; 189 for(ibin=0;ibin<xbins;ibin++) { 190 //joffset = ibin + jbin * parent::m_axes 191 offset = joffset; 192 for(kbin=0;kbin<zbins;kbin++) { 193 //offset = joffset + kbin * parent::m_ 194 _entries += parent::m_bin_entries[offs 195 offset += zoffset; 196 } 197 joffset++; 198 } 199 return _entries; 200 } 201 202 TN bin_entries_z(int aK) const { 203 if(!parent::m_dimension) return 0; 204 bn_t kbin; 205 if(!parent::m_axes[2].in_range_to_absolute 206 bn_t ibin,jbin; 207 TO offset; 208 bn_t xbins = parent::m_axes[0].bins()+2; 209 bn_t ybins = parent::m_axes[1].bins()+2; 210 TO yoffset = parent::m_axes[1].m_offset; 211 TO zoffset = parent::m_axes[2].m_offset; 212 TO koffset = kbin * zoffset; 213 TN _entries = 0; 214 for(ibin=0;ibin<xbins;ibin++) { 215 //koffset = ibin + kbin * parent::m_axes 216 offset = koffset; 217 for(jbin=0;jbin<ybins;jbin++) { 218 //offset = koffset + jbin * parent::m_ 219 _entries += parent::m_bin_entries[offs 220 offset += yoffset; 221 } 222 koffset++; 223 } 224 return _entries; 225 } 226 227 TW bin_height_x(int aI) const { 228 //to slow : return get_ith_axis_bin_height 229 if(!parent::m_dimension) return 0; 230 bn_t ibin; 231 if(!parent::m_axes[0].in_range_to_absolute 232 bn_t ybins = parent::m_axes[1].bins()+2; 233 bn_t zbins = parent::m_axes[2].bins()+2; 234 TO yoffset = parent::m_axes[1].m_offset; 235 TO zoffset = parent::m_axes[2].m_offset; 236 TO joffset = ibin; 237 TW sw = 0; 238 for(bn_t jbin=0;jbin<ybins;jbin++) { 239 //joffset = ibin + jbin * parent::m_axes 240 TO offset = joffset; 241 for(bn_t kbin=0;kbin<zbins;kbin++) { 242 //offset = joffset + kbin * parent::m_ 243 sw += this->get_bin_height(offset); 244 offset += zoffset; 245 } 246 joffset += yoffset; 247 } 248 return sw; 249 } 250 251 TW bin_height_y(int aJ) const { 252 if(!parent::m_dimension) return 0; 253 bn_t jbin; 254 if(!parent::m_axes[1].in_range_to_absolute 255 bn_t xbins = parent::m_axes[0].bins()+2; 256 bn_t zbins = parent::m_axes[2].bins()+2; 257 TO yoffset = parent::m_axes[1].m_offset; 258 TO zoffset = parent::m_axes[2].m_offset; 259 TO joffset = jbin * yoffset; 260 TW sw = 0; 261 for(bn_t ibin=0;ibin<xbins;ibin++) { 262 //joffset = ibin + jbin * parent::m_axes 263 TO offset = joffset; 264 for(bn_t kbin=0;kbin<zbins;kbin++) { 265 //offset = joffset + kbin * parent::m_ 266 sw += this->get_bin_height(offset); 267 offset += zoffset; 268 } 269 joffset++; 270 } 271 return sw; 272 } 273 274 TW bin_height_z(int aK) const { 275 if(!parent::m_dimension) return 0; 276 bn_t kbin; 277 if(!parent::m_axes[2].in_range_to_absolute 278 bn_t xbins = parent::m_axes[0].bins()+2; 279 bn_t ybins = parent::m_axes[1].bins()+2; 280 TO yoffset = parent::m_axes[1].m_offset; 281 TO zoffset = parent::m_axes[2].m_offset; 282 TO koffset = kbin * zoffset; 283 TW sw = 0; 284 for(bn_t ibin=0;ibin<xbins;ibin++) { 285 //koffset = ibin + kbin * parent::m_axes 286 TO offset = koffset; 287 for(bn_t jbin=0;jbin<ybins;jbin++) { 288 //offset = koffset + jbin * parent::m_ 289 sw += this->get_bin_height(offset); 290 offset += yoffset; 291 } 292 koffset++; 293 } 294 return sw; 295 } 296 297 TC Sxyw() const {return parent::m_in_range_p 298 TC Syzw() const {return parent::m_in_range_p 299 TC Szxw() const {return parent::m_in_range_p 300 public: 301 //NOTE : print is a Python keyword. 302 void hprint(std::ostream& a_out) { 303 // A la HPRINT. 304 a_out << parent::dimension() << parent::ti 305 a_out 306 << " * ENTRIES = " << parent::all_entrie 307 308 } 309 public: 310 b3(const std::string& a_title, 311 bn_t aXnumber,TC aXmin,TC aXmax, 312 bn_t aYnumber,TC aYmin,TC aYmax, 313 bn_t aZnumber,TC aZmin,TC aZmax) 314 { 315 parent::m_title = a_title; 316 std::vector<bn_t> nbins; 317 nbins.push_back(aXnumber); 318 nbins.push_back(aYnumber); 319 nbins.push_back(aZnumber); 320 std::vector<TC> mins; 321 mins.push_back(aXmin); 322 mins.push_back(aYmin); 323 mins.push_back(aZmin); 324 std::vector<TC> maxs; 325 maxs.push_back(aXmax); 326 maxs.push_back(aYmax); 327 maxs.push_back(aZmax); 328 parent::configure(3,nbins,mins,maxs); 329 } 330 331 b3(const std::string& a_title, 332 const std::vector<TC>& a_edges_x, 333 const std::vector<TC>& a_edges_y, 334 const std::vector<TC>& a_edges_z) 335 { 336 parent::m_title = a_title; 337 std::vector< std::vector<TC> > edges(3); 338 edges[0] = a_edges_x; 339 edges[1] = a_edges_y; 340 edges[2] = a_edges_z; 341 parent::configure(3,edges); 342 } 343 344 virtual ~b3(){} 345 protected: 346 b3(const b3& a_from):parent(a_from) {} 347 b3& operator=(const b3& a_from){parent::oper 348 public: 349 bool configure(bn_t aXnumber,TC aXmin,TC aXm 350 bn_t aYnumber,TC aYmin,TC aYm 351 bn_t aZnumber,TC aZmin,TC aZm 352 std::vector<bn_t> nbins; 353 nbins.push_back(aXnumber); 354 nbins.push_back(aYnumber); 355 nbins.push_back(aZnumber); 356 std::vector<TC> mins; 357 mins.push_back(aXmin); 358 mins.push_back(aYmin); 359 mins.push_back(aZmin); 360 std::vector<TC> maxs; 361 maxs.push_back(aXmax); 362 maxs.push_back(aYmax); 363 maxs.push_back(aZmax); 364 return parent::configure(3,nbins,mins,maxs 365 } 366 367 bool configure(const std::vector<TC>& a_edge 368 const std::vector<TC>& a_edge 369 const std::vector<TC>& a_edge 370 std::vector< std::vector<TC> > edges(3); 371 edges[0] = a_edges_x; 372 edges[1] = a_edges_y; 373 edges[2] = a_edges_z; 374 return parent::configure(3,edges); 375 } 376 377 protected: 378 bool _find_offset(int aI,int aJ,int aK,TO& a 379 if(parent::m_dimension!=3) {a_offset=0;ret 380 bn_t ibin,jbin,kbin; 381 if(!parent::m_axes[0].in_range_to_absolute 382 if(!parent::m_axes[1].in_range_to_absolute 383 if(!parent::m_axes[2].in_range_to_absolute 384 a_offset = ibin + jbin * parent::m_axes[1] 385 return true; 386 } 387 }; 388 389 }} 390 391 #endif 392 393 394 395