Warning, /include/Geant4/tools/histo/c3d 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_c3d
0005 #define tools_histo_c3d
0006
0007 #include "base_cloud"
0008
0009 #include "../mnmx"
0010
0011 #include "h3d"
0012
0013 namespace tools {
0014 namespace histo {
0015
0016 class c3d : public base_cloud {
0017 public:
0018 static const std::string& s_class() {
0019 static const std::string s_v("tools::histo::c3d");
0020 return s_v;
0021 }
0022 public:
0023 bool set_title(const std::string&);
0024 unsigned int dimension() const {return 3;}
0025 bool reset();
0026 unsigned int entries() const;
0027 public:
0028 double sum_of_weights() const;
0029 bool convert_to_histogram();
0030 bool is_converted() const;
0031 bool scale(double);
0032 public:
0033 bool fill(double,double,double,double = 1);
0034 double lower_edge_x() const;
0035 double upper_edge_x() const;
0036 double lower_edge_y() const;
0037 double upper_edge_y() const;
0038 double lower_edge_z() const;
0039 double upper_edge_z() const;
0040 double value_x(unsigned int) const;
0041 double value_y(unsigned int) const;
0042 double value_z(unsigned int) const;
0043 double weight(unsigned int) const;
0044 double mean_x() const;
0045 double mean_y() const;
0046 double mean_z() const;
0047 double rms_x() const;
0048 double rms_y() const;
0049 double rms_z() const;
0050 bool convert(unsigned int,double,double,
0051 unsigned int,double,double,
0052 unsigned int,double,double);
0053 bool convert(const std::vector<double>&,
0054 const std::vector<double>&,
0055 const std::vector<double>&);
0056 const histo::h3d& histogram() const;
0057 bool fill_histogram(histo::h3d& a_histo) const {
0058 size_t number = m_xs.size();
0059 for(size_t index=0;index<number;index++) {
0060 if(!a_histo.fill(m_xs[index],m_ys[index],m_zs[index],m_ws[index])) return false;
0061 }
0062 return true;
0063 }
0064 bool set_conversion_parameters(unsigned int,double,double,
0065 unsigned int,double,double,
0066 unsigned int,double,double);
0067
0068 bool set_histogram(h3d* a_histo){ //we take ownership of a_histo.
0069 reset();
0070 m_histo = a_histo;
0071 return true;
0072 }
0073 public:
0074 c3d();
0075 c3d(const std::string&,int=base_cloud::UNLIMITED());
0076 virtual ~c3d(){delete m_histo;}
0077 public:
0078 c3d(const c3d& a_from)
0079 :base_cloud(a_from)
0080 ,m_xs(a_from.m_xs)
0081 ,m_ys(a_from.m_ys)
0082 ,m_zs(a_from.m_zs)
0083 ,m_lower_x(a_from.m_lower_x)
0084 ,m_upper_x(a_from.m_upper_x)
0085 ,m_lower_y(a_from.m_lower_y)
0086 ,m_upper_y(a_from.m_upper_y)
0087 ,m_lower_z(a_from.m_lower_z)
0088 ,m_upper_z(a_from.m_upper_z)
0089 ,m_Sxw(a_from.m_Sxw)
0090 ,m_Sx2w(a_from.m_Sx2w)
0091 ,m_Syw(a_from.m_Syw)
0092 ,m_Sy2w(a_from.m_Sy2w)
0093 ,m_Szw(a_from.m_Szw)
0094 ,m_Sz2w(a_from.m_Sz2w)
0095 ,m_cnv_x_num(a_from.m_cnv_x_num)
0096 ,m_cnv_x_min(a_from.m_cnv_x_min)
0097 ,m_cnv_x_max(a_from.m_cnv_x_max)
0098 ,m_cnv_y_num(a_from.m_cnv_y_num)
0099 ,m_cnv_y_min(a_from.m_cnv_y_min)
0100 ,m_cnv_y_max(a_from.m_cnv_y_max)
0101 ,m_cnv_z_num(a_from.m_cnv_z_num)
0102 ,m_cnv_z_min(a_from.m_cnv_z_min)
0103 ,m_cnv_z_max(a_from.m_cnv_z_max)
0104 ,m_histo(0)
0105 {
0106 if(a_from.m_histo) {
0107 m_histo = new histo::h3d(*a_from.m_histo);
0108 }
0109 }
0110
0111 c3d& operator=(const c3d& a_from) {
0112 base_cloud::operator=(a_from);
0113 if(&a_from==this) return *this;
0114 m_xs = a_from.m_xs;
0115 m_ys = a_from.m_ys;
0116 m_zs = a_from.m_zs;
0117 m_lower_x = a_from.m_lower_x;
0118 m_upper_x = a_from.m_upper_x;
0119 m_lower_y = a_from.m_lower_y;
0120 m_upper_y = a_from.m_upper_y;
0121 m_lower_z = a_from.m_lower_z;
0122 m_upper_z = a_from.m_upper_z;
0123 m_Sxw = a_from.m_Sxw;
0124 m_Sx2w = a_from.m_Sx2w;
0125 m_Syw = a_from.m_Syw;
0126 m_Sy2w = a_from.m_Sy2w;
0127 m_Szw = a_from.m_Szw;
0128 m_Sz2w = a_from.m_Sz2w;
0129 m_cnv_x_num = a_from.m_cnv_x_num;
0130 m_cnv_x_min = a_from.m_cnv_x_min;
0131 m_cnv_x_max = a_from.m_cnv_x_max;
0132 m_cnv_y_num = a_from.m_cnv_y_num;
0133 m_cnv_y_min = a_from.m_cnv_y_min;
0134 m_cnv_y_max = a_from.m_cnv_y_max;
0135 m_cnv_z_num = a_from.m_cnv_z_num;
0136 m_cnv_z_min = a_from.m_cnv_z_min;
0137 m_cnv_z_max = a_from.m_cnv_z_max;
0138 delete m_histo;
0139 m_histo = 0;
0140 if(a_from.m_histo) {
0141 m_histo = new histo::h3d(*a_from.m_histo);
0142 }
0143 return *this;
0144 }
0145
0146 protected:
0147 void clear();
0148 protected:
0149 std::vector<double> m_xs;
0150 std::vector<double> m_ys;
0151 std::vector<double> m_zs;
0152 double m_lower_x;
0153 double m_upper_x;
0154 double m_lower_y;
0155 double m_upper_y;
0156 double m_lower_z;
0157 double m_upper_z;
0158 double m_Sxw;
0159 double m_Sx2w;
0160 double m_Syw;
0161 double m_Sy2w;
0162 double m_Szw;
0163 double m_Sz2w;
0164 //
0165 unsigned int m_cnv_x_num;
0166 double m_cnv_x_min;
0167 double m_cnv_x_max;
0168 unsigned int m_cnv_y_num;
0169 double m_cnv_y_min;
0170 double m_cnv_y_max;
0171 unsigned int m_cnv_z_num;
0172 double m_cnv_z_min;
0173 double m_cnv_z_max;
0174 histo::h3d* m_histo;
0175 };
0176
0177 }}
0178
0179 namespace tools {
0180 namespace histo {
0181
0182 inline
0183 c3d::c3d()
0184 :base_cloud(UNLIMITED())
0185 ,m_lower_x(0)
0186 ,m_upper_x(0)
0187 ,m_lower_y(0)
0188 ,m_upper_y(0)
0189 ,m_lower_z(0)
0190 ,m_upper_z(0)
0191 ,m_Sxw(0)
0192 ,m_Sx2w(0)
0193 ,m_Syw(0)
0194 ,m_Sy2w(0)
0195 ,m_Szw(0)
0196 ,m_Sz2w(0)
0197 ,m_cnv_x_num(0)
0198 ,m_cnv_x_min(0)
0199 ,m_cnv_x_max(0)
0200 ,m_cnv_y_num(0)
0201 ,m_cnv_y_min(0)
0202 ,m_cnv_y_max(0)
0203 ,m_cnv_z_num(0)
0204 ,m_cnv_z_min(0)
0205 ,m_cnv_z_max(0)
0206 ,m_histo(0)
0207 {}
0208
0209 inline
0210 c3d::c3d(const std::string& a_title,int aLimit)
0211 :base_cloud(aLimit)
0212 ,m_lower_x(0)
0213 ,m_upper_x(0)
0214 ,m_lower_y(0)
0215 ,m_upper_y(0)
0216 ,m_lower_z(0)
0217 ,m_upper_z(0)
0218 ,m_Sxw(0)
0219 ,m_Sx2w(0)
0220 ,m_Syw(0)
0221 ,m_Sy2w(0)
0222 ,m_Szw(0)
0223 ,m_Sz2w(0)
0224 ,m_cnv_x_num(0)
0225 ,m_cnv_x_min(0)
0226 ,m_cnv_x_max(0)
0227 ,m_cnv_y_num(0)
0228 ,m_cnv_y_min(0)
0229 ,m_cnv_y_max(0)
0230 ,m_cnv_z_num(0)
0231 ,m_cnv_z_min(0)
0232 ,m_cnv_z_max(0)
0233 ,m_histo(0)
0234 {
0235 set_title(a_title);
0236 }
0237
0238 inline
0239 bool c3d::is_converted() const {return m_histo ? true : false;}
0240
0241 inline
0242 void c3d::clear(){
0243 m_lower_x = 0;
0244 m_upper_x = 0;
0245 m_lower_y = 0;
0246 m_upper_y = 0;
0247 m_lower_z = 0;
0248 m_upper_z = 0;
0249 m_Sw = 0;
0250 m_Sxw = 0;
0251 m_Sx2w = 0;
0252 m_Syw = 0;
0253 m_Sy2w = 0;
0254 m_Szw = 0;
0255 m_Sz2w = 0;
0256 m_xs.clear();
0257 m_ys.clear();
0258 m_zs.clear();
0259 m_ws.clear();
0260 }
0261
0262 inline
0263 bool c3d::convert(
0264 unsigned int a_bins_x,double a_lower_edge_x,double a_upper_edge_x
0265 ,unsigned int a_bins_y,double a_lower_edge_y,double a_upper_edge_y
0266 ,unsigned int a_bins_z,double a_lower_edge_z,double a_upper_edge_z
0267 ) {
0268 if(m_histo) return true; // Done.
0269 m_histo = new histo::h3d(base_cloud::title(),
0270 a_bins_x,a_lower_edge_x,a_upper_edge_x,
0271 a_bins_y,a_lower_edge_y,a_upper_edge_y,
0272 a_bins_z,a_lower_edge_z,a_upper_edge_z);
0273 if(!m_histo) return false;
0274 bool status = fill_histogram(*m_histo);
0275 clear();
0276 return status;
0277 }
0278
0279 inline
0280 bool c3d::convert_to_histogram(){
0281 if( (m_cnv_x_num<=0) || (m_cnv_x_max<=m_cnv_x_min) ||
0282 (m_cnv_y_num<=0) || (m_cnv_y_max<=m_cnv_y_min) ||
0283 (m_cnv_z_num<=0) || (m_cnv_z_max<=m_cnv_z_min) ) {
0284 double dx = 0.01 * (upper_edge_x() - lower_edge_x())/BINS();
0285 double dy = 0.01 * (upper_edge_y() - lower_edge_y())/BINS();
0286 double dz = 0.01 * (upper_edge_z() - lower_edge_z())/BINS();
0287 return convert(BINS(),lower_edge_x(),upper_edge_x()+dx,
0288 BINS(),lower_edge_y(),upper_edge_y()+dy,
0289 BINS(),lower_edge_z(),upper_edge_z()+dz);
0290 } else {
0291 return convert(m_cnv_x_num,m_cnv_x_min,m_cnv_x_max,
0292 m_cnv_y_num,m_cnv_y_min,m_cnv_y_max,
0293 m_cnv_z_num,m_cnv_z_min,m_cnv_z_max);
0294 }
0295 }
0296
0297 inline
0298 bool c3d::set_title(const std::string& a_title){
0299 m_title = a_title;
0300 if(m_histo) m_histo->set_title(a_title);
0301 return true;
0302 }
0303
0304 inline
0305 bool c3d::scale(double a_scale) {
0306 if(m_histo) {
0307 return m_histo->scale(a_scale);
0308 } else {
0309 size_t number = m_ws.size();
0310 for(size_t index=0;index<number;index++) m_ws[index] *= a_scale;
0311 m_Sw *= a_scale;
0312 m_Sxw *= a_scale;
0313 m_Sx2w *= a_scale;
0314 m_Syw *= a_scale;
0315 m_Sy2w *= a_scale;
0316 m_Szw *= a_scale;
0317 m_Sz2w *= a_scale;
0318 return true;
0319 }
0320 }
0321
0322 inline
0323 bool c3d::set_conversion_parameters(
0324 unsigned int aCnvXnumber,double aCnvXmin,double aCnvXmax
0325 ,unsigned int aCnvYnumber,double aCnvYmin,double aCnvYmax
0326 ,unsigned int aCnvZnumber,double aCnvZmin,double aCnvZmax
0327 ){
0328 m_cnv_x_num = aCnvXnumber;
0329 m_cnv_x_min = aCnvXmin;
0330 m_cnv_x_max = aCnvXmax;
0331 m_cnv_y_num = aCnvYnumber;
0332 m_cnv_y_min = aCnvYmin;
0333 m_cnv_y_max = aCnvYmax;
0334 m_cnv_z_num = aCnvZnumber;
0335 m_cnv_z_min = aCnvZmin;
0336 m_cnv_z_max = aCnvZmax;
0337 return true;
0338 }
0339
0340
0341 inline
0342 const h3d& c3d::histogram() const {
0343 if(!m_histo) const_cast<c3d&>(*this).convert_to_histogram();
0344 return *m_histo;
0345 }
0346
0347 inline
0348 bool c3d::reset() {
0349 clear();
0350 delete m_histo;
0351 m_histo = 0;
0352 return true;
0353 }
0354
0355 inline
0356 bool c3d::fill(double aX,double aY,double aZ,double aW){
0357 if(!m_histo && (m_limit!=UNLIMITED()) && ((int)m_xs.size()>=m_limit)){
0358 convert_to_histogram();
0359 }
0360
0361 if(m_histo) {
0362 return m_histo->fill(aX,aY,aZ,aW);
0363 } else {
0364 if(m_xs.size()) {
0365 m_lower_x = mn<double>(aX,m_lower_x);
0366 m_upper_x = mx<double>(aX,m_upper_x);
0367 } else {
0368 m_lower_x = aX;
0369 m_upper_x = aX;
0370 }
0371 if(m_ys.size()) {
0372 m_lower_y = mn<double>(aY,m_lower_y);
0373 m_upper_y = mx<double>(aY,m_upper_y);
0374 } else {
0375 m_lower_y = aY;
0376 m_upper_y = aY;
0377 }
0378 if(m_zs.size()) {
0379 m_lower_z = mn<double>(aZ,m_lower_z);
0380 m_upper_z = mx<double>(aZ,m_upper_z);
0381 } else {
0382 m_lower_z = aZ;
0383 m_upper_z = aZ;
0384 }
0385 m_xs.push_back(aX);
0386 m_ys.push_back(aY);
0387 m_zs.push_back(aZ);
0388 m_ws.push_back(aW);
0389 m_Sw += aW;
0390 double xw = aX * aW;
0391 m_Sxw += xw;
0392 m_Sx2w += aX * xw;
0393 double yw = aY * aW;
0394 m_Syw += yw;
0395 m_Sy2w += aY * yw;
0396 double zw = aZ * aW;
0397 m_Szw += zw;
0398 m_Sz2w += aZ * zw;
0399 return true;
0400 }
0401 }
0402
0403 inline
0404 bool c3d::convert(
0405 const std::vector<double>& a_edges_x
0406 ,const std::vector<double>& a_edges_y
0407 ,const std::vector<double>& a_edges_z
0408 ) {
0409 if(m_histo) return true;
0410 m_histo = new histo::h3d(base_cloud::title(),
0411 a_edges_x,a_edges_y,a_edges_z);
0412 if(!m_histo) return false;
0413 bool status = fill_histogram(*m_histo);
0414 clear();
0415 return status;
0416 }
0417
0418 inline
0419 double c3d::sum_of_weights() const {
0420 return (m_histo ? m_histo->sum_bin_heights() : m_Sw);
0421 }
0422
0423 inline
0424 unsigned int c3d::entries() const {
0425 return m_histo ? m_histo->all_entries() : (unsigned int)m_ws.size();
0426 }
0427
0428 inline
0429 double c3d::lower_edge_x() const {
0430 return m_histo ? m_histo->axis_x().lower_edge() : m_lower_x;
0431 }
0432 inline
0433 double c3d::lower_edge_y() const {
0434 return m_histo ? m_histo->axis_y().lower_edge() : m_lower_y;
0435 }
0436 inline
0437 double c3d::lower_edge_z() const {
0438 return m_histo ? m_histo->axis_z().lower_edge() : m_lower_z;
0439 }
0440 inline
0441 double c3d::upper_edge_x() const {
0442 return m_histo ? m_histo->axis_x().upper_edge() : m_upper_x;
0443 }
0444 inline
0445 double c3d::upper_edge_y() const {
0446 return m_histo ? m_histo->axis_y().upper_edge() : m_upper_y;
0447 }
0448 inline
0449 double c3d::upper_edge_z() const {
0450 return m_histo ? m_histo->axis_z().upper_edge() : m_upper_z;
0451 }
0452 inline
0453 double c3d::value_x(unsigned int a_index) const {
0454 return m_histo ? 0 : m_xs[a_index];
0455 }
0456 inline
0457 double c3d::value_y(unsigned int a_index) const {
0458 return m_histo ? 0 : m_ys[a_index];
0459 }
0460 inline
0461 double c3d::value_z(unsigned int a_index) const {
0462 return m_histo ? 0 : m_zs[a_index];
0463 }
0464 inline
0465 double c3d::weight(unsigned int a_index) const {
0466 return m_histo ? 0 : m_ws[a_index];
0467 }
0468 inline
0469 double c3d::mean_x() const {
0470 return m_histo ? m_histo->mean_x() : (m_Sw?m_Sxw/m_Sw:0);
0471 }
0472 inline
0473 double c3d::mean_y() const {
0474 return m_histo ? m_histo->mean_y() : (m_Sw?m_Syw/m_Sw:0);
0475 }
0476 inline
0477 double c3d::mean_z() const {
0478 return m_histo ? m_histo->mean_z() : (m_Sw?m_Szw/m_Sw:0);
0479 }
0480 inline
0481 double c3d::rms_x() const {
0482 double rms = 0; //FIXME nan.
0483 if(m_histo) {
0484 rms = m_histo->rms_x();
0485 } else {
0486 if(m_Sw==0) {
0487 } else {
0488 double mean = m_Sxw / m_Sw;
0489 rms = ::sqrt(::fabs( (m_Sx2w / m_Sw) - mean * mean));
0490 }
0491 }
0492 return rms;
0493 }
0494 inline
0495 double c3d::rms_y() const {
0496 double rms = 0; //FIXME nan.
0497 if(m_histo) {
0498 rms = m_histo->rms_y();
0499 } else {
0500 if(m_Sw==0) {
0501 } else {
0502 double mean = m_Syw / m_Sw;
0503 rms = ::sqrt(::fabs( (m_Sy2w / m_Sw) - mean * mean));
0504 }
0505 }
0506 return rms;
0507 }
0508 inline
0509 double c3d::rms_z() const {
0510 double rms = 0; //FIXME nan.
0511 if(m_histo) {
0512 rms = m_histo->rms_z();
0513 } else {
0514 if(m_Sw==0) {
0515 } else {
0516 double mean = m_Szw / m_Sw;
0517 rms = ::sqrt(::fabs( (m_Sz2w / m_Sw) - mean * mean));
0518 }
0519 }
0520 return rms;
0521 }
0522
0523 }}
0524
0525 #endif
