Page MenuHomeHEPForge
Files F7879809

No OneTemporary
Actions

View Options

diff --git a/include/HEJ/detail/Unweighter_impl.hh b/include/HEJ/detail/Unweighter_impl.hh
index feeb97c..a8a44eb 100644
--- a/include/HEJ/detail/Unweighter_impl.hh
+++ b/include/HEJ/detail/Unweighter_impl.hh
@@ -1,147 +1,147 @@
/** \file
* \brief Unweighter Class Implementation for template functions
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
#include "HEJ/Unweighter.hh"
#include <limits>
#include <utility>
#include <cstdint>
#include "HEJ/Event.hh"
#include "HEJ/RNG.hh"
namespace HEJ {
namespace detail {
// primitive histogram class
// we don't want to do much with this,
// so try to avoid an external dependence like YODA or ROOT
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) {
- // needs intermediate cast since negative double -> unsigned
- // is *undefined behaviour*
- return static_cast<int>((pos - min_)/width_);
+ int const idx = std::abs((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);
}
}
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);
- // Set range minimaly higher to include min/max inside range
- const double min_lwt = 0.99999*std::log(extremal_awts.first);
- const double max_lwt = 1.00001*std::log(extremal_awts.second);
+ 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

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 9:07 PM (1 d, 20 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3806201
Default Alt Text
(4 KB)

Event Timeline