// Copyright (C) 2010, Guy Barrand. All rights reserved. // See the file tools.license for terms. #ifndef tools_histo_b3 #define tools_histo_b3 #include "base_histo" #include namespace tools { namespace histo { template class b3 : public base_histo { typedef base_histo parent; public: typedef axis axis_t; typedef typename base_histo::bn_t bn_t; protected: enum {AxisX=0,AxisY=1,AxisZ=2}; public: virtual TH bin_error(int,int,int) const = 0; //for print public: void update_fast_getters() { m_in_range_entries = 0; m_in_range_Sw = 0; m_in_range_Sxw = 0; m_in_range_Syw = 0; m_in_range_Szw = 0; m_in_range_Sx2w = 0; m_in_range_Sy2w = 0; m_in_range_Sz2w = 0; bn_t ibin,jbin,kbin,joffset,offset; bn_t xbins = parent::m_axes[0].bins(); bn_t ybins = parent::m_axes[1].bins(); bn_t zbins = parent::m_axes[2].bins(); bn_t yoffset = parent::m_axes[1].m_offset; bn_t zoffset = parent::m_axes[2].m_offset; for(ibin=1;ibin<=xbins;ibin++) { joffset = ibin + yoffset; for(jbin=1;jbin<=ybins;jbin++) { //joffset = ibin + jbin * parent::m_axes[1].m_offset; offset = joffset + zoffset; for(kbin=1;kbin<=zbins;kbin++) { //offset = joffset + kbin * parent::m_axes[2].m_offset; m_in_range_entries += parent::m_bin_entries[offset]; m_in_range_Sw += parent::m_bin_Sw[offset]; m_in_range_Sxw += parent::m_bin_Sxw[offset][0]; m_in_range_Syw += parent::m_bin_Sxw[offset][1]; m_in_range_Szw += parent::m_bin_Sxw[offset][2]; m_in_range_Sx2w += parent::m_bin_Sx2w[offset][0]; m_in_range_Sy2w += parent::m_bin_Sx2w[offset][1]; m_in_range_Sz2w += parent::m_bin_Sx2w[offset][2]; offset += zoffset; } joffset += yoffset; } } } // Partition : int coord_to_index_x(TC aCoord) const { return axis_x().coord_to_index(aCoord); } int coord_to_index_y(TC aCoord) const { return axis_y().coord_to_index(aCoord); } int coord_to_index_z(TC aCoord) const { return axis_z().coord_to_index(aCoord); } TC mean_x() const { if(m_in_range_Sw==0) return 0; return m_in_range_Sxw/m_in_range_Sw; } TC mean_y() const { if(m_in_range_Sw==0) return 0; return m_in_range_Syw/m_in_range_Sw; } TC mean_z() const { if(m_in_range_Sw==0) return 0; return m_in_range_Szw/m_in_range_Sw; } TC rms_x() const { if(m_in_range_Sw==0) return 0; TC mean = m_in_range_Sxw/m_in_range_Sw; return ::sqrt(::fabs((m_in_range_Sx2w / m_in_range_Sw) - mean * mean)); } TC rms_y() const { if(m_in_range_Sw==0) return 0; TC mean = m_in_range_Syw/m_in_range_Sw; return ::sqrt(::fabs((m_in_range_Sy2w / m_in_range_Sw) - mean * mean)); } TC rms_z() const { if(m_in_range_Sw==0) return 0; TC mean = m_in_range_Szw/m_in_range_Sw; return ::sqrt(::fabs((m_in_range_Sz2w / m_in_range_Sw) - mean * mean)); } // bins : TN bin_entries(int aI,int aJ,int aK) const { if(parent::m_bin_number==0) return 0; bn_t ibin,jbin,kbin; if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) return 0; if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) return 0; if(!parent::m_axes[2].in_range_to_absolute_index(aK,kbin)) return 0; bn_t offset = ibin + jbin * parent::m_axes[1].m_offset + kbin * parent::m_axes[2].m_offset; return parent::m_bin_entries[offset]; } TH bin_height(int aI,int aJ,int aK) const { if(parent::m_bin_number==0) return 0; bn_t ibin,jbin,kbin; if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) return 0; if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) return 0; if(!parent::m_axes[2].in_range_to_absolute_index(aK,kbin)) return 0; bn_t offset = ibin + jbin * parent::m_axes[1].m_offset + kbin * parent::m_axes[2].m_offset; return this->get_bin_height(offset); } TC bin_center_x(int aI) const {return parent::m_axes[0].bin_center(aI);} TC bin_center_y(int aJ) const {return parent::m_axes[1].bin_center(aJ);} TC bin_center_z(int aK) const {return parent::m_axes[2].bin_center(aK);} TC bin_mean_x(int aI,int aJ,int aK) const { if(parent::m_bin_number==0) return 0; bn_t ibin,jbin,kbin; if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) return 0; if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) return 0; if(!parent::m_axes[2].in_range_to_absolute_index(aK,kbin)) return 0; bn_t offset = ibin + jbin * parent::m_axes[1].m_offset + kbin * parent::m_axes[2].m_offset; TW sw = parent::m_bin_Sw[offset]; if(sw==0) return 0; return parent::m_bin_Sxw[offset][AxisX]/sw; } TC bin_mean_y(int aI,int aJ,int aK) const { if(parent::m_bin_number==0) return 0; bn_t ibin,jbin,kbin; if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) return 0; if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) return 0; if(!parent::m_axes[2].in_range_to_absolute_index(aK,kbin)) return 0; bn_t offset = ibin + jbin * parent::m_axes[1].m_offset + kbin * parent::m_axes[2].m_offset; TW sw = parent::m_bin_Sw[offset]; if(sw==0) return 0; return parent::m_bin_Sxw[offset][AxisY]/sw; } TC bin_mean_z(int aI,int aJ,int aK) const { if(parent::m_bin_number==0) return 0; bn_t ibin,jbin,kbin; if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) return 0; if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) return 0; if(!parent::m_axes[2].in_range_to_absolute_index(aK,kbin)) return 0; bn_t offset = ibin + jbin * parent::m_axes[1].m_offset + kbin * parent::m_axes[2].m_offset; TW sw = parent::m_bin_Sw[offset]; if(sw==0) return 0; return parent::m_bin_Sxw[offset][AxisZ]/sw; } TC bin_rms_x(int aI,int aJ,int aK) const { if(parent::m_bin_number==0) return 0; bn_t ibin,jbin,kbin; if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) return 0; if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) return 0; if(!parent::m_axes[2].in_range_to_absolute_index(aK,kbin)) return 0; bn_t offset = ibin + jbin * parent::m_axes[1].m_offset + kbin * parent::m_axes[2].m_offset; TW sw = parent::m_bin_Sw[offset]; if(sw==0) return 0; TC sxw = parent::m_bin_Sxw[offset][AxisX]; TC sx2w = parent::m_bin_Sx2w[offset][AxisX]; TC mean = sxw/sw; return ::sqrt(::fabs((sx2w / sw) - mean * mean)); } TC bin_rms_y(int aI,int aJ,int aK) const { if(parent::m_bin_number==0) return 0; bn_t ibin,jbin,kbin; if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) return 0; if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) return 0; if(!parent::m_axes[2].in_range_to_absolute_index(aK,kbin)) return 0; bn_t offset = ibin + jbin * parent::m_axes[1].m_offset + kbin * parent::m_axes[2].m_offset; TW sw = parent::m_bin_Sw[offset]; if(sw==0) return 0; TC sxw = parent::m_bin_Sxw[offset][AxisY]; TC sx2w = parent::m_bin_Sx2w[offset][AxisY]; TC mean = sxw/sw; return ::sqrt(::fabs((sx2w / sw) - mean * mean)); } TC bin_rms_z(int aI,int aJ,int aK) const { if(parent::m_bin_number==0) return 0; bn_t ibin,jbin,kbin; if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) return 0; if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) return 0; if(!parent::m_axes[2].in_range_to_absolute_index(aK,kbin)) return 0; bn_t offset = ibin + jbin * parent::m_axes[1].m_offset + kbin * parent::m_axes[2].m_offset; TW sw = parent::m_bin_Sw[offset]; if(sw==0) return 0; TC sxw = parent::m_bin_Sxw[offset][AxisZ]; TC sx2w = parent::m_bin_Sx2w[offset][AxisZ]; TC mean = sxw/sw; return ::sqrt(::fabs((sx2w / sw) - mean * mean)); } // Axes : const axis_t& axis_x() const {return parent::m_axes[0];} const axis_t& axis_y() const {return parent::m_axes[1];} const axis_t& axis_z() const {return parent::m_axes[2];} axis_t& axis_x() {return parent::m_axes[0];} //touchy axis_t& axis_y() {return parent::m_axes[1];} //touchy axis_t& axis_z() {return parent::m_axes[2];} //touchy // Projection : TN bin_entries_x(int aI) const { if(!parent::m_dimension) return 0; bn_t ibin; if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) return 0; bn_t jbin,kbin,offset; bn_t ybins = parent::m_axes[1].bins()+2; bn_t zbins = parent::m_axes[2].bins()+2; bn_t yoffset = parent::m_axes[1].m_offset; bn_t zoffset = parent::m_axes[2].m_offset; bn_t joffset = ibin; TN _entries = 0; for(jbin=0;jbinget_bin_height(offset); offset += zoffset; } joffset += yoffset; } return sw; } TW bin_height_y(int aJ) const { if(!parent::m_dimension) return 0; bn_t jbin; if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) return 0; bn_t xbins = parent::m_axes[0].bins()+2; bn_t zbins = parent::m_axes[2].bins()+2; bn_t yoffset = parent::m_axes[1].m_offset; bn_t zoffset = parent::m_axes[2].m_offset; bn_t joffset = jbin * yoffset; TW sw = 0; for(bn_t ibin=0;ibinget_bin_height(offset); offset += zoffset; } joffset++; } return sw; } TW bin_height_z(int aK) const { if(!parent::m_dimension) return 0; bn_t kbin; if(!parent::m_axes[2].in_range_to_absolute_index(aK,kbin)) return 0; bn_t xbins = parent::m_axes[0].bins()+2; bn_t ybins = parent::m_axes[1].bins()+2; bn_t yoffset = parent::m_axes[1].m_offset; bn_t zoffset = parent::m_axes[2].m_offset; bn_t koffset = kbin * zoffset; TW sw = 0; for(bn_t ibin=0;ibinget_bin_height(offset); offset += yoffset; } koffset++; } return sw; } public: //NOTE : print is a Python keyword. void hprint(std::ostream& a_out) { // A la HPRINT. a_out << parent::dimension() << parent::title() << std::endl; a_out << " * ENTRIES = " << parent::all_entries() << std::endl; } public: b3(const std::string& a_title, bn_t aXnumber,TC aXmin,TC aXmax, bn_t aYnumber,TC aYmin,TC aYmax, bn_t aZnumber,TC aZmin,TC aZmax) :m_in_range_entries(0) ,m_in_range_Sw(0) ,m_in_range_Sxw(0) ,m_in_range_Syw(0) ,m_in_range_Szw(0) ,m_in_range_Sx2w(0) ,m_in_range_Sy2w(0) ,m_in_range_Sz2w(0) { parent::m_title = a_title; std::vector nbins; nbins.push_back(aXnumber); nbins.push_back(aYnumber); nbins.push_back(aZnumber); std::vector mins; mins.push_back(aXmin); mins.push_back(aYmin); mins.push_back(aZmin); std::vector maxs; maxs.push_back(aXmax); maxs.push_back(aYmax); maxs.push_back(aZmax); parent::configure(3,nbins,mins,maxs); } b3(const std::string& a_title, const std::vector& aEdgesX, const std::vector& aEdgesY, const std::vector& aEdgesZ) :m_in_range_entries(0) ,m_in_range_Sw(0) ,m_in_range_Sxw(0) ,m_in_range_Syw(0) ,m_in_range_Szw(0) ,m_in_range_Sx2w(0) ,m_in_range_Sy2w(0) ,m_in_range_Sz2w(0) { parent::m_title = a_title; std::vector< std::vector > edges(3); edges[0] = aEdgesX; edges[1] = aEdgesY; edges[2] = aEdgesZ; parent::configure(3,edges); } virtual ~b3(){} protected: b3(const b3& a_from) : parent(a_from) ,m_in_range_entries(a_from.m_in_range_entries) ,m_in_range_Sw(a_from.m_in_range_Sw) ,m_in_range_Sxw(a_from.m_in_range_Sxw) ,m_in_range_Syw(a_from.m_in_range_Syw) ,m_in_range_Szw(a_from.m_in_range_Szw) ,m_in_range_Sx2w(a_from.m_in_range_Sx2w) ,m_in_range_Sy2w(a_from.m_in_range_Sy2w) ,m_in_range_Sz2w(a_from.m_in_range_Sz2w) { update_fast_getters(); } b3& operator=(const b3& a_from){ parent::operator=(a_from); m_in_range_entries = a_from.m_in_range_entries; m_in_range_Sw = a_from.m_in_range_Sw; m_in_range_Sxw = a_from.m_in_range_Sxw; m_in_range_Syw = a_from.m_in_range_Syw; m_in_range_Szw = a_from.m_in_range_Szw; m_in_range_Sx2w = a_from.m_in_range_Sx2w; m_in_range_Sy2w = a_from.m_in_range_Sy2w; m_in_range_Sz2w = a_from.m_in_range_Sz2w; update_fast_getters(); return *this; } public: bool configure(bn_t aXnumber,TC aXmin,TC aXmax, bn_t aYnumber,TC aYmin,TC aYmax, bn_t aZnumber,TC aZmin,TC aZmax){ m_in_range_entries = 0; m_in_range_Sw = 0; m_in_range_Sxw = 0; m_in_range_Syw = 0; m_in_range_Szw = 0; m_in_range_Sx2w = 0; m_in_range_Sy2w = 0; m_in_range_Sz2w = 0; std::vector nbins; nbins.push_back(aXnumber); nbins.push_back(aYnumber); nbins.push_back(aZnumber); std::vector mins; mins.push_back(aXmin); mins.push_back(aYmin); mins.push_back(aZmin); std::vector maxs; maxs.push_back(aXmax); maxs.push_back(aYmax); maxs.push_back(aZmax); return parent::configure(3,nbins,mins,maxs); } bool configure(const std::vector& aEdgesX, const std::vector& aEdgesY, const std::vector& aEdgesZ){ m_in_range_entries = 0; m_in_range_Sw = 0; m_in_range_Sxw = 0; m_in_range_Syw = 0; m_in_range_Szw = 0; m_in_range_Sx2w = 0; m_in_range_Sy2w = 0; m_in_range_Sz2w = 0; std::vector< std::vector > edges(3); edges[0] = aEdgesX; edges[1] = aEdgesY; edges[2] = aEdgesZ; return parent::configure(3,edges); } protected: TN m_in_range_entries; TW m_in_range_Sw; TC m_in_range_Sxw; TC m_in_range_Syw; TC m_in_range_Szw; TC m_in_range_Sx2w; TC m_in_range_Sy2w; TC m_in_range_Sz2w; }; }} #endif