Geant4 Cross Reference |
1 // Copyright (C) 2010, Guy Barrand. All rights 2 // See the file tools.license for terms. 3 4 #ifndef tools_histo_c3d 5 #define tools_histo_c3d 6 7 #include "base_cloud" 8 9 #include "../mnmx" 10 11 #include "h3d" 12 13 namespace tools { 14 namespace histo { 15 16 class c3d : public base_cloud { 17 public: 18 static const std::string& s_class() { 19 static const std::string s_v("tools::histo 20 return s_v; 21 } 22 public: 23 bool set_title(const std::string&); 24 unsigned int dimension() const {return 3;} 25 bool reset(); 26 unsigned int entries() const; 27 public: 28 double sum_of_weights() const; 29 bool convert_to_histogram(); 30 bool is_converted() const; 31 bool scale(double); 32 public: 33 bool fill(double,double,double,double = 1); 34 double lower_edge_x() const; 35 double upper_edge_x() const; 36 double lower_edge_y() const; 37 double upper_edge_y() const; 38 double lower_edge_z() const; 39 double upper_edge_z() const; 40 double value_x(unsigned int) const; 41 double value_y(unsigned int) const; 42 double value_z(unsigned int) const; 43 double weight(unsigned int) const; 44 double mean_x() const; 45 double mean_y() const; 46 double mean_z() const; 47 double rms_x() const; 48 double rms_y() const; 49 double rms_z() const; 50 bool convert(unsigned int,double,double, 51 unsigned int,double,double, 52 unsigned int,double,double); 53 bool convert(const std::vector<double>&, 54 const std::vector<doubl 55 const std::vector<doubl 56 const histo::h3d& histogram() const; 57 bool fill_histogram(histo::h3d& a_histo) con 58 size_t number = m_xs.size(); 59 for(size_t index=0;index<number;index++) { 60 if(!a_histo.fill(m_xs[index],m_ys[index] 61 } 62 return true; 63 } 64 bool set_conversion_parameters(unsigned int, 65 unsigned int, 66 unsigned int, 67 68 bool set_histogram(h3d* a_histo){ //we take 69 reset(); 70 m_histo = a_histo; 71 return true; 72 } 73 public: 74 c3d(); 75 c3d(const std::string&,int=base_cloud::UNLIM 76 virtual ~c3d(){delete m_histo;} 77 public: 78 c3d(const c3d& a_from) 79 :base_cloud(a_from) 80 ,m_xs(a_from.m_xs) 81 ,m_ys(a_from.m_ys) 82 ,m_zs(a_from.m_zs) 83 ,m_lower_x(a_from.m_lower_x) 84 ,m_upper_x(a_from.m_upper_x) 85 ,m_lower_y(a_from.m_lower_y) 86 ,m_upper_y(a_from.m_upper_y) 87 ,m_lower_z(a_from.m_lower_z) 88 ,m_upper_z(a_from.m_upper_z) 89 ,m_Sxw(a_from.m_Sxw) 90 ,m_Sx2w(a_from.m_Sx2w) 91 ,m_Syw(a_from.m_Syw) 92 ,m_Sy2w(a_from.m_Sy2w) 93 ,m_Szw(a_from.m_Szw) 94 ,m_Sz2w(a_from.m_Sz2w) 95 ,m_cnv_x_num(a_from.m_cnv_x_num) 96 ,m_cnv_x_min(a_from.m_cnv_x_min) 97 ,m_cnv_x_max(a_from.m_cnv_x_max) 98 ,m_cnv_y_num(a_from.m_cnv_y_num) 99 ,m_cnv_y_min(a_from.m_cnv_y_min) 100 ,m_cnv_y_max(a_from.m_cnv_y_max) 101 ,m_cnv_z_num(a_from.m_cnv_z_num) 102 ,m_cnv_z_min(a_from.m_cnv_z_min) 103 ,m_cnv_z_max(a_from.m_cnv_z_max) 104 ,m_histo(0) 105 { 106 if(a_from.m_histo) { 107 m_histo = new histo::h3d(*a_from.m_histo 108 } 109 } 110 111 c3d& operator=(const c3d& a_from) { 112 base_cloud::operator=(a_from); 113 if(&a_from==this) return *this; 114 m_xs = a_from.m_xs; 115 m_ys = a_from.m_ys; 116 m_zs = a_from.m_zs; 117 m_lower_x = a_from.m_lower_x; 118 m_upper_x = a_from.m_upper_x; 119 m_lower_y = a_from.m_lower_y; 120 m_upper_y = a_from.m_upper_y; 121 m_lower_z = a_from.m_lower_z; 122 m_upper_z = a_from.m_upper_z; 123 m_Sxw = a_from.m_Sxw; 124 m_Sx2w = a_from.m_Sx2w; 125 m_Syw = a_from.m_Syw; 126 m_Sy2w = a_from.m_Sy2w; 127 m_Szw = a_from.m_Szw; 128 m_Sz2w = a_from.m_Sz2w; 129 m_cnv_x_num = a_from.m_cnv_x_num; 130 m_cnv_x_min = a_from.m_cnv_x_min; 131 m_cnv_x_max = a_from.m_cnv_x_max; 132 m_cnv_y_num = a_from.m_cnv_y_num; 133 m_cnv_y_min = a_from.m_cnv_y_min; 134 m_cnv_y_max = a_from.m_cnv_y_max; 135 m_cnv_z_num = a_from.m_cnv_z_num; 136 m_cnv_z_min = a_from.m_cnv_z_min; 137 m_cnv_z_max = a_from.m_cnv_z_max; 138 delete m_histo; 139 m_histo = 0; 140 if(a_from.m_histo) { 141 m_histo = new histo::h3d(*a_from.m_histo 142 } 143 return *this; 144 } 145 146 protected: 147 void clear(); 148 protected: 149 std::vector<double> m_xs; 150 std::vector<double> m_ys; 151 std::vector<double> m_zs; 152 double m_lower_x; 153 double m_upper_x; 154 double m_lower_y; 155 double m_upper_y; 156 double m_lower_z; 157 double m_upper_z; 158 double m_Sxw; 159 double m_Sx2w; 160 double m_Syw; 161 double m_Sy2w; 162 double m_Szw; 163 double m_Sz2w; 164 // 165 unsigned int m_cnv_x_num; 166 double m_cnv_x_min; 167 double m_cnv_x_max; 168 unsigned int m_cnv_y_num; 169 double m_cnv_y_min; 170 double m_cnv_y_max; 171 unsigned int m_cnv_z_num; 172 double m_cnv_z_min; 173 double m_cnv_z_max; 174 histo::h3d* m_histo; 175 }; 176 177 }} 178 179 namespace tools { 180 namespace histo { 181 182 inline 183 c3d::c3d() 184 :base_cloud(UNLIMITED()) 185 ,m_lower_x(0) 186 ,m_upper_x(0) 187 ,m_lower_y(0) 188 ,m_upper_y(0) 189 ,m_lower_z(0) 190 ,m_upper_z(0) 191 ,m_Sxw(0) 192 ,m_Sx2w(0) 193 ,m_Syw(0) 194 ,m_Sy2w(0) 195 ,m_Szw(0) 196 ,m_Sz2w(0) 197 ,m_cnv_x_num(0) 198 ,m_cnv_x_min(0) 199 ,m_cnv_x_max(0) 200 ,m_cnv_y_num(0) 201 ,m_cnv_y_min(0) 202 ,m_cnv_y_max(0) 203 ,m_cnv_z_num(0) 204 ,m_cnv_z_min(0) 205 ,m_cnv_z_max(0) 206 ,m_histo(0) 207 {} 208 209 inline 210 c3d::c3d(const std::string& a_title,int aLimit 211 :base_cloud(aLimit) 212 ,m_lower_x(0) 213 ,m_upper_x(0) 214 ,m_lower_y(0) 215 ,m_upper_y(0) 216 ,m_lower_z(0) 217 ,m_upper_z(0) 218 ,m_Sxw(0) 219 ,m_Sx2w(0) 220 ,m_Syw(0) 221 ,m_Sy2w(0) 222 ,m_Szw(0) 223 ,m_Sz2w(0) 224 ,m_cnv_x_num(0) 225 ,m_cnv_x_min(0) 226 ,m_cnv_x_max(0) 227 ,m_cnv_y_num(0) 228 ,m_cnv_y_min(0) 229 ,m_cnv_y_max(0) 230 ,m_cnv_z_num(0) 231 ,m_cnv_z_min(0) 232 ,m_cnv_z_max(0) 233 ,m_histo(0) 234 { 235 set_title(a_title); 236 } 237 238 inline 239 bool c3d::is_converted() const {return m_histo 240 241 inline 242 void c3d::clear(){ 243 m_lower_x = 0; 244 m_upper_x = 0; 245 m_lower_y = 0; 246 m_upper_y = 0; 247 m_lower_z = 0; 248 m_upper_z = 0; 249 m_Sw = 0; 250 m_Sxw = 0; 251 m_Sx2w = 0; 252 m_Syw = 0; 253 m_Sy2w = 0; 254 m_Szw = 0; 255 m_Sz2w = 0; 256 m_xs.clear(); 257 m_ys.clear(); 258 m_zs.clear(); 259 m_ws.clear(); 260 } 261 262 inline 263 bool c3d::convert( 264 unsigned int a_bins_x,double a_lower_edge_x,d 265 ,unsigned int a_bins_y,double a_lower_edge_y,d 266 ,unsigned int a_bins_z,double a_lower_edge_z,d 267 ) { 268 if(m_histo) return true; // Done. 269 m_histo = new histo::h3d(base_cloud::title() 270 a_bins_x,a_lower_ed 271 a_bins_y,a_lower_ed 272 a_bins_z,a_lower_ed 273 if(!m_histo) return false; 274 bool status = fill_histogram(*m_histo); 275 clear(); 276 return status; 277 } 278 279 inline 280 bool c3d::convert_to_histogram(){ 281 if( (m_cnv_x_num<=0) || (m_cnv_x_max<=m_cnv_ 282 (m_cnv_y_num<=0) || (m_cnv_y_max<=m_cnv_ 283 (m_cnv_z_num<=0) || (m_cnv_z_max<=m_cnv_ 284 double dx = 0.01 * (upper_edge_x() - lower 285 double dy = 0.01 * (upper_edge_y() - lower 286 double dz = 0.01 * (upper_edge_z() - lower 287 return convert(BINS(),lower_edge_x(),upper 288 BINS(),lower_edge_y(),upper 289 BINS(),lower_edge_z(),upper 290 } else { 291 return convert(m_cnv_x_num,m_cnv_x_min,m_c 292 m_cnv_y_num,m_cnv_y_min,m_c 293 m_cnv_z_num,m_cnv_z_min,m_c 294 } 295 } 296 297 inline 298 bool c3d::set_title(const std::string& a_title 299 m_title = a_title; 300 if(m_histo) m_histo->set_title(a_title); 301 return true; 302 } 303 304 inline 305 bool c3d::scale(double a_scale) { 306 if(m_histo) { 307 return m_histo->scale(a_scale); 308 } else { 309 size_t number = m_ws.size(); 310 for(size_t index=0;index<number;index++) m 311 m_Sw *= a_scale; 312 m_Sxw *= a_scale; 313 m_Sx2w *= a_scale; 314 m_Syw *= a_scale; 315 m_Sy2w *= a_scale; 316 m_Szw *= a_scale; 317 m_Sz2w *= a_scale; 318 return true; 319 } 320 } 321 322 inline 323 bool c3d::set_conversion_parameters( 324 unsigned int aCnvXnumber,double aCnvXmin,doub 325 ,unsigned int aCnvYnumber,double aCnvYmin,doub 326 ,unsigned int aCnvZnumber,double aCnvZmin,doub 327 ){ 328 m_cnv_x_num = aCnvXnumber; 329 m_cnv_x_min = aCnvXmin; 330 m_cnv_x_max = aCnvXmax; 331 m_cnv_y_num = aCnvYnumber; 332 m_cnv_y_min = aCnvYmin; 333 m_cnv_y_max = aCnvYmax; 334 m_cnv_z_num = aCnvZnumber; 335 m_cnv_z_min = aCnvZmin; 336 m_cnv_z_max = aCnvZmax; 337 return true; 338 } 339 340 341 inline 342 const h3d& c3d::histogram() const { 343 if(!m_histo) const_cast<c3d&>(*this).convert 344 return *m_histo; 345 } 346 347 inline 348 bool c3d::reset() { 349 clear(); 350 delete m_histo; 351 m_histo = 0; 352 return true; 353 } 354 355 inline 356 bool c3d::fill(double aX,double aY,double aZ,d 357 if(!m_histo && (m_limit!=UNLIMITED()) && ((i 358 convert_to_histogram(); 359 } 360 361 if(m_histo) { 362 return m_histo->fill(aX,aY,aZ,aW); 363 } else { 364 if(m_xs.size()) { 365 m_lower_x = mn<double>(aX,m_lower_x); 366 m_upper_x = mx<double>(aX,m_upper_x); 367 } else { 368 m_lower_x = aX; 369 m_upper_x = aX; 370 } 371 if(m_ys.size()) { 372 m_lower_y = mn<double>(aY,m_lower_y); 373 m_upper_y = mx<double>(aY,m_upper_y); 374 } else { 375 m_lower_y = aY; 376 m_upper_y = aY; 377 } 378 if(m_zs.size()) { 379 m_lower_z = mn<double>(aZ,m_lower_z); 380 m_upper_z = mx<double>(aZ,m_upper_z); 381 } else { 382 m_lower_z = aZ; 383 m_upper_z = aZ; 384 } 385 m_xs.push_back(aX); 386 m_ys.push_back(aY); 387 m_zs.push_back(aZ); 388 m_ws.push_back(aW); 389 m_Sw += aW; 390 double xw = aX * aW; 391 m_Sxw += xw; 392 m_Sx2w += aX * xw; 393 double yw = aY * aW; 394 m_Syw += yw; 395 m_Sy2w += aY * yw; 396 double zw = aZ * aW; 397 m_Szw += zw; 398 m_Sz2w += aZ * zw; 399 return true; 400 } 401 } 402 403 inline 404 bool c3d::convert( 405 const std::vector<double>& a_edges_x 406 ,const std::vector<double>& a_edges_y 407 ,const std::vector<double>& a_edges_z 408 ) { 409 if(m_histo) return true; 410 m_histo = new histo::h3d(base_cloud::title() 411 a_edges_x,a_edges_y 412 if(!m_histo) return false; 413 bool status = fill_histogram(*m_histo); 414 clear(); 415 return status; 416 } 417 418 inline 419 double c3d::sum_of_weights() const { 420 return (m_histo ? m_histo->sum_bin_heights() 421 } 422 423 inline 424 unsigned int c3d::entries() const { 425 return m_histo ? m_histo->all_entries() : (u 426 } 427 428 inline 429 double c3d::lower_edge_x() const { 430 return m_histo ? m_histo->axis_x().lower_edg 431 } 432 inline 433 double c3d::lower_edge_y() const { 434 return m_histo ? m_histo->axis_y().lower_edg 435 } 436 inline 437 double c3d::lower_edge_z() const { 438 return m_histo ? m_histo->axis_z().lower_edg 439 } 440 inline 441 double c3d::upper_edge_x() const { 442 return m_histo ? m_histo->axis_x().upper_edg 443 } 444 inline 445 double c3d::upper_edge_y() const { 446 return m_histo ? m_histo->axis_y().upper_edg 447 } 448 inline 449 double c3d::upper_edge_z() const { 450 return m_histo ? m_histo->axis_z().upper_edg 451 } 452 inline 453 double c3d::value_x(unsigned int a_index) cons 454 return m_histo ? 0 : m_xs[a_index]; 455 } 456 inline 457 double c3d::value_y(unsigned int a_index) cons 458 return m_histo ? 0 : m_ys[a_index]; 459 } 460 inline 461 double c3d::value_z(unsigned int a_index) cons 462 return m_histo ? 0 : m_zs[a_index]; 463 } 464 inline 465 double c3d::weight(unsigned int a_index) const 466 return m_histo ? 0 : m_ws[a_index]; 467 } 468 inline 469 double c3d::mean_x() const { 470 return m_histo ? m_histo->mean_x() : (m_Sw?m 471 } 472 inline 473 double c3d::mean_y() const { 474 return m_histo ? m_histo->mean_y() : (m_Sw?m 475 } 476 inline 477 double c3d::mean_z() const { 478 return m_histo ? m_histo->mean_z() : (m_Sw?m 479 } 480 inline 481 double c3d::rms_x() const { 482 double rms = 0; //FIXME nan. 483 if(m_histo) { 484 rms = m_histo->rms_x(); 485 } else { 486 if(m_Sw==0) { 487 } else { 488 double mean = m_Sxw / m_Sw; 489 rms = ::sqrt(::fabs( (m_Sx2w / m_Sw) - m 490 } 491 } 492 return rms; 493 } 494 inline 495 double c3d::rms_y() const { 496 double rms = 0; //FIXME nan. 497 if(m_histo) { 498 rms = m_histo->rms_y(); 499 } else { 500 if(m_Sw==0) { 501 } else { 502 double mean = m_Syw / m_Sw; 503 rms = ::sqrt(::fabs( (m_Sy2w / m_Sw) - m 504 } 505 } 506 return rms; 507 } 508 inline 509 double c3d::rms_z() const { 510 double rms = 0; //FIXME nan. 511 if(m_histo) { 512 rms = m_histo->rms_z(); 513 } else { 514 if(m_Sw==0) { 515 } else { 516 double mean = m_Szw / m_Sw; 517 rms = ::sqrt(::fabs( (m_Sz2w / m_Sw) - m 518 } 519 } 520 return rms; 521 } 522 523 }} 524 525 #endif