1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
|
/*********************************************
* A simple driver for reweighting management
* S. Carrazza - June 2017
********************************************/
#include <iostream>
#include <string>
#include <map>
#include <vector>
#include <stdexcept>
#include <cmath>
#include <fstream>
#include <iomanip>
using std::string;
using std::vector;
using std::map;
using std::cout;
using std::endl;
using std::scientific;
using std::setw;
// export fortran methods
extern "C" void bookupeqbins_(const char *, const double&, const double&, const double&, int);
extern "C" void filld_(const char*, const double&, const double*, int);
/**
* @brief The reweight class
* This class performs the management of the reweighting
* output and elaboration.
*/
class reweight
{
public:
/**
* @brief Allocate class and create file.
*/
reweight()
{
_output.open("reweight.dat");
}
/**
* Closes the output file
*/
~reweight()
{
_output.close();
}
/**
* Books analysis histograms (alla fortran)
* and stores histogram bins.
*/
void book(string const& key, double const& step,
double const& min, double const& max)
{
// general argument tests
if (min >= max)
throw std::runtime_error("reweight::book: min is smaller than max for key " + key);
if (std::floor( (max-min)/step ) != (max-min)/step)
throw std::runtime_error("reweight::book: step is not compatible with min/max for key " + key);
// book histograms as usual
bookupeqbins_(key.c_str(), step, min, max, key.size());
// populate dict
if (_bins.find(key) == _bins.end())
{
const int N = (max-min)/step;
_bins[key] = linspace(min, max, N);
_output << scientific << key << "\t" << min << "\t" << max << "\t" << N << endl;
}
else
throw std::runtime_error("reweight::book: key " + key + "already exists.");
}
/**
* Take weights and dump to file.
*/
void fill(string const& key, double const& obs, const double* weight)
{
// fill the fortran histogram
filld_(key.c_str(), obs, weight, key.size());
_binid[key] = std::lower_bound(_bins[key].begin(), _bins[key].end(), obs)-_bins[key].begin();
}
/**
* Store event.
*/
void store(const double* weight, double const& l)
{
static bool header = true;
if (header)
{
for (int i = 1; i <= 7; i++)
_output << "W" << i << "\t";
_output << "log";
for (map<string, vector<double> >::iterator it = _bins.begin(); it != _bins.end(); ++it)
_output << "\t" << it->first;
_output << endl;
header = false;
}
for (int i = 0; i < 7; i++)
_output << scientific << weight[i] << "\t";
_output << l;
for (map<string, vector<double> >::iterator it = _bins.begin(); it != _bins.end(); ++it)
_output << "\t" << _binid[it->first];
_output << endl;
}
/**
* Create linear spaced vector of points.
*/
vector<double> linspace(double const& min, double const& max, int const& size)
{
const double step = std::fabs(max-min)/size;
vector<double> res(size+1, 0);
for (size_t i = 0; i < res.size(); i++)
res[i] = min + i*step;
return res;
}
private:
map<string, vector<double> > _bins;
map<string, int> _binid;
std::ofstream _output;
};
// The singleton which controls the reweight instance
reweight& get_instance()
{
static reweight rw;
return rw;
}
// Fortran interface
extern "C"
{
/**
* Insert histogram to the analysis.
*/
void reweight_book_(const char* key, double const& step,
double const& min, double const& max, int len)
{
string skey(key, len);
get_instance().book(skey, step, min, max);
}
/**
* Fill histogram with weights.
*/
void reweight_fill_(const char* key, double const& obs,
const double* weight, int len)
{
string skey(key, len);
get_instance().fill(skey, obs, weight);
}
/**
* Store the log term to file
*/
void reweight_store_(const double* weight, double const& l)
{
get_instance().store(weight, l);
}
}
|