Geant4 Cross Reference |
1 // Copyright (C) 2010, Guy Barrand. All rights 2 // See the file tools.license for terms. 3 4 #ifndef tools_histo_p1 5 #define tools_histo_p1 6 7 #include "b1" 8 #include "profile_data" 9 10 namespace tools { 11 namespace histo { 12 13 //TC is for a coordinate. 14 //TO is for an offset/index to identify/find a 15 //TN is for a number of entries. 16 //TW is for a weight. 17 //TH is for a height. Should be the same as TV 18 //TV is for a value (in general same as TC). 19 20 template <class TC,class TO,class TN,class TW, 21 class p1 : public b1<TC,TO,TN,TW,TH> { 22 typedef b1<TC,TO,TN,TW,TH> parent; 23 public: 24 typedef profile_data<TC,TO,TN,TW,TV> pd_t; 25 typedef typename parent::bn_t bn_t; 26 typedef typename parent::axis_t axis_t; 27 typedef std::vector<TV> vs_t; 28 protected: 29 virtual TH get_bin_height(TO a_offset) const 30 return (parent::m_bin_Sw[a_offset] ? (m_bi 31 } 32 public: 33 bool equals(const p1& a_from,const TW& a_pre 34 if(!parent::equals(a_from,a_prec,a_fabs)) 35 if(m_cut_v!=a_from.m_cut_v) return false; 36 if(!numbers_are_equal(m_min_v,a_from.m_min 37 if(!numbers_are_equal(m_max_v,a_from.m_max 38 if(!vectors_are_equal(m_bin_Svw,a_from.m_b 39 if(!vectors_are_equal(m_bin_Sv2w,a_from.m_ 40 return true; 41 } 42 bool equals_TH(const p1& a_from,const TW& a_ 43 if(!parent::equals_TH(a_from,a_prec,a_fabs 44 if(m_cut_v!=a_from.m_cut_v) return false; 45 if(!numbers_are_equal(m_min_v,a_from.m_min 46 if(!numbers_are_equal(m_max_v,a_from.m_max 47 if(!vectors_are_equal(m_bin_Svw,a_from.m_b 48 if(!vectors_are_equal(m_bin_Sv2w,a_from.m_ 49 return true; 50 } 51 52 virtual TH bin_error(int aI) const { //TH sh 53 TO offset; 54 if(!parent::_find_offset(aI,offset)) retur 55 56 //FIXME Is it correct ? 57 // TProfile::GetBinError with kERRORMEAN m 58 // Stat_t cont = fArray[bin]; 59 // Stat_t sum = parent::m_bin_entries.fA 60 // Stat_t err2 = fSumw2.fArray[bin]; 61 // if (sum == 0) return 0; 62 // Stat_t eprim; 63 // Stat_t contsum = cont/sum; 64 // Stat_t eprim2 = TMath::Abs(err2/sum - 65 // eprim = TMath::Sqrt(eprim2); 66 // ... ??? 67 // if (fErrorMode == kERRORMEAN) return e 68 69 TW sw = parent::m_bin_Sw[offset]; //R 70 if(sw==0) return 0; 71 TV svw = m_bin_Svw[offset]; //ROOT cont 72 TV sv2w = m_bin_Sv2w[offset]; //ROOT err2 73 TV _mean = (svw / sw); //ROOT conts 74 TV _rms = ::sqrt(::fabs((sv2w/sw) - _mean 75 // rms = get_bin_rms_value. 76 return _rms/::sqrt(sw); //ROOT kERRORMEAN 77 } 78 79 public: 80 bool multiply(TW a_factor){ 81 if(!parent::base_multiply(a_factor)) retur 82 for(bn_t ibin=0;ibin<parent::m_bin_number; 83 m_bin_Svw[ibin] *= a_factor; 84 } 85 return true; 86 } 87 bool scale(TW a_factor) {return multiply(a_f 88 89 TV bin_Svw(int aI) const { 90 TO offset; 91 if(!parent::_find_offset(aI,offset)) retur 92 return m_bin_Svw[offset]; 93 } 94 TV bin_Sv2w(int aI) const { 95 TO offset; 96 if(!parent::_find_offset(aI,offset)) retur 97 return m_bin_Sv2w[offset]; 98 } 99 100 bool reset() { 101 parent::base_reset(); 102 for(bn_t ibin=0;ibin<parent::m_bin_number; 103 m_bin_Svw[ibin] = 0; 104 m_bin_Sv2w[ibin] = 0; 105 } 106 return true; 107 } 108 109 void copy_from_data(const pd_t& a_from) { 110 parent::base_from_data(a_from); 111 m_bin_Svw = a_from.m_bin_Svw; 112 m_bin_Sv2w = a_from.m_bin_Sv2w; 113 m_cut_v = a_from.m_cut_v; 114 m_min_v = a_from.m_min_v; 115 m_max_v = a_from.m_max_v; 116 } 117 pd_t get_histo_data() const { 118 pd_t hd(parent::dac()); 119 hd.m_is_profile = true; 120 hd.m_bin_Svw = m_bin_Svw; 121 hd.m_bin_Sv2w = m_bin_Sv2w; 122 hd.m_cut_v = m_cut_v; 123 hd.m_min_v = m_min_v; 124 hd.m_max_v = m_max_v; 125 return hd; 126 } 127 128 bool fill(TC aX,TV aV,TW aWeight = 1) { 129 if(parent::m_dimension!=1) return false; 130 131 if(m_cut_v) { 132 if( (aV<m_min_v) || (aV>=m_max_v) ) { 133 return true; 134 } 135 } 136 137 bn_t ibin; 138 if(!parent::m_axes[0].coord_to_absolute_in 139 TO offset = ibin; 140 141 parent::m_bin_entries[offset]++; 142 parent::m_bin_Sw[offset] += aWeight; 143 parent::m_bin_Sw2[offset] += aWeight * aWe 144 145 TC xw = aX * aWeight; 146 TC x2w = aX * xw; 147 parent::m_bin_Sxw[offset][0] += xw; 148 parent::m_bin_Sx2w[offset][0] += x2w; 149 150 bool inRange = true; 151 if(ibin==0) inRange = false; 152 else if(ibin==(parent::m_axes[0].m_number_ 153 154 parent::m_all_entries++; 155 if(inRange) { 156 // fast getters : 157 parent::m_in_range_entries++; 158 parent::m_in_range_Sw += aWeight; 159 parent::m_in_range_Sw2 += aWeight*aWeigh 160 161 parent::m_in_range_Sxw[0] += xw; 162 parent::m_in_range_Sx2w[0] += x2w; 163 } 164 165 // Profile part : 166 TV vw = aV * aWeight; 167 m_bin_Svw[offset] += vw; 168 m_bin_Sv2w[offset] += aV * vw; 169 170 return true; 171 } 172 173 bool set_bin_content(bn_t a_ibin,TN a_entrie 174 if(parent::m_dimension!=1) return false; 175 if(a_ibin>(parent::m_axes[0].m_number_of_b 176 177 bool inRange = true; 178 if(a_ibin==0) inRange = false; 179 else if(a_ibin==(parent::m_axes[0].m_numbe 180 181 TO offset = a_ibin; 182 183 parent::m_all_entries -= parent::m_bin_ent 184 if(inRange) { 185 parent::m_in_range_entries -= parent::m_ 186 parent::m_in_range_Sw -= parent::m_bin_S 187 parent::m_in_range_Sw2 -= parent::m_bin_ 188 parent::m_in_range_Sxw[0] -= parent::m_b 189 parent::m_in_range_Sx2w[0] -= parent::m_ 190 } 191 192 parent::m_bin_entries[offset] = a_entries; 193 parent::m_bin_Sw[offset] = a_Sw; 194 parent::m_bin_Sw2[offset] = a_Sw2; 195 196 parent::m_bin_Sxw[offset][0] = a_Sxw; 197 parent::m_bin_Sx2w[offset][0] = a_Sx2w; 198 199 parent::m_all_entries += a_entries; 200 if(inRange) { 201 parent::m_in_range_entries += a_entries; 202 parent::m_in_range_Sw += a_Sw; 203 parent::m_in_range_Sw2 += a_Sw2; 204 205 parent::m_in_range_Sxw[0] += a_Sxw; 206 parent::m_in_range_Sx2w[0] += a_Sx2w; 207 } 208 209 // Profile part : 210 m_bin_Svw[offset] = a_Svw; 211 m_bin_Sv2w[offset] = a_Sv2w; 212 213 return true; 214 } 215 216 bool get_bin_content(bn_t a_ibin,TN& a_entri 217 if(parent::m_dimension!=1) { 218 a_entries = 0;a_Sw = 0;a_Sw2 = 0;a_Sxw = 219 return false; 220 } 221 if(a_ibin>(parent::m_axes[0].m_number_of_b 222 a_entries = 0;a_Sw = 0;a_Sw2 = 0;a_Sxw = 223 return false; 224 } 225 226 TO offset = a_ibin; 227 228 a_entries = parent::m_bin_entries[offset]; 229 a_Sw = parent::m_bin_Sw[offset]; 230 a_Sw2 = parent::m_bin_Sw2[offset]; 231 232 a_Sxw = parent::m_bin_Sxw[offset][0]; 233 a_Sx2w = parent::m_bin_Sx2w[offset][0]; 234 235 // Profile part : 236 a_Svw = m_bin_Svw[offset]; 237 a_Sv2w = m_bin_Sv2w[offset]; 238 239 return true; 240 } 241 242 TV bin_rms_value(int aI) const { 243 TO offset; 244 if(!parent::_find_offset(aI,offset)) retur 245 TW sw = parent::m_bin_Sw[offset]; 246 if(sw==0) return 0; 247 TV svw = m_bin_Svw[offset]; 248 TV sv2w = m_bin_Sv2w[offset]; 249 TV _mean = (svw / sw); 250 return ::sqrt(::fabs((sv2w / sw) - _mean * 251 } 252 253 bool add(const p1& a_histo){ 254 parent::base_add(a_histo); 255 for(bn_t ibin=0;ibin<parent::m_bin_number; 256 m_bin_Svw[ibin] += a_histo.m_bin_Svw[ibi 257 m_bin_Sv2w[ibin] += a_histo.m_bin_Sv2w[i 258 } 259 return true; 260 } 261 bool subtract(const p1& a_histo){ 262 parent::base_subtract(a_histo); 263 for(bn_t ibin=0;ibin<parent::m_bin_number; 264 m_bin_Svw[ibin] -= a_histo.m_bin_Svw[ibi 265 m_bin_Sv2w[ibin] -= a_histo.m_bin_Sv2w[i 266 } 267 return true; 268 } 269 270 bool gather_bins(unsigned int a_factor) { // 271 if(!a_factor) return false; 272 273 // actual bin number must be a multiple of 274 275 const axis_t& _axis = parent::axis(); 276 277 bn_t n = _axis.bins(); 278 if(!n) return false; 279 280 bn_t new_n = n/a_factor; 281 if(a_factor*new_n!=n) return false; 282 283 p1* new_h = 0; 284 if(_axis.is_fixed_binning()) { 285 new_h = new p1(parent::m_title,new_n,_ax 286 } else { 287 const std::vector<TC>& _edges = _axis.ed 288 std::vector<TC> new_edges(new_n+1); 289 for(bn_t ibin=0;ibin<new_n;ibin++) { 290 new_edges[ibin] = _edges[ibin*a_factor 291 } 292 new_edges[new_n] = _edges[n]; //upper ed 293 new_h = new p1(parent::m_title,new_edges 294 } 295 if(!new_h) return false; 296 297 new_h->m_cut_v = m_cut_v; 298 new_h->m_min_v = m_min_v; 299 new_h->m_max_v = m_max_v; 300 301 bn_t offset,new_offset,offac; 302 for(bn_t ibin=0;ibin<new_n;ibin++) { 303 new_offset = ibin+1; 304 offset = a_factor*ibin+1; 305 for(unsigned int ifac=0;ifac<a_factor;if 306 offac = offset+ifac; 307 new_h->m_bin_entries[new_offset] += pa 308 new_h->m_bin_Sw[new_offset] += parent: 309 new_h->m_bin_Sw2[new_offset] += parent 310 new_h->m_bin_Sxw[new_offset][0] += par 311 new_h->m_bin_Sx2w[new_offset][0] += pa 312 313 new_h->m_bin_Svw[new_offset] += m_bin_ 314 new_h->m_bin_Sv2w[new_offset] += m_bin 315 } 316 } 317 318 //underflow : 319 new_offset = 0; 320 offac = 0; 321 new_h->m_bin_entries[new_offset] = parent: 322 new_h->m_bin_Sw[new_offset] = parent::m_bi 323 new_h->m_bin_Sw2[new_offset] = parent::m_b 324 new_h->m_bin_Sxw[new_offset][0] = parent:: 325 new_h->m_bin_Sx2w[new_offset][0] = parent: 326 new_h->m_bin_Svw[new_offset] = m_bin_Svw[o 327 new_h->m_bin_Sv2w[new_offset] = m_bin_Sv2w 328 329 //overflow : 330 new_offset = new_n+1; 331 offac = n+1; 332 new_h->m_bin_entries[new_offset] = parent: 333 new_h->m_bin_Sw[new_offset] = parent::m_bi 334 new_h->m_bin_Sw2[new_offset] = parent::m_b 335 new_h->m_bin_Sxw[new_offset][0] = parent:: 336 new_h->m_bin_Sx2w[new_offset][0] = parent: 337 new_h->m_bin_Svw[new_offset] = m_bin_Svw[o 338 new_h->m_bin_Sv2w[new_offset] = m_bin_Sv2w 339 340 *this = *new_h; 341 return true; 342 } 343 344 bool cut_v() const {return m_cut_v;} 345 TV min_v() const {return m_min_v;} 346 TV max_v() const {return m_max_v;} 347 348 public: 349 p1(const std::string& a_title,bn_t aXnumber, 350 :parent(a_title,aXnumber,aXmin,aXmax) 351 ,m_cut_v(false) 352 ,m_min_v(0) 353 ,m_max_v(0) 354 { 355 m_bin_Svw.resize(parent::m_bin_number,0); 356 m_bin_Sv2w.resize(parent::m_bin_number,0); 357 } 358 359 p1(const std::string& a_title,bn_t aXnumber, 360 :parent(a_title,aXnumber,aXmin,aXmax) 361 ,m_cut_v(true) 362 ,m_min_v(aVmin) 363 ,m_max_v(aVmax) 364 { 365 m_bin_Svw.resize(parent::m_bin_number,0); 366 m_bin_Sv2w.resize(parent::m_bin_number,0); 367 } 368 369 p1(const std::string& a_title,const std::vec 370 :parent(a_title,a_edges) 371 ,m_cut_v(false) 372 ,m_min_v(0) 373 ,m_max_v(0) 374 { 375 m_bin_Svw.resize(parent::m_bin_number,0); 376 m_bin_Sv2w.resize(parent::m_bin_number,0); 377 } 378 379 p1(const std::string& a_title,const std::vec 380 :parent(a_title,a_edges) 381 ,m_cut_v(true) 382 ,m_min_v(aVmin) 383 ,m_max_v(aVmax) 384 { 385 m_bin_Svw.resize(parent::m_bin_number,0); 386 m_bin_Sv2w.resize(parent::m_bin_number,0); 387 } 388 389 virtual ~p1(){} 390 public: 391 p1(const p1& a_from) 392 :parent(a_from) 393 ,m_cut_v(a_from.m_cut_v) 394 ,m_min_v(a_from.m_min_v) 395 ,m_max_v(a_from.m_max_v) 396 ,m_bin_Svw(a_from.m_bin_Svw) 397 ,m_bin_Sv2w(a_from.m_bin_Sv2w) 398 {} 399 p1& operator=(const p1& a_from){ 400 parent::operator=(a_from); 401 m_cut_v = a_from.m_cut_v; 402 m_min_v = a_from.m_min_v; 403 m_max_v = a_from.m_max_v; 404 m_bin_Svw = a_from.m_bin_Svw; 405 m_bin_Sv2w = a_from.m_bin_Sv2w; 406 return *this; 407 } 408 public: 409 bool configure(bn_t aXnumber,TC aXmin,TC aXm 410 if(!parent::configure(aXnumber,aXmin,aXmax 411 m_bin_Svw.clear(); 412 m_bin_Sv2w.clear(); 413 m_bin_Svw.resize(parent::m_bin_number,0); 414 m_bin_Sv2w.resize(parent::m_bin_number,0); 415 m_cut_v = false; 416 m_min_v = 0; 417 m_max_v = 0; 418 return true; 419 } 420 bool configure(const std::vector<TC>& a_edge 421 if(!parent::configure(a_edges)) return fal 422 m_bin_Svw.clear(); 423 m_bin_Sv2w.clear(); 424 m_bin_Svw.resize(parent::m_bin_number,0); 425 m_bin_Sv2w.resize(parent::m_bin_number,0); 426 m_cut_v = false; 427 m_min_v = 0; 428 m_max_v = 0; 429 return true; 430 } 431 bool configure(bn_t aXnumber,TC aXmin,TC aXm 432 if(!parent::configure(aXnumber,aXmin,aXmax 433 m_bin_Svw.clear(); 434 m_bin_Sv2w.clear(); 435 m_bin_Svw.resize(parent::m_bin_number,0); 436 m_bin_Sv2w.resize(parent::m_bin_number,0); 437 m_cut_v = true; 438 m_min_v = aVmin; 439 m_max_v = aVmax; 440 return true; 441 } 442 bool configure(const std::vector<TC>& a_edge 443 if(!parent::configure(a_edges)) return fal 444 m_bin_Svw.clear(); 445 m_bin_Sv2w.clear(); 446 m_bin_Svw.resize(parent::m_bin_number,0); 447 m_bin_Sv2w.resize(parent::m_bin_number,0); 448 m_cut_v = true; 449 m_min_v = aVmin; 450 m_max_v = aVmax; 451 return true; 452 } 453 public: 454 const vs_t& bins_sum_vw() const {return m_bi 455 const vs_t& bins_sum_v2w() const {return m_b 456 457 TW get_Svw() const { 458 TW sw = 0; 459 for(TO ibin=0;ibin<parent::m_bin_number;ib 460 if(!histo::is_out(parent::m_axes,ibin)) 461 sw += m_bin_Svw[ibin]; 462 } 463 } 464 return sw; 465 } 466 TW get_Sv2w() const { 467 TW sw = 0; 468 for(TO ibin=0;ibin<parent::m_bin_number;ib 469 if(!histo::is_out(parent::m_axes,ibin)) 470 sw += m_bin_Sv2w[ibin]; 471 } 472 } 473 return sw; 474 } 475 protected: 476 bool m_cut_v; 477 TV m_min_v; 478 TV m_max_v; 479 vs_t m_bin_Svw; 480 vs_t m_bin_Sv2w; 481 }; 482 483 }} 484 485 #endif 486