Warning, /include/Geant4/tools/histo/p1 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_p1
0005 #define tools_histo_p1
0006
0007 #include "b1"
0008 #include "profile_data"
0009
0010 namespace tools {
0011 namespace histo {
0012
0013 //TC is for a coordinate.
0014 //TO is for an offset/index to identify/find a bind.
0015 //TN is for a number of entries.
0016 //TW is for a weight.
0017 //TH is for a height. Should be the same as TV.
0018 //TV is for a value (in general same as TC).
0019
0020 template <class TC,class TO,class TN,class TW,class TH,class TV>
0021 class p1 : public b1<TC,TO,TN,TW,TH> {
0022 typedef b1<TC,TO,TN,TW,TH> parent;
0023 public:
0024 typedef profile_data<TC,TO,TN,TW,TV> pd_t;
0025 typedef typename parent::bn_t bn_t;
0026 typedef typename parent::axis_t axis_t;
0027 typedef std::vector<TV> vs_t;
0028 protected:
0029 virtual TH get_bin_height(TO a_offset) const {
0030 return (parent::m_bin_Sw[a_offset] ? (m_bin_Svw[a_offset]/parent::m_bin_Sw[a_offset]):0);
0031 }
0032 public:
0033 bool equals(const p1& a_from,const TW& a_prec,TW(*a_fabs)(TW)) const {
0034 if(!parent::equals(a_from,a_prec,a_fabs)) return false;
0035 if(m_cut_v!=a_from.m_cut_v) return false;
0036 if(!numbers_are_equal(m_min_v,a_from.m_min_v,a_prec,a_fabs)) return false;
0037 if(!numbers_are_equal(m_max_v,a_from.m_max_v,a_prec,a_fabs)) return false;
0038 if(!vectors_are_equal(m_bin_Svw,a_from.m_bin_Svw,a_prec,a_fabs)) return false;
0039 if(!vectors_are_equal(m_bin_Sv2w,a_from.m_bin_Sv2w,a_prec,a_fabs)) return false;
0040 return true;
0041 }
0042 bool equals_TH(const p1& a_from,const TW& a_prec,TW(*a_fabs)(TW)) const {
0043 if(!parent::equals_TH(a_from,a_prec,a_fabs,false)) return false;
0044 if(m_cut_v!=a_from.m_cut_v) return false;
0045 if(!numbers_are_equal(m_min_v,a_from.m_min_v,a_prec,a_fabs)) return false;
0046 if(!numbers_are_equal(m_max_v,a_from.m_max_v,a_prec,a_fabs)) return false;
0047 if(!vectors_are_equal(m_bin_Svw,a_from.m_bin_Svw,a_prec,a_fabs)) return false;
0048 if(!vectors_are_equal(m_bin_Sv2w,a_from.m_bin_Sv2w,a_prec,a_fabs)) return false;
0049 return true;
0050 }
0051
0052 virtual TH bin_error(int aI) const { //TH should be the same as TV
0053 TO offset;
0054 if(!parent::_find_offset(aI,offset)) return 0;
0055
0056 //FIXME Is it correct ?
0057 // TProfile::GetBinError with kERRORMEAN mode does :
0058 // Stat_t cont = fArray[bin]; //Svw (see TProfile::Fill)
0059 // Stat_t sum = parent::m_bin_entries.fArray[bin]; //Sw
0060 // Stat_t err2 = fSumw2.fArray[bin]; //Sv2w
0061 // if (sum == 0) return 0;
0062 // Stat_t eprim;
0063 // Stat_t contsum = cont/sum;
0064 // Stat_t eprim2 = TMath::Abs(err2/sum - contsum*contsum);
0065 // eprim = TMath::Sqrt(eprim2);
0066 // ... ???
0067 // if (fErrorMode == kERRORMEAN) return eprim/TMath::Sqrt(sum);
0068
0069 TW sw = parent::m_bin_Sw[offset]; //ROOT sum
0070 if(sw==0) return 0;
0071 TV svw = m_bin_Svw[offset]; //ROOT cont
0072 TV sv2w = m_bin_Sv2w[offset]; //ROOT err2
0073 TV _mean = (svw / sw); //ROOT contsum
0074 TV _rms = ::sqrt(::fabs((sv2w/sw) - _mean * _mean)); //ROOT eprim
0075 // rms = get_bin_rms_value.
0076 return _rms/::sqrt(sw); //ROOT kERRORMEAN mode returned value
0077 }
0078
0079 public:
0080 bool multiply(TW a_factor){
0081 if(!parent::base_multiply(a_factor)) return false;
0082 for(bn_t ibin=0;ibin<parent::m_bin_number;ibin++) {
0083 m_bin_Svw[ibin] *= a_factor;
0084 }
0085 return true;
0086 }
0087 bool scale(TW a_factor) {return multiply(a_factor);}
0088
0089 TV bin_Svw(int aI) const {
0090 TO offset;
0091 if(!parent::_find_offset(aI,offset)) return 0;
0092 return m_bin_Svw[offset];
0093 }
0094 TV bin_Sv2w(int aI) const {
0095 TO offset;
0096 if(!parent::_find_offset(aI,offset)) return 0;
0097 return m_bin_Sv2w[offset];
0098 }
0099
0100 bool reset() {
0101 parent::base_reset();
0102 for(bn_t ibin=0;ibin<parent::m_bin_number;ibin++) {
0103 m_bin_Svw[ibin] = 0;
0104 m_bin_Sv2w[ibin] = 0;
0105 }
0106 return true;
0107 }
0108
0109 void copy_from_data(const pd_t& a_from) {
0110 parent::base_from_data(a_from);
0111 m_bin_Svw = a_from.m_bin_Svw;
0112 m_bin_Sv2w = a_from.m_bin_Sv2w;
0113 m_cut_v = a_from.m_cut_v;
0114 m_min_v = a_from.m_min_v;
0115 m_max_v = a_from.m_max_v;
0116 }
0117 pd_t get_histo_data() const {
0118 pd_t hd(parent::dac());
0119 hd.m_is_profile = true;
0120 hd.m_bin_Svw = m_bin_Svw;
0121 hd.m_bin_Sv2w = m_bin_Sv2w;
0122 hd.m_cut_v = m_cut_v;
0123 hd.m_min_v = m_min_v;
0124 hd.m_max_v = m_max_v;
0125 return hd;
0126 }
0127
0128 bool fill(TC aX,TV aV,TW aWeight = 1) {
0129 if(parent::m_dimension!=1) return false;
0130
0131 if(m_cut_v) {
0132 if( (aV<m_min_v) || (aV>=m_max_v) ) {
0133 return true;
0134 }
0135 }
0136
0137 bn_t ibin;
0138 if(!parent::m_axes[0].coord_to_absolute_index(aX,ibin)) return false;
0139 TO offset = ibin;
0140
0141 parent::m_bin_entries[offset]++;
0142 parent::m_bin_Sw[offset] += aWeight;
0143 parent::m_bin_Sw2[offset] += aWeight * aWeight;
0144
0145 TC xw = aX * aWeight;
0146 TC x2w = aX * xw;
0147 parent::m_bin_Sxw[offset][0] += xw;
0148 parent::m_bin_Sx2w[offset][0] += x2w;
0149
0150 bool inRange = true;
0151 if(ibin==0) inRange = false;
0152 else if(ibin==(parent::m_axes[0].m_number_of_bins+1)) inRange = false;
0153
0154 parent::m_all_entries++;
0155 if(inRange) {
0156 // fast getters :
0157 parent::m_in_range_entries++;
0158 parent::m_in_range_Sw += aWeight;
0159 parent::m_in_range_Sw2 += aWeight*aWeight;
0160
0161 parent::m_in_range_Sxw[0] += xw;
0162 parent::m_in_range_Sx2w[0] += x2w;
0163 }
0164
0165 // Profile part :
0166 TV vw = aV * aWeight;
0167 m_bin_Svw[offset] += vw;
0168 m_bin_Sv2w[offset] += aV * vw;
0169
0170 return true;
0171 }
0172
0173 bool set_bin_content(bn_t a_ibin,TN a_entries,TW a_Sw,TW a_Sw2,TC a_Sxw,TC a_Sx2w,TC a_Svw,TC a_Sv2w) {
0174 if(parent::m_dimension!=1) return false;
0175 if(a_ibin>(parent::m_axes[0].m_number_of_bins+1)) return false;
0176
0177 bool inRange = true;
0178 if(a_ibin==0) inRange = false;
0179 else if(a_ibin==(parent::m_axes[0].m_number_of_bins+1)) inRange = false;
0180
0181 TO offset = a_ibin;
0182
0183 parent::m_all_entries -= parent::m_bin_entries[offset];
0184 if(inRange) {
0185 parent::m_in_range_entries -= parent::m_bin_entries[offset];
0186 parent::m_in_range_Sw -= parent::m_bin_Sw[offset];
0187 parent::m_in_range_Sw2 -= parent::m_bin_Sw2[offset];
0188 parent::m_in_range_Sxw[0] -= parent::m_bin_Sxw[offset][0];
0189 parent::m_in_range_Sx2w[0] -= parent::m_bin_Sx2w[offset][0];
0190 }
0191
0192 parent::m_bin_entries[offset] = a_entries;
0193 parent::m_bin_Sw[offset] = a_Sw;
0194 parent::m_bin_Sw2[offset] = a_Sw2;
0195
0196 parent::m_bin_Sxw[offset][0] = a_Sxw;
0197 parent::m_bin_Sx2w[offset][0] = a_Sx2w;
0198
0199 parent::m_all_entries += a_entries;
0200 if(inRange) {
0201 parent::m_in_range_entries += a_entries;
0202 parent::m_in_range_Sw += a_Sw;
0203 parent::m_in_range_Sw2 += a_Sw2;
0204
0205 parent::m_in_range_Sxw[0] += a_Sxw;
0206 parent::m_in_range_Sx2w[0] += a_Sx2w;
0207 }
0208
0209 // Profile part :
0210 m_bin_Svw[offset] = a_Svw;
0211 m_bin_Sv2w[offset] = a_Sv2w;
0212
0213 return true;
0214 }
0215
0216 bool get_bin_content(bn_t a_ibin,TN& a_entries,TW& a_Sw,TW& a_Sw2,TC& a_Sxw,TC& a_Sx2w,TC& a_Svw,TC& a_Sv2w) {
0217 if(parent::m_dimension!=1) {
0218 a_entries = 0;a_Sw = 0;a_Sw2 = 0;a_Sxw = 0;a_Sx2w = 0;
0219 return false;
0220 }
0221 if(a_ibin>(parent::m_axes[0].m_number_of_bins+1)) {
0222 a_entries = 0;a_Sw = 0;a_Sw2 = 0;a_Sxw = 0;a_Sx2w = 0;
0223 return false;
0224 }
0225
0226 TO offset = a_ibin;
0227
0228 a_entries = parent::m_bin_entries[offset];
0229 a_Sw = parent::m_bin_Sw[offset];
0230 a_Sw2 = parent::m_bin_Sw2[offset];
0231
0232 a_Sxw = parent::m_bin_Sxw[offset][0];
0233 a_Sx2w = parent::m_bin_Sx2w[offset][0];
0234
0235 // Profile part :
0236 a_Svw = m_bin_Svw[offset];
0237 a_Sv2w = m_bin_Sv2w[offset];
0238
0239 return true;
0240 }
0241
0242 TV bin_rms_value(int aI) const {
0243 TO offset;
0244 if(!parent::_find_offset(aI,offset)) return 0;
0245 TW sw = parent::m_bin_Sw[offset];
0246 if(sw==0) return 0;
0247 TV svw = m_bin_Svw[offset];
0248 TV sv2w = m_bin_Sv2w[offset];
0249 TV _mean = (svw / sw);
0250 return ::sqrt(::fabs((sv2w / sw) - _mean * _mean));
0251 }
0252
0253 bool add(const p1& a_histo){
0254 parent::base_add(a_histo);
0255 for(bn_t ibin=0;ibin<parent::m_bin_number;ibin++) {
0256 m_bin_Svw[ibin] += a_histo.m_bin_Svw[ibin];
0257 m_bin_Sv2w[ibin] += a_histo.m_bin_Sv2w[ibin];
0258 }
0259 return true;
0260 }
0261 bool subtract(const p1& a_histo){
0262 parent::base_subtract(a_histo);
0263 for(bn_t ibin=0;ibin<parent::m_bin_number;ibin++) {
0264 m_bin_Svw[ibin] -= a_histo.m_bin_Svw[ibin];
0265 m_bin_Sv2w[ibin] -= a_histo.m_bin_Sv2w[ibin];
0266 }
0267 return true;
0268 }
0269
0270 bool gather_bins(unsigned int a_factor) { //for exa 2,3.
0271 if(!a_factor) return false;
0272
0273 // actual bin number must be a multiple of a_factor.
0274
0275 const axis_t& _axis = parent::axis();
0276
0277 bn_t n = _axis.bins();
0278 if(!n) return false;
0279
0280 bn_t new_n = n/a_factor;
0281 if(a_factor*new_n!=n) return false;
0282
0283 p1* new_h = 0;
0284 if(_axis.is_fixed_binning()) {
0285 new_h = new p1(parent::m_title,new_n,_axis.lower_edge(),_axis.upper_edge());
0286 } else {
0287 const std::vector<TC>& _edges = _axis.edges();
0288 std::vector<TC> new_edges(new_n+1);
0289 for(bn_t ibin=0;ibin<new_n;ibin++) {
0290 new_edges[ibin] = _edges[ibin*a_factor];
0291 }
0292 new_edges[new_n] = _edges[n]; //upper edge.
0293 new_h = new p1(parent::m_title,new_edges);
0294 }
0295 if(!new_h) return false;
0296
0297 new_h->m_cut_v = m_cut_v;
0298 new_h->m_min_v = m_min_v;
0299 new_h->m_max_v = m_max_v;
0300
0301 bn_t offset,new_offset,offac;
0302 for(bn_t ibin=0;ibin<new_n;ibin++) {
0303 new_offset = ibin+1;
0304 offset = a_factor*ibin+1;
0305 for(unsigned int ifac=0;ifac<a_factor;ifac++) {
0306 offac = offset+ifac;
0307 new_h->m_bin_entries[new_offset] += parent::m_bin_entries[offac];
0308 new_h->m_bin_Sw[new_offset] += parent::m_bin_Sw[offac];
0309 new_h->m_bin_Sw2[new_offset] += parent::m_bin_Sw2[offac];
0310 new_h->m_bin_Sxw[new_offset][0] += parent::m_bin_Sxw[offac][0];
0311 new_h->m_bin_Sx2w[new_offset][0] += parent::m_bin_Sx2w[offac][0];
0312
0313 new_h->m_bin_Svw[new_offset] += m_bin_Svw[offac];
0314 new_h->m_bin_Sv2w[new_offset] += m_bin_Sv2w[offac];
0315 }
0316 }
0317
0318 //underflow :
0319 new_offset = 0;
0320 offac = 0;
0321 new_h->m_bin_entries[new_offset] = parent::m_bin_entries[offac];
0322 new_h->m_bin_Sw[new_offset] = parent::m_bin_Sw[offac];
0323 new_h->m_bin_Sw2[new_offset] = parent::m_bin_Sw2[offac];
0324 new_h->m_bin_Sxw[new_offset][0] = parent::m_bin_Sxw[offac][0];
0325 new_h->m_bin_Sx2w[new_offset][0] = parent::m_bin_Sx2w[offac][0];
0326 new_h->m_bin_Svw[new_offset] = m_bin_Svw[offac];
0327 new_h->m_bin_Sv2w[new_offset] = m_bin_Sv2w[offac];
0328
0329 //overflow :
0330 new_offset = new_n+1;
0331 offac = n+1;
0332 new_h->m_bin_entries[new_offset] = parent::m_bin_entries[offac];
0333 new_h->m_bin_Sw[new_offset] = parent::m_bin_Sw[offac];
0334 new_h->m_bin_Sw2[new_offset] = parent::m_bin_Sw2[offac];
0335 new_h->m_bin_Sxw[new_offset][0] = parent::m_bin_Sxw[offac][0];
0336 new_h->m_bin_Sx2w[new_offset][0] = parent::m_bin_Sx2w[offac][0];
0337 new_h->m_bin_Svw[new_offset] = m_bin_Svw[offac];
0338 new_h->m_bin_Sv2w[new_offset] = m_bin_Sv2w[offac];
0339
0340 *this = *new_h;
0341 return true;
0342 }
0343
0344 bool cut_v() const {return m_cut_v;}
0345 TV min_v() const {return m_min_v;}
0346 TV max_v() const {return m_max_v;}
0347
0348 public:
0349 p1(const std::string& a_title,bn_t aXnumber,TC aXmin,TC aXmax)
0350 :parent(a_title,aXnumber,aXmin,aXmax)
0351 ,m_cut_v(false)
0352 ,m_min_v(0)
0353 ,m_max_v(0)
0354 {
0355 m_bin_Svw.resize(parent::m_bin_number,0);
0356 m_bin_Sv2w.resize(parent::m_bin_number,0);
0357 }
0358
0359 p1(const std::string& a_title,bn_t aXnumber,TC aXmin,TC aXmax,TV aVmin,TV aVmax)
0360 :parent(a_title,aXnumber,aXmin,aXmax)
0361 ,m_cut_v(true)
0362 ,m_min_v(aVmin)
0363 ,m_max_v(aVmax)
0364 {
0365 m_bin_Svw.resize(parent::m_bin_number,0);
0366 m_bin_Sv2w.resize(parent::m_bin_number,0);
0367 }
0368
0369 p1(const std::string& a_title,const std::vector<TC>& a_edges)
0370 :parent(a_title,a_edges)
0371 ,m_cut_v(false)
0372 ,m_min_v(0)
0373 ,m_max_v(0)
0374 {
0375 m_bin_Svw.resize(parent::m_bin_number,0);
0376 m_bin_Sv2w.resize(parent::m_bin_number,0);
0377 }
0378
0379 p1(const std::string& a_title,const std::vector<TC>& a_edges,TV aVmin,TV aVmax)
0380 :parent(a_title,a_edges)
0381 ,m_cut_v(true)
0382 ,m_min_v(aVmin)
0383 ,m_max_v(aVmax)
0384 {
0385 m_bin_Svw.resize(parent::m_bin_number,0);
0386 m_bin_Sv2w.resize(parent::m_bin_number,0);
0387 }
0388
0389 virtual ~p1(){}
0390 public:
0391 p1(const p1& a_from)
0392 :parent(a_from)
0393 ,m_cut_v(a_from.m_cut_v)
0394 ,m_min_v(a_from.m_min_v)
0395 ,m_max_v(a_from.m_max_v)
0396 ,m_bin_Svw(a_from.m_bin_Svw)
0397 ,m_bin_Sv2w(a_from.m_bin_Sv2w)
0398 {}
0399 p1& operator=(const p1& a_from){
0400 parent::operator=(a_from);
0401 m_cut_v = a_from.m_cut_v;
0402 m_min_v = a_from.m_min_v;
0403 m_max_v = a_from.m_max_v;
0404 m_bin_Svw = a_from.m_bin_Svw;
0405 m_bin_Sv2w = a_from.m_bin_Sv2w;
0406 return *this;
0407 }
0408 public:
0409 bool configure(bn_t aXnumber,TC aXmin,TC aXmax){
0410 if(!parent::configure(aXnumber,aXmin,aXmax)) return false;
0411 m_bin_Svw.clear();
0412 m_bin_Sv2w.clear();
0413 m_bin_Svw.resize(parent::m_bin_number,0);
0414 m_bin_Sv2w.resize(parent::m_bin_number,0);
0415 m_cut_v = false;
0416 m_min_v = 0;
0417 m_max_v = 0;
0418 return true;
0419 }
0420 bool configure(const std::vector<TC>& a_edges) {
0421 if(!parent::configure(a_edges)) return false;
0422 m_bin_Svw.clear();
0423 m_bin_Sv2w.clear();
0424 m_bin_Svw.resize(parent::m_bin_number,0);
0425 m_bin_Sv2w.resize(parent::m_bin_number,0);
0426 m_cut_v = false;
0427 m_min_v = 0;
0428 m_max_v = 0;
0429 return true;
0430 }
0431 bool configure(bn_t aXnumber,TC aXmin,TC aXmax,TV aVmin,TV aVmax){
0432 if(!parent::configure(aXnumber,aXmin,aXmax)) return false;
0433 m_bin_Svw.clear();
0434 m_bin_Sv2w.clear();
0435 m_bin_Svw.resize(parent::m_bin_number,0);
0436 m_bin_Sv2w.resize(parent::m_bin_number,0);
0437 m_cut_v = true;
0438 m_min_v = aVmin;
0439 m_max_v = aVmax;
0440 return true;
0441 }
0442 bool configure(const std::vector<TC>& a_edges,TV aVmin,TV aVmax) {
0443 if(!parent::configure(a_edges)) return false;
0444 m_bin_Svw.clear();
0445 m_bin_Sv2w.clear();
0446 m_bin_Svw.resize(parent::m_bin_number,0);
0447 m_bin_Sv2w.resize(parent::m_bin_number,0);
0448 m_cut_v = true;
0449 m_min_v = aVmin;
0450 m_max_v = aVmax;
0451 return true;
0452 }
0453 public:
0454 const vs_t& bins_sum_vw() const {return m_bin_Svw;}
0455 const vs_t& bins_sum_v2w() const {return m_bin_Sv2w;}
0456
0457 TW get_Svw() const {
0458 TW sw = 0;
0459 for(TO ibin=0;ibin<parent::m_bin_number;ibin++) {
0460 if(!histo::is_out(parent::m_axes,ibin)) {
0461 sw += m_bin_Svw[ibin];
0462 }
0463 }
0464 return sw;
0465 }
0466 TW get_Sv2w() const {
0467 TW sw = 0;
0468 for(TO ibin=0;ibin<parent::m_bin_number;ibin++) {
0469 if(!histo::is_out(parent::m_axes,ibin)) {
0470 sw += m_bin_Sv2w[ibin];
0471 }
0472 }
0473 return sw;
0474 }
0475 protected:
0476 bool m_cut_v;
0477 TV m_min_v;
0478 TV m_max_v;
0479 vs_t m_bin_Svw;
0480 vs_t m_bin_Sv2w;
0481 };
0482
0483 }}
0484
0485 #endif
0486
