Warning, /include/Geant4/tools/histo/b3 is written in an unsupported language. File is not indexed.
0001 // Copyright (C) 2010, Guy Barrand. All rights reserved.
0002 // See the file tools.license for terms.
0003
0004 #ifndef tools_histo_b3
0005 #define tools_histo_b3
0006
0007 #include "base_histo"
0008
0009 #include <ostream>
0010
0011 namespace tools {
0012 namespace histo {
0013
0014 template <class TC,class TO,class TN,class TW,class TH>
0015 class b3 : public base_histo<TC,TO,TN,TW,TH> {
0016 typedef base_histo<TC,TO,TN,TW,TH> parent;
0017 public:
0018 typedef base_histo<TC,TO,TN,TW,TH> base_histo_t;
0019 typedef typename parent::axis_t axis_t;
0020 typedef typename parent::bn_t bn_t;
0021 protected:
0022 enum {AxisX=0,AxisY=1,AxisZ=2};
0023 public:
0024 virtual TH bin_error(int,int,int) const = 0; //for print
0025 public:
0026 // Partition :
0027 int coord_to_index_x(TC aCoord) const {
0028 return axis_x().coord_to_index(aCoord);
0029 }
0030 int coord_to_index_y(TC aCoord) const {
0031 return axis_y().coord_to_index(aCoord);
0032 }
0033 int coord_to_index_z(TC aCoord) const {
0034 return axis_z().coord_to_index(aCoord);
0035 }
0036
0037 TC mean_x() const {
0038 if(parent::m_in_range_Sw==0) return 0;
0039 return parent::m_in_range_Sxw[0]/parent::m_in_range_Sw;
0040 }
0041
0042 TC mean_y() const {
0043 if(parent::m_in_range_Sw==0) return 0;
0044 return parent::m_in_range_Sxw[1]/parent::m_in_range_Sw;
0045 }
0046
0047 TC mean_z() const {
0048 if(parent::m_in_range_Sw==0) return 0;
0049 return parent::m_in_range_Sxw[2]/parent::m_in_range_Sw;
0050 }
0051
0052 TC rms_x() const {
0053 if(parent::m_in_range_Sw==0) return 0;
0054 TC mean = parent::m_in_range_Sxw[0]/parent::m_in_range_Sw;
0055 return ::sqrt(::fabs((parent::m_in_range_Sx2w[0] / parent::m_in_range_Sw) - mean * mean));
0056 }
0057
0058 TC rms_y() const {
0059 if(parent::m_in_range_Sw==0) return 0;
0060 TC mean = parent::m_in_range_Sxw[1]/parent::m_in_range_Sw;
0061 return ::sqrt(::fabs((parent::m_in_range_Sx2w[1] / parent::m_in_range_Sw) - mean * mean));
0062 }
0063
0064 TC rms_z() const {
0065 if(parent::m_in_range_Sw==0) return 0;
0066 TC mean = parent::m_in_range_Sxw[2]/parent::m_in_range_Sw;
0067 return ::sqrt(::fabs((parent::m_in_range_Sx2w[2] / parent::m_in_range_Sw) - mean * mean));
0068 }
0069
0070 // bins :
0071 TN bin_entries(int aI,int aJ,int aK) const {
0072 TO offset;
0073 if(!_find_offset(aI,aJ,aK,offset)) return 0;
0074 return parent::m_bin_entries[offset];
0075 }
0076
0077 TH bin_height(int aI,int aJ,int aK) const {
0078 TO offset;
0079 if(!_find_offset(aI,aJ,aK,offset)) return 0;
0080 return this->get_bin_height(offset);
0081 }
0082
0083 TC bin_center_x(int aI) const {return parent::m_axes[0].bin_center(aI);}
0084 TC bin_center_y(int aJ) const {return parent::m_axes[1].bin_center(aJ);}
0085 TC bin_center_z(int aK) const {return parent::m_axes[2].bin_center(aK);}
0086
0087 TC bin_mean_x(int aI,int aJ,int aK) const {
0088 TO offset;
0089 if(!_find_offset(aI,aJ,aK,offset)) return 0;
0090 TW sw = parent::m_bin_Sw[offset];
0091 if(sw==0) return 0;
0092 return parent::m_bin_Sxw[offset][AxisX]/sw;
0093 }
0094
0095 TC bin_mean_y(int aI,int aJ,int aK) const {
0096 TO offset;
0097 if(!_find_offset(aI,aJ,aK,offset)) return 0;
0098 TW sw = parent::m_bin_Sw[offset];
0099 if(sw==0) return 0;
0100 return parent::m_bin_Sxw[offset][AxisY]/sw;
0101 }
0102
0103 TC bin_mean_z(int aI,int aJ,int aK) const {
0104 TO offset;
0105 if(!_find_offset(aI,aJ,aK,offset)) return 0;
0106 TW sw = parent::m_bin_Sw[offset];
0107 if(sw==0) return 0;
0108 return parent::m_bin_Sxw[offset][AxisZ]/sw;
0109 }
0110
0111 TC bin_rms_x(int aI,int aJ,int aK) const {
0112 TO offset;
0113 if(!_find_offset(aI,aJ,aK,offset)) return 0;
0114 TW sw = parent::m_bin_Sw[offset];
0115 if(sw==0) return 0;
0116 TC sxw = parent::m_bin_Sxw[offset][AxisX];
0117 TC sx2w = parent::m_bin_Sx2w[offset][AxisX];
0118 TC mean = sxw/sw;
0119 return ::sqrt(::fabs((sx2w / sw) - mean * mean));
0120 }
0121
0122 TC bin_rms_y(int aI,int aJ,int aK) const {
0123 TO offset;
0124 if(!_find_offset(aI,aJ,aK,offset)) return 0;
0125 TW sw = parent::m_bin_Sw[offset];
0126 if(sw==0) return 0;
0127 TC sxw = parent::m_bin_Sxw[offset][AxisY];
0128 TC sx2w = parent::m_bin_Sx2w[offset][AxisY];
0129 TC mean = sxw/sw;
0130 return ::sqrt(::fabs((sx2w / sw) - mean * mean));
0131 }
0132
0133 TC bin_rms_z(int aI,int aJ,int aK) const {
0134 TO offset;
0135 if(!_find_offset(aI,aJ,aK,offset)) return 0;
0136 TW sw = parent::m_bin_Sw[offset];
0137 if(sw==0) return 0;
0138 TC sxw = parent::m_bin_Sxw[offset][AxisZ];
0139 TC sx2w = parent::m_bin_Sx2w[offset][AxisZ];
0140 TC mean = sxw/sw;
0141 return ::sqrt(::fabs((sx2w / sw) - mean * mean));
0142 }
0143
0144 // Axes :
0145 const axis_t& axis_x() const {return parent::m_axes[0];}
0146 const axis_t& axis_y() const {return parent::m_axes[1];}
0147 const axis_t& axis_z() const {return parent::m_axes[2];}
0148 axis_t& axis_x() {return parent::m_axes[0];} //touchy
0149 axis_t& axis_y() {return parent::m_axes[1];} //touchy
0150 axis_t& axis_z() {return parent::m_axes[2];} //touchy
0151
0152 // Projection :
0153 TN bin_entries_x(int aI) const {
0154 if(!parent::m_dimension) return 0;
0155 bn_t ibin;
0156 if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) return 0;
0157 bn_t jbin,kbin,offset;
0158 bn_t ybins = parent::m_axes[1].bins()+2;
0159 bn_t zbins = parent::m_axes[2].bins()+2;
0160 TO yoffset = parent::m_axes[1].m_offset;
0161 TO zoffset = parent::m_axes[2].m_offset;
0162 TO joffset = ibin;
0163 TN _entries = 0;
0164 for(jbin=0;jbin<ybins;jbin++) {
0165 offset = joffset;
0166 for(kbin=0;kbin<zbins;kbin++) {
0167 _entries += parent::m_bin_entries[offset];
0168 offset += zoffset;
0169 }
0170 joffset += yoffset;
0171 }
0172 return _entries;
0173 }
0174
0175 TN bin_entries_y(int aJ) const {
0176 if(!parent::m_dimension) return 0;
0177 bn_t jbin;
0178 if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) return 0;
0179 bn_t ibin,kbin;
0180 TO offset;
0181 bn_t xbins = parent::m_axes[0].bins()+2;
0182 bn_t zbins = parent::m_axes[2].bins()+2;
0183 TO yoffset = parent::m_axes[1].m_offset;
0184 TO zoffset = parent::m_axes[2].m_offset;
0185 TO joffset = jbin * yoffset;
0186 TN _entries = 0;
0187 for(ibin=0;ibin<xbins;ibin++) {
0188 offset = joffset;
0189 for(kbin=0;kbin<zbins;kbin++) {
0190 _entries += parent::m_bin_entries[offset];
0191 offset += zoffset;
0192 }
0193 joffset++;
0194 }
0195 return _entries;
0196 }
0197
0198 TN bin_entries_z(int aK) const {
0199 if(!parent::m_dimension) return 0;
0200 bn_t kbin;
0201 if(!parent::m_axes[2].in_range_to_absolute_index(aK,kbin)) return 0;
0202 bn_t ibin,jbin;
0203 TO offset;
0204 bn_t xbins = parent::m_axes[0].bins()+2;
0205 bn_t ybins = parent::m_axes[1].bins()+2;
0206 TO yoffset = parent::m_axes[1].m_offset;
0207 TO zoffset = parent::m_axes[2].m_offset;
0208 TO koffset = kbin * zoffset;
0209 TN _entries = 0;
0210 for(ibin=0;ibin<xbins;ibin++) {
0211 offset = koffset;
0212 for(jbin=0;jbin<ybins;jbin++) {
0213 _entries += parent::m_bin_entries[offset];
0214 offset += yoffset;
0215 }
0216 koffset++;
0217 }
0218 return _entries;
0219 }
0220
0221 TW bin_height_x(int aI) const {
0222 //to slow : return get_ith_axis_bin_height(0,aI);
0223 if(!parent::m_dimension) return 0;
0224 bn_t ibin;
0225 if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) return 0;
0226 bn_t ybins = parent::m_axes[1].bins()+2;
0227 bn_t zbins = parent::m_axes[2].bins()+2;
0228 TO yoffset = parent::m_axes[1].m_offset;
0229 TO zoffset = parent::m_axes[2].m_offset;
0230 TO joffset = ibin;
0231 TW sw = 0;
0232 for(bn_t jbin=0;jbin<ybins;jbin++) {
0233 TO offset = joffset;
0234 for(bn_t kbin=0;kbin<zbins;kbin++) {
0235 sw += this->get_bin_height(offset);
0236 offset += zoffset;
0237 }
0238 joffset += yoffset;
0239 }
0240 return sw;
0241 }
0242
0243 TW bin_height_y(int aJ) const {
0244 if(!parent::m_dimension) return 0;
0245 bn_t jbin;
0246 if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) return 0;
0247 bn_t xbins = parent::m_axes[0].bins()+2;
0248 bn_t zbins = parent::m_axes[2].bins()+2;
0249 TO yoffset = parent::m_axes[1].m_offset;
0250 TO zoffset = parent::m_axes[2].m_offset;
0251 TO joffset = jbin * yoffset;
0252 TW sw = 0;
0253 for(bn_t ibin=0;ibin<xbins;ibin++) {
0254 TO offset = joffset;
0255 for(bn_t kbin=0;kbin<zbins;kbin++) {
0256 sw += this->get_bin_height(offset);
0257 offset += zoffset;
0258 }
0259 joffset++;
0260 }
0261 return sw;
0262 }
0263
0264 TW bin_height_z(int aK) const {
0265 if(!parent::m_dimension) return 0;
0266 bn_t kbin;
0267 if(!parent::m_axes[2].in_range_to_absolute_index(aK,kbin)) return 0;
0268 bn_t xbins = parent::m_axes[0].bins()+2;
0269 bn_t ybins = parent::m_axes[1].bins()+2;
0270 TO yoffset = parent::m_axes[1].m_offset;
0271 TO zoffset = parent::m_axes[2].m_offset;
0272 TO koffset = kbin * zoffset;
0273 TW sw = 0;
0274 for(bn_t ibin=0;ibin<xbins;ibin++) {
0275 TO offset = koffset;
0276 for(bn_t jbin=0;jbin<ybins;jbin++) {
0277 sw += this->get_bin_height(offset);
0278 offset += yoffset;
0279 }
0280 koffset++;
0281 }
0282 return sw;
0283 }
0284
0285 TC Sxyw() const {return parent::m_in_range_plane_Sxyw[0];}
0286 TC Syzw() const {return parent::m_in_range_plane_Sxyw[1];}
0287 TC Szxw() const {return parent::m_in_range_plane_Sxyw[2];}
0288 public:
0289 //NOTE : print is a Python keyword.
0290 void hprint(std::ostream& a_out) {
0291 // A la HPRINT.
0292 a_out << parent::dimension() << parent::title() << std::endl;
0293 a_out
0294 << " * ENTRIES = " << parent::all_entries() << std::endl;
0295
0296 }
0297 public:
0298 b3(const std::string& a_title,
0299 bn_t aXnumber,TC aXmin,TC aXmax,
0300 bn_t aYnumber,TC aYmin,TC aYmax,
0301 bn_t aZnumber,TC aZmin,TC aZmax)
0302 {
0303 parent::m_title = a_title;
0304 std::vector<bn_t> nbins;
0305 nbins.push_back(aXnumber);
0306 nbins.push_back(aYnumber);
0307 nbins.push_back(aZnumber);
0308 std::vector<TC> mins;
0309 mins.push_back(aXmin);
0310 mins.push_back(aYmin);
0311 mins.push_back(aZmin);
0312 std::vector<TC> maxs;
0313 maxs.push_back(aXmax);
0314 maxs.push_back(aYmax);
0315 maxs.push_back(aZmax);
0316 parent::configure(3,nbins,mins,maxs);
0317 }
0318
0319 b3(const std::string& a_title,
0320 const std::vector<TC>& a_edges_x,
0321 const std::vector<TC>& a_edges_y,
0322 const std::vector<TC>& a_edges_z)
0323 {
0324 parent::m_title = a_title;
0325 std::vector< std::vector<TC> > edges(3);
0326 edges[0] = a_edges_x;
0327 edges[1] = a_edges_y;
0328 edges[2] = a_edges_z;
0329 parent::configure(3,edges);
0330 }
0331
0332 virtual ~b3(){}
0333 protected:
0334 b3(const b3& a_from):parent(a_from) {}
0335 b3& operator=(const b3& a_from){parent::operator=(a_from);return *this;}
0336 public:
0337 bool configure(bn_t aXnumber,TC aXmin,TC aXmax,
0338 bn_t aYnumber,TC aYmin,TC aYmax,
0339 bn_t aZnumber,TC aZmin,TC aZmax){
0340 std::vector<bn_t> nbins;
0341 nbins.push_back(aXnumber);
0342 nbins.push_back(aYnumber);
0343 nbins.push_back(aZnumber);
0344 std::vector<TC> mins;
0345 mins.push_back(aXmin);
0346 mins.push_back(aYmin);
0347 mins.push_back(aZmin);
0348 std::vector<TC> maxs;
0349 maxs.push_back(aXmax);
0350 maxs.push_back(aYmax);
0351 maxs.push_back(aZmax);
0352 return parent::configure(3,nbins,mins,maxs);
0353 }
0354
0355 bool configure(const std::vector<TC>& a_edges_x,
0356 const std::vector<TC>& a_edges_y,
0357 const std::vector<TC>& a_edges_z){
0358 std::vector< std::vector<TC> > edges(3);
0359 edges[0] = a_edges_x;
0360 edges[1] = a_edges_y;
0361 edges[2] = a_edges_z;
0362 return parent::configure(3,edges);
0363 }
0364
0365 protected:
0366 bool _find_offset(int aI,int aJ,int aK,TO& a_offset) const {
0367 if(parent::m_dimension!=3) {a_offset=0;return false;}
0368 bn_t ibin,jbin,kbin;
0369 if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) {a_offset=0;return false;}
0370 if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) {a_offset=0;return false;}
0371 if(!parent::m_axes[2].in_range_to_absolute_index(aK,kbin)) {a_offset=0;return false;}
0372 a_offset = ibin + jbin * parent::m_axes[1].m_offset + kbin * parent::m_axes[2].m_offset;
0373 return true;
0374 }
0375 };
0376
0377 }}
0378
0379 #endif
0380
0381
0382
0383
