inlib
1.2.0
|
00001 // Copyright (C) 2010, Guy Barrand. All rights reserved. 00002 // See the file inlib.license for terms. 00003 00004 #ifndef inlib_histo_b3 00005 #define inlib_histo_b3 00006 00007 #include "base_histo" 00008 00009 #include <ostream> 00010 00011 namespace inlib { 00012 namespace histo { 00013 00014 template <class TC,class TN,class TW,class TH> 00015 class b3 : public base_histo<TC,TN,TW,TH> { 00016 typedef base_histo<TC,TN,TW,TH> parent; 00017 typedef axis<TC> axis_t; 00018 public: 00019 typedef typename base_histo<TC,TN,TW,TH>::bn_t bn_t; 00020 protected: 00021 enum {AxisX=0,AxisY=1,AxisZ=2}; 00022 public: 00023 virtual TH bin_error(int,int,int) const = 0; //for print 00024 public: 00025 void update_fast_getters() { 00026 m_in_range_entries = 0; 00027 m_in_range_Sw = 0; 00028 m_in_range_Sxw = 0; 00029 m_in_range_Syw = 0; 00030 m_in_range_Szw = 0; 00031 m_in_range_Sx2w = 0; 00032 m_in_range_Sy2w = 0; 00033 m_in_range_Sz2w = 0; 00034 bn_t ibin,jbin,kbin,joffset,offset; 00035 bn_t xbins = parent::m_axes[0].bins(); 00036 bn_t ybins = parent::m_axes[1].bins(); 00037 bn_t zbins = parent::m_axes[2].bins(); 00038 bn_t yoffset = parent::m_axes[1].m_offset; 00039 bn_t zoffset = parent::m_axes[2].m_offset; 00040 for(ibin=1;ibin<=xbins;ibin++) { 00041 joffset = ibin + yoffset; 00042 for(jbin=1;jbin<=ybins;jbin++) { 00043 //joffset = ibin + jbin * parent::m_axes[1].m_offset; 00044 offset = joffset + zoffset; 00045 for(kbin=1;kbin<=zbins;kbin++) { 00046 //offset = joffset + kbin * parent::m_axes[2].m_offset; 00047 00048 m_in_range_entries += parent::m_bin_entries[offset]; 00049 m_in_range_Sw += parent::m_bin_Sw[offset]; 00050 m_in_range_Sxw += parent::m_bin_Sxw[offset][0]; 00051 m_in_range_Syw += parent::m_bin_Sxw[offset][1]; 00052 m_in_range_Szw += parent::m_bin_Sxw[offset][2]; 00053 m_in_range_Sx2w += parent::m_bin_Sx2w[offset][0]; 00054 m_in_range_Sy2w += parent::m_bin_Sx2w[offset][1]; 00055 m_in_range_Sz2w += parent::m_bin_Sx2w[offset][2]; 00056 00057 offset += zoffset; 00058 } 00059 joffset += yoffset; 00060 } 00061 } 00062 } 00063 00064 // Partition : 00065 int coord_to_index_x(TC aCoord) const { 00066 return axis_x().coord_to_index(aCoord); 00067 } 00068 int coord_to_index_y(TC aCoord) const { 00069 return axis_y().coord_to_index(aCoord); 00070 } 00071 int coord_to_index_z(TC aCoord) const { 00072 return axis_z().coord_to_index(aCoord); 00073 } 00074 00075 TC mean_x() const { 00076 if(m_in_range_Sw==0) return 0; 00077 return m_in_range_Sxw/m_in_range_Sw; 00078 } 00079 00080 TC mean_y() const { 00081 if(m_in_range_Sw==0) return 0; 00082 return m_in_range_Syw/m_in_range_Sw; 00083 } 00084 00085 TC mean_z() const { 00086 if(m_in_range_Sw==0) return 0; 00087 return m_in_range_Szw/m_in_range_Sw; 00088 } 00089 00090 TC rms_x() const { 00091 if(m_in_range_Sw==0) return 0; 00092 TC mean = m_in_range_Sxw/m_in_range_Sw; 00093 return ::sqrt(::fabs((m_in_range_Sx2w / m_in_range_Sw) - mean * mean)); 00094 } 00095 00096 TC rms_y() const { 00097 if(m_in_range_Sw==0) return 0; 00098 TC mean = m_in_range_Syw/m_in_range_Sw; 00099 return ::sqrt(::fabs((m_in_range_Sy2w / m_in_range_Sw) - mean * mean)); 00100 } 00101 00102 TC rms_z() const { 00103 if(m_in_range_Sw==0) return 0; 00104 TC mean = m_in_range_Szw/m_in_range_Sw; 00105 return ::sqrt(::fabs((m_in_range_Sz2w / m_in_range_Sw) - mean * mean)); 00106 } 00107 00108 // bins : 00109 TN bin_entries(int aI,int aJ,int aK) const { 00110 if(parent::m_bin_number==0) return 0; 00111 bn_t ibin,jbin,kbin; 00112 if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) return 0; 00113 if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) return 0; 00114 if(!parent::m_axes[2].in_range_to_absolute_index(aK,kbin)) return 0; 00115 bn_t offset = ibin + jbin * parent::m_axes[1].m_offset + kbin * parent::m_axes[2].m_offset; 00116 return parent::m_bin_entries[offset]; 00117 } 00118 00119 TH bin_height(int aI,int aJ,int aK) const { 00120 if(parent::m_bin_number==0) return 0; 00121 bn_t ibin,jbin,kbin; 00122 if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) return 0; 00123 if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) return 0; 00124 if(!parent::m_axes[2].in_range_to_absolute_index(aK,kbin)) return 0; 00125 bn_t offset = ibin + jbin * parent::m_axes[1].m_offset + kbin * parent::m_axes[2].m_offset; 00126 return this->get_bin_height(offset); 00127 } 00128 00129 TC bin_center_x(int aI) const {return parent::m_axes[0].bin_center(aI);} 00130 TC bin_center_y(int aJ) const {return parent::m_axes[1].bin_center(aJ);} 00131 TC bin_center_z(int aK) const {return parent::m_axes[2].bin_center(aK);} 00132 00133 TC bin_mean_x(int aI,int aJ,int aK) const { 00134 if(parent::m_bin_number==0) return 0; 00135 bn_t ibin,jbin,kbin; 00136 if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) return 0; 00137 if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) return 0; 00138 if(!parent::m_axes[2].in_range_to_absolute_index(aK,kbin)) return 0; 00139 bn_t offset = ibin + jbin * parent::m_axes[1].m_offset + kbin * parent::m_axes[2].m_offset; 00140 TW sw = parent::m_bin_Sw[offset]; 00141 if(sw==0) return 0; 00142 return parent::m_bin_Sxw[offset][AxisX]/sw; 00143 } 00144 00145 TC bin_mean_y(int aI,int aJ,int aK) const { 00146 if(parent::m_bin_number==0) return 0; 00147 bn_t ibin,jbin,kbin; 00148 if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) return 0; 00149 if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) return 0; 00150 if(!parent::m_axes[2].in_range_to_absolute_index(aK,kbin)) return 0; 00151 bn_t offset = ibin + jbin * parent::m_axes[1].m_offset + kbin * parent::m_axes[2].m_offset; 00152 TW sw = parent::m_bin_Sw[offset]; 00153 if(sw==0) return 0; 00154 return parent::m_bin_Sxw[offset][AxisY]/sw; 00155 } 00156 00157 TC bin_mean_z(int aI,int aJ,int aK) const { 00158 if(parent::m_bin_number==0) return 0; 00159 bn_t ibin,jbin,kbin; 00160 if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) return 0; 00161 if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) return 0; 00162 if(!parent::m_axes[2].in_range_to_absolute_index(aK,kbin)) return 0; 00163 bn_t offset = ibin + jbin * parent::m_axes[1].m_offset + kbin * parent::m_axes[2].m_offset; 00164 TW sw = parent::m_bin_Sw[offset]; 00165 if(sw==0) return 0; 00166 return parent::m_bin_Sxw[offset][AxisZ]/sw; 00167 } 00168 00169 TC bin_rms_x(int aI,int aJ,int aK) const { 00170 if(parent::m_bin_number==0) return 0; 00171 bn_t ibin,jbin,kbin; 00172 if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) return 0; 00173 if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) return 0; 00174 if(!parent::m_axes[2].in_range_to_absolute_index(aK,kbin)) return 0; 00175 bn_t offset = ibin + jbin * parent::m_axes[1].m_offset + kbin * parent::m_axes[2].m_offset; 00176 TW sw = parent::m_bin_Sw[offset]; 00177 if(sw==0) return 0; 00178 TC sxw = parent::m_bin_Sxw[offset][AxisX]; 00179 TC sx2w = parent::m_bin_Sx2w[offset][AxisX]; 00180 TC mean = sxw/sw; 00181 return ::sqrt(::fabs((sx2w / sw) - mean * mean)); 00182 } 00183 00184 TC bin_rms_y(int aI,int aJ,int aK) const { 00185 if(parent::m_bin_number==0) return 0; 00186 bn_t ibin,jbin,kbin; 00187 if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) return 0; 00188 if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) return 0; 00189 if(!parent::m_axes[2].in_range_to_absolute_index(aK,kbin)) return 0; 00190 bn_t offset = ibin + jbin * parent::m_axes[1].m_offset + kbin * parent::m_axes[2].m_offset; 00191 TW sw = parent::m_bin_Sw[offset]; 00192 if(sw==0) return 0; 00193 TC sxw = parent::m_bin_Sxw[offset][AxisY]; 00194 TC sx2w = parent::m_bin_Sx2w[offset][AxisY]; 00195 TC mean = sxw/sw; 00196 return ::sqrt(::fabs((sx2w / sw) - mean * mean)); 00197 } 00198 00199 TC bin_rms_z(int aI,int aJ,int aK) const { 00200 if(parent::m_bin_number==0) return 0; 00201 bn_t ibin,jbin,kbin; 00202 if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) return 0; 00203 if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) return 0; 00204 if(!parent::m_axes[2].in_range_to_absolute_index(aK,kbin)) return 0; 00205 bn_t offset = ibin + jbin * parent::m_axes[1].m_offset + kbin * parent::m_axes[2].m_offset; 00206 TW sw = parent::m_bin_Sw[offset]; 00207 if(sw==0) return 0; 00208 TC sxw = parent::m_bin_Sxw[offset][AxisZ]; 00209 TC sx2w = parent::m_bin_Sx2w[offset][AxisZ]; 00210 TC mean = sxw/sw; 00211 return ::sqrt(::fabs((sx2w / sw) - mean * mean)); 00212 } 00213 00214 // Axes : 00215 const axis_t& axis_x() const {return parent::m_axes[0];} 00216 const axis_t& axis_y() const {return parent::m_axes[1];} 00217 const axis_t& axis_z() const {return parent::m_axes[2];} 00218 axis_t& axis_x() {return parent::m_axes[0];} //touchy 00219 axis_t& axis_y() {return parent::m_axes[1];} //touchy 00220 axis_t& axis_z() {return parent::m_axes[2];} //touchy 00221 00222 // Projection : 00223 TN bin_entries_x(int aI) const { 00224 if(!parent::m_dimension) return 0; 00225 bn_t ibin; 00226 if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) return 0; 00227 bn_t jbin,kbin,offset; 00228 bn_t ybins = parent::m_axes[1].bins()+2; 00229 bn_t zbins = parent::m_axes[2].bins()+2; 00230 bn_t yoffset = parent::m_axes[1].m_offset; 00231 bn_t zoffset = parent::m_axes[2].m_offset; 00232 bn_t joffset = ibin; 00233 TN entries = 0; 00234 for(jbin=0;jbin<ybins;jbin++) { 00235 //joffset = ibin + jbin * parent::m_axes[1].m_offset; 00236 offset = joffset; 00237 for(kbin=0;kbin<zbins;kbin++) { 00238 //offset = joffset + kbin * parent::m_axes[2].m_offset; 00239 entries += parent::m_bin_entries[offset]; 00240 offset += zoffset; 00241 } 00242 joffset += yoffset; 00243 } 00244 return entries; 00245 } 00246 00247 TN bin_entries_y(int aJ) const { 00248 if(!parent::m_dimension) return 0; 00249 bn_t jbin; 00250 if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) return 0; 00251 bn_t ibin,kbin,offset; 00252 bn_t xbins = parent::m_axes[0].bins()+2; 00253 bn_t zbins = parent::m_axes[2].bins()+2; 00254 bn_t yoffset = parent::m_axes[1].m_offset; 00255 bn_t zoffset = parent::m_axes[2].m_offset; 00256 bn_t joffset = jbin * yoffset; 00257 TN entries = 0; 00258 for(ibin=0;ibin<xbins;ibin++) { 00259 //joffset = ibin + jbin * parent::m_axes[1].m_offset; 00260 offset = joffset; 00261 for(kbin=0;kbin<zbins;kbin++) { 00262 //offset = joffset + kbin * parent::m_axes[2].m_offset; 00263 entries += parent::m_bin_entries[offset]; 00264 offset += zoffset; 00265 } 00266 joffset++; 00267 } 00268 return entries; 00269 } 00270 00271 TN bin_entries_z(int aK) const { 00272 if(!parent::m_dimension) return 0; 00273 bn_t kbin; 00274 if(!parent::m_axes[2].in_range_to_absolute_index(aK,kbin)) return 0; 00275 bn_t ibin,jbin,offset; 00276 bn_t xbins = parent::m_axes[0].bins()+2; 00277 bn_t ybins = parent::m_axes[1].bins()+2; 00278 bn_t yoffset = parent::m_axes[1].m_offset; 00279 bn_t zoffset = parent::m_axes[2].m_offset; 00280 bn_t koffset = kbin * zoffset; 00281 TN entries = 0; 00282 for(ibin=0;ibin<xbins;ibin++) { 00283 //koffset = ibin + kbin * parent::m_axes[2].m_offset; 00284 offset = koffset; 00285 for(jbin=0;jbin<ybins;jbin++) { 00286 //offset = koffset + jbin * parent::m_axes[1].m_offset; 00287 entries += parent::m_bin_entries[offset]; 00288 offset += yoffset; 00289 } 00290 koffset++; 00291 } 00292 return entries; 00293 } 00294 00295 TW bin_height_x(int aI) const { 00296 //to slow : return get_ith_axis_bin_height(0,aI); 00297 if(!parent::m_dimension) return 0; 00298 bn_t ibin; 00299 if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) return 0; 00300 bn_t ybins = parent::m_axes[1].bins()+2; 00301 bn_t zbins = parent::m_axes[2].bins()+2; 00302 bn_t yoffset = parent::m_axes[1].m_offset; 00303 bn_t zoffset = parent::m_axes[2].m_offset; 00304 bn_t joffset = ibin; 00305 TW sw = 0; 00306 for(bn_t jbin=0;jbin<ybins;jbin++) { 00307 //joffset = ibin + jbin * parent::m_axes[1].m_offset; 00308 bn_t offset = joffset; 00309 for(bn_t kbin=0;kbin<zbins;kbin++) { 00310 //offset = joffset + kbin * parent::m_axes[2].m_offset; 00311 sw += this->get_bin_height(offset); 00312 offset += zoffset; 00313 } 00314 joffset += yoffset; 00315 } 00316 return sw; 00317 } 00318 00319 TW bin_height_y(int aJ) const { 00320 if(!parent::m_dimension) return 0; 00321 bn_t jbin; 00322 if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) return 0; 00323 bn_t xbins = parent::m_axes[0].bins()+2; 00324 bn_t zbins = parent::m_axes[2].bins()+2; 00325 bn_t yoffset = parent::m_axes[1].m_offset; 00326 bn_t zoffset = parent::m_axes[2].m_offset; 00327 bn_t joffset = jbin * yoffset; 00328 TW sw = 0; 00329 for(bn_t ibin=0;ibin<xbins;ibin++) { 00330 //joffset = ibin + jbin * parent::m_axes[1].m_offset; 00331 bn_t offset = joffset; 00332 for(bn_t kbin=0;kbin<zbins;kbin++) { 00333 //offset = joffset + kbin * parent::m_axes[2].m_offset; 00334 sw += this->get_bin_height(offset); 00335 offset += zoffset; 00336 } 00337 joffset++; 00338 } 00339 return sw; 00340 } 00341 00342 TW bin_height_z(int aK) const { 00343 if(!parent::m_dimension) return 0; 00344 bn_t kbin; 00345 if(!parent::m_axes[2].in_range_to_absolute_index(aK,kbin)) return 0; 00346 bn_t xbins = parent::m_axes[0].bins()+2; 00347 bn_t ybins = parent::m_axes[1].bins()+2; 00348 bn_t yoffset = parent::m_axes[1].m_offset; 00349 bn_t zoffset = parent::m_axes[2].m_offset; 00350 bn_t koffset = kbin * zoffset; 00351 TW sw = 0; 00352 for(bn_t ibin=0;ibin<xbins;ibin++) { 00353 //koffset = ibin + kbin * parent::m_axes[2].m_offset; 00354 bn_t offset = koffset; 00355 for(bn_t jbin=0;jbin<ybins;jbin++) { 00356 //offset = koffset + jbin * parent::m_axes[1].m_offset; 00357 sw += this->get_bin_height(offset); 00358 offset += yoffset; 00359 } 00360 koffset++; 00361 } 00362 return sw; 00363 } 00364 public: 00365 b3(const std::string& a_title, 00366 bn_t aXnumber,TC aXmin,TC aXmax, 00367 bn_t aYnumber,TC aYmin,TC aYmax, 00368 bn_t aZnumber,TC aZmin,TC aZmax) 00369 :m_in_range_entries(0) 00370 ,m_in_range_Sw(0) 00371 ,m_in_range_Sxw(0) 00372 ,m_in_range_Syw(0) 00373 ,m_in_range_Szw(0) 00374 ,m_in_range_Sx2w(0) 00375 ,m_in_range_Sy2w(0) 00376 ,m_in_range_Sz2w(0) 00377 { 00378 parent::m_title = a_title; 00379 std::vector<bn_t> ns; 00380 ns.push_back(aXnumber); 00381 ns.push_back(aYnumber); 00382 ns.push_back(aZnumber); 00383 std::vector<TC> mins; 00384 mins.push_back(aXmin); 00385 mins.push_back(aYmin); 00386 mins.push_back(aZmin); 00387 std::vector<TC> maxs; 00388 maxs.push_back(aXmax); 00389 maxs.push_back(aYmax); 00390 maxs.push_back(aZmax); 00391 parent::configure(3,ns,mins,maxs); 00392 } 00393 00394 b3(const std::string& a_title, 00395 const std::vector<TC>& aEdgesX, 00396 const std::vector<TC>& aEdgesY, 00397 const std::vector<TC>& aEdgesZ) 00398 :m_in_range_entries(0) 00399 ,m_in_range_Sw(0) 00400 ,m_in_range_Sxw(0) 00401 ,m_in_range_Syw(0) 00402 ,m_in_range_Szw(0) 00403 ,m_in_range_Sx2w(0) 00404 ,m_in_range_Sy2w(0) 00405 ,m_in_range_Sz2w(0) 00406 { 00407 parent::m_title = a_title; 00408 std::vector< std::vector<TC> > edges(3); 00409 edges[0] = aEdgesX; 00410 edges[1] = aEdgesY; 00411 edges[2] = aEdgesZ; 00412 parent::configure(3,edges); 00413 } 00414 00415 virtual ~b3(){} 00416 protected: 00417 b3(const b3& a_from) 00418 : parent(a_from) 00419 ,m_in_range_entries(a_from.m_in_range_entries) 00420 ,m_in_range_Sw(a_from.m_in_range_Sw) 00421 ,m_in_range_Sxw(a_from.m_in_range_Sxw) 00422 ,m_in_range_Syw(a_from.m_in_range_Syw) 00423 ,m_in_range_Szw(a_from.m_in_range_Szw) 00424 ,m_in_range_Sx2w(a_from.m_in_range_Sx2w) 00425 ,m_in_range_Sy2w(a_from.m_in_range_Sy2w) 00426 ,m_in_range_Sz2w(a_from.m_in_range_Sz2w) 00427 { 00428 update_fast_getters(); 00429 } 00430 b3& operator=(const b3& a_from){ 00431 parent::operator=(a_from); 00432 m_in_range_entries = a_from.m_in_range_entries; 00433 m_in_range_Sw = a_from.m_in_range_Sw; 00434 m_in_range_Sxw = a_from.m_in_range_Sxw; 00435 m_in_range_Syw = a_from.m_in_range_Syw; 00436 m_in_range_Szw = a_from.m_in_range_Szw; 00437 m_in_range_Sx2w = a_from.m_in_range_Sx2w; 00438 m_in_range_Sy2w = a_from.m_in_range_Sy2w; 00439 m_in_range_Sz2w = a_from.m_in_range_Sz2w; 00440 update_fast_getters(); 00441 return *this; 00442 } 00443 protected: 00444 TN m_in_range_entries; 00445 TW m_in_range_Sw; 00446 TC m_in_range_Sxw; 00447 TC m_in_range_Syw; 00448 TC m_in_range_Szw; 00449 TC m_in_range_Sx2w; 00450 TC m_in_range_Sy2w; 00451 TC m_in_range_Sz2w; 00452 }; 00453 00454 }} 00455 00456 #endif 00457 00458 00459 00460