Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/externals/g4tools/include/tools/histo/c3d

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  1 // Copyright (C) 2010, Guy Barrand. All rights reserved.
  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::c3d");
 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<double>&,
 55                        const std::vector<double>&);
 56   const histo::h3d& histogram() const;
 57   bool fill_histogram(histo::h3d& a_histo) const {
 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],m_zs[index],m_ws[index])) return false;
 61     }
 62     return true;
 63   }
 64   bool set_conversion_parameters(unsigned int,double,double,
 65                                  unsigned int,double,double,
 66                                  unsigned int,double,double);
 67 
 68   bool set_histogram(h3d* a_histo){ //we take ownership of a_histo.
 69     reset();
 70     m_histo = a_histo;
 71     return true;
 72   }
 73 public:
 74   c3d();
 75   c3d(const std::string&,int=base_cloud::UNLIMITED());
 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 ? true : false;}
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,double a_upper_edge_x
265 ,unsigned int a_bins_y,double a_lower_edge_y,double a_upper_edge_y
266 ,unsigned int a_bins_z,double a_lower_edge_z,double a_upper_edge_z
267 ) {
268   if(m_histo) return true; // Done.
269   m_histo = new histo::h3d(base_cloud::title(),
270                            a_bins_x,a_lower_edge_x,a_upper_edge_x,
271                            a_bins_y,a_lower_edge_y,a_upper_edge_y,
272                            a_bins_z,a_lower_edge_z,a_upper_edge_z);
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_x_min) ||
282       (m_cnv_y_num<=0) || (m_cnv_y_max<=m_cnv_y_min) ||
283       (m_cnv_z_num<=0) || (m_cnv_z_max<=m_cnv_z_min) ) {
284     double dx = 0.01 * (upper_edge_x() - lower_edge_x())/BINS();
285     double dy = 0.01 * (upper_edge_y() - lower_edge_y())/BINS();
286     double dz = 0.01 * (upper_edge_z() - lower_edge_z())/BINS();
287     return convert(BINS(),lower_edge_x(),upper_edge_x()+dx,
288                    BINS(),lower_edge_y(),upper_edge_y()+dy,
289                    BINS(),lower_edge_z(),upper_edge_z()+dz);
290   } else {
291     return convert(m_cnv_x_num,m_cnv_x_min,m_cnv_x_max,
292                    m_cnv_y_num,m_cnv_y_min,m_cnv_y_max,
293                    m_cnv_z_num,m_cnv_z_min,m_cnv_z_max);
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_ws[index] *= a_scale;
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,double aCnvXmax
325 ,unsigned int aCnvYnumber,double aCnvYmin,double aCnvYmax
326 ,unsigned int aCnvZnumber,double aCnvZmin,double aCnvZmax
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_to_histogram();
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,double aW){
357   if(!m_histo && (m_limit!=UNLIMITED()) && ((int)m_xs.size()>=m_limit)){
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,a_edges_z);
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() : m_Sw);
421 }
422 
423 inline
424 unsigned int c3d::entries() const {
425   return m_histo ? m_histo->all_entries() : (unsigned int)m_ws.size();
426 }
427 
428 inline
429 double c3d::lower_edge_x() const {
430   return m_histo ? m_histo->axis_x().lower_edge() : m_lower_x;
431 }
432 inline
433 double c3d::lower_edge_y() const {
434   return m_histo ? m_histo->axis_y().lower_edge() : m_lower_y;
435 }
436 inline
437 double c3d::lower_edge_z() const {
438   return m_histo ? m_histo->axis_z().lower_edge() : m_lower_z;
439 }
440 inline
441 double c3d::upper_edge_x() const {
442   return m_histo ? m_histo->axis_x().upper_edge() : m_upper_x;
443 }
444 inline
445 double c3d::upper_edge_y() const {
446   return m_histo ? m_histo->axis_y().upper_edge() : m_upper_y;
447 }
448 inline
449 double c3d::upper_edge_z() const {
450   return m_histo ? m_histo->axis_z().upper_edge() : m_upper_z;
451 }
452 inline
453 double c3d::value_x(unsigned int a_index) const {
454   return m_histo ? 0 : m_xs[a_index];
455 }
456 inline
457 double c3d::value_y(unsigned int a_index) const {
458   return m_histo ? 0 : m_ys[a_index];
459 }
460 inline
461 double c3d::value_z(unsigned int a_index) const {
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_Sxw/m_Sw:0);
471 }
472 inline
473 double c3d::mean_y() const {
474   return m_histo ? m_histo->mean_y() : (m_Sw?m_Syw/m_Sw:0);
475 }
476 inline
477 double c3d::mean_z() const {
478   return m_histo ? m_histo->mean_z() : (m_Sw?m_Szw/m_Sw:0);
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) - mean * mean));
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) - mean * mean));
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) - mean * mean));
518     }
519   }
520   return rms;
521 }
522 
523 }}
524 
525 #endif