hej is hosted by Hepforge, IPPP Durham
HEJ 2.1.4
High energy resummation for hadron colliders
Loading...
Searching...
No Matches
Unweighter_impl.hh
Go to the documentation of this file.
1
8#pragma once
9
10#include <cmath>
11#include <cstdint>
12#include <limits>
13#include <utility>
14
15#include "HEJ/Event.hh"
16#include "HEJ/RNG.hh"
17
18namespace HEJ {
19 namespace detail {
23 class Histogram {
24 public:
25 Histogram(size_t nbins, double min, double max):
26 bins_(nbins),
27 min_{min},
28 width_{(max-min)/nbins}
29 {}
30
31 size_t bin_idx(double pos) {
32 int const idx = std::floor((pos - min_)/width_);
33 if(idx < 0) return 0;
34 return std::min(static_cast<size_t>(idx), bins_.size()-1);
35 }
36
37 void fill(double pos, double val) {
38 size_t const bin = bin_idx(pos);
39 if(bin >= bins_.size()) return;
40 bins_[bin] += val;
41 }
42
43 double operator[](size_t i) const {
44 return bins_[i];
45 }
46 double & operator[](size_t i) {
47 return bins_[i];
48 }
49
50 std::pair<double, double> peak_bin() const {
51 const auto max = std::max_element(begin(bins_), end(bins_));
52 if(max == end(bins_)) {
53 static constexpr double nan = std::numeric_limits<double>::quiet_NaN();
54 return std::make_pair(nan, nan);
55 }
56 const size_t nbin = std::distance(begin(bins_), max);
57 return std::make_pair(min_ + (nbin + 0.5)*width_, *max);
58 }
59
60 private:
61 std::vector<double> bins_;
62 double min_, width_;
63 };
64
66 template<class ConstIt>
67 std::pair<double, double> awt_range(ConstIt begin, ConstIt end) {
68 // can't use std::minmax_element
69 // because boost filter iterators are not assignable
70 assert(begin != end);
71 double min = std::numeric_limits<double>::infinity();
72 double max = 0.;
73 for(auto it = begin; it != end; ++it) {
74 const double awt = std::abs(it->central().weight);
75 if(awt == 0) continue;
76 if(awt < min) {
77 min = awt;
78 } else if (awt > max) {
79 max = awt;
80 }
81 }
82 return std::make_pair(min, max);
83 }
84 } // namespace detail
85
86 template<class ConstIt>
87 void Unweighter::set_cut_to_peakwt(ConstIt begin, ConstIt end, double max_dev){
88 if(begin == end) {
89 throw std::invalid_argument{"Empty range of events"};
90 }
91 double err = 0.;
92 double awt_sum = 0.;
93 const auto extremal_awts = detail::awt_range(begin, end);
94 const double min_lwt = std::log(extremal_awts.first);
95 const double max_lwt = std::log(extremal_awts.second);
96 const size_t nevents = std::distance(begin, end);
97 // heuristic value for number of bins
98 const size_t nbins = std::sqrt(nevents);
99 detail::Histogram wtwt{nbins, min_lwt, max_lwt};
100 detail::Histogram nwt{nbins, min_lwt, max_lwt};
101
102 for(auto it=begin; it != end; ++it){
103 const double awt = std::abs(it->central().weight);
104 nwt.fill(std::log(awt), 1.);
105 }
106
107 const auto cut = nevents/nbins;
108 for(; begin != end; ++begin){
109 const double awt = std::abs(begin->central().weight);
110 const double log_wt = std::log(awt);
111 const double bin_idx = nwt.bin_idx(log_wt);
112 assert(bin_idx<nbins);
113 assert(bin_idx>=0);
114 // prune low statistics bins with a heuristic cut
115 if(nwt[bin_idx] >= cut){
116 wtwt.fill(log_wt, awt);
117 const double tmp = awt*log_wt;
118 err += tmp*tmp;
119 awt_sum += awt;
120 }
121 }
122
123 const auto peak = wtwt.peak_bin();
124 err = std::sqrt(err)/awt_sum;
125 set_cut(std::exp(peak.first + max_dev*err));
126 }
127
128 template<class ConstIt>
129 void Unweighter::set_cut_to_maxwt(ConstIt begin, ConstIt end){
130 set_cut(-1);
131 for(; begin != end; ++begin){
132 const double awt = std::abs(begin->central().weight);
133 if( awt > get_cut())
134 set_cut(awt);
135 }
136 }
137 template<class Iterator>
139 Iterator begin, Iterator end, RNG & ran
140 ) const {
141 if(get_cut() < 0) return end;
142 return std::remove_if(begin, end,
143 [&](auto & ev) -> bool {
144 return this->discard(ran, ev);
145 });
146 }
147} // namespace HEJ
Declares the Event class and helpers.
Interface for pseudorandom number generators.
void set_cut_to_maxwt(std::vector< Event > const &events)
Set cut as max(weight) of events.
Definition: Unweighter.hh:35
void set_cut(double max_weight)
Explicitly set cut.
Definition: Unweighter.hh:31
void set_cut_to_peakwt(std::vector< Event > const &events, double max_dev)
Estimate some reasonable cut for partial unweighting.
Definition: Unweighter.hh:47
optional< Event > unweight(Event ev, RNG &ran) const
Unweight one event, returns original event if weight > get_cut()
double get_cut() const
Returns current value of the cut.
Definition: Unweighter.hh:55
Definition: Unweighter_impl.hh:23
Histogram(size_t nbins, double min, double max)
Definition: Unweighter_impl.hh:25
void fill(double pos, double val)
Definition: Unweighter_impl.hh:37
double & operator[](size_t i)
Definition: Unweighter_impl.hh:46
size_t bin_idx(double pos)
Definition: Unweighter_impl.hh:31
double operator[](size_t i) const
Definition: Unweighter_impl.hh:43
std::pair< double, double > peak_bin() const
Definition: Unweighter_impl.hh:50
std::pair< double, double > awt_range(ConstIt begin, ConstIt end)
Get minimal and maximal absolute weight.
Definition: Unweighter_impl.hh:67
Main HEJ 2 Namespace.
Definition: mainpage.dox:1
@ discard
Definition: StatusCode.hh:18
Interface for random number generator.
Definition: RNG.hh:19