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 //joffset = ibin + jbin * parent::m_axes[1].m_offset;
0166 offset = joffset;
0167 for(kbin=0;kbin<zbins;kbin++) {
0168 //offset = joffset + kbin * parent::m_axes[2].m_offset;
0169 _entries += parent::m_bin_entries[offset];
0170 offset += zoffset;
0171 }
0172 joffset += yoffset;
0173 }
0174 return _entries;
0175 }
0176
0177 TN bin_entries_y(int aJ) const {
0178 if(!parent::m_dimension) return 0;
0179 bn_t jbin;
0180 if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) return 0;
0181 bn_t ibin,kbin;
0182 TO offset;
0183 bn_t xbins = parent::m_axes[0].bins()+2;
0184 bn_t zbins = parent::m_axes[2].bins()+2;
0185 TO yoffset = parent::m_axes[1].m_offset;
0186 TO zoffset = parent::m_axes[2].m_offset;
0187 TO joffset = jbin * yoffset;
0188 TN _entries = 0;
0189 for(ibin=0;ibin<xbins;ibin++) {
0190 //joffset = ibin + jbin * parent::m_axes[1].m_offset;
0191 offset = joffset;
0192 for(kbin=0;kbin<zbins;kbin++) {
0193 //offset = joffset + kbin * parent::m_axes[2].m_offset;
0194 _entries += parent::m_bin_entries[offset];
0195 offset += zoffset;
0196 }
0197 joffset++;
0198 }
0199 return _entries;
0200 }
0201
0202 TN bin_entries_z(int aK) const {
0203 if(!parent::m_dimension) return 0;
0204 bn_t kbin;
0205 if(!parent::m_axes[2].in_range_to_absolute_index(aK,kbin)) return 0;
0206 bn_t ibin,jbin;
0207 TO offset;
0208 bn_t xbins = parent::m_axes[0].bins()+2;
0209 bn_t ybins = parent::m_axes[1].bins()+2;
0210 TO yoffset = parent::m_axes[1].m_offset;
0211 TO zoffset = parent::m_axes[2].m_offset;
0212 TO koffset = kbin * zoffset;
0213 TN _entries = 0;
0214 for(ibin=0;ibin<xbins;ibin++) {
0215 //koffset = ibin + kbin * parent::m_axes[2].m_offset;
0216 offset = koffset;
0217 for(jbin=0;jbin<ybins;jbin++) {
0218 //offset = koffset + jbin * parent::m_axes[1].m_offset;
0219 _entries += parent::m_bin_entries[offset];
0220 offset += yoffset;
0221 }
0222 koffset++;
0223 }
0224 return _entries;
0225 }
0226
0227 TW bin_height_x(int aI) const {
0228 //to slow : return get_ith_axis_bin_height(0,aI);
0229 if(!parent::m_dimension) return 0;
0230 bn_t ibin;
0231 if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) return 0;
0232 bn_t ybins = parent::m_axes[1].bins()+2;
0233 bn_t zbins = parent::m_axes[2].bins()+2;
0234 TO yoffset = parent::m_axes[1].m_offset;
0235 TO zoffset = parent::m_axes[2].m_offset;
0236 TO joffset = ibin;
0237 TW sw = 0;
0238 for(bn_t jbin=0;jbin<ybins;jbin++) {
0239 //joffset = ibin + jbin * parent::m_axes[1].m_offset;
0240 TO offset = joffset;
0241 for(bn_t kbin=0;kbin<zbins;kbin++) {
0242 //offset = joffset + kbin * parent::m_axes[2].m_offset;
0243 sw += this->get_bin_height(offset);
0244 offset += zoffset;
0245 }
0246 joffset += yoffset;
0247 }
0248 return sw;
0249 }
0250
0251 TW bin_height_y(int aJ) const {
0252 if(!parent::m_dimension) return 0;
0253 bn_t jbin;
0254 if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) return 0;
0255 bn_t xbins = parent::m_axes[0].bins()+2;
0256 bn_t zbins = parent::m_axes[2].bins()+2;
0257 TO yoffset = parent::m_axes[1].m_offset;
0258 TO zoffset = parent::m_axes[2].m_offset;
0259 TO joffset = jbin * yoffset;
0260 TW sw = 0;
0261 for(bn_t ibin=0;ibin<xbins;ibin++) {
0262 //joffset = ibin + jbin * parent::m_axes[1].m_offset;
0263 TO offset = joffset;
0264 for(bn_t kbin=0;kbin<zbins;kbin++) {
0265 //offset = joffset + kbin * parent::m_axes[2].m_offset;
0266 sw += this->get_bin_height(offset);
0267 offset += zoffset;
0268 }
0269 joffset++;
0270 }
0271 return sw;
0272 }
0273
0274 TW bin_height_z(int aK) const {
0275 if(!parent::m_dimension) return 0;
0276 bn_t kbin;
0277 if(!parent::m_axes[2].in_range_to_absolute_index(aK,kbin)) return 0;
0278 bn_t xbins = parent::m_axes[0].bins()+2;
0279 bn_t ybins = parent::m_axes[1].bins()+2;
0280 TO yoffset = parent::m_axes[1].m_offset;
0281 TO zoffset = parent::m_axes[2].m_offset;
0282 TO koffset = kbin * zoffset;
0283 TW sw = 0;
0284 for(bn_t ibin=0;ibin<xbins;ibin++) {
0285 //koffset = ibin + kbin * parent::m_axes[2].m_offset;
0286 TO offset = koffset;
0287 for(bn_t jbin=0;jbin<ybins;jbin++) {
0288 //offset = koffset + jbin * parent::m_axes[1].m_offset;
0289 sw += this->get_bin_height(offset);
0290 offset += yoffset;
0291 }
0292 koffset++;
0293 }
0294 return sw;
0295 }
0296
0297 TC Sxyw() const {return parent::m_in_range_plane_Sxyw[0];}
0298 TC Syzw() const {return parent::m_in_range_plane_Sxyw[1];}
0299 TC Szxw() const {return parent::m_in_range_plane_Sxyw[2];}
0300 public:
0301 //NOTE : print is a Python keyword.
0302 void hprint(std::ostream& a_out) {
0303 // A la HPRINT.
0304 a_out << parent::dimension() << parent::title() << std::endl;
0305 a_out
0306 << " * ENTRIES = " << parent::all_entries() << std::endl;
0307
0308 }
0309 public:
0310 b3(const std::string& a_title,
0311 bn_t aXnumber,TC aXmin,TC aXmax,
0312 bn_t aYnumber,TC aYmin,TC aYmax,
0313 bn_t aZnumber,TC aZmin,TC aZmax)
0314 {
0315 parent::m_title = a_title;
0316 std::vector<bn_t> nbins;
0317 nbins.push_back(aXnumber);
0318 nbins.push_back(aYnumber);
0319 nbins.push_back(aZnumber);
0320 std::vector<TC> mins;
0321 mins.push_back(aXmin);
0322 mins.push_back(aYmin);
0323 mins.push_back(aZmin);
0324 std::vector<TC> maxs;
0325 maxs.push_back(aXmax);
0326 maxs.push_back(aYmax);
0327 maxs.push_back(aZmax);
0328 parent::configure(3,nbins,mins,maxs);
0329 }
0330
0331 b3(const std::string& a_title,
0332 const std::vector<TC>& a_edges_x,
0333 const std::vector<TC>& a_edges_y,
0334 const std::vector<TC>& a_edges_z)
0335 {
0336 parent::m_title = a_title;
0337 std::vector< std::vector<TC> > edges(3);
0338 edges[0] = a_edges_x;
0339 edges[1] = a_edges_y;
0340 edges[2] = a_edges_z;
0341 parent::configure(3,edges);
0342 }
0343
0344 virtual ~b3(){}
0345 protected:
0346 b3(const b3& a_from):parent(a_from) {}
0347 b3& operator=(const b3& a_from){parent::operator=(a_from);return *this;}
0348 public:
0349 bool configure(bn_t aXnumber,TC aXmin,TC aXmax,
0350 bn_t aYnumber,TC aYmin,TC aYmax,
0351 bn_t aZnumber,TC aZmin,TC aZmax){
0352 std::vector<bn_t> nbins;
0353 nbins.push_back(aXnumber);
0354 nbins.push_back(aYnumber);
0355 nbins.push_back(aZnumber);
0356 std::vector<TC> mins;
0357 mins.push_back(aXmin);
0358 mins.push_back(aYmin);
0359 mins.push_back(aZmin);
0360 std::vector<TC> maxs;
0361 maxs.push_back(aXmax);
0362 maxs.push_back(aYmax);
0363 maxs.push_back(aZmax);
0364 return parent::configure(3,nbins,mins,maxs);
0365 }
0366
0367 bool configure(const std::vector<TC>& a_edges_x,
0368 const std::vector<TC>& a_edges_y,
0369 const std::vector<TC>& a_edges_z){
0370 std::vector< std::vector<TC> > edges(3);
0371 edges[0] = a_edges_x;
0372 edges[1] = a_edges_y;
0373 edges[2] = a_edges_z;
0374 return parent::configure(3,edges);
0375 }
0376
0377 protected:
0378 bool _find_offset(int aI,int aJ,int aK,TO& a_offset) const {
0379 if(parent::m_dimension!=3) {a_offset=0;return false;}
0380 bn_t ibin,jbin,kbin;
0381 if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) {a_offset=0;return false;}
0382 if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) {a_offset=0;return false;}
0383 if(!parent::m_axes[2].in_range_to_absolute_index(aK,kbin)) {a_offset=0;return false;}
0384 a_offset = ibin + jbin * parent::m_axes[1].m_offset + kbin * parent::m_axes[2].m_offset;
0385 return true;
0386 }
0387 };
0388
0389 }}
0390
0391 #endif
0392
0393
0394
0395
