inlib
1.2.0
|
00001 // Copyright (C) 2010, Guy Barrand. All rights reserved. 00002 // See the file inlib.license for terms. 00003 00004 #ifndef inlib_histo_p1 00005 #define inlib_histo_p1 00006 00007 #include "b1" 00008 #include "profile_data" 00009 00010 namespace inlib { 00011 namespace histo { 00012 00013 //TC is for a coordinate. 00014 //TW is for a weight. 00015 //TH is for a height. Should be the same as TV. 00016 //TV is for a value (in general same as TC). 00017 00018 template <class TC,class TN,class TW,class TH,class TV> 00019 class p1 : public b1<TC,TN,TW,TH> { 00020 typedef b1<TC,TN,TW,TH> parent; 00021 typedef typename base_histo<TC,TN,TW,TH>::bn_t bn_t; 00022 protected: 00023 virtual TH get_bin_height(int a_offset) const { 00024 return (parent::m_bin_Sw[a_offset] ? (m_bin_Svw[a_offset]/parent::m_bin_Sw[a_offset]):0); 00025 } 00026 public: 00027 00028 virtual TH bin_error(int aI) const { //TH should be the same as TV 00029 if(parent::m_bin_number==0) return 0; 00030 bn_t offset; 00031 if(!parent::m_axes[0].in_range_to_absolute_index(aI,offset)) return 0; 00032 00033 //FIXME Is it correct ? 00034 // TProfile::GetBinError with kERRORMEAN mode does : 00035 // Stat_t cont = fArray[bin]; //Svw (see TProfile::Fill) 00036 // Stat_t sum = parent::m_bin_entries.fArray[bin]; //Sw 00037 // Stat_t err2 = fSumw2.fArray[bin]; //Sv2w 00038 // if (sum == 0) return 0; 00039 // Stat_t eprim; 00040 // Stat_t contsum = cont/sum; 00041 // Stat_t eprim2 = TMath::Abs(err2/sum - contsum*contsum); 00042 // eprim = TMath::Sqrt(eprim2); 00043 // ... ??? 00044 // if (fErrorMode == kERRORMEAN) return eprim/TMath::Sqrt(sum); 00045 00046 TW sw = parent::m_bin_Sw[offset]; //ROOT sum 00047 if(sw==0) return 0; 00048 TV svw = m_bin_Svw[offset]; //ROOT cont 00049 TV sv2w = m_bin_Sv2w[offset]; //ROOT err2 00050 TV mean = (svw / sw); //ROOT contsum 00051 TV rms = ::sqrt(::fabs((sv2w/sw) - mean * mean)); //ROOT eprim 00052 // rms = get_bin_rms_value. 00053 return rms/::sqrt(sw); //ROOT kERRORMEAN mode returned value 00054 } 00055 00056 public: 00057 bool multiply(TW aFactor){ 00058 if(!parent::base_multiply(aFactor)) return false; 00059 for(bn_t ibin=0;ibin<parent::m_bin_number;ibin++) { 00060 m_bin_Svw[ibin] *= aFactor; 00061 } 00062 this->update_fast_getters(); 00063 return true; 00064 } 00065 bool scale(TW aFactor) {return multiply(aFactor);} 00066 00067 TV bin_Svw(int aI) const { 00068 if(parent::m_bin_number==0) return 0; 00069 bn_t offset; 00070 if(!parent::m_axes[0].in_range_to_absolute_index(aI,offset)) return 0; 00071 return parent::m_bin_Svw(offset); 00072 } 00073 TV bin_Sv2w(int aI) const { 00074 if(parent::m_bin_number==0) return 0; 00075 bn_t offset; 00076 if(!parent::m_axes[0].in_range_to_absolute_index(aI,offset)) return 0; 00077 return parent::m_bin_Sv2w(offset); 00078 } 00079 00080 bool reset() { 00081 parent::base_reset(); 00082 for(bn_t ibin=0;ibin<parent::m_bin_number;ibin++) { 00083 m_bin_Svw[ibin] = 0; 00084 m_bin_Sv2w[ibin] = 0; 00085 } 00086 this->update_fast_getters(); 00087 return true; 00088 } 00089 00090 void copy_from_data(const profile_data<TC,TN,TW,TV>& a_from) { 00091 parent::base_from_data(a_from); 00092 m_bin_Svw = a_from.m_bin_Svw; 00093 m_bin_Sv2w = a_from.m_bin_Sv2w; 00094 m_cut_v = a_from.m_cut_v; 00095 m_min_v = a_from.m_min_v; 00096 m_max_v = a_from.m_max_v; 00097 } 00098 profile_data<TC,TN,TW,TV> get_histo_data() const { 00099 profile_data<TC,TN,TW,TV> hd(parent::base_get_data()); 00100 hd.m_is_profile = true; 00101 hd.m_bin_Svw = m_bin_Svw; 00102 hd.m_bin_Sv2w = m_bin_Sv2w; 00103 hd.m_cut_v = m_cut_v; 00104 hd.m_min_v = m_min_v; 00105 hd.m_max_v = m_max_v; 00106 return hd; 00107 } 00108 00109 bool fill(TC aX,TV aV,TW aWeight = 1) { 00110 //m_coords[0] = aX; 00111 //return fill_bin(m_coords,aV,aWeight); 00112 00113 if(!parent::m_dimension) return false; 00114 if(m_cut_v) { 00115 if( (aV<m_min_v) || (aV>=m_max_v) ) { 00116 return true; 00117 } 00118 } 00119 00120 bn_t offset; 00121 if(!parent::m_axes[0].coord_to_absolute_index(aX,offset)) return false; 00122 00123 parent::m_bin_entries[offset]++; 00124 parent::m_bin_Sw[offset] += aWeight; 00125 parent::m_bin_Sw2[offset] += aWeight * aWeight; 00126 00127 TC xw = aX * aWeight; 00128 TC x2w = aX * xw; 00129 parent::m_bin_Sxw[offset][0] += xw; 00130 parent::m_bin_Sx2w[offset][0] += x2w; 00131 00132 // Profile part : 00133 TV vw = aV * aWeight; 00134 m_bin_Svw[offset] += vw; 00135 m_bin_Sv2w[offset] += aV * vw; 00136 00137 return true; 00138 } 00139 00140 TV bin_rms_value(int aI) const { 00141 if(parent::m_bin_number==0) return 0; 00142 bn_t offset; 00143 if(!parent::m_axes[0].in_range_to_absolute_index(aI,offset)) return 0; 00144 00145 TW sw = parent::m_bin_Sw[offset]; 00146 if(sw==0) return 0; 00147 TV svw = m_bin_Svw[offset]; 00148 TV sv2w = m_bin_Sv2w[offset]; 00149 TV mean = (svw / sw); 00150 return ::sqrt(::fabs((sv2w / sw) - mean * mean)); 00151 } 00152 00153 bool add(const p1& a_histo){ 00154 parent::base_add(a_histo); 00155 for(bn_t ibin=0;ibin<parent::m_bin_number;ibin++) { 00156 m_bin_Svw[ibin] += a_histo.m_bin_Svw[ibin]; 00157 m_bin_Sv2w[ibin] += a_histo.m_bin_Sv2w[ibin]; 00158 } 00159 this->update_fast_getters(); 00160 return true; 00161 } 00162 bool subtract(const p1& a_histo){ 00163 parent::base_subtract(a_histo); 00164 for(bn_t ibin=0;ibin<parent::m_bin_number;ibin++) { 00165 m_bin_Svw[ibin] -= a_histo.m_bin_Svw[ibin]; 00166 m_bin_Sv2w[ibin] -= a_histo.m_bin_Sv2w[ibin]; 00167 } 00168 this->update_fast_getters(); 00169 return true; 00170 } 00171 00172 bool gather_bins(unsigned int a_factor) { //for exa 2,3. 00173 if(!a_factor) return false; 00174 00175 // actual bin number must be a multiple of a_factor. 00176 typedef typename b1<TC,TN,TW,TH>::bn_t bn_t; 00177 00178 const histo::axis<TC>& _axis = this->axis(); 00179 00180 bn_t n = _axis.bins(); 00181 if(!n) return false; 00182 00183 bn_t new_n = n/a_factor; 00184 if(a_factor*new_n!=n) return false; 00185 00186 p1* new_h = 0; 00187 if(_axis.is_fixed_binning()) { 00188 new_h = new p1(parent::m_title, 00189 new_n,_axis.lower_edge(),_axis.upper_edge()); 00190 } else { 00191 const std::vector<TC>& _edges = _axis.edges(); 00192 std::vector<TC> new_edges(new_n+1); 00193 for(bn_t ibin=0;ibin<new_n;ibin++) { 00194 new_edges[ibin] = _edges[ibin*a_factor]; 00195 } 00196 new_edges[new_n] = _edges[n]; //upper edge. 00197 new_h = new p1(parent::m_title,new_edges); 00198 } 00199 if(!new_h) return false; 00200 00201 new_h->m_cut_v = m_cut_v; 00202 new_h->m_min_v = m_min_v; 00203 new_h->m_max_v = m_max_v; 00204 00205 bn_t offset,new_offset,offac; 00206 for(bn_t ibin=0;ibin<new_n;ibin++) { 00207 new_offset = ibin+1; 00208 offset = a_factor*ibin+1; 00209 for(unsigned int ifac=0;ifac<a_factor;ifac++) { 00210 offac = offset+ifac; 00211 new_h->m_bin_entries[new_offset] += parent::m_bin_entries[offac]; 00212 new_h->m_bin_Sw[new_offset] += parent::m_bin_Sw[offac]; 00213 new_h->m_bin_Sw2[new_offset] += parent::m_bin_Sw2[offac]; 00214 new_h->m_bin_Sxw[new_offset][0] += parent::m_bin_Sxw[offac][0]; 00215 new_h->m_bin_Sx2w[new_offset][0] += parent::m_bin_Sx2w[offac][0]; 00216 00217 new_h->m_bin_Svw[new_offset] += m_bin_Svw[offac]; 00218 new_h->m_bin_Sv2w[new_offset] += m_bin_Sv2w[offac]; 00219 } 00220 } 00221 00222 //underflow : 00223 new_offset = 0; 00224 offac = 0; 00225 new_h->m_bin_entries[new_offset] = parent::m_bin_entries[offac]; 00226 new_h->m_bin_Sw[new_offset] = parent::m_bin_Sw[offac]; 00227 new_h->m_bin_Sw2[new_offset] = parent::m_bin_Sw2[offac]; 00228 new_h->m_bin_Sxw[new_offset][0] = parent::m_bin_Sxw[offac][0]; 00229 new_h->m_bin_Sx2w[new_offset][0] = parent::m_bin_Sx2w[offac][0]; 00230 new_h->m_bin_Svw[new_offset] = m_bin_Svw[offac]; 00231 new_h->m_bin_Sv2w[new_offset] = m_bin_Sv2w[offac]; 00232 00233 //overflow : 00234 new_offset = new_n+1; 00235 offac = n+1; 00236 new_h->m_bin_entries[new_offset] = parent::m_bin_entries[offac]; 00237 new_h->m_bin_Sw[new_offset] = parent::m_bin_Sw[offac]; 00238 new_h->m_bin_Sw2[new_offset] = parent::m_bin_Sw2[offac]; 00239 new_h->m_bin_Sxw[new_offset][0] = parent::m_bin_Sxw[offac][0]; 00240 new_h->m_bin_Sx2w[new_offset][0] = parent::m_bin_Sx2w[offac][0]; 00241 new_h->m_bin_Svw[new_offset] = m_bin_Svw[offac]; 00242 new_h->m_bin_Sv2w[new_offset] = m_bin_Sv2w[offac]; 00243 00244 *this = *new_h; 00245 return true; 00246 } 00247 00248 bool cut_v() const {return m_cut_v;} 00249 TV min_v() const {return m_min_v;} 00250 TV max_v() const {return m_max_v;} 00251 00252 public: 00253 p1(const std::string& a_title, 00254 bn_t aXnumber,TC aXmin,TC aXmax) 00255 : parent(a_title,aXnumber,aXmin,aXmax) 00256 ,m_cut_v(false) 00257 ,m_min_v(0) 00258 ,m_max_v(0) 00259 { 00260 m_bin_Svw.resize(parent::m_bin_number,0); 00261 m_bin_Sv2w.resize(parent::m_bin_number,0); 00262 } 00263 00264 p1(const std::string& a_title, 00265 bn_t aXnumber,TC aXmin,TC aXmax, 00266 TV aVmin,TV aVmax) 00267 : parent(a_title,aXnumber,aXmin,aXmax) 00268 ,m_cut_v(true) 00269 ,m_min_v(aVmin) 00270 ,m_max_v(aVmax) 00271 { 00272 m_bin_Svw.resize(parent::m_bin_number,0); 00273 m_bin_Sv2w.resize(parent::m_bin_number,0); 00274 } 00275 00276 p1(const std::string& a_title, 00277 const std::vector<TC>& aEdges) 00278 : parent(a_title,aEdges) 00279 ,m_cut_v(false) 00280 ,m_min_v(0) 00281 ,m_max_v(0) 00282 { 00283 m_bin_Svw.resize(parent::m_bin_number,0); 00284 m_bin_Sv2w.resize(parent::m_bin_number,0); 00285 } 00286 00287 p1(const std::string& a_title, 00288 const std::vector<TC>& aEdges, 00289 TV aVmin,TV aVmax) 00290 : parent(a_title,aEdges) 00291 ,m_cut_v(true) 00292 ,m_min_v(aVmin) 00293 ,m_max_v(aVmax) 00294 { 00295 m_bin_Svw.resize(parent::m_bin_number,0); 00296 m_bin_Sv2w.resize(parent::m_bin_number,0); 00297 } 00298 00299 virtual ~p1(){} 00300 public: 00301 p1(const p1& a_from) 00302 : parent(a_from) 00303 ,m_cut_v(a_from.m_cut_v) 00304 ,m_min_v(a_from.m_min_v) 00305 ,m_max_v(a_from.m_max_v) 00306 ,m_bin_Svw(a_from.m_bin_Svw) 00307 ,m_bin_Sv2w(a_from.m_bin_Sv2w) 00308 {} 00309 p1& operator=(const p1& a_from){ 00310 parent::operator=(a_from); 00311 m_cut_v = a_from.m_cut_v; 00312 m_min_v = a_from.m_min_v; 00313 m_max_v = a_from.m_max_v; 00314 m_bin_Svw = a_from.m_bin_Svw; 00315 m_bin_Sv2w = a_from.m_bin_Sv2w; 00316 return *this; 00317 } 00318 public: 00319 const std::vector<TV>& bins_sum_vw() const {return m_bin_Svw;} 00320 const std::vector<TV>& bins_sum_v2w() const {return m_bin_Sv2w;} 00321 protected: 00322 bool m_cut_v; 00323 TV m_min_v; 00324 TV m_max_v; 00325 std::vector<TV> m_bin_Svw; 00326 std::vector<TV> m_bin_Sv2w; 00327 }; 00328 00329 }} 00330 00331 #endif 00332