Warning, /include/Geant4/tools/histo/c2d 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_c2d
0005 #define tools_histo_c2d
0006
0007 #include "base_cloud"
0008
0009 #include "../mnmx"
0010
0011 #include "h2d"
0012
0013 namespace tools {
0014 namespace histo {
0015
0016 class c2d : public base_cloud {
0017 public:
0018 static const std::string& s_class() {
0019 static const std::string s_v("tools::histo::c2d");
0020 return s_v;
0021 }
0022 public:
0023 bool set_title(const std::string&);
0024 unsigned int dimension() const {return 2;}
0025 bool reset();
0026 unsigned int entries() const { return m_histo ? m_histo->all_entries() : (unsigned int)m_ws.size();}
0027 public:
0028 double sum_of_weights() const { return (m_histo ? m_histo->sum_bin_heights() : m_Sw);}
0029 bool convert_to_histogram();
0030 bool is_converted() const {return m_histo ? true : false;}
0031 bool scale(double);
0032 public:
0033 bool fill(double,double,double = 1);
0034 double lower_edge_x() const {return m_histo ? m_histo->axis_x().lower_edge() : m_lower_x;}
0035 double lower_edge_y() const { return m_histo ? m_histo->axis_y().lower_edge() : m_lower_y;}
0036 double upper_edge_x() const { return m_histo ? m_histo->axis_x().upper_edge() : m_upper_x;}
0037 double upper_edge_y() const { return m_histo ? m_histo->axis_y().upper_edge() : m_upper_y;}
0038 double value_x(unsigned int a_index) const { return m_histo ? 0 : m_xs[a_index];}
0039 double value_y(unsigned int a_index) const { return m_histo ? 0 : m_ys[a_index];}
0040 double weight(unsigned int a_index) const { return m_histo ? 0 : m_ws[a_index];}
0041 double mean_x() const {return m_histo ? m_histo->mean_x() : (m_Sw?m_Sxw/m_Sw:0);}
0042 double mean_y() const {return m_histo ? m_histo->mean_y() : (m_Sw?m_Syw/m_Sw:0);}
0043 double rms_x() const;
0044 double rms_y() const;
0045 bool convert(unsigned int,double,double,unsigned int,double,double);
0046 bool convert(const std::vector<double>&,const std::vector<double>&);
0047 const histo::h2d& histogram() const;
0048 template <class HISTO>
0049 bool fill_histogram(HISTO& a_histo) const {
0050 size_t number = m_xs.size();
0051 for(size_t index=0;index<number;index++) {
0052 if(!a_histo.fill(m_xs[index],m_ys[index],m_ws[index])) return false;
0053 }
0054 return true;
0055 }
0056 bool set_conversion_parameters(unsigned int,double,double,unsigned int,double,double);
0057
0058 bool set_histogram(h2d* a_histo){ //we take ownership of a_histo.
0059 reset();
0060 m_histo = a_histo;
0061 return true;
0062 }
0063 public:
0064 c2d();
0065 c2d(const std::string&,int=base_cloud::UNLIMITED());
0066 virtual ~c2d(){delete m_histo;}
0067 public:
0068 c2d(const c2d& a_from)
0069 :base_cloud(a_from)
0070 ,m_xs(a_from.m_xs)
0071 ,m_ys(a_from.m_ys)
0072 ,m_lower_x(a_from.m_lower_x)
0073 ,m_upper_x(a_from.m_upper_x)
0074 ,m_lower_y(a_from.m_lower_y)
0075 ,m_upper_y(a_from.m_upper_y)
0076 ,m_Sxw(a_from.m_Sxw)
0077 ,m_Sx2w(a_from.m_Sx2w)
0078 ,m_Syw(a_from.m_Syw)
0079 ,m_Sy2w(a_from.m_Sy2w)
0080 ,m_cnv_x_num(a_from.m_cnv_x_num)
0081 ,m_cnv_x_min(a_from.m_cnv_x_min)
0082 ,m_cnv_x_max(a_from.m_cnv_x_max)
0083 ,m_cnv_y_num(a_from.m_cnv_y_num)
0084 ,m_cnv_y_min(a_from.m_cnv_y_min)
0085 ,m_cnv_y_max(a_from.m_cnv_y_max)
0086 ,m_histo(0)
0087 {
0088 if(a_from.m_histo) {
0089 m_histo = new histo::h2d(*a_from.m_histo);
0090 }
0091 }
0092
0093 c2d& operator=(const c2d& a_from) {
0094 base_cloud::operator=(a_from);
0095 if(&a_from==this) return *this;
0096 m_xs = a_from.m_xs;
0097 m_ys = a_from.m_ys;
0098 m_lower_x = a_from.m_lower_x;
0099 m_upper_x = a_from.m_upper_x;
0100 m_lower_y = a_from.m_lower_y;
0101 m_upper_y = a_from.m_upper_y;
0102 m_Sxw = a_from.m_Sxw;
0103 m_Sx2w = a_from.m_Sx2w;
0104 m_Syw = a_from.m_Syw;
0105 m_Sy2w = a_from.m_Sy2w;
0106 m_cnv_x_num = a_from.m_cnv_x_num;
0107 m_cnv_x_min = a_from.m_cnv_x_min;
0108 m_cnv_x_max = a_from.m_cnv_x_max;
0109 m_cnv_y_num = a_from.m_cnv_y_num;
0110 m_cnv_y_min = a_from.m_cnv_y_min;
0111 m_cnv_y_max = a_from.m_cnv_y_max;
0112 delete m_histo;
0113 m_histo = 0;
0114 if(a_from.m_histo) {
0115 m_histo = new histo::h2d(*a_from.m_histo);
0116 }
0117 return *this;
0118 }
0119 public: //AIDA API
0120 double lowerEdgeX() const {return lower_edge_x();}
0121 double lowerEdgeY() const {return lower_edge_y();}
0122 double upperEdgeX() const {return upper_edge_x();}
0123 double upperEdgeY() const {return upper_edge_y();}
0124 template <class HISTO>
0125 bool fillHistogram(HISTO& a_histo) const {return fill_histogram<HISTO>(a_histo);}
0126 protected:
0127 void clear();
0128 protected:
0129 std::vector<double> m_xs;
0130 std::vector<double> m_ys;
0131 double m_lower_x;
0132 double m_upper_x;
0133 double m_lower_y;
0134 double m_upper_y;
0135 double m_Sxw;
0136 double m_Sx2w;
0137 double m_Syw;
0138 double m_Sy2w;
0139 //
0140 unsigned int m_cnv_x_num;
0141 double m_cnv_x_min;
0142 double m_cnv_x_max;
0143 unsigned int m_cnv_y_num;
0144 double m_cnv_y_min;
0145 double m_cnv_y_max;
0146 histo::h2d* m_histo;
0147 };
0148
0149 }}
0150
0151 namespace tools {
0152 namespace histo {
0153
0154 inline
0155 c2d::c2d()
0156 :base_cloud(UNLIMITED())
0157 ,m_lower_x(0)
0158 ,m_upper_x(0)
0159 ,m_lower_y(0)
0160 ,m_upper_y(0)
0161 ,m_Sxw(0)
0162 ,m_Sx2w(0)
0163 ,m_Syw(0)
0164 ,m_Sy2w(0)
0165 ,m_cnv_x_num(0)
0166 ,m_cnv_x_min(0)
0167 ,m_cnv_x_max(0)
0168 ,m_cnv_y_num(0)
0169 ,m_cnv_y_min(0)
0170 ,m_cnv_y_max(0)
0171 ,m_histo(0)
0172 {}
0173
0174 inline
0175 c2d::c2d(const std::string& a_title,int aLimit)
0176 :base_cloud(aLimit)
0177 ,m_lower_x(0)
0178 ,m_upper_x(0)
0179 ,m_lower_y(0)
0180 ,m_upper_y(0)
0181 ,m_Sxw(0)
0182 ,m_Sx2w(0)
0183 ,m_Syw(0)
0184 ,m_Sy2w(0)
0185 ,m_cnv_x_num(0)
0186 ,m_cnv_x_min(0)
0187 ,m_cnv_x_max(0)
0188 ,m_cnv_y_num(0)
0189 ,m_cnv_y_min(0)
0190 ,m_cnv_y_max(0)
0191 ,m_histo(0)
0192 {
0193 set_title(a_title);
0194 }
0195
0196 inline
0197 void c2d::clear(){
0198 m_lower_x = 0;
0199 m_upper_x = 0;
0200 m_lower_y = 0;
0201 m_upper_y = 0;
0202 m_Sw = 0;
0203 m_Sxw = 0;
0204 m_Sx2w = 0;
0205 m_Syw = 0;
0206 m_Sy2w = 0;
0207 m_xs.clear();
0208 m_ys.clear();
0209 m_ws.clear();
0210 }
0211
0212 inline
0213 bool c2d::convert(
0214 unsigned int a_bins_x,double a_lower_edge_x,double a_upper_edge_x
0215 ,unsigned int a_bins_y,double a_lower_edge_y,double a_upper_edge_y
0216 ) {
0217 if(m_histo) return true; // Done.
0218 m_histo = new histo::h2d(base_cloud::title(),
0219 a_bins_x,a_lower_edge_x,a_upper_edge_x,
0220 a_bins_y,a_lower_edge_y,a_upper_edge_y);
0221 if(!m_histo) return false;
0222 bool status = fill_histogram(*m_histo);
0223 clear();
0224 return status;
0225 }
0226
0227 inline
0228 bool c2d::convert_to_histogram(){
0229 if( (m_cnv_x_num<=0) || (m_cnv_x_max<=m_cnv_x_min) ||
0230 (m_cnv_y_num<=0) || (m_cnv_y_max<=m_cnv_y_min) ) {
0231 double dx = 0.01 * (upper_edge_x() - lower_edge_x())/BINS();
0232 double dy = 0.01 * (upper_edge_y() - lower_edge_y())/BINS();
0233 return convert(BINS(),lower_edge_x(),upper_edge_x()+dx,
0234 BINS(),lower_edge_y(),upper_edge_y()+dy);
0235 } else {
0236 return convert(m_cnv_x_num,m_cnv_x_min,m_cnv_x_max,
0237 m_cnv_y_num,m_cnv_y_min,m_cnv_y_max);
0238 }
0239 }
0240
0241
0242 inline
0243 bool c2d::set_title(const std::string& a_title){
0244 m_title = a_title;
0245 if(m_histo) m_histo->set_title(a_title);
0246 return true;
0247 }
0248
0249 inline
0250 bool c2d::scale(double a_scale) {
0251 if(m_histo) {
0252 return m_histo->scale(a_scale);
0253 } else {
0254 size_t number = m_ws.size();
0255 for(size_t index=0;index<number;index++) m_ws[index] *= a_scale;
0256 m_Sw *= a_scale;
0257 m_Sxw *= a_scale;
0258 m_Sx2w *= a_scale;
0259 m_Syw *= a_scale;
0260 m_Sy2w *= a_scale;
0261 return true;
0262 }
0263 }
0264
0265 inline
0266 bool c2d::reset() {
0267 clear();
0268 delete m_histo;
0269 m_histo = 0;
0270 return true;
0271 }
0272
0273 inline
0274 bool c2d::fill(double aX,double aY,double aW){
0275 if(!m_histo && (m_limit!=UNLIMITED()) && ((int)m_xs.size()>=m_limit)){
0276 convert_to_histogram();
0277 }
0278
0279 if(m_histo) {
0280 return m_histo->fill(aX,aY,aW);
0281 } else {
0282 if(m_xs.size()) {
0283 m_lower_x = mn<double>(aX,m_lower_x);
0284 m_upper_x = mx<double>(aX,m_upper_x);
0285 } else {
0286 m_lower_x = aX;
0287 m_upper_x = aX;
0288 }
0289 if(m_ys.size()) {
0290 m_lower_y = mn<double>(aY,m_lower_y);
0291 m_upper_y = mx<double>(aY,m_upper_y);
0292 } else {
0293 m_lower_y = aY;
0294 m_upper_y = aY;
0295 }
0296 m_xs.push_back(aX);
0297 m_ys.push_back(aY);
0298 m_ws.push_back(aW);
0299 m_Sw += aW;
0300 double xw = aX * aW;
0301 m_Sxw += xw;
0302 m_Sx2w += aX * xw;
0303 double yw = aY * aW;
0304 m_Syw += yw;
0305 m_Sy2w += aY * yw;
0306 return true;
0307 }
0308 }
0309
0310 inline
0311 bool c2d::convert(const std::vector<double>& a_edges_x,const std::vector<double>& a_edges_y) {
0312 if(m_histo) return true;
0313 m_histo = new histo::h2d(base_cloud::title(),
0314 a_edges_x,a_edges_y);
0315 if(!m_histo) return false;
0316 bool status = fill_histogram(*m_histo);
0317 clear();
0318 return status;
0319 }
0320
0321 inline
0322 bool c2d::set_conversion_parameters(
0323 unsigned int aCnvXnumber,double aCnvXmin,double aCnvXmax
0324 ,unsigned int aCnvYnumber,double aCnvYmin,double aCnvYmax
0325 ){
0326 m_cnv_x_num = aCnvXnumber;
0327 m_cnv_x_min = aCnvXmin;
0328 m_cnv_x_max = aCnvXmax;
0329 m_cnv_y_num = aCnvYnumber;
0330 m_cnv_y_min = aCnvYmin;
0331 m_cnv_y_max = aCnvYmax;
0332 return true;
0333 }
0334
0335 inline
0336 const h2d& c2d::histogram() const {
0337 if(!m_histo) const_cast<c2d&>(*this).convert_to_histogram();
0338 return *m_histo;
0339 }
0340
0341 inline
0342 double c2d::rms_x() const {
0343 double rms = 0; //FIXME nan.
0344 if(m_histo) {
0345 rms = m_histo->rms_x();
0346 } else {
0347 if(m_Sw==0) {
0348 } else {
0349 double mean = m_Sxw / m_Sw;
0350 rms = ::sqrt(::fabs( (m_Sx2w / m_Sw) - mean * mean));
0351 }
0352 }
0353 return rms;
0354 }
0355 inline
0356 double c2d::rms_y() const {
0357 double rms = 0; //FIXME nan.
0358 if(m_histo) {
0359 rms = m_histo->rms_y();
0360 } else {
0361 if(m_Sw==0) {
0362 } else {
0363 double mean = m_Syw / m_Sw;
0364 rms = ::sqrt(::fabs( (m_Sy2w / m_Sw) - mean * mean));
0365 }
0366 }
0367 return rms;
0368 }
0369
0370 }}
0371
0372 #endif
