Analysis/LWH/Histogram1D.h

0001 // -*- C++ -*-
0002 //
0003 // Histogram1D.h is a part of ThePEG - Toolkit for HEP Event Generation
0004 // Copyright (C) 1999-2019 Leif Lonnblad
0005 //
0006 // ThePEG is licenced under version 3 of the GPL, see COPYING for details.
0007 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
0008 //
0009 #ifndef LWH_Histogram1D_H
0010 #define LWH_Histogram1D_H
0011 //
0012 // This is the declaration of the Histogram1D class.
0013 //
0014
0015 #include "AIHistogram1D.h"
0016 #include "ManagedObject.h"
0017 #include "Axis.h"
0018 #include "VariAxis.h"
0019 #include <vector>
0020 #include <stdexcept>
0021
0022 namespace LWH {
0023
0024 using namespace AIDA;
0025
0026 /**
0027  * User level interface to 1D Histogram.
0028  */
0029 class Histogram1D: public IHistogram1D, public ManagedObject {
0030
0031 public:
0032
0033   /** HistFactory is a friend. */
0034   friend class HistogramFactory;
0035
0036 public:
0037
0038   /**
0039    * Standard constructor.
0040    */
0041   Histogram1D(int n, double lo, double up)
0042     : fax(new Axis(n, lo, up)), vax(0),
0043       sum(n + 2), sumw(n + 2), sumw2(n + 2), sumxw(n + 2), sumx2w(n + 2) {
0044     ax = fax;
0045   }
0046
0047   /**
0048    * Standard constructor for variable bin width.
0049    */
0050   Histogram1D(const std::vector<double> & edges)
0051     : fax(0), vax(new VariAxis(edges)),
0052       sum(edges.size() + 1), sumw(edges.size() + 1), sumw2(edges.size() + 1),
0053       sumxw(edges.size() + 1), sumx2w(edges.size() + 1) {
0054     ax = vax;
0055   }
0056
0057   /**
0058    * Copy constructor.
0059    */
0060   Histogram1D(const Histogram1D & h)
0061     : IBaseHistogram(h), IHistogram(h), IHistogram1D(h), ManagedObject(h),
0062       fax(0), vax(0), sum(h.sum), sumw(h.sumw), sumw2(h.sumw2),
0063       sumxw(h.sumxw), sumx2w(h.sumx2w) {
0064     const VariAxis * hvax = dynamic_cast<const VariAxis *>(h.ax);
0065     if ( vax ) ax = vax = new VariAxis(*hvax);
0066     else ax = fax = new Axis(dynamic_cast<const Axis &>(*h.ax));
0067 }
0068
0069   /// Destructor.
0070   virtual ~Histogram1D() {
0071     delete ax;
0072   }
0073
0074   /**
0075    * Get the Histogram's title.
0076    * @return The Histogram's title.
0077    */
0078   std::string title() const {
0079     return theTitle;
0080   }
0081
0082   /**
0083    * Get the Histogram's title.
0084    * @return The Histogram's title.
0085    */
0086   std::string name() const {
0087     return theTitle;
0088   }
0089
0090   /**
0091    * Set the histogram title.
0092    * @param title The title.
0093    * @return false If title cannot be changed.
0094    */
0095   bool setTitle(const std::string & title) {
0096     theTitle = title;
0097     return true;
0098   }
0099
0100   /**
0101    * Not implemented in LWH. will throw an exception.
0102    */
0103   IAnnotation & annotation() {
0104     throw std::runtime_error("LWH cannot handle annotations");
0105   }
0106
0107   /**
0108    * Not implemented in LWH. will throw an exception.
0109    */
0110   const IAnnotation & annotation() const {
0111     throw std::runtime_error("LWH cannot handle annotations");
0112   }
0113
0114   /**
0115    * Get the Histogram's dimension.
0116    * @return The Histogram's dimension.
0117    */
0118   int dimension() const {
0119     return 1;
0120   }
0121
0122   /**
0123    * Reset the Histogram; as if just created.
0124    * @return false If something goes wrong.
0125    */
0126   bool reset() {
0127     sum = std::vector<int>(ax->bins() + 2);
0128     sumw = std::vector<double>(ax->bins() + 2);
0129     sumxw = std::vector<double>(ax->bins() + 2);
0130     sumx2w = std::vector<double>(ax->bins() + 2);
0131     sumw2 = std::vector<double>(ax->bins() + 2);
0132     return true;
0133   }
0134
0135   /**
0136    * Get the number of in-range entries in the Histogram.
0137    * @return The number of in-range entries.
0138    *
0139    */
0140   int entries() const {
0141     int si = 0;
0142     for ( int i = 2; i < ax->bins() + 2; ++i ) si += sum[i];
0143     return si;
0144   }
0145
0146   /**
0147    * Sum of the entries in all the IHistogram's bins,
0148    * i.e in-range bins, UNDERFLOW and OVERFLOW.
0149    * This is equivalent to the number of times the
0150    * method fill was invoked.
0151    * @return The sum of all the entries.
0152    */
0153   int allEntries() const {
0154     return entries() + extraEntries();
0155   }
0156
0157   /**
0158    * Number of entries in the UNDERFLOW and OVERFLOW bins.
0159    * @return The number of entries outside the range of the IHistogram.
0160    */
0161   int extraEntries() const {
0162     return sum[0] + sum[1];
0163   }
0164
0165   /**
0166    * Number of equivalent entries,
0167    * i.e. <tt>SUM[ weight ] ^ 2 / SUM[ weight^2 ]</tt>
0168    * @return The number of equivalent entries.
0169    */
0170   double equivalentBinEntries() const {
0171     double sw = 0.0;
0172     double sw2 = 0.0;
0173     for ( int i = 2; i < ax->bins() + 2; ++i ) {
0174       sw += sumw[i];
0175       sw2 += sumw2[i];
0176     }
0177     return sw2/(sw*sw);
0178   }
0179
0180   /**
0181    * Sum of in-range bin heights in the IHistogram,
0182    * UNDERFLOW and OVERFLOW bins are excluded.
0183    * @return The sum of the in-range bins heights.
0184    *
0185    */
0186   double sumBinHeights() const {
0187     double sw = 0.0;
0188     for ( int i = 2; i < ax->bins() + 2; ++i ) sw += sumw[i];
0189     return sw;
0190   }
0191
0192   /**
0193    * Sum of the heights of all the IHistogram's bins,
0194    * i.e in-range bins, UNDERFLOW and OVERFLOW.
0195    * @return The sum of all the bins heights.
0196    */
0197   double sumAllBinHeights() const {
0198     return sumBinHeights() + sumExtraBinHeights();
0199   }
0200
0201   /**
0202    * Sum of heights in the UNDERFLOW and OVERFLOW bins.
0203    * @return The sum of the heights of the out-of-range bins.
0204    */
0205   double sumExtraBinHeights() const {
0206     return sumw[0] + sumw[1];
0207   }
0208
0209   /**
0210    * Minimum height of the in-range bins,
0211    * i.e. not considering the UNDERFLOW and OVERFLOW bins.
0212    * @return The minimum height among the in-range bins.
0213    */
0214   double minBinHeight() const {
0215     double minw = sumw[2];
0216     for ( int i = 3; i < ax->bins() + 2; ++i ) minw = std::min(minw, sumw[i]);
0217     return minw;
0218   }
0219
0220   /**
0221    * Maximum height of the in-range bins,
0222    * i.e. not considering the UNDERFLOW and OVERFLOW bins.
0223    * @return The maximum height among the in-range bins.
0224    */
0225   double maxBinHeight() const{
0226     double maxw = sumw[2];
0227     for ( int i = 3; i < ax->bins() + 2; ++i ) maxw = std::max(maxw, sumw[i]);
0228     return maxw;
0229   }
0230
0231   /**
0232    * Fill the IHistogram1D with a value and the
0233    * corresponding weight.
0234    * @param x      The value to be filled in.
0235    * @param weight The corresponding weight (by default 1).
0236    * @return false If the weight is <0 or >1 (?).
0237    */
0238   bool fill(double x, double weight = 1.) {
0239     int i = ax->coordToIndex(x) + 2;
0240     ++sum[i];
0241     sumw[i] += weight;
0242     sumxw[i] += x*weight;
0243     sumx2w[i] += x*x*weight;
0244     sumw2[i] += weight*weight;
0245     return weight >= 0 && weight <= 1;
0246   }
0247
0248   /**
0249    * The weighted mean of a bin.
0250    * @param index The bin number (0...N-1) or OVERFLOW or UNDERFLOW.
0251    * @return      The mean of the corresponding bin.
0252    */
0253   double binMean(int index) const {
0254     int i = index + 2;
0255     return sumw[i] != 0.0? sumxw[i]/sumw[i]:
0256       ( vax? vax->binMidPoint(index): fax->binMidPoint(index) );
0257   };
0258
0259   /**
0260    * The weighted RMS of a bin.
0261    * @param index The bin number (0...N-1) or OVERFLOW or UNDERFLOW.
0262    * @return      The RMS of the corresponding bin.
0263    */
0264   double binRms(int index) const {
0265     int i = index + 2;
0266     return sumw[i] == 0.0 || sum[i] < 2? ax->binWidth(index):
0267       std::sqrt(std::max(sumw[i]*sumx2w[i] - sumxw[i]*sumxw[i], 0.0))/sumw[i];
0268   };
0269
0270   /**
0271    * Number of entries in the corresponding bin (ie the number of
0272    * times fill was called for this bin).
0273    * @param index The bin number (0...N-1) or OVERFLOW or UNDERFLOW.
0274    * @return      The number of entries in the corresponding bin.
0275    */
0276   int binEntries(int index) const {
0277     return sum[index + 2];
0278   }
0279
0280   /**
0281    * Total height of the corresponding bin (ie the sum of the weights
0282    * in this bin).
0283    * @param index The bin number (0...N-1) or OVERFLOW or UNDERFLOW.
0284    * @return      The height of the corresponding bin.
0285    */
0286   double binHeight(int index) const {
0287     return sumw[index + 2];
0288   }
0289
0290   /**
0291    * The error of a given bin.
0292    * @param index The bin number (0...N-1) or OVERFLOW or UNDERFLOW.
0293    * @return      The error on the corresponding bin.
0294    *
0295    */
0296   double binError(int index) const {
0297     return std::sqrt(sumw2[index + 2]);
0298   }
0299
0300   /**
0301    * The mean of the whole IHistogram1D.
0302    * @return The mean of the IHistogram1D.
0303    */
0304   double mean() const {
0305     double s = 0.0;
0306     double sx = 0.0;
0307     for ( int i = 2; i < ax->bins() + 2; ++i ) {
0308       s += sumw[i];
0309       sx += sumxw[i];
0310     }
0311     return s != 0.0? sx/s: 0.0;
0312   }
0313
0314   /**
0315    * The RMS of the whole IHistogram1D.
0316    * @return The RMS if the IHistogram1D.
0317    */
0318   double rms() const {
0319     double s = 0.0;
0320     double sx = 0.0;
0321     double sx2 = 0.0;
0322     for ( int i = 2; i < ax->bins() + 2; ++i ) {
0323       s += sumw[i];
0324       sx += sumxw[i];
0325       sx2 += sumx2w[i];
0326     }
0327     return s != 0.0? std::sqrt(std::max(s*sx2 - sx*sx, 0.0))/s:
0328       ax->upperEdge() - ax->lowerEdge();
0329   }
0330
0331   /**
0332    * Get the x axis of the IHistogram1D.
0333    * @return The x coordinate IAxis.
0334    */
0335   const IAxis & axis() const {
0336     return *ax;
0337   }
0338
0339   /**
0340    * Get the bin number corresponding to a given coordinate along the
0341    * x axis.  This is a convenience method, equivalent to
0342    * <tt>axis().coordToIndex(coord)</tt>.
0343    * @param coord The coordinalte along the x axis.
0344    * @return      The corresponding bin number.
0345    */
0346   int coordToIndex(double coord) const {
0347     return ax->coordToIndex(coord);
0348   }
0349
0350   /**
0351    * Add to this Histogram1D the contents of another IHistogram1D.
0352    * @param h The Histogram1D to be added to this IHistogram1D.
0353    * @return false If the IHistogram1Ds binnings are incompatible.
0354    */
0355   bool add(const Histogram1D & h) {
0356     if ( ax->upperEdge() != h.ax->upperEdge() ||
0357      ax->lowerEdge() != h.ax->lowerEdge() ||
0358      ax->bins() != h.ax->bins() ) return false;
0359     for ( int i = 0; i < ax->bins() + 2; ++i ) {
0360       sum[i] += h.sum[i];
0361       sumw[i] += h.sumw[i];
0362       sumxw[i] += h.sumxw[i];
0363       sumx2w[i] += h.sumx2w[i];
0364       sumw2[i] += h.sumw2[i];
0365     }
0366     return true;
0367   }
0368
0369   /**
0370    * Add to this IHistogram1D the contents of another IHistogram1D.
0371    * @param hist The IHistogram1D to be added to this IHistogram1D.
0372    * @return false If the IHistogram1Ds binnings are incompatible.
0373    */
0374   bool add(const IHistogram1D & hist) {
0375     return add(dynamic_cast<const Histogram1D &>(hist));
0376   }
0377
0378   /**
0379    * Scale the contents of this histogram with the given factor.
0380    * @param s the scaling factor to use.
0381    */
0382   bool scale(double s) {
0383     for ( int i = 0; i < ax->bins() + 2; ++i ) {
0384       sumw[i] *= s;
0385       sumxw[i] *= s;
0386       sumx2w[i] *= s;
0387       sumw2[i] *= s*s;
0388     }
0389     return true;
0390   }
0391
0392   /**
0393    * Scale the given histogram so that the integral over all bins
0394    * (including overflow) gives \a intg. This function also corrects
0395    * for the bin-widths, which means that it should only be run once
0396    * for each histogram. Further rescaling must be done with the
0397    * scale(double) function.
0398    */
0399   void normalize(double intg) {
0400     double oldintg = sumAllBinHeights();
0401     if ( oldintg == 0.0 ) return;
0402     for ( int i = 0; i < ax->bins() + 2; ++i ) {
0403       double fac = intg/oldintg;
0404       if ( i >= 2 ) fac /= (ax->binUpperEdge(i - 2) - ax->binLowerEdge(i - 2));
0405       sumw[i] *= fac;
0406       sumxw[i] *= fac;
0407       sumx2w[i] *= fac;
0408       sumw2[i] *= fac*fac;
0409     }
0410   }
0411
0412   /**
0413    * Return the integral over the histogram bins assuming it has been
0414    * normalize()d.
0415    */
0416   double integral() const {
0417     double intg = sumw[0] + sumw[1];
0418     for ( int i = 2; i < ax->bins() + 2; ++i )
0419       intg += sumw[i]*(ax->binUpperEdge(i - 2) - ax->binLowerEdge(i - 2));
0420     return intg;
0421   }
0422
0423   /**
0424    * Not implemented in LWH.
0425    * @return null pointer always.
0426    */
0427   void * cast(const std::string &) const {
0428     return 0;
0429   }
0430
0431   /**
0432    * Write out the histogram in the AIDA xml format.
0433    */
0434   bool writeXML(std::ostream & os, std::string path, std::string name) {
0435     os << "  <histogram1d name=\"" << name
0436        << "\"\n    title=\"" << title()
0437        << "\" path=\"" << path
0438        << "\">\n    <axis max=\"" << ax->upperEdge()
0439        << "\" numberOfBins=\"" << ax->bins()
0440        << "\" min=\"" << ax->lowerEdge()
0441        << "\" direction=\"x\"";
0442     if ( vax ) {
0443       os << ">\n";
0444       for ( int i = 0, N = ax->bins() - 1; i < N; ++i )
0445     os << "      <binBorder value=\"" << ax->binUpperEdge(i) << "\"/>\n";
0446       os << "    </axis>\n";
0447     } else {
0448       os << "/>\n";
0449     }
0450     os << "    <statistics entries=\"" << entries()
0451        << "\">\n      <statistic mean=\"" << mean()
0452        << "\" direction=\"x\"\n        rms=\"" << rms()
0453        << "\"/>\n    </statistics>\n    <data1d>\n";
0454     for ( int i = 0; i < ax->bins() + 2; ++i ) if ( sum[i] ) {
0455       os << "      <bin1d binNum=\"";
0456       if ( i == 0 ) os << "UNDERFLOW";
0457       else if ( i == 1 ) os << "OVERFLOW";
0458       else os << i - 2;
0459       os << "\" entries=\"" << sum[i]
0460      << "\" height=\"" << sumw[i]
0461      << "\"\n        error=\"" << std::sqrt(sumw2[i])
0462      << "\" error2=\"" << sumw2[i]
0463      << "\"\n        weightedMean=\"" << binMean(i - 2)
0464      << "\" weightedRms=\"" << binRms(i - 2)
0465      << "\"/>\n";
0466     }
0467     os << "    </data1d>\n  </histogram1d>" << std::endl;
0468     return true;
0469   }
0470
0471
0472   /**
0473    * Write out the histogram in a flat text file suitable for
0474    * eg. gnuplot to read. The coloums are layed out as 'x w w2 n'.
0475    */
0476   bool writeFLAT(std::ostream & os, std::string path, std::string name) {
0477     os << "# " << path << "/" << name << " " << ax->lowerEdge()
0478        << " " << ax->bins() << " " << ax->upperEdge()
0479        << " \"" << title() << " \"" << std::endl;
0480     for ( int i = 2; i < ax->bins() + 2; ++i )
0481       os << 0.5*(ax->binLowerEdge(i - 2) + ax->binUpperEdge(i - 2)) << " "
0482      << sumw[i] << " " << sqrt(sumw2[i]) << " " << sum[i] << std::endl;
0483     os << std::endl;
0484     return true;
0485   }
0486
0487 private:
0488
0489   /** The title */
0490   std::string theTitle;
0491
0492   /** The axis. */
0493   IAxis * ax;
0494
0495   /** Pointer (possibly null) to a axis with fixed bin width. */
0496   Axis * fax;
0497
0498   /** Pointer (possibly null) to a axis with fixed bin width. */
0499   VariAxis * vax;
0500
0501   /** The counts. */
0502   std::vector<int> sum;
0503
0504   /** The weights. */
0505   std::vector<double> sumw;
0506
0507   /** The squared weights. */
0508   std::vector<double> sumw2;
0509
0510   /** The weighted x-values. */
0511   std::vector<double> sumxw;
0512
0513   /** The weighted x-square-values. */
0514   std::vector<double> sumx2w;
0515
0516 };
0517
0518 }
0519
0520 #endif /* LWH_Histogram1D_H */