Rivet/Tools/HistoGroup.hh

0001 // -*- C++ -*-
0002 #ifndef RIVET_HISTOGROUP_HH
0003 #define RIVET_HISTOGROUP_HH
0004
0005 #include "Rivet/Config/RivetCommon.hh"
0006 #include "Rivet/Tools/Logging.hh"
0007 #include "Rivet/Tools/RivetYODA.hh"
0008
0009 namespace Rivet {
0010
0011   namespace {
0012
0013     /// @brief Dummy fill adapter for the HistoGroup
0014     ///
0015     /// The group bin is already being determined in the fill method.
0016     template <size_t FillDim, typename T, typename AxisT>
0017     typename YODA::FillableStorage<FillDim, T, AxisT>::FillAdapterT
0018     groupAdapter = [](auto& /* ptr */, auto&& /* coords */, double /* weight */, double /* fraction */) { };
0019
0020   }
0021
0022   template <typename GroupAxisT, typename... AxisT>
0023   class HistoGroup : public YODA::FillableStorage<sizeof...(AxisT)+1, BinnedHistoPtr<AxisT...>, GroupAxisT>,
0024                      public YODA::Fillable {
0025   public:
0026
0027     using BaseT = YODA::FillableStorage<sizeof...(AxisT)+1, BinnedHistoPtr<AxisT...>, GroupAxisT>;
0028     using BinContentT = BinnedHistoPtr<AxisT...>;
0029     using FillDim = std::integral_constant<size_t, sizeof...(AxisT)+1>;
0030
0031     HistoGroup() : BaseT(groupAdapter<FillDim{}, BinContentT, GroupAxisT>) { }
0032
0033     HistoGroup(const std::vector<GroupAxisT>& edges)
0034          : BaseT(YODA::Axis<GroupAxisT>(edges), groupAdapter<FillDim{}, BinContentT, GroupAxisT>) { }
0035
0036     HistoGroup(std::initializer_list<GroupAxisT>&& edges)
0037          : BaseT(YODA::Axis<GroupAxisT>(std::move(edges)), groupAdapter<FillDim{}, BinContentT, GroupAxisT>) { }
0038
0039     int fill(GroupAxisT tmpCoord, AxisT... coords, const double weight = 1.0, const double fraction = 1.0) {
0040       auto& bin = BaseT::binAt(tmpCoord);
0041       if (!bin.raw())  return -1; // nullptr if bin not booked
0042       return bin->fill({coords...}, weight, fraction);
0043     }
0044
0045     /// @name Reset methods
0046     /// @{
0047
0048     /// @brief Reset the histogram.
0049     ///
0050     /// Keep the binning but set all bin contents and related quantities to zero
0051     void reset() noexcept { BaseT::reset(); }
0052
0053     /// @}
0054
0055     /// @name Whole histo data
0056     /// @{
0057
0058     /// @brief Fill dimension of the object (number of conent axes + temprary axis)
0059     size_t fillDim() const noexcept { return FillDim{}; }
0060
0061     /// @brief Total dimension of the object (number of fill axes + filled value)
0062     size_t dim() const noexcept { return fillDim()+1; }
0063
0064     /// @brief Get the number of fills (fractional fills are possible).
0065     double numEntries(const bool includeOverflows=true) const noexcept {
0066       double n = 0;
0067       for (const auto& b : BaseT::bins(includeOverflows)) {
0068         if (!b.get())  continue;
0069         n += b->numEntries(includeOverflows);
0070       }
0071       return n;
0072     }
0073
0074     /// @brief Get the effective number of fills.
0075     double effNumEntries(const bool includeOverflows=true) const noexcept {
0076       double n = 0;
0077       for (const auto& b : BaseT::bins(includeOverflows)) {
0078         if (!b.get())  continue;
0079         n += b->effNumEntries(includeOverflows);
0080       }
0081       return n;
0082     }
0083
0084     /// @brief Calculates sum of weights in histo group.
0085     double sumW(const bool includeOverflows=true) const noexcept {
0086       double sumw = 0;
0087       for (const auto& b : BaseT::bins(includeOverflows)) {
0088         if (!b.get())  continue;
0089         sumw += b->sumW(includeOverflows);
0090       }
0091       return sumw;
0092     }
0093
0094     /// @brief Calculates sum of squared weights in histo group.
0095     double sumW2(const bool includeOverflows=true) const noexcept {
0096       double sumw2 = 0;
0097       for (const auto& b : BaseT::bins(includeOverflows)) {
0098         if (!b.get())  continue;
0099         sumw2 += b->sumW2(includeOverflows);
0100       }
0101       return sumw2;
0102     }
0103
0104     /// @brief Return the vector of sum of weights for each histo in the group.
0105     vector<double> sumWGroup(const bool includeOverflows=true) const noexcept {
0106       vector<double> rtn;
0107       rtn.reserve(BaseT::numBins(true));
0108       for (const auto& b : BaseT::bins(true)) {
0109         if (!b.get()) {
0110           rtn.push_back(0.0);
0111           continue;
0112         }
0113         rtn.push_back( b->sumW(includeOverflows) );
0114       }
0115       return rtn;
0116     }
0117
0118     /// @brief Get the total volume of the histogram group.
0119     double integral(const bool includeOverflows=true) const noexcept {
0120       return sumW(includeOverflows);
0121     }
0122
0123     /// @brief Get the total volume error of the histogram group.
0124     double integralError(const bool includeOverflows=true) const noexcept {
0125       return sqrt(sumW2(includeOverflows));
0126     }
0127
0128     /// @}
0129
0130     /// @name Transformations
0131     /// @{
0132
0133     void divByGroupWidth() const noexcept {
0134       for (auto& bin : BaseT::bins(true, true)) {
0135         if (!bin.get())  continue;
0136         const double bw = bin.dVol();
0137         if (bw)  bin->scaleW(1.0/bw);
0138       }
0139     }
0140
0141     /// @brief Rescale as if all fill weights had been different by factor @a scalefactor.
0142     void scaleW(const double scalefactor) noexcept {
0143       for (auto& bin : BaseT::bins(true, true)) {
0144         if (!bin.get())  continue;
0145         bin->scaleW(scalefactor);
0146       }
0147     }
0148
0149     /// @brief Rescale as if all fill weights had been different by factor @a scalefactor along dimension @a i.
0150     void scale(const size_t i, const double scalefactor) noexcept {
0151       for (auto& bin : BaseT::bins(true, true)) {
0152         if (!bin.get())  continue;
0153         bin->scale(i, scalefactor);
0154       }
0155     }
0156
0157     /// @brief Normalize the (visible) histo "volume" to the @a normto value.
0158     ///
0159     /// @note Consider each sub-histogram separately when calculating the integral.
0160     void normalize(const double normto=1.0, const bool includeOverflows=true) {
0161       for (auto& bin : BaseT::bins(true, true)) {
0162         if (!bin.get())  continue;
0163         if (bin->integral(includeOverflows) != 0.)  bin->normalize(normto, includeOverflows);
0164       }
0165     }
0166
0167     /// @brief Normalize the (visible) histo "volume" to the @a normto value.
0168     ///
0169     /// @note Treat the histogram group as an entity when calculating the integral.
0170     void normalizeGroup(const double normto=1.0, const bool includeOverflows=true) {
0171       const double oldintegral = integral(includeOverflows);
0172       if (oldintegral == 0) {
0173         MSG_DEBUG("Attempted to normalize a histogram group with null area; skipping.");
0174         return;
0175       }
0176       scaleW(normto / oldintegral);
0177     }
0178
0179     /// @}
0180
0181     protected:
0182
0183     /// Get a Log object based on the name() property of the calling analysis object.
0184     Log& getLog() const {
0185       string logname = "Rivet.HistoGroup";
0186       return Log::getLog(logname);
0187     }
0188
0189   };
0190
0191   /// @name Convenient aliases and dimension-specific short-hands
0192   /// @{
0193
0194   template <typename GroupAxisT, typename... AxisT>
0195   using HistoGroupPtr = std::shared_ptr<HistoGroup<GroupAxisT, AxisT...>>;
0196   using Histo1DGroupPtr = HistoGroupPtr<double,double>;
0197   using Histo2DGroupPtr = HistoGroupPtr<double,double,double>;
0198
0199   /// @}
0200
0201 }
0202
0203 #endif