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 // rms = get_bin_rms_value.
0068 return rms/::sqrt(sw); //ROOT kERRORMEAN mode returned value
0069 }
0070
0071 public:
0072 bool multiply(TW a_factor){
0073 if(!parent::base_multiply(a_factor)) return false;
0074 for(bn_t ibin=0;ibin<parent::m_bin_number;ibin++) {
0075 m_bin_Svw[ibin] *= a_factor;
0076 }
0077 return true;
0078 }
0079 bool scale(TW a_factor) {return multiply(a_factor);}
0080
0081 TV bin_Svw(int aI,int aJ) const {
0082 TO offset;
0083 if(!parent::_find_offset(aI,aJ,offset)) return 0;
0084 return m_bin_Svw[offset];
0085 }
0086 TV bin_Sv2w(int aI,int aJ) const {
0087 TO offset;
0088 if(!parent::_find_offset(aI,aJ,offset)) return 0;
0089 return m_bin_Sv2w[offset];
0090 }
0091
0092 bool reset() {
0093 parent::base_reset();
0094 for(bn_t ibin=0;ibin<parent::m_bin_number;ibin++) {
0095 m_bin_Svw[ibin] = 0;
0096 m_bin_Sv2w[ibin] = 0;
0097 }
0098 return true;
0099 }
0100
0101 void copy_from_data(const pd_t& a_from) {
0102 parent::base_from_data(a_from);
0103 m_bin_Svw = a_from.m_bin_Svw;
0104 m_bin_Sv2w = a_from.m_bin_Sv2w;
0105 m_cut_v = a_from.m_cut_v;
0106 m_min_v = a_from.m_min_v;
0107 m_max_v = a_from.m_max_v;
0108 }
0109 pd_t get_histo_data() const {
0110 pd_t hd(parent::dac());
0111 hd.m_is_profile = true;
0112 hd.m_bin_Svw = m_bin_Svw;
0113 hd.m_bin_Sv2w = m_bin_Sv2w;
0114 hd.m_cut_v = m_cut_v;
0115 hd.m_min_v = m_min_v;
0116 hd.m_max_v = m_max_v;
0117 return hd;
0118 }
0119
0120 bool fill(TC aX,TC aY,TV aV,TW aWeight = 1) {
0121 if(parent::m_dimension!=2) return false;
0122
0123 if(m_cut_v) {
0124 if( (aV<m_min_v) || (aV>=m_max_v) ) {
0125 return true;
0126 }
0127 }
0128
0129 bn_t ibin,jbin;
0130 if(!parent::m_axes[0].coord_to_absolute_index(aX,ibin)) return false;
0131 if(!parent::m_axes[1].coord_to_absolute_index(aY,jbin)) return false;
0132 bn_t offset = ibin + jbin * parent::m_axes[1].m_offset;
0133
0134 parent::m_bin_entries[offset]++;
0135 parent::m_bin_Sw[offset] += aWeight;
0136 parent::m_bin_Sw2[offset] += aWeight * aWeight;
0137
0138 TC xw = aX * aWeight;
0139 TC x2w = aX * xw;
0140 parent::m_bin_Sxw[offset][0] += xw;
0141 parent::m_bin_Sx2w[offset][0] += x2w;
0142
0143 TC yw = aY * aWeight;
0144 TC y2w = aY * yw;
0145 parent::m_bin_Sxw[offset][1] += yw;
0146 parent::m_bin_Sx2w[offset][1] += y2w;
0147
0148 bool inRange = true;
0149 if(ibin==0) inRange = false;
0150 else if(ibin==(parent::m_axes[0].m_number_of_bins+1)) inRange = false;
0151
0152 if(jbin==0) inRange = false;
0153 else if(jbin==(parent::m_axes[1].m_number_of_bins+1)) inRange = false;
0154
0155 parent::m_all_entries++;
0156 if(inRange) {
0157 parent::m_in_range_plane_Sxyw[0] += aX * aY * aWeight;
0158
0159 // fast getters :
0160 parent::m_in_range_entries++;
0161 parent::m_in_range_Sw += aWeight;
0162 parent::m_in_range_Sw2 += aWeight*aWeight;
0163
0164 parent::m_in_range_Sxw[0] += xw;
0165 parent::m_in_range_Sx2w[0] += x2w;
0166
0167 parent::m_in_range_Sxw[1] += yw;
0168 parent::m_in_range_Sx2w[1] += y2w;
0169 }
0170
0171 // Profile part :
0172 TV vw = aV * aWeight;
0173 m_bin_Svw[offset] += vw;
0174 m_bin_Sv2w[offset] += aV * vw;
0175
0176 return true;
0177 }
0178
0179 TV bin_rms_value(int aI,int aJ) const {
0180 TO offset;
0181 if(!parent::_find_offset(aI,aJ,offset)) return 0;
0182
0183 TW sw = parent::m_bin_Sw[offset];
0184 if(sw==0) return 0;
0185 TV svw = m_bin_Svw[offset];
0186 TV sv2w = m_bin_Sv2w[offset];
0187 TV mean = (svw / sw);
0188 return ::sqrt(::fabs((sv2w / sw) - mean * mean));
0189 }
0190
0191 bool add(const p2& a_histo){
0192 parent::base_add(a_histo);
0193 for(bn_t ibin=0;ibin<parent::m_bin_number;ibin++) {
0194 m_bin_Svw[ibin] += a_histo.m_bin_Svw[ibin];
0195 m_bin_Sv2w[ibin] += a_histo.m_bin_Sv2w[ibin];
0196 }
0197 return true;
0198 }
0199 bool subtract(const p2& a_histo){
0200 parent::base_subtract(a_histo);
0201 for(bn_t ibin=0;ibin<parent::m_bin_number;ibin++) {
0202 m_bin_Svw[ibin] -= a_histo.m_bin_Svw[ibin];
0203 m_bin_Sv2w[ibin] -= a_histo.m_bin_Sv2w[ibin];
0204 }
0205 return true;
0206 }
0207
0208 bool cut_v() const {return m_cut_v;}
0209 TV min_v() const {return m_min_v;}
0210 TV max_v() const {return m_max_v;}
0211 public:
0212 p2(const std::string& a_title,
0213 bn_t aXnumber,TC aXmin,TC aXmax,
0214 bn_t aYnumber,TC aYmin,TC aYmax)
0215 :parent(a_title,aXnumber,aXmin,aXmax,aYnumber,aYmin,aYmax)
0216 ,m_cut_v(false)
0217 ,m_min_v(0)
0218 ,m_max_v(0)
0219 {
0220 m_bin_Svw.resize(parent::m_bin_number,0);
0221 m_bin_Sv2w.resize(parent::m_bin_number,0);
0222 }
0223
0224 p2(const std::string& a_title,
0225 bn_t aXnumber,TC aXmin,TC aXmax,
0226 bn_t aYnumber,TC aYmin,TC aYmax,
0227 TV aVmin,TV aVmax)
0228 :parent(a_title,aXnumber,aXmin,aXmax,aYnumber,aYmin,aYmax)
0229 ,m_cut_v(true)
0230 ,m_min_v(aVmin)
0231 ,m_max_v(aVmax)
0232 {
0233 m_bin_Svw.resize(parent::m_bin_number,0);
0234 m_bin_Sv2w.resize(parent::m_bin_number,0);
0235 }
0236
0237 p2(const std::string& a_title,
0238 const std::vector<TC>& a_edges_x,
0239 const std::vector<TC>& a_edges_y)
0240 :parent(a_title,a_edges_x,a_edges_y)
0241 ,m_cut_v(false)
0242 ,m_min_v(0)
0243 ,m_max_v(0)
0244 {
0245 m_bin_Svw.resize(parent::m_bin_number,0);
0246 m_bin_Sv2w.resize(parent::m_bin_number,0);
0247 }
0248
0249 p2(const std::string& a_title,
0250 const std::vector<TC>& a_edges_x,
0251 const std::vector<TC>& a_edges_y,
0252 TV aVmin,TV aVmax)
0253 :parent(a_title,a_edges_x,a_edges_y)
0254 ,m_cut_v(true)
0255 ,m_min_v(aVmin)
0256 ,m_max_v(aVmax)
0257 {
0258 m_bin_Svw.resize(parent::m_bin_number,0);
0259 m_bin_Sv2w.resize(parent::m_bin_number,0);
0260 }
0261
0262 virtual ~p2(){}
0263 public:
0264 p2(const p2& a_from)
0265 :parent(a_from)
0266 ,m_cut_v(a_from.m_cut_v)
0267 ,m_min_v(a_from.m_min_v)
0268 ,m_max_v(a_from.m_max_v)
0269 ,m_bin_Svw(a_from.m_bin_Svw)
0270 ,m_bin_Sv2w(a_from.m_bin_Sv2w)
0271 {}
0272 p2& operator=(const p2& a_from){
0273 parent::operator=(a_from);
0274 m_cut_v = a_from.m_cut_v;
0275 m_min_v = a_from.m_min_v;
0276 m_max_v = a_from.m_max_v;
0277 m_bin_Svw = a_from.m_bin_Svw;
0278 m_bin_Sv2w = a_from.m_bin_Sv2w;
0279 return *this;
0280 }
0281 public:
0282 bool configure(bn_t aXnumber,TC aXmin,TC aXmax,bn_t aYnumber,TC aYmin,TC aYmax){
0283 if(!parent::configure(aXnumber,aXmin,aXmax,aYnumber,aYmin,aYmax)) return false;
0284 m_bin_Svw.clear();
0285 m_bin_Sv2w.clear();
0286 m_bin_Svw.resize(parent::m_bin_number,0);
0287 m_bin_Sv2w.resize(parent::m_bin_number,0);
0288 m_cut_v = false;
0289 m_min_v = 0;
0290 m_max_v = 0;
0291 return true;
0292 }
0293 bool configure(const std::vector<TC>& a_edges_x,const std::vector<TC>& a_edges_y) {
0294 if(!parent::configure(a_edges_x,a_edges_y)) return false;
0295 m_bin_Svw.clear();
0296 m_bin_Sv2w.clear();
0297 m_bin_Svw.resize(parent::m_bin_number,0);
0298 m_bin_Sv2w.resize(parent::m_bin_number,0);
0299 m_cut_v = false;
0300 m_min_v = 0;
0301 m_max_v = 0;
0302 return true;
0303 }
0304 bool configure(bn_t aXnumber,TC aXmin,TC aXmax,bn_t aYnumber,TC aYmin,TC aYmax,TV aVmin,TV aVmax){
0305 if(!parent::configure(aXnumber,aXmin,aXmax,aYnumber,aYmin,aYmax)) return false;
0306 m_bin_Svw.clear();
0307 m_bin_Sv2w.clear();
0308 m_bin_Svw.resize(parent::m_bin_number,0);
0309 m_bin_Sv2w.resize(parent::m_bin_number,0);
0310 m_cut_v = true;
0311 m_min_v = aVmin;
0312 m_max_v = aVmax;
0313 return true;
0314 }
0315 bool configure(const std::vector<TC>& a_edges_x,const std::vector<TC>& a_edges_y,TV aVmin,TV aVmax) {
0316 if(!parent::configure(a_edges_x,a_edges_y)) return false;
0317 m_bin_Svw.clear();
0318 m_bin_Sv2w.clear();
0319 m_bin_Svw.resize(parent::m_bin_number,0);
0320 m_bin_Sv2w.resize(parent::m_bin_number,0);
0321 m_cut_v = true;
0322 m_min_v = aVmin;
0323 m_max_v = aVmax;
0324 return true;
0325 }
0326 public:
0327 const vs_t& bins_sum_vw() const {return m_bin_Svw;}
0328 const vs_t& bins_sum_v2w() const {return m_bin_Sv2w;}
0329
0330 TW get_Svw() const {
0331 TW sw = 0;
0332 for(TO ibin=0;ibin<parent::m_bin_number;ibin++) {
0333 if(!histo::is_out(parent::m_axes,ibin)) {
0334 sw += m_bin_Svw[ibin];
0335 }
0336 }
0337 return sw;
0338 }
0339 TW get_Sv2w() const {
0340 TW sw = 0;
0341 for(TO ibin=0;ibin<parent::m_bin_number;ibin++) {
0342 if(!histo::is_out(parent::m_axes,ibin)) {
0343 sw += m_bin_Sv2w[ibin];
0344 }
0345 }
0346 return sw;
0347 }
0348 protected:
0349 bool m_cut_v;
0350 TV m_min_v;
0351 TV m_max_v;
0352 vs_t m_bin_Svw;
0353 vs_t m_bin_Sv2w;
0354 };
0355
0356 }}
0357
0358 #endif
0359
