Warning, /include/Geant4/tools/histo/p2 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_p2
0005 #define tools_histo_p2
0006
0007 #include "b2"
0008 #include "profile_data"
0009
0010 namespace tools {
0011 namespace histo {
0012
0013 template <class TC,class TO,class TN,class TW,class TH,class TV>
0014 class p2 : public b2<TC,TO,TN,TW,TH> {
0015 typedef b2<TC,TO,TN,TW,TH> parent;
0016 public:
0017 typedef profile_data<TC,TO,TN,TW,TV> pd_t;
0018 typedef typename parent::bn_t bn_t;
0019 typedef std::vector<TV> vs_t;
0020 protected:
0021 virtual TH get_bin_height(TO a_offset) const {
0022 return (parent::m_bin_Sw[a_offset] ? (m_bin_Svw[a_offset]/parent::m_bin_Sw[a_offset]):0);
0023 }
0024 public:
0025 bool equals(const p2& a_from,const TW& a_prec,TW(*a_fabs)(TW)) const {
0026 if(!parent::equals(a_from,a_prec,a_fabs)) return false;
0027 if(m_cut_v!=a_from.m_cut_v) return false;
0028 if(!numbers_are_equal(m_min_v,a_from.m_min_v,a_prec,a_fabs)) return false;
0029 if(!numbers_are_equal(m_max_v,a_from.m_max_v,a_prec,a_fabs)) return false;
0030 if(!vectors_are_equal(m_bin_Svw,a_from.m_bin_Svw,a_prec,a_fabs)) return false;
0031 if(!vectors_are_equal(m_bin_Sv2w,a_from.m_bin_Sv2w,a_prec,a_fabs)) return false;
0032 return true;
0033 }
0034 bool equals_TH(const p2& a_from,const TW& a_prec,TW(*a_fabs)(TW)) const {
0035 if(!parent::equals_TH(a_from,a_prec,a_fabs,false)) return false;
0036 if(m_cut_v!=a_from.m_cut_v) return false;
0037 if(!numbers_are_equal(m_min_v,a_from.m_min_v,a_prec,a_fabs)) return false;
0038 if(!numbers_are_equal(m_max_v,a_from.m_max_v,a_prec,a_fabs)) return false;
0039 if(!vectors_are_equal(m_bin_Svw,a_from.m_bin_Svw,a_prec,a_fabs)) return false;
0040 if(!vectors_are_equal(m_bin_Sv2w,a_from.m_bin_Sv2w,a_prec,a_fabs)) return false;
0041 return true;
0042 }
0043
0044 virtual TH bin_error(int aI,int aJ) const { //TH should be the same as TV
0045 TO offset;
0046 if(!parent::_find_offset(aI,aJ,offset)) return 0;
0047
0048 //FIXME Is it correct ?
0049 // TProfile::GetBinError with kERRORMEAN mode does :
0050 // Stat_t cont = fArray[bin]; //Svw (see TProfile::Fill)
0051 // Stat_t sum = parent::m_bin_entries.fArray[bin]; //Sw
0052 // Stat_t err2 = fSumw2.fArray[bin]; //Sv2w
0053 // if (sum == 0) return 0;
0054 // Stat_t eprim;
0055 // Stat_t contsum = cont/sum;
0056 // Stat_t eprim2 = TMath::Abs(err2/sum - contsum*contsum);
0057 // eprim = TMath::Sqrt(eprim2);
0058 // ... ???
0059 // if (fErrorMode == kERRORMEAN) return eprim/TMath::Sqrt(sum);
0060
0061 TW sw = parent::m_bin_Sw[offset]; //ROOT sum
0062 if(sw==0) return 0;
0063 TV svw = m_bin_Svw[offset]; //ROOT cont
0064 TV sv2w = m_bin_Sv2w[offset]; //ROOT err2
0065 TV mean = (svw / sw); //ROOT contsum
0066 TV rms = ::sqrt(::fabs((sv2w/sw) - mean * mean)); //ROOT eprim
0067 return rms/::sqrt(sw); //ROOT kERRORMEAN mode returned value
0068 }
0069
0070 public:
0071 bool multiply(TW a_factor){
0072 if(!parent::base_multiply(a_factor)) return false;
0073 for(bn_t ibin=0;ibin<parent::m_bin_number;ibin++) {
0074 m_bin_Svw[ibin] *= a_factor;
0075 }
0076 return true;
0077 }
0078 bool scale(TW a_factor) {return multiply(a_factor);}
0079
0080 TV bin_Svw(int aI,int aJ) const {
0081 TO offset;
0082 if(!parent::_find_offset(aI,aJ,offset)) return 0;
0083 return m_bin_Svw[offset];
0084 }
0085 TV bin_Sv2w(int aI,int aJ) const {
0086 TO offset;
0087 if(!parent::_find_offset(aI,aJ,offset)) return 0;
0088 return m_bin_Sv2w[offset];
0089 }
0090
0091 bool reset() {
0092 parent::base_reset();
0093 for(bn_t ibin=0;ibin<parent::m_bin_number;ibin++) {
0094 m_bin_Svw[ibin] = 0;
0095 m_bin_Sv2w[ibin] = 0;
0096 }
0097 return true;
0098 }
0099
0100 void copy_from_data(const pd_t& a_from) {
0101 parent::base_from_data(a_from);
0102 m_bin_Svw = a_from.m_bin_Svw;
0103 m_bin_Sv2w = a_from.m_bin_Sv2w;
0104 m_cut_v = a_from.m_cut_v;
0105 m_min_v = a_from.m_min_v;
0106 m_max_v = a_from.m_max_v;
0107 }
0108 pd_t get_histo_data() const {
0109 pd_t hd(parent::dac());
0110 hd.m_is_profile = true;
0111 hd.m_bin_Svw = m_bin_Svw;
0112 hd.m_bin_Sv2w = m_bin_Sv2w;
0113 hd.m_cut_v = m_cut_v;
0114 hd.m_min_v = m_min_v;
0115 hd.m_max_v = m_max_v;
0116 return hd;
0117 }
0118
0119 bool fill(TC aX,TC aY,TV aV,TW aWeight = 1) {
0120 if(parent::m_dimension!=2) return false;
0121
0122 if(m_cut_v) {
0123 if( (aV<m_min_v) || (aV>=m_max_v) ) {
0124 return true;
0125 }
0126 }
0127
0128 bn_t ibin,jbin;
0129 if(!parent::m_axes[0].coord_to_absolute_index(aX,ibin)) return false;
0130 if(!parent::m_axes[1].coord_to_absolute_index(aY,jbin)) return false;
0131 bn_t offset = ibin + jbin * parent::m_axes[1].m_offset;
0132
0133 parent::m_bin_entries[offset]++;
0134 parent::m_bin_Sw[offset] += aWeight;
0135 parent::m_bin_Sw2[offset] += aWeight * aWeight;
0136
0137 TC xw = aX * aWeight;
0138 TC x2w = aX * xw;
0139 parent::m_bin_Sxw[offset][0] += xw;
0140 parent::m_bin_Sx2w[offset][0] += x2w;
0141
0142 TC yw = aY * aWeight;
0143 TC y2w = aY * yw;
0144 parent::m_bin_Sxw[offset][1] += yw;
0145 parent::m_bin_Sx2w[offset][1] += y2w;
0146
0147 bool inRange = true;
0148 if(ibin==0) inRange = false;
0149 else if(ibin==(parent::m_axes[0].m_number_of_bins+1)) inRange = false;
0150
0151 if(jbin==0) inRange = false;
0152 else if(jbin==(parent::m_axes[1].m_number_of_bins+1)) inRange = false;
0153
0154 parent::m_all_entries++;
0155 if(inRange) {
0156 parent::m_in_range_plane_Sxyw[0] += aX * aY * aWeight;
0157
0158 // fast getters :
0159 parent::m_in_range_entries++;
0160 parent::m_in_range_Sw += aWeight;
0161 parent::m_in_range_Sw2 += aWeight*aWeight;
0162
0163 parent::m_in_range_Sxw[0] += xw;
0164 parent::m_in_range_Sx2w[0] += x2w;
0165
0166 parent::m_in_range_Sxw[1] += yw;
0167 parent::m_in_range_Sx2w[1] += y2w;
0168 }
0169
0170 // Profile part :
0171 TV vw = aV * aWeight;
0172 m_bin_Svw[offset] += vw;
0173 m_bin_Sv2w[offset] += aV * vw;
0174
0175 return true;
0176 }
0177
0178 TV bin_rms_value(int aI,int aJ) const {
0179 TO offset;
0180 if(!parent::_find_offset(aI,aJ,offset)) return 0;
0181
0182 TW sw = parent::m_bin_Sw[offset];
0183 if(sw==0) return 0;
0184 TV svw = m_bin_Svw[offset];
0185 TV sv2w = m_bin_Sv2w[offset];
0186 TV mean = (svw / sw);
0187 return ::sqrt(::fabs((sv2w / sw) - mean * mean));
0188 }
0189
0190 bool add(const p2& a_histo){
0191 parent::base_add(a_histo);
0192 for(bn_t ibin=0;ibin<parent::m_bin_number;ibin++) {
0193 m_bin_Svw[ibin] += a_histo.m_bin_Svw[ibin];
0194 m_bin_Sv2w[ibin] += a_histo.m_bin_Sv2w[ibin];
0195 }
0196 return true;
0197 }
0198 bool subtract(const p2& a_histo){
0199 parent::base_subtract(a_histo);
0200 for(bn_t ibin=0;ibin<parent::m_bin_number;ibin++) {
0201 m_bin_Svw[ibin] -= a_histo.m_bin_Svw[ibin];
0202 m_bin_Sv2w[ibin] -= a_histo.m_bin_Sv2w[ibin];
0203 }
0204 return true;
0205 }
0206
0207 bool cut_v() const {return m_cut_v;}
0208 TV min_v() const {return m_min_v;}
0209 TV max_v() const {return m_max_v;}
0210 public:
0211 p2(const std::string& a_title,
0212 bn_t aXnumber,TC aXmin,TC aXmax,
0213 bn_t aYnumber,TC aYmin,TC aYmax)
0214 :parent(a_title,aXnumber,aXmin,aXmax,aYnumber,aYmin,aYmax)
0215 ,m_cut_v(false)
0216 ,m_min_v(0)
0217 ,m_max_v(0)
0218 {
0219 m_bin_Svw.resize(parent::m_bin_number,0);
0220 m_bin_Sv2w.resize(parent::m_bin_number,0);
0221 }
0222
0223 p2(const std::string& a_title,
0224 bn_t aXnumber,TC aXmin,TC aXmax,
0225 bn_t aYnumber,TC aYmin,TC aYmax,
0226 TV aVmin,TV aVmax)
0227 :parent(a_title,aXnumber,aXmin,aXmax,aYnumber,aYmin,aYmax)
0228 ,m_cut_v(true)
0229 ,m_min_v(aVmin)
0230 ,m_max_v(aVmax)
0231 {
0232 m_bin_Svw.resize(parent::m_bin_number,0);
0233 m_bin_Sv2w.resize(parent::m_bin_number,0);
0234 }
0235
0236 p2(const std::string& a_title,
0237 const std::vector<TC>& a_edges_x,
0238 const std::vector<TC>& a_edges_y)
0239 :parent(a_title,a_edges_x,a_edges_y)
0240 ,m_cut_v(false)
0241 ,m_min_v(0)
0242 ,m_max_v(0)
0243 {
0244 m_bin_Svw.resize(parent::m_bin_number,0);
0245 m_bin_Sv2w.resize(parent::m_bin_number,0);
0246 }
0247
0248 p2(const std::string& a_title,
0249 const std::vector<TC>& a_edges_x,
0250 const std::vector<TC>& a_edges_y,
0251 TV aVmin,TV aVmax)
0252 :parent(a_title,a_edges_x,a_edges_y)
0253 ,m_cut_v(true)
0254 ,m_min_v(aVmin)
0255 ,m_max_v(aVmax)
0256 {
0257 m_bin_Svw.resize(parent::m_bin_number,0);
0258 m_bin_Sv2w.resize(parent::m_bin_number,0);
0259 }
0260
0261 virtual ~p2(){}
0262 public:
0263 p2(const p2& a_from)
0264 :parent(a_from)
0265 ,m_cut_v(a_from.m_cut_v)
0266 ,m_min_v(a_from.m_min_v)
0267 ,m_max_v(a_from.m_max_v)
0268 ,m_bin_Svw(a_from.m_bin_Svw)
0269 ,m_bin_Sv2w(a_from.m_bin_Sv2w)
0270 {}
0271 p2& operator=(const p2& a_from){
0272 parent::operator=(a_from);
0273 m_cut_v = a_from.m_cut_v;
0274 m_min_v = a_from.m_min_v;
0275 m_max_v = a_from.m_max_v;
0276 m_bin_Svw = a_from.m_bin_Svw;
0277 m_bin_Sv2w = a_from.m_bin_Sv2w;
0278 return *this;
0279 }
0280 public:
0281 bool configure(bn_t aXnumber,TC aXmin,TC aXmax,bn_t aYnumber,TC aYmin,TC aYmax){
0282 if(!parent::configure(aXnumber,aXmin,aXmax,aYnumber,aYmin,aYmax)) return false;
0283 m_bin_Svw.clear();
0284 m_bin_Sv2w.clear();
0285 m_bin_Svw.resize(parent::m_bin_number,0);
0286 m_bin_Sv2w.resize(parent::m_bin_number,0);
0287 m_cut_v = false;
0288 m_min_v = 0;
0289 m_max_v = 0;
0290 return true;
0291 }
0292 bool configure(const std::vector<TC>& a_edges_x,const std::vector<TC>& a_edges_y) {
0293 if(!parent::configure(a_edges_x,a_edges_y)) return false;
0294 m_bin_Svw.clear();
0295 m_bin_Sv2w.clear();
0296 m_bin_Svw.resize(parent::m_bin_number,0);
0297 m_bin_Sv2w.resize(parent::m_bin_number,0);
0298 m_cut_v = false;
0299 m_min_v = 0;
0300 m_max_v = 0;
0301 return true;
0302 }
0303 bool configure(bn_t aXnumber,TC aXmin,TC aXmax,bn_t aYnumber,TC aYmin,TC aYmax,TV aVmin,TV aVmax){
0304 if(!parent::configure(aXnumber,aXmin,aXmax,aYnumber,aYmin,aYmax)) return false;
0305 m_bin_Svw.clear();
0306 m_bin_Sv2w.clear();
0307 m_bin_Svw.resize(parent::m_bin_number,0);
0308 m_bin_Sv2w.resize(parent::m_bin_number,0);
0309 m_cut_v = true;
0310 m_min_v = aVmin;
0311 m_max_v = aVmax;
0312 return true;
0313 }
0314 bool configure(const std::vector<TC>& a_edges_x,const std::vector<TC>& a_edges_y,TV aVmin,TV aVmax) {
0315 if(!parent::configure(a_edges_x,a_edges_y)) return false;
0316 m_bin_Svw.clear();
0317 m_bin_Sv2w.clear();
0318 m_bin_Svw.resize(parent::m_bin_number,0);
0319 m_bin_Sv2w.resize(parent::m_bin_number,0);
0320 m_cut_v = true;
0321 m_min_v = aVmin;
0322 m_max_v = aVmax;
0323 return true;
0324 }
0325 public:
0326 const vs_t& bins_sum_vw() const {return m_bin_Svw;}
0327 const vs_t& bins_sum_v2w() const {return m_bin_Sv2w;}
0328
0329 TW get_Svw() const {
0330 TW sw = 0;
0331 for(TO ibin=0;ibin<parent::m_bin_number;ibin++) {
0332 if(!histo::is_out(parent::m_axes,ibin)) {
0333 sw += m_bin_Svw[ibin];
0334 }
0335 }
0336 return sw;
0337 }
0338 TW get_Sv2w() const {
0339 TW sw = 0;
0340 for(TO ibin=0;ibin<parent::m_bin_number;ibin++) {
0341 if(!histo::is_out(parent::m_axes,ibin)) {
0342 sw += m_bin_Sv2w[ibin];
0343 }
0344 }
0345 return sw;
0346 }
0347 protected:
0348 bool m_cut_v;
0349 TV m_min_v;
0350 TV m_max_v;
0351 vs_t m_bin_Svw;
0352 vs_t m_bin_Sv2w;
0353 };
0354
0355 }}
0356
0357 #endif
0358
