Geant4 Cross Reference

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

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_base_histo
  5 #define tools_histo_base_histo
  6 
  7 #ifdef TOOLS_MEM
  8 #include "../mem"
  9 #endif
 10 
 11 #include "histo_data"
 12 
 13 #include <cmath>
 14 #include <map> //for annotations
 15 #include <ostream>
 16 
 17 namespace tools {
 18 namespace histo {
 19 
 20 //TC is for a coordinate.
 21 //TO is for an offset used to identify a bin.
 22 //TN is for a number of entries.
 23 //TW is for a weight.
 24 //TH is for a height.
 25 
 26 template <class TC,class TO,class TN,class TW,class TH>
 27 class base_histo : protected histo_data<TC,TO,TN,TW> {
 28   typedef histo_data<TC,TO,TN,TW> parent;
 29 private:
 30   static const std::string& s_class() {
 31     static const std::string s_v("tools::histo::base_histo");
 32     return s_v;
 33   }
 34 public:
 35   typedef histo_data<TC,TO,TN,TW> hd_t;
 36   typedef axis<TC,TO> axis_t;
 37   typedef typename axis_t::bn_t bn_t;
 38   typedef unsigned int dim_t;
 39   typedef TC coordinate_t;
 40   typedef TO offset_t;
 41   typedef TN num_entries_t;
 42   typedef TW weight_t;
 43   typedef TH height_t;
 44 protected:
 45   virtual TH get_bin_height(TO) const = 0;  //histo/profile
 46 protected:
 47   void base_from_data(const hd_t& a_from) {parent::operator=(a_from);}
 48 #ifdef tools_histo_base_histo //for backward compatibility with tools
 49   hd_t base_get_data() const {
 50     hd_t hd;
 51     hd = *this;
 52     return hd;
 53   }
 54 #endif
 55 public:
 56   const hd_t& dac() const {return *this;} //data accessor.
 57 protected:
 58   base_histo():parent() {
 59 #ifdef TOOLS_MEM
 60     mem::increment(s_class().c_str());
 61 #endif
 62   }
 63 protected:
 64   virtual ~base_histo() {
 65 #ifdef TOOLS_MEM
 66     mem::decrement(s_class().c_str());
 67 #endif
 68   }
 69 protected:
 70   base_histo(const base_histo& a_from):parent(a_from) {
 71 #ifdef TOOLS_MEM
 72     mem::increment(s_class().c_str());
 73 #endif
 74   }
 75 
 76   base_histo& operator=(const base_histo& a_from) {
 77     if(&a_from==this) return *this;
 78     parent::operator=(a_from);
 79     return *this;
 80   }
 81 
 82 public:
 83   bool equals(const base_histo& a_from,const TW& a_prec,TW(*a_fabs)(TW)) const {
 84     return parent::equals(a_from,a_prec,a_fabs);
 85   }
 86 
 87   const std::string& title() const {return parent::m_title;}
 88   bool set_title(const std::string& a_title){parent::m_title = a_title;return true;}
 89   dim_t dimension() const {return parent::m_dimension;}
 90   dim_t number_of_planes() const {return dim_planes(parent::m_dimension);}
 91 
 92   TN entries() const {
 93     return parent::m_in_range_entries; //not set if reading a TH from a CERN-ROOT file.
 94   }
 95   TN all_entries() const {
 96     return parent::m_all_entries; //works also is reading histo from a CERN-ROOT file.
 97   }
 98 
 99   TN extra_entries() const {
100     return parent::m_all_entries-parent::m_in_range_entries; //works also is reading histo from a CERN-ROOT file.
101   }
102   TW equivalent_bin_entries() const {
103     TW sw = 0;
104     TW sw2 = 0;
105     for(TO ibin=0;ibin<parent::m_bin_number;ibin++) {
106       if(!histo::is_out(parent::m_axes,ibin)) {
107         sw += parent::m_bin_Sw[ibin];
108         sw2 += parent::m_bin_Sw2[ibin];
109       }
110     }
111     if(sw2==0) return 0;
112     return (sw * sw)/sw2;
113   }
114   TH sum_bin_heights() const {
115     TH sh = 0;
116     for(TO ibin=0;ibin<parent::m_bin_number;ibin++) {
117       if(!histo::is_out(parent::m_axes,ibin)) {
118         sh += get_bin_height(ibin);
119       }
120     }
121     return sh;
122   }
123   TH sum_all_bin_heights() const {
124     TH sh = 0;
125     for(TO ibin=0;ibin<parent::m_bin_number;ibin++) {
126       sh += get_bin_height(ibin);
127     }
128     return sh;
129   }
130 
131   TH sum_extra_bin_heights() const {
132     TH sh = 0;
133     for(TO ibin=0;ibin<parent::m_bin_number;ibin++) {
134       if(histo::is_out(parent::m_axes,ibin)) {
135         sh += get_bin_height(ibin);
136       }
137     }
138     return sh;
139   }
140 
141   TH min_bin_height() const {
142     TH value = 0;
143     bool first = true;
144     for(TO ibin=0;ibin<parent::m_bin_number;ibin++) {
145       if(!histo::is_out(parent::m_axes,ibin)) {
146         TH vbin = get_bin_height(ibin);
147         if(first) {
148           first = false;
149           value = vbin;
150         } else {
151           if(vbin<=value) value = vbin;
152         }
153       }
154     }
155     return value;
156   }
157 
158   TH max_bin_height() const {
159     TH value = 0;
160     bool first = true;
161     for(TO ibin=0;ibin<parent::m_bin_number;ibin++) {
162       if(!histo::is_out(parent::m_axes,ibin)) {
163         TH vbin = get_bin_height(ibin);
164         if(first) {
165           first = false;
166           value = vbin;
167         } else {
168           if(vbin>=value) value = vbin;
169         }
170       }
171     }
172     return value;
173   }
174 
175   bool min_bin_height_with_entries(TH& a_value) const {
176     TH value = 0;
177     bool first = true;
178     for(TO ibin=0;ibin<parent::m_bin_number;ibin++) {
179       if(!histo::is_out(parent::m_axes,ibin) && (parent::m_bin_entries[ibin]>0) ) {
180         TH vbin = get_bin_height(ibin);
181         if(first) {
182           first = false;
183           value = vbin;
184         } else {
185           if(vbin<=value) value = vbin;
186         }
187       }
188     }
189     a_value = value;
190     return first?false:true; //return true if at least one bin with entries processed.
191   }
192 
193   bool max_bin_height_with_entries(TH& a_value) const {
194     TH value = 0;
195     bool first = true;
196     for(TO ibin=0;ibin<parent::m_bin_number;ibin++) {
197       if(!histo::is_out(parent::m_axes,ibin) && (parent::m_bin_entries[ibin]>0) ) {
198         TH vbin = get_bin_height(ibin);
199         if(first) {
200           first = false;
201           value = vbin;
202         } else {
203           if(vbin>=value) value = vbin;
204         }
205       }
206     }
207     a_value = value;
208     return first?false:true; //return true if at least one bin with entries processed.
209   }
210 
211   bool has_entries_per_bin() const { //to detect histos coming from TH streaming out of a root file.
212     // it assumes that update_fast_getters() had been applied.
213     if(parent::m_in_range_entries) return true;
214     // may be a from-root histo :
215     if(parent::m_in_range_Sw) return false;
216     // no in range entries and weight :
217     return true; //for exa not filled = ok.
218   }
219 
220 public: //histo_data
221   bool get_ith_axis_Sxw(dim_t a_axis,TC& a_value) const {
222     a_value = 0;
223     if(a_axis>=parent::m_dimension) return false;
224     for(TO ibin=0;ibin<parent::m_bin_number;ibin++) {
225       if(!histo::is_out(parent::m_axes,ibin)) {
226         a_value += parent::m_bin_Sxw[ibin][a_axis];
227       }
228     }
229     return true;
230   }
231 
232   bool get_ith_axis_Sx2w(dim_t a_axis,TC& a_value) const {
233     a_value = 0;
234     if(a_axis>=parent::m_dimension) return false;
235     for(TO ibin=0;ibin<parent::m_bin_number;ibin++) {
236       if(!histo::is_out(parent::m_axes,ibin)) {
237         a_value += parent::m_bin_Sx2w[ibin][a_axis];
238       }
239     }
240     return true;
241   }
242 
243   TW get_in_range_Sw() const {return parent::m_in_range_Sw;}   //for CERN-ROOT file writing.
244   TW get_in_range_Sw2() const {return parent::m_in_range_Sw2;} //for CERN-ROOT file writing.
245 
246   void get_Sw_Sw2(TW& a_sw,TW& a_sw2) const {
247     a_sw = 0;
248     a_sw2 = 0;
249     for(TO ibin=0;ibin<parent::m_bin_number;ibin++) {
250       if(!histo::is_out(parent::m_axes,ibin)) {
251         a_sw += parent::m_bin_Sw[ibin];
252         a_sw2 += parent::m_bin_Sw2[ibin];
253       }
254     }
255   }
256   void get_all_Sw_Sw2(TW& a_sw,TW& a_sw2) const {
257     a_sw = 0;
258     a_sw2 = 0;
259     for(TO ibin=0;ibin<parent::m_bin_number;ibin++) {
260       a_sw += parent::m_bin_Sw[ibin];
261       a_sw2 += parent::m_bin_Sw2[ibin];
262     }
263   }
264 
265 /*
266   TW get_all_Sw() const {
267     TW sw = 0;
268     for(TO ibin=0;ibin<m_bin_number;ibin++) sw += m_bin_Sw[ibin];
269     return sw;
270   }
271   TN get_all_entries() const {
272     TN number = 0;
273     for(TO ibin=0;ibin<m_bin_number;ibin++) {
274       number += m_bin_entries[ibin];
275     }
276     return number;
277   }
278 
279   // for tools/wroot/streamers :
280   TN get_entries() const {
281     TN number = 0;
282     for(TO ibin=0;ibin<m_bin_number;ibin++) {
283       if(!histo::is_out(m_axes,ibin)) {
284         number += m_bin_entries[ibin];
285       }
286     }
287     return number;
288   }
289 */
290 protected:
291   enum {AxisX=0,AxisY=1,AxisZ=2};
292 
293   bool configure(dim_t a_dim,
294                  const std::vector<bn_t>& aNumbers,
295                  const std::vector<TC>& aMins,
296                  const std::vector<TC>& aMaxs) {
297     // Clear :
298     parent::m_bin_entries.clear();
299     parent::m_bin_Sw.clear();
300     parent::m_bin_Sw2.clear();
301     parent::m_bin_Sxw.clear();
302     parent::m_bin_Sx2w.clear();
303     parent::m_in_range_Sxw.clear();
304     parent::m_in_range_Sx2w.clear();
305     parent::m_axes.clear();
306     parent::m_in_range_plane_Sxyw.clear();
307     parent::m_annotations.clear();
308 
309     parent::m_bin_number = 0;
310     parent::m_dimension = 0;
311     parent::m_all_entries = 0;
312     parent::m_in_range_entries = 0;
313     parent::m_in_range_Sw = 0;
314     parent::m_in_range_Sw2 = 0;
315     parent::m_in_range_Sxw.resize(a_dim,0);
316     parent::m_in_range_Sx2w.resize(a_dim,0);
317 
318     // Some checks :
319     if(!a_dim) return false;
320     parent::m_axes.resize(a_dim);
321     // Setup axes :
322     for(dim_t iaxis=0;iaxis<a_dim;iaxis++) {
323       if(!parent::m_axes[iaxis].configure(aNumbers[iaxis],aMins[iaxis],aMaxs[iaxis])) {
324         // do not do :
325         //   m_axes.clear()
326         // so that :
327         //   b1::axis(),b2::axis_[x,y]()
328         // do not crash in case of a bad booking.
329         //m_axes.clear();
330         return false;
331       }
332     }
333 
334     parent::m_dimension = a_dim;
335 
336     base_allocate(); //set m_bin_number.
337 
338     return true;
339   }
340 
341   bool configure(dim_t a_dim,const std::vector< std::vector<TC> >& a_edges) {
342     // Clear :
343     parent::m_bin_entries.clear();
344     parent::m_bin_Sw.clear();
345     parent::m_bin_Sw2.clear();
346     parent::m_bin_Sxw.clear();
347     parent::m_bin_Sx2w.clear();
348     parent::m_in_range_Sxw.clear();
349     parent::m_in_range_Sx2w.clear();
350     parent::m_axes.clear();
351     parent::m_in_range_plane_Sxyw.clear();
352     parent::m_annotations.clear();
353 
354     parent::m_bin_number = 0;
355     parent::m_dimension = 0;
356     parent::m_all_entries = 0;
357     parent::m_in_range_entries = 0;
358     parent::m_in_range_Sw = 0;
359     parent::m_in_range_Sw2 = 0;
360     parent::m_in_range_Sxw.resize(a_dim,0);
361     parent::m_in_range_Sx2w.resize(a_dim,0);
362 
363     // Some checks :
364     if(!a_dim) return false;
365     parent::m_axes.resize(a_dim);
366     // Setup axes :
367     for(dim_t iaxis=0;iaxis<a_dim;iaxis++) {
368       if(!parent::m_axes[iaxis].configure(a_edges[iaxis])) {
369         //m_axes.clear();
370         return false;
371       }
372     }
373 
374     parent::m_dimension = a_dim;
375 
376     base_allocate(); //set m_bin_number.
377 
378     return true;
379   }
380 
381   void base_reset() {
382     // Reset content (different of clear that deallocate all internal things).
383     for(TO ibin=0;ibin<parent::m_bin_number;ibin++) {
384       parent::m_bin_entries[ibin] = 0;
385       parent::m_bin_Sw[ibin] = 0;
386       parent::m_bin_Sw2[ibin] = 0;
387       for(dim_t iaxis=0;iaxis<parent::m_dimension;iaxis++) {
388         parent::m_bin_Sxw[ibin][iaxis] = 0;
389         parent::m_bin_Sx2w[ibin][iaxis] = 0;
390       }
391     }
392     parent::m_in_range_plane_Sxyw.assign(dim_planes(parent::m_dimension),0);
393     //profile not done here.
394     parent::reset_fast_getters();
395   }
396 
397 protected:
398   void base_allocate() {
399     dim_t iaxis;
400     // Add two bins for the [under,out]flow data.
401     TO n_bin = 1;
402     for(iaxis=0;iaxis<parent::m_dimension;iaxis++) {
403       n_bin *= (parent::m_axes[iaxis].bins() + 2);
404     }
405 
406     parent::m_bin_entries.resize(n_bin,0);
407     parent::m_bin_Sw.resize(n_bin,0);
408     parent::m_bin_Sw2.resize(n_bin,0);
409 
410     std::vector<TC> empty;
411     empty.resize(parent::m_dimension,0);
412     parent::m_bin_Sxw.resize(n_bin,empty);
413     parent::m_bin_Sx2w.resize(n_bin,empty);
414 
415     parent::m_bin_number = n_bin; // All bins : [in-range, underflow, outflow] bins.
416 
417     parent::m_axes[0].m_offset = 1;
418     for(iaxis=1;iaxis<parent::m_dimension;iaxis++) {
419       parent::m_axes[iaxis].m_offset = parent::m_axes[iaxis-1].m_offset * (parent::m_axes[iaxis-1].bins()+2);
420     }
421 
422     parent::m_in_range_plane_Sxyw.resize(dim_planes(parent::m_dimension),0);
423   }
424 
425 public:
426   // to access data from methods :
427   const std::vector<TN>& bins_entries() const {return parent::m_bin_entries;}
428   const std::vector<TW>& bins_sum_w() const {return parent::m_bin_Sw;}
429   const std::vector<TW>& bins_sum_w2() const {return parent::m_bin_Sw2;}
430   const std::vector< std::vector<TC> >& bins_sum_xw() const {return parent::m_bin_Sxw;}
431   const std::vector< std::vector<TC> >& bins_sum_x2w() const {return parent::m_bin_Sx2w;}
432   const std::vector<TC>& in_range_planes_xyw() const {return parent::m_in_range_plane_Sxyw;}
433 
434 public:
435   const axis_t& get_axis(int a_index) const {return parent::m_axes[a_index];}
436   offset_t get_bins() const {return parent::m_bin_number;}
437   const std::string& get_title() const {return parent::m_title;}
438   dim_t get_dimension() const {return parent::m_dimension;}
439   bool is_valid() const {return (parent::m_dimension?true:false);}
440 
441 public: //annotations :
442   typedef std::map<std::string,std::string> annotations_t;
443   const annotations_t& annotations() const {return parent::m_annotations;}
444   annotations_t annotations() {return parent::m_annotations;}
445 
446   void add_annotation(const std::string& a_key,const std::string& a_value) {
447     parent::m_annotations[a_key] = a_value; //override if a_key already exists.
448   }
449   bool annotation(const std::string& a_key,std::string& a_value) const {
450     annotations_t::const_iterator it = parent::m_annotations.find(a_key);
451     if(it==parent::m_annotations.end()) {a_value.clear();return false;}
452     a_value = (*it).second;
453     return true;
454   }
455 
456   void set_annotations(const annotations_t& a_annotations) {parent::m_annotations = a_annotations;}
457 
458   void hprint_annotations(std::ostream& a_out) {
459     a_out << " * ANNOTATIONS :" << std::endl;
460     annotations_t::const_iterator it;
461     for(it=parent::m_annotations.begin();it!=parent::m_annotations.end();++it) {
462     //out << " * (" << index << ") "
463       a_out << " *  " << (*it).first << " = " << (*it).second << std::endl;
464     }
465   }
466 
467 protected:
468   bool is_compatible(const base_histo& a_histo){
469     if(parent::m_dimension!=a_histo.m_dimension) return false;
470     for(dim_t iaxis=0;iaxis<parent::m_dimension;iaxis++) {
471       if(!parent::m_axes[iaxis].is_compatible(a_histo.m_axes[iaxis])) return false;
472     }
473     return true;
474   }
475 
476   void base_add(const base_histo& a_histo){
477     // The only histogram operation that makes sense.
478     for(TO ibin=0;ibin<parent::m_bin_number;ibin++) {
479       parent::m_bin_entries[ibin] += a_histo.m_bin_entries[ibin];
480       parent::m_bin_Sw[ibin] += a_histo.m_bin_Sw[ibin];
481       parent::m_bin_Sw2[ibin] += a_histo.m_bin_Sw2[ibin];
482       for(dim_t iaxis=0;iaxis<parent::m_dimension;iaxis++) {
483         parent::m_bin_Sxw[ibin][iaxis] += a_histo.m_bin_Sxw[ibin][iaxis];
484         parent::m_bin_Sx2w[ibin][iaxis] += a_histo.m_bin_Sx2w[ibin][iaxis];
485       }
486     }
487    {size_t nplane = parent::m_in_range_plane_Sxyw.size();
488     for(size_t iplane=0;iplane<nplane;iplane++)
489       parent::m_in_range_plane_Sxyw[iplane] += a_histo.m_in_range_plane_Sxyw[iplane];}
490     parent::update_fast_getters();
491   }
492 
493   void base_subtract(const base_histo& a_histo) {
494     //ill-defined operation. We keep that because of the "ill-defined past".
495     // We build a new histo with one entry in each bin.
496     for(TO ibin=0;ibin<parent::m_bin_number;ibin++) {
497       parent::m_bin_entries[ibin] = 1;
498 
499       parent::m_bin_Sw[ibin] -= a_histo.m_bin_Sw[ibin];
500       // Yes, it is a += in the below.
501       parent::m_bin_Sw2[ibin] += a_histo.m_bin_Sw2[ibin];
502       for(dim_t iaxis=0;iaxis<parent::m_dimension;iaxis++) {
503         parent::m_bin_Sxw[ibin][iaxis] -= a_histo.m_bin_Sxw[ibin][iaxis];
504         parent::m_bin_Sx2w[ibin][iaxis] -= a_histo.m_bin_Sx2w[ibin][iaxis];
505       }
506     }
507 
508  //{for(dim_t iplane=0;iplane<nplane;iplane++) parent::m_in_range_plane_Sxyw[iplane] ??? a_histo.m_in_range_plane_Sxyw[iplane];}
509 
510     parent::update_fast_getters();
511   }
512 
513   bool base_multiply(const base_histo& a_histo) {
514     //ill-defined operation. We keep that because of the "ill-defined past".
515 
516     // We build a new histo with one entry in each bin of weight :
517     //   this.w * a_histo.w
518     // The current histo is overriden with this new histo.
519     // The m_bin_Sw2 computation is consistent with FreeHEP and CERN-ROOT.
520 
521     if(!is_compatible(a_histo)) return false;
522 
523     std::vector<int> is(parent::m_dimension);
524     for(TO ibin=0;ibin<parent::m_bin_number;ibin++) {
525       TW swa = parent::m_bin_Sw[ibin];
526       TW sw2a = parent::m_bin_Sw2[ibin];
527       TW swb = a_histo.m_bin_Sw[ibin];
528       TW sw2b = a_histo.m_bin_Sw2[ibin];
529       TW sw = swa * swb;
530       parent::m_bin_entries[ibin] = 1;
531       parent::m_bin_Sw[ibin] = sw;
532       parent::m_bin_Sw2[ibin] = sw2a * swb * swb + sw2b * swa * swa;
533       histo::get_indices(parent::m_axes,ibin,is);
534       for(dim_t iaxis=0;iaxis<parent::m_dimension;iaxis++) {
535         TC x = parent::m_axes[iaxis].bin_center(is[iaxis]);
536         parent::m_bin_Sxw[ibin][iaxis] = x * sw;
537         parent::m_bin_Sx2w[ibin][iaxis] = x * x * sw;
538       }
539     }
540 
541  //{for(dim_t iplane=0;iplane<nplane;iplane++) parent::m_in_range_plane_Sxyw[iplane] ??? a_histo.m_in_range_plane_Sxyw[iplane];}
542 
543     parent::update_fast_getters();
544     return true;
545   }
546 
547   bool base_divide(const base_histo& a_histo) {
548     //ill-defined operation. We keep that because of the "ill-defined past".
549 
550     // We build a new histo with one entry in each bin of weight :
551     //   this.w / a_histo.w
552     // The current histo is overriden with this new histo.
553     // The m_bin_Sw2 computation is consistent with FreeHEP and ROOT.
554 
555     if(!is_compatible(a_histo)) return false;
556 
557     std::vector<int> is(parent::m_dimension);
558     for(TO ibin=0;ibin<parent::m_bin_number;ibin++) {
559       histo::get_indices(parent::m_axes,ibin,is);
560       TW swa = parent::m_bin_Sw[ibin];
561       TW swb = a_histo.m_bin_Sw[ibin];
562       TW sw2a = parent::m_bin_Sw2[ibin];
563       TW sw2b = a_histo.m_bin_Sw2[ibin];
564       if(swb!=0) {
565         parent::m_bin_entries[ibin] = 1;
566         TW sw = swa / swb;
567         parent::m_bin_Sw[ibin] = sw;
568         TW swb2 = swb * swb;
569         parent::m_bin_Sw2[ibin] = sw2a / swb2 + sw2b * swa * swa /(swb2*swb2);
570         for(dim_t iaxis=0;iaxis<parent::m_dimension;iaxis++) {
571           TC x = parent::m_axes[iaxis].bin_center(is[iaxis]);
572           parent::m_bin_Sxw[ibin][iaxis] = x * sw;
573           parent::m_bin_Sx2w[ibin][iaxis] = x * x * sw;
574         }
575       } else {
576         parent::m_bin_entries[ibin] = 0;
577         parent::m_bin_Sw[ibin] = 0;
578         parent::m_bin_Sw2[ibin] = 0;
579         for(dim_t iaxis=0;iaxis<parent::m_dimension;iaxis++) {
580           parent::m_bin_Sxw[ibin][iaxis] = 0;
581           parent::m_bin_Sx2w[ibin][iaxis] = 0;
582         }
583       }
584     }
585 
586  //{for(dim_t iplane=0;iplane<nplane;iplane++) parent::m_in_range_plane_Sxyw[iplane] ??? a_histo.m_in_range_plane_Sxyw[iplane];}
587 
588     parent::update_fast_getters();
589     return true;
590   }
591 
592   bool base_multiply(TW a_factor) {
593     if(a_factor<0) return false;
594     TW factor2 = a_factor * a_factor;
595     for(TO ibin=0;ibin<parent::m_bin_number;ibin++) {
596       parent::m_bin_Sw[ibin] *= a_factor;
597       parent::m_bin_Sw2[ibin] *= factor2;
598       for(dim_t iaxis=0;iaxis<parent::m_dimension;iaxis++) {
599         parent::m_bin_Sxw[ibin][iaxis] *= a_factor;
600         parent::m_bin_Sx2w[ibin][iaxis] *= a_factor;
601       }
602     }
603    {size_t nplane = parent::m_in_range_plane_Sxyw.size();
604     for(size_t iplane=0;iplane<nplane;iplane++) parent::m_in_range_plane_Sxyw[iplane] *= a_factor;}
605     parent::update_fast_getters();
606     return true;
607   }
608 
609   bool get_ith_axis_mean(dim_t a_axis,TC& a_value) const {
610     a_value = 0;
611     if(a_axis>=parent::m_dimension) return false;
612     TW sw = 0;
613     TC sxw = 0;
614     for(TO ibin=0;ibin<parent::m_bin_number;ibin++) {
615       if(!histo::is_out(parent::m_axes,ibin)) {
616         sw += parent::m_bin_Sw[ibin];
617         sxw += parent::m_bin_Sxw[ibin][a_axis];
618       }
619     }
620     if(sw==0) return false;
621     a_value = sxw/sw;
622     return true;
623   }
624 
625   bool get_ith_axis_rms(dim_t a_axis,TC& a_value) const {
626     a_value = 0;
627     if(a_axis>=parent::m_dimension) return false;
628     TW sw = 0;
629     TC sxw = 0;
630     TC sx2w = 0;
631     for(TO ibin=0;ibin<parent::m_bin_number;ibin++) {
632       if(!histo::is_out(parent::m_axes,ibin)) {
633         sw += parent::m_bin_Sw[ibin];
634         sxw += parent::m_bin_Sxw[ibin][a_axis];
635         sx2w += parent::m_bin_Sx2w[ibin][a_axis];
636       }
637     }
638     if(sw==0) return false;
639     TC mean = sxw/sw;
640     a_value = ::sqrt(::fabs((sx2w / sw) - mean * mean));
641     return true;
642   }
643 
644   TN get_bin_entries(const std::vector<int>& aIs) const {
645     if(parent::m_bin_number==0) return 0;
646     TO offset;
647     if(!histo::get_offset(parent::m_axes,aIs,offset)) return 0;
648     return parent::m_bin_entries[offset];
649   }
650 };
651 
652 // predefined annotation keys :
653 inline const std::string& key_axis_x_title() {
654   static const std::string s_v("axis_x.title");
655   return s_v;
656 }
657 inline const std::string& key_axis_y_title() {
658   static const std::string s_v("axis_y.title");
659   return s_v;
660 }
661 inline const std::string& key_axis_z_title() {
662   static const std::string s_v("axis_z.title");
663   return s_v;
664 }
665 
666 }}
667 
668 #endif