Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-07-30 02:33:14

0001 #ifndef NPSTAT_HISTOND_HH_
0002 #define NPSTAT_HISTOND_HH_
0003 
0004 /*! 
0005 // \file HistoND.h
0006 //
0007 // \brief Arbitrary-dimensional histogram template
0008 //
0009 // Author: I. Volobouev
0010 //
0011 // July 2010
0012 */
0013 
0014 #include "JetMETCorrections/InterpolationTables/interface/ArrayND.h"
0015 #include "JetMETCorrections/InterpolationTables/interface/HistoAxis.h"
0016 
0017 namespace npstat {
0018   /**
0019     // (Almost) arbitrary-dimensional histogram with binning determined
0020     // by the second template parameter (typically HistoAxis or NUHistoAxis).
0021     // The dimensionality must not exceed CHAR_BIT*sizeof(unsigned long)-1
0022     // which is normally 31/63 on 32/64-bit systems.
0023     //
0024     // The template parameter class (Numeric) must be such that it can be
0025     // used as the template parameter of ArrayND class. For a typical usage
0026     // pattern, Numeric should also support operator += between itself and
0027     // the weights with which the histogram is filled (see, however, the
0028     // description of the "dispatch" method which is not subject to
0029     // this recommendation).
0030     //
0031     // If the "fillC" method is used to accumulate the data then the weights
0032     // must support multiplication by a double, and then it must be possible
0033     // to use the "+=" operator to add such a product to Numeric.
0034     //
0035     // Note that there are no methods which would allow the user to examine
0036     // the bin contents of the histogram using bin numbers. This is
0037     // intentional: almost always such examinations are performed in a loop
0038     // over indices, and it is more efficient to grab a reference to the 
0039     // underlying array using the "binContents()" method and then examine
0040     // that array directly.
0041     */
0042   template <typename Numeric, class Axis = HistoAxis>
0043   class HistoND {
0044     template <typename Num2, class Axis2>
0045     friend class HistoND;
0046 
0047   public:
0048     typedef Numeric value_type;
0049     typedef Axis axis_type;
0050 
0051     enum RebinType { SAMPLE = 0, SUM, AVERAGE };
0052 
0053     /** Main constructor for arbitrary-dimensional histograms */
0054     explicit HistoND(const std::vector<Axis>& axes,
0055                      const char* title = nullptr,
0056                      const char* accumulatedDataLabel = nullptr);
0057 
0058     /** Convenience constructor for 1-d histograms */
0059     explicit HistoND(const Axis& xAxis, const char* title = nullptr, const char* accumulatedDataLabel = nullptr);
0060 
0061     /** Convenience constructor for 2-d histograms */
0062     HistoND(const Axis& xAxis,
0063             const Axis& yAxis,
0064             const char* title = nullptr,
0065             const char* accumulatedDataLabel = nullptr);
0066 
0067     /** Convenience constructor for 3-d histograms */
0068     HistoND(const Axis& xAxis,
0069             const Axis& yAxis,
0070             const Axis& zAxis,
0071             const char* title = nullptr,
0072             const char* accumulatedDataLabel = nullptr);
0073 
0074     /** Convenience constructor for 4-d histograms */
0075     HistoND(const Axis& xAxis,
0076             const Axis& yAxis,
0077             const Axis& zAxis,
0078             const Axis& tAxis,
0079             const char* title = nullptr,
0080             const char* accumulatedDataLabel = nullptr);
0081 
0082     /** Convenience constructor for 5-d histograms */
0083     HistoND(const Axis& xAxis,
0084             const Axis& yAxis,
0085             const Axis& zAxis,
0086             const Axis& tAxis,
0087             const Axis& vAxis,
0088             const char* title = nullptr,
0089             const char* accumulatedDataLabel = nullptr);
0090 
0091     /**
0092         // Simple constructor for uniformly binned histograms without
0093         // axis labels. Sequence size returned by the size() method of
0094         // both "shape" and "boundingBox" arguments must be the same.
0095         */
0096     HistoND(const ArrayShape& shape,
0097             const BoxND<double>& boundingBox,
0098             const char* title = nullptr,
0099             const char* accumulatedDataLabel = nullptr);
0100 
0101     /**
0102         // Converting constructor. The functor will be applied to all bins
0103         // of the argument histogram to fill the bins of the constructed
0104         // histogram. If the title and data label are not provided, they
0105         // will be cleared.
0106         */
0107     template <typename Num2, class Functor>
0108     HistoND(const HistoND<Num2, Axis>& h,
0109             const Functor& f,
0110             const char* title = nullptr,
0111             const char* accumulatedDataLabel = nullptr);
0112 
0113     /**
0114         // A slicing constructor. The new histogram will be created by
0115         // slicing another histogram. See the description of the slicing
0116         // constructor in the "ArrayND" class for the meaning of arguments
0117         // "indices" and "nIndices". The data of the newly created histogram
0118         // is cleared.
0119         */
0120     template <typename Num2>
0121     HistoND(const HistoND<Num2, Axis>& h, const unsigned* indices, unsigned nIndices, const char* title = nullptr);
0122 
0123     /**
0124         // A constructor that inserts a new axis into a histogram
0125         // (as if the argument histogram was a slice of the new histogram).
0126         // The "newAxisNumber" argument specifies the number of the
0127         // new axis in the axis sequence of the constructed histogram.
0128         // If the "newAxisNumber" exceeds the number of axes of the
0129         // argument histogram, the new axis will become last. The data
0130         // of the newly created histogram is cleared.
0131         */
0132     template <typename Num2>
0133     HistoND(const HistoND<Num2, Axis>& h, const Axis& newAxis, unsigned newAxisNumber, const char* title = nullptr);
0134 
0135     /**
0136         // Create a rebinned histogram with the same axis coverage.
0137         // Note that not all such operations will be meaningful if the
0138         // bin contents do not belong to one of the floating point types.
0139         // The "newBinCounts" argument specifies the new number of bins
0140         // along each axis. The length of this array (provided by the
0141         // "lenNewBinCounts" argument) should be equal to the input
0142         // histogram dimensionality.
0143         //
0144         // The "shifts" argument can be meaningfully specified with the
0145         // "rType" argument set to "SAMPLE". These shifts will be added
0146         // to the bin centers of the created histogram when the bin contents
0147         // are looked up in the input histogram. This can be useful in case
0148         // the bin center lookup without shifts would fall exactly on the
0149         // bin edge. Naturally, the length of the "shifts" array should be
0150         // equal to the input histogram dimensionality.
0151         */
0152     template <typename Num2>
0153     HistoND(const HistoND<Num2, Axis>& h,
0154             RebinType rType,
0155             const unsigned* newBinCounts,
0156             unsigned lenNewBinCounts,
0157             const double* shifts = nullptr,
0158             const char* title = nullptr);
0159 
0160     /** Copy constructor */
0161     HistoND(const HistoND&);
0162 
0163     /** 
0164         // Assignment operator. Works even when the binning of the two
0165         // histograms is not compatible.
0166         */
0167     HistoND& operator=(const HistoND&);
0168 
0169     /** Default constructor not implemented **/
0170     HistoND() = delete;
0171 
0172     /** Histogram dimensionality */
0173     inline unsigned dim() const { return dim_; }
0174 
0175     /** Histogram title */
0176     inline const std::string& title() const { return title_; }
0177 
0178     /** Label associated with accumulated data */
0179     inline const std::string& accumulatedDataLabel() const { return accumulatedDataLabel_; }
0180 
0181     /** Retrive a reference to the array of bin contents */
0182     inline const ArrayND<Numeric>& binContents() const { return data_; }
0183 
0184     /** Retrive a reference to the array of overflows */
0185     inline const ArrayND<Numeric>& overflows() const { return overflow_; }
0186 
0187     /** Inspect histogram axes */
0188     inline const std::vector<Axis>& axes() const { return axes_; }
0189 
0190     /** Inspect a histogram axis for the given dimension */
0191     inline const Axis& axis(const unsigned i) const { return axes_.at(i); }
0192 
0193     /** Total number of bins */
0194     inline unsigned long nBins() const { return data_.length(); }
0195 
0196     /** Total number of fills performed */
0197     inline unsigned long nFillsTotal() const { return fillCount_; }
0198 
0199     /** Total number of fills which fell inside the histogram range */
0200     inline unsigned long nFillsInRange() const { return fillCount_ - overCount_; }
0201 
0202     /** Total number of fills which fell outside the histogram range */
0203     inline unsigned long nFillsOver() const { return overCount_; }
0204 
0205     /** 
0206         // This method returns "true" if the method isUniform()
0207         // of each histogram axis returns "true" 
0208         */
0209     bool isUniformlyBinned() const;
0210 
0211     /** Modify the histogram title */
0212     inline void setTitle(const char* newtitle) {
0213       title_ = newtitle ? newtitle : "";
0214       ++modCount_;
0215     }
0216 
0217     /** Modify the label associated with accumulated data */
0218     inline void setAccumulatedDataLabel(const char* newlabel) {
0219       accumulatedDataLabel_ = newlabel ? newlabel : "";
0220       ++modCount_;
0221     }
0222 
0223     /** Modify the label for the histogram axis with the given number */
0224     inline void setAxisLabel(const unsigned axisNum, const char* newlabel) {
0225       axes_.at(axisNum).setLabel(newlabel);
0226       ++modCount_;
0227     }
0228 
0229     /**
0230         // This method returns width/area/volume/etc. of a single bin.
0231         // 1.0 is returned for a dimensionless histogram.
0232         */
0233     double binVolume(unsigned long binNumber = 0) const;
0234 
0235     /**
0236         // Position of the bin center. Length of the "coords" array
0237         // (filled on return) should be equal to the dimensionality
0238         // of the histogram.
0239         */
0240     void binCenter(unsigned long binNumber, double* coords, unsigned lenCoords) const;
0241 
0242     /**
0243         // Convenience function which fills out a vector of bin centers
0244         // in the same order as the linear order of binContents().
0245         // The class "Point" must have a subscript operator, default
0246         // constructor, copy constructor, and the size() method (use,
0247         // for example, std::array).
0248         */
0249     template <class Point>
0250     void allBinCenters(std::vector<Point>* centers) const;
0251 
0252     /** Bounding box for the given bin */
0253     void binBox(unsigned long binNumber, BoxND<double>* box) const;
0254 
0255     /** Bounding box for the whole histogram */
0256     BoxND<double> boundingBox() const;
0257 
0258     /**
0259         // Volume of the histogram bounding box (this direct call is faster
0260         // than calling boundingBox().volume() ). This function returns 1.0
0261         // for 0-dim histogram, axis interval length for 1-d histogram, etc.
0262         */
0263     double volume() const;
0264 
0265     /** Integral of the histogram */
0266     double integral() const;
0267 
0268     /** Clear the histogram contents (both bins and overflows) */
0269     void clear();
0270 
0271     /** This method clears the bin contents but not overflows */
0272     void clearBinContents();
0273 
0274     /** This method clears overflows but not the bin contents */
0275     void clearOverflows();
0276 
0277     /** Comparison for equality */
0278     bool operator==(const HistoND&) const;
0279 
0280     /** Logical negation of operator== */
0281     bool operator!=(const HistoND&) const;
0282 
0283     /** 
0284         // Check data for equality (both bin contents and overflows).
0285         // Do not compare axes, labels, fill counts, etc.
0286         */
0287     bool isSameData(const HistoND&) const;
0288 
0289     /**
0290         // Fill function for histograms of arbitrary dimensionality.
0291         // The length of the "coords" array should be equal to the
0292         // histogram dimensionality. The Numeric type must have the "+="
0293         // operator defined with the Num2 type on the right side.
0294         */
0295     template <typename Num2>
0296     void fill(const double* coords, unsigned coordLength, const Num2& weight);
0297 
0298     //@{
0299     /**
0300         // Convenience "fill" method for histograms of corresponding
0301         // dimensionality
0302         */
0303     template <typename Num2>
0304     void fill(const Num2& weight);
0305 
0306     template <typename Num2>
0307     void fill(double x0, const Num2& weight);
0308 
0309     template <typename Num2>
0310     void fill(double x0, double x1, const Num2& weight);
0311 
0312     template <typename Num2>
0313     void fill(double x0, double x1, double x2, const Num2& weight);
0314 
0315     template <typename Num2>
0316     void fill(double x0, double x1, double x2, double x3, const Num2& weight);
0317 
0318     template <typename Num2>
0319     void fill(double x0, double x1, double x2, double x3, double x4, const Num2& weight);
0320 
0321     template <typename Num2>
0322     void fill(double x0, double x1, double x2, double x3, double x4, double x5, const Num2& weight);
0323 
0324     template <typename Num2>
0325     void fill(double x0, double x1, double x2, double x3, double x4, double x5, double x6, const Num2& weight);
0326 
0327     template <typename Num2>
0328     void fill(
0329         double x0, double x1, double x2, double x3, double x4, double x5, double x6, double x7, const Num2& weight);
0330 
0331     template <typename Num2>
0332     void fill(double x0,
0333               double x1,
0334               double x2,
0335               double x3,
0336               double x4,
0337               double x5,
0338               double x6,
0339               double x7,
0340               double x8,
0341               const Num2& weight);
0342 
0343     template <typename Num2>
0344     void fill(double x0,
0345               double x1,
0346               double x2,
0347               double x3,
0348               double x4,
0349               double x5,
0350               double x6,
0351               double x7,
0352               double x8,
0353               double x9,
0354               const Num2& weight);
0355     //@}
0356 
0357     /**
0358         // Location-based dispatch method. The provided binary functor
0359         // will be called with the approprite histogram bin value as the
0360         // first argument and the weight as the second (functor return value
0361         // is ignored). This allows for a very general use of the histogram
0362         // binning functionality. For example, with a proper functor, the
0363         // histogram bins can be filled with pointers to an arbitrary class
0364         // (this is the only way to use classes which do not have default
0365         // constructors as bin contents) and the functor can be used to
0366         // dispatch class methods. Depending on the exact nature of the
0367         // functor, multiple things might be modified as the result of this
0368         // call: the bin value, the weight, and the functor internal state.
0369         */
0370     template <typename Num2, class Functor>
0371     void dispatch(const double* coords, unsigned coordLength, Num2& weight, Functor& f);
0372 
0373     //@{
0374     /**
0375         // Convenience "dispatch" method for histograms of corresponding
0376         // dimensionality
0377         */
0378     template <typename Num2, class Functor>
0379     void dispatch(Num2& weight, Functor& f);
0380 
0381     template <typename Num2, class Functor>
0382     void dispatch(double x0, Num2& weight, Functor& f);
0383 
0384     template <typename Num2, class Functor>
0385     void dispatch(double x0, double x1, Num2& weight, Functor& f);
0386 
0387     template <typename Num2, class Functor>
0388     void dispatch(double x0, double x1, double x2, Num2& weight, Functor& f);
0389 
0390     template <typename Num2, class Functor>
0391     void dispatch(double x0, double x1, double x2, double x3, Num2& weight, Functor& f);
0392 
0393     template <typename Num2, class Functor>
0394     void dispatch(double x0, double x1, double x2, double x3, double x4, Num2& weight, Functor& f);
0395 
0396     template <typename Num2, class Functor>
0397     void dispatch(double x0, double x1, double x2, double x3, double x4, double x5, Num2& weight, Functor& f);
0398 
0399     template <typename Num2, class Functor>
0400     void dispatch(double x0, double x1, double x2, double x3, double x4, double x5, double x6, Num2& weight, Functor& f);
0401 
0402     template <typename Num2, class Functor>
0403     void dispatch(double x0,
0404                   double x1,
0405                   double x2,
0406                   double x3,
0407                   double x4,
0408                   double x5,
0409                   double x6,
0410                   double x7,
0411                   Num2& weight,
0412                   Functor& f);
0413 
0414     template <typename Num2, class Functor>
0415     void dispatch(double x0,
0416                   double x1,
0417                   double x2,
0418                   double x3,
0419                   double x4,
0420                   double x5,
0421                   double x6,
0422                   double x7,
0423                   double x8,
0424                   Num2& weight,
0425                   Functor& f);
0426 
0427     template <typename Num2, class Functor>
0428     void dispatch(double x0,
0429                   double x1,
0430                   double x2,
0431                   double x3,
0432                   double x4,
0433                   double x5,
0434                   double x6,
0435                   double x7,
0436                   double x8,
0437                   double x9,
0438                   Num2& weight,
0439                   Functor& f);
0440     //@}
0441 
0442     /**
0443         // The "examine" functions allow the user to access bin contents
0444         // when bins are addressed by their coordinates. Use "binContents()"
0445         // to access the data by bin numbers. Overflow bins will be accessed
0446         // if the given coordinates fall outside the histogram range.
0447         */
0448     const Numeric& examine(const double* coords, unsigned coordLength) const;
0449 
0450     //@{
0451     /**
0452         // Convenience "examine" method for histograms of corresponding
0453         // dimensionality
0454         */
0455     const Numeric& examine() const;
0456 
0457     const Numeric& examine(double x0) const;
0458 
0459     const Numeric& examine(double x0, double x1) const;
0460 
0461     const Numeric& examine(double x0, double x1, double x2) const;
0462 
0463     const Numeric& examine(double x0, double x1, double x2, double x3) const;
0464 
0465     const Numeric& examine(double x0, double x1, double x2, double x3, double x4) const;
0466 
0467     const Numeric& examine(double x0, double x1, double x2, double x3, double x4, double x5) const;
0468 
0469     const Numeric& examine(double x0, double x1, double x2, double x3, double x4, double x5, double x6) const;
0470 
0471     const Numeric& examine(double x0, double x1, double x2, double x3, double x4, double x5, double x6, double x7) const;
0472 
0473     const Numeric& examine(
0474         double x0, double x1, double x2, double x3, double x4, double x5, double x6, double x7, double x8) const;
0475 
0476     const Numeric& examine(
0477         double x0, double x1, double x2, double x3, double x4, double x5, double x6, double x7, double x8, double x9)
0478         const;
0479     //@}
0480 
0481     /**
0482         // The "closestBin" functions are similar to the "examine" functions
0483         // but always return a valid bin and never overflow. This can be
0484         // useful for implementing lookup tables with constant extrapolation
0485         // outside of the histogram range.
0486         */
0487     const Numeric& closestBin(const double* coords, unsigned coordLength) const;
0488 
0489     //@{
0490     /**
0491         // Convenience "closestBin" method for histograms of corresponding
0492         // dimensionality
0493         */
0494     const Numeric& closestBin() const;
0495 
0496     const Numeric& closestBin(double x0) const;
0497 
0498     const Numeric& closestBin(double x0, double x1) const;
0499 
0500     const Numeric& closestBin(double x0, double x1, double x2) const;
0501 
0502     const Numeric& closestBin(double x0, double x1, double x2, double x3) const;
0503 
0504     const Numeric& closestBin(double x0, double x1, double x2, double x3, double x4) const;
0505 
0506     const Numeric& closestBin(double x0, double x1, double x2, double x3, double x4, double x5) const;
0507 
0508     const Numeric& closestBin(double x0, double x1, double x2, double x3, double x4, double x5, double x6) const;
0509 
0510     const Numeric& closestBin(
0511         double x0, double x1, double x2, double x3, double x4, double x5, double x6, double x7) const;
0512 
0513     const Numeric& closestBin(
0514         double x0, double x1, double x2, double x3, double x4, double x5, double x6, double x7, double x8) const;
0515 
0516     const Numeric& closestBin(
0517         double x0, double x1, double x2, double x3, double x4, double x5, double x6, double x7, double x8, double x9)
0518         const;
0519     //@}
0520 
0521     /**
0522         // The "fillC" functions are similar to the "fill" methods but
0523         // they preserve the centroid of the deposit. Note that, if the
0524         // histogram dimensionality is high, "fillC" works significantly
0525         // slower than the corresponding "fill". Also note that there
0526         // must be at least 2 bins in each dimension in order for this
0527         // function to work.
0528         //
0529         // A word of caution. What is added to the bins is the input weight
0530         // multiplied by another weight calculated using the bin proximity.
0531         // If the input weight is just 1 (which happens quite often in
0532         // practice), the product of the weights is normally less than 1.
0533         // If the histogram template parameter is one of the integer types,
0534         // operator += will convert this product to 0 before adding it to
0535         // the bin! Therefore, it is best to use "fillC" only with floating
0536         // point template parameters (float, double, etc).
0537         // 
0538         // Currently, the "fillC" methods work sensibly only in the case
0539         // the binning is uniform (i.e., the second template parameter is
0540         // HistoAxis rather than, let say, NUHistoAxis). They typically
0541         // will not even compile if the binning is not uniform.
0542         */
0543     template <typename Num2>
0544     void fillC(const double* coords, unsigned coordLength, const Num2& weight);
0545 
0546     //@{
0547     /**
0548         // Convenience "fillC" method for histograms of corresponding
0549         // dimensionality
0550         */
0551     template <typename Num2>
0552     void fillC(const Num2& weight);
0553 
0554     template <typename Num2>
0555     void fillC(double x0, const Num2& weight);
0556 
0557     template <typename Num2>
0558     void fillC(double x0, double x1, const Num2& weight);
0559 
0560     template <typename Num2>
0561     void fillC(double x0, double x1, double x2, const Num2& weight);
0562 
0563     template <typename Num2>
0564     void fillC(double x0, double x1, double x2, double x3, const Num2& weight);
0565 
0566     template <typename Num2>
0567     void fillC(double x0, double x1, double x2, double x3, double x4, const Num2& weight);
0568 
0569     template <typename Num2>
0570     void fillC(double x0, double x1, double x2, double x3, double x4, double x5, const Num2& weight);
0571 
0572     template <typename Num2>
0573     void fillC(double x0, double x1, double x2, double x3, double x4, double x5, double x6, const Num2& weight);
0574 
0575     template <typename Num2>
0576     void fillC(
0577         double x0, double x1, double x2, double x3, double x4, double x5, double x6, double x7, const Num2& weight);
0578 
0579     template <typename Num2>
0580     void fillC(double x0,
0581                double x1,
0582                double x2,
0583                double x3,
0584                double x4,
0585                double x5,
0586                double x6,
0587                double x7,
0588                double x8,
0589                const Num2& weight);
0590 
0591     template <typename Num2>
0592     void fillC(double x0,
0593                double x1,
0594                double x2,
0595                double x3,
0596                double x4,
0597                double x5,
0598                double x6,
0599                double x7,
0600                double x8,
0601                double x9,
0602                const Num2& weight);
0603     //@}
0604 
0605     /**
0606         // Fill from another histogram. Compatibility of axis limits
0607         // will not be checked, but compatibility of array shapes will be.
0608         */
0609     template <typename Num2>
0610     HistoND& operator+=(const HistoND<Num2, Axis>& r);
0611 
0612     /**
0613         // Subtract contents of another histogram. Equivalent to multiplying
0614         // the contents of the other histogram by -1 and then adding them.
0615         // One of the consequences of this approach is that, for histograms
0616         // "a" and "b", the sequence of operations "a += b; a -= b;" does not
0617         // leave histogram "a" unchanged: although its bin contents will
0618         // remain the same (up to round-off errors), the fill counts will
0619         // increase by twice the fill counts of "b".
0620         */
0621     template <typename Num2>
0622     HistoND& operator-=(const HistoND<Num2, Axis>& r);
0623 
0624     //@{
0625     /** Method to set contents of individual bins (no bounds checking) */
0626     template <typename Num2>
0627     void setBin(const unsigned* index, unsigned indexLen, const Num2& v);
0628 
0629     template <typename Num2>
0630     void setBin(const Num2& v);
0631 
0632     template <typename Num2>
0633     void setBin(unsigned i0, const Num2& v);
0634 
0635     template <typename Num2>
0636     void setBin(unsigned i0, unsigned i1, const Num2& v);
0637 
0638     template <typename Num2>
0639     void setBin(unsigned i0, unsigned i1, unsigned i2, const Num2& v);
0640 
0641     template <typename Num2>
0642     void setBin(unsigned i0, unsigned i1, unsigned i2, unsigned i3, const Num2& v);
0643 
0644     template <typename Num2>
0645     void setBin(unsigned i0, unsigned i1, unsigned i2, unsigned i3, unsigned i4, const Num2& v);
0646 
0647     template <typename Num2>
0648     void setBin(unsigned i0, unsigned i1, unsigned i2, unsigned i3, unsigned i4, unsigned i5, const Num2& v);
0649 
0650     template <typename Num2>
0651     void setBin(
0652         unsigned i0, unsigned i1, unsigned i2, unsigned i3, unsigned i4, unsigned i5, unsigned i6, const Num2& v);
0653 
0654     template <typename Num2>
0655     void setBin(unsigned i0,
0656                 unsigned i1,
0657                 unsigned i2,
0658                 unsigned i3,
0659                 unsigned i4,
0660                 unsigned i5,
0661                 unsigned i6,
0662                 unsigned i7,
0663                 const Num2& v);
0664 
0665     template <typename Num2>
0666     void setBin(unsigned i0,
0667                 unsigned i1,
0668                 unsigned i2,
0669                 unsigned i3,
0670                 unsigned i4,
0671                 unsigned i5,
0672                 unsigned i6,
0673                 unsigned i7,
0674                 unsigned i8,
0675                 const Num2& v);
0676 
0677     template <typename Num2>
0678     void setBin(unsigned i0,
0679                 unsigned i1,
0680                 unsigned i2,
0681                 unsigned i3,
0682                 unsigned i4,
0683                 unsigned i5,
0684                 unsigned i6,
0685                 unsigned i7,
0686                 unsigned i8,
0687                 unsigned i9,
0688                 const Num2& v);
0689 
0690     template <typename Num2>
0691     inline void setLinearBin(const unsigned long index, const Num2& v) {
0692       data_.linearValue(index) = v;
0693       ++modCount_;
0694     }
0695     //@}
0696 
0697     //@{
0698     /** Method to set contents of individual bins with bounds checking */
0699     template <typename Num2>
0700     void setBinAt(const unsigned* index, unsigned indexLen, const Num2& v);
0701 
0702     template <typename Num2>
0703     void setBinAt(const Num2& v);
0704 
0705     template <typename Num2>
0706     void setBinAt(unsigned i0, const Num2& v);
0707 
0708     template <typename Num2>
0709     void setBinAt(unsigned i0, unsigned i1, const Num2& v);
0710 
0711     template <typename Num2>
0712     void setBinAt(unsigned i0, unsigned i1, unsigned i2, const Num2& v);
0713 
0714     template <typename Num2>
0715     void setBinAt(unsigned i0, unsigned i1, unsigned i2, unsigned i3, const Num2& v);
0716 
0717     template <typename Num2>
0718     void setBinAt(unsigned i0, unsigned i1, unsigned i2, unsigned i3, unsigned i4, const Num2& v);
0719 
0720     template <typename Num2>
0721     void setBinAt(unsigned i0, unsigned i1, unsigned i2, unsigned i3, unsigned i4, unsigned i5, const Num2& v);
0722 
0723     template <typename Num2>
0724     void setBinAt(
0725         unsigned i0, unsigned i1, unsigned i2, unsigned i3, unsigned i4, unsigned i5, unsigned i6, const Num2& v);
0726 
0727     template <typename Num2>
0728     void setBinAt(unsigned i0,
0729                   unsigned i1,
0730                   unsigned i2,
0731                   unsigned i3,
0732                   unsigned i4,
0733                   unsigned i5,
0734                   unsigned i6,
0735                   unsigned i7,
0736                   const Num2& v);
0737 
0738     template <typename Num2>
0739     void setBinAt(unsigned i0,
0740                   unsigned i1,
0741                   unsigned i2,
0742                   unsigned i3,
0743                   unsigned i4,
0744                   unsigned i5,
0745                   unsigned i6,
0746                   unsigned i7,
0747                   unsigned i8,
0748                   const Num2& v);
0749 
0750     template <typename Num2>
0751     void setBinAt(unsigned i0,
0752                   unsigned i1,
0753                   unsigned i2,
0754                   unsigned i3,
0755                   unsigned i4,
0756                   unsigned i5,
0757                   unsigned i6,
0758                   unsigned i7,
0759                   unsigned i8,
0760                   unsigned i9,
0761                   const Num2& v);
0762 
0763     template <typename Num2>
0764     inline void setLinearBinAt(const unsigned long index, const Num2& v) {
0765       data_.linearValueAt(index) = v;
0766       ++modCount_;
0767     }
0768     //@}
0769 
0770     /** This method sets all bin contents in one fell swoop */
0771     template <typename Num2>
0772     void setBinContents(const Num2* data, unsigned long dataLength, bool clearOverflows = true);
0773 
0774     /** This method sets all overflows in one fell swoop */
0775     template <typename Num2>
0776     void setOverflows(const Num2* data, unsigned long dataLength);
0777 
0778     /**
0779         // Setting bin contents to some constant value.
0780         // The Num2 type should allow automatic conversion to Numeric.
0781         */
0782     template <typename Num2>
0783     inline void setBinsToConst(const Num2& value) {
0784       data_.constFill(value);
0785       ++modCount_;
0786     }
0787 
0788     /**
0789         // Setting overflows to some constant value.
0790         // The Num2 type should allow automatic conversion to Numeric.
0791         */
0792     template <typename Num2>
0793     inline void setOverflowsToConst(const Num2& value) {
0794       overflow_.constFill(value);
0795       ++modCount_;
0796     }
0797 
0798     /**
0799         // This member function instructs the histogram to recalculate
0800         // the number of fills from data. It may be useful to call this
0801         // function after "setBinContents" in case the contents are filled
0802         // from another histogram.
0803         */
0804     void recalculateNFillsFromData();
0805 
0806     //@{
0807     /**
0808         // This method is intended for data format conversion
0809         // programs only, not for typical histogramming use
0810         */
0811     inline void setNFillsTotal(const unsigned long i) {
0812       fillCount_ = i;
0813       ++modCount_;
0814     }
0815     inline void setNFillsOver(const unsigned long i) {
0816       overCount_ = i;
0817       ++modCount_;
0818     }
0819     //@}
0820 
0821     /** In-place multiplication by a scalar (scaling) */
0822     template <typename Num2>
0823     HistoND& operator*=(const Num2& r);
0824 
0825     /** In-place division by a scalar */
0826     template <typename Num2>
0827     HistoND& operator/=(const Num2& r);
0828 
0829     //@{
0830     /** Multiplication by a value which is different for every bin */
0831     template <typename Num2>
0832     void scaleBinContents(const Num2* data, unsigned long dataLength);
0833 
0834     template <typename Num2>
0835     void scaleOverflows(const Num2* data, unsigned long dataLength);
0836     //@}
0837 
0838     //@{
0839     /**
0840         // In-place addition of a scalar to all bins. Equivalent to calling
0841         // the "fill" function with the same weight once for every bin.
0842         */
0843     template <typename Num2>
0844     void addToBinContents(const Num2& weight);
0845 
0846     template <typename Num2>
0847     void addToOverflows(const Num2& weight);
0848     //@}
0849 
0850     //@{
0851     /**
0852         // In-place addition of an array. Equivalent to calling the "fill"
0853         // function once for every bin with the weight taken from the
0854         // corresponding array element.
0855         */
0856     template <typename Num2>
0857     void addToBinContents(const Num2* data, unsigned long dataLength);
0858 
0859     template <typename Num2>
0860     void addToOverflows(const Num2* data, unsigned long dataLength);
0861     //@}
0862 
0863     /**
0864         // Add contents of all bins inside the given box to the accumulator.
0865         // Note that Numeric type must support multiplication by a double
0866         // in order for this function to work (it calculates the overlap
0867         // fraction of each bin with the box and multiplies bin content
0868         // by that fraction for subsequent accumulation). The operation
0869         // Acc += Numeric must be defined.
0870         */
0871     template <typename Acc>
0872     void accumulateBinsInBox(const BoxND<double>& box, Acc* acc, bool calculateAverage = false) const;
0873 
0874     //@{
0875     /**
0876         // Code for projecting one histogram onto another. For now,
0877         // this is done for bin contents only, not for overflows.
0878         // The projection should be created in advance from this
0879         // histogram with the aid of the slicing constructor. The indices
0880         // used in that constructor should be provided here as well.
0881         //
0882         // Note that you might want to recalculate the number of fills
0883         // from data after performing all projections needed.
0884         */
0885     template <typename Num2, typename Num3>
0886     void addToProjection(HistoND<Num2, Axis>* projection,
0887                          AbsArrayProjector<Numeric, Num3>& projector,
0888                          const unsigned* projectedIndices,
0889                          unsigned nProjectedIndices) const;
0890 
0891     template <typename Num2, typename Num3>
0892     void addToProjection(HistoND<Num2, Axis>* projection,
0893                          AbsVisitor<Numeric, Num3>& projector,
0894                          const unsigned* projectedIndices,
0895                          unsigned nProjectedIndices) const;
0896     //@}
0897 
0898     /** Transpose the histogram axes and bin contents */
0899     HistoND transpose(unsigned axisNum1, unsigned axisNum2) const;
0900 
0901     /**
0902         // This method returns the number of modifications
0903         // performed on the histogram since its creation. This number
0904         // is always increasing during the lifetime of the histogram
0905         // object. Its main property is as follows: if the method
0906         // "getModCount" returns the same number twice, there should
0907         // be no changes in the histogram object (so that a drawing
0908         // program does not need to redraw the histogram image).
0909         //
0910         // This number is pure transient, it is not serialized and
0911         // does not participate in histogram comparisons for equality.
0912         */
0913     inline unsigned long getModCount() const { return modCount_; }
0914 
0915     /**
0916         // Indicate that the histogram contents have changed. Should
0917         // be used by any code which directly modifies histogram bins
0918         // (after using const_cast on the relevant reference).
0919         */
0920     inline void incrModCount() { ++modCount_; }
0921 
0922     //@{
0923     /** Method related to "geners" I/O */
0924     inline gs::ClassId classId() const { return gs::ClassId(*this); }
0925     bool write(std::ostream& of) const;
0926     //@}
0927 
0928     static const char* classname();
0929     static inline unsigned version() { return 1; }
0930     static HistoND* read(const gs::ClassId& id, std::istream& in);
0931 
0932   private:
0933     // Special constructor which speeds up the "transpose" operation.
0934     // Does not do full error checking (some of it is done in transpose).
0935     HistoND(const HistoND& r, unsigned ax1, unsigned ax2);
0936 
0937     template <typename Num2>
0938     void fillPreservingCentroid(const Num2& weight);
0939 
0940     template <typename Acc>
0941     void accumulateBinsLoop(unsigned level,
0942                             const BoxND<double>& box,
0943                             unsigned* idx,
0944                             Acc* accumulator,
0945                             double overlapFraction,
0946                             long double* wsum) const;
0947     std::string title_;
0948     std::string accumulatedDataLabel_;
0949     ArrayND<Numeric> data_;
0950     ArrayND<Numeric> overflow_;
0951     std::vector<Axis> axes_;
0952     mutable std::vector<double> weightBuf_;
0953     mutable std::vector<unsigned> indexBuf_;
0954     unsigned long fillCount_;
0955     unsigned long overCount_;
0956     unsigned long modCount_;
0957     unsigned dim_;
0958   };
0959 
0960   /**
0961     // Reset negative histogram bins to zero and then divide histogram
0962     // bin contents by the histogram integral. If the "knownNonNegative"
0963     // argument is true, it will be assumed that there are no negative
0964     // bins, and their explicit reset is unnecessary.
0965     //
0966     // This function will throw npstat::NpstatRuntimeError in case the histogram
0967     // is empty after all negative bins are reset.
0968     //
0969     // This function is not a member of the HistoND class itself because
0970     // these operations do not necessarily make sense for all bin types.
0971     // Making such operation a member would make creation of HistoND
0972     // scripting API (e.g., for python) more difficult.
0973     */
0974   template <typename Histo>
0975   void convertHistoToDensity(Histo* histogram, bool knownNonNegative = false);
0976 
0977   /**
0978     // Generate a density scanning map for subsequent use with
0979     // the "DensityScanND" template. Naturally, only histograms
0980     // with uniform binning can be used here.
0981     */
0982   template <typename Histo>
0983   std::vector<LinearMapper1d> densityScanHistoMap(const Histo& histo);
0984 
0985   /**
0986     // Generate a density scanning map for subsequent use with the
0987     // "DensityScanND" template when a density is to be convolved with
0988     // the histogram data. Only histograms with uniform binning
0989     // can be used here.
0990     //
0991     // The "doubleDataRange" should be set "true" in case the data
0992     // will be mirrored (or just empty range added) to avoid circular
0993     // spilling after convolution.
0994     */
0995   template <typename Histo>
0996   std::vector<CircularMapper1d> convolutionHistoMap(const Histo& histo, bool doubleDataRange);
0997 }  // namespace npstat
0998 
0999 #include <cassert>
1000 #include "JetMETCorrections/InterpolationTables/interface/NpstatException.h"
1001 #include <sstream>
1002 #include <climits>
1003 #include <algorithm>
1004 
1005 #include "Alignment/Geners/interface/CPP11_auto_ptr.hh"
1006 #include "Alignment/Geners/interface/binaryIO.hh"
1007 
1008 namespace npstat {
1009   namespace Private {
1010     template <class Axis>
1011     ArrayShape makeHistoShape(const std::vector<Axis>& axes) {
1012       const unsigned n = axes.size();
1013       ArrayShape result;
1014       result.reserve(n);
1015       for (unsigned i = 0; i < n; ++i)
1016         result.push_back(axes[i].nBins());
1017       return result;
1018     }
1019 
1020     template <class Axis>
1021     ArrayShape makeHistoShape(const Axis& xAxis) {
1022       ArrayShape result;
1023       result.reserve(1U);
1024       result.push_back(xAxis.nBins());
1025       return result;
1026     }
1027 
1028     template <class Axis>
1029     ArrayShape makeHistoShape(const Axis& xAxis, const Axis& yAxis) {
1030       ArrayShape result;
1031       result.reserve(2U);
1032       result.push_back(xAxis.nBins());
1033       result.push_back(yAxis.nBins());
1034       return result;
1035     }
1036 
1037     template <class Axis>
1038     ArrayShape makeHistoShape(const Axis& xAxis, const Axis& yAxis, const Axis& zAxis) {
1039       ArrayShape result;
1040       result.reserve(3U);
1041       result.push_back(xAxis.nBins());
1042       result.push_back(yAxis.nBins());
1043       result.push_back(zAxis.nBins());
1044       return result;
1045     }
1046 
1047     template <class Axis>
1048     ArrayShape makeHistoShape(const Axis& xAxis, const Axis& yAxis, const Axis& zAxis, const Axis& tAxis) {
1049       ArrayShape result;
1050       result.reserve(4U);
1051       result.push_back(xAxis.nBins());
1052       result.push_back(yAxis.nBins());
1053       result.push_back(zAxis.nBins());
1054       result.push_back(tAxis.nBins());
1055       return result;
1056     }
1057 
1058     template <class Axis>
1059     ArrayShape makeHistoShape(
1060         const Axis& xAxis, const Axis& yAxis, const Axis& zAxis, const Axis& tAxis, const Axis& vAxis) {
1061       ArrayShape result;
1062       result.reserve(5U);
1063       result.push_back(xAxis.nBins());
1064       result.push_back(yAxis.nBins());
1065       result.push_back(zAxis.nBins());
1066       result.push_back(tAxis.nBins());
1067       result.push_back(vAxis.nBins());
1068       return result;
1069     }
1070 
1071     template <class Axis>
1072     std::vector<Axis> rebinAxes(const std::vector<Axis>& axes, const unsigned* newBins, const unsigned lenNewBins) {
1073       const unsigned dim = axes.size();
1074       if (lenNewBins != dim)
1075         throw npstat::NpstatInvalidArgument(
1076             "In npstat::Private::rebinAxes: invalid length "
1077             "of the new bins array");
1078       assert(newBins);
1079       std::vector<Axis> newAxes;
1080       newAxes.reserve(dim);
1081       for (unsigned i = 0; i < dim; ++i)
1082         newAxes.push_back(axes[i].rebin(newBins[i]));
1083       return newAxes;
1084     }
1085 
1086     template <class Axis>
1087     std::vector<Axis> axesOfASlice(const std::vector<Axis>& axes,
1088                                    const unsigned* fixedIndices,
1089                                    const unsigned nFixedIndices) {
1090       const unsigned dim = axes.size();
1091       std::vector<Axis> newAxes;
1092       if (nFixedIndices == 0U)
1093         throw npstat::NpstatInvalidArgument(
1094             "In npstat::Private::axesOfASlice: "
1095             "at least one fixed index must be specified");
1096       if (nFixedIndices > dim)
1097         throw npstat::NpstatInvalidArgument("In npstat::Private::axesOfASlice: too many fixed indices");
1098       assert(fixedIndices);
1099       for (unsigned i = 0; i < nFixedIndices; ++i)
1100         if (fixedIndices[i] >= dim)
1101           throw npstat::NpstatInvalidArgument("In npstat::Private::axesOfASlice: fixed index out of range");
1102       newAxes.reserve(dim - nFixedIndices);
1103       for (unsigned i = 0; i < dim; ++i) {
1104         bool fixed = false;
1105         for (unsigned j = 0; j < nFixedIndices; ++j)
1106           if (fixedIndices[j] == i) {
1107             fixed = true;
1108             break;
1109           }
1110         if (!fixed)
1111           newAxes.push_back(axes[i]);
1112       }
1113       if (newAxes.size() != dim - nFixedIndices)
1114         throw npstat::NpstatInvalidArgument("In npstat::Private::axesOfASlice: duplicate fixed index");
1115       return newAxes;
1116     }
1117 
1118     template <class Axis>
1119     ArrayShape shapeOfASlice(const std::vector<Axis>& axes,
1120                              const unsigned* fixedIndices,
1121                              const unsigned nFixedIndices) {
1122       const unsigned dim = axes.size();
1123       if (nFixedIndices == 0U)
1124         throw npstat::NpstatInvalidArgument(
1125             "In npstat::Private::shapeOfASlice: "
1126             "at least one fixed index must be specified");
1127       if (nFixedIndices > dim)
1128         throw npstat::NpstatInvalidArgument("In npstat::Private::shapeOfASlice: too many fixed indices");
1129       assert(fixedIndices);
1130 
1131       // Check that the fixed indices are within range
1132       for (unsigned j = 0; j < nFixedIndices; ++j)
1133         if (fixedIndices[j] >= dim)
1134           throw npstat::NpstatInvalidArgument("In npstat::Private::shapeOfASlice: fixed index out of range");
1135 
1136       // Build the shape for the slice
1137       ArrayShape sh;
1138       if (nFixedIndices < dim)
1139         sh.reserve(dim - nFixedIndices);
1140       for (unsigned i = 0; i < dim; ++i) {
1141         bool fixed = false;
1142         for (unsigned j = 0; j < nFixedIndices; ++j)
1143           if (fixedIndices[j] == i) {
1144             fixed = true;
1145             break;
1146           }
1147         if (!fixed)
1148           sh.push_back(axes[i].nBins());
1149       }
1150       if (sh.size() != dim - nFixedIndices)
1151         throw npstat::NpstatInvalidArgument("In npstat::Private::shapeOfASlice: duplicate fixed index");
1152       return sh;
1153     }
1154 
1155     template <class Axis>
1156     std::vector<Axis> addAxis(const std::vector<Axis>& axes, const Axis& newAxis, const unsigned newAxisNumber) {
1157       const unsigned dim = axes.size();
1158       std::vector<Axis> newAxes;
1159       newAxes.reserve(dim + 1U);
1160       unsigned iadd = 0;
1161       for (unsigned i = 0; i < dim; ++i) {
1162         if (newAxisNumber == i)
1163           newAxes.push_back(newAxis);
1164         else
1165           newAxes.push_back(axes[iadd++]);
1166       }
1167       if (iadd == dim)
1168         newAxes.push_back(newAxis);
1169       else
1170         newAxes.push_back(axes[iadd]);
1171       return newAxes;
1172     }
1173 
1174     template <class Axis>
1175     ArrayShape shapeWithExtraAxis(const std::vector<Axis>& axes, const Axis& newAxis, const unsigned newAxisNumber) {
1176       const unsigned dim = axes.size();
1177       ArrayShape result;
1178       result.reserve(dim + 1U);
1179       unsigned iadd = 0;
1180       for (unsigned i = 0; i < dim; ++i) {
1181         if (newAxisNumber == i)
1182           result.push_back(newAxis.nBins());
1183         else
1184           result.push_back(axes[iadd++].nBins());
1185       }
1186       if (iadd == dim)
1187         result.push_back(newAxis.nBins());
1188       else
1189         result.push_back(axes[iadd].nBins());
1190       return result;
1191     }
1192 
1193     inline void h_badargs(const char* method) {
1194       std::ostringstream os;
1195       os << "In npstat::HistoND::" << method << ": number of arguments"
1196          << " is incompatible with histogram dimensionality";
1197       throw npstat::NpstatInvalidArgument(os.str());
1198     }
1199   }  // namespace Private
1200 
1201   template <typename Numeric, class Axis>
1202   template <typename Acc>
1203   void HistoND<Numeric, Axis>::accumulateBinsLoop(const unsigned level,
1204                                                   const BoxND<double>& box,
1205                                                   unsigned* idx,
1206                                                   Acc* accumulator,
1207                                                   const double overlapFraction,
1208                                                   long double* wsum) const {
1209     const Interval<double>& boxSide(box[level]);
1210     const Axis& axis(axes_[level]);
1211     const unsigned nbins = axis.nBins();
1212     const bool lastLevel = level == dim_ - 1U;
1213     for (unsigned i = 0; i < nbins; ++i) {
1214       const double over = overlapFraction * axis.binInterval(i).overlapFraction(boxSide);
1215       if (over > 0.0) {
1216         idx[level] = i;
1217         if (lastLevel) {
1218           *accumulator += over * data_.value(idx, dim_);
1219           *wsum += over;
1220         } else
1221           accumulateBinsLoop(level + 1U, box, idx, accumulator, over, wsum);
1222       }
1223     }
1224   }
1225 
1226   template <typename Numeric, class Axis>
1227   template <typename Acc>
1228   void HistoND<Numeric, Axis>::accumulateBinsInBox(const BoxND<double>& box,
1229                                                    Acc* accumulator,
1230                                                    const bool calculateAverage) const {
1231     if (box.size() != dim_)
1232       throw npstat::NpstatInvalidArgument(
1233           "In npstat::HistoND::accumulateBinsInBox: "
1234           "incompatible box dimensionality");
1235     assert(accumulator);
1236     if (dim_) {
1237       long double wsum = 0.0L;
1238       for (unsigned i = 0; i < dim_; ++i)
1239         indexBuf_[i] = 0U;
1240       accumulateBinsLoop(0U, box, &indexBuf_[0], accumulator, 1.0, &wsum);
1241       if (calculateAverage && wsum > 0.0L)
1242         *accumulator *= static_cast<double>(1.0L / wsum);
1243     } else
1244       *accumulator += 1.0 * data_();
1245   }
1246 
1247   template <typename Numeric, class Axis>
1248   inline void HistoND<Numeric, Axis>::clearBinContents() {
1249     data_.clear();
1250     fillCount_ = 0UL;
1251     ++modCount_;
1252   }
1253 
1254   template <typename Numeric, class Axis>
1255   inline void HistoND<Numeric, Axis>::clearOverflows() {
1256     overflow_.clear();
1257     overCount_ = 0UL;
1258     ++modCount_;
1259   }
1260 
1261   template <typename Numeric, class Axis>
1262   inline void HistoND<Numeric, Axis>::clear() {
1263     clearBinContents();
1264     clearOverflows();
1265     ++modCount_;
1266   }
1267 
1268   template <typename Numeric, class Axis>
1269   HistoND<Numeric, Axis>::HistoND(const std::vector<Axis>& axesIn, const char* title, const char* label)
1270       : title_(title ? title : ""),
1271         accumulatedDataLabel_(label ? label : ""),
1272         data_(Private::makeHistoShape(axesIn)),
1273         overflow_(ArrayShape(axesIn.size(), 3U)),
1274         axes_(axesIn),
1275         weightBuf_(axesIn.size()),
1276         indexBuf_(2U * axesIn.size()),
1277         modCount_(0UL),
1278         dim_(axesIn.size()) {
1279     if (dim_ >= CHAR_BIT * sizeof(unsigned long))
1280       throw npstat::NpstatInvalidArgument(
1281           "In npstat::HistoND constructor: requested histogram "
1282           "dimensionality is not supported (too large)");
1283     clear();
1284   }
1285 
1286   template <typename Numeric, class Axis>
1287   HistoND<Numeric, Axis>::HistoND(const Axis& xAxis, const char* title, const char* label)
1288       : title_(title ? title : ""),
1289         accumulatedDataLabel_(label ? label : ""),
1290         data_(Private::makeHistoShape(xAxis)),
1291         overflow_(ArrayShape(1U, 3U)),
1292         weightBuf_(1U),
1293         indexBuf_(2U * 1U),
1294         modCount_(0UL),
1295         dim_(1U) {
1296     axes_.reserve(dim_);
1297     axes_.push_back(xAxis);
1298     clear();
1299   }
1300 
1301   template <typename Numeric, class Axis>
1302   HistoND<Numeric, Axis>::HistoND(const Axis& xAxis, const Axis& yAxis, const char* title, const char* label)
1303       : title_(title ? title : ""),
1304         accumulatedDataLabel_(label ? label : ""),
1305         data_(Private::makeHistoShape(xAxis, yAxis)),
1306         overflow_(ArrayShape(2U, 3U)),
1307         weightBuf_(2U),
1308         indexBuf_(2U * 2U),
1309         modCount_(0UL),
1310         dim_(2U) {
1311     axes_.reserve(dim_);
1312     axes_.push_back(xAxis);
1313     axes_.push_back(yAxis);
1314     clear();
1315   }
1316 
1317   template <typename Numeric, class Axis>
1318   HistoND<Numeric, Axis>::HistoND(
1319       const Axis& xAxis, const Axis& yAxis, const Axis& zAxis, const char* title, const char* label)
1320       : title_(title ? title : ""),
1321         accumulatedDataLabel_(label ? label : ""),
1322         data_(Private::makeHistoShape(xAxis, yAxis, zAxis)),
1323         overflow_(ArrayShape(3U, 3U)),
1324         weightBuf_(3U),
1325         indexBuf_(2U * 3U),
1326         modCount_(0UL),
1327         dim_(3U) {
1328     axes_.reserve(dim_);
1329     axes_.push_back(xAxis);
1330     axes_.push_back(yAxis);
1331     axes_.push_back(zAxis);
1332     clear();
1333   }
1334 
1335   template <typename Numeric, class Axis>
1336   HistoND<Numeric, Axis>::HistoND(
1337       const Axis& xAxis, const Axis& yAxis, const Axis& zAxis, const Axis& tAxis, const char* title, const char* label)
1338       : title_(title ? title : ""),
1339         accumulatedDataLabel_(label ? label : ""),
1340         data_(Private::makeHistoShape(xAxis, yAxis, zAxis, tAxis)),
1341         overflow_(ArrayShape(4U, 3U)),
1342         weightBuf_(4U),
1343         indexBuf_(2U * 4U),
1344         modCount_(0UL),
1345         dim_(4U) {
1346     axes_.reserve(dim_);
1347     axes_.push_back(xAxis);
1348     axes_.push_back(yAxis);
1349     axes_.push_back(zAxis);
1350     axes_.push_back(tAxis);
1351     clear();
1352   }
1353 
1354   template <typename Numeric, class Axis>
1355   HistoND<Numeric, Axis>::HistoND(const Axis& xAxis,
1356                                   const Axis& yAxis,
1357                                   const Axis& zAxis,
1358                                   const Axis& tAxis,
1359                                   const Axis& vAxis,
1360                                   const char* title,
1361                                   const char* label)
1362       : title_(title ? title : ""),
1363         accumulatedDataLabel_(label ? label : ""),
1364         data_(Private::makeHistoShape(xAxis, yAxis, zAxis, tAxis, vAxis)),
1365         overflow_(ArrayShape(5U, 3U)),
1366         weightBuf_(5U),
1367         indexBuf_(2U * 5U),
1368         modCount_(0UL),
1369         dim_(5U) {
1370     axes_.reserve(dim_);
1371     axes_.push_back(xAxis);
1372     axes_.push_back(yAxis);
1373     axes_.push_back(zAxis);
1374     axes_.push_back(tAxis);
1375     axes_.push_back(vAxis);
1376     clear();
1377   }
1378 
1379   template <typename Numeric, class Axis>
1380   HistoND<Numeric, Axis>::HistoND(const ArrayShape& shape,
1381                                   const BoxND<double>& boundingBox,
1382                                   const char* title,
1383                                   const char* label)
1384       : title_(title ? title : ""),
1385         accumulatedDataLabel_(label ? label : ""),
1386         data_(shape),
1387         overflow_(ArrayShape(shape.size(), 3U)),
1388         weightBuf_(shape.size()),
1389         indexBuf_(2U * shape.size()),
1390         modCount_(0UL),
1391         dim_(shape.size()) {
1392     if (boundingBox.size() != dim_)
1393       throw npstat::NpstatInvalidArgument(
1394           "In npstat::HistoND constructor: "
1395           "incompatible bounding box dimensionality");
1396     if (dim_ >= CHAR_BIT * sizeof(unsigned long))
1397       throw npstat::NpstatInvalidArgument(
1398           "In npstat::HistoND constructor: requested histogram "
1399           "dimensionality is not supported (too large)");
1400     axes_.reserve(dim_);
1401     for (unsigned i = 0; i < dim_; ++i)
1402       axes_.push_back(Axis(shape[i], boundingBox[i].min(), boundingBox[i].max()));
1403     clear();
1404   }
1405 
1406   template <typename Numeric, class Axis>
1407   template <typename Num2, class Functor>
1408   HistoND<Numeric, Axis>::HistoND(const HistoND<Num2, Axis>& r, const Functor& f, const char* title, const char* label)
1409       : title_(title ? title : ""),
1410         accumulatedDataLabel_(label ? label : ""),
1411         data_(r.data_, f),
1412         overflow_(r.overflow_, f),
1413         axes_(r.axes_),
1414         weightBuf_(r.dim_),
1415         indexBuf_(2U * r.dim_),
1416         fillCount_(r.fillCount_),
1417         overCount_(r.overCount_),
1418         modCount_(0UL),
1419         dim_(r.dim_) {}
1420 
1421   template <typename Numeric, class Axis>
1422   template <typename Num2>
1423   HistoND<Numeric, Axis>::HistoND(const HistoND<Num2, Axis>& h,
1424                                   const unsigned* indices,
1425                                   const unsigned nIndices,
1426                                   const char* title)
1427       : title_(title ? title : ""),
1428         accumulatedDataLabel_(h.accumulatedDataLabel_),
1429         data_(Private::shapeOfASlice(h.axes_, indices, nIndices)),
1430         overflow_(ArrayShape(data_.rank(), 3U)),
1431         axes_(Private::axesOfASlice(h.axes_, indices, nIndices)),
1432         weightBuf_(data_.rank()),
1433         indexBuf_(2U * data_.rank()),
1434         modCount_(0UL),
1435         dim_(data_.rank()) {
1436     clear();
1437   }
1438 
1439   template <typename Numeric, class Axis>
1440   template <typename Num2>
1441   HistoND<Numeric, Axis>::HistoND(const HistoND<Num2, Axis>& h,
1442                                   const Axis& newAxis,
1443                                   const unsigned newAxisNumber,
1444                                   const char* title)
1445       : title_(title ? title : ""),
1446         accumulatedDataLabel_(h.accumulatedDataLabel_),
1447         data_(Private::shapeWithExtraAxis(h.axes_, newAxis, newAxisNumber)),
1448         overflow_(data_.rank(), 3U),
1449         axes_(Private::addAxis(h.axes_, newAxis, newAxisNumber)),
1450         weightBuf_(data_.rank()),
1451         indexBuf_(2U * data_.rank()),
1452         modCount_(0UL),
1453         dim_(data_.rank()) {
1454     if (dim_ >= CHAR_BIT * sizeof(unsigned long))
1455       throw npstat::NpstatInvalidArgument(
1456           "In npstat::HistoND constructor: requested histogram "
1457           "dimensionality is not supported (too large)");
1458     clear();
1459   }
1460 
1461   template <typename Numeric, class Axis>
1462   template <typename Num2>
1463   HistoND<Numeric, Axis>::HistoND(const HistoND<Num2, Axis>& h,
1464                                   const RebinType rType,
1465                                   const unsigned* newBinCounts,
1466                                   const unsigned lenNewBinCounts,
1467                                   const double* shifts,
1468                                   const char* title)
1469       : title_(title ? title : h.title_.c_str()),
1470         accumulatedDataLabel_(h.accumulatedDataLabel_),
1471         data_(newBinCounts, lenNewBinCounts),
1472         overflow_(h.overflow_),
1473         axes_(Private::rebinAxes(h.axes_, newBinCounts, lenNewBinCounts)),
1474         weightBuf_(h.dim_),
1475         indexBuf_(2U * h.dim_),
1476         fillCount_(h.fillCount_),
1477         overCount_(h.overCount_),
1478         modCount_(0UL),
1479         dim_(h.dim_) {
1480     const unsigned long newBins = data_.length();
1481     const Axis* ax = &axes_[0];
1482     unsigned* ubuf = &indexBuf_[0];
1483 
1484     // Fill out the bins of the new histogram
1485     if (rType == SAMPLE) {
1486       double* buf = &weightBuf_[0];
1487       for (unsigned long ibin = 0; ibin < newBins; ++ibin) {
1488         data_.convertLinearIndex(ibin, ubuf, dim_);
1489         if (shifts)
1490           for (unsigned i = 0; i < dim_; ++i)
1491             buf[i] = ax[i].binCenter(ubuf[i]) + shifts[i];
1492         else
1493           for (unsigned i = 0; i < dim_; ++i)
1494             buf[i] = ax[i].binCenter(ubuf[i]);
1495         data_.linearValue(ibin) = h.examine(buf, dim_);
1496       }
1497     } else {
1498       const Numeric zero = Numeric();
1499       BoxND<double> binLimits(dim_);
1500       for (unsigned long ibin = 0; ibin < newBins; ++ibin) {
1501         data_.convertLinearIndex(ibin, ubuf, dim_);
1502         for (unsigned i = 0; i < dim_; ++i)
1503           binLimits[i] = ax[i].binInterval(ubuf[i]);
1504         Numeric& thisBin(data_.linearValue(ibin));
1505         thisBin = zero;
1506         h.accumulateBinsInBox(binLimits, &thisBin, rType == AVERAGE);
1507       }
1508     }
1509   }
1510 
1511   template <typename Numeric, class Axis>
1512   bool HistoND<Numeric, Axis>::isUniformlyBinned() const {
1513     for (unsigned i = 0; i < dim_; ++i)
1514       if (!axes_[i].isUniform())
1515         return false;
1516     return true;
1517   }
1518 
1519   template <typename Numeric, class Axis>
1520   double HistoND<Numeric, Axis>::integral() const {
1521     typedef typename PreciseType<Numeric>::type Precise;
1522 
1523     if (dim_ == 0U)
1524       return 0.0;
1525     if (isUniformlyBinned()) {
1526       Precise sum = data_.template sum<Precise>();
1527       return static_cast<double>(sum) * binVolume();
1528     } else {
1529       Precise sum = Precise();
1530       const Numeric* data = data_.data();
1531       const unsigned long len = data_.length();
1532       for (unsigned long i = 0; i < len; ++i)
1533         sum += data[i] * binVolume(i);
1534       return static_cast<double>(sum);
1535     }
1536   }
1537 
1538   template <typename Numeric, class Axis>
1539   BoxND<double> HistoND<Numeric, Axis>::boundingBox() const {
1540     BoxND<double> box;
1541     if (dim_) {
1542       box.reserve(dim_);
1543       const Axis* ax = &axes_[0];
1544       for (unsigned i = 0; i < dim_; ++i)
1545         box.push_back(ax[i].interval());
1546     }
1547     return box;
1548   }
1549 
1550   template <typename Numeric, class Axis>
1551   void HistoND<Numeric, Axis>::binCenter(const unsigned long binNumber,
1552                                          double* coords,
1553                                          const unsigned lenCoords) const {
1554     if (dim_ != lenCoords)
1555       throw npstat::NpstatInvalidArgument(
1556           "In npstat::HistoND::binCenter: "
1557           "incompatible input point dimensionality");
1558     if (dim_) {
1559       assert(coords);
1560       data_.convertLinearIndex(binNumber, &indexBuf_[0], dim_);
1561       const Axis* ax = &axes_[0];
1562       for (unsigned i = 0; i < dim_; ++i)
1563         coords[i] = ax[i].binCenter(indexBuf_[i]);
1564     }
1565   }
1566 
1567   template <typename Numeric, class Axis>
1568   template <class Point>
1569   void HistoND<Numeric, Axis>::allBinCenters(std::vector<Point>* centers) const {
1570     assert(centers);
1571     centers->clear();
1572     const unsigned long len = data_.length();
1573     centers->reserve(len);
1574     unsigned* ibuf = &indexBuf_[0];
1575     const Axis* ax = &axes_[0];
1576     Point center;
1577     if (center.size() < dim_)
1578       throw npstat::NpstatInvalidArgument(
1579           "In npstat::HistoND::allBinCenters: "
1580           "incompatible point dimensionality (too small)");
1581     typename Point::value_type* cdat = &center[0];
1582 
1583     for (unsigned long i = 0; i < len; ++i) {
1584       data_.convertLinearIndex(i, ibuf, dim_);
1585       for (unsigned idim = 0; idim < dim_; ++idim)
1586         cdat[idim] = ax[idim].binCenter(ibuf[idim]);
1587       centers->push_back(center);
1588     }
1589   }
1590 
1591   template <typename Numeric, class Axis>
1592   void HistoND<Numeric, Axis>::binBox(const unsigned long binNumber, BoxND<double>* box) const {
1593     assert(box);
1594     box->clear();
1595     if (dim_) {
1596       box->reserve(dim_);
1597       data_.convertLinearIndex(binNumber, &indexBuf_[0], dim_);
1598       const Axis* ax = &axes_[0];
1599       for (unsigned i = 0; i < dim_; ++i)
1600         box->push_back(ax[i].binInterval(indexBuf_[i]));
1601     }
1602   }
1603 
1604   template <typename Numeric, class Axis>
1605   inline bool HistoND<Numeric, Axis>::isSameData(const HistoND& r) const {
1606     return dim_ == r.dim_ && overflow_ == r.overflow_ && data_ == r.data_;
1607   }
1608 
1609   template <typename Numeric, class Axis>
1610   inline bool HistoND<Numeric, Axis>::operator==(const HistoND& r) const {
1611     return dim_ == r.dim_ && fillCount_ == r.fillCount_ && overCount_ == r.overCount_ && title_ == r.title_ &&
1612            accumulatedDataLabel_ == r.accumulatedDataLabel_ && axes_ == r.axes_ && overflow_ == r.overflow_ &&
1613            data_ == r.data_;
1614   }
1615 
1616   template <typename Numeric, class Axis>
1617   inline bool HistoND<Numeric, Axis>::operator!=(const HistoND& r) const {
1618     return !(*this == r);
1619   }
1620 
1621   template <typename Numeric, class Axis>
1622   double HistoND<Numeric, Axis>::binVolume(const unsigned long binNumber) const {
1623     double v = 1.0;
1624     if (dim_) {
1625       data_.convertLinearIndex(binNumber, &indexBuf_[0], dim_);
1626       const Axis* ax = &axes_[0];
1627       for (unsigned i = 0; i < dim_; ++i)
1628         v *= ax[i].binWidth(indexBuf_[i]);
1629     }
1630     return v;
1631   }
1632 
1633   template <typename Numeric, class Axis>
1634   double HistoND<Numeric, Axis>::volume() const {
1635     double v = 1.0;
1636     if (dim_) {
1637       const Axis* ax = &axes_[0];
1638       for (unsigned i = 0; i < dim_; ++i)
1639         v *= (ax[i].max() - ax[i].min());
1640     }
1641     return v;
1642   }
1643 
1644   template <typename Numeric, class Axis>
1645   template <typename Num2>
1646   void HistoND<Numeric, Axis>::fill(const double* coords, const unsigned coordLength, const Num2& w) {
1647     if (coordLength != dim_)
1648       throw npstat::NpstatInvalidArgument(
1649           "In npstat::HistoND::fill: "
1650           "incompatible input point dimensionality");
1651     if (coordLength) {
1652       assert(coords);
1653       unsigned* idx = &indexBuf_[0];
1654       unsigned* over = idx + dim_;
1655       const Axis* ax = &axes_[0];
1656       unsigned overflown = 0U;
1657       for (unsigned i = 0; i < dim_; ++i) {
1658         over[i] = ax[i].overflowIndex(coords[i], idx + i);
1659         overflown |= (over[i] - 1U);
1660       }
1661       if (overflown) {
1662         overflow_.value(over, dim_) += w;
1663         ++overCount_;
1664       } else
1665         data_.value(idx, dim_) += w;
1666     } else
1667       data_() += w;
1668     ++fillCount_;
1669     ++modCount_;
1670   }
1671 
1672   template <typename Numeric, class Axis>
1673   template <typename Num2, class Functor>
1674   void HistoND<Numeric, Axis>::dispatch(const double* coords, const unsigned coordLength, Num2& w, Functor& f) {
1675     if (coordLength != dim_)
1676       throw npstat::NpstatInvalidArgument(
1677           "In npstat::HistoND::dispatch: "
1678           "incompatible input point dimensionality");
1679     if (coordLength) {
1680       assert(coords);
1681       unsigned* idx = &indexBuf_[0];
1682       unsigned* over = idx + dim_;
1683       const Axis* ax = &axes_[0];
1684       unsigned overflown = 0U;
1685       for (unsigned i = 0; i < dim_; ++i) {
1686         over[i] = ax[i].overflowIndex(coords[i], idx + i);
1687         overflown |= (over[i] - 1U);
1688       }
1689       if (overflown)
1690         f(overflow_.value(over, dim_), w);
1691       else
1692         f(data_.value(idx, dim_), w);
1693     } else
1694       f(data_(), w);
1695     ++modCount_;
1696   }
1697 
1698   template <typename Numeric, class Axis>
1699   const Numeric& HistoND<Numeric, Axis>::examine(const double* coords, const unsigned coordLength) const {
1700     if (coordLength != dim_)
1701       throw npstat::NpstatInvalidArgument(
1702           "In npstat::HistoND::examine: "
1703           "incompatible input point dimensionality");
1704     if (coordLength) {
1705       assert(coords);
1706       unsigned* idx = &indexBuf_[0];
1707       unsigned* over = idx + dim_;
1708       const Axis* ax = &axes_[0];
1709       unsigned overflown = 0U;
1710       for (unsigned i = 0; i < dim_; ++i) {
1711         over[i] = ax[i].overflowIndex(coords[i], idx + i);
1712         overflown |= (over[i] - 1U);
1713       }
1714       if (overflown)
1715         return overflow_.value(over, dim_);
1716       else
1717         return data_.value(idx, dim_);
1718     } else
1719       return data_();
1720   }
1721 
1722   template <typename Numeric, class Axis>
1723   const Numeric& HistoND<Numeric, Axis>::closestBin(const double* coords, const unsigned coordLength) const {
1724     if (coordLength != dim_)
1725       throw npstat::NpstatInvalidArgument(
1726           "In npstat::HistoND::closestBin: "
1727           "incompatible input point dimensionality");
1728     if (coordLength) {
1729       assert(coords);
1730       unsigned* idx = &indexBuf_[0];
1731       const Axis* ax = &axes_[0];
1732       for (unsigned i = 0; i < dim_; ++i)
1733         idx[i] = ax[i].closestValidBin(coords[i]);
1734       return data_.value(idx, dim_);
1735     } else
1736       return data_();
1737   }
1738 
1739   template <typename Numeric, class Axis>
1740   template <typename Num2>
1741   void HistoND<Numeric, Axis>::fillPreservingCentroid(const Num2& value) {
1742     const double* weights = &weightBuf_[0];
1743     const unsigned* cell = &indexBuf_[0];
1744     const unsigned long* strides = data_.strides();
1745     const unsigned long maxcycle = 1UL << dim_;
1746     for (unsigned long icycle = 0; icycle < maxcycle; ++icycle) {
1747       double w = 1.0;
1748       unsigned long icell = 0UL;
1749       for (unsigned i = 0; i < dim_; ++i) {
1750         if (icycle & (1UL << i)) {
1751           w *= (1.0 - weights[i]);
1752           icell += strides[i] * (cell[i] + 1U);
1753         } else {
1754           w *= weights[i];
1755           icell += strides[i] * cell[i];
1756         }
1757       }
1758       data_.linearValue(icell) += (value * w);
1759     }
1760   }
1761 
1762   template <typename Numeric, class Axis>
1763   template <typename Num2>
1764   void HistoND<Numeric, Axis>::fillC(const double* coords, const unsigned coordLength, const Num2& w) {
1765     if (coordLength != dim_)
1766       throw npstat::NpstatInvalidArgument(
1767           "In npstat::HistoND::fillC: "
1768           "incompatible input point dimensionality");
1769     if (coordLength) {
1770       assert(coords);
1771       double* wg = &weightBuf_[0];
1772       unsigned* idx = &indexBuf_[0];
1773       unsigned* over = idx + dim_;
1774       const Axis* ax = &axes_[0];
1775       unsigned overflown = 0U;
1776       for (unsigned i = 0; i < dim_; ++i) {
1777         over[i] = ax[i].overflowIndexWeighted(coords[i], idx + i, wg + i);
1778         overflown |= (over[i] - 1U);
1779       }
1780       if (overflown) {
1781         overflow_.value(over, dim_) += w;
1782         ++overCount_;
1783       } else
1784         fillPreservingCentroid(w);
1785     } else
1786       data_() += w;
1787     ++fillCount_;
1788     ++modCount_;
1789   }
1790 
1791   template <typename Numeric, class Axis>
1792   template <typename Num2>
1793   inline void HistoND<Numeric, Axis>::fill(const Num2& w) {
1794     if (dim_)
1795       Private::h_badargs("fill");
1796     data_() += w;
1797     ++fillCount_;
1798     ++modCount_;
1799   }
1800 
1801   template <typename Numeric, class Axis>
1802   template <typename Num2, class Functor>
1803   inline void HistoND<Numeric, Axis>::dispatch(Num2& w, Functor& f) {
1804     if (dim_)
1805       Private::h_badargs("dispatch");
1806     f(data_(), w);
1807     ++modCount_;
1808   }
1809 
1810   template <typename Numeric, class Axis>
1811   template <typename Num2>
1812   inline void HistoND<Numeric, Axis>::fillC(const Num2& w) {
1813     if (dim_)
1814       Private::h_badargs("fillC");
1815     data_() += w;
1816     ++fillCount_;
1817     ++modCount_;
1818   }
1819 
1820   template <typename Numeric, class Axis>
1821   inline const Numeric& HistoND<Numeric, Axis>::examine() const {
1822     if (dim_)
1823       Private::h_badargs("examine");
1824     return data_();
1825   }
1826 
1827   template <typename Numeric, class Axis>
1828   inline const Numeric& HistoND<Numeric, Axis>::closestBin() const {
1829     if (dim_)
1830       Private::h_badargs("closestBin");
1831     return data_();
1832   }
1833 
1834   template <typename Numeric, class Axis>
1835   template <typename Num2>
1836   void HistoND<Numeric, Axis>::fill(const double x0, const Num2& w) {
1837     if (dim_ != 1U)
1838       Private::h_badargs("fill");
1839     unsigned i0 = 0;
1840     const unsigned ov0 = axes_[0].overflowIndex(x0, &i0);
1841     if (ov0 == 1U)
1842       data_(i0) += w;
1843     else {
1844       overflow_(ov0) += w;
1845       ++overCount_;
1846     }
1847     ++fillCount_;
1848     ++modCount_;
1849   }
1850 
1851   template <typename Numeric, class Axis>
1852   template <typename Num2, class Functor>
1853   void HistoND<Numeric, Axis>::dispatch(const double x0, Num2& w, Functor& f) {
1854     if (dim_ != 1U)
1855       Private::h_badargs("dispatch");
1856     unsigned i0 = 0;
1857     const unsigned ov0 = axes_[0].overflowIndex(x0, &i0);
1858     if (ov0 == 1U)
1859       f(data_(i0), w);
1860     else
1861       f(overflow_(ov0), w);
1862     ++modCount_;
1863   }
1864 
1865   template <typename Numeric, class Axis>
1866   template <typename Num2>
1867   void HistoND<Numeric, Axis>::fillC(const double x0, const Num2& w) {
1868     if (dim_ != 1U)
1869       Private::h_badargs("fillC");
1870     double* wg = &weightBuf_[0];
1871     unsigned* idx = &indexBuf_[0];
1872     const unsigned ov0 = axes_[0].overflowIndexWeighted(x0, idx, wg);
1873     if (ov0 == 1U)
1874       fillPreservingCentroid(w);
1875     else {
1876       overflow_(ov0) += w;
1877       ++overCount_;
1878     }
1879     ++fillCount_;
1880     ++modCount_;
1881   }
1882 
1883   template <typename Numeric, class Axis>
1884   inline const Numeric& HistoND<Numeric, Axis>::examine(const double x0) const {
1885     if (dim_ != 1U)
1886       Private::h_badargs("examine");
1887     unsigned i0 = 0;
1888     const unsigned ov0 = axes_[0].overflowIndex(x0, &i0);
1889     if (ov0 == 1U)
1890       return data_(i0);
1891     else
1892       return overflow_(ov0);
1893   }
1894 
1895   template <typename Numeric, class Axis>
1896   inline const Numeric& HistoND<Numeric, Axis>::closestBin(const double x0) const {
1897     if (dim_ != 1U)
1898       Private::h_badargs("closestBin");
1899     const unsigned i0 = axes_[0].closestValidBin(x0);
1900     return data_(i0);
1901   }
1902 
1903   template <typename Numeric, class Axis>
1904   template <typename Num2>
1905   void HistoND<Numeric, Axis>::fill(const double x0, const double x1, const Num2& w) {
1906     if (dim_ != 2U)
1907       Private::h_badargs("fill");
1908     unsigned i0 = 0, i1 = 0;
1909     const Axis* ax = &axes_[0];
1910     const unsigned o0 = ax[0].overflowIndex(x0, &i0);
1911     const unsigned o1 = ax[1].overflowIndex(x1, &i1);
1912     if (o0 == 1U && o1 == 1U)
1913       data_(i0, i1) += w;
1914     else {
1915       overflow_(o0, o1) += w;
1916       ++overCount_;
1917     }
1918     ++fillCount_;
1919     ++modCount_;
1920   }
1921 
1922   template <typename Numeric, class Axis>
1923   template <typename Num2, class Functor>
1924   void HistoND<Numeric, Axis>::dispatch(const double x0, const double x1, Num2& w, Functor& f) {
1925     if (dim_ != 2U)
1926       Private::h_badargs("dispatch");
1927     unsigned i0 = 0, i1 = 0;
1928     const Axis* ax = &axes_[0];
1929     const unsigned o0 = ax[0].overflowIndex(x0, &i0);
1930     const unsigned o1 = ax[1].overflowIndex(x1, &i1);
1931     if (o0 == 1U && o1 == 1U)
1932       f(data_(i0, i1), w);
1933     else
1934       f(overflow_(o0, o1), 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, const Num2& w) {
1941     if (dim_ != 2U)
1942       Private::h_badargs("fillC");
1943     double* wg = &weightBuf_[0];
1944     unsigned* idx = &indexBuf_[0];
1945     const Axis* ax = &axes_[0];
1946     const unsigned o0 = ax[0].overflowIndexWeighted(x0, idx + 0, wg + 0);
1947     const unsigned o1 = ax[1].overflowIndexWeighted(x1, idx + 1, wg + 1);
1948     if (o0 == 1U && o1 == 1U)
1949       fillPreservingCentroid(w);
1950     else {
1951       overflow_(o0, o1) += w;
1952       ++overCount_;
1953     }
1954     ++fillCount_;
1955     ++modCount_;
1956   }
1957 
1958   template <typename Numeric, class Axis>
1959   const Numeric& HistoND<Numeric, Axis>::examine(const double x0, const double x1) const {
1960     if (dim_ != 2U)
1961       Private::h_badargs("examine");
1962     unsigned i0 = 0, i1 = 0;
1963     const Axis* ax = &axes_[0];
1964     const unsigned o0 = ax[0].overflowIndex(x0, &i0);
1965     const unsigned o1 = ax[1].overflowIndex(x1, &i1);
1966     if (o0 == 1U && o1 == 1U)
1967       return data_(i0, i1);
1968     else
1969       return overflow_(o0, o1);
1970   }
1971 
1972   template <typename Numeric, class Axis>
1973   const Numeric& HistoND<Numeric, Axis>::closestBin(const double x0, const double x1) const {
1974     if (dim_ != 2U)
1975       Private::h_badargs("closestBin");
1976     const Axis* ax = &axes_[0];
1977     const unsigned i0 = ax[0].closestValidBin(x0);
1978     const unsigned i1 = ax[1].closestValidBin(x1);
1979     return data_(i0, i1);
1980   }
1981 
1982   template <typename Numeric, class Axis>
1983   template <typename Num2>
1984   void HistoND<Numeric, Axis>::fill(const double x0, const double x1, const double x2, const Num2& w) {
1985     if (dim_ != 3U)
1986       Private::h_badargs("fill");
1987     unsigned i0 = 0, i1 = 0, i2 = 0;
1988     const Axis* ax = &axes_[0];
1989     const unsigned o0 = ax[0].overflowIndex(x0, &i0);
1990     const unsigned o1 = ax[1].overflowIndex(x1, &i1);
1991     const unsigned o2 = ax[2].overflowIndex(x2, &i2);
1992     if (o0 == 1U && o1 == 1U && o2 == 1U)
1993       data_(i0, i1, i2) += w;
1994     else {
1995       overflow_(o0, o1, o2) += w;
1996       ++overCount_;
1997     }
1998     ++fillCount_;
1999     ++modCount_;
2000   }
2001 
2002   template <typename Numeric, class Axis>
2003   template <typename Num2, class Functor>
2004   void HistoND<Numeric, Axis>::dispatch(const double x0, const double x1, const double x2, Num2& w, Functor& f) {
2005     if (dim_ != 3U)
2006       Private::h_badargs("dispatch");
2007     unsigned i0 = 0, i1 = 0, i2 = 0;
2008     const Axis* ax = &axes_[0];
2009     const unsigned o0 = ax[0].overflowIndex(x0, &i0);
2010     const unsigned o1 = ax[1].overflowIndex(x1, &i1);
2011     const unsigned o2 = ax[2].overflowIndex(x2, &i2);
2012     if (o0 == 1U && o1 == 1U && o2 == 1U)
2013       f(data_(i0, i1, i2), w);
2014     else
2015       f(overflow_(o0, o1, o2), w);
2016     ++modCount_;
2017   }
2018 
2019   template <typename Numeric, class Axis>
2020   template <typename Num2>
2021   void HistoND<Numeric, Axis>::fillC(const double x0, const double x1, const double x2, const Num2& w) {
2022     if (dim_ != 3U)
2023       Private::h_badargs("fillC");
2024     double* wg = &weightBuf_[0];
2025     unsigned* idx = &indexBuf_[0];
2026     const Axis* ax = &axes_[0];
2027     const unsigned o0 = ax[0].overflowIndexWeighted(x0, idx + 0, wg + 0);
2028     const unsigned o1 = ax[1].overflowIndexWeighted(x1, idx + 1, wg + 1);
2029     const unsigned o2 = ax[2].overflowIndexWeighted(x2, idx + 2, wg + 2);
2030     if (o0 == 1U && o1 == 1U && o2 == 1U)
2031       fillPreservingCentroid(w);
2032     else {
2033       overflow_(o0, o1, o2) += w;
2034       ++overCount_;
2035     }
2036     ++fillCount_;
2037     ++modCount_;
2038   }
2039 
2040   template <typename Numeric, class Axis>
2041   const Numeric& HistoND<Numeric, Axis>::examine(const double x0, const double x1, const double x2) const {
2042     if (dim_ != 3U)
2043       Private::h_badargs("examine");
2044     unsigned i0 = 0, i1 = 0, i2 = 0;
2045     const Axis* ax = &axes_[0];
2046     const unsigned o0 = ax[0].overflowIndex(x0, &i0);
2047     const unsigned o1 = ax[1].overflowIndex(x1, &i1);
2048     const unsigned o2 = ax[2].overflowIndex(x2, &i2);
2049     if (o0 == 1U && o1 == 1U && o2 == 1U)
2050       return data_(i0, i1, i2);
2051     else
2052       return overflow_(o0, o1, o2);
2053   }
2054 
2055   template <typename Numeric, class Axis>
2056   const Numeric& HistoND<Numeric, Axis>::closestBin(const double x0, const double x1, const double x2) const {
2057     if (dim_ != 3U)
2058       Private::h_badargs("closestBin");
2059     const Axis* ax = &axes_[0];
2060     const unsigned i0 = ax[0].closestValidBin(x0);
2061     const unsigned i1 = ax[1].closestValidBin(x1);
2062     const unsigned i2 = ax[2].closestValidBin(x2);
2063     return data_(i0, i1, i2);
2064   }
2065 
2066   template <typename Numeric, class Axis>
2067   template <typename Num2>
2068   void HistoND<Numeric, Axis>::fill(const double x0, const double x1, const double x2, const double x3, const Num2& w) {
2069     if (dim_ != 4U)
2070       Private::h_badargs("fill");
2071     unsigned i0 = 0, i1 = 0, i2 = 0, i3 = 0;
2072     const Axis* ax = &axes_[0];
2073     const unsigned o0 = ax[0].overflowIndex(x0, &i0);
2074     const unsigned o1 = ax[1].overflowIndex(x1, &i1);
2075     const unsigned o2 = ax[2].overflowIndex(x2, &i2);
2076     const unsigned o3 = ax[3].overflowIndex(x3, &i3);
2077     if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U)
2078       data_(i0, i1, i2, i3) += w;
2079     else {
2080       overflow_(o0, o1, o2, o3) += w;
2081       ++overCount_;
2082     }
2083     ++fillCount_;
2084     ++modCount_;
2085   }
2086 
2087   template <typename Numeric, class Axis>
2088   template <typename Num2, class Functor>
2089   void HistoND<Numeric, Axis>::dispatch(
2090       const double x0, const double x1, const double x2, const double x3, Num2& w, Functor& f) {
2091     if (dim_ != 4U)
2092       Private::h_badargs("dispatch");
2093     unsigned i0 = 0, i1 = 0, i2 = 0, i3 = 0;
2094     const Axis* ax = &axes_[0];
2095     const unsigned o0 = ax[0].overflowIndex(x0, &i0);
2096     const unsigned o1 = ax[1].overflowIndex(x1, &i1);
2097     const unsigned o2 = ax[2].overflowIndex(x2, &i2);
2098     const unsigned o3 = ax[3].overflowIndex(x3, &i3);
2099     if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U)
2100       f(data_(i0, i1, i2, i3), w);
2101     else
2102       f(overflow_(o0, o1, o2, o3), w);
2103     ++modCount_;
2104   }
2105 
2106   template <typename Numeric, class Axis>
2107   template <typename Num2>
2108   void HistoND<Numeric, Axis>::fillC(const double x0, const double x1, const double x2, const double x3, const Num2& w) {
2109     if (dim_ != 4U)
2110       Private::h_badargs("fillC");
2111     double* wg = &weightBuf_[0];
2112     unsigned* idx = &indexBuf_[0];
2113     const Axis* ax = &axes_[0];
2114     const unsigned o0 = ax[0].overflowIndexWeighted(x0, idx + 0, wg + 0);
2115     const unsigned o1 = ax[1].overflowIndexWeighted(x1, idx + 1, wg + 1);
2116     const unsigned o2 = ax[2].overflowIndexWeighted(x2, idx + 2, wg + 2);
2117     const unsigned o3 = ax[3].overflowIndexWeighted(x3, idx + 3, wg + 3);
2118     if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U)
2119       fillPreservingCentroid(w);
2120     else {
2121       overflow_(o0, o1, o2, o3) += w;
2122       ++overCount_;
2123     }
2124     ++fillCount_;
2125     ++modCount_;
2126   }
2127 
2128   template <typename Numeric, class Axis>
2129   const Numeric& HistoND<Numeric, Axis>::examine(const double x0,
2130                                                  const double x1,
2131                                                  const double x2,
2132                                                  const double x3) const {
2133     if (dim_ != 4U)
2134       Private::h_badargs("examine");
2135     unsigned i0 = 0, i1 = 0, i2 = 0, i3 = 0;
2136     const Axis* ax = &axes_[0];
2137     const unsigned o0 = ax[0].overflowIndex(x0, &i0);
2138     const unsigned o1 = ax[1].overflowIndex(x1, &i1);
2139     const unsigned o2 = ax[2].overflowIndex(x2, &i2);
2140     const unsigned o3 = ax[3].overflowIndex(x3, &i3);
2141     if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U)
2142       return data_(i0, i1, i2, i3);
2143     else
2144       return overflow_(o0, o1, o2, o3);
2145   }
2146 
2147   template <typename Numeric, class Axis>
2148   const Numeric& HistoND<Numeric, Axis>::closestBin(const double x0,
2149                                                     const double x1,
2150                                                     const double x2,
2151                                                     const double x3) const {
2152     if (dim_ != 4U)
2153       Private::h_badargs("closestBin");
2154     const Axis* ax = &axes_[0];
2155     const unsigned i0 = ax[0].closestValidBin(x0);
2156     const unsigned i1 = ax[1].closestValidBin(x1);
2157     const unsigned i2 = ax[2].closestValidBin(x2);
2158     const unsigned i3 = ax[3].closestValidBin(x3);
2159     return data_(i0, i1, i2, i3);
2160   }
2161 
2162   template <typename Numeric, class Axis>
2163   template <typename Num2>
2164   void HistoND<Numeric, Axis>::fill(
2165       const double x0, const double x1, const double x2, const double x3, const double x4, const Num2& w) {
2166     if (dim_ != 5U)
2167       Private::h_badargs("fill");
2168     unsigned i0 = 0, i1 = 0, i2 = 0, i3 = 0, i4 = 0;
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       data_(i0, i1, i2, i3, i4) += w;
2177     else {
2178       overflow_(o0, o1, o2, o3, o4) += w;
2179       ++overCount_;
2180     }
2181     ++fillCount_;
2182     ++modCount_;
2183   }
2184 
2185   template <typename Numeric, class Axis>
2186   template <typename Num2, class Functor>
2187   void HistoND<Numeric, Axis>::dispatch(
2188       const double x0, const double x1, const double x2, const double x3, const double x4, Num2& w, Functor& f) {
2189     if (dim_ != 5U)
2190       Private::h_badargs("dispatch");
2191     unsigned i0 = 0, i1 = 0, i2 = 0, i3 = 0, i4 = 0;
2192     const Axis* ax = &axes_[0];
2193     const unsigned o0 = ax[0].overflowIndex(x0, &i0);
2194     const unsigned o1 = ax[1].overflowIndex(x1, &i1);
2195     const unsigned o2 = ax[2].overflowIndex(x2, &i2);
2196     const unsigned o3 = ax[3].overflowIndex(x3, &i3);
2197     const unsigned o4 = ax[4].overflowIndex(x4, &i4);
2198     if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U && o4 == 1U)
2199       f(data_(i0, i1, i2, i3, i4), w);
2200     else
2201       f(overflow_(o0, o1, o2, o3, o4), w);
2202     ++modCount_;
2203   }
2204 
2205   template <typename Numeric, class Axis>
2206   template <typename Num2>
2207   void HistoND<Numeric, Axis>::fillC(
2208       const double x0, const double x1, const double x2, const double x3, const double x4, const Num2& w) {
2209     if (dim_ != 5U)
2210       Private::h_badargs("fillC");
2211     double* wg = &weightBuf_[0];
2212     unsigned* idx = &indexBuf_[0];
2213     const Axis* ax = &axes_[0];
2214     const unsigned o0 = ax[0].overflowIndexWeighted(x0, idx + 0, wg + 0);
2215     const unsigned o1 = ax[1].overflowIndexWeighted(x1, idx + 1, wg + 1);
2216     const unsigned o2 = ax[2].overflowIndexWeighted(x2, idx + 2, wg + 2);
2217     const unsigned o3 = ax[3].overflowIndexWeighted(x3, idx + 3, wg + 3);
2218     const unsigned o4 = ax[4].overflowIndexWeighted(x4, idx + 4, wg + 4);
2219     if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U && o4 == 1U)
2220       fillPreservingCentroid(w);
2221     else {
2222       overflow_(o0, o1, o2, o3, o4) += w;
2223       ++overCount_;
2224     }
2225     ++fillCount_;
2226     ++modCount_;
2227   }
2228 
2229   template <typename Numeric, class Axis>
2230   const Numeric& HistoND<Numeric, Axis>::examine(
2231       const double x0, const double x1, const double x2, const double x3, const double x4) const {
2232     if (dim_ != 5U)
2233       Private::h_badargs("examine");
2234     unsigned i0 = 0, i1 = 0, i2 = 0, i3 = 0, i4 = 0;
2235     const Axis* ax = &axes_[0];
2236     const unsigned o0 = ax[0].overflowIndex(x0, &i0);
2237     const unsigned o1 = ax[1].overflowIndex(x1, &i1);
2238     const unsigned o2 = ax[2].overflowIndex(x2, &i2);
2239     const unsigned o3 = ax[3].overflowIndex(x3, &i3);
2240     const unsigned o4 = ax[4].overflowIndex(x4, &i4);
2241     if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U && o4 == 1U)
2242       return data_(i0, i1, i2, i3, i4);
2243     else
2244       return overflow_(o0, o1, o2, o3, o4);
2245   }
2246 
2247   template <typename Numeric, class Axis>
2248   const Numeric& HistoND<Numeric, Axis>::closestBin(
2249       const double x0, const double x1, const double x2, const double x3, const double x4) const {
2250     if (dim_ != 5U)
2251       Private::h_badargs("closestBin");
2252     const Axis* ax = &axes_[0];
2253     const unsigned i0 = ax[0].closestValidBin(x0);
2254     const unsigned i1 = ax[1].closestValidBin(x1);
2255     const unsigned i2 = ax[2].closestValidBin(x2);
2256     const unsigned i3 = ax[3].closestValidBin(x3);
2257     const unsigned i4 = ax[4].closestValidBin(x4);
2258     return data_(i0, i1, i2, i3, i4);
2259   }
2260 
2261   template <typename Numeric, class Axis>
2262   template <typename Num2>
2263   void HistoND<Numeric, Axis>::fill(const double x0,
2264                                     const double x1,
2265                                     const double x2,
2266                                     const double x3,
2267                                     const double x4,
2268                                     const double x5,
2269                                     const Num2& w) {
2270     if (dim_ != 6U)
2271       Private::h_badargs("fill");
2272     unsigned i0 = 0, i1 = 0, i2 = 0, i3 = 0, i4 = 0, i5 = 0;
2273     const Axis* ax = &axes_[0];
2274     const unsigned o0 = ax[0].overflowIndex(x0, &i0);
2275     const unsigned o1 = ax[1].overflowIndex(x1, &i1);
2276     const unsigned o2 = ax[2].overflowIndex(x2, &i2);
2277     const unsigned o3 = ax[3].overflowIndex(x3, &i3);
2278     const unsigned o4 = ax[4].overflowIndex(x4, &i4);
2279     const unsigned o5 = ax[5].overflowIndex(x5, &i5);
2280     if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U && o4 == 1U && o5 == 1U)
2281       data_(i0, i1, i2, i3, i4, i5) += w;
2282     else {
2283       overflow_(o0, o1, o2, o3, o4, o5) += w;
2284       ++overCount_;
2285     }
2286     ++fillCount_;
2287     ++modCount_;
2288   }
2289 
2290   template <typename Numeric, class Axis>
2291   template <typename Num2, class Functor>
2292   void HistoND<Numeric, Axis>::dispatch(const double x0,
2293                                         const double x1,
2294                                         const double x2,
2295                                         const double x3,
2296                                         const double x4,
2297                                         const double x5,
2298                                         Num2& w,
2299                                         Functor& f) {
2300     if (dim_ != 6U)
2301       Private::h_badargs("dispatch");
2302     unsigned i0 = 0, i1 = 0, i2 = 0, i3 = 0, i4 = 0, i5 = 0;
2303     const Axis* ax = &axes_[0];
2304     const unsigned o0 = ax[0].overflowIndex(x0, &i0);
2305     const unsigned o1 = ax[1].overflowIndex(x1, &i1);
2306     const unsigned o2 = ax[2].overflowIndex(x2, &i2);
2307     const unsigned o3 = ax[3].overflowIndex(x3, &i3);
2308     const unsigned o4 = ax[4].overflowIndex(x4, &i4);
2309     const unsigned o5 = ax[5].overflowIndex(x5, &i5);
2310     if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U && o4 == 1U && o5 == 1U)
2311       f(data_(i0, i1, i2, i3, i4, i5), w);
2312     else
2313       f(overflow_(o0, o1, o2, o3, o4, o5), w);
2314     ++modCount_;
2315   }
2316 
2317   template <typename Numeric, class Axis>
2318   template <typename Num2>
2319   void HistoND<Numeric, Axis>::fillC(const double x0,
2320                                      const double x1,
2321                                      const double x2,
2322                                      const double x3,
2323                                      const double x4,
2324                                      const double x5,
2325                                      const Num2& w) {
2326     if (dim_ != 6U)
2327       Private::h_badargs("fillC");
2328     double* wg = &weightBuf_[0];
2329     unsigned* idx = &indexBuf_[0];
2330     const Axis* ax = &axes_[0];
2331     const unsigned o0 = ax[0].overflowIndexWeighted(x0, idx + 0, wg + 0);
2332     const unsigned o1 = ax[1].overflowIndexWeighted(x1, idx + 1, wg + 1);
2333     const unsigned o2 = ax[2].overflowIndexWeighted(x2, idx + 2, wg + 2);
2334     const unsigned o3 = ax[3].overflowIndexWeighted(x3, idx + 3, wg + 3);
2335     const unsigned o4 = ax[4].overflowIndexWeighted(x4, idx + 4, wg + 4);
2336     const unsigned o5 = ax[5].overflowIndexWeighted(x5, idx + 5, wg + 5);
2337     if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U && o4 == 1U && o5 == 1U)
2338       fillPreservingCentroid(w);
2339     else {
2340       overflow_(o0, o1, o2, o3, o4, o5) += w;
2341       ++overCount_;
2342     }
2343     ++fillCount_;
2344     ++modCount_;
2345   }
2346 
2347   template <typename Numeric, class Axis>
2348   const Numeric& HistoND<Numeric, Axis>::examine(
2349       const double x0, const double x1, const double x2, const double x3, const double x4, const double x5) const {
2350     if (dim_ != 6U)
2351       Private::h_badargs("examine");
2352     unsigned i0 = 0, i1 = 0, i2 = 0, i3 = 0, i4 = 0, i5 = 0;
2353     const Axis* ax = &axes_[0];
2354     const unsigned o0 = ax[0].overflowIndex(x0, &i0);
2355     const unsigned o1 = ax[1].overflowIndex(x1, &i1);
2356     const unsigned o2 = ax[2].overflowIndex(x2, &i2);
2357     const unsigned o3 = ax[3].overflowIndex(x3, &i3);
2358     const unsigned o4 = ax[4].overflowIndex(x4, &i4);
2359     const unsigned o5 = ax[5].overflowIndex(x5, &i5);
2360     if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U && o4 == 1U && o5 == 1U)
2361       return data_(i0, i1, i2, i3, i4, i5);
2362     else
2363       return overflow_(o0, o1, o2, o3, o4, o5);
2364   }
2365 
2366   template <typename Numeric, class Axis>
2367   const Numeric& HistoND<Numeric, Axis>::closestBin(
2368       const double x0, const double x1, const double x2, const double x3, const double x4, const double x5) const {
2369     if (dim_ != 6U)
2370       Private::h_badargs("closestBin");
2371     const Axis* ax = &axes_[0];
2372     const unsigned i0 = ax[0].closestValidBin(x0);
2373     const unsigned i1 = ax[1].closestValidBin(x1);
2374     const unsigned i2 = ax[2].closestValidBin(x2);
2375     const unsigned i3 = ax[3].closestValidBin(x3);
2376     const unsigned i4 = ax[4].closestValidBin(x4);
2377     const unsigned i5 = ax[5].closestValidBin(x5);
2378     return data_(i0, i1, i2, i3, i4, i5);
2379   }
2380 
2381   template <typename Numeric, class Axis>
2382   template <typename Num2>
2383   void HistoND<Numeric, Axis>::fill(const double x0,
2384                                     const double x1,
2385                                     const double x2,
2386                                     const double x3,
2387                                     const double x4,
2388                                     const double x5,
2389                                     const double x6,
2390                                     const Num2& w) {
2391     if (dim_ != 7U)
2392       Private::h_badargs("fill");
2393     unsigned i0 = 0, i1 = 0, i2 = 0, i3 = 0, i4 = 0, i5 = 0, i6 = 0;
2394     const Axis* ax = &axes_[0];
2395     const unsigned o0 = ax[0].overflowIndex(x0, &i0);
2396     const unsigned o1 = ax[1].overflowIndex(x1, &i1);
2397     const unsigned o2 = ax[2].overflowIndex(x2, &i2);
2398     const unsigned o3 = ax[3].overflowIndex(x3, &i3);
2399     const unsigned o4 = ax[4].overflowIndex(x4, &i4);
2400     const unsigned o5 = ax[5].overflowIndex(x5, &i5);
2401     const unsigned o6 = ax[6].overflowIndex(x6, &i6);
2402     if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U && o4 == 1U && o5 == 1U && o6 == 1U)
2403       data_(i0, i1, i2, i3, i4, i5, i6) += w;
2404     else {
2405       overflow_(o0, o1, o2, o3, o4, o5, o6) += w;
2406       ++overCount_;
2407     }
2408     ++fillCount_;
2409     ++modCount_;
2410   }
2411 
2412   template <typename Numeric, class Axis>
2413   template <typename Num2, class Functor>
2414   void HistoND<Numeric, Axis>::dispatch(const double x0,
2415                                         const double x1,
2416                                         const double x2,
2417                                         const double x3,
2418                                         const double x4,
2419                                         const double x5,
2420                                         const double x6,
2421                                         Num2& w,
2422                                         Functor& f) {
2423     if (dim_ != 7U)
2424       Private::h_badargs("dispatch");
2425     unsigned i0 = 0, i1 = 0, i2 = 0, i3 = 0, i4 = 0, i5 = 0, i6 = 0;
2426     const Axis* ax = &axes_[0];
2427     const unsigned o0 = ax[0].overflowIndex(x0, &i0);
2428     const unsigned o1 = ax[1].overflowIndex(x1, &i1);
2429     const unsigned o2 = ax[2].overflowIndex(x2, &i2);
2430     const unsigned o3 = ax[3].overflowIndex(x3, &i3);
2431     const unsigned o4 = ax[4].overflowIndex(x4, &i4);
2432     const unsigned o5 = ax[5].overflowIndex(x5, &i5);
2433     const unsigned o6 = ax[6].overflowIndex(x6, &i6);
2434     if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U && o4 == 1U && o5 == 1U && o6 == 1U)
2435       f(data_(i0, i1, i2, i3, i4, i5, i6), w);
2436     else
2437       f(overflow_(o0, o1, o2, o3, o4, o5, o6), w);
2438     ++modCount_;
2439   }
2440 
2441   template <typename Numeric, class Axis>
2442   template <typename Num2>
2443   void HistoND<Numeric, Axis>::fillC(const double x0,
2444                                      const double x1,
2445                                      const double x2,
2446                                      const double x3,
2447                                      const double x4,
2448                                      const double x5,
2449                                      const double x6,
2450                                      const Num2& w) {
2451     if (dim_ != 7U)
2452       Private::h_badargs("fillC");
2453     double* wg = &weightBuf_[0];
2454     unsigned* idx = &indexBuf_[0];
2455     const Axis* ax = &axes_[0];
2456     const unsigned o0 = ax[0].overflowIndexWeighted(x0, idx + 0, wg + 0);
2457     const unsigned o1 = ax[1].overflowIndexWeighted(x1, idx + 1, wg + 1);
2458     const unsigned o2 = ax[2].overflowIndexWeighted(x2, idx + 2, wg + 2);
2459     const unsigned o3 = ax[3].overflowIndexWeighted(x3, idx + 3, wg + 3);
2460     const unsigned o4 = ax[4].overflowIndexWeighted(x4, idx + 4, wg + 4);
2461     const unsigned o5 = ax[5].overflowIndexWeighted(x5, idx + 5, wg + 5);
2462     const unsigned o6 = ax[6].overflowIndexWeighted(x6, idx + 6, wg + 6);
2463     if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U && o4 == 1U && o5 == 1U && o6 == 1U)
2464       fillPreservingCentroid(w);
2465     else {
2466       overflow_(o0, o1, o2, o3, o4, o5, o6) += w;
2467       ++overCount_;
2468     }
2469     ++fillCount_;
2470     ++modCount_;
2471   }
2472 
2473   template <typename Numeric, class Axis>
2474   const Numeric& HistoND<Numeric, Axis>::examine(const double x0,
2475                                                  const double x1,
2476                                                  const double x2,
2477                                                  const double x3,
2478                                                  const double x4,
2479                                                  const double x5,
2480                                                  const double x6) const {
2481     if (dim_ != 7U)
2482       Private::h_badargs("examine");
2483     unsigned i0 = 0, i1 = 0, i2 = 0, i3 = 0, i4 = 0, i5 = 0, i6 = 0;
2484     const Axis* ax = &axes_[0];
2485     const unsigned o0 = ax[0].overflowIndex(x0, &i0);
2486     const unsigned o1 = ax[1].overflowIndex(x1, &i1);
2487     const unsigned o2 = ax[2].overflowIndex(x2, &i2);
2488     const unsigned o3 = ax[3].overflowIndex(x3, &i3);
2489     const unsigned o4 = ax[4].overflowIndex(x4, &i4);
2490     const unsigned o5 = ax[5].overflowIndex(x5, &i5);
2491     const unsigned o6 = ax[6].overflowIndex(x6, &i6);
2492     if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U && o4 == 1U && o5 == 1U && o6 == 1U)
2493       return data_(i0, i1, i2, i3, i4, i5, i6);
2494     else
2495       return overflow_(o0, o1, o2, o3, o4, o5, o6);
2496   }
2497 
2498   template <typename Numeric, class Axis>
2499   const Numeric& HistoND<Numeric, Axis>::closestBin(const double x0,
2500                                                     const double x1,
2501                                                     const double x2,
2502                                                     const double x3,
2503                                                     const double x4,
2504                                                     const double x5,
2505                                                     const double x6) const {
2506     if (dim_ != 7U)
2507       Private::h_badargs("closestBin");
2508     const Axis* ax = &axes_[0];
2509     const unsigned i0 = ax[0].closestValidBin(x0);
2510     const unsigned i1 = ax[1].closestValidBin(x1);
2511     const unsigned i2 = ax[2].closestValidBin(x2);
2512     const unsigned i3 = ax[3].closestValidBin(x3);
2513     const unsigned i4 = ax[4].closestValidBin(x4);
2514     const unsigned i5 = ax[5].closestValidBin(x5);
2515     const unsigned i6 = ax[6].closestValidBin(x6);
2516     return data_(i0, i1, i2, i3, i4, i5, i6);
2517   }
2518 
2519   template <typename Numeric, class Axis>
2520   template <typename Num2>
2521   void HistoND<Numeric, Axis>::fill(const double x0,
2522                                     const double x1,
2523                                     const double x2,
2524                                     const double x3,
2525                                     const double x4,
2526                                     const double x5,
2527                                     const double x6,
2528                                     const double x7,
2529                                     const Num2& w) {
2530     if (dim_ != 8U)
2531       Private::h_badargs("fill");
2532     unsigned i0 = 0, i1 = 0, i2 = 0, i3 = 0, i4 = 0, i5 = 0, i6 = 0, i7 = 0;
2533     const Axis* ax = &axes_[0];
2534     const unsigned o0 = ax[0].overflowIndex(x0, &i0);
2535     const unsigned o1 = ax[1].overflowIndex(x1, &i1);
2536     const unsigned o2 = ax[2].overflowIndex(x2, &i2);
2537     const unsigned o3 = ax[3].overflowIndex(x3, &i3);
2538     const unsigned o4 = ax[4].overflowIndex(x4, &i4);
2539     const unsigned o5 = ax[5].overflowIndex(x5, &i5);
2540     const unsigned o6 = ax[6].overflowIndex(x6, &i6);
2541     const unsigned o7 = ax[7].overflowIndex(x7, &i7);
2542     if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U && o4 == 1U && o5 == 1U && o6 == 1U && o7 == 1U)
2543       data_(i0, i1, i2, i3, i4, i5, i6, i7) += w;
2544     else {
2545       overflow_(o0, o1, o2, o3, o4, o5, o6, o7) += w;
2546       ++overCount_;
2547     }
2548     ++fillCount_;
2549     ++modCount_;
2550   }
2551 
2552   template <typename Numeric, class Axis>
2553   template <typename Num2, class Functor>
2554   void HistoND<Numeric, Axis>::dispatch(const double x0,
2555                                         const double x1,
2556                                         const double x2,
2557                                         const double x3,
2558                                         const double x4,
2559                                         const double x5,
2560                                         const double x6,
2561                                         const double x7,
2562                                         Num2& w,
2563                                         Functor& f) {
2564     if (dim_ != 8U)
2565       Private::h_badargs("dispatch");
2566     unsigned i0 = 0, i1 = 0, i2 = 0, i3 = 0, i4 = 0, i5 = 0, i6 = 0, i7 = 0;
2567     const Axis* ax = &axes_[0];
2568     const unsigned o0 = ax[0].overflowIndex(x0, &i0);
2569     const unsigned o1 = ax[1].overflowIndex(x1, &i1);
2570     const unsigned o2 = ax[2].overflowIndex(x2, &i2);
2571     const unsigned o3 = ax[3].overflowIndex(x3, &i3);
2572     const unsigned o4 = ax[4].overflowIndex(x4, &i4);
2573     const unsigned o5 = ax[5].overflowIndex(x5, &i5);
2574     const unsigned o6 = ax[6].overflowIndex(x6, &i6);
2575     const unsigned o7 = ax[7].overflowIndex(x7, &i7);
2576     if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U && o4 == 1U && o5 == 1U && o6 == 1U && o7 == 1U)
2577       f(data_(i0, i1, i2, i3, i4, i5, i6, i7), w);
2578     else
2579       f(overflow_(o0, o1, o2, o3, o4, o5, o6, o7), w);
2580     ++modCount_;
2581   }
2582 
2583   template <typename Numeric, class Axis>
2584   template <typename Num2>
2585   void HistoND<Numeric, Axis>::fillC(const double x0,
2586                                      const double x1,
2587                                      const double x2,
2588                                      const double x3,
2589                                      const double x4,
2590                                      const double x5,
2591                                      const double x6,
2592                                      const double x7,
2593                                      const Num2& w) {
2594     if (dim_ != 8U)
2595       Private::h_badargs("fillC");
2596     double* wg = &weightBuf_[0];
2597     unsigned* idx = &indexBuf_[0];
2598     const Axis* ax = &axes_[0];
2599     const unsigned o0 = ax[0].overflowIndexWeighted(x0, idx + 0, wg + 0);
2600     const unsigned o1 = ax[1].overflowIndexWeighted(x1, idx + 1, wg + 1);
2601     const unsigned o2 = ax[2].overflowIndexWeighted(x2, idx + 2, wg + 2);
2602     const unsigned o3 = ax[3].overflowIndexWeighted(x3, idx + 3, wg + 3);
2603     const unsigned o4 = ax[4].overflowIndexWeighted(x4, idx + 4, wg + 4);
2604     const unsigned o5 = ax[5].overflowIndexWeighted(x5, idx + 5, wg + 5);
2605     const unsigned o6 = ax[6].overflowIndexWeighted(x6, idx + 6, wg + 6);
2606     const unsigned o7 = ax[7].overflowIndexWeighted(x7, idx + 7, wg + 7);
2607     if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U && o4 == 1U && o5 == 1U && o6 == 1U && o7 == 1U)
2608       fillPreservingCentroid(w);
2609     else {
2610       overflow_(o0, o1, o2, o3, o4, o5, o6, o7) += w;
2611       ++overCount_;
2612     }
2613     ++fillCount_;
2614     ++modCount_;
2615   }
2616 
2617   template <typename Numeric, class Axis>
2618   const Numeric& HistoND<Numeric, Axis>::examine(const double x0,
2619                                                  const double x1,
2620                                                  const double x2,
2621                                                  const double x3,
2622                                                  const double x4,
2623                                                  const double x5,
2624                                                  const double x6,
2625                                                  const double x7) const {
2626     if (dim_ != 8U)
2627       Private::h_badargs("examine");
2628     unsigned i0 = 0, i1 = 0, i2 = 0, i3 = 0, i4 = 0, i5 = 0, i6 = 0, i7 = 0;
2629     const Axis* ax = &axes_[0];
2630     const unsigned o0 = ax[0].overflowIndex(x0, &i0);
2631     const unsigned o1 = ax[1].overflowIndex(x1, &i1);
2632     const unsigned o2 = ax[2].overflowIndex(x2, &i2);
2633     const unsigned o3 = ax[3].overflowIndex(x3, &i3);
2634     const unsigned o4 = ax[4].overflowIndex(x4, &i4);
2635     const unsigned o5 = ax[5].overflowIndex(x5, &i5);
2636     const unsigned o6 = ax[6].overflowIndex(x6, &i6);
2637     const unsigned o7 = ax[7].overflowIndex(x7, &i7);
2638     if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U && o4 == 1U && o5 == 1U && o6 == 1U && o7 == 1U)
2639       return data_(i0, i1, i2, i3, i4, i5, i6, i7);
2640     else
2641       return overflow_(o0, o1, o2, o3, o4, o5, o6, o7);
2642   }
2643 
2644   template <typename Numeric, class Axis>
2645   const Numeric& HistoND<Numeric, Axis>::closestBin(const double x0,
2646                                                     const double x1,
2647                                                     const double x2,
2648                                                     const double x3,
2649                                                     const double x4,
2650                                                     const double x5,
2651                                                     const double x6,
2652                                                     const double x7) const {
2653     if (dim_ != 8U)
2654       Private::h_badargs("closestBin");
2655     const Axis* ax = &axes_[0];
2656     const unsigned i0 = ax[0].closestValidBin(x0);
2657     const unsigned i1 = ax[1].closestValidBin(x1);
2658     const unsigned i2 = ax[2].closestValidBin(x2);
2659     const unsigned i3 = ax[3].closestValidBin(x3);
2660     const unsigned i4 = ax[4].closestValidBin(x4);
2661     const unsigned i5 = ax[5].closestValidBin(x5);
2662     const unsigned i6 = ax[6].closestValidBin(x6);
2663     const unsigned i7 = ax[7].closestValidBin(x7);
2664     return data_(i0, i1, i2, i3, i4, i5, i6, i7);
2665   }
2666 
2667   template <typename Numeric, class Axis>
2668   template <typename Num2>
2669   void HistoND<Numeric, Axis>::fill(const double x0,
2670                                     const double x1,
2671                                     const double x2,
2672                                     const double x3,
2673                                     const double x4,
2674                                     const double x5,
2675                                     const double x6,
2676                                     const double x7,
2677                                     const double x8,
2678                                     const Num2& w) {
2679     if (dim_ != 9U)
2680       Private::h_badargs("fill");
2681     unsigned i0 = 0, i1 = 0, i2 = 0, i3 = 0, i4 = 0, i5 = 0, i6 = 0, i7 = 0, i8 = 0;
2682     const Axis* ax = &axes_[0];
2683     const unsigned o0 = ax[0].overflowIndex(x0, &i0);
2684     const unsigned o1 = ax[1].overflowIndex(x1, &i1);
2685     const unsigned o2 = ax[2].overflowIndex(x2, &i2);
2686     const unsigned o3 = ax[3].overflowIndex(x3, &i3);
2687     const unsigned o4 = ax[4].overflowIndex(x4, &i4);
2688     const unsigned o5 = ax[5].overflowIndex(x5, &i5);
2689     const unsigned o6 = ax[6].overflowIndex(x6, &i6);
2690     const unsigned o7 = ax[7].overflowIndex(x7, &i7);
2691     const unsigned o8 = ax[8].overflowIndex(x8, &i8);
2692     if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U && o4 == 1U && o5 == 1U && o6 == 1U && o7 == 1U && o8 == 1U)
2693       data_(i0, i1, i2, i3, i4, i5, i6, i7, i8) += w;
2694     else {
2695       overflow_(o0, o1, o2, o3, o4, o5, o6, o7, o8) += w;
2696       ++overCount_;
2697     }
2698     ++fillCount_;
2699     ++modCount_;
2700   }
2701 
2702   template <typename Numeric, class Axis>
2703   template <typename Num2, class Functor>
2704   void HistoND<Numeric, Axis>::dispatch(const double x0,
2705                                         const double x1,
2706                                         const double x2,
2707                                         const double x3,
2708                                         const double x4,
2709                                         const double x5,
2710                                         const double x6,
2711                                         const double x7,
2712                                         const double x8,
2713                                         Num2& w,
2714                                         Functor& f) {
2715     if (dim_ != 9U)
2716       Private::h_badargs("dispatch");
2717     unsigned i0 = 0, i1 = 0, i2 = 0, i3 = 0, i4 = 0, i5 = 0, i6 = 0, i7 = 0, i8 = 0;
2718     const Axis* ax = &axes_[0];
2719     const unsigned o0 = ax[0].overflowIndex(x0, &i0);
2720     const unsigned o1 = ax[1].overflowIndex(x1, &i1);
2721     const unsigned o2 = ax[2].overflowIndex(x2, &i2);
2722     const unsigned o3 = ax[3].overflowIndex(x3, &i3);
2723     const unsigned o4 = ax[4].overflowIndex(x4, &i4);
2724     const unsigned o5 = ax[5].overflowIndex(x5, &i5);
2725     const unsigned o6 = ax[6].overflowIndex(x6, &i6);
2726     const unsigned o7 = ax[7].overflowIndex(x7, &i7);
2727     const unsigned o8 = ax[8].overflowIndex(x8, &i8);
2728     if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U && o4 == 1U && o5 == 1U && o6 == 1U && o7 == 1U && o8 == 1U)
2729       f(data_(i0, i1, i2, i3, i4, i5, i6, i7, i8), w);
2730     else
2731       f(overflow_(o0, o1, o2, o3, o4, o5, o6, o7, o8), w);
2732     ++modCount_;
2733   }
2734 
2735   template <typename Numeric, class Axis>
2736   template <typename Num2>
2737   void HistoND<Numeric, Axis>::fillC(const double x0,
2738                                      const double x1,
2739                                      const double x2,
2740                                      const double x3,
2741                                      const double x4,
2742                                      const double x5,
2743                                      const double x6,
2744                                      const double x7,
2745                                      const double x8,
2746                                      const Num2& w) {
2747     if (dim_ != 9U)
2748       Private::h_badargs("fillC");
2749     double* wg = &weightBuf_[0];
2750     unsigned* idx = &indexBuf_[0];
2751     const Axis* ax = &axes_[0];
2752     const unsigned o0 = ax[0].overflowIndexWeighted(x0, idx + 0, wg + 0);
2753     const unsigned o1 = ax[1].overflowIndexWeighted(x1, idx + 1, wg + 1);
2754     const unsigned o2 = ax[2].overflowIndexWeighted(x2, idx + 2, wg + 2);
2755     const unsigned o3 = ax[3].overflowIndexWeighted(x3, idx + 3, wg + 3);
2756     const unsigned o4 = ax[4].overflowIndexWeighted(x4, idx + 4, wg + 4);
2757     const unsigned o5 = ax[5].overflowIndexWeighted(x5, idx + 5, wg + 5);
2758     const unsigned o6 = ax[6].overflowIndexWeighted(x6, idx + 6, wg + 6);
2759     const unsigned o7 = ax[7].overflowIndexWeighted(x7, idx + 7, wg + 7);
2760     const unsigned o8 = ax[8].overflowIndexWeighted(x8, idx + 8, wg + 8);
2761     if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U && o4 == 1U && o5 == 1U && o6 == 1U && o7 == 1U && o8 == 1U)
2762       fillPreservingCentroid(w);
2763     else {
2764       overflow_(o0, o1, o2, o3, o4, o5, o6, o7, o8) += w;
2765       ++overCount_;
2766     }
2767     ++fillCount_;
2768     ++modCount_;
2769   }
2770 
2771   template <typename Numeric, class Axis>
2772   const Numeric& HistoND<Numeric, Axis>::examine(const double x0,
2773                                                  const double x1,
2774                                                  const double x2,
2775                                                  const double x3,
2776                                                  const double x4,
2777                                                  const double x5,
2778                                                  const double x6,
2779                                                  const double x7,
2780                                                  const double x8) const {
2781     if (dim_ != 9U)
2782       Private::h_badargs("examine");
2783     unsigned i0 = 0, i1 = 0, i2 = 0, i3 = 0, i4 = 0, i5 = 0, i6 = 0, i7 = 0, i8 = 0;
2784     const Axis* ax = &axes_[0];
2785     const unsigned o0 = ax[0].overflowIndex(x0, &i0);
2786     const unsigned o1 = ax[1].overflowIndex(x1, &i1);
2787     const unsigned o2 = ax[2].overflowIndex(x2, &i2);
2788     const unsigned o3 = ax[3].overflowIndex(x3, &i3);
2789     const unsigned o4 = ax[4].overflowIndex(x4, &i4);
2790     const unsigned o5 = ax[5].overflowIndex(x5, &i5);
2791     const unsigned o6 = ax[6].overflowIndex(x6, &i6);
2792     const unsigned o7 = ax[7].overflowIndex(x7, &i7);
2793     const unsigned o8 = ax[8].overflowIndex(x8, &i8);
2794     if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U && o4 == 1U && o5 == 1U && o6 == 1U && o7 == 1U && o8 == 1U)
2795       return data_(i0, i1, i2, i3, i4, i5, i6, i7, i8);
2796     else
2797       return overflow_(o0, o1, o2, o3, o4, o5, o6, o7, o8);
2798   }
2799 
2800   template <typename Numeric, class Axis>
2801   const Numeric& HistoND<Numeric, Axis>::closestBin(const double x0,
2802                                                     const double x1,
2803                                                     const double x2,
2804                                                     const double x3,
2805                                                     const double x4,
2806                                                     const double x5,
2807                                                     const double x6,
2808                                                     const double x7,
2809                                                     const double x8) const {
2810     if (dim_ != 9U)
2811       Private::h_badargs("closestBin");
2812     const Axis* ax = &axes_[0];
2813     const unsigned i0 = ax[0].closestValidBin(x0);
2814     const unsigned i1 = ax[1].closestValidBin(x1);
2815     const unsigned i2 = ax[2].closestValidBin(x2);
2816     const unsigned i3 = ax[3].closestValidBin(x3);
2817     const unsigned i4 = ax[4].closestValidBin(x4);
2818     const unsigned i5 = ax[5].closestValidBin(x5);
2819     const unsigned i6 = ax[6].closestValidBin(x6);
2820     const unsigned i7 = ax[7].closestValidBin(x7);
2821     const unsigned i8 = ax[8].closestValidBin(x8);
2822     return data_(i0, i1, i2, i3, i4, i5, i6, i7, i8);
2823   }
2824 
2825   template <typename Numeric, class Axis>
2826   template <typename Num2>
2827   void HistoND<Numeric, Axis>::fill(const double x0,
2828                                     const double x1,
2829                                     const double x2,
2830                                     const double x3,
2831                                     const double x4,
2832                                     const double x5,
2833                                     const double x6,
2834                                     const double x7,
2835                                     const double x8,
2836                                     const double x9,
2837                                     const Num2& w) {
2838     if (dim_ != 10U)
2839       Private::h_badargs("fill");
2840     unsigned i0 = 0, i1 = 0, i2 = 0, i3 = 0, i4 = 0, i5 = 0, i6 = 0, i7 = 0, i8 = 0, i9 = 0;
2841     const Axis* ax = &axes_[0];
2842     const unsigned o0 = ax[0].overflowIndex(x0, &i0);
2843     const unsigned o1 = ax[1].overflowIndex(x1, &i1);
2844     const unsigned o2 = ax[2].overflowIndex(x2, &i2);
2845     const unsigned o3 = ax[3].overflowIndex(x3, &i3);
2846     const unsigned o4 = ax[4].overflowIndex(x4, &i4);
2847     const unsigned o5 = ax[5].overflowIndex(x5, &i5);
2848     const unsigned o6 = ax[6].overflowIndex(x6, &i6);
2849     const unsigned o7 = ax[7].overflowIndex(x7, &i7);
2850     const unsigned o8 = ax[8].overflowIndex(x8, &i8);
2851     const unsigned o9 = ax[9].overflowIndex(x9, &i9);
2852     if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U && o4 == 1U && o5 == 1U && o6 == 1U && o7 == 1U && o8 == 1U &&
2853         o9 == 1U)
2854       data_(i0, i1, i2, i3, i4, i5, i6, i7, i8, i9) += w;
2855     else {
2856       overflow_(o0, o1, o2, o3, o4, o5, o6, o7, o8, o9) += w;
2857       ++overCount_;
2858     }
2859     ++fillCount_;
2860     ++modCount_;
2861   }
2862 
2863   template <typename Numeric, class Axis>
2864   template <typename Num2, class Functor>
2865   void HistoND<Numeric, Axis>::dispatch(const double x0,
2866                                         const double x1,
2867                                         const double x2,
2868                                         const double x3,
2869                                         const double x4,
2870                                         const double x5,
2871                                         const double x6,
2872                                         const double x7,
2873                                         const double x8,
2874                                         const double x9,
2875                                         Num2& w,
2876                                         Functor& f) {
2877     if (dim_ != 10U)
2878       Private::h_badargs("dispatch");
2879     unsigned i0 = 0, i1 = 0, i2 = 0, i3 = 0, i4 = 0, i5 = 0, i6 = 0, i7 = 0, i8 = 0, i9 = 0;
2880     const Axis* ax = &axes_[0];
2881     const unsigned o0 = ax[0].overflowIndex(x0, &i0);
2882     const unsigned o1 = ax[1].overflowIndex(x1, &i1);
2883     const unsigned o2 = ax[2].overflowIndex(x2, &i2);
2884     const unsigned o3 = ax[3].overflowIndex(x3, &i3);
2885     const unsigned o4 = ax[4].overflowIndex(x4, &i4);
2886     const unsigned o5 = ax[5].overflowIndex(x5, &i5);
2887     const unsigned o6 = ax[6].overflowIndex(x6, &i6);
2888     const unsigned o7 = ax[7].overflowIndex(x7, &i7);
2889     const unsigned o8 = ax[8].overflowIndex(x8, &i8);
2890     const unsigned o9 = ax[9].overflowIndex(x9, &i9);
2891     if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U && o4 == 1U && o5 == 1U && o6 == 1U && o7 == 1U && o8 == 1U &&
2892         o9 == 1U)
2893       f(data_(i0, i1, i2, i3, i4, i5, i6, i7, i8, i9), w);
2894     else
2895       f(overflow_(o0, o1, o2, o3, o4, o5, o6, o7, o8, o9), w);
2896     ++modCount_;
2897   }
2898 
2899   template <typename Numeric, class Axis>
2900   template <typename Num2>
2901   void HistoND<Numeric, Axis>::fillC(const double x0,
2902                                      const double x1,
2903                                      const double x2,
2904                                      const double x3,
2905                                      const double x4,
2906                                      const double x5,
2907                                      const double x6,
2908                                      const double x7,
2909                                      const double x8,
2910                                      const double x9,
2911                                      const Num2& w) {
2912     if (dim_ != 10U)
2913       Private::h_badargs("fillC");
2914     double* wg = &weightBuf_[0];
2915     unsigned* idx = &indexBuf_[0];
2916     const Axis* ax = &axes_[0];
2917     const unsigned o0 = ax[0].overflowIndexWeighted(x0, idx + 0, wg + 0);
2918     const unsigned o1 = ax[1].overflowIndexWeighted(x1, idx + 1, wg + 1);
2919     const unsigned o2 = ax[2].overflowIndexWeighted(x2, idx + 2, wg + 2);
2920     const unsigned o3 = ax[3].overflowIndexWeighted(x3, idx + 3, wg + 3);
2921     const unsigned o4 = ax[4].overflowIndexWeighted(x4, idx + 4, wg + 4);
2922     const unsigned o5 = ax[5].overflowIndexWeighted(x5, idx + 5, wg + 5);
2923     const unsigned o6 = ax[6].overflowIndexWeighted(x6, idx + 6, wg + 6);
2924     const unsigned o7 = ax[7].overflowIndexWeighted(x7, idx + 7, wg + 7);
2925     const unsigned o8 = ax[8].overflowIndexWeighted(x8, idx + 8, wg + 8);
2926     const unsigned o9 = ax[9].overflowIndexWeighted(x9, idx + 9, wg + 9);
2927     if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U && o4 == 1U && o5 == 1U && o6 == 1U && o7 == 1U && o8 == 1U &&
2928         o9 == 1U)
2929       fillPreservingCentroid(w);
2930     else {
2931       overflow_(o0, o1, o2, o3, o4, o5, o6, o7, o8, o9) += w;
2932       ++overCount_;
2933     }
2934     ++fillCount_;
2935     ++modCount_;
2936   }
2937 
2938   template <typename Numeric, class Axis>
2939   const Numeric& HistoND<Numeric, Axis>::examine(const double x0,
2940                                                  const double x1,
2941                                                  const double x2,
2942                                                  const double x3,
2943                                                  const double x4,
2944                                                  const double x5,
2945                                                  const double x6,
2946                                                  const double x7,
2947                                                  const double x8,
2948                                                  const double x9) const {
2949     if (dim_ != 10U)
2950       Private::h_badargs("examine");
2951     unsigned i0 = 0, i1 = 0, i2 = 0, i3 = 0, i4 = 0, i5 = 0, i6 = 0, i7 = 0, i8 = 0, i9 = 0;
2952     const Axis* ax = &axes_[0];
2953     const unsigned o0 = ax[0].overflowIndex(x0, &i0);
2954     const unsigned o1 = ax[1].overflowIndex(x1, &i1);
2955     const unsigned o2 = ax[2].overflowIndex(x2, &i2);
2956     const unsigned o3 = ax[3].overflowIndex(x3, &i3);
2957     const unsigned o4 = ax[4].overflowIndex(x4, &i4);
2958     const unsigned o5 = ax[5].overflowIndex(x5, &i5);
2959     const unsigned o6 = ax[6].overflowIndex(x6, &i6);
2960     const unsigned o7 = ax[7].overflowIndex(x7, &i7);
2961     const unsigned o8 = ax[8].overflowIndex(x8, &i8);
2962     const unsigned o9 = ax[9].overflowIndex(x9, &i9);
2963     if (o0 == 1U && o1 == 1U && o2 == 1U && o3 == 1U && o4 == 1U && o5 == 1U && o6 == 1U && o7 == 1U && o8 == 1U &&
2964         o9 == 1U)
2965       return data_(i0, i1, i2, i3, i4, i5, i6, i7, i8, i9);
2966     else
2967       return overflow_(o0, o1, o2, o3, o4, o5, o6, o7, o8, o9);
2968   }
2969 
2970   template <typename Numeric, class Axis>
2971   const Numeric& HistoND<Numeric, Axis>::closestBin(const double x0,
2972                                                     const double x1,
2973                                                     const double x2,
2974                                                     const double x3,
2975                                                     const double x4,
2976                                                     const double x5,
2977                                                     const double x6,
2978                                                     const double x7,
2979                                                     const double x8,
2980                                                     const double x9) const {
2981     if (dim_ != 10U)
2982       Private::h_badargs("closestBin");
2983     const Axis* ax = &axes_[0];
2984     const unsigned i0 = ax[0].closestValidBin(x0);
2985     const unsigned i1 = ax[1].closestValidBin(x1);
2986     const unsigned i2 = ax[2].closestValidBin(x2);
2987     const unsigned i3 = ax[3].closestValidBin(x3);
2988     const unsigned i4 = ax[4].closestValidBin(x4);
2989     const unsigned i5 = ax[5].closestValidBin(x5);
2990     const unsigned i6 = ax[6].closestValidBin(x6);
2991     const unsigned i7 = ax[7].closestValidBin(x7);
2992     const unsigned i8 = ax[8].closestValidBin(x8);
2993     const unsigned i9 = ax[9].closestValidBin(x9);
2994     return data_(i0, i1, i2, i3, i4, i5, i6, i7, i8, i9);
2995   }
2996 
2997   template <typename Numeric, class Axis>
2998   template <typename Num2>
2999   inline void HistoND<Numeric, Axis>::setBin(const unsigned* index, const unsigned indexLen, const Num2& v) {
3000     data_.value(index, indexLen) = v;
3001     ++modCount_;
3002   }
3003 
3004   template <typename Numeric, class Axis>
3005   template <typename Num2>
3006   inline void HistoND<Numeric, Axis>::setBinAt(const unsigned* index, const unsigned indexLen, const Num2& v) {
3007     data_.valueAt(index, indexLen) = v;
3008     ++modCount_;
3009   }
3010 
3011   template <typename Numeric, class Axis>
3012   template <typename Num2>
3013   inline void HistoND<Numeric, Axis>::setBin(const Num2& v) {
3014     data_() = v;
3015     ++modCount_;
3016   }
3017 
3018   template <typename Numeric, class Axis>
3019   template <typename Num2>
3020   inline void HistoND<Numeric, Axis>::setBinAt(const Num2& v) {
3021     data_.at() = v;
3022     ++modCount_;
3023   }
3024 
3025   template <typename Numeric, class Axis>
3026   template <typename Num2>
3027   inline void HistoND<Numeric, Axis>::setBin(const unsigned i0, const Num2& v) {
3028     data_(i0) = v;
3029     ++modCount_;
3030   }
3031 
3032   template <typename Numeric, class Axis>
3033   template <typename Num2>
3034   inline void HistoND<Numeric, Axis>::setBinAt(const unsigned i0, const Num2& v) {
3035     data_.at(i0) = v;
3036     ++modCount_;
3037   }
3038 
3039   template <typename Numeric, class Axis>
3040   template <typename Num2>
3041   inline void HistoND<Numeric, Axis>::setBin(const unsigned i0, const unsigned i1, const Num2& v) {
3042     data_(i0, i1) = v;
3043     ++modCount_;
3044   }
3045 
3046   template <typename Numeric, class Axis>
3047   template <typename Num2>
3048   inline void HistoND<Numeric, Axis>::setBinAt(const unsigned i0, const unsigned i1, const Num2& v) {
3049     data_.at(i0, i1) = v;
3050     ++modCount_;
3051   }
3052 
3053   template <typename Numeric, class Axis>
3054   template <typename Num2>
3055   inline void HistoND<Numeric, Axis>::setBin(const unsigned i0, const unsigned i1, const unsigned i2, const Num2& v) {
3056     data_(i0, i1, i2) = v;
3057     ++modCount_;
3058   }
3059 
3060   template <typename Numeric, class Axis>
3061   template <typename Num2>
3062   inline void HistoND<Numeric, Axis>::setBin(
3063       const unsigned i0, const unsigned i1, const unsigned i2, const unsigned i3, const Num2& v) {
3064     data_(i0, i1, i2, i3) = v;
3065     ++modCount_;
3066   }
3067 
3068   template <typename Numeric, class Axis>
3069   template <typename Num2>
3070   inline void HistoND<Numeric, Axis>::setBin(
3071       const unsigned i0, const unsigned i1, const unsigned i2, const unsigned i3, const unsigned i4, const Num2& v) {
3072     data_(i0, i1, i2, i3, i4) = v;
3073     ++modCount_;
3074   }
3075 
3076   template <typename Numeric, class Axis>
3077   template <typename Num2>
3078   inline void HistoND<Numeric, Axis>::setBin(const unsigned i0,
3079                                              const unsigned i1,
3080                                              const unsigned i2,
3081                                              const unsigned i3,
3082                                              const unsigned i4,
3083                                              const unsigned i5,
3084                                              const Num2& v) {
3085     data_(i0, i1, i2, i3, i4, i5) = v;
3086     ++modCount_;
3087   }
3088 
3089   template <typename Numeric, class Axis>
3090   template <typename Num2>
3091   inline void HistoND<Numeric, Axis>::setBin(const unsigned i0,
3092                                              const unsigned i1,
3093                                              const unsigned i2,
3094                                              const unsigned i3,
3095                                              const unsigned i4,
3096                                              const unsigned i5,
3097                                              const unsigned i6,
3098                                              const Num2& v) {
3099     data_(i0, i1, i2, i3, i4, i5, i6) = v;
3100     ++modCount_;
3101   }
3102 
3103   template <typename Numeric, class Axis>
3104   template <typename Num2>
3105   inline void HistoND<Numeric, Axis>::setBin(const unsigned i0,
3106                                              const unsigned i1,
3107                                              const unsigned i2,
3108                                              const unsigned i3,
3109                                              const unsigned i4,
3110                                              const unsigned i5,
3111                                              const unsigned i6,
3112                                              const unsigned i7,
3113                                              const Num2& v) {
3114     data_(i0, i1, i2, i3, i4, i5, i6, i7) = v;
3115     ++modCount_;
3116   }
3117 
3118   template <typename Numeric, class Axis>
3119   template <typename Num2>
3120   inline void HistoND<Numeric, Axis>::setBin(const unsigned i0,
3121                                              const unsigned i1,
3122                                              const unsigned i2,
3123                                              const unsigned i3,
3124                                              const unsigned i4,
3125                                              const unsigned i5,
3126                                              const unsigned i6,
3127                                              const unsigned i7,
3128                                              const unsigned i8,
3129                                              const Num2& v) {
3130     data_(i0, i1, i2, i3, i4, i5, i6, i7, i8) = v;
3131     ++modCount_;
3132   }
3133 
3134   template <typename Numeric, class Axis>
3135   template <typename Num2>
3136   inline void HistoND<Numeric, Axis>::setBin(const unsigned i0,
3137                                              const unsigned i1,
3138                                              const unsigned i2,
3139                                              const unsigned i3,
3140                                              const unsigned i4,
3141                                              const unsigned i5,
3142                                              const unsigned i6,
3143                                              const unsigned i7,
3144                                              const unsigned i8,
3145                                              const unsigned i9,
3146                                              const Num2& v) {
3147     data_(i0, i1, i2, i3, i4, i5, i6, i7, i8, i9) = 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, const unsigned i1, const unsigned i2, const Num2& v) {
3154     data_.at(i0, i1, i2) = v;
3155     ++modCount_;
3156   }
3157 
3158   template <typename Numeric, class Axis>
3159   template <typename Num2>
3160   inline void HistoND<Numeric, Axis>::setBinAt(
3161       const unsigned i0, const unsigned i1, const unsigned i2, const unsigned i3, const Num2& v) {
3162     data_.at(i0, i1, i2, i3) = v;
3163     ++modCount_;
3164   }
3165 
3166   template <typename Numeric, class Axis>
3167   template <typename Num2>
3168   inline void HistoND<Numeric, Axis>::setBinAt(
3169       const unsigned i0, const unsigned i1, const unsigned i2, const unsigned i3, const unsigned i4, const Num2& v) {
3170     data_.at(i0, i1, i2, i3, i4) = v;
3171     ++modCount_;
3172   }
3173 
3174   template <typename Numeric, class Axis>
3175   template <typename Num2>
3176   inline void HistoND<Numeric, Axis>::setBinAt(const unsigned i0,
3177                                                const unsigned i1,
3178                                                const unsigned i2,
3179                                                const unsigned i3,
3180                                                const unsigned i4,
3181                                                const unsigned i5,
3182                                                const Num2& v) {
3183     data_.at(i0, i1, i2, i3, i4, i5) = v;
3184     ++modCount_;
3185   }
3186 
3187   template <typename Numeric, class Axis>
3188   template <typename Num2>
3189   inline void HistoND<Numeric, Axis>::setBinAt(const unsigned i0,
3190                                                const unsigned i1,
3191                                                const unsigned i2,
3192                                                const unsigned i3,
3193                                                const unsigned i4,
3194                                                const unsigned i5,
3195                                                const unsigned i6,
3196                                                const Num2& v) {
3197     data_.at(i0, i1, i2, i3, i4, i5, i6) = v;
3198     ++modCount_;
3199   }
3200 
3201   template <typename Numeric, class Axis>
3202   template <typename Num2>
3203   inline void HistoND<Numeric, Axis>::setBinAt(const unsigned i0,
3204                                                const unsigned i1,
3205                                                const unsigned i2,
3206                                                const unsigned i3,
3207                                                const unsigned i4,
3208                                                const unsigned i5,
3209                                                const unsigned i6,
3210                                                const unsigned i7,
3211                                                const Num2& v) {
3212     data_.at(i0, i1, i2, i3, i4, i5, i6, i7) = v;
3213     ++modCount_;
3214   }
3215 
3216   template <typename Numeric, class Axis>
3217   template <typename Num2>
3218   inline void HistoND<Numeric, Axis>::setBinAt(const unsigned i0,
3219                                                const unsigned i1,
3220                                                const unsigned i2,
3221                                                const unsigned i3,
3222                                                const unsigned i4,
3223                                                const unsigned i5,
3224                                                const unsigned i6,
3225                                                const unsigned i7,
3226                                                const unsigned i8,
3227                                                const Num2& v) {
3228     data_.at(i0, i1, i2, i3, i4, i5, i6, i7, i8) = v;
3229     ++modCount_;
3230   }
3231 
3232   template <typename Numeric, class Axis>
3233   template <typename Num2>
3234   inline void HistoND<Numeric, Axis>::setBinAt(const unsigned i0,
3235                                                const unsigned i1,
3236                                                const unsigned i2,
3237                                                const unsigned i3,
3238                                                const unsigned i4,
3239                                                const unsigned i5,
3240                                                const unsigned i6,
3241                                                const unsigned i7,
3242                                                const unsigned i8,
3243                                                const unsigned i9,
3244                                                const Num2& v) {
3245     data_.at(i0, i1, i2, i3, i4, i5, i6, i7, i8, i9) = v;
3246     ++modCount_;
3247   }
3248 
3249   template <typename Numeric, class Axis>
3250   template <typename Num2>
3251   inline HistoND<Numeric, Axis>& HistoND<Numeric, Axis>::operator+=(const HistoND<Num2, Axis>& r) {
3252     data_ += r.data_;
3253     overflow_ += r.overflow_;
3254     fillCount_ += r.fillCount_;
3255     overCount_ += r.overCount_;
3256     ++modCount_;
3257     return *this;
3258   }
3259 
3260   template <typename Numeric, class Axis>
3261   template <typename Num2>
3262   inline HistoND<Numeric, Axis>& HistoND<Numeric, Axis>::operator-=(const HistoND<Num2, Axis>& r) {
3263     data_ -= r.data_;
3264     overflow_ -= r.overflow_;
3265 
3266     // Subtraction does not make much sense for fill counts.
3267     // We will assume that what we want should be equivalent
3268     // to the in-place multiplication of the other histogram
3269     // by -1 and then adding.
3270     //
3271     fillCount_ += r.fillCount_;
3272     overCount_ += r.overCount_;
3273 
3274     ++modCount_;
3275     return *this;
3276   }
3277 
3278   template <typename Numeric, class Axis>
3279   template <typename Num2>
3280   inline HistoND<Numeric, Axis>& HistoND<Numeric, Axis>::operator*=(const Num2& r) {
3281     data_ *= r;
3282     overflow_ *= r;
3283     ++modCount_;
3284     return *this;
3285   }
3286 
3287   template <typename Numeric, class Axis>
3288   template <typename Num2>
3289   inline HistoND<Numeric, Axis>& HistoND<Numeric, Axis>::operator/=(const Num2& r) {
3290     data_ /= r;
3291     overflow_ /= r;
3292     ++modCount_;
3293     return *this;
3294   }
3295 
3296   template <typename Numeric, class Axis>
3297   HistoND<Numeric, Axis>::HistoND(const HistoND& r, const unsigned ax1, const unsigned ax2)
3298       : title_(r.title_),
3299         accumulatedDataLabel_(r.accumulatedDataLabel_),
3300         data_(r.data_.transpose(ax1, ax2)),
3301         overflow_(r.overflow_.transpose(ax1, ax2)),
3302         axes_(r.axes_),
3303         weightBuf_(r.weightBuf_),
3304         indexBuf_(r.indexBuf_),
3305         fillCount_(r.fillCount_),
3306         overCount_(r.overCount_),
3307         modCount_(0UL),
3308         dim_(r.dim_) {
3309     std::swap(axes_[ax1], axes_[ax2]);
3310   }
3311 
3312   template <typename Numeric, class Axis>
3313   HistoND<Numeric, Axis>::HistoND(const HistoND& r)
3314       : title_(r.title_),
3315         accumulatedDataLabel_(r.accumulatedDataLabel_),
3316         data_(r.data_),
3317         overflow_(r.overflow_),
3318         axes_(r.axes_),
3319         weightBuf_(r.weightBuf_),
3320         indexBuf_(r.indexBuf_),
3321         fillCount_(r.fillCount_),
3322         overCount_(r.overCount_),
3323         modCount_(0UL),
3324         dim_(r.dim_) {}
3325 
3326   template <typename Numeric, class Axis>
3327   HistoND<Numeric, Axis>& HistoND<Numeric, Axis>::operator=(const HistoND& r) {
3328     if (&r != this) {
3329       title_ = r.title_;
3330       accumulatedDataLabel_ = r.accumulatedDataLabel_;
3331       data_.uninitialize();
3332       data_ = r.data_;
3333       overflow_.uninitialize();
3334       overflow_ = r.overflow_;
3335       axes_ = r.axes_;
3336       weightBuf_ = r.weightBuf_;
3337       indexBuf_ = r.indexBuf_;
3338       fillCount_ = r.fillCount_;
3339       overCount_ = r.overCount_;
3340       dim_ = r.dim_;
3341       ++modCount_;
3342     }
3343     return *this;
3344   }
3345 
3346   template <typename Numeric, class Axis>
3347   HistoND<Numeric, Axis> HistoND<Numeric, Axis>::transpose(const unsigned axisNum1, const unsigned axisNum2) const {
3348     if (axisNum1 >= dim_ || axisNum2 >= dim_)
3349       throw npstat::NpstatOutOfRange(
3350           "In npstat::HistoND::transpose: "
3351           "axis number is out of range");
3352     if (axisNum1 == axisNum2)
3353       // Just make a copy
3354       return *this;
3355     else
3356       return HistoND(*this, axisNum1, axisNum2);
3357   }
3358 
3359   template <typename Numeric, class Axis>
3360   template <typename Num2>
3361   void HistoND<Numeric, Axis>::addToBinContents(const Num2& weight) {
3362     const unsigned long nDat = data_.length();
3363     Numeric* data = const_cast<Numeric*>(data_.data());
3364     for (unsigned long i = 0; i < nDat; ++i)
3365       data[i] += weight;
3366     fillCount_ += nDat;
3367     ++modCount_;
3368   }
3369 
3370   template <typename Numeric, class Axis>
3371   template <typename Num2>
3372   void HistoND<Numeric, Axis>::addToOverflows(const Num2& weight) {
3373     const unsigned long nOver = overflow_.length();
3374     Numeric* data = const_cast<Numeric*>(overflow_.data());
3375     for (unsigned long i = 0; i < nOver; ++i)
3376       data[i] += weight;
3377     overCount_ += nOver;
3378     fillCount_ += nOver;
3379     ++modCount_;
3380   }
3381 
3382   template <typename Numeric, class Axis>
3383   template <typename Num2>
3384   void HistoND<Numeric, Axis>::addToBinContents(const Num2* data, const unsigned long dataLength) {
3385     if (dataLength != data_.length())
3386       throw npstat::NpstatInvalidArgument("In npstat::HistoND::addToBinContents: incompatible data length");
3387     assert(data);
3388     Numeric* dat = const_cast<Numeric*>(data_.data());
3389     for (unsigned long i = 0; i < dataLength; ++i)
3390       dat[i] += data[i];
3391     fillCount_ += dataLength;
3392     ++modCount_;
3393   }
3394 
3395   template <typename Numeric, class Axis>
3396   template <typename Num2>
3397   void HistoND<Numeric, Axis>::addToOverflows(const Num2* data, const unsigned long dataLength) {
3398     if (dataLength != overflow_.length())
3399       throw npstat::NpstatInvalidArgument("In npstat::HistoND::addToOverflows: incompatible data length");
3400     assert(data);
3401     Numeric* dat = const_cast<Numeric*>(overflow_.data());
3402     for (unsigned long i = 0; i < dataLength; ++i)
3403       dat[i] += data[i];
3404     overCount_ += dataLength;
3405     fillCount_ += dataLength;
3406     ++modCount_;
3407   }
3408 
3409   template <typename Numeric, class Axis>
3410   template <typename Num2>
3411   void HistoND<Numeric, Axis>::scaleBinContents(const Num2* data, const unsigned long dataLength) {
3412     if (dataLength != data_.length())
3413       throw npstat::NpstatInvalidArgument("In npstat::HistoND::scaleBinContents: incompatible data length");
3414     assert(data);
3415     Numeric* dat = const_cast<Numeric*>(data_.data());
3416     for (unsigned long i = 0; i < dataLength; ++i)
3417       dat[i] *= data[i];
3418     ++modCount_;
3419   }
3420 
3421   template <typename Numeric, class Axis>
3422   template <typename Num2>
3423   void HistoND<Numeric, Axis>::scaleOverflows(const Num2* data, const unsigned long dataLength) {
3424     if (dataLength != overflow_.length())
3425       throw npstat::NpstatInvalidArgument("In npstat::HistoND::scaleOverflows: incompatible data length");
3426     assert(data);
3427     Numeric* dat = const_cast<Numeric*>(overflow_.data());
3428     for (unsigned long i = 0; i < dataLength; ++i)
3429       dat[i] *= data[i];
3430     ++modCount_;
3431   }
3432 
3433   template <typename Numeric, class Axis>
3434   template <typename Num2>
3435   inline void HistoND<Numeric, Axis>::setBinContents(const Num2* data,
3436                                                      const unsigned long dataLength,
3437                                                      const bool clearOverflowsNow) {
3438     data_.setData(data, dataLength);
3439     if (clearOverflowsNow)
3440       clearOverflows();
3441     ++modCount_;
3442   }
3443 
3444   template <typename Numeric, class Axis>
3445   template <typename Num2>
3446   inline void HistoND<Numeric, Axis>::setOverflows(const Num2* data, const unsigned long dataLength) {
3447     overflow_.setData(data, dataLength);
3448     ++modCount_;
3449   }
3450 
3451   template <typename Numeric, class Axis>
3452   inline void HistoND<Numeric, Axis>::recalculateNFillsFromData() {
3453     const long double nOver = overflow_.template sum<long double>();
3454     const long double nData = data_.template sum<long double>();
3455     overCount_ = static_cast<unsigned long>(nOver);
3456     fillCount_ = static_cast<unsigned long>(nData + nOver);
3457     ++modCount_;
3458   }
3459 
3460   template <typename Numeric, class Axis>
3461   template <typename Num2, typename Num3>
3462   inline void HistoND<Numeric, Axis>::addToProjection(HistoND<Num2, Axis>* projection,
3463                                                       AbsArrayProjector<Numeric, Num3>& projector,
3464                                                       const unsigned* projectedIndices,
3465                                                       const unsigned nProjectedIndices) const {
3466     assert(projection);
3467     data_.addToProjection(&projection->data_, projector, projectedIndices, nProjectedIndices);
3468     projection->fillCount_ += projection->nBins();
3469     projection->modCount_++;
3470   }
3471 
3472   template <typename Numeric, class Axis>
3473   template <typename Num2, typename Num3>
3474   inline void HistoND<Numeric, Axis>::addToProjection(HistoND<Num2, Axis>* projection,
3475                                                       AbsVisitor<Numeric, Num3>& projector,
3476                                                       const unsigned* projectedIndices,
3477                                                       const unsigned nProjectedIndices) const {
3478     assert(projection);
3479     data_.addToProjection(&projection->data_, projector, projectedIndices, nProjectedIndices);
3480     projection->fillCount_ += projection->nBins();
3481     projection->modCount_++;
3482   }
3483 
3484   template <typename Numeric, class Axis>
3485   const char* HistoND<Numeric, Axis>::classname() {
3486     static const std::string myClass(gs::template_class_name<Numeric, Axis>("npstat::HistoND"));
3487     return myClass.c_str();
3488   }
3489 
3490   template <typename Numeric, class Axis>
3491   bool HistoND<Numeric, Axis>::write(std::ostream& of) const {
3492     gs::write_pod(of, title_);
3493     gs::write_pod(of, accumulatedDataLabel_);
3494     gs::write_pod(of, fillCount_);
3495     gs::write_pod(of, overCount_);
3496 
3497     return !of.fail() && gs::write_obj_vector(of, axes_) && data_.classId().write(of) && data_.write(of) &&
3498            overflow_.write(of);
3499   }
3500 
3501   template <typename Numeric, class Axis>
3502   HistoND<Numeric, Axis>* HistoND<Numeric, Axis>::read(const gs::ClassId& id, std::istream& in) {
3503     static const gs::ClassId current(gs::ClassId::makeId<HistoND<Numeric, Axis> >());
3504     current.ensureSameId(id);
3505 
3506     std::string title;
3507     gs::read_pod(in, &title);
3508 
3509     std::string accumulatedDataLabel;
3510     gs::read_pod(in, &accumulatedDataLabel);
3511 
3512     unsigned long fillCount = 0, overCount = 0;
3513     gs::read_pod(in, &fillCount);
3514     gs::read_pod(in, &overCount);
3515     if (in.fail())
3516       throw gs::IOReadFailure("In npstat::HistoND::read: input stream failure");
3517 
3518     std::vector<Axis> axes;
3519     gs::read_heap_obj_vector_as_placed(in, &axes);
3520     gs::ClassId ida(in, 1);
3521     ArrayND<Numeric> data, over;
3522     ArrayND<Numeric>::restore(ida, in, &data);
3523     ArrayND<Numeric>::restore(ida, in, &over);
3524     CPP11_auto_ptr<HistoND<Numeric, Axis> > result(
3525         new HistoND<Numeric, Axis>(axes, title.c_str(), accumulatedDataLabel.c_str()));
3526     result->data_ = data;
3527     result->overflow_ = over;
3528     result->fillCount_ = fillCount;
3529     result->overCount_ = overCount;
3530     return result.release();
3531   }
3532 
3533   template <typename Histo>
3534   void convertHistoToDensity(Histo* h, const bool knownNonNegative) {
3535     assert(h);
3536     if (!knownNonNegative)
3537       (const_cast<ArrayND<typename Histo::value_type>&>(h->binContents())).makeNonNegative();
3538     const double integ = h->integral();
3539     *h /= integ;
3540   }
3541 
3542   template <typename Histo>
3543   std::vector<LinearMapper1d> densityScanHistoMap(const Histo& histo) {
3544     std::vector<LinearMapper1d> result;
3545     const unsigned d = histo.dim();
3546     result.reserve(d);
3547     for (unsigned i = 0; i < d; ++i) {
3548       const LinearMapper1d& m = histo.axis(i).binNumberMapper(false);
3549       result.push_back(m.inverse());
3550     }
3551     return result;
3552   }
3553 
3554   template <typename Histo>
3555   std::vector<CircularMapper1d> convolutionHistoMap(const Histo& histo, const bool doubleRange) {
3556     std::vector<CircularMapper1d> result;
3557     const unsigned d = histo.dim();
3558     result.reserve(d);
3559     for (unsigned i = 0; i < d; ++i)
3560       result.push_back(histo.axis(i).kernelScanMapper(doubleRange));
3561     return result;
3562   }
3563 }  // namespace npstat
3564 
3565 #endif  // NPSTAT_HISTOND_HH_