hej
is hosted by
Hepforge
,
IPPP Durham
HEJ
2.1.4
High energy resummation for hadron colliders
Loading...
Searching...
No Matches
include
HEJ
detail
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
18
namespace
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>
138
Iterator
Unweighter::unweight
(
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
Event.hh
Declares the Event class and helpers.
RNG.hh
Interface for pseudorandom number generators.
HEJ::Unweighter::set_cut_to_maxwt
void set_cut_to_maxwt(std::vector< Event > const &events)
Set cut as max(weight) of events.
Definition:
Unweighter.hh:35
HEJ::Unweighter::set_cut
void set_cut(double max_weight)
Explicitly set cut.
Definition:
Unweighter.hh:31
HEJ::Unweighter::set_cut_to_peakwt
void set_cut_to_peakwt(std::vector< Event > const &events, double max_dev)
Estimate some reasonable cut for partial unweighting.
Definition:
Unweighter.hh:47
HEJ::Unweighter::unweight
optional< Event > unweight(Event ev, RNG &ran) const
Unweight one event, returns original event if weight > get_cut()
HEJ::Unweighter::get_cut
double get_cut() const
Returns current value of the cut.
Definition:
Unweighter.hh:55
HEJ::detail::Histogram
Definition:
Unweighter_impl.hh:23
HEJ::detail::Histogram::Histogram
Histogram(size_t nbins, double min, double max)
Definition:
Unweighter_impl.hh:25
HEJ::detail::Histogram::fill
void fill(double pos, double val)
Definition:
Unweighter_impl.hh:37
HEJ::detail::Histogram::operator[]
double & operator[](size_t i)
Definition:
Unweighter_impl.hh:46
HEJ::detail::Histogram::bin_idx
size_t bin_idx(double pos)
Definition:
Unweighter_impl.hh:31
HEJ::detail::Histogram::operator[]
double operator[](size_t i) const
Definition:
Unweighter_impl.hh:43
HEJ::detail::Histogram::peak_bin
std::pair< double, double > peak_bin() const
Definition:
Unweighter_impl.hh:50
HEJ::detail::awt_range
std::pair< double, double > awt_range(ConstIt begin, ConstIt end)
Get minimal and maximal absolute weight.
Definition:
Unweighter_impl.hh:67
HEJ
Main HEJ 2 Namespace.
Definition:
mainpage.dox:1
HEJ::discard
@ discard
Definition:
StatusCode.hh:18
HEJ::RNG
Interface for random number generator.
Definition:
RNG.hh:19
Generated by
1.9.5