Geant4 Cross Reference |
1 // Copyright (C) 2010, Guy Barrand. All rights 2 // See the file tools.license for terms. 3 4 #ifndef tools_histo_c2d 5 #define tools_histo_c2d 6 7 #include "base_cloud" 8 9 #include "../mnmx" 10 11 #include "h2d" 12 13 namespace tools { 14 namespace histo { 15 16 class c2d : 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 2;} 25 bool reset(); 26 unsigned int entries() const { return m_hist 27 public: 28 double sum_of_weights() const { return (m_hi 29 bool convert_to_histogram(); 30 bool is_converted() const {return m_histo ? 31 bool scale(double); 32 public: 33 bool fill(double,double,double = 1); 34 double lower_edge_x() const {return m_histo 35 double lower_edge_y() const { return m_histo 36 double upper_edge_x() const { return m_histo 37 double upper_edge_y() const { return m_histo 38 double value_x(unsigned int a_index) const { 39 double value_y(unsigned int a_index) const { 40 double weight(unsigned int a_index) const { 41 double mean_x() const {return m_histo ? m_hi 42 double mean_y() const {return m_histo ? m_hi 43 double rms_x() const; 44 double rms_y() const; 45 bool convert(unsigned int,double,double,unsi 46 bool convert(const std::vector<double>&,cons 47 const histo::h2d& histogram() const; 48 template <class HISTO> 49 bool fill_histogram(HISTO& a_histo) const { 50 size_t number = m_xs.size(); 51 for(size_t index=0;index<number;index++) { 52 if(!a_histo.fill(m_xs[index],m_ys[index] 53 } 54 return true; 55 } 56 bool set_conversion_parameters(unsigned int, 57 58 bool set_histogram(h2d* a_histo){ //we take 59 reset(); 60 m_histo = a_histo; 61 return true; 62 } 63 public: 64 c2d(); 65 c2d(const std::string&,int=base_cloud::UNLIM 66 virtual ~c2d(){delete m_histo;} 67 public: 68 c2d(const c2d& a_from) 69 :base_cloud(a_from) 70 ,m_xs(a_from.m_xs) 71 ,m_ys(a_from.m_ys) 72 ,m_lower_x(a_from.m_lower_x) 73 ,m_upper_x(a_from.m_upper_x) 74 ,m_lower_y(a_from.m_lower_y) 75 ,m_upper_y(a_from.m_upper_y) 76 ,m_Sxw(a_from.m_Sxw) 77 ,m_Sx2w(a_from.m_Sx2w) 78 ,m_Syw(a_from.m_Syw) 79 ,m_Sy2w(a_from.m_Sy2w) 80 ,m_cnv_x_num(a_from.m_cnv_x_num) 81 ,m_cnv_x_min(a_from.m_cnv_x_min) 82 ,m_cnv_x_max(a_from.m_cnv_x_max) 83 ,m_cnv_y_num(a_from.m_cnv_y_num) 84 ,m_cnv_y_min(a_from.m_cnv_y_min) 85 ,m_cnv_y_max(a_from.m_cnv_y_max) 86 ,m_histo(0) 87 { 88 if(a_from.m_histo) { 89 m_histo = new histo::h2d(*a_from.m_histo 90 } 91 } 92 93 c2d& operator=(const c2d& a_from) { 94 base_cloud::operator=(a_from); 95 if(&a_from==this) return *this; 96 m_xs = a_from.m_xs; 97 m_ys = a_from.m_ys; 98 m_lower_x = a_from.m_lower_x; 99 m_upper_x = a_from.m_upper_x; 100 m_lower_y = a_from.m_lower_y; 101 m_upper_y = a_from.m_upper_y; 102 m_Sxw = a_from.m_Sxw; 103 m_Sx2w = a_from.m_Sx2w; 104 m_Syw = a_from.m_Syw; 105 m_Sy2w = a_from.m_Sy2w; 106 m_cnv_x_num = a_from.m_cnv_x_num; 107 m_cnv_x_min = a_from.m_cnv_x_min; 108 m_cnv_x_max = a_from.m_cnv_x_max; 109 m_cnv_y_num = a_from.m_cnv_y_num; 110 m_cnv_y_min = a_from.m_cnv_y_min; 111 m_cnv_y_max = a_from.m_cnv_y_max; 112 delete m_histo; 113 m_histo = 0; 114 if(a_from.m_histo) { 115 m_histo = new histo::h2d(*a_from.m_histo 116 } 117 return *this; 118 } 119 public: //AIDA API 120 double lowerEdgeX() const {return lower_edge 121 double lowerEdgeY() const {return lower_edge 122 double upperEdgeX() const {return upper_edge 123 double upperEdgeY() const {return upper_edge 124 template <class HISTO> 125 bool fillHistogram(HISTO& a_histo) const {re 126 protected: 127 void clear(); 128 protected: 129 std::vector<double> m_xs; 130 std::vector<double> m_ys; 131 double m_lower_x; 132 double m_upper_x; 133 double m_lower_y; 134 double m_upper_y; 135 double m_Sxw; 136 double m_Sx2w; 137 double m_Syw; 138 double m_Sy2w; 139 // 140 unsigned int m_cnv_x_num; 141 double m_cnv_x_min; 142 double m_cnv_x_max; 143 unsigned int m_cnv_y_num; 144 double m_cnv_y_min; 145 double m_cnv_y_max; 146 histo::h2d* m_histo; 147 }; 148 149 }} 150 151 namespace tools { 152 namespace histo { 153 154 inline 155 c2d::c2d() 156 :base_cloud(UNLIMITED()) 157 ,m_lower_x(0) 158 ,m_upper_x(0) 159 ,m_lower_y(0) 160 ,m_upper_y(0) 161 ,m_Sxw(0) 162 ,m_Sx2w(0) 163 ,m_Syw(0) 164 ,m_Sy2w(0) 165 ,m_cnv_x_num(0) 166 ,m_cnv_x_min(0) 167 ,m_cnv_x_max(0) 168 ,m_cnv_y_num(0) 169 ,m_cnv_y_min(0) 170 ,m_cnv_y_max(0) 171 ,m_histo(0) 172 {} 173 174 inline 175 c2d::c2d(const std::string& a_title,int aLimit 176 :base_cloud(aLimit) 177 ,m_lower_x(0) 178 ,m_upper_x(0) 179 ,m_lower_y(0) 180 ,m_upper_y(0) 181 ,m_Sxw(0) 182 ,m_Sx2w(0) 183 ,m_Syw(0) 184 ,m_Sy2w(0) 185 ,m_cnv_x_num(0) 186 ,m_cnv_x_min(0) 187 ,m_cnv_x_max(0) 188 ,m_cnv_y_num(0) 189 ,m_cnv_y_min(0) 190 ,m_cnv_y_max(0) 191 ,m_histo(0) 192 { 193 set_title(a_title); 194 } 195 196 inline 197 void c2d::clear(){ 198 m_lower_x = 0; 199 m_upper_x = 0; 200 m_lower_y = 0; 201 m_upper_y = 0; 202 m_Sw = 0; 203 m_Sxw = 0; 204 m_Sx2w = 0; 205 m_Syw = 0; 206 m_Sy2w = 0; 207 m_xs.clear(); 208 m_ys.clear(); 209 m_ws.clear(); 210 } 211 212 inline 213 bool c2d::convert( 214 unsigned int a_bins_x,double a_lower_edge_x,d 215 ,unsigned int a_bins_y,double a_lower_edge_y,d 216 ) { 217 if(m_histo) return true; // Done. 218 m_histo = new histo::h2d(base_cloud::title() 219 a_bins_x,a_lower_ed 220 a_bins_y,a_lower_ed 221 if(!m_histo) return false; 222 bool status = fill_histogram(*m_histo); 223 clear(); 224 return status; 225 } 226 227 inline 228 bool c2d::convert_to_histogram(){ 229 if( (m_cnv_x_num<=0) || (m_cnv_x_max<=m_cnv_ 230 (m_cnv_y_num<=0) || (m_cnv_y_max<=m_cnv_ 231 double dx = 0.01 * (upper_edge_x() - lower 232 double dy = 0.01 * (upper_edge_y() - lower 233 return convert(BINS(),lower_edge_x(),upper 234 BINS(),lower_edge_y(),upper 235 } else { 236 return convert(m_cnv_x_num,m_cnv_x_min,m_c 237 m_cnv_y_num,m_cnv_y_min,m_c 238 } 239 } 240 241 242 inline 243 bool c2d::set_title(const std::string& a_title 244 m_title = a_title; 245 if(m_histo) m_histo->set_title(a_title); 246 return true; 247 } 248 249 inline 250 bool c2d::scale(double a_scale) { 251 if(m_histo) { 252 return m_histo->scale(a_scale); 253 } else { 254 size_t number = m_ws.size(); 255 for(size_t index=0;index<number;index++) m 256 m_Sw *= a_scale; 257 m_Sxw *= a_scale; 258 m_Sx2w *= a_scale; 259 m_Syw *= a_scale; 260 m_Sy2w *= a_scale; 261 return true; 262 } 263 } 264 265 inline 266 bool c2d::reset() { 267 clear(); 268 delete m_histo; 269 m_histo = 0; 270 return true; 271 } 272 273 inline 274 bool c2d::fill(double aX,double aY,double aW){ 275 if(!m_histo && (m_limit!=UNLIMITED()) && ((i 276 convert_to_histogram(); 277 } 278 279 if(m_histo) { 280 return m_histo->fill(aX,aY,aW); 281 } else { 282 if(m_xs.size()) { 283 m_lower_x = mn<double>(aX,m_lower_x); 284 m_upper_x = mx<double>(aX,m_upper_x); 285 } else { 286 m_lower_x = aX; 287 m_upper_x = aX; 288 } 289 if(m_ys.size()) { 290 m_lower_y = mn<double>(aY,m_lower_y); 291 m_upper_y = mx<double>(aY,m_upper_y); 292 } else { 293 m_lower_y = aY; 294 m_upper_y = aY; 295 } 296 m_xs.push_back(aX); 297 m_ys.push_back(aY); 298 m_ws.push_back(aW); 299 m_Sw += aW; 300 double xw = aX * aW; 301 m_Sxw += xw; 302 m_Sx2w += aX * xw; 303 double yw = aY * aW; 304 m_Syw += yw; 305 m_Sy2w += aY * yw; 306 return true; 307 } 308 } 309 310 inline 311 bool c2d::convert(const std::vector<double>& a 312 if(m_histo) return true; 313 m_histo = new histo::h2d(base_cloud::title() 314 a_edges_x,a_edges_y 315 if(!m_histo) return false; 316 bool status = fill_histogram(*m_histo); 317 clear(); 318 return status; 319 } 320 321 inline 322 bool c2d::set_conversion_parameters( 323 unsigned int aCnvXnumber,double aCnvXmin,doub 324 ,unsigned int aCnvYnumber,double aCnvYmin,doub 325 ){ 326 m_cnv_x_num = aCnvXnumber; 327 m_cnv_x_min = aCnvXmin; 328 m_cnv_x_max = aCnvXmax; 329 m_cnv_y_num = aCnvYnumber; 330 m_cnv_y_min = aCnvYmin; 331 m_cnv_y_max = aCnvYmax; 332 return true; 333 } 334 335 inline 336 const h2d& c2d::histogram() const { 337 if(!m_histo) const_cast<c2d&>(*this).convert 338 return *m_histo; 339 } 340 341 inline 342 double c2d::rms_x() const { 343 double rms = 0; //FIXME nan. 344 if(m_histo) { 345 rms = m_histo->rms_x(); 346 } else { 347 if(m_Sw==0) { 348 } else { 349 double mean = m_Sxw / m_Sw; 350 rms = ::sqrt(::fabs( (m_Sx2w / m_Sw) - m 351 } 352 } 353 return rms; 354 } 355 inline 356 double c2d::rms_y() const { 357 double rms = 0; //FIXME nan. 358 if(m_histo) { 359 rms = m_histo->rms_y(); 360 } else { 361 if(m_Sw==0) { 362 } else { 363 double mean = m_Syw / m_Sw; 364 rms = ::sqrt(::fabs( (m_Sy2w / m_Sw) - m 365 } 366 } 367 return rms; 368 } 369 370 }} 371 372 #endif