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_b2 00005 #define inlib_histo_b2 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 b2 : public base_histo<TC,TN,TW,TH> { 00016 typedef base_histo<TC,TN,TW,TH> parent; 00017 typedef axis<TC> axis_t; 00018 protected: 00019 enum {AxisX=0,AxisY=1}; 00020 public: 00021 typedef typename base_histo<TC,TN,TW,TH>::bn_t bn_t; 00022 public: 00023 virtual TH bin_error(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_Sx2w = 0; 00031 m_in_range_Sy2w = 0; 00032 bn_t ibin,jbin,offset; 00033 bn_t xbins = parent::m_axes[0].bins(); 00034 bn_t ybins = parent::m_axes[1].bins(); 00035 bn_t yoffset = parent::m_axes[1].m_offset; 00036 for(ibin=1;ibin<=xbins;ibin++) { 00037 offset = ibin + yoffset; 00038 for(jbin=1;jbin<=ybins;jbin++) { 00039 //offset = ibin + jbin * m_axes[1].m_offset; 00040 00041 m_in_range_entries += parent::m_bin_entries[offset]; 00042 m_in_range_Sw += parent::m_bin_Sw[offset]; 00043 m_in_range_Sxw += parent::m_bin_Sxw[offset][0]; 00044 m_in_range_Syw += parent::m_bin_Sxw[offset][1]; 00045 m_in_range_Sx2w += parent::m_bin_Sx2w[offset][0]; 00046 m_in_range_Sy2w += parent::m_bin_Sx2w[offset][1]; 00047 00048 offset += yoffset; 00049 } 00050 } 00051 } 00052 00053 // Partition : 00054 TC mean_x() const { 00055 if(m_in_range_Sw==0) return 0; 00056 return m_in_range_Sxw/m_in_range_Sw; 00057 } 00058 00059 TC mean_y() const { 00060 if(m_in_range_Sw==0) return 0; 00061 return m_in_range_Syw/m_in_range_Sw; 00062 } 00063 00064 TC rms_x() const { 00065 if(m_in_range_Sw==0) return 0; 00066 TC mean = m_in_range_Sxw/m_in_range_Sw; 00067 return ::sqrt(::fabs((m_in_range_Sx2w / m_in_range_Sw) - mean * mean)); 00068 } 00069 00070 TC rms_y() const { 00071 if(m_in_range_Sw==0) return 0; 00072 TC mean = m_in_range_Syw/m_in_range_Sw; 00073 return ::sqrt(::fabs((m_in_range_Sy2w / m_in_range_Sw) - mean * mean)); 00074 } 00075 00076 int coord_to_index_x(TC aCoord) const { 00077 return axis_x().coord_to_index(aCoord); 00078 } 00079 int coord_to_index_y(TC aCoord) const { 00080 return axis_y().coord_to_index(aCoord); 00081 } 00082 00083 // bins : 00084 TN bin_entries(int aI,int aJ) const { 00085 if(parent::m_bin_number==0) return 0; 00086 bn_t ibin,jbin; 00087 if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) return 0; 00088 if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) return 0; 00089 bn_t offset = ibin + jbin * parent::m_axes[1].m_offset; 00090 return parent::m_bin_entries[offset]; 00091 } 00092 00093 TW bin_Sw(int aI,int aJ) const { 00094 if(parent::m_bin_number==0) return 0; 00095 bn_t ibin,jbin; 00096 if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) return 0; 00097 if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) return 0; 00098 bn_t offset = ibin + jbin * parent::m_axes[1].m_offset; 00099 return parent::m_bin_Sw(offset); 00100 } 00101 00102 TW bin_Sw2(int aI,int aJ) const { 00103 if(parent::m_bin_number==0) return 0; 00104 bn_t ibin,jbin; 00105 if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) return 0; 00106 if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) return 0; 00107 bn_t offset = ibin + jbin * parent::m_axes[1].m_offset; 00108 return parent::m_bin_Sw2(offset); 00109 } 00110 TC bin_Sxw(int aI,int aJ) const { 00111 if(parent::m_bin_number==0) return 0; 00112 bn_t ibin,jbin; 00113 if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) return 0; 00114 if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) return 0; 00115 bn_t offset = ibin + jbin * parent::m_axes[1].m_offset; 00116 return parent::m_bin_Sxw[offset][AxisX]; 00117 } 00118 TC bin_Sx2w(int aI,int aJ) const { 00119 if(parent::m_bin_number==0) return 0; 00120 bn_t ibin,jbin; 00121 if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) return 0; 00122 if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) return 0; 00123 bn_t offset = ibin + jbin * parent::m_axes[1].m_offset; 00124 return parent::m_bin_Sx2w[offset][AxisX]; 00125 } 00126 TC bin_Syw(int aI,int aJ) const { 00127 if(parent::m_bin_number==0) return 0; 00128 bn_t ibin,jbin; 00129 if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) return 0; 00130 if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) return 0; 00131 bn_t offset = ibin + jbin * parent::m_axes[1].m_offset; 00132 return parent::m_bin_Sxw[offset][AxisY]; 00133 } 00134 TC bin_Sy2w(int aI,int aJ) const { 00135 if(parent::m_bin_number==0) return 0; 00136 bn_t ibin,jbin; 00137 if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) return 0; 00138 if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) return 0; 00139 bn_t offset = ibin + jbin * parent::m_axes[1].m_offset; 00140 return parent::m_bin_Sx2w[offset][AxisY]; 00141 } 00142 00143 TH bin_height(int aI,int aJ) const { 00144 if(parent::m_bin_number==0) return 0; 00145 bn_t ibin,jbin; 00146 if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) return 0; 00147 if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) return 0; 00148 bn_t offset = ibin + jbin * parent::m_axes[1].m_offset; 00149 return this->get_bin_height(offset); 00150 } 00151 00152 TC bin_center_x(int aI) const { 00153 return parent::m_axes[0].bin_center(aI); 00154 } 00155 TC bin_center_y(int aJ) const { 00156 return parent::m_axes[1].bin_center(aJ); 00157 } 00158 00159 TC bin_mean_x(int aI,int aJ) const { 00160 if(parent::m_bin_number==0) return 0; 00161 bn_t ibin,jbin; 00162 if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) return 0; 00163 if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) return 0; 00164 bn_t offset = ibin + jbin * parent::m_axes[1].m_offset; 00165 TW sw = parent::m_bin_Sw[offset]; 00166 if(sw==0) return 0; 00167 return parent::m_bin_Sxw[offset][AxisX]/sw; 00168 } 00169 00170 TC bin_mean_y(int aI,int aJ) const { 00171 if(parent::m_bin_number==0) return 0; 00172 bn_t ibin,jbin; 00173 if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) return 0; 00174 if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) return 0; 00175 bn_t offset = ibin + jbin * parent::m_axes[1].m_offset; 00176 TW sw = parent::m_bin_Sw[offset]; 00177 if(sw==0) return 0; 00178 return parent::m_bin_Sxw[offset][AxisY]/sw; 00179 } 00180 00181 TC bin_rms_x(int aI,int aJ) const { 00182 if(parent::m_bin_number==0) return 0; 00183 bn_t ibin,jbin; 00184 if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) return 0; 00185 if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) return 0; 00186 bn_t offset = ibin + jbin * parent::m_axes[1].m_offset; 00187 TW sw = parent::m_bin_Sw[offset]; 00188 if(sw==0) return 0; 00189 TC sxw = parent::m_bin_Sxw[offset][AxisX]; 00190 TC sx2w = parent::m_bin_Sx2w[offset][AxisX]; 00191 TC mean = sxw/sw; 00192 return ::sqrt(::fabs((sx2w / sw) - mean * mean)); 00193 } 00194 00195 TC bin_rms_y(int aI,int aJ) const { 00196 if(parent::m_bin_number==0) return 0; 00197 bn_t ibin,jbin; 00198 if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) return 0; 00199 if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) return 0; 00200 bn_t offset = ibin + jbin * parent::m_axes[1].m_offset; 00201 TW sw = parent::m_bin_Sw[offset]; 00202 if(sw==0) return 0; 00203 TC sxw = parent::m_bin_Sxw[offset][AxisY]; 00204 TC sx2w = parent::m_bin_Sx2w[offset][AxisY]; 00205 TC mean = sxw/sw; 00206 return ::sqrt(::fabs((sx2w / sw) - mean * mean)); 00207 } 00208 00209 // Axes : 00210 const axis<TC>& axis_x() const {return parent::m_axes[0];} 00211 const axis<TC>& axis_y() const {return parent::m_axes[1];} 00212 axis<TC>& axis_x() {return parent::m_axes[0];} //touchy 00213 axis<TC>& axis_y() {return parent::m_axes[1];} //touchy 00214 // Projection : 00215 00216 TN bin_entries_x(int aI) const { 00217 if(!parent::m_dimension) return 0; 00218 bn_t ibin; 00219 if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) return 0; 00220 bn_t ybins = parent::m_axes[1].bins()+2; 00221 bn_t offset; 00222 TN entries = 0; 00223 for(bn_t jbin=0;jbin<ybins;jbin++) { 00224 offset = ibin + jbin * parent::m_axes[1].m_offset; 00225 entries += parent::m_bin_entries[offset]; 00226 } 00227 return entries; 00228 } 00229 00230 TW bin_height_x(int aI) const { 00231 if(!parent::m_dimension) return 0; 00232 bn_t ibin; 00233 if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) return 0; 00234 bn_t ybins = parent::m_axes[1].bins()+2; 00235 bn_t offset; 00236 TW sw = 0; 00237 for(bn_t jbin=0;jbin<ybins;jbin++) { 00238 offset = ibin + jbin * parent::m_axes[1].m_offset; 00239 sw += this->get_bin_height(offset); 00240 } 00241 return sw; 00242 } 00243 00244 TN bin_entries_y(int aJ) const { 00245 if(!parent::m_dimension) return 0; 00246 bn_t jbin; 00247 if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) return 0; 00248 bn_t xbins = parent::m_axes[0].bins()+2; 00249 bn_t offset; 00250 TN entries = 0; 00251 for(bn_t ibin=0;ibin<xbins;ibin++) { 00252 offset = ibin + jbin * parent::m_axes[1].m_offset; 00253 entries += parent::m_bin_entries[offset]; 00254 } 00255 return entries; 00256 } 00257 00258 TW bin_height_y(int aJ) const { 00259 if(!parent::m_dimension) return 0; 00260 bn_t jbin; 00261 if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) return 0; 00262 bn_t xbins = parent::m_axes[0].bins()+2; 00263 bn_t offset; 00264 TW sw = 0; 00265 for(bn_t ibin=0;ibin<xbins;ibin++) { 00266 offset = ibin + jbin * parent::m_axes[1].m_offset; 00267 sw += this->get_bin_height(offset); 00268 } 00269 return sw; 00270 } 00271 00272 public: 00273 //NOTE : print is a Python keyword. 00274 void hprint(std::ostream& a_out) { 00275 // A la HPRINT. 00276 a_out << parent::dimension() << parent::title() << std::endl; 00277 a_out 00278 << " * ENTRIES = " << parent::all_entries() << std::endl; 00279 00280 // 6 | 7 | 8 00281 // ----------- 00282 // 3 | 4 | 5 00283 // ----------- 00284 // 0 | 1 | 2 00285 00286 TW height_0 = bin_height(axis_t::UNDERFLOW_BIN, 00287 axis_t::UNDERFLOW_BIN); 00288 TW height_2 = bin_height(axis_t::OVERFLOW_BIN, 00289 axis_t::UNDERFLOW_BIN); 00290 TW height_6 = bin_height(axis_t::UNDERFLOW_BIN, 00291 axis_t::OVERFLOW_BIN); 00292 TW height_8 = bin_height(axis_t::OVERFLOW_BIN, 00293 axis_t::OVERFLOW_BIN); 00294 00295 bn_t i,j; 00296 00297 TW height_1 = 0; 00298 TW height_7 = 0; 00299 for(i=0;i<axis_x().bins();i++){ 00300 height_1 += bin_height(i,axis_t::UNDERFLOW_BIN); 00301 height_7 += bin_height(i,axis_t::OVERFLOW_BIN); 00302 } 00303 00304 TW height_3 = 0; 00305 TW height_5 = 0; 00306 for(j=0;j<axis_y().bins();j++){ 00307 height_3 += bin_height(axis_t::UNDERFLOW_BIN,j); 00308 height_5 += bin_height(axis_t::OVERFLOW_BIN,j); 00309 } 00310 00311 TW height_4 = 0; 00312 for(i=0;i<axis_x().bins();i++){ 00313 for(j=0;j<axis_y().bins();j++){ 00314 height_4 += bin_height(i,j); 00315 } 00316 } 00317 00318 a_out 00319 << " " << height_6 << " " << height_7 << " " << height_8 << std::endl; 00320 a_out 00321 << " " << height_3 << " " << height_4 << " " << height_5 << std::endl; 00322 a_out 00323 << " " << height_0 << " " << height_1 << " " << height_2 << std::endl; 00324 00325 // Some bins : 00326 bn_t xbins = axis_x().bins(); 00327 bn_t ybins = axis_y().bins(); 00328 a_out 00329 << " * ENTRIES[0,0] = " 00330 << bin_entries(0,0) 00331 << " * HEIGHT[0,0] = " 00332 << bin_height(0,0) 00333 << " * ERROR[0,0] = " 00334 << bin_error(0,0) 00335 << std::endl; 00336 a_out 00337 << " * ENTRIES[N/2,N/2] = " 00338 << bin_entries(xbins/2,ybins/2) 00339 << " * HEIGHT[N/2,N/2] = " 00340 << bin_height(xbins/2,ybins/2) 00341 << " * ERROR[N/2,N/2] = " 00342 << bin_error(xbins/2,ybins/2) 00343 << std::endl; 00344 a_out 00345 << " * ENTRIES[N-1,N-1] = " 00346 << bin_entries(xbins-1,ybins-1) 00347 << " * HEIGHT[N-1,N-1] = " 00348 << bin_height(xbins-1,ybins-1) 00349 << " * ERROR[N-1,N-1] = " 00350 << bin_error(xbins-1,ybins-1) 00351 << std::endl; 00352 } 00353 protected: 00354 b2(const std::string& a_title, 00355 bn_t aXnumber,TC aXmin,TC aXmax, 00356 bn_t aYnumber,TC aYmin,TC aYmax) 00357 :m_in_range_entries(0) 00358 ,m_in_range_Sw(0) 00359 ,m_in_range_Sxw(0) 00360 ,m_in_range_Syw(0) 00361 ,m_in_range_Sx2w(0) 00362 ,m_in_range_Sy2w(0) 00363 { 00364 parent::m_title = a_title; 00365 std::vector<bn_t> ns; 00366 ns.push_back(aXnumber); 00367 ns.push_back(aYnumber); 00368 std::vector<TC> mins; 00369 mins.push_back(aXmin); 00370 mins.push_back(aYmin); 00371 std::vector<TC> maxs; 00372 maxs.push_back(aXmax); 00373 maxs.push_back(aYmax); 00374 parent::configure(2,ns,mins,maxs); 00375 } 00376 00377 b2(const std::string& a_title, 00378 const std::vector<TC>& aEdgesX, 00379 const std::vector<TC>& aEdgesY) 00380 :m_in_range_entries(0) 00381 ,m_in_range_Sw(0) 00382 ,m_in_range_Sxw(0) 00383 ,m_in_range_Syw(0) 00384 ,m_in_range_Sx2w(0) 00385 ,m_in_range_Sy2w(0) 00386 { 00387 parent::m_title = a_title; 00388 std::vector< std::vector<TC> > edges(2); 00389 edges[0] = aEdgesX; 00390 edges[1] = aEdgesY; 00391 parent::configure(2,edges); 00392 } 00393 00394 virtual ~b2(){} 00395 protected: 00396 b2(const b2& a_from) 00397 : parent(a_from) 00398 ,m_in_range_entries(a_from.m_in_range_entries) 00399 ,m_in_range_Sw(a_from.m_in_range_Sw) 00400 ,m_in_range_Sxw(a_from.m_in_range_Sxw) 00401 ,m_in_range_Syw(a_from.m_in_range_Syw) 00402 ,m_in_range_Sx2w(a_from.m_in_range_Sx2w) 00403 ,m_in_range_Sy2w(a_from.m_in_range_Sy2w) 00404 { 00405 update_fast_getters(); 00406 } 00407 00408 b2& operator=(const b2& a_from) { 00409 parent::operator=(a_from); 00410 m_in_range_entries = a_from.m_in_range_entries; 00411 m_in_range_Sw = a_from.m_in_range_Sw; 00412 m_in_range_Sxw = a_from.m_in_range_Sxw; 00413 m_in_range_Syw = a_from.m_in_range_Syw; 00414 m_in_range_Sx2w = a_from.m_in_range_Sx2w; 00415 m_in_range_Sy2w = a_from.m_in_range_Sy2w; 00416 update_fast_getters(); 00417 return *this; 00418 } 00419 00420 protected: 00421 TN m_in_range_entries; 00422 TW m_in_range_Sw; 00423 TC m_in_range_Sxw; 00424 TC m_in_range_Syw; 00425 TC m_in_range_Sx2w; 00426 TC m_in_range_Sy2w; 00427 }; 00428 00429 }} 00430 00431 #endif 00432 00433 00434 00435