CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
HistoND.h
Go to the documentation of this file.
1 #ifndef NPSTAT_HISTOND_HH_
2 #define NPSTAT_HISTOND_HH_
3 
14 #include "JetMETCorrections/InterpolationTables/interface/ArrayND.h"
15 #include "JetMETCorrections/InterpolationTables/interface/HistoAxis.h"
16 
17 namespace npstat {
42  template <typename Numeric, class Axis=HistoAxis>
43  class HistoND
44  {
45  template <typename Num2, class Axis2> friend class HistoND;
46 
47  public:
48  typedef Numeric value_type;
49  typedef Axis axis_type;
50 
52  HistoND(const std::vector<Axis>& axes, const char* title=0,
53  const char* accumulatedDataLabel=0);
54 
56  HistoND(const Axis& xAxis, const char* title=0,
57  const char* accumulatedDataLabel=0);
58 
60  HistoND(const Axis& xAxis, const Axis& yAxis,
61  const char* title=0, const char* accumulatedDataLabel=0);
62 
64  HistoND(const Axis& xAxis, const Axis& yAxis, const Axis& zAxis,
65  const char* title=0, const char* accumulatedDataLabel=0);
66 
68  HistoND(const Axis& xAxis, const Axis& yAxis,
69  const Axis& zAxis, const Axis& tAxis,
70  const char* title=0, const char* accumulatedDataLabel=0);
71 
73  HistoND(const Axis& xAxis, const Axis& yAxis,
74  const Axis& zAxis, const Axis& tAxis, const Axis& vAxis,
75  const char* title=0, const char* accumulatedDataLabel=0);
76 
82  HistoND(const ArrayShape& shape, const BoxND<double>& boundingBox,
83  const char* title=0, const char* accumulatedDataLabel=0);
84 
91  template <typename Num2, class Functor>
92  HistoND(const HistoND<Num2,Axis>& h, const Functor& f,
93  const char* title=0, const char* accumulatedDataLabel=0);
94 
102  template <typename Num2>
103  HistoND(const HistoND<Num2,Axis>& h, const unsigned *indices,
104  unsigned nIndices, const char* title=0);
105 
115  template <typename Num2>
116  HistoND(const HistoND<Num2,Axis>& h, const Axis& newAxis,
117  unsigned newAxisNumber, const char* title=0);
118 
120  HistoND(const HistoND&);
121 
126  HistoND& operator=(const HistoND&);
127 
129  inline unsigned dim() const {return dim_;}
130 
132  inline const std::string& title() const {return title_;}
133 
135  inline const std::string& accumulatedDataLabel() const
136  {return accumulatedDataLabel_;}
137 
139  inline const ArrayND<Numeric>& binContents() const {return data_;}
140 
142  inline const ArrayND<Numeric>& overflows() const {return overflow_;}
143 
145  inline const std::vector<Axis>& axes() const {return axes_;}
146 
148  inline const Axis& axis(const unsigned i) const
149  {return axes_.at(i);}
150 
152  inline unsigned long nBins() const {return data_.length();}
153 
155  inline unsigned long nFillsTotal() const {return fillCount_;}
156 
158  inline unsigned long nFillsInRange() const
159  {return fillCount_ - overCount_;}
160 
162  inline unsigned long nFillsOver() const {return overCount_;}
163 
168  bool isUniformlyBinned() const;
169 
171  inline void setTitle(const char* newtitle)
172  {title_ = newtitle ? newtitle : ""; ++modCount_;}
173 
175  inline void setAccumulatedDataLabel(const char* newlabel)
176  {accumulatedDataLabel_ = newlabel ? newlabel : ""; ++modCount_;}
177 
179  inline void setAxisLabel(const unsigned axisNum, const char* newlabel)
180  {axes_.at(axisNum).setLabel(newlabel); ++modCount_;}
181 
186  double binVolume(unsigned long binNumber=0) const;
187 
193  void binCenter(unsigned long binNumber,
194  double* coords, unsigned lenCoords) const;
195 
203  template <class Point>
204  void allBinCenters(std::vector<Point>* centers) const;
205 
207  void binBox(unsigned long binNumber, BoxND<double>* box) const;
208 
210  BoxND<double> boundingBox() const;
211 
217  double volume() const;
218 
220  double integral() const;
221 
223  void clear();
224 
226  void clearBinContents();
227 
229  void clearOverflows();
230 
232  bool operator==(const HistoND&) const;
233 
235  bool operator!=(const HistoND&) const;
236 
241  bool isSameData(const HistoND&) const;
242 
249  template <typename Num2>
250  void fill(const double* coords, unsigned coordLength,
251  const Num2& weight);
252 
254 
258  template <typename Num2>
259  void fill(const Num2& weight);
260 
261  template <typename Num2>
262  void fill(double x0, const Num2& weight);
263 
264  template <typename Num2>
265  void fill(double x0, double x1, const Num2& weight);
266 
267  template <typename Num2>
268  void fill(double x0, double x1, double x2, const Num2& weight);
269 
270  template <typename Num2>
271  void fill(double x0, double x1, double x2, double x3,
272  const Num2& weight);
273 
274  template <typename Num2>
275  void fill(double x0, double x1, double x2, double x3, double x4,
276  const Num2& weight);
277 
278  template <typename Num2>
279  void fill(double x0, double x1, double x2, double x3, double x4,
280  double x5, const Num2& weight);
281 
282  template <typename Num2>
283  void fill(double x0, double x1, double x2, double x3, double x4,
284  double x5, double x6, const Num2& weight);
285 
286  template <typename Num2>
287  void fill(double x0, double x1, double x2, double x3, double x4,
288  double x5, double x6, double x7, const Num2& weight);
289 
290  template <typename Num2>
291  void fill(double x0, double x1, double x2, double x3, double x4,
292  double x5, double x6, double x7, double x8,
293  const Num2& weight);
294 
295  template <typename Num2>
296  void fill(double x0, double x1, double x2, double x3, double x4,
297  double x5, double x6, double x7, double x8, double x9,
298  const Num2& weight);
300 
314  template <typename Num2, class Functor>
315  void dispatch(const double* coords, unsigned coordLength,
316  Num2& weight, Functor& f);
317 
319 
323  template <typename Num2, class Functor>
324  void dispatch(Num2& weight, Functor& f);
325 
326  template <typename Num2, class Functor>
327  void dispatch(double x0, Num2& weight, Functor& f);
328 
329  template <typename Num2, class Functor>
330  void dispatch(double x0, double x1, Num2& weight, Functor& f);
331 
332  template <typename Num2, class Functor>
333  void dispatch(double x0, double x1, double x2, Num2& weight,
334  Functor& f);
335 
336  template <typename Num2, class Functor>
337  void dispatch(double x0, double x1, double x2, double x3,
338  Num2& weight, Functor& f);
339 
340  template <typename Num2, class Functor>
341  void dispatch(double x0, double x1, double x2, double x3, double x4,
342  Num2& weight, Functor& f);
343 
344  template <typename Num2, class Functor>
345  void dispatch(double x0, double x1, double x2, double x3, double x4,
346  double x5, Num2& weight, Functor& f);
347 
348  template <typename Num2, class Functor>
349  void dispatch(double x0, double x1, double x2, double x3, double x4,
350  double x5, double x6, Num2& weight, Functor& f);
351 
352  template <typename Num2, class Functor>
353  void dispatch(double x0, double x1, double x2, double x3, double x4,
354  double x5, double x6, double x7, Num2& weight,
355  Functor& f);
356 
357  template <typename Num2, class Functor>
358  void dispatch(double x0, double x1, double x2, double x3, double x4,
359  double x5, double x6, double x7, double x8,
360  Num2& weight, Functor& f);
361 
362  template <typename Num2, class Functor>
363  void dispatch(double x0, double x1, double x2, double x3, double x4,
364  double x5, double x6, double x7, double x8, double x9,
365  Num2& weight, Functor& f);
367 
374  const Numeric& examine(const double* coords,
375  unsigned coordLength) const;
376 
378 
382  const Numeric& examine() const;
383 
384  const Numeric& examine(double x0) const;
385 
386  const Numeric& examine(double x0, double x1) const;
387 
388  const Numeric& examine(double x0, double x1, double x2) const;
389 
390  const Numeric& examine(double x0, double x1, double x2,
391  double x3) const;
392 
393  const Numeric& examine(double x0, double x1, double x2, double x3,
394  double x4) const;
395 
396  const Numeric& examine(double x0, double x1, double x2, double x3,
397  double x4, double x5) const;
398 
399  const Numeric& examine(double x0, double x1, double x2, double x3,
400  double x4, double x5, double x6) const;
401 
402  const Numeric& examine(double x0, double x1, double x2, double x3,
403  double x4, double x5, double x6,
404  double x7) const;
405 
406  const Numeric& examine(double x0, double x1, double x2, double x3,
407  double x4, double x5, double x6, double x7,
408  double x8) const;
409 
410  const Numeric& examine(double x0, double x1, double x2, double x3,
411  double x4, double x5, double x6, double x7,
412  double x8, double x9) const;
414 
421  const Numeric& closestBin(const double* coords,
422  unsigned coordLength) const;
423 
425 
429  const Numeric& closestBin() const;
430 
431  const Numeric& closestBin(double x0) const;
432 
433  const Numeric& closestBin(double x0, double x1) const;
434 
435  const Numeric& closestBin(double x0, double x1, double x2) const;
436 
437  const Numeric& closestBin(double x0, double x1, double x2,
438  double x3) const;
439 
440  const Numeric& closestBin(double x0, double x1, double x2, double x3,
441  double x4) const;
442 
443  const Numeric& closestBin(double x0, double x1, double x2, double x3,
444  double x4, double x5) const;
445 
446  const Numeric& closestBin(double x0, double x1, double x2, double x3,
447  double x4, double x5, double x6) const;
448 
449  const Numeric& closestBin(double x0, double x1, double x2, double x3,
450  double x4, double x5, double x6,
451  double x7) const;
452 
453  const Numeric& closestBin(double x0, double x1, double x2, double x3,
454  double x4, double x5, double x6, double x7,
455  double x8) const;
456 
457  const Numeric& closestBin(double x0, double x1, double x2, double x3,
458  double x4, double x5, double x6, double x7,
459  double x8, double x9) const;
461 
484  template <typename Num2>
485  void fillC(const double* coords, unsigned coordLength,
486  const Num2& weight);
487 
489 
493  template <typename Num2>
494  void fillC(const Num2& weight);
495 
496  template <typename Num2>
497  void fillC(double x0, const Num2& weight);
498 
499  template <typename Num2>
500  void fillC(double x0, double x1, const Num2& weight);
501 
502  template <typename Num2>
503  void fillC(double x0, double x1, double x2, const Num2& weight);
504 
505  template <typename Num2>
506  void fillC(double x0, double x1, double x2, double x3,
507  const Num2& weight);
508 
509  template <typename Num2>
510  void fillC(double x0, double x1, double x2, double x3, double x4,
511  const Num2& weight);
512 
513  template <typename Num2>
514  void fillC(double x0, double x1, double x2, double x3, double x4,
515  double x5, const Num2& weight);
516 
517  template <typename Num2>
518  void fillC(double x0, double x1, double x2, double x3, double x4,
519  double x5, double x6, const Num2& weight);
520 
521  template <typename Num2>
522  void fillC(double x0, double x1, double x2, double x3, double x4,
523  double x5, double x6, double x7, const Num2& weight);
524 
525  template <typename Num2>
526  void fillC(double x0, double x1, double x2, double x3, double x4,
527  double x5, double x6, double x7, double x8,
528  const Num2& weight);
529 
530  template <typename Num2>
531  void fillC(double x0, double x1, double x2, double x3, double x4,
532  double x5, double x6, double x7, double x8, double x9,
533  const Num2& weight);
535 
540  template <typename Num2>
541  HistoND& operator+=(const HistoND<Num2,Axis>& r);
542 
552  template <typename Num2>
553  HistoND& operator-=(const HistoND<Num2,Axis>& r);
554 
556 
557  template <typename Num2>
558  void setBin(const unsigned *index, unsigned indexLen, const Num2& v);
559 
560  template <typename Num2>
561  void setBin(const Num2& v);
562 
563  template <typename Num2>
564  void setBin(unsigned i0, const Num2& v);
565 
566  template <typename Num2>
567  void setBin(unsigned i0, unsigned i1, const Num2& v);
568 
569  template <typename Num2>
570  void setBin(unsigned i0, unsigned i1, unsigned i2, const Num2& v);
571 
572  template <typename Num2>
573  void setBin(unsigned i0, unsigned i1, unsigned i2, unsigned i3,
574  const Num2& v);
575 
576  template <typename Num2>
577  void setBin(unsigned i0, unsigned i1, unsigned i2, unsigned i3,
578  unsigned i4, const Num2& v);
579 
580  template <typename Num2>
581  void setBin(unsigned i0, unsigned i1, unsigned i2, unsigned i3,
582  unsigned i4, unsigned i5, const Num2& v);
583 
584  template <typename Num2>
585  void setBin(unsigned i0, unsigned i1, unsigned i2, unsigned i3,
586  unsigned i4, unsigned i5, unsigned i6, const Num2& v);
587 
588  template <typename Num2>
589  void setBin(unsigned i0, unsigned i1, unsigned i2, unsigned i3,
590  unsigned i4, unsigned i5, unsigned i6, unsigned i7,
591  const Num2& v);
592 
593  template <typename Num2>
594  void setBin(unsigned i0, unsigned i1, unsigned i2, unsigned i3,
595  unsigned i4, unsigned i5, unsigned i6, unsigned i7,
596  unsigned i8, const Num2& v);
597 
598  template <typename Num2>
599  void setBin(unsigned i0, unsigned i1, unsigned i2, unsigned i3,
600  unsigned i4, unsigned i5, unsigned i6, unsigned i7,
601  unsigned i8, unsigned i9, const Num2& v);
602 
603  template <typename Num2>
604  inline void setLinearBin(const unsigned long index, const Num2& v)
605  {data_.linearValue(index) = v; ++modCount_;}
607 
609 
610  template <typename Num2>
611  void setBinAt(const unsigned *index, unsigned indexLen, const Num2& v);
612 
613  template <typename Num2>
614  void setBinAt(const Num2& v);
615 
616  template <typename Num2>
617  void setBinAt(unsigned i0, const Num2& v);
618 
619  template <typename Num2>
620  void setBinAt(unsigned i0, unsigned i1, const Num2& v);
621 
622  template <typename Num2>
623  void setBinAt(unsigned i0, unsigned i1, unsigned i2, const Num2& v);
624 
625  template <typename Num2>
626  void setBinAt(unsigned i0, unsigned i1, unsigned i2, unsigned i3,
627  const Num2& v);
628 
629  template <typename Num2>
630  void setBinAt(unsigned i0, unsigned i1, unsigned i2, unsigned i3,
631  unsigned i4, const Num2& v);
632 
633  template <typename Num2>
634  void setBinAt(unsigned i0, unsigned i1, unsigned i2, unsigned i3,
635  unsigned i4, unsigned i5, const Num2& v);
636 
637  template <typename Num2>
638  void setBinAt(unsigned i0, unsigned i1, unsigned i2, unsigned i3,
639  unsigned i4, unsigned i5, unsigned i6, const Num2& v);
640 
641  template <typename Num2>
642  void setBinAt(unsigned i0, unsigned i1, unsigned i2, unsigned i3,
643  unsigned i4, unsigned i5, unsigned i6, unsigned i7,
644  const Num2& v);
645 
646  template <typename Num2>
647  void setBinAt(unsigned i0, unsigned i1, unsigned i2, unsigned i3,
648  unsigned i4, unsigned i5, unsigned i6, unsigned i7,
649  unsigned i8, const Num2& v);
650 
651  template <typename Num2>
652  void setBinAt(unsigned i0, unsigned i1, unsigned i2, unsigned i3,
653  unsigned i4, unsigned i5, unsigned i6, unsigned i7,
654  unsigned i8, unsigned i9, const Num2& v);
655 
656  template <typename Num2>
657  inline void setLinearBinAt(const unsigned long index, const Num2& v)
658  {data_.linearValueAt(index) = v; ++modCount_;}
660 
662  template <typename Num2>
663  void setBinContents(const Num2* data, unsigned long dataLength,
664  bool clearOverflows=true);
665 
667  template <typename Num2>
668  void setOverflows(const Num2* data, unsigned long dataLength);
669 
674  template <typename Num2>
675  inline void setBinsToConst(const Num2& value)
676  {data_.constFill(value); ++modCount_;}
677 
682  template <typename Num2>
683  inline void setOverflowsToConst(const Num2& value)
684  {overflow_.constFill(value); ++modCount_;}
685 
692  void recalculateNFillsFromData();
693 
695 
699  inline void setNFillsTotal(const unsigned long i)
700  {fillCount_ = i; ++modCount_;}
701  inline void setNFillsOver(const unsigned long i)
702  {overCount_ = i; ++modCount_;}
704 
706  template <typename Num2>
707  HistoND& operator*=(const Num2& r);
708 
710  template <typename Num2>
711  HistoND& operator/=(const Num2& r);
712 
714 
715  template <typename Num2>
716  void scaleBinContents(const Num2* data, unsigned long dataLength);
717 
718  template <typename Num2>
719  void scaleOverflows(const Num2* data, unsigned long dataLength);
721 
723 
727  template <typename Num2>
728  void addToBinContents(const Num2& weight);
729 
730  template <typename Num2>
731  void addToOverflows(const Num2& weight);
733 
735 
740  template <typename Num2>
741  void addToBinContents(const Num2* data, unsigned long dataLength);
742 
743  template <typename Num2>
744  void addToOverflows(const Num2* data, unsigned long dataLength);
746 
755  template <typename Acc>
756  void accumulateBinsInBox(const BoxND<double>& box, Acc* acc) const;
757 
759 
769  template <typename Num2, typename Num3>
770  void addToProjection(HistoND<Num2,Axis>* projection,
771  AbsArrayProjector<Numeric,Num3>& projector,
772  const unsigned *projectedIndices,
773  unsigned nProjectedIndices) const;
774 
775  template <typename Num2, typename Num3>
776  void addToProjection(HistoND<Num2,Axis>* projection,
777  AbsVisitor<Numeric,Num3>& projector,
778  const unsigned *projectedIndices,
779  unsigned nProjectedIndices) const;
781 
783  HistoND transpose(unsigned axisNum1, unsigned axisNum2) const;
784 
797  inline unsigned long getModCount() const {return modCount_;}
798 
804  inline void incrModCount() {++modCount_;}
805 
807 
808  inline gs::ClassId classId() const {return gs::ClassId(*this);}
809  bool write(std::ostream& of) const;
811 
812  static const char* classname();
813  static inline unsigned version() {return 1;}
814  static HistoND* read(const gs::ClassId& id, std::istream& in);
815 
816  private:
817  HistoND();
818 
819  // Special constructor which speeds up the "transpose" operation.
820  // Does not do full error checking (some of it is done in transpose).
821  HistoND(const HistoND& r, unsigned ax1, unsigned ax2);
822 
823  template <typename Num2>
824  void fillPreservingCentroid(const Num2& weight);
825 
826  template <typename Acc>
827  void accumulateBinsLoop(unsigned level, const BoxND<double>& box,
828  unsigned* idx, Acc* accumulator,
829  double overlapFraction) const;
830  std::string title_;
831  std::string accumulatedDataLabel_;
832  ArrayND<Numeric> data_;
833  ArrayND<Numeric> overflow_;
834  std::vector<Axis> axes_;
835  mutable std::vector<double> weightBuf_;
836  mutable std::vector<unsigned> indexBuf_;
837  unsigned long fillCount_;
838  unsigned long overCount_;
839  unsigned long modCount_;
840  unsigned dim_;
841 
842  };
843 
858  template <typename Histo>
859  void convertHistoToDensity(Histo* histogram, bool knownNonNegative=false);
860 
866  template <typename Histo>
867  std::vector<LinearMapper1d> densityScanHistoMap(const Histo& histo);
868 
879  template <typename Histo>
880  std::vector<CircularMapper1d> convolutionHistoMap(const Histo& histo,
881  bool doubleDataRange);
882 }
883 
884 #include <cassert>
885 #include "JetMETCorrections/InterpolationTables/interface/NpstatException.h"
886 #include <sstream>
887 #include <climits>
888 #include <algorithm>
889 
890 #include "Alignment/Geners/interface/CPP11_auto_ptr.hh"
891 #include "Alignment/Geners/interface/binaryIO.hh"
892 
893 namespace npstat {
894  namespace Private {
895  template <class Axis>
896  ArrayShape makeHistoShape(const std::vector<Axis>& axes)
897  {
898  const unsigned n = axes.size();
899  ArrayShape result;
900  result.reserve(n);
901  for (unsigned i=0; i<n; ++i)
902  result.push_back(axes[i].nBins());
903  return result;
904  }
905 
906  template <class Axis>
907  ArrayShape makeHistoShape(const Axis& xAxis)
908  {
909  ArrayShape result;
910  result.reserve(1U);
911  result.push_back(xAxis.nBins());
912  return result;
913  }
914 
915  template <class Axis>
916  ArrayShape makeHistoShape(const Axis& xAxis, const Axis& yAxis)
917  {
918  ArrayShape result;
919  result.reserve(2U);
920  result.push_back(xAxis.nBins());
921  result.push_back(yAxis.nBins());
922  return result;
923  }
924 
925  template <class Axis>
926  ArrayShape makeHistoShape(const Axis& xAxis,
927  const Axis& yAxis,
928  const Axis& zAxis)
929  {
930  ArrayShape result;
931  result.reserve(3U);
932  result.push_back(xAxis.nBins());
933  result.push_back(yAxis.nBins());
934  result.push_back(zAxis.nBins());
935  return result;
936  }
937 
938  template <class Axis>
939  ArrayShape makeHistoShape(const Axis& xAxis, const Axis& yAxis,
940  const Axis& zAxis, const Axis& tAxis)
941  {
942  ArrayShape result;
943  result.reserve(4U);
944  result.push_back(xAxis.nBins());
945  result.push_back(yAxis.nBins());
946  result.push_back(zAxis.nBins());
947  result.push_back(tAxis.nBins());
948  return result;
949  }
950 
951  template <class Axis>
952  ArrayShape makeHistoShape(const Axis& xAxis, const Axis& yAxis,
953  const Axis& zAxis, const Axis& tAxis,
954  const Axis& vAxis)
955  {
956  ArrayShape result;
957  result.reserve(5U);
958  result.push_back(xAxis.nBins());
959  result.push_back(yAxis.nBins());
960  result.push_back(zAxis.nBins());
961  result.push_back(tAxis.nBins());
962  result.push_back(vAxis.nBins());
963  return result;
964  }
965 
966  template <class Axis>
967  std::vector<Axis> axesOfASlice(const std::vector<Axis>& axes,
968  const unsigned *fixedIndices,
969  const unsigned nFixedIndices)
970  {
971  const unsigned dim = axes.size();
972  std::vector<Axis> newAxes;
973  if (nFixedIndices == 0U) throw npstat::NpstatInvalidArgument(
974  "In npstat::Private::axesOfASlice: "
975  "at least one fixed index must be specified");
976  if (nFixedIndices > dim) throw npstat::NpstatInvalidArgument(
977  "In npstat::Private::axesOfASlice: too many fixed indices");
978  assert(fixedIndices);
979  for (unsigned i=0; i<nFixedIndices; ++i)
980  if (fixedIndices[i] >= dim) throw npstat::NpstatInvalidArgument(
981  "In npstat::Private::axesOfASlice: fixed index out of range");
982  newAxes.reserve(dim - nFixedIndices);
983  for (unsigned i=0; i<dim; ++i)
984  {
985  bool fixed = false;
986  for (unsigned j=0; j<nFixedIndices; ++j)
987  if (fixedIndices[j] == i)
988  {
989  fixed = true;
990  break;
991  }
992  if (!fixed)
993  newAxes.push_back(axes[i]);
994  }
995  if (newAxes.size() != dim - nFixedIndices)
996  throw npstat::NpstatInvalidArgument(
997  "In npstat::Private::axesOfASlice: duplicate fixed index");
998  return newAxes;
999  }
1000 
1001  template <class Axis>
1002  ArrayShape shapeOfASlice(const std::vector<Axis>& axes,
1003  const unsigned *fixedIndices,
1004  const unsigned nFixedIndices)
1005  {
1006  const unsigned dim = axes.size();
1007  if (nFixedIndices == 0U) throw npstat::NpstatInvalidArgument(
1008  "In npstat::Private::shapeOfASlice: "
1009  "at least one fixed index must be specified");
1010  if (nFixedIndices > dim) throw npstat::NpstatInvalidArgument(
1011  "In npstat::Private::shapeOfASlice: too many fixed indices");
1012  assert(fixedIndices);
1013 
1014  // Check that the fixed indices are within range
1015  for (unsigned j=0; j<nFixedIndices; ++j)
1016  if (fixedIndices[j] >= dim) throw npstat::NpstatInvalidArgument(
1017  "In npstat::Private::shapeOfASlice: fixed index out of range");
1018 
1019  // Build the shape for the slice
1020  ArrayShape sh;
1021  if (nFixedIndices < dim)
1022  sh.reserve(dim - nFixedIndices);
1023  for (unsigned i=0; i<dim; ++i)
1024  {
1025  bool fixed = false;
1026  for (unsigned j=0; j<nFixedIndices; ++j)
1027  if (fixedIndices[j] == i)
1028  {
1029  fixed = true;
1030  break;
1031  }
1032  if (!fixed)
1033  sh.push_back(axes[i].nBins());
1034  }
1035  if (sh.size() != dim - nFixedIndices)
1036  throw npstat::NpstatInvalidArgument(
1037  "In npstat::Private::shapeOfASlice: duplicate fixed index");
1038  return sh;
1039  }
1040 
1041  template <class Axis>
1042  std::vector<Axis> addAxis(const std::vector<Axis>& axes,
1043  const Axis& newAxis,
1044  const unsigned newAxisNumber)
1045  {
1046  const unsigned dim = axes.size();
1047  std::vector<Axis> newAxes;
1048  newAxes.reserve(dim + 1U);
1049  unsigned iadd = 0;
1050  for (unsigned i=0; i<dim; ++i)
1051  {
1052  if (newAxisNumber == i)
1053  newAxes.push_back(newAxis);
1054  else
1055  newAxes.push_back(axes[iadd++]);
1056  }
1057  if (iadd == dim)
1058  newAxes.push_back(newAxis);
1059  else
1060  newAxes.push_back(axes[iadd]);
1061  return newAxes;
1062  }
1063 
1064  template <class Axis>
1065  ArrayShape shapeWithExtraAxis(const std::vector<Axis>& axes,
1066  const Axis& newAxis,
1067  const unsigned newAxisNumber)
1068  {
1069  const unsigned dim = axes.size();
1070  ArrayShape result;
1071  result.reserve(dim + 1U);
1072  unsigned iadd = 0;
1073  for (unsigned i=0; i<dim; ++i)
1074  {
1075  if (newAxisNumber == i)
1076  result.push_back(newAxis.nBins());
1077  else
1078  result.push_back(axes[iadd++].nBins());
1079  }
1080  if (iadd == dim)
1081  result.push_back(newAxis.nBins());
1082  else
1083  result.push_back(axes[iadd].nBins());
1084  return result;
1085  }
1086 
1087  inline void h_badargs(const char* method)
1088  {
1089  std::ostringstream os;
1090  os << "In npstat::HistoND::" << method << ": number of arguments"
1091  << " is incompatible with histogram dimensionality";
1092  throw npstat::NpstatInvalidArgument(os.str());
1093  }
1094  }
1095 
1096  template <typename Numeric, class Axis>
1097  template <typename Acc>
1098  void HistoND<Numeric,Axis>::accumulateBinsLoop(
1099  const unsigned level, const BoxND<double>& box,
1100  unsigned* idx, Acc* accumulator, const double overlapFraction) const
1101  {
1102  const Interval<double>& boxSide(box[level]);
1103  const Axis& axis(axes_[level]);
1104  const unsigned nbins = axis.nBins();
1105  const bool lastLevel = level == dim_ - 1U;
1106  for (unsigned i=0; i<nbins; ++i)
1107  {
1108  const double over = overlapFraction*
1109  axis.binInterval(i).overlapFraction(boxSide);
1110  if (over > 0.0)
1111  {
1112  idx[level] = i;
1113  if (lastLevel)
1114  *accumulator += over*data_.value(idx, dim_);
1115  else
1116  accumulateBinsLoop(level+1U, box, idx, accumulator, over);
1117  }
1118  }
1119  }
1120 
1121  template <typename Numeric, class Axis>
1122  template <typename Acc>
1123  void HistoND<Numeric,Axis>::accumulateBinsInBox(
1124  const BoxND<double>& box, Acc* accumulator) const
1125  {
1126  if (box.size() != dim_) throw npstat::NpstatInvalidArgument(
1127  "In npstat::HistoND::accumulateBinsInBox: "
1128  "incompatible box dimensionality");
1129  assert(accumulator);
1130  if (dim_)
1131  {
1132  for (unsigned i=0; i<dim_; ++i)
1133  indexBuf_[i] = 0U;
1134  accumulateBinsLoop(0U, box, &indexBuf_[0], accumulator, 1.0);
1135  }
1136  else
1137  *accumulator += 1.0*data_();
1138  }
1139 
1140  template <typename Numeric, class Axis>
1141  inline void HistoND<Numeric,Axis>::clearBinContents()
1142  {
1143  data_.clear();
1144  fillCount_ = 0UL;
1145  ++modCount_;
1146  }
1147 
1148  template <typename Numeric, class Axis>
1149  inline void HistoND<Numeric,Axis>::clearOverflows()
1150  {
1151  overflow_.clear();
1152  overCount_ = 0UL;
1153  ++modCount_;
1154  }
1155 
1156  template <typename Numeric, class Axis>
1157  inline void HistoND<Numeric,Axis>::clear()
1158  {
1159  clearBinContents();
1160  clearOverflows();
1161  ++modCount_;
1162  }
1163 
1164  template <typename Numeric, class Axis>
1165  HistoND<Numeric,Axis>::HistoND(const std::vector<Axis>& axesIn,
1166  const char* title, const char* label)
1167  : title_(title ? title : ""),
1168  accumulatedDataLabel_(label ? label : ""),
1169  data_(Private::makeHistoShape(axesIn)),
1170  overflow_(ArrayShape(axesIn.size(), 3U)),
1171  axes_(axesIn),
1172  weightBuf_(axesIn.size()),
1173  indexBuf_(2U*axesIn.size()),
1174  modCount_(0UL),
1175  dim_(axesIn.size())
1176  {
1177  if (dim_ >= CHAR_BIT*sizeof(unsigned long))
1178  throw npstat::NpstatInvalidArgument(
1179  "In npstat::HistoND constructor: requested histogram "
1180  "dimensionality is not supported (too large)");
1181  clear();
1182  }
1183 
1184  template <typename Numeric, class Axis>
1185  HistoND<Numeric,Axis>::HistoND(const Axis& xAxis,
1186  const char* title, const char* label)
1187  : title_(title ? title : ""),
1188  accumulatedDataLabel_(label ? label : ""),
1189  data_(Private::makeHistoShape(xAxis)),
1190  overflow_(ArrayShape(1U, 3U)),
1191  weightBuf_(1U),
1192  indexBuf_(2U*1U),
1193  modCount_(0UL),
1194  dim_(1U)
1195  {
1196  axes_.reserve(dim_);
1197  axes_.push_back(xAxis);
1198  clear();
1199  }
1200 
1201  template <typename Numeric, class Axis>
1202  HistoND<Numeric,Axis>::HistoND(const Axis& xAxis, const Axis& yAxis,
1203  const char* title, const char* label)
1204  : title_(title ? title : ""),
1205  accumulatedDataLabel_(label ? label : ""),
1206  data_(Private::makeHistoShape(xAxis, yAxis)),
1207  overflow_(ArrayShape(2U, 3U)),
1208  weightBuf_(2U),
1209  indexBuf_(2U*2U),
1210  modCount_(0UL),
1211  dim_(2U)
1212  {
1213  axes_.reserve(dim_);
1214  axes_.push_back(xAxis);
1215  axes_.push_back(yAxis);
1216  clear();
1217  }
1218 
1219  template <typename Numeric, class Axis>
1220  HistoND<Numeric,Axis>::HistoND(const Axis& xAxis, const Axis& yAxis,
1221  const Axis& zAxis, const char* title,
1222  const char* label)
1223  : title_(title ? title : ""),
1224  accumulatedDataLabel_(label ? label : ""),
1225  data_(Private::makeHistoShape(xAxis, yAxis, zAxis)),
1226  overflow_(ArrayShape(3U, 3U)),
1227  weightBuf_(3U),
1228  indexBuf_(2U*3U),
1229  modCount_(0UL),
1230  dim_(3U)
1231  {
1232  axes_.reserve(dim_);
1233  axes_.push_back(xAxis);
1234  axes_.push_back(yAxis);
1235  axes_.push_back(zAxis);
1236  clear();
1237  }
1238 
1239  template <typename Numeric, class Axis>
1240  HistoND<Numeric,Axis>::HistoND(const Axis& xAxis, const Axis& yAxis,
1241  const Axis& zAxis, const Axis& tAxis,
1242  const char* title, const char* label)
1243  : title_(title ? title : ""),
1244  accumulatedDataLabel_(label ? label : ""),
1245  data_(Private::makeHistoShape(xAxis, yAxis, zAxis, tAxis)),
1246  overflow_(ArrayShape(4U, 3U)),
1247  weightBuf_(4U),
1248  indexBuf_(2U*4U),
1249  modCount_(0UL),
1250  dim_(4U)
1251  {
1252  axes_.reserve(dim_);
1253  axes_.push_back(xAxis);
1254  axes_.push_back(yAxis);
1255  axes_.push_back(zAxis);
1256  axes_.push_back(tAxis);
1257  clear();
1258  }
1259 
1260  template <typename Numeric, class Axis>
1261  HistoND<Numeric,Axis>::HistoND(const Axis& xAxis, const Axis& yAxis,
1262  const Axis& zAxis, const Axis& tAxis,
1263  const Axis& vAxis,
1264  const char* title, const char* label)
1265  : title_(title ? title : ""),
1266  accumulatedDataLabel_(label ? label : ""),
1267  data_(Private::makeHistoShape(xAxis, yAxis, zAxis, tAxis, vAxis)),
1268  overflow_(ArrayShape(5U, 3U)),
1269  weightBuf_(5U),
1270  indexBuf_(2U*5U),
1271  modCount_(0UL),
1272  dim_(5U)
1273  {
1274  axes_.reserve(dim_);
1275  axes_.push_back(xAxis);
1276  axes_.push_back(yAxis);
1277  axes_.push_back(zAxis);
1278  axes_.push_back(tAxis);
1279  axes_.push_back(vAxis);
1280  clear();
1281  }
1282 
1283  template <typename Numeric, class Axis>
1284  HistoND<Numeric,Axis>::HistoND(const ArrayShape& shape,
1285  const BoxND<double>& boundingBox,
1286  const char* title, const char* label)
1287  : title_(title ? title : ""),
1288  accumulatedDataLabel_(label ? label : ""),
1289  data_(shape),
1290  overflow_(ArrayShape(shape.size(), 3U)),
1291  weightBuf_(shape.size()),
1292  indexBuf_(2U*shape.size()),
1293  modCount_(0UL),
1294  dim_(shape.size())
1295  {
1296  if (boundingBox.size() != dim_) throw npstat::NpstatInvalidArgument(
1297  "In npstat::HistoND constructor: "
1298  "incompatible bounding box dimensionality");
1299  if (dim_ >= CHAR_BIT*sizeof(unsigned long))
1300  throw npstat::NpstatInvalidArgument(
1301  "In npstat::HistoND constructor: requested histogram "
1302  "dimensionality is not supported (too large)");
1303  axes_.reserve(dim_);
1304  for (unsigned i=0; i<dim_; ++i)
1305  axes_.push_back(Axis(shape[i],
1306  boundingBox[i].min(),
1307  boundingBox[i].max()));
1308  clear();
1309  }
1310 
1311  template <typename Numeric, class Axis>
1312  template <typename Num2, class Functor>
1313  HistoND<Numeric,Axis>::HistoND(
1314  const HistoND<Num2,Axis>& r, const Functor& f,
1315  const char* title, const char* label)
1316  : title_(title ? title : ""),
1317  accumulatedDataLabel_(label ? label : ""),
1318  data_(r.data_, f),
1319  overflow_(r.overflow_, f),
1320  axes_(r.axes_),
1321  weightBuf_(r.dim_),
1322  indexBuf_(2U*r.dim_),
1323  fillCount_(r.fillCount_),
1324  overCount_(r.overCount_),
1325  modCount_(0UL),
1326  dim_(r.dim_)
1327  {
1328  }
1329 
1330  template <typename Numeric, class Axis>
1331  template <typename Num2>
1332  HistoND<Numeric,Axis>::HistoND(
1333  const HistoND<Num2,Axis>& h, const unsigned *indices,
1334  const unsigned nIndices, const char* title)
1335  : title_(title ? title : ""),
1336  accumulatedDataLabel_(h.accumulatedDataLabel_),
1337  data_(Private::shapeOfASlice(h.axes_, indices, nIndices)),
1338  overflow_(ArrayShape(data_.rank(), 3U)),
1339  axes_(Private::axesOfASlice(h.axes_, indices, nIndices)),
1340  weightBuf_(data_.rank()),
1341  indexBuf_(2U*data_.rank()),
1342  modCount_(0UL),
1343  dim_(data_.rank())
1344  {
1345  clear();
1346  }
1347 
1348  template <typename Numeric, class Axis>
1349  template <typename Num2>
1350  HistoND<Numeric,Axis>::HistoND(
1351  const HistoND<Num2,Axis>& h, const Axis& newAxis,
1352  const unsigned newAxisNumber, const char* title)
1353  : title_(title ? title : ""),
1354  accumulatedDataLabel_(h.accumulatedDataLabel_),
1355  data_(Private::shapeWithExtraAxis(h.axes_, newAxis, newAxisNumber)),
1356  overflow_(data_.rank(), 3U),
1357  axes_(Private::addAxis(h.axes_, newAxis, newAxisNumber)),
1358  weightBuf_(data_.rank()),
1359  indexBuf_(2U*data_.rank()),
1360  modCount_(0UL),
1361  dim_(data_.rank())
1362  {
1363  if (dim_ >= CHAR_BIT*sizeof(unsigned long))
1364  throw npstat::NpstatInvalidArgument(
1365  "In npstat::HistoND constructor: requested histogram "
1366  "dimensionality is not supported (too large)");
1367  clear();
1368  }
1369 
1370  template <typename Numeric, class Axis>
1371  bool HistoND<Numeric,Axis>::isUniformlyBinned() const
1372  {
1373  for (unsigned i=0; i<dim_; ++i)
1374  if (!axes_[i].isUniform())
1375  return false;
1376  return true;
1377  }
1378 
1379  template <typename Numeric, class Axis>
1380  double HistoND<Numeric,Axis>::integral() const
1381  {
1382  typedef typename PreciseType<Numeric>::type Precise;
1383 
1384  if (dim_ == 0U)
1385  return 0.0;
1386  if (isUniformlyBinned())
1387  {
1388  Precise sum = data_.template sum<Precise>();
1389  return static_cast<double>(sum)*binVolume();
1390  }
1391  else
1392  {
1393  Precise sum = Precise();
1394  const Numeric* data = data_.data();
1395  const unsigned long len = data_.length();
1396  for (unsigned long i=0; i<len; ++i)
1397  sum += data[i]*binVolume(i);
1398  return static_cast<double>(sum);
1399  }
1400  }
1401 
1402  template <typename Numeric, class Axis>
1403  BoxND<double> HistoND<Numeric,Axis>::boundingBox() const
1404  {
1405  BoxND<double> box;
1406  if (dim_)
1407  {
1408  box.reserve(dim_);
1409  const Axis* ax = &axes_[0];
1410  for (unsigned i=0; i<dim_; ++i)
1411  box.push_back(ax[i].interval());
1412  }
1413  return box;
1414  }
1415 
1416  template <typename Numeric, class Axis>
1417  void HistoND<Numeric,Axis>::binCenter(
1418  const unsigned long binNumber,
1419  double* coords, const unsigned lenCoords) const
1420  {
1421  if (dim_ != lenCoords) throw npstat::NpstatInvalidArgument(
1422  "In npstat::HistoND::binCenter: "
1423  "incompatible input point dimensionality");
1424  if (dim_)
1425  {
1426  assert(coords);
1427  data_.convertLinearIndex(binNumber, &indexBuf_[0], dim_);
1428  const Axis* ax = &axes_[0];
1429  for (unsigned i=0; i<dim_; ++i)
1430  coords[i] = ax[i].binCenter(indexBuf_[i]);
1431  }
1432  }
1433 
1434  template <typename Numeric, class Axis>
1435  template <class Point>
1436  void HistoND<Numeric,Axis>::allBinCenters(
1437  std::vector<Point>* centers) const
1438  {
1439  assert(centers);
1440  centers->clear();
1441  const unsigned long len = data_.length();
1442  centers->reserve(len);
1443  unsigned* ibuf = &indexBuf_[0];
1444  const Axis* ax = &axes_[0];
1445  Point center;
1446  if (center.size() < dim_) throw npstat::NpstatInvalidArgument(
1447  "In npstat::HistoND::allBinCenters: "
1448  "incompatible point dimensionality (too small)");
1449  typename Point::value_type* cdat = &center[0];
1450 
1451  for (unsigned long i=0; i<len; ++i)
1452  {
1453  data_.convertLinearIndex(i, ibuf, dim_);
1454  for (unsigned idim=0; idim<dim_; ++idim)
1455  cdat[idim] = ax[idim].binCenter(ibuf[idim]);
1456  centers->push_back(center);
1457  }
1458  }
1459 
1460  template <typename Numeric, class Axis>
1461  void HistoND<Numeric,Axis>::binBox(const unsigned long binNumber,
1462  BoxND<double>* box) const
1463  {
1464  assert(box);
1465  box->clear();
1466  if (dim_)
1467  {
1468  box->reserve(dim_);
1469  data_.convertLinearIndex(binNumber, &indexBuf_[0], dim_);
1470  const Axis* ax = &axes_[0];
1471  for (unsigned i=0; i<dim_; ++i)
1472  box->push_back(ax[i].binInterval(indexBuf_[i]));
1473  }
1474  }
1475 
1476  template <typename Numeric, class Axis>
1477  inline bool HistoND<Numeric,Axis>::isSameData(const HistoND& r) const
1478  {
1479  return dim_ == r.dim_ &&
1480  overflow_ == r.overflow_ &&
1481  data_ == r.data_;
1482  }
1483 
1484  template <typename Numeric, class Axis>
1485  inline bool HistoND<Numeric,Axis>::operator==(const HistoND& r) const
1486  {
1487  return dim_ == r.dim_ &&
1488  fillCount_ == r.fillCount_ &&
1489  overCount_ == r.overCount_ &&
1490  title_ == r.title_ &&
1491  accumulatedDataLabel_ == r.accumulatedDataLabel_ &&
1492  axes_ == r.axes_ &&
1493  overflow_ == r.overflow_ &&
1494  data_ == r.data_;
1495  }
1496 
1497  template <typename Numeric, class Axis>
1498  inline bool HistoND<Numeric,Axis>::operator!=(const HistoND& r) const
1499  {
1500  return !(*this == r);
1501  }
1502 
1503  template <typename Numeric, class Axis>
1504  double HistoND<Numeric,Axis>::binVolume(
1505  const unsigned long binNumber) const
1506  {
1507  double v = 1.0;
1508  if (dim_)
1509  {
1510  data_.convertLinearIndex(binNumber, &indexBuf_[0], dim_);
1511  const Axis* ax = &axes_[0];
1512  for (unsigned i=0; i<dim_; ++i)
1513  v *= ax[i].binWidth(indexBuf_[i]);
1514  }
1515  return v;
1516  }
1517 
1518  template <typename Numeric, class Axis>
1519  double HistoND<Numeric,Axis>::volume() const
1520  {
1521  double v = 1.0;
1522  if (dim_)
1523  {
1524  const Axis* ax = &axes_[0];
1525  for (unsigned i=0; i<dim_; ++i)
1526  v *= (ax[i].max() - ax[i].min());
1527  }
1528  return v;
1529  }
1530 
1531  template <typename Numeric, class Axis>
1532  template <typename Num2>
1533  void HistoND<Numeric,Axis>::fill(
1534  const double* coords, const unsigned coordLength, const Num2& w)
1535  {
1536  if (coordLength != dim_) throw npstat::NpstatInvalidArgument(
1537  "In npstat::HistoND::fill: "
1538  "incompatible input point dimensionality");
1539  if (coordLength)
1540  {
1541  assert(coords);
1542  unsigned* idx = &indexBuf_[0];
1543  unsigned* over = idx + dim_;
1544  const Axis* ax = &axes_[0];
1545  unsigned overflown = 0U;
1546  for (unsigned i=0; i<dim_; ++i)
1547  {
1548  over[i] = ax[i].overflowIndex(coords[i], idx + i);
1549  overflown |= (over[i] - 1U);
1550  }
1551  if (overflown)
1552  {
1553  overflow_.value(over, dim_) += w;
1554  ++overCount_;
1555  }
1556  else
1557  data_.value(idx, dim_) += w;
1558  }
1559  else
1560  data_() += w;
1561  ++fillCount_; ++modCount_;
1562  }
1563 
1564  template <typename Numeric, class Axis>
1565  template <typename Num2, class Functor>
1566  void HistoND<Numeric,Axis>::dispatch(
1567  const double* coords, const unsigned coordLength, Num2& w, Functor& f)
1568  {
1569  if (coordLength != dim_) throw npstat::NpstatInvalidArgument(
1570  "In npstat::HistoND::dispatch: "
1571  "incompatible input point dimensionality");
1572  if (coordLength)
1573  {
1574  assert(coords);
1575  unsigned* idx = &indexBuf_[0];
1576  unsigned* over = idx + dim_;
1577  const Axis* ax = &axes_[0];
1578  unsigned overflown = 0U;
1579  for (unsigned i=0; i<dim_; ++i)
1580  {
1581  over[i] = ax[i].overflowIndex(coords[i], idx + i);
1582  overflown |= (over[i] - 1U);
1583  }
1584  if (overflown)
1585  f(overflow_.value(over, dim_), w);
1586  else
1587  f(data_.value(idx, dim_), w);
1588  }
1589  else
1590  f(data_(), w);
1591  ++modCount_;
1592  }
1593 
1594  template <typename Numeric, class Axis>
1595  const Numeric& HistoND<Numeric,Axis>::examine(
1596  const double* coords, const unsigned coordLength) const
1597  {
1598  if (coordLength != dim_) throw npstat::NpstatInvalidArgument(
1599  "In npstat::HistoND::examine: "
1600  "incompatible input point dimensionality");
1601  if (coordLength)
1602  {
1603  assert(coords);
1604  unsigned* idx = &indexBuf_[0];
1605  unsigned* over = idx + dim_;
1606  const Axis* ax = &axes_[0];
1607  unsigned overflown = 0U;
1608  for (unsigned i=0; i<dim_; ++i)
1609  {
1610  over[i] = ax[i].overflowIndex(coords[i], idx + i);
1611  overflown |= (over[i] - 1U);
1612  }
1613  if (overflown)
1614  return overflow_.value(over, dim_);
1615  else
1616  return data_.value(idx, dim_);
1617  }
1618  else
1619  return data_();
1620  }
1621 
1622  template <typename Numeric, class Axis>
1623  const Numeric& HistoND<Numeric,Axis>::closestBin(
1624  const double* coords, const unsigned coordLength) const
1625  {
1626  if (coordLength != dim_) throw npstat::NpstatInvalidArgument(
1627  "In npstat::HistoND::closestBin: "
1628  "incompatible input point dimensionality");
1629  if (coordLength)
1630  {
1631  assert(coords);
1632  unsigned* idx = &indexBuf_[0];
1633  const Axis* ax = &axes_[0];
1634  for (unsigned i=0; i<dim_; ++i)
1635  idx[i] = ax[i].closestValidBin(coords[i]);
1636  return data_.value(idx, dim_);
1637  }
1638  else
1639  return data_();
1640  }
1641 
1642  template <typename Numeric, class Axis>
1643  template <typename Num2>
1644  void HistoND<Numeric,Axis>::fillPreservingCentroid(const Num2& value)
1645  {
1646  const double* weights = &weightBuf_[0];
1647  const unsigned* cell = &indexBuf_[0];
1648  const unsigned long* strides = data_.strides();
1649  const unsigned long maxcycle = 1UL << dim_;
1650  for (unsigned long icycle=0; icycle<maxcycle; ++icycle)
1651  {
1652  double w = 1.0;
1653  unsigned long icell = 0UL;
1654  for (unsigned i=0; i<dim_; ++i)
1655  {
1656  if (icycle & (1UL << i))
1657  {
1658  w *= (1.0 - weights[i]);
1659  icell += strides[i]*(cell[i] + 1U);
1660  }
1661  else
1662  {
1663  w *= weights[i];
1664  icell += strides[i]*cell[i];
1665  }
1666  }
1667  data_.linearValue(icell) += (value * w);
1668  }
1669  }
1670 
1671  template <typename Numeric, class Axis>
1672  template <typename Num2>
1673  void HistoND<Numeric,Axis>::fillC(
1674  const double* coords, const unsigned coordLength, const Num2& w)
1675  {
1676  if (coordLength != dim_) throw npstat::NpstatInvalidArgument(
1677  "In npstat::HistoND::fillC: "
1678  "incompatible input point dimensionality");
1679  if (coordLength)
1680  {
1681  assert(coords);
1682  double* wg = &weightBuf_[0];
1683  unsigned* idx = &indexBuf_[0];
1684  unsigned* over = idx + dim_;
1685  const Axis* ax = &axes_[0];
1686  unsigned overflown = 0U;
1687  for (unsigned i=0; i<dim_; ++i)
1688  {
1689  over[i] = ax[i].overflowIndexWeighted(coords[i], idx+i, wg+i);
1690  overflown |= (over[i] - 1U);
1691  }
1692  if (overflown)
1693  {
1694  overflow_.value(over, dim_) += w;
1695  ++overCount_;
1696  }
1697  else
1698  fillPreservingCentroid(w);
1699  }
1700  else
1701  data_() += w;
1702  ++fillCount_; ++modCount_;
1703  }
1704 
1705  template <typename Numeric, class Axis>
1706  template <typename Num2>
1707  inline void HistoND<Numeric,Axis>::fill(const Num2& w)
1708  {
1709  if (dim_) Private::h_badargs("fill");
1710  data_() += w;
1711  ++fillCount_; ++modCount_;
1712  }
1713 
1714  template <typename Numeric, class Axis>
1715  template <typename Num2, class Functor>
1716  inline void HistoND<Numeric,Axis>::dispatch(Num2& w, Functor& f)
1717  {
1718  if (dim_) Private::h_badargs("dispatch");
1719  f(data_(), w);
1720  ++modCount_;
1721  }
1722 
1723  template <typename Numeric, class Axis>
1724  template <typename Num2>
1725  inline void HistoND<Numeric,Axis>::fillC(const Num2& w)
1726  {
1727  if (dim_) Private::h_badargs("fillC");
1728  data_() += w;
1729  ++fillCount_; ++modCount_;
1730  }
1731 
1732  template <typename Numeric, class Axis>
1733  inline const Numeric& HistoND<Numeric,Axis>::examine() const
1734  {
1735  if (dim_) Private::h_badargs("examine");
1736  return data_();
1737  }
1738 
1739  template <typename Numeric, class Axis>
1740  inline const Numeric& HistoND<Numeric,Axis>::closestBin() const
1741  {
1742  if (dim_) Private::h_badargs("closestBin");
1743  return data_();
1744  }
1745 
1746  template <typename Numeric, class Axis>
1747  template <typename Num2>
1748  void HistoND<Numeric,Axis>::fill(const double x0, const Num2& w)
1749  {
1750  if (dim_ != 1U) Private::h_badargs("fill");
1751  unsigned i0;
1752  const unsigned ov0 = axes_[0].overflowIndex(x0, &i0);
1753  if (ov0 == 1U)
1754  data_(i0) += w;
1755  else
1756  {
1757  overflow_(ov0) += w;
1758  ++overCount_;
1759  }
1760  ++fillCount_; ++modCount_;
1761  }
1762 
1763  template <typename Numeric, class Axis>
1764  template <typename Num2, class Functor>
1765  void HistoND<Numeric,Axis>::dispatch(const double x0, Num2& w, Functor& f)
1766  {
1767  if (dim_ != 1U) Private::h_badargs("dispatch");
1768  unsigned i0;
1769  const unsigned ov0 = axes_[0].overflowIndex(x0, &i0);
1770  if (ov0 == 1U)
1771  f(data_(i0), w);
1772  else
1773  f(overflow_(ov0), w);
1774  ++modCount_;
1775  }
1776 
1777  template <typename Numeric, class Axis>
1778  template <typename Num2>
1779  void HistoND<Numeric,Axis>::fillC(const double x0, const Num2& w)
1780  {
1781  if (dim_ != 1U) Private::h_badargs("fillC");
1782  double* wg = &weightBuf_[0];
1783  unsigned* idx = &indexBuf_[0];
1784  const unsigned ov0 = axes_[0].overflowIndexWeighted(x0, idx, wg);
1785  if (ov0 == 1U)
1786  fillPreservingCentroid(w);
1787  else
1788  {
1789  overflow_(ov0) += w;
1790  ++overCount_;
1791  }
1792  ++fillCount_; ++modCount_;
1793  }
1794 
1795  template <typename Numeric, class Axis>
1796  inline const Numeric& HistoND<Numeric,Axis>::examine(const double x0) const
1797  {
1798  if (dim_ != 1U) Private::h_badargs("examine");
1799  unsigned i0;
1800  const unsigned ov0 = axes_[0].overflowIndex(x0, &i0);
1801  if (ov0 == 1U)
1802  return data_(i0);
1803  else
1804  return overflow_(ov0);
1805  }
1806 
1807  template <typename Numeric, class Axis>
1808  inline const Numeric& HistoND<Numeric,Axis>::closestBin(const double x0) const
1809  {
1810  if (dim_ != 1U) Private::h_badargs("closestBin");
1811  const unsigned i0 = axes_[0].closestValidBin(x0);
1812  return data_(i0);
1813  }
1814 
1815  template <typename Numeric, class Axis>
1816  template <typename Num2>
1817  void HistoND<Numeric,Axis>::fill(const double x0, const double x1,
1818  const Num2& w)
1819  {
1820  if (dim_ != 2U) Private::h_badargs("fill");
1821  unsigned i0, i1;
1822  const Axis* ax = &axes_[0];
1823  const unsigned o0 = ax[0].overflowIndex(x0, &i0);
1824  const unsigned o1 = ax[1].overflowIndex(x1, &i1);
1825  if (o0 == 1U && o1 == 1U)
1826  data_(i0, i1) += w;
1827  else
1828  {
1829  overflow_(o0, o1) += w;
1830  ++overCount_;
1831  }
1832  ++fillCount_; ++modCount_;
1833  }
1834 
1835  template <typename Numeric, class Axis>
1836  template <typename Num2, class Functor>
1837  void HistoND<Numeric,Axis>::dispatch(const double x0, const double x1,
1838  Num2& w, Functor& f)
1839  {
1840  if (dim_ != 2U) Private::h_badargs("dispatch");
1841  unsigned i0, i1;
1842  const Axis* ax = &axes_[0];
1843  const unsigned o0 = ax[0].overflowIndex(x0, &i0);
1844  const unsigned o1 = ax[1].overflowIndex(x1, &i1);
1845  if (o0 == 1U && o1 == 1U)
1846  f(data_(i0, i1), w);
1847  else
1848  f(overflow_(o0, o1), w);
1849  ++modCount_;
1850  }
1851 
1852  template <typename Numeric, class Axis>
1853  template <typename Num2>
1854  void HistoND<Numeric,Axis>::fillC(const double x0, const double x1,
1855  const Num2& w)
1856  {
1857  if (dim_ != 2U) Private::h_badargs("fillC");
1858  double* wg = &weightBuf_[0];
1859  unsigned* idx = &indexBuf_[0];
1860  const Axis* ax = &axes_[0];
1861  const unsigned o0 = ax[0].overflowIndexWeighted(x0, idx+0, wg+0);
1862  const unsigned o1 = ax[1].overflowIndexWeighted(x1, idx+1, wg+1);
1863  if (o0 == 1U && o1 == 1U)
1864  fillPreservingCentroid(w);
1865  else
1866  {
1867  overflow_(o0, o1) += w;
1868  ++overCount_;
1869  }
1870  ++fillCount_; ++modCount_;
1871  }
1872 
1873  template <typename Numeric, class Axis>
1874  const Numeric& HistoND<Numeric,Axis>::examine(const double x0,
1875  const double x1) const
1876  {
1877  if (dim_ != 2U) Private::h_badargs("examine");
1878  unsigned i0, i1;
1879  const Axis* ax = &axes_[0];
1880  const unsigned o0 = ax[0].overflowIndex(x0, &i0);
1881  const unsigned o1 = ax[1].overflowIndex(x1, &i1);
1882  if (o0 == 1U && o1 == 1U)
1883  return data_(i0, i1);
1884  else
1885  return overflow_(o0, o1);
1886  }
1887 
1888  template <typename Numeric, class Axis>
1889  const Numeric& HistoND<Numeric,Axis>::closestBin(const double x0,
1890  const double x1) const
1891  {
1892  if (dim_ != 2U) Private::h_badargs("closestBin");
1893  const Axis* ax = &axes_[0];
1894  const unsigned i0 = ax[0].closestValidBin(x0);
1895  const unsigned i1 = ax[1].closestValidBin(x1);
1896  return data_(i0, i1);
1897  }
1898 
1899  template <typename Numeric, class Axis>
1900  template <typename Num2>
1901  void HistoND<Numeric,Axis>::fill(const double x0, const double x1,
1902  const double x2, const Num2& w)
1903  {
1904  if (dim_ != 3U) Private::h_badargs("fill");
1905  unsigned i0, i1, i2;
1906  const Axis* ax = &axes_[0];
1907  const unsigned o0 = ax[0].overflowIndex(x0, &i0);
1908  const unsigned o1 = ax[1].overflowIndex(x1, &i1);
1909  const unsigned o2 = ax[2].overflowIndex(x2, &i2);
1910  if (o0 == 1U && o1 == 1U && o2 == 1U)
1911  data_(i0, i1, i2) += w;
1912  else
1913  {
1914  overflow_(o0, o1, o2) += w;
1915  ++overCount_;
1916  }
1917  ++fillCount_; ++modCount_;
1918  }
1919 
1920  template <typename Numeric, class Axis>
1921  template <typename Num2, class Functor>
1922  void HistoND<Numeric,Axis>::dispatch(const double x0, const double x1,
1923  const double x2, Num2& w, Functor& f)
1924  {
1925  if (dim_ != 3U) Private::h_badargs("dispatch");
1926  unsigned i0, i1, i2;
1927  const Axis* ax = &axes_[0];
1928  const unsigned o0 = ax[0].overflowIndex(x0, &i0);
1929  const unsigned o1 = ax[1].overflowIndex(x1, &i1);
1930  const unsigned o2 = ax[2].overflowIndex(x2, &i2);
1931  if (o0 == 1U && o1 == 1U && o2 == 1U)
1932  f(data_(i0, i1, i2), w);
1933  else
1934  f(overflow_(o0, o1, o2), w);
1935  ++modCount_;
1936  }
1937 
1938  template <typename Numeric, class Axis>
1939  template <typename Num2>
1940  void HistoND<Numeric,Axis>::fillC(const double x0, const double x1,
1941  const double x2, const Num2& w)
1942  {
1943  if (dim_ != 3U) Private::h_badargs("fillC");
1944  double* wg = &weightBuf_[0];
1945  unsigned* idx = &indexBuf_[0];
1946  const Axis* ax = &axes_[0];
1947  const unsigned o0 = ax[0].overflowIndexWeighted(x0, idx+0, wg+0);
1948  const unsigned o1 = ax[1].overflowIndexWeighted(x1, idx+1, wg+1);
1949  const unsigned o2 = ax[2].overflowIndexWeighted(x2, idx+2, wg+2);
1950  if (o0 == 1U && o1 == 1U && o2 == 1U)
1951  fillPreservingCentroid(w);
1952  else
1953  {
1954  overflow_(o0, o1, o2) += w;
1955  ++overCount_;
1956  }
1957  ++fillCount_; ++modCount_;
1958  }
1959 
1960  template <typename Numeric, class Axis>
1961  const Numeric& HistoND<Numeric,Axis>::examine(const double x0,
1962  const double x1,
1963  const double x2) const
1964  {
1965  if (dim_ != 3U) Private::h_badargs("examine");
1966  unsigned i0, i1, i2;
1967  const Axis* ax = &axes_[0];
1968  const unsigned o0 = ax[0].overflowIndex(x0, &i0);
1969  const unsigned o1 = ax[1].overflowIndex(x1, &i1);
1970  const unsigned o2 = ax[2].overflowIndex(x2, &i2);
1971  if (o0 == 1U && o1 == 1U && o2 == 1U)
1972  return data_(i0, i1, i2);
1973  else
1974  return overflow_(o0, o1, o2);
1975  }
1976 
1977  template <typename Numeric, class Axis>
1978  const Numeric& HistoND<Numeric,Axis>::closestBin(const double x0,
1979  const double x1,
1980  const double x2) const
1981  {
1982  if (dim_ != 3U) Private::h_badargs("closestBin");
1983  const Axis* ax = &axes_[0];
1984  const unsigned i0 = ax[0].closestValidBin(x0);
1985  const unsigned i1 = ax[1].closestValidBin(x1);
1986  const unsigned i2 = ax[2].closestValidBin(x2);
1987  return data_(i0, i1, i2);
1988  }
1989 
1990  template <typename Numeric, class Axis>
1991  template <typename Num2>
1992  void HistoND<Numeric,Axis>::fill(const double x0, const double x1,
1993  const double x2, const double x3,
1994  const Num2& w)
1995  {
1996  if (dim_ != 4U) Private::h_badargs("fill");
1997  unsigned i0, i1, i2, i3;
1998  const Axis* ax = &axes_[0];
1999  const unsigned o0 = ax[0].overflowIndex(x0, &i0);
2000  const unsigned o1 = ax[1].overflowIndex(x1, &i1);
2001  const unsigned o2 = ax[2].overflowIndex(x2, &i2);
2002  const unsigned o3 = ax[3].overflowIndex(x3, &i3);
2003  if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U)
2004  data_(i0, i1, i2, i3) += w;
2005  else
2006  {
2007  overflow_(o0, o1, o2, o3) += w;
2008  ++overCount_;
2009  }
2010  ++fillCount_; ++modCount_;
2011  }
2012 
2013  template <typename Numeric, class Axis>
2014  template <typename Num2, class Functor>
2015  void HistoND<Numeric,Axis>::dispatch(const double x0, const double x1,
2016  const double x2, const double x3,
2017  Num2& w, Functor& f)
2018  {
2019  if (dim_ != 4U) Private::h_badargs("dispatch");
2020  unsigned i0, i1, i2, i3;
2021  const Axis* ax = &axes_[0];
2022  const unsigned o0 = ax[0].overflowIndex(x0, &i0);
2023  const unsigned o1 = ax[1].overflowIndex(x1, &i1);
2024  const unsigned o2 = ax[2].overflowIndex(x2, &i2);
2025  const unsigned o3 = ax[3].overflowIndex(x3, &i3);
2026  if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U)
2027  f(data_(i0, i1, i2, i3), w);
2028  else
2029  f(overflow_(o0, o1, o2, o3), w);
2030  ++modCount_;
2031  }
2032 
2033  template <typename Numeric, class Axis>
2034  template <typename Num2>
2035  void HistoND<Numeric,Axis>::fillC(const double x0, const double x1,
2036  const double x2, const double x3,
2037  const Num2& w)
2038  {
2039  if (dim_ != 4U) Private::h_badargs("fillC");
2040  double* wg = &weightBuf_[0];
2041  unsigned* idx = &indexBuf_[0];
2042  const Axis* ax = &axes_[0];
2043  const unsigned o0 = ax[0].overflowIndexWeighted(x0, idx+0, wg+0);
2044  const unsigned o1 = ax[1].overflowIndexWeighted(x1, idx+1, wg+1);
2045  const unsigned o2 = ax[2].overflowIndexWeighted(x2, idx+2, wg+2);
2046  const unsigned o3 = ax[3].overflowIndexWeighted(x3, idx+3, wg+3);
2047  if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U)
2048  fillPreservingCentroid(w);
2049  else
2050  {
2051  overflow_(o0, o1, o2, o3) += w;
2052  ++overCount_;
2053  }
2054  ++fillCount_; ++modCount_;
2055  }
2056 
2057  template <typename Numeric, class Axis>
2058  const Numeric& HistoND<Numeric,Axis>::examine(const double x0,
2059  const double x1,
2060  const double x2,
2061  const double x3) const
2062  {
2063  if (dim_ != 4U) Private::h_badargs("examine");
2064  unsigned i0, i1, i2, i3;
2065  const Axis* ax = &axes_[0];
2066  const unsigned o0 = ax[0].overflowIndex(x0, &i0);
2067  const unsigned o1 = ax[1].overflowIndex(x1, &i1);
2068  const unsigned o2 = ax[2].overflowIndex(x2, &i2);
2069  const unsigned o3 = ax[3].overflowIndex(x3, &i3);
2070  if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U)
2071  return data_(i0, i1, i2, i3);
2072  else
2073  return overflow_(o0, o1, o2, o3);
2074  }
2075 
2076  template <typename Numeric, class Axis>
2077  const Numeric& HistoND<Numeric,Axis>::closestBin(const double x0,
2078  const double x1,
2079  const double x2,
2080  const double x3) const
2081  {
2082  if (dim_ != 4U) Private::h_badargs("closestBin");
2083  const Axis* ax = &axes_[0];
2084  const unsigned i0 = ax[0].closestValidBin(x0);
2085  const unsigned i1 = ax[1].closestValidBin(x1);
2086  const unsigned i2 = ax[2].closestValidBin(x2);
2087  const unsigned i3 = ax[3].closestValidBin(x3);
2088  return data_(i0, i1, i2, i3);
2089  }
2090 
2091  template <typename Numeric, class Axis>
2092  template <typename Num2>
2093  void HistoND<Numeric,Axis>::fill(const double x0, const double x1,
2094  const double x2, const double x3,
2095  const double x4, const Num2& w)
2096  {
2097  if (dim_ != 5U) Private::h_badargs("fill");
2098  unsigned i0, i1, i2, i3, i4;
2099  const Axis* ax = &axes_[0];
2100  const unsigned o0 = ax[0].overflowIndex(x0, &i0);
2101  const unsigned o1 = ax[1].overflowIndex(x1, &i1);
2102  const unsigned o2 = ax[2].overflowIndex(x2, &i2);
2103  const unsigned o3 = ax[3].overflowIndex(x3, &i3);
2104  const unsigned o4 = ax[4].overflowIndex(x4, &i4);
2105  if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U && o4 == 1U)
2106  data_(i0, i1, i2, i3, i4) += w;
2107  else
2108  {
2109  overflow_(o0, o1, o2, o3, o4) += w;
2110  ++overCount_;
2111  }
2112  ++fillCount_; ++modCount_;
2113  }
2114 
2115  template <typename Numeric, class Axis>
2116  template <typename Num2, class Functor>
2117  void HistoND<Numeric,Axis>::dispatch(const double x0, const double x1,
2118  const double x2, const double x3,
2119  const double x4, Num2& w, Functor& f)
2120  {
2121  if (dim_ != 5U) Private::h_badargs("dispatch");
2122  unsigned i0, i1, i2, i3, i4;
2123  const Axis* ax = &axes_[0];
2124  const unsigned o0 = ax[0].overflowIndex(x0, &i0);
2125  const unsigned o1 = ax[1].overflowIndex(x1, &i1);
2126  const unsigned o2 = ax[2].overflowIndex(x2, &i2);
2127  const unsigned o3 = ax[3].overflowIndex(x3, &i3);
2128  const unsigned o4 = ax[4].overflowIndex(x4, &i4);
2129  if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U && o4 == 1U)
2130  f(data_(i0, i1, i2, i3, i4), w);
2131  else
2132  f(overflow_(o0, o1, o2, o3, o4), w);
2133  ++modCount_;
2134  }
2135 
2136  template <typename Numeric, class Axis>
2137  template <typename Num2>
2138  void HistoND<Numeric,Axis>::fillC(const double x0, const double x1,
2139  const double x2, const double x3,
2140  const double x4, const Num2& w)
2141  {
2142  if (dim_ != 5U) Private::h_badargs("fillC");
2143  double* wg = &weightBuf_[0];
2144  unsigned* idx = &indexBuf_[0];
2145  const Axis* ax = &axes_[0];
2146  const unsigned o0 = ax[0].overflowIndexWeighted(x0, idx+0, wg+0);
2147  const unsigned o1 = ax[1].overflowIndexWeighted(x1, idx+1, wg+1);
2148  const unsigned o2 = ax[2].overflowIndexWeighted(x2, idx+2, wg+2);
2149  const unsigned o3 = ax[3].overflowIndexWeighted(x3, idx+3, wg+3);
2150  const unsigned o4 = ax[4].overflowIndexWeighted(x4, idx+4, wg+4);
2151  if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U && o4 == 1U)
2152  fillPreservingCentroid(w);
2153  else
2154  {
2155  overflow_(o0, o1, o2, o3, o4) += w;
2156  ++overCount_;
2157  }
2158  ++fillCount_; ++modCount_;
2159  }
2160 
2161  template <typename Numeric, class Axis>
2162  const Numeric& HistoND<Numeric,Axis>::examine(
2163  const double x0, const double x1,
2164  const double x2, const double x3,
2165  const double x4) const
2166  {
2167  if (dim_ != 5U) Private::h_badargs("examine");
2168  unsigned i0, i1, i2, i3, i4;
2169  const Axis* ax = &axes_[0];
2170  const unsigned o0 = ax[0].overflowIndex(x0, &i0);
2171  const unsigned o1 = ax[1].overflowIndex(x1, &i1);
2172  const unsigned o2 = ax[2].overflowIndex(x2, &i2);
2173  const unsigned o3 = ax[3].overflowIndex(x3, &i3);
2174  const unsigned o4 = ax[4].overflowIndex(x4, &i4);
2175  if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U && o4 == 1U)
2176  return data_(i0, i1, i2, i3, i4);
2177  else
2178  return overflow_(o0, o1, o2, o3, o4);
2179  }
2180 
2181  template <typename Numeric, class Axis>
2182  const Numeric& HistoND<Numeric,Axis>::closestBin(const double x0,
2183  const double x1,
2184  const double x2,
2185  const double x3,
2186  const double x4) const
2187  {
2188  if (dim_ != 5U) Private::h_badargs("closestBin");
2189  const Axis* ax = &axes_[0];
2190  const unsigned i0 = ax[0].closestValidBin(x0);
2191  const unsigned i1 = ax[1].closestValidBin(x1);
2192  const unsigned i2 = ax[2].closestValidBin(x2);
2193  const unsigned i3 = ax[3].closestValidBin(x3);
2194  const unsigned i4 = ax[4].closestValidBin(x4);
2195  return data_(i0, i1, i2, i3, i4);
2196  }
2197 
2198  template <typename Numeric, class Axis>
2199  template <typename Num2>
2200  void HistoND<Numeric,Axis>::fill(const double x0, const double x1,
2201  const double x2, const double x3,
2202  const double x4, const double x5,
2203  const Num2& w)
2204  {
2205  if (dim_ != 6U) Private::h_badargs("fill");
2206  unsigned i0, i1, i2, i3, i4, i5;
2207  const Axis* ax = &axes_[0];
2208  const unsigned o0 = ax[0].overflowIndex(x0, &i0);
2209  const unsigned o1 = ax[1].overflowIndex(x1, &i1);
2210  const unsigned o2 = ax[2].overflowIndex(x2, &i2);
2211  const unsigned o3 = ax[3].overflowIndex(x3, &i3);
2212  const unsigned o4 = ax[4].overflowIndex(x4, &i4);
2213  const unsigned o5 = ax[5].overflowIndex(x5, &i5);
2214  if (o0 == 1U && o1 == 1U && o2 == 1U &&
2215  o3 == 1U && o4 == 1U && o5 == 1U)
2216  data_(i0, i1, i2, i3, i4, i5) += w;
2217  else
2218  {
2219  overflow_(o0, o1, o2, o3, o4, o5) += w;
2220  ++overCount_;
2221  }
2222  ++fillCount_; ++modCount_;
2223  }
2224 
2225  template <typename Numeric, class Axis>
2226  template <typename Num2, class Functor>
2227  void HistoND<Numeric,Axis>::dispatch(const double x0, const double x1,
2228  const double x2, const double x3,
2229  const double x4, const double x5,
2230  Num2& w, Functor& f)
2231  {
2232  if (dim_ != 6U) Private::h_badargs("dispatch");
2233  unsigned i0, i1, i2, i3, i4, i5;
2234  const Axis* ax = &axes_[0];
2235  const unsigned o0 = ax[0].overflowIndex(x0, &i0);
2236  const unsigned o1 = ax[1].overflowIndex(x1, &i1);
2237  const unsigned o2 = ax[2].overflowIndex(x2, &i2);
2238  const unsigned o3 = ax[3].overflowIndex(x3, &i3);
2239  const unsigned o4 = ax[4].overflowIndex(x4, &i4);
2240  const unsigned o5 = ax[5].overflowIndex(x5, &i5);
2241  if (o0 == 1U && o1 == 1U && o2 == 1U &&
2242  o3 == 1U && o4 == 1U && o5 == 1U)
2243  f(data_(i0, i1, i2, i3, i4, i5), w);
2244  else
2245  f(overflow_(o0, o1, o2, o3, o4, o5), w);
2246  ++modCount_;
2247  }
2248 
2249  template <typename Numeric, class Axis>
2250  template <typename Num2>
2251  void HistoND<Numeric,Axis>::fillC(const double x0, const double x1,
2252  const double x2, const double x3,
2253  const double x4, const double x5,
2254  const Num2& w)
2255  {
2256  if (dim_ != 6U) Private::h_badargs("fillC");
2257  double* wg = &weightBuf_[0];
2258  unsigned* idx = &indexBuf_[0];
2259  const Axis* ax = &axes_[0];
2260  const unsigned o0 = ax[0].overflowIndexWeighted(x0, idx+0, wg+0);
2261  const unsigned o1 = ax[1].overflowIndexWeighted(x1, idx+1, wg+1);
2262  const unsigned o2 = ax[2].overflowIndexWeighted(x2, idx+2, wg+2);
2263  const unsigned o3 = ax[3].overflowIndexWeighted(x3, idx+3, wg+3);
2264  const unsigned o4 = ax[4].overflowIndexWeighted(x4, idx+4, wg+4);
2265  const unsigned o5 = ax[5].overflowIndexWeighted(x5, idx+5, wg+5);
2266  if (o0 == 1U && o1 == 1U && o2 == 1U &&
2267  o3 == 1U && o4 == 1U && o5 == 1U)
2268  fillPreservingCentroid(w);
2269  else
2270  {
2271  overflow_(o0, o1, o2, o3, o4, o5) += w;
2272  ++overCount_;
2273  }
2274  ++fillCount_; ++modCount_;
2275  }
2276 
2277  template <typename Numeric, class Axis>
2278  const Numeric& HistoND<Numeric,Axis>::examine(const double x0,
2279  const double x1,
2280  const double x2,
2281  const double x3,
2282  const double x4,
2283  const double x5) const
2284  {
2285  if (dim_ != 6U) Private::h_badargs("examine");
2286  unsigned i0, i1, i2, i3, i4, i5;
2287  const Axis* ax = &axes_[0];
2288  const unsigned o0 = ax[0].overflowIndex(x0, &i0);
2289  const unsigned o1 = ax[1].overflowIndex(x1, &i1);
2290  const unsigned o2 = ax[2].overflowIndex(x2, &i2);
2291  const unsigned o3 = ax[3].overflowIndex(x3, &i3);
2292  const unsigned o4 = ax[4].overflowIndex(x4, &i4);
2293  const unsigned o5 = ax[5].overflowIndex(x5, &i5);
2294  if (o0 == 1U && o1 == 1U && o2 == 1U &&
2295  o3 == 1U && o4 == 1U && o5 == 1U)
2296  return data_(i0, i1, i2, i3, i4, i5);
2297  else
2298  return overflow_(o0, o1, o2, o3, o4, o5);
2299  }
2300 
2301  template <typename Numeric, class Axis>
2302  const Numeric& HistoND<Numeric,Axis>::closestBin(const double x0,
2303  const double x1,
2304  const double x2,
2305  const double x3,
2306  const double x4,
2307  const double x5) const
2308  {
2309  if (dim_ != 6U) Private::h_badargs("closestBin");
2310  const Axis* ax = &axes_[0];
2311  const unsigned i0 = ax[0].closestValidBin(x0);
2312  const unsigned i1 = ax[1].closestValidBin(x1);
2313  const unsigned i2 = ax[2].closestValidBin(x2);
2314  const unsigned i3 = ax[3].closestValidBin(x3);
2315  const unsigned i4 = ax[4].closestValidBin(x4);
2316  const unsigned i5 = ax[5].closestValidBin(x5);
2317  return data_(i0, i1, i2, i3, i4, i5);
2318  }
2319 
2320  template <typename Numeric, class Axis>
2321  template <typename Num2>
2322  void HistoND<Numeric,Axis>::fill(const double x0, const double x1,
2323  const double x2, const double x3,
2324  const double x4, const double x5,
2325  const double x6, const Num2& w)
2326  {
2327  if (dim_ != 7U) Private::h_badargs("fill");
2328  unsigned i0, i1, i2, i3, i4, i5, i6;
2329  const Axis* ax = &axes_[0];
2330  const unsigned o0 = ax[0].overflowIndex(x0, &i0);
2331  const unsigned o1 = ax[1].overflowIndex(x1, &i1);
2332  const unsigned o2 = ax[2].overflowIndex(x2, &i2);
2333  const unsigned o3 = ax[3].overflowIndex(x3, &i3);
2334  const unsigned o4 = ax[4].overflowIndex(x4, &i4);
2335  const unsigned o5 = ax[5].overflowIndex(x5, &i5);
2336  const unsigned o6 = ax[6].overflowIndex(x6, &i6);
2337  if (o0 == 1U && o1 == 1U && o2 == 1U &&
2338  o3 == 1U && o4 == 1U && o5 == 1U && o6 == 1U)
2339  data_(i0, i1, i2, i3, i4, i5, i6) += w;
2340  else
2341  {
2342  overflow_(o0, o1, o2, o3, o4, o5, o6) += w;
2343  ++overCount_;
2344  }
2345  ++fillCount_; ++modCount_;
2346  }
2347 
2348  template <typename Numeric, class Axis>
2349  template <typename Num2, class Functor>
2350  void HistoND<Numeric,Axis>::dispatch(const double x0, const double x1,
2351  const double x2, const double x3,
2352  const double x4, const double x5,
2353  const double x6, Num2& w, Functor& f)
2354  {
2355  if (dim_ != 7U) Private::h_badargs("dispatch");
2356  unsigned i0, i1, i2, i3, i4, i5, i6;
2357  const Axis* ax = &axes_[0];
2358  const unsigned o0 = ax[0].overflowIndex(x0, &i0);
2359  const unsigned o1 = ax[1].overflowIndex(x1, &i1);
2360  const unsigned o2 = ax[2].overflowIndex(x2, &i2);
2361  const unsigned o3 = ax[3].overflowIndex(x3, &i3);
2362  const unsigned o4 = ax[4].overflowIndex(x4, &i4);
2363  const unsigned o5 = ax[5].overflowIndex(x5, &i5);
2364  const unsigned o6 = ax[6].overflowIndex(x6, &i6);
2365  if (o0 == 1U && o1 == 1U && o2 == 1U &&
2366  o3 == 1U && o4 == 1U && o5 == 1U && o6 == 1U)
2367  f(data_(i0, i1, i2, i3, i4, i5, i6), w);
2368  else
2369  f(overflow_(o0, o1, o2, o3, o4, o5, o6), w);
2370  ++modCount_;
2371  }
2372 
2373  template <typename Numeric, class Axis>
2374  template <typename Num2>
2375  void HistoND<Numeric,Axis>::fillC(const double x0, const double x1,
2376  const double x2, const double x3,
2377  const double x4, const double x5,
2378  const double x6, const Num2& w)
2379  {
2380  if (dim_ != 7U) Private::h_badargs("fillC");
2381  double* wg = &weightBuf_[0];
2382  unsigned* idx = &indexBuf_[0];
2383  const Axis* ax = &axes_[0];
2384  const unsigned o0 = ax[0].overflowIndexWeighted(x0, idx+0, wg+0);
2385  const unsigned o1 = ax[1].overflowIndexWeighted(x1, idx+1, wg+1);
2386  const unsigned o2 = ax[2].overflowIndexWeighted(x2, idx+2, wg+2);
2387  const unsigned o3 = ax[3].overflowIndexWeighted(x3, idx+3, wg+3);
2388  const unsigned o4 = ax[4].overflowIndexWeighted(x4, idx+4, wg+4);
2389  const unsigned o5 = ax[5].overflowIndexWeighted(x5, idx+5, wg+5);
2390  const unsigned o6 = ax[6].overflowIndexWeighted(x6, idx+6, wg+6);
2391  if (o0 == 1U && o1 == 1U && o2 == 1U &&
2392  o3 == 1U && o4 == 1U && o5 == 1U && o6 == 1U)
2393  fillPreservingCentroid(w);
2394  else
2395  {
2396  overflow_(o0, o1, o2, o3, o4, o5, o6) += w;
2397  ++overCount_;
2398  }
2399  ++fillCount_; ++modCount_;
2400  }
2401 
2402  template <typename Numeric, class Axis>
2403  const Numeric& HistoND<Numeric,Axis>::examine(
2404  const double x0, const double x1,
2405  const double x2, const double x3,
2406  const double x4, const double x5,
2407  const double x6) const
2408  {
2409  if (dim_ != 7U) Private::h_badargs("examine");
2410  unsigned i0, i1, i2, i3, i4, i5, i6;
2411  const Axis* ax = &axes_[0];
2412  const unsigned o0 = ax[0].overflowIndex(x0, &i0);
2413  const unsigned o1 = ax[1].overflowIndex(x1, &i1);
2414  const unsigned o2 = ax[2].overflowIndex(x2, &i2);
2415  const unsigned o3 = ax[3].overflowIndex(x3, &i3);
2416  const unsigned o4 = ax[4].overflowIndex(x4, &i4);
2417  const unsigned o5 = ax[5].overflowIndex(x5, &i5);
2418  const unsigned o6 = ax[6].overflowIndex(x6, &i6);
2419  if (o0 == 1U && o1 == 1U && o2 == 1U &&
2420  o3 == 1U && o4 == 1U && o5 == 1U && o6 == 1U)
2421  return data_(i0, i1, i2, i3, i4, i5, i6);
2422  else
2423  return overflow_(o0, o1, o2, o3, o4, o5, o6);
2424  }
2425 
2426  template <typename Numeric, class Axis>
2427  const Numeric& HistoND<Numeric,Axis>::closestBin(const double x0,
2428  const double x1,
2429  const double x2,
2430  const double x3,
2431  const double x4,
2432  const double x5,
2433  const double x6) const
2434  {
2435  if (dim_ != 7U) Private::h_badargs("closestBin");
2436  const Axis* ax = &axes_[0];
2437  const unsigned i0 = ax[0].closestValidBin(x0);
2438  const unsigned i1 = ax[1].closestValidBin(x1);
2439  const unsigned i2 = ax[2].closestValidBin(x2);
2440  const unsigned i3 = ax[3].closestValidBin(x3);
2441  const unsigned i4 = ax[4].closestValidBin(x4);
2442  const unsigned i5 = ax[5].closestValidBin(x5);
2443  const unsigned i6 = ax[6].closestValidBin(x6);
2444  return data_(i0, i1, i2, i3, i4, i5, i6);
2445  }
2446 
2447  template <typename Numeric, class Axis>
2448  template <typename Num2>
2449  void HistoND<Numeric,Axis>::fill(const double x0, const double x1,
2450  const double x2, const double x3,
2451  const double x4, const double x5,
2452  const double x6, const double x7,
2453  const Num2& w)
2454  {
2455  if (dim_ != 8U) Private::h_badargs("fill");
2456  unsigned i0, i1, i2, i3, i4, i5, i6, i7;
2457  const Axis* ax = &axes_[0];
2458  const unsigned o0 = ax[0].overflowIndex(x0, &i0);
2459  const unsigned o1 = ax[1].overflowIndex(x1, &i1);
2460  const unsigned o2 = ax[2].overflowIndex(x2, &i2);
2461  const unsigned o3 = ax[3].overflowIndex(x3, &i3);
2462  const unsigned o4 = ax[4].overflowIndex(x4, &i4);
2463  const unsigned o5 = ax[5].overflowIndex(x5, &i5);
2464  const unsigned o6 = ax[6].overflowIndex(x6, &i6);
2465  const unsigned o7 = ax[7].overflowIndex(x7, &i7);
2466  if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U &&
2467  o4 == 1U && o5 == 1U && o6 == 1U && o7 == 1U)
2468  data_(i0, i1, i2, i3, i4, i5, i6, i7) += w;
2469  else
2470  {
2471  overflow_(o0, o1, o2, o3, o4, o5, o6, o7) += w;
2472  ++overCount_;
2473  }
2474  ++fillCount_; ++modCount_;
2475  }
2476 
2477  template <typename Numeric, class Axis>
2478  template <typename Num2, class Functor>
2479  void HistoND<Numeric,Axis>::dispatch(const double x0, const double x1,
2480  const double x2, const double x3,
2481  const double x4, const double x5,
2482  const double x6, const double x7,
2483  Num2& w, Functor& f)
2484  {
2485  if (dim_ != 8U) Private::h_badargs("dispatch");
2486  unsigned i0, i1, i2, i3, i4, i5, i6, i7;
2487  const Axis* ax = &axes_[0];
2488  const unsigned o0 = ax[0].overflowIndex(x0, &i0);
2489  const unsigned o1 = ax[1].overflowIndex(x1, &i1);
2490  const unsigned o2 = ax[2].overflowIndex(x2, &i2);
2491  const unsigned o3 = ax[3].overflowIndex(x3, &i3);
2492  const unsigned o4 = ax[4].overflowIndex(x4, &i4);
2493  const unsigned o5 = ax[5].overflowIndex(x5, &i5);
2494  const unsigned o6 = ax[6].overflowIndex(x6, &i6);
2495  const unsigned o7 = ax[7].overflowIndex(x7, &i7);
2496  if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U &&
2497  o4 == 1U && o5 == 1U && o6 == 1U && o7 == 1U)
2498  f(data_(i0, i1, i2, i3, i4, i5, i6, i7), w);
2499  else
2500  f(overflow_(o0, o1, o2, o3, o4, o5, o6, o7), w);
2501  ++modCount_;
2502  }
2503 
2504  template <typename Numeric, class Axis>
2505  template <typename Num2>
2506  void HistoND<Numeric,Axis>::fillC(const double x0, const double x1,
2507  const double x2, const double x3,
2508  const double x4, const double x5,
2509  const double x6, const double x7,
2510  const Num2& w)
2511  {
2512  if (dim_ != 8U) Private::h_badargs("fillC");
2513  double* wg = &weightBuf_[0];
2514  unsigned* idx = &indexBuf_[0];
2515  const Axis* ax = &axes_[0];
2516  const unsigned o0 = ax[0].overflowIndexWeighted(x0, idx+0, wg+0);
2517  const unsigned o1 = ax[1].overflowIndexWeighted(x1, idx+1, wg+1);
2518  const unsigned o2 = ax[2].overflowIndexWeighted(x2, idx+2, wg+2);
2519  const unsigned o3 = ax[3].overflowIndexWeighted(x3, idx+3, wg+3);
2520  const unsigned o4 = ax[4].overflowIndexWeighted(x4, idx+4, wg+4);
2521  const unsigned o5 = ax[5].overflowIndexWeighted(x5, idx+5, wg+5);
2522  const unsigned o6 = ax[6].overflowIndexWeighted(x6, idx+6, wg+6);
2523  const unsigned o7 = ax[7].overflowIndexWeighted(x7, idx+7, wg+7);
2524  if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U &&
2525  o4 == 1U && o5 == 1U && o6 == 1U && o7 == 1U)
2526  fillPreservingCentroid(w);
2527  else
2528  {
2529  overflow_(o0, o1, o2, o3, o4, o5, o6, o7) += w;
2530  ++overCount_;
2531  }
2532  ++fillCount_; ++modCount_;
2533  }
2534 
2535  template <typename Numeric, class Axis>
2536  const Numeric& HistoND<Numeric,Axis>::examine(
2537  const double x0, const double x1,
2538  const double x2, const double x3,
2539  const double x4, const double x5,
2540  const double x6,
2541  const double x7) const
2542  {
2543  if (dim_ != 8U) Private::h_badargs("examine");
2544  unsigned i0, i1, i2, i3, i4, i5, i6, i7;
2545  const Axis* ax = &axes_[0];
2546  const unsigned o0 = ax[0].overflowIndex(x0, &i0);
2547  const unsigned o1 = ax[1].overflowIndex(x1, &i1);
2548  const unsigned o2 = ax[2].overflowIndex(x2, &i2);
2549  const unsigned o3 = ax[3].overflowIndex(x3, &i3);
2550  const unsigned o4 = ax[4].overflowIndex(x4, &i4);
2551  const unsigned o5 = ax[5].overflowIndex(x5, &i5);
2552  const unsigned o6 = ax[6].overflowIndex(x6, &i6);
2553  const unsigned o7 = ax[7].overflowIndex(x7, &i7);
2554  if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U &&
2555  o4 == 1U && o5 == 1U && o6 == 1U && o7 == 1U)
2556  return data_(i0, i1, i2, i3, i4, i5, i6, i7);
2557  else
2558  return overflow_(o0, o1, o2, o3, o4, o5, o6, o7);
2559  }
2560 
2561  template <typename Numeric, class Axis>
2562  const Numeric& HistoND<Numeric,Axis>::closestBin(const double x0,
2563  const double x1,
2564  const double x2,
2565  const double x3,
2566  const double x4,
2567  const double x5,
2568  const double x6,
2569  const double x7) const
2570  {
2571  if (dim_ != 8U) Private::h_badargs("closestBin");
2572  const Axis* ax = &axes_[0];
2573  const unsigned i0 = ax[0].closestValidBin(x0);
2574  const unsigned i1 = ax[1].closestValidBin(x1);
2575  const unsigned i2 = ax[2].closestValidBin(x2);
2576  const unsigned i3 = ax[3].closestValidBin(x3);
2577  const unsigned i4 = ax[4].closestValidBin(x4);
2578  const unsigned i5 = ax[5].closestValidBin(x5);
2579  const unsigned i6 = ax[6].closestValidBin(x6);
2580  const unsigned i7 = ax[7].closestValidBin(x7);
2581  return data_(i0, i1, i2, i3, i4, i5, i6, i7);
2582  }
2583 
2584  template <typename Numeric, class Axis>
2585  template <typename Num2>
2586  void HistoND<Numeric,Axis>::fill(const double x0, const double x1,
2587  const double x2, const double x3,
2588  const double x4, const double x5,
2589  const double x6, const double x7,
2590  const double x8, const Num2& w)
2591  {
2592  if (dim_ != 9U) Private::h_badargs("fill");
2593  unsigned i0, i1, i2, i3, i4, i5, i6, i7, i8;
2594  const Axis* ax = &axes_[0];
2595  const unsigned o0 = ax[0].overflowIndex(x0, &i0);
2596  const unsigned o1 = ax[1].overflowIndex(x1, &i1);
2597  const unsigned o2 = ax[2].overflowIndex(x2, &i2);
2598  const unsigned o3 = ax[3].overflowIndex(x3, &i3);
2599  const unsigned o4 = ax[4].overflowIndex(x4, &i4);
2600  const unsigned o5 = ax[5].overflowIndex(x5, &i5);
2601  const unsigned o6 = ax[6].overflowIndex(x6, &i6);
2602  const unsigned o7 = ax[7].overflowIndex(x7, &i7);
2603  const unsigned o8 = ax[8].overflowIndex(x8, &i8);
2604  if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U && o4 == 1U &&
2605  o5 == 1U && o6 == 1U && o7 == 1U && o8 == 1U)
2606  data_(i0, i1, i2, i3, i4, i5, i6, i7, i8) += w;
2607  else
2608  {
2609  overflow_(o0, o1, o2, o3, o4, o5, o6, o7, o8) += w;
2610  ++overCount_;
2611  }
2612  ++fillCount_; ++modCount_;
2613  }
2614 
2615  template <typename Numeric, class Axis>
2616  template <typename Num2, class Functor>
2617  void HistoND<Numeric,Axis>::dispatch(const double x0, const double x1,
2618  const double x2, const double x3,
2619  const double x4, const double x5,
2620  const double x6, const double x7,
2621  const double x8, Num2& w, Functor& f)
2622  {
2623  if (dim_ != 9U) Private::h_badargs("dispatch");
2624  unsigned i0, i1, i2, i3, i4, i5, i6, i7, i8;
2625  const Axis* ax = &axes_[0];
2626  const unsigned o0 = ax[0].overflowIndex(x0, &i0);
2627  const unsigned o1 = ax[1].overflowIndex(x1, &i1);
2628  const unsigned o2 = ax[2].overflowIndex(x2, &i2);
2629  const unsigned o3 = ax[3].overflowIndex(x3, &i3);
2630  const unsigned o4 = ax[4].overflowIndex(x4, &i4);
2631  const unsigned o5 = ax[5].overflowIndex(x5, &i5);
2632  const unsigned o6 = ax[6].overflowIndex(x6, &i6);
2633  const unsigned o7 = ax[7].overflowIndex(x7, &i7);
2634  const unsigned o8 = ax[8].overflowIndex(x8, &i8);
2635  if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U && o4 == 1U &&
2636  o5 == 1U && o6 == 1U && o7 == 1U && o8 == 1U)
2637  f(data_(i0, i1, i2, i3, i4, i5, i6, i7, i8), w);
2638  else
2639  f(overflow_(o0, o1, o2, o3, o4, o5, o6, o7, o8), w);
2640  ++modCount_;
2641  }
2642 
2643  template <typename Numeric, class Axis>
2644  template <typename Num2>
2645  void HistoND<Numeric,Axis>::fillC(const double x0, const double x1,
2646  const double x2, const double x3,
2647  const double x4, const double x5,
2648  const double x6, const double x7,
2649  const double x8, const Num2& w)
2650  {
2651  if (dim_ != 9U) Private::h_badargs("fillC");
2652  double* wg = &weightBuf_[0];
2653  unsigned* idx = &indexBuf_[0];
2654  const Axis* ax = &axes_[0];
2655  const unsigned o0 = ax[0].overflowIndexWeighted(x0, idx+0, wg+0);
2656  const unsigned o1 = ax[1].overflowIndexWeighted(x1, idx+1, wg+1);
2657  const unsigned o2 = ax[2].overflowIndexWeighted(x2, idx+2, wg+2);
2658  const unsigned o3 = ax[3].overflowIndexWeighted(x3, idx+3, wg+3);
2659  const unsigned o4 = ax[4].overflowIndexWeighted(x4, idx+4, wg+4);
2660  const unsigned o5 = ax[5].overflowIndexWeighted(x5, idx+5, wg+5);
2661  const unsigned o6 = ax[6].overflowIndexWeighted(x6, idx+6, wg+6);
2662  const unsigned o7 = ax[7].overflowIndexWeighted(x7, idx+7, wg+7);
2663  const unsigned o8 = ax[8].overflowIndexWeighted(x8, idx+8, wg+8);
2664  if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U && o4 == 1U &&
2665  o5 == 1U && o6 == 1U && o7 == 1U && o8 == 1U)
2666  fillPreservingCentroid(w);
2667  else
2668  {
2669  overflow_(o0, o1, o2, o3, o4, o5, o6, o7, o8) += w;
2670  ++overCount_;
2671  }
2672  ++fillCount_; ++modCount_;
2673  }
2674 
2675  template <typename Numeric, class Axis>
2676  const Numeric& HistoND<Numeric,Axis>::examine(
2677  const double x0, const double x1,
2678  const double x2, const double x3,
2679  const double x4, const double x5,
2680  const double x6, const double x7,
2681  const double x8) const
2682  {
2683  if (dim_ != 9U) Private::h_badargs("examine");
2684  unsigned i0, i1, i2, i3, i4, i5, i6, i7, i8;
2685  const Axis* ax = &axes_[0];
2686  const unsigned o0 = ax[0].overflowIndex(x0, &i0);
2687  const unsigned o1 = ax[1].overflowIndex(x1, &i1);
2688  const unsigned o2 = ax[2].overflowIndex(x2, &i2);
2689  const unsigned o3 = ax[3].overflowIndex(x3, &i3);
2690  const unsigned o4 = ax[4].overflowIndex(x4, &i4);
2691  const unsigned o5 = ax[5].overflowIndex(x5, &i5);
2692  const unsigned o6 = ax[6].overflowIndex(x6, &i6);
2693  const unsigned o7 = ax[7].overflowIndex(x7, &i7);
2694  const unsigned o8 = ax[8].overflowIndex(x8, &i8);
2695  if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U && o4 == 1U &&
2696  o5 == 1U && o6 == 1U && o7 == 1U && o8 == 1U)
2697  return data_(i0, i1, i2, i3, i4, i5, i6, i7, i8);
2698  else
2699  return overflow_(o0, o1, o2, o3, o4, o5, o6, o7, o8);
2700  }
2701 
2702  template <typename Numeric, class Axis>
2703  const Numeric& HistoND<Numeric,Axis>::closestBin(const double x0,
2704  const double x1,
2705  const double x2,
2706  const double x3,
2707  const double x4,
2708  const double x5,
2709  const double x6,
2710  const double x7,
2711  const double x8) const
2712  {
2713  if (dim_ != 9U) Private::h_badargs("closestBin");
2714  const Axis* ax = &axes_[0];
2715  const unsigned i0 = ax[0].closestValidBin(x0);
2716  const unsigned i1 = ax[1].closestValidBin(x1);
2717  const unsigned i2 = ax[2].closestValidBin(x2);
2718  const unsigned i3 = ax[3].closestValidBin(x3);
2719  const unsigned i4 = ax[4].closestValidBin(x4);
2720  const unsigned i5 = ax[5].closestValidBin(x5);
2721  const unsigned i6 = ax[6].closestValidBin(x6);
2722  const unsigned i7 = ax[7].closestValidBin(x7);
2723  const unsigned i8 = ax[8].closestValidBin(x8);
2724  return data_(i0, i1, i2, i3, i4, i5, i6, i7, i8);
2725  }
2726 
2727  template <typename Numeric, class Axis>
2728  template <typename Num2>
2729  void HistoND<Numeric,Axis>::fill(const double x0, const double x1,
2730  const double x2, const double x3,
2731  const double x4, const double x5,
2732  const double x6, const double x7,
2733  const double x8, const double x9,
2734  const Num2& w)
2735  {
2736  if (dim_ != 10U) Private::h_badargs("fill");
2737  unsigned i0, i1, i2, i3, i4, i5, i6, i7, i8, i9;
2738  const Axis* ax = &axes_[0];
2739  const unsigned o0 = ax[0].overflowIndex(x0, &i0);
2740  const unsigned o1 = ax[1].overflowIndex(x1, &i1);
2741  const unsigned o2 = ax[2].overflowIndex(x2, &i2);
2742  const unsigned o3 = ax[3].overflowIndex(x3, &i3);
2743  const unsigned o4 = ax[4].overflowIndex(x4, &i4);
2744  const unsigned o5 = ax[5].overflowIndex(x5, &i5);
2745  const unsigned o6 = ax[6].overflowIndex(x6, &i6);
2746  const unsigned o7 = ax[7].overflowIndex(x7, &i7);
2747  const unsigned o8 = ax[8].overflowIndex(x8, &i8);
2748  const unsigned o9 = ax[9].overflowIndex(x9, &i9);
2749  if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U && o4 == 1U &&
2750  o5 == 1U && o6 == 1U && o7 == 1U && o8 == 1U && o9 == 1U)
2751  data_(i0, i1, i2, i3, i4, i5, i6, i7, i8, i9) += w;
2752  else
2753  {
2754  overflow_(o0, o1, o2, o3, o4, o5, o6, o7, o8, o9) += w;
2755  ++overCount_;
2756  }
2757  ++fillCount_; ++modCount_;
2758  }
2759 
2760  template <typename Numeric, class Axis>
2761  template <typename Num2, class Functor>
2762  void HistoND<Numeric,Axis>::dispatch(const double x0, const double x1,
2763  const double x2, const double x3,
2764  const double x4, const double x5,
2765  const double x6, const double x7,
2766  const double x8, const double x9,
2767  Num2& w, Functor& f)
2768  {
2769  if (dim_ != 10U) Private::h_badargs("dispatch");
2770  unsigned i0, i1, i2, i3, i4, i5, i6, i7, i8, i9;
2771  const Axis* ax = &axes_[0];
2772  const unsigned o0 = ax[0].overflowIndex(x0, &i0);
2773  const unsigned o1 = ax[1].overflowIndex(x1, &i1);
2774  const unsigned o2 = ax[2].overflowIndex(x2, &i2);
2775  const unsigned o3 = ax[3].overflowIndex(x3, &i3);
2776  const unsigned o4 = ax[4].overflowIndex(x4, &i4);
2777  const unsigned o5 = ax[5].overflowIndex(x5, &i5);
2778  const unsigned o6 = ax[6].overflowIndex(x6, &i6);
2779  const unsigned o7 = ax[7].overflowIndex(x7, &i7);
2780  const unsigned o8 = ax[8].overflowIndex(x8, &i8);
2781  const unsigned o9 = ax[9].overflowIndex(x9, &i9);
2782  if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U && o4 == 1U &&
2783  o5 == 1U && o6 == 1U && o7 == 1U && o8 == 1U && o9 == 1U)
2784  f(data_(i0, i1, i2, i3, i4, i5, i6, i7, i8, i9), w);
2785  else
2786  f(overflow_(o0, o1, o2, o3, o4, o5, o6, o7, o8, o9), w);
2787  ++modCount_;
2788  }
2789 
2790  template <typename Numeric, class Axis>
2791  template <typename Num2>
2792  void HistoND<Numeric,Axis>::fillC(const double x0, const double x1,
2793  const double x2, const double x3,
2794  const double x4, const double x5,
2795  const double x6, const double x7,
2796  const double x8, const double x9,
2797  const Num2& w)
2798  {
2799  if (dim_ != 10U) Private::h_badargs("fillC");
2800  double* wg = &weightBuf_[0];
2801  unsigned* idx = &indexBuf_[0];
2802  const Axis* ax = &axes_[0];
2803  const unsigned o0 = ax[0].overflowIndexWeighted(x0, idx+0, wg+0);
2804  const unsigned o1 = ax[1].overflowIndexWeighted(x1, idx+1, wg+1);
2805  const unsigned o2 = ax[2].overflowIndexWeighted(x2, idx+2, wg+2);
2806  const unsigned o3 = ax[3].overflowIndexWeighted(x3, idx+3, wg+3);
2807  const unsigned o4 = ax[4].overflowIndexWeighted(x4, idx+4, wg+4);
2808  const unsigned o5 = ax[5].overflowIndexWeighted(x5, idx+5, wg+5);
2809  const unsigned o6 = ax[6].overflowIndexWeighted(x6, idx+6, wg+6);
2810  const unsigned o7 = ax[7].overflowIndexWeighted(x7, idx+7, wg+7);
2811  const unsigned o8 = ax[8].overflowIndexWeighted(x8, idx+8, wg+8);
2812  const unsigned o9 = ax[9].overflowIndexWeighted(x9, idx+9, wg+9);
2813  if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U && o4 == 1U &&
2814  o5 == 1U && o6 == 1U && o7 == 1U && o8 == 1U && o9 == 1U)
2815  fillPreservingCentroid(w);
2816  else
2817  {
2818  overflow_(o0, o1, o2, o3, o4, o5, o6, o7, o8, o9) += w;
2819  ++overCount_;
2820  }
2821  ++fillCount_; ++modCount_;
2822  }
2823 
2824  template <typename Numeric, class Axis>
2825  const Numeric& HistoND<Numeric,Axis>::examine(
2826  const double x0, const double x1,
2827  const double x2, const double x3,
2828  const double x4, const double x5,
2829  const double x6, const double x7,
2830  const double x8,
2831  const double x9) const
2832  {
2833  if (dim_ != 10U) Private::h_badargs("examine");
2834  unsigned i0, i1, i2, i3, i4, i5, i6, i7, i8, i9;
2835  const Axis* ax = &axes_[0];
2836  const unsigned o0 = ax[0].overflowIndex(x0, &i0);
2837  const unsigned o1 = ax[1].overflowIndex(x1, &i1);
2838  const unsigned o2 = ax[2].overflowIndex(x2, &i2);
2839  const unsigned o3 = ax[3].overflowIndex(x3, &i3);
2840  const unsigned o4 = ax[4].overflowIndex(x4, &i4);
2841  const unsigned o5 = ax[5].overflowIndex(x5, &i5);
2842  const unsigned o6 = ax[6].overflowIndex(x6, &i6);
2843  const unsigned o7 = ax[7].overflowIndex(x7, &i7);
2844  const unsigned o8 = ax[8].overflowIndex(x8, &i8);
2845  const unsigned o9 = ax[9].overflowIndex(x9, &i9);
2846  if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U && o4 == 1U &&
2847  o5 == 1U && o6 == 1U && o7 == 1U && o8 == 1U && o9 == 1U)
2848  return data_(i0, i1, i2, i3, i4, i5, i6, i7, i8, i9);
2849  else
2850  return overflow_(o0, o1, o2, o3, o4, o5, o6, o7, o8, o9);
2851  }
2852 
2853  template <typename Numeric, class Axis>
2854  const Numeric& HistoND<Numeric,Axis>::closestBin(const double x0,
2855  const double x1,
2856  const double x2,
2857  const double x3,
2858  const double x4,
2859  const double x5,
2860  const double x6,
2861  const double x7,
2862  const double x8,
2863  const double x9) const
2864  {
2865  if (dim_ != 10U) Private::h_badargs("closestBin");
2866  const Axis* ax = &axes_[0];
2867  const unsigned i0 = ax[0].closestValidBin(x0);
2868  const unsigned i1 = ax[1].closestValidBin(x1);
2869  const unsigned i2 = ax[2].closestValidBin(x2);
2870  const unsigned i3 = ax[3].closestValidBin(x3);
2871  const unsigned i4 = ax[4].closestValidBin(x4);
2872  const unsigned i5 = ax[5].closestValidBin(x5);
2873  const unsigned i6 = ax[6].closestValidBin(x6);
2874  const unsigned i7 = ax[7].closestValidBin(x7);
2875  const unsigned i8 = ax[8].closestValidBin(x8);
2876  const unsigned i9 = ax[9].closestValidBin(x9);
2877  return data_(i0, i1, i2, i3, i4, i5, i6, i7, i8, i9);
2878  }
2879 
2880  template <typename Numeric, class Axis>
2881  template <typename Num2>
2882  inline void HistoND<Numeric,Axis>::setBin(const unsigned *index,
2883  const unsigned indexLen,
2884  const Num2& v)
2885  {
2886  data_.value(index, indexLen) = v;
2887  ++modCount_;
2888  }
2889 
2890  template <typename Numeric, class Axis>
2891  template <typename Num2>
2892  inline void HistoND<Numeric,Axis>::setBinAt(const unsigned *index,
2893  const unsigned indexLen,
2894  const Num2& v)
2895  {
2896  data_.valueAt(index, indexLen) = v;
2897  ++modCount_;
2898  }
2899 
2900  template <typename Numeric, class Axis>
2901  template <typename Num2>
2902  inline void HistoND<Numeric,Axis>::setBin(const Num2& v)
2903  {
2904  data_() = v;
2905  ++modCount_;
2906  }
2907 
2908  template <typename Numeric, class Axis>
2909  template <typename Num2>
2910  inline void HistoND<Numeric,Axis>::setBinAt(const Num2& v)
2911  {
2912  data_.at() = v;
2913  ++modCount_;
2914  }
2915 
2916  template <typename Numeric, class Axis>
2917  template <typename Num2>
2918  inline void HistoND<Numeric,Axis>::setBin(
2919  const unsigned i0, const Num2& v)
2920  {
2921  data_(i0) = v;
2922  ++modCount_;
2923  }
2924 
2925  template <typename Numeric, class Axis>
2926  template <typename Num2>
2927  inline void HistoND<Numeric,Axis>::setBinAt(
2928  const unsigned i0, const Num2& v)
2929  {
2930  data_.at(i0) = v;
2931  ++modCount_;
2932  }
2933 
2934  template <typename Numeric, class Axis>
2935  template <typename Num2>
2936  inline void HistoND<Numeric,Axis>::setBin(const unsigned i0,
2937  const unsigned i1,
2938  const Num2& v)
2939  {
2940  data_(i0, i1) = v;
2941  ++modCount_;
2942  }
2943 
2944  template <typename Numeric, class Axis>
2945  template <typename Num2>
2946  inline void HistoND<Numeric,Axis>::setBinAt(const unsigned i0,
2947  const unsigned i1,
2948  const Num2& v)
2949  {
2950  data_.at(i0, i1) = v;
2951  ++modCount_;
2952  }
2953 
2954  template <typename Numeric, class Axis>
2955  template <typename Num2>
2956  inline void HistoND<Numeric,Axis>::setBin(const unsigned i0,
2957  const unsigned i1,
2958  const unsigned i2,
2959  const Num2& v)
2960  {
2961  data_(i0, i1, i2) = v;
2962  ++modCount_;
2963  }
2964 
2965  template <typename Numeric, class Axis>
2966  template <typename Num2>
2967  inline void HistoND<Numeric,Axis>::setBin(const unsigned i0,
2968  const unsigned i1,
2969  const unsigned i2,
2970  const unsigned i3,
2971  const Num2& v)
2972  {
2973  data_(i0, i1, i2, i3) = v;
2974  ++modCount_;
2975  }
2976 
2977  template <typename Numeric, class Axis>
2978  template <typename Num2>
2979  inline void HistoND<Numeric,Axis>::setBin(const unsigned i0,
2980  const unsigned i1,
2981  const unsigned i2,
2982  const unsigned i3,
2983  const unsigned i4,
2984  const Num2& v)
2985  {
2986  data_(i0, i1, i2, i3, i4) = v;
2987  ++modCount_;
2988  }
2989 
2990  template <typename Numeric, class Axis>
2991  template <typename Num2>
2992  inline void HistoND<Numeric,Axis>::setBin(const unsigned i0,
2993  const unsigned i1,
2994  const unsigned i2,
2995  const unsigned i3,
2996  const unsigned i4,
2997  const unsigned i5,
2998  const Num2& v)
2999  {
3000  data_(i0, i1, i2, i3, i4, i5) = v;
3001  ++modCount_;
3002  }
3003 
3004  template <typename Numeric, class Axis>
3005  template <typename Num2>
3006  inline void HistoND<Numeric,Axis>::setBin(const unsigned i0,
3007  const unsigned i1,
3008  const unsigned i2,
3009  const unsigned i3,
3010  const unsigned i4,
3011  const unsigned i5,
3012  const unsigned i6,
3013  const Num2& v)
3014  {
3015  data_(i0, i1, i2, i3, i4, i5, i6) = v;
3016  ++modCount_;
3017  }
3018 
3019  template <typename Numeric, class Axis>
3020  template <typename Num2>
3021  inline void HistoND<Numeric,Axis>::setBin(const unsigned i0,
3022  const unsigned i1,
3023  const unsigned i2,
3024  const unsigned i3,
3025  const unsigned i4,
3026  const unsigned i5,
3027  const unsigned i6,
3028  const unsigned i7,
3029  const Num2& v)
3030  {
3031  data_(i0, i1, i2, i3, i4, i5, i6, i7) = v;
3032  ++modCount_;
3033  }
3034 
3035  template <typename Numeric, class Axis>
3036  template <typename Num2>
3037  inline void HistoND<Numeric,Axis>::setBin(const unsigned i0,
3038  const unsigned i1,
3039  const unsigned i2,
3040  const unsigned i3,
3041  const unsigned i4,
3042  const unsigned i5,
3043  const unsigned i6,
3044  const unsigned i7,
3045  const unsigned i8,
3046  const Num2& v)
3047  {
3048  data_(i0, i1, i2, i3, i4, i5, i6, i7, i8) = v;
3049  ++modCount_;
3050  }
3051 
3052  template <typename Numeric, class Axis>
3053  template <typename Num2>
3054  inline void HistoND<Numeric,Axis>::setBin(const unsigned i0,
3055  const unsigned i1,
3056  const unsigned i2,
3057  const unsigned i3,
3058  const unsigned i4,
3059  const unsigned i5,
3060  const unsigned i6,
3061  const unsigned i7,
3062  const unsigned i8,
3063  const unsigned i9,
3064  const Num2& v)
3065  {
3066  data_(i0, i1, i2, i3, i4, i5, i6, i7, i8, i9) = v;
3067  ++modCount_;
3068  }
3069 
3070  template <typename Numeric, class Axis>
3071  template <typename Num2>
3072  inline void HistoND<Numeric,Axis>::setBinAt(const unsigned i0,
3073  const unsigned i1,
3074  const unsigned i2,
3075  const Num2& v)
3076  {
3077  data_.at(i0, i1, i2) = v;
3078  ++modCount_;
3079  }
3080 
3081  template <typename Numeric, class Axis>
3082  template <typename Num2>
3083  inline void HistoND<Numeric,Axis>::setBinAt(const unsigned i0,
3084  const unsigned i1,
3085  const unsigned i2,
3086  const unsigned i3,
3087  const Num2& v)
3088  {
3089  data_.at(i0, i1, i2, i3) = v;
3090  ++modCount_;
3091  }
3092 
3093  template <typename Numeric, class Axis>
3094  template <typename Num2>
3095  inline void HistoND<Numeric,Axis>::setBinAt(const unsigned i0,
3096  const unsigned i1,
3097  const unsigned i2,
3098  const unsigned i3,
3099  const unsigned i4,
3100  const Num2& v)
3101  {
3102  data_.at(i0, i1, i2, i3, i4) = v;
3103  ++modCount_;
3104  }
3105 
3106  template <typename Numeric, class Axis>
3107  template <typename Num2>
3108  inline void HistoND<Numeric,Axis>::setBinAt(const unsigned i0,
3109  const unsigned i1,
3110  const unsigned i2,
3111  const unsigned i3,
3112  const unsigned i4,
3113  const unsigned i5,
3114  const Num2& v)
3115  {
3116  data_.at(i0, i1, i2, i3, i4, i5) = v;
3117  ++modCount_;
3118  }
3119 
3120  template <typename Numeric, class Axis>
3121  template <typename Num2>
3122  inline void HistoND<Numeric,Axis>::setBinAt(const unsigned i0,
3123  const unsigned i1,
3124  const unsigned i2,
3125  const unsigned i3,
3126  const unsigned i4,
3127  const unsigned i5,
3128  const unsigned i6,
3129  const Num2& v)
3130  {
3131  data_.at(i0, i1, i2, i3, i4, i5, i6) = v;
3132  ++modCount_;
3133  }
3134 
3135  template <typename Numeric, class Axis>
3136  template <typename Num2>
3137  inline void HistoND<Numeric,Axis>::setBinAt(const unsigned i0,
3138  const unsigned i1,
3139  const unsigned i2,
3140  const unsigned i3,
3141  const unsigned i4,
3142  const unsigned i5,
3143  const unsigned i6,
3144  const unsigned i7,
3145  const Num2& v)
3146  {
3147  data_.at(i0, i1, i2, i3, i4, i5, i6, i7) = v;
3148  ++modCount_;
3149  }
3150 
3151  template <typename Numeric, class Axis>
3152  template <typename Num2>
3153  inline void HistoND<Numeric,Axis>::setBinAt(const unsigned i0,
3154  const unsigned i1,
3155  const unsigned i2,
3156  const unsigned i3,
3157  const unsigned i4,
3158  const unsigned i5,
3159  const unsigned i6,
3160  const unsigned i7,
3161  const unsigned i8,
3162  const Num2& v)
3163  {
3164  data_.at(i0, i1, i2, i3, i4, i5, i6, i7, i8) = v;
3165  ++modCount_;
3166  }
3167 
3168  template <typename Numeric, class Axis>
3169  template <typename Num2>
3170  inline void HistoND<Numeric,Axis>::setBinAt(const unsigned i0,
3171  const unsigned i1,
3172  const unsigned i2,
3173  const unsigned i3,
3174  const unsigned i4,
3175  const unsigned i5,
3176  const unsigned i6,
3177  const unsigned i7,
3178  const unsigned i8,
3179  const unsigned i9,
3180  const Num2& v)
3181  {
3182  data_.at(i0, i1, i2, i3, i4, i5, i6, i7, i8, i9) = v;
3183  ++modCount_;
3184  }
3185 
3186  template <typename Numeric, class Axis>
3187  template <typename Num2>
3188  inline HistoND<Numeric,Axis>& HistoND<Numeric,Axis>::operator+=(
3189  const HistoND<Num2,Axis>& r)
3190  {
3191  data_ += r.data_;
3192  overflow_ += r.overflow_;
3193  fillCount_ += r.fillCount_;
3194  overCount_ += r.overCount_;
3195  ++modCount_;
3196  return *this;
3197  }
3198 
3199  template <typename Numeric, class Axis>
3200  template <typename Num2>
3201  inline HistoND<Numeric,Axis>& HistoND<Numeric,Axis>::operator-=(
3202  const HistoND<Num2,Axis>& r)
3203  {
3204  data_ -= r.data_;
3205  overflow_ -= r.overflow_;
3206 
3207  // Subtraction does not make much sense for fill counts.
3208  // We will assume that what we want should be equivalent
3209  // to the in-place multiplication of the other histogram
3210  // by -1 and then adding.
3211  //
3212  fillCount_ += r.fillCount_;
3213  overCount_ += r.overCount_;
3214 
3215  ++modCount_;
3216  return *this;
3217  }
3218 
3219  template <typename Numeric, class Axis>
3220  template <typename Num2>
3221  inline HistoND<Numeric,Axis>&
3222  HistoND<Numeric,Axis>::operator*=(const Num2& r)
3223  {
3224  data_ *= r;
3225  overflow_ *= r;
3226  ++modCount_;
3227  return *this;
3228  }
3229 
3230  template <typename Numeric, class Axis>
3231  template <typename Num2>
3232  inline HistoND<Numeric,Axis>&
3233  HistoND<Numeric,Axis>::operator/=(const Num2& r)
3234  {
3235  data_ /= r;
3236  overflow_ /= r;
3237  ++modCount_;
3238  return *this;
3239  }
3240 
3241  template <typename Numeric, class Axis>
3242  HistoND<Numeric,Axis>::HistoND(const HistoND& r,
3243  const unsigned ax1,
3244  const unsigned ax2)
3245  : title_(r.title_),
3246  accumulatedDataLabel_(r.accumulatedDataLabel_),
3247  data_(r.data_.transpose(ax1, ax2)),
3248  overflow_(r.overflow_.transpose(ax1, ax2)),
3249  axes_(r.axes_),
3250  weightBuf_(r.weightBuf_),
3251  indexBuf_(r.indexBuf_),
3252  fillCount_(r.fillCount_),
3253  overCount_(r.overCount_),
3254  modCount_(0UL),
3255  dim_(r.dim_)
3256  {
3257  std::swap(axes_[ax1], axes_[ax2]);
3258  }
3259 
3260  template <typename Numeric, class Axis>
3261  HistoND<Numeric,Axis>::HistoND(const HistoND& r)
3262  : title_(r.title_),
3263  accumulatedDataLabel_(r.accumulatedDataLabel_),
3264  data_(r.data_),
3265  overflow_(r.overflow_),
3266  axes_(r.axes_),
3267  weightBuf_(r.weightBuf_),
3268  indexBuf_(r.indexBuf_),
3269  fillCount_(r.fillCount_),
3270  overCount_(r.overCount_),
3271  modCount_(0UL),
3272  dim_(r.dim_)
3273  {
3274  }
3275 
3276  template <typename Numeric, class Axis>
3277  HistoND<Numeric,Axis>& HistoND<Numeric,Axis>::operator=(const HistoND& r)
3278  {
3279  if (&r != this)
3280  {
3281  title_ = r.title_;
3282  accumulatedDataLabel_ = r.accumulatedDataLabel_;
3283  data_.uninitialize();
3284  data_ = r.data_;
3285  overflow_.uninitialize();
3286  overflow_ = r.overflow_;
3287  axes_ = r.axes_;
3288  weightBuf_ = r.weightBuf_;
3289  indexBuf_ = r.indexBuf_;
3290  fillCount_ = r.fillCount_;
3291  overCount_ = r.overCount_;
3292  dim_ = r.dim_;
3293  ++modCount_;
3294  }
3295  return *this;
3296  }
3297 
3298  template <typename Numeric, class Axis>
3299  HistoND<Numeric,Axis> HistoND<Numeric,Axis>::transpose(
3300  const unsigned axisNum1, const unsigned axisNum2) const
3301  {
3302  if (axisNum1 >= dim_ || axisNum2 >= dim_)
3303  throw npstat::NpstatOutOfRange("In npstat::HistoND::transpose: "
3304  "axis number is out of range");
3305  if (axisNum1 == axisNum2)
3306  // Just make a copy
3307  return *this;
3308  else
3309  return HistoND(*this, axisNum1, axisNum2);
3310  }
3311 
3312  template <typename Numeric, class Axis>
3313  template <typename Num2>
3314  void HistoND<Numeric,Axis>::addToBinContents(const Num2& weight)
3315  {
3316  const unsigned long nDat = data_.length();
3317  Numeric* data = const_cast<Numeric*>(data_.data());
3318  for (unsigned long i=0; i<nDat; ++i)
3319  data[i] += weight;
3320  fillCount_ += nDat;
3321  ++modCount_;
3322  }
3323 
3324  template <typename Numeric, class Axis>
3325  template <typename Num2>
3326  void HistoND<Numeric,Axis>::addToOverflows(const Num2& weight)
3327  {
3328  const unsigned long nOver = overflow_.length();
3329  Numeric* data = const_cast<Numeric*>(overflow_.data());
3330  for (unsigned long i=0; i<nOver; ++i)
3331  data[i] += weight;
3332  overCount_ += nOver;
3333  fillCount_ += nOver;
3334  ++modCount_;
3335  }
3336 
3337  template <typename Numeric, class Axis>
3338  template <typename Num2>
3339  void HistoND<Numeric,Axis>::addToBinContents(
3340  const Num2* data, const unsigned long dataLength)
3341  {
3342  if (dataLength != data_.length()) throw npstat::NpstatInvalidArgument(
3343  "In npstat::HistoND::addToBinContents: incompatible data length");
3344  assert(data);
3345  Numeric* dat = const_cast<Numeric*>(data_.data());
3346  for (unsigned long i=0; i<dataLength; ++i)
3347  dat[i] += data[i];
3348  fillCount_ += dataLength;
3349  ++modCount_;
3350  }
3351 
3352  template <typename Numeric, class Axis>
3353  template <typename Num2>
3354  void HistoND<Numeric,Axis>::addToOverflows(
3355  const Num2* data, const unsigned long dataLength)
3356  {
3357  if (dataLength != overflow_.length()) throw npstat::NpstatInvalidArgument(
3358  "In npstat::HistoND::addToOverflows: incompatible data length");
3359  assert(data);
3360  Numeric* dat = const_cast<Numeric*>(overflow_.data());
3361  for (unsigned long i=0; i<dataLength; ++i)
3362  dat[i] += data[i];
3363  overCount_ += dataLength;
3364  fillCount_ += dataLength;
3365  ++modCount_;
3366  }
3367 
3368  template <typename Numeric, class Axis>
3369  template <typename Num2>
3370  void HistoND<Numeric,Axis>::scaleBinContents(
3371  const Num2* data, const unsigned long dataLength)
3372  {
3373  if (dataLength != data_.length()) throw npstat::NpstatInvalidArgument(
3374  "In npstat::HistoND::scaleBinContents: incompatible data length");
3375  assert(data);
3376  Numeric* dat = const_cast<Numeric*>(data_.data());
3377  for (unsigned long i=0; i<dataLength; ++i)
3378  dat[i] *= data[i];
3379  ++modCount_;
3380  }
3381 
3382  template <typename Numeric, class Axis>
3383  template <typename Num2>
3384  void HistoND<Numeric,Axis>::scaleOverflows(
3385  const Num2* data, const unsigned long dataLength)
3386  {
3387  if (dataLength != overflow_.length()) throw npstat::NpstatInvalidArgument(
3388  "In npstat::HistoND::scaleOverflows: incompatible data length");
3389  assert(data);
3390  Numeric* dat = const_cast<Numeric*>(overflow_.data());
3391  for (unsigned long i=0; i<dataLength; ++i)
3392  dat[i] *= data[i];
3393  ++modCount_;
3394  }
3395 
3396  template <typename Numeric, class Axis>
3397  template <typename Num2>
3398  inline void HistoND<Numeric,Axis>::setBinContents(
3399  const Num2* data, const unsigned long dataLength,
3400  const bool clearOverflowsNow)
3401  {
3402  data_.setData(data, dataLength);
3403  if (clearOverflowsNow)
3404  clearOverflows();
3405  ++modCount_;
3406  }
3407 
3408  template <typename Numeric, class Axis>
3409  template <typename Num2>
3410  inline void HistoND<Numeric,Axis>::setOverflows(
3411  const Num2* data, const unsigned long dataLength)
3412  {
3413  overflow_.setData(data, dataLength);
3414  ++modCount_;
3415  }
3416 
3417  template <typename Numeric, class Axis>
3418  inline void HistoND<Numeric,Axis>::recalculateNFillsFromData()
3419  {
3420  const long double nOver = overflow_.template sum<long double>();
3421  const long double nData = data_.template sum<long double>();
3422  overCount_ = static_cast<unsigned long>(nOver);
3423  fillCount_ = static_cast<unsigned long>(nData + nOver);
3424  ++modCount_;
3425  }
3426 
3427  template <typename Numeric, class Axis>
3428  template <typename Num2, typename Num3>
3429  inline void HistoND<Numeric,Axis>::addToProjection(
3430  HistoND<Num2,Axis>* projection,
3431  AbsArrayProjector<Numeric,Num3>& projector,
3432  const unsigned *projectedIndices,
3433  const unsigned nProjectedIndices) const
3434  {
3435  assert(projection);
3436  data_.addToProjection(&projection->data_, projector,
3437  projectedIndices, nProjectedIndices);
3438  projection->fillCount_ += projection->nBins();
3439  projection->modCount_++;
3440  }
3441 
3442  template <typename Numeric, class Axis>
3443  template <typename Num2, typename Num3>
3444  inline void HistoND<Numeric,Axis>::addToProjection(
3445  HistoND<Num2,Axis>* projection,
3446  AbsVisitor<Numeric,Num3>& projector,
3447  const unsigned *projectedIndices,
3448  const unsigned nProjectedIndices) const
3449  {
3450  assert(projection);
3451  data_.addToProjection(&projection->data_, projector,
3452  projectedIndices, nProjectedIndices);
3453  projection->fillCount_ += projection->nBins();
3454  projection->modCount_++;
3455  }
3456 
3457  template <typename Numeric, class Axis>
3458  const char* HistoND<Numeric,Axis>::classname()
3459  {
3460  static const std::string myClass(gs::template_class_name<Numeric,Axis>(
3461  "npstat::HistoND"));
3462  return myClass.c_str();
3463  }
3464 
3465  template<typename Numeric, class Axis>
3466  bool HistoND<Numeric,Axis>::write(std::ostream& of) const
3467  {
3468  gs::write_pod(of, title_);
3469  gs::write_pod(of, accumulatedDataLabel_);
3470  gs::write_pod(of, fillCount_);
3471  gs::write_pod(of, overCount_);
3472 
3473  return !of.fail() &&
3474  gs::write_obj_vector(of, axes_) &&
3475  data_.classId().write(of) &&
3476  data_.write(of) &&
3477  overflow_.write(of);
3478  }
3479 
3480  template<typename Numeric, class Axis>
3481  HistoND<Numeric,Axis>* HistoND<Numeric,Axis>::read(const gs::ClassId& id,
3482  std::istream& in)
3483  {
3484  static const gs::ClassId current(
3485  gs::ClassId::makeId<HistoND<Numeric,Axis> >());
3486  current.ensureSameId(id);
3487 
3488  std::string title;
3489  gs::read_pod(in, &title);
3490 
3491  std::string accumulatedDataLabel;
3492  gs::read_pod(in, &accumulatedDataLabel);
3493 
3494  unsigned long fillCount = 0, overCount = 0;
3495  gs::read_pod(in, &fillCount);
3496  gs::read_pod(in, &overCount);
3497  if (in.fail()) throw gs::IOReadFailure(
3498  "In npstat::HistoND::read: input stream failure");
3499 
3500  std::vector<Axis> axes;
3501  gs::read_heap_obj_vector_as_placed(in, &axes);
3502  gs::ClassId ida(in, 1);
3503  ArrayND<Numeric> data, over;
3504  ArrayND<Numeric>::restore(ida, in, &data);
3505  ArrayND<Numeric>::restore(ida, in, &over);
3506  CPP11_auto_ptr<HistoND<Numeric,Axis> > result(
3507  new HistoND<Numeric,Axis>(axes, title.c_str(),
3508  accumulatedDataLabel.c_str()));
3509  result->data_ = data;
3510  result->overflow_ = over;
3511  result->fillCount_ = fillCount;
3512  result->overCount_ = overCount;
3513  return result.release();
3514  }
3515 
3516  template <typename Histo>
3517  void convertHistoToDensity(Histo* h, const bool knownNonNegative)
3518  {
3519  assert(h);
3520  if (!knownNonNegative)
3521  (const_cast<ArrayND<typename
3522  Histo::value_type>&>(h->binContents())).makeNonNegative();
3523  const double integ = h->integral();
3524  *h /= integ;
3525  }
3526 
3527  template <typename Histo>
3528  std::vector<LinearMapper1d> densityScanHistoMap(const Histo& histo)
3529  {
3530  std::vector<LinearMapper1d> result;
3531  const unsigned d = histo.dim();
3532  result.reserve(d);
3533  for (unsigned i=0; i<d; ++i)
3534  {
3535  const LinearMapper1d& m = histo.axis(i).binNumberMapper(false);
3536  result.push_back(m.inverse());
3537  }
3538  return result;
3539  }
3540 
3541  template <typename Histo>
3542  std::vector<CircularMapper1d> convolutionHistoMap(
3543  const Histo& histo, const bool doubleRange)
3544  {
3545  std::vector<CircularMapper1d> result;
3546  const unsigned d = histo.dim();
3547  result.reserve(d);
3548  for (unsigned i=0; i<d; ++i)
3549  result.push_back(histo.axis(i).kernelScanMapper(doubleRange));
3550  return result;
3551  }
3552 }
3553 
3554 
3555 #endif // NPSTAT_HISTOND_HH_
3556 
i
int i
Definition: DBlmapReader.cc:9
npstat::HistoND::setBinAt
void setBinAt(const unsigned *index, unsigned indexLen, const Num2 &v)
Definition: HistoND.h:2892
npstat::HistoND::version
static unsigned version()
Definition: HistoND.h:813
npstat::HistoND::operator==
bool operator==(const HistoND &) const
Definition: HistoND.h:1485
npstat::HistoND::fillCount_
unsigned long fillCount_
Definition: HistoND.h:837
npstat::HistoND::dim_
unsigned dim_
Definition: HistoND.h:840
npstat::HistoND::getModCount
unsigned long getModCount() const
Definition: HistoND.h:797
npstat::HistoND::write
bool write(std::ostream &of) const
Definition: HistoND.h:3466
npstat::HistoND::operator+=
HistoND & operator+=(const HistoND< Num2, Axis > &r)
npstat::convolutionHistoMap
std::vector< CircularMapper1d > convolutionHistoMap(const Histo &histo, bool doubleDataRange)
Definition: HistoND.h:3542
npstat::HistoND::addToProjection
void addToProjection(HistoND< Num2, Axis > *projection, AbsArrayProjector< Numeric, Num3 > &projector, const unsigned *projectedIndices, unsigned nProjectedIndices) const
Definition: HistoND.h:3429
npstat::HistoND::accumulateBinsInBox
void accumulateBinsInBox(const BoxND< double > &box, Acc *acc) const
Definition: HistoND.h:1123
npstat::HistoND::accumulatedDataLabel
const std::string & accumulatedDataLabel() const
Definition: HistoND.h:135
npstat::HistoND::clear
void clear()
Definition: HistoND.h:1157
npstat::HistoND::modCount_
unsigned long modCount_
Definition: HistoND.h:839
npstat::HistoND::examine
const Numeric & examine() const
Definition: HistoND.h:1733
npstat::HistoND::isSameData
bool isSameData(const HistoND &) const
Definition: HistoND.h:1477
npstat::HistoND::nFillsTotal
unsigned long nFillsTotal() const
Definition: HistoND.h:155
npstat::HistoND::binCenter
void binCenter(unsigned long binNumber, double *coords, unsigned lenCoords) const
Definition: HistoND.h:1417
npstat::HistoND::axis
const Axis & axis(const unsigned i) const
Definition: HistoND.h:148
PFRecoTauDiscriminationAgainstElectronMVA3GBR_cfi.method
tuple method
Definition: PFRecoTauDiscriminationAgainstElectronMVA3GBR_cfi.py:15
npstat::HistoND::operator=
HistoND & operator=(const HistoND &)
Definition: HistoND.h:3277
npstat::HistoND::dispatch
void dispatch(const double *coords, unsigned coordLength, Num2 &weight, Functor &f)
Definition: HistoND.h:1566
npstat::ArrayShape
std::vector< unsigned > ArrayShape
Definition: ArrayShape.h:21
npstat::HistoND::operator*=
HistoND & operator*=(const Num2 &r)
npstat::ArrayND::linearValue
Numeric & linearValue(unsigned long index)
Definition: ArrayND.h:3139
npstat::ArrayND::length
unsigned long length() const
Definition: ArrayND.h:231
npstat::HistoND::value_type
Numeric value_type
Definition: HistoND.h:48
npstat::HistoND::setLinearBin
void setLinearBin(const unsigned long index, const Num2 &v)
Definition: HistoND.h:604
npstat::HistoND::fill
void fill(const double *coords, unsigned coordLength, const Num2 &weight)
Definition: HistoND.h:1533
npstat::HistoND::title
const std::string & title() const
Definition: HistoND.h:132
npstat::HistoND::title_
std::string title_
Definition: HistoND.h:830
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
create_public_lumi_plots.ax
tuple ax
Definition: create_public_lumi_plots.py:1091
npstat::ArrayND::linearValueAt
Numeric & linearValueAt(unsigned long index)
Definition: ArrayND.h:3153
npstat::HistoND::recalculateNFillsFromData
void recalculateNFillsFromData()
Definition: HistoND.h:3418
npstat::HistoND::scaleBinContents
void scaleBinContents(const Num2 *data, unsigned long dataLength)
Definition: HistoND.h:3370
npstat::LinearMapper1d::inverse
LinearMapper1d inverse() const
Definition: LinearMapper1d.h:54
npstat::HistoND::accumulatedDataLabel_
std::string accumulatedDataLabel_
Definition: HistoND.h:831
npstat::HistoND::volume
double volume() const
Definition: HistoND.h:1519
npstat::convertHistoToDensity
void convertHistoToDensity(Histo *histogram, bool knownNonNegative=false)
Definition: HistoND.h:3517
npstat::ArrayND::constFill
ArrayND & constFill(Numeric c)
Definition: ArrayND.h:4584
npstat::HistoND::setNFillsOver
void setNFillsOver(const unsigned long i)
Definition: HistoND.h:701
npstat::HistoND::binVolume
double binVolume(unsigned long binNumber=0) const
Definition: HistoND.h:1504
npstat::HistoND::classId
gs::ClassId classId() const
Definition: HistoND.h:808
npstat::HistoND::axis_type
Axis axis_type
Definition: HistoND.h:49
npstat::HistoND::binContents
const ArrayND< Numeric > & binContents() const
Definition: HistoND.h:139
NpstatException.h
Exceptions for the npstat namespace.
npstat::Private::h_badargs
void h_badargs(const char *method)
Definition: HistoND.h:1087
npstat::HistoND::setAccumulatedDataLabel
void setAccumulatedDataLabel(const char *newlabel)
Definition: HistoND.h:175
npstat::HistoND::setOverflowsToConst
void setOverflowsToConst(const Num2 &value)
Definition: HistoND.h:683
npstat::HistoND::addToBinContents
void addToBinContents(const Num2 &weight)
Definition: HistoND.h:3314
max
const T & max(const T &a, const T &b)
Definition: MaterialBudgetTrackerHistos.cc:4
npstat::HistoND::accumulateBinsLoop
void accumulateBinsLoop(unsigned level, const BoxND< double > &box, unsigned *idx, Acc *accumulator, double overlapFraction) const
Definition: HistoND.h:1098
std::swap
void swap(edm::DataFrameContainer &lhs, edm::DataFrameContainer &rhs)
Definition: DataFrameContainer.h:223
npstat::HistoND::setBin
void setBin(const unsigned *index, unsigned indexLen, const Num2 &v)
Definition: HistoND.h:2882
npstat::HistoND::classname
static const char * classname()
Definition: HistoND.h:3458
query.result
tuple result
Definition: query.py:137
npstat::HistoND::addToOverflows
void addToOverflows(const Num2 &weight)
Definition: HistoND.h:3326
npstat::HistoND::setNFillsTotal
void setNFillsTotal(const unsigned long i)
Definition: HistoND.h:699
npstat::HistoND::data_
ArrayND< Numeric > data_
Definition: HistoND.h:832
Point
math::XYZPoint Point
Definition: TrackerDpgAnalysis.cc:107
npstat::HistoND::fillC
void fillC(const double *coords, unsigned coordLength, const Num2 &weight)
Definition: HistoND.h:1673
j
int j
Definition: DBlmapReader.cc:9
npstat::PreciseType::type
PreciseTypeHelper< T, gs::IOIsNumber< T >::value >::type type
Definition: PreciseType.h:39
f
double f[11][100]
Definition: MuScleFitUtils.cc:78
npstat::HistoND::isUniformlyBinned
bool isUniformlyBinned() const
Definition: HistoND.h:1371
npstat::Private::makeHistoShape
ArrayShape makeHistoShape(const std::vector< Axis > &axes)
Definition: HistoND.h:896
npstat::HistoND::nFillsOver
unsigned long nFillsOver() const
Definition: HistoND.h:162
npstat::HistoND::integral
double integral() const
Definition: HistoND.h:1380
cond::ecalcond::value_type
Container::value_type value_type
Definition: EcalChannelStatusPyWrapper.cc:33
npstat::HistoND::weightBuf_
std::vector< double > weightBuf_
Definition: HistoND.h:835
npstat::HistoND::axes
const std::vector< Axis > & axes() const
Definition: HistoND.h:145
npstat::HistoND::overCount_
unsigned long overCount_
Definition: HistoND.h:838
npstat::HistoND::allBinCenters
void allBinCenters(std::vector< Point > *centers) const
Definition: HistoND.h:1436
npstat::ArrayND::restore
static void restore(const gs::ClassId &id, std::istream &in, ArrayND *array)
Definition: ArrayND.h:5870
npstat::HistoND::indexBuf_
std::vector< unsigned > indexBuf_
Definition: HistoND.h:836
h
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
npstat::HistoND::clearOverflows
void clearOverflows()
Definition: HistoND.h:1149
customizeTrackingMonitorSeedNumber.idx
tuple idx
DEBUGGING if hasattr(process,&quot;trackMonIterativeTracking2012&quot;): print &quot;trackMonIterativeTracking2012 D...
Definition: customizeTrackingMonitorSeedNumber.py:15
npstat::HistoND::nBins
unsigned long nBins() const
Definition: HistoND.h:152
npstat::Private::addAxis
std::vector< Axis > addAxis(const std::vector< Axis > &axes, const Axis &newAxis, const unsigned newAxisNumber)
Definition: HistoND.h:1042
npstat::HistoND::closestBin
const Numeric & closestBin() const
Definition: HistoND.h:1740
npstat::HistoND::boundingBox
BoxND< double > boundingBox() const
Definition: HistoND.h:1403
npstat::HistoND::operator-=
HistoND & operator-=(const HistoND< Num2, Axis > &r)
create_public_pileup_plots.weights
list weights
Definition: create_public_pileup_plots.py:205
npstat::HistoND::setTitle
void setTitle(const char *newtitle)
Definition: HistoND.h:171
npstat::HistoND::incrModCount
void incrModCount()
Definition: HistoND.h:804
npstat::HistoND::operator/=
HistoND & operator/=(const Num2 &r)
npstat::HistoND::setOverflows
void setOverflows(const Num2 *data, unsigned long dataLength)
Definition: HistoND.h:3410
npstat::HistoND::binBox
void binBox(unsigned long binNumber, BoxND< double > *box) const
Definition: HistoND.h:1461
npstat::HistoND::fillPreservingCentroid
void fillPreservingCentroid(const Num2 &weight)
Definition: HistoND.h:1644
data
char data[epos_bytes_allocation]
Definition: EPOS_Wrapper.h:82
npstat::HistoND::nFillsInRange
unsigned long nFillsInRange() const
Definition: HistoND.h:158
npstat::Private::axesOfASlice
std::vector< Axis > axesOfASlice(const std::vector< Axis > &axes, const unsigned *fixedIndices, const unsigned nFixedIndices)
Definition: HistoND.h:967
npstat::HistoND::setBinContents
void setBinContents(const Num2 *data, unsigned long dataLength, bool clearOverflows=true)
Definition: HistoND.h:3398
npstat::HistoND::setAxisLabel
void setAxisLabel(const unsigned axisNum, const char *newlabel)
Definition: HistoND.h:179
npstat::HistoND::setBinsToConst
void setBinsToConst(const Num2 &value)
Definition: HistoND.h:675
npstat::HistoND::read
static HistoND * read(const gs::ClassId &id, std::istream &in)
Definition: HistoND.h:3481
npstat::HistoND::operator!=
bool operator!=(const HistoND &) const
Definition: HistoND.h:1498
npstat::HistoND::setLinearBinAt
void setLinearBinAt(const unsigned long index, const Num2 &v)
Definition: HistoND.h:657
npstat::HistoND::axes_
std::vector< Axis > axes_
Definition: HistoND.h:834
testEve_cfg.level
tuple level
Definition: testEve_cfg.py:34
npstat::Private::shapeOfASlice
ArrayShape shapeOfASlice(const std::vector< Axis > &axes, const unsigned *fixedIndices, const unsigned nFixedIndices)
Definition: HistoND.h:1002
w
T w() const
Definition: extBasic3DVector.h:245
histoStyle.weight
int weight
Definition: histoStyle.py:50
npstat::HistoND::overflow_
ArrayND< Numeric > overflow_
Definition: HistoND.h:833
npstat::HistoND::clearBinContents
void clearBinContents()
Definition: HistoND.h:1141
npstat::HistoND::overflows
const ArrayND< Numeric > & overflows() const
Definition: HistoND.h:142
npstat::densityScanHistoMap
std::vector< LinearMapper1d > densityScanHistoMap(const Histo &histo)
Definition: HistoND.h:3528
npstat::HistoND::transpose
HistoND transpose(unsigned axisNum1, unsigned axisNum2) const
Definition: HistoND.h:3299
findQualityFiles.size
tuple size
Write out results.
Definition: findQualityFiles.py:442
HistoAxis.h
Histogram axis with equidistant bins.
npstat::HistoND::scaleOverflows
void scaleOverflows(const Num2 *data, unsigned long dataLength)
Definition: HistoND.h:3384
npstat::HistoND::dim
unsigned dim() const
Definition: HistoND.h:129
ArrayND.h
Arbitrary-dimensional array template.
npstat::Private::shapeWithExtraAxis
ArrayShape shapeWithExtraAxis(const std::vector< Axis > &axes, const Axis &newAxis, const unsigned newAxisNumber)
Definition: HistoND.h:1065
Generated for CMSSW Reference Manual by doxygen 1.8.5