Warning, /include/Geant4/tools/histo/h3 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_h3
0005 #define tools_histo_h3
0006
0007 #include "b3"
0008
0009 namespace tools {
0010 namespace histo {
0011
0012 template <class TC,class TO,class TN,class TW,class TH>
0013 class h3 : public b3<TC,TO,TN,TW,TH> {
0014 typedef b3<TC,TO,TN,TW,TH> parent;
0015 public:
0016 typedef histo_data<TC,TO,TN,TW> hd_t;
0017 typedef typename parent::bn_t bn_t;
0018 protected:
0019 virtual TH get_bin_height(TO a_offset) const { //TH should be the same as TW
0020 return parent::m_bin_Sw[a_offset];
0021 }
0022 public:
0023
0024 virtual TH bin_error(int aI,int aJ,int aK) const {
0025 TO offset;
0026 if(!parent::_find_offset(aI,aJ,aK,offset)) return 0;
0027 return ::sqrt(parent::m_bin_Sw2[offset]);
0028 }
0029
0030 public:
0031 bool multiply(TW a_factor){return parent::base_multiply(a_factor);}
0032 bool scale(TW a_factor) {return multiply(a_factor);}
0033
0034 void copy_from_data(const hd_t& a_from) {parent::base_from_data(a_from);}
0035 hd_t get_histo_data() const {return *this;} //deprecated. Keep it for g4tools.
0036
0037 bool reset() {
0038 parent::base_reset();
0039 return true;
0040 }
0041
0042 bool fill(TC aX,TC aY,TC aZ,TW aWeight = 1) {
0043 if(parent::m_dimension!=3) return false;
0044
0045 bn_t ibin,jbin,kbin;
0046 if(!parent::m_axes[0].coord_to_absolute_index(aX,ibin)) return false;
0047 if(!parent::m_axes[1].coord_to_absolute_index(aY,jbin)) return false;
0048 if(!parent::m_axes[2].coord_to_absolute_index(aZ,kbin)) return false;
0049 TO offset = ibin + jbin * parent::m_axes[1].m_offset + kbin * parent::m_axes[2].m_offset;
0050
0051 parent::m_bin_entries[offset]++;
0052 parent::m_bin_Sw[offset] += aWeight;
0053 parent::m_bin_Sw2[offset] += aWeight * aWeight;
0054
0055 TC xw = aX * aWeight;
0056 TC x2w = aX * xw;
0057 parent::m_bin_Sxw[offset][0] += xw;
0058 parent::m_bin_Sx2w[offset][0] += x2w;
0059
0060 TC yw = aY * aWeight;
0061 TC y2w = aY * yw;
0062 parent::m_bin_Sxw[offset][1] += yw;
0063 parent::m_bin_Sx2w[offset][1] += y2w;
0064
0065 TC zw = aZ * aWeight;
0066 TC z2w = aZ * zw;
0067 parent::m_bin_Sxw[offset][2] += zw;
0068 parent::m_bin_Sx2w[offset][2] += z2w;
0069
0070 bool inRange = true;
0071 if(ibin==0) inRange = false;
0072 else if(ibin==(parent::m_axes[0].m_number_of_bins+1)) inRange = false;
0073
0074 if(jbin==0) inRange = false;
0075 else if(jbin==(parent::m_axes[1].m_number_of_bins+1)) inRange = false;
0076
0077 if(kbin==0) inRange = false;
0078 else if(kbin==(parent::m_axes[2].m_number_of_bins+1)) inRange = false;
0079
0080 parent::m_all_entries++;
0081 if(inRange) {
0082 parent::m_in_range_plane_Sxyw[0] += aX * aY * aWeight;
0083 parent::m_in_range_plane_Sxyw[1] += aY * aZ * aWeight;
0084 parent::m_in_range_plane_Sxyw[2] += aZ * aX * aWeight;
0085
0086 // fast getters :
0087 parent::m_in_range_entries++;
0088 parent::m_in_range_Sw += aWeight;
0089 parent::m_in_range_Sw2 += aWeight*aWeight;
0090
0091 parent::m_in_range_Sxw[0] += xw;
0092 parent::m_in_range_Sx2w[0] += x2w;
0093
0094 parent::m_in_range_Sxw[1] += yw;
0095 parent::m_in_range_Sx2w[1] += y2w;
0096
0097 parent::m_in_range_Sxw[2] += zw;
0098 parent::m_in_range_Sx2w[2] += z2w;
0099 }
0100
0101 return true;
0102 }
0103
0104 bool set_bin_content(bn_t a_ibin,bn_t a_jbin,bn_t a_kbin,
0105 TN a_entries,TW a_Sw,TW a_Sw2,
0106 TC a_Sxw,TC a_Sx2w,TC a_Syw,TC a_Sy2w,TC a_Szw,TC a_Sz2w) {
0107 if(parent::m_dimension!=3) return false;
0108 if(a_ibin>(parent::m_axes[0].m_number_of_bins+1)) return false;
0109 if(a_jbin>(parent::m_axes[1].m_number_of_bins+1)) return false;
0110 if(a_kbin>(parent::m_axes[2].m_number_of_bins+1)) return false;
0111
0112 bool inRange = true;
0113 if(a_ibin==0) inRange = false;
0114 else if(a_ibin==(parent::m_axes[0].m_number_of_bins+1)) inRange = false;
0115 if(a_jbin==0) inRange = false;
0116 else if(a_jbin==(parent::m_axes[1].m_number_of_bins+1)) inRange = false;
0117 if(a_kbin==0) inRange = false;
0118 else if(a_kbin==(parent::m_axes[2].m_number_of_bins+1)) inRange = false;
0119
0120 TO offset = a_ibin + a_jbin * parent::m_axes[1].m_offset + a_kbin * parent::m_axes[2].m_offset;
0121
0122 parent::m_all_entries -= parent::m_bin_entries[offset];
0123 if(inRange) {
0124 parent::m_in_range_entries -= parent::m_bin_entries[offset];
0125 parent::m_in_range_Sw -= parent::m_bin_Sw[offset];
0126 parent::m_in_range_Sw2 -= parent::m_bin_Sw2[offset];
0127 parent::m_in_range_Sxw[0] -= parent::m_bin_Sxw[offset][0];
0128 parent::m_in_range_Sx2w[0] -= parent::m_bin_Sx2w[offset][0];
0129 parent::m_in_range_Sxw[1] -= parent::m_bin_Sxw[offset][1];
0130 parent::m_in_range_Sx2w[1] -= parent::m_bin_Sx2w[offset][1];
0131 parent::m_in_range_Sxw[2] -= parent::m_bin_Sxw[offset][2];
0132 parent::m_in_range_Sx2w[2] -= parent::m_bin_Sx2w[offset][2];
0133 }
0134
0135 parent::m_bin_entries[offset] = a_entries;
0136 parent::m_bin_Sw[offset] = a_Sw;
0137 parent::m_bin_Sw2[offset] = a_Sw2;
0138
0139 parent::m_bin_Sxw[offset][0] = a_Sxw;
0140 parent::m_bin_Sx2w[offset][0] = a_Sx2w;
0141 parent::m_bin_Sxw[offset][1] = a_Syw;
0142 parent::m_bin_Sx2w[offset][1] = a_Sy2w;
0143 parent::m_bin_Sxw[offset][2] = a_Szw;
0144 parent::m_bin_Sx2w[offset][2] = a_Sz2w;
0145
0146 parent::m_all_entries += a_entries;
0147 if(inRange) {
0148 //parent::m_in_range_plane_Sxyw[0,1,2] ??? ill-defined.
0149
0150 parent::m_in_range_entries += a_entries;
0151 parent::m_in_range_Sw += a_Sw;
0152 parent::m_in_range_Sw2 += a_Sw2;
0153
0154 parent::m_in_range_Sxw[0] += a_Sxw;
0155 parent::m_in_range_Sx2w[0] += a_Sx2w;
0156 parent::m_in_range_Sxw[1] += a_Syw;
0157 parent::m_in_range_Sx2w[1] += a_Sy2w;
0158 parent::m_in_range_Sxw[2] += a_Szw;
0159 parent::m_in_range_Sx2w[2] += a_Sz2w;
0160 }
0161
0162 return true;
0163 }
0164
0165 bool get_bin_content(bn_t a_ibin,bn_t a_jbin,bn_t a_kbin,
0166 TN& a_entries,TW& a_Sw,TW& a_Sw2,
0167 TC& a_Sxw,TC& a_Sx2w,
0168 TC& a_Syw,TC& a_Sy2w,
0169 TC& a_Szw,TC& a_Sz2w) {
0170 if(parent::m_dimension!=3) {
0171 a_entries = 0;a_Sw = 0;a_Sw2 = 0;
0172 a_Sxw = 0;a_Sx2w = 0;
0173 a_Syw = 0;a_Sy2w = 0;
0174 a_Szw = 0;a_Sz2w = 0;
0175 return false;
0176 }
0177 if(a_ibin>(parent::m_axes[0].m_number_of_bins+1)) {
0178 a_entries = 0;a_Sw = 0;a_Sw2 = 0;
0179 a_Sxw = 0;a_Sx2w = 0;
0180 a_Syw = 0;a_Sy2w = 0;
0181 a_Szw = 0;a_Sz2w = 0;
0182 return false;
0183 }
0184 if(a_jbin>(parent::m_axes[1].m_number_of_bins+1)) {
0185 a_entries = 0;a_Sw = 0;a_Sw2 = 0;
0186 a_Sxw = 0;a_Sx2w = 0;
0187 a_Syw = 0;a_Sy2w = 0;
0188 a_Szw = 0;a_Sz2w = 0;
0189 return false;
0190 }
0191 if(a_kbin>(parent::m_axes[2].m_number_of_bins+1)) {
0192 a_entries = 0;a_Sw = 0;a_Sw2 = 0;
0193 a_Sxw = 0;a_Sx2w = 0;
0194 a_Syw = 0;a_Sy2w = 0;
0195 a_Szw = 0;a_Sz2w = 0;
0196 return false;
0197 }
0198
0199 TO offset = a_ibin + a_jbin * parent::m_axes[1].m_offset + a_kbin * parent::m_axes[2].m_offset;
0200
0201 a_entries = parent::m_bin_entries[offset];
0202 a_Sw = parent::m_bin_Sw[offset];
0203 a_Sw2 = parent::m_bin_Sw2[offset];
0204
0205 a_Sxw = parent::m_bin_Sxw[offset][0];
0206 a_Sx2w = parent::m_bin_Sx2w[offset][0];
0207 a_Syw = parent::m_bin_Sxw[offset][1];
0208 a_Sy2w = parent::m_bin_Sx2w[offset][1];
0209 a_Szw = parent::m_bin_Sxw[offset][2];
0210 a_Sz2w = parent::m_bin_Sx2w[offset][2];
0211
0212 return true;
0213 }
0214
0215 bool add(const h3& a_histo){
0216 parent::base_add(a_histo);
0217 return true;
0218 }
0219 bool subtract(const h3& a_histo){
0220 parent::base_subtract(a_histo);
0221 return true;
0222 }
0223
0224 bool multiply(const h3& a_histo) {
0225 return parent::base_multiply(a_histo);
0226 }
0227
0228 bool divide(const h3& a_histo) {
0229 return parent::base_divide(a_histo);
0230 }
0231
0232 bool equals_TH(const h3& a_from,const TW& a_prec,TW(*a_fabs)(TW)) const {
0233 if(!parent::equals_TH(a_from,a_prec,a_fabs,true)) return false;
0234 return true;
0235 }
0236
0237 void not_a_profile() const {}
0238
0239 public: //CERN-ROOT API
0240 bool Fill(TC aX,TC aY,TC aZ,TW aWeight = 1) {return fill(aX,aY,aZ,aWeight);}
0241
0242 public:
0243 /*
0244 // Slices :
0245 h2d* slice_xy(int aKbeg,int aKend) const;
0246 h2d* slice_yz(int aIbeg,int aIend) const;
0247 h2d* slice_xz(int aJbeg,int aJend) const;
0248 */
0249 public:
0250 h3(const std::string& a_title,
0251 bn_t aXnumber,TC aXmin,TC aXmax,
0252 bn_t aYnumber,TC aYmin,TC aYmax,
0253 bn_t aZnumber,TC aZmin,TC aZmax)
0254 :parent(a_title,aXnumber,aXmin,aXmax,
0255 aYnumber,aYmin,aYmax,
0256 aZnumber,aZmin,aZmax)
0257 {}
0258
0259 h3(const std::string& a_title,
0260 const std::vector<TC>& a_edges_x,
0261 const std::vector<TC>& a_edges_y,
0262 const std::vector<TC>& a_edges_z)
0263 :parent(a_title,a_edges_x,a_edges_y,a_edges_z)
0264 {}
0265
0266 virtual ~h3(){}
0267 public:
0268 h3(const h3& a_from): parent(a_from){}
0269 h3& operator=(const h3& a_from){
0270 parent::operator=(a_from);
0271 return *this;
0272 }
0273 };
0274
0275 }}
0276
0277 #endif
0278
0279
0280
0281
