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_p2 00005 #define inlib_histo_p2 00006 00007 #include "b2" 00008 #include "profile_data" 00009 00010 namespace inlib { 00011 namespace histo { 00012 00013 template <class TC,class TN,class TW,class TH,class TV> 00014 class p2 : public b2<TC,TN,TW,TH> { 00015 typedef b2<TC,TN,TW,TH> parent; 00016 typedef typename base_histo<TC,TN,TW,TH>::bn_t bn_t; 00017 protected: 00018 virtual TH get_bin_height(int a_offset) const { 00019 return (parent::m_bin_Sw[a_offset] ? (m_bin_Svw[a_offset]/parent::m_bin_Sw[a_offset]):0); 00020 } 00021 public: 00022 00023 virtual TH bin_error(int aI,int aJ) const { //TH should be the same as TV 00024 if(parent::m_bin_number==0) return 0; 00025 bn_t ibin; 00026 if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) return 0; 00027 bn_t jbin; 00028 if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) return 0; 00029 bn_t offset = ibin + jbin * parent::m_axes[1].m_offset; 00030 00031 //FIXME Is it correct ? 00032 // TProfile::GetBinError with kERRORMEAN mode does : 00033 // Stat_t cont = fArray[bin]; //Svw (see TProfile::Fill) 00034 // Stat_t sum = parent::m_bin_entries.fArray[bin]; //Sw 00035 // Stat_t err2 = fSumw2.fArray[bin]; //Sv2w 00036 // if (sum == 0) return 0; 00037 // Stat_t eprim; 00038 // Stat_t contsum = cont/sum; 00039 // Stat_t eprim2 = TMath::Abs(err2/sum - contsum*contsum); 00040 // eprim = TMath::Sqrt(eprim2); 00041 // ... ??? 00042 // if (fErrorMode == kERRORMEAN) return eprim/TMath::Sqrt(sum); 00043 00044 TW sw = parent::m_bin_Sw[offset]; //ROOT sum 00045 if(sw==0) return 0; 00046 TV svw = m_bin_Svw[offset]; //ROOT cont 00047 TV sv2w = m_bin_Sv2w[offset]; //ROOT err2 00048 TV mean = (svw / sw); //ROOT contsum 00049 TV rms = ::sqrt(::fabs((sv2w/sw) - mean * mean)); //ROOT eprim 00050 // rms = get_bin_rms_value. 00051 return rms/::sqrt(sw); //ROOT kERRORMEAN mode returned value 00052 } 00053 00054 public: 00055 bool multiply(TW aFactor){ 00056 if(!parent::base_multiply(aFactor)) return false; 00057 for(bn_t ibin=0;ibin<parent::m_bin_number;ibin++) { 00058 m_bin_Svw[ibin] *= aFactor; 00059 } 00060 this->update_fast_getters(); 00061 return true; 00062 } 00063 bool scale(TW aFactor) {return multiply(aFactor);} 00064 00065 TV bin_Svw(int aI,int aJ) const { 00066 if(parent::m_bin_number==0) return 0; 00067 bn_t ibin; 00068 if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) return 0; 00069 bn_t jbin; 00070 if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) return 0; 00071 bn_t offset = ibin + jbin * parent::m_axes[1].m_offset; 00072 return m_bin_Svw[offset]; 00073 } 00074 TV bin_Sv2w(int aI,int aJ) const { 00075 if(parent::m_bin_number==0) return 0; 00076 bn_t ibin; 00077 if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) return 0; 00078 bn_t jbin; 00079 if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) return 0; 00080 bn_t offset = ibin + jbin * parent::m_axes[1].m_offset; 00081 return m_bin_Sv2w[offset]; 00082 } 00083 00084 bool reset() { 00085 parent::base_reset(); 00086 for(bn_t ibin=0;ibin<parent::m_bin_number;ibin++) { 00087 m_bin_Svw[ibin] = 0; 00088 m_bin_Sv2w[ibin] = 0; 00089 } 00090 this->update_fast_getters(); 00091 return true; 00092 } 00093 00094 void copy_from_data(const profile_data<TC,TN,TW,TV>& a_from) { 00095 parent::base_from_data(a_from); 00096 m_bin_Svw = a_from.m_bin_Svw; 00097 m_bin_Sv2w = a_from.m_bin_Sv2w; 00098 m_cut_v = a_from.m_cut_v; 00099 m_min_v = a_from.m_min_v; 00100 m_max_v = a_from.m_max_v; 00101 } 00102 profile_data<TC,TN,TW,TV> get_histo_data() const { 00103 profile_data<TC,TN,TW,TV> hd(parent::base_get_data()); 00104 hd.m_is_profile = true; 00105 hd.m_bin_Svw = m_bin_Svw; 00106 hd.m_bin_Sv2w = m_bin_Sv2w; 00107 hd.m_cut_v = m_cut_v; 00108 hd.m_min_v = m_min_v; 00109 hd.m_max_v = m_max_v; 00110 return hd; 00111 } 00112 00113 bool fill(TC aX,TC aY,TV aV,TW aWeight = 1) { 00114 //m_coords[0] = aX; 00115 //m_coords[1] = aY; 00116 //return fill_bin(m_coords,aV,aWeight); 00117 00118 if(m_cut_v) { 00119 if( (aV<m_min_v) || (aV>=m_max_v) ) { 00120 return true; 00121 } 00122 } 00123 00124 if(parent::m_dimension<=0) return false; 00125 00126 bn_t ibin,jbin; 00127 if(!parent::m_axes[0].coord_to_absolute_index(aX,ibin)) return false; 00128 if(!parent::m_axes[1].coord_to_absolute_index(aY,jbin)) return false; 00129 bn_t offset = ibin + jbin * parent::m_axes[1].m_offset; 00130 00131 parent::m_bin_entries[offset]++; 00132 parent::m_bin_Sw[offset] += aWeight; 00133 parent::m_bin_Sw2[offset] += aWeight * aWeight; 00134 00135 TC xw = aX * aWeight; 00136 TC x2w = aX * xw; 00137 parent::m_bin_Sxw[offset][0] += xw; 00138 parent::m_bin_Sx2w[offset][0] += x2w; 00139 00140 TC yw = aY * aWeight; 00141 TC y2w = aY * yw; 00142 parent::m_bin_Sxw[offset][1] += yw; 00143 parent::m_bin_Sx2w[offset][1] += y2w; 00144 00145 bool inRange = true; 00146 if(ibin==0) inRange = false; 00147 else if(ibin==(parent::m_axes[0].m_number_of_bins+1)) inRange = false; 00148 00149 if(jbin==0) inRange = false; 00150 else if(jbin==(parent::m_axes[1].m_number_of_bins+1)) inRange = false; 00151 00152 if(inRange) { 00153 parent::m_in_range_entries++; 00154 parent::m_in_range_Sw += aWeight; 00155 00156 parent::m_in_range_Sxw += xw; 00157 parent::m_in_range_Sx2w += x2w; 00158 00159 parent::m_in_range_Syw += yw; 00160 parent::m_in_range_Sy2w += y2w; 00161 } 00162 00163 // Profile part : 00164 TV vw = aV * aWeight; 00165 m_bin_Svw[offset] += vw; 00166 m_bin_Sv2w[offset] += aV * vw; 00167 00168 return true; 00169 } 00170 00171 TV bin_rms_value(int aI,int aJ) const { 00172 if(parent::m_bin_number==0) return 0; 00173 bn_t ibin; 00174 if(!parent::m_axes[0].in_range_to_absolute_index(aI,ibin)) return 0; 00175 bn_t jbin; 00176 if(!parent::m_axes[1].in_range_to_absolute_index(aJ,jbin)) return 0; 00177 bn_t offset = ibin + jbin * parent::m_axes[1].m_offset; 00178 00179 TW sw = parent::m_bin_Sw[offset]; 00180 if(sw==0) return 0; 00181 TV svw = m_bin_Svw[offset]; 00182 TV sv2w = m_bin_Sv2w[offset]; 00183 TV mean = (svw / sw); 00184 return ::sqrt(::fabs((sv2w / sw) - mean * mean)); 00185 } 00186 00187 bool add(const p2& a_histo){ 00188 base_add(a_histo); 00189 for(bn_t ibin=0;ibin<parent::m_bin_number;ibin++) { 00190 m_bin_Svw[ibin] += a_histo.m_bin_Svw[ibin]; 00191 m_bin_Sv2w[ibin] += a_histo.m_bin_Sv2w[ibin]; 00192 } 00193 this->update_fast_getters(); 00194 return true; 00195 } 00196 bool subtract(const p2& a_histo){ 00197 base_subtract(a_histo); 00198 for(bn_t ibin=0;ibin<parent::m_bin_number;ibin++) { 00199 m_bin_Svw[ibin] -= a_histo.m_bin_Svw[ibin]; 00200 m_bin_Sv2w[ibin] -= a_histo.m_bin_Sv2w[ibin]; 00201 } 00202 this->update_fast_getters(); 00203 return true; 00204 } 00205 00206 bool cut_v() const {return m_cut_v;} 00207 TV min_v() const {return m_min_v;} 00208 TV max_v() const {return m_max_v;} 00209 public: 00210 p2(const std::string& a_title, 00211 bn_t aXnumber,TC aXmin,TC aXmax, 00212 bn_t aYnumber,TC aYmin,TC aYmax) 00213 : parent(a_title,aXnumber,aXmin,aXmax,aYnumber,aYmin,aYmax) 00214 ,m_cut_v(false) 00215 ,m_min_v(0) 00216 ,m_max_v(0) 00217 { 00218 m_bin_Svw.resize(parent::m_bin_number,0); 00219 m_bin_Sv2w.resize(parent::m_bin_number,0); 00220 } 00221 00222 p2(const std::string& a_title, 00223 bn_t aXnumber,TC aXmin,TC aXmax, 00224 bn_t aYnumber,TC aYmin,TC aYmax, 00225 TV aVmin,TV aVmax) 00226 : parent(a_title,aXnumber,aXmin,aXmax,aYnumber,aYmin,aYmax) 00227 ,m_cut_v(true) 00228 ,m_min_v(aVmin) 00229 ,m_max_v(aVmax) 00230 { 00231 m_bin_Svw.resize(parent::m_bin_number,0); 00232 m_bin_Sv2w.resize(parent::m_bin_number,0); 00233 } 00234 00235 p2(const std::string& a_title, 00236 const std::vector<TC>& aEdgesX, 00237 const std::vector<TC>& aEdgesY) 00238 : parent(a_title,aEdgesX,aEdgesY) 00239 ,m_cut_v(false) 00240 ,m_min_v(0) 00241 ,m_max_v(0) 00242 { 00243 m_bin_Svw.resize(parent::m_bin_number,0); 00244 m_bin_Sv2w.resize(parent::m_bin_number,0); 00245 } 00246 00247 p2(const std::string& a_title, 00248 const std::vector<TC>& aEdgesX, 00249 const std::vector<TC>& aEdgesY, 00250 TV aVmin,TV aVmax) 00251 : parent(a_title,aEdgesX,aEdgesY) 00252 ,m_cut_v(true) 00253 ,m_min_v(aVmin) 00254 ,m_max_v(aVmax) 00255 { 00256 m_bin_Svw.resize(parent::m_bin_number,0); 00257 m_bin_Sv2w.resize(parent::m_bin_number,0); 00258 } 00259 00260 virtual ~p2(){} 00261 public: 00262 p2(const p2& a_from) 00263 : parent(a_from) 00264 ,m_cut_v(a_from.m_cut_v) 00265 ,m_min_v(a_from.m_min_v) 00266 ,m_max_v(a_from.m_max_v) 00267 ,m_bin_Svw(a_from.m_bin_Svw) 00268 ,m_bin_Sv2w(a_from.m_bin_Sv2w) 00269 {} 00270 p2& operator=(const p2& a_from){ 00271 parent::operator=(a_from); 00272 m_cut_v = a_from.m_cut_v; 00273 m_min_v = a_from.m_min_v; 00274 m_max_v = a_from.m_max_v; 00275 m_bin_Svw = a_from.m_bin_Svw; 00276 m_bin_Sv2w = a_from.m_bin_Sv2w; 00277 return *this; 00278 } 00279 public: 00280 const std::vector<TV>& bins_sum_vw() const {return m_bin_Svw;} 00281 const std::vector<TV>& bins_sum_v2w() const {return m_bin_Sv2w;} 00282 protected: 00283 bool m_cut_v; 00284 TV m_min_v; 00285 TV m_max_v; 00286 std::vector<TV> m_bin_Svw; 00287 std::vector<TV> m_bin_Sv2w; 00288 }; 00289 00290 }} 00291 00292 #endif 00293