Unweighter_impl.hh

Go to the documentation of this file.

 
#pragma once
 
#include <cmath>
#include <cstdint>
#include <limits>
#include <utility>
 
#include "HEJ/Event.hh"
#include "HEJ/RNG.hh"
 
namespace HEJ {
  namespace detail {
    class Histogram {
    public:
      Histogram(size_t nbins, double min, double max):
        bins_(nbins),
        min_{min},
        width_{(max-min)/nbins}
      {}
 
      size_t bin_idx(double pos) {
        int const idx = std::floor((pos - min_)/width_);
        if(idx < 0) return 0;
        return std::min(static_cast<size_t>(idx), bins_.size()-1);
      }
 
      void fill(double pos, double val) {
        size_t const bin = bin_idx(pos);
        if(bin >= bins_.size()) return;
        bins_[bin] += val;
      }
 
      double operator[](size_t i) const {
        return bins_[i];
      }
      double & operator[](size_t i) {
        return bins_[i];
      }
 
      std::pair<double, double> peak_bin() const {
        const auto max = std::max_element(begin(bins_), end(bins_));
        if(max == end(bins_)) {
          static constexpr double nan = std::numeric_limits<double>::quiet_NaN();
          return std::make_pair(nan, nan);
        }
        const size_t nbin = std::distance(begin(bins_), max);
        return std::make_pair(min_ + (nbin + 0.5)*width_, *max);
      }
 
    private:
      std::vector<double> bins_;
      double min_, width_;
    };
 
    template<class ConstIt>
    std::pair<double, double> awt_range(ConstIt begin, ConstIt end) {
      // can't use std::minmax_element
      // because boost filter iterators are not assignable
      assert(begin != end);
      double min = std::numeric_limits<double>::infinity();
      double max = 0.;
      for(auto it = begin; it != end; ++it) {
        const double awt = std::abs(it->central().weight);
        if(awt == 0) continue;
        if(awt < min) {
          min = awt;
        } else if (awt > max) {
          max = awt;
        }
      }
      return std::make_pair(min, max);
    }
  } // namespace detail
 
  template<class ConstIt>
  void Unweighter::set_cut_to_peakwt(ConstIt begin, ConstIt end, double max_dev){
    if(begin == end) {
      throw std::invalid_argument{"Empty range of events"};
    }
    double err = 0.;
    double awt_sum = 0.;
    const auto extremal_awts = detail::awt_range(begin, end);
    const double min_lwt = std::log(extremal_awts.first);
    const double max_lwt = std::log(extremal_awts.second);
    const size_t nevents = std::distance(begin, end);
    // heuristic value for number of bins
    const size_t nbins = std::sqrt(nevents);
    detail::Histogram wtwt{nbins, min_lwt, max_lwt};
    detail::Histogram nwt{nbins, min_lwt, max_lwt};
 
    for(auto it=begin; it != end; ++it){
      const double awt = std::abs(it->central().weight);
      nwt.fill(std::log(awt), 1.);
    }
 
    const auto cut = nevents/nbins;
    for(; begin != end; ++begin){
      const double awt = std::abs(begin->central().weight);
      const double log_wt = std::log(awt);
      const double bin_idx = nwt.bin_idx(log_wt);
      assert(bin_idx<nbins);
      assert(bin_idx>=0);
      // prune low statistics bins with a heuristic cut
      if(nwt[bin_idx] >= cut){
        wtwt.fill(log_wt, awt);
        const double tmp = awt*log_wt;
        err += tmp*tmp;
        awt_sum += awt;
      }
    }
 
    const auto peak = wtwt.peak_bin();
    err = std::sqrt(err)/awt_sum;
    set_cut(std::exp(peak.first + max_dev*err));
  }
 
  template<class ConstIt>
  void Unweighter::set_cut_to_maxwt(ConstIt begin, ConstIt end){
    set_cut(-1);
    for(; begin != end; ++begin){
      const double awt = std::abs(begin->central().weight);
      if( awt > get_cut())
        set_cut(awt);
    }
  }
  template<class Iterator>
  Iterator Unweighter::unweight(
    Iterator begin, Iterator end, RNG & ran
  ) const {
    if(get_cut() < 0) return end;
    return std::remove_if(begin, end,
                [&](auto & ev) -> bool {
                    return this->discard(ran, ev);
                });
  }
} // namespace HEJ