|
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
|
#ifndef ALPHA_T_H
#define ALPHA_T_H
#include <cmath>
#include <vector>
#include <numeric>
#include <functional>
#include <algorithm>
struct alpha_T {
template <class LVs>
double operator()(const LVs& p4s) const {
typedef typename LVs::value_type LorentzV;
if (p4s.size() < 2)
return 0;
std::vector<double> pTs;
transform(p4s.begin(), p4s.end(), back_inserter(pTs), std::mem_fn(&LorentzV::Pt));
const double DsumPT = minimum_deltaSumPT(pTs);
const double sumPT = accumulate(pTs.begin(), pTs.end(), double(0));
const LorentzV sumP4 = accumulate(p4s.begin(), p4s.end(), LorentzV());
return 0.5 * (sumPT - DsumPT) / sqrt(sumPT * sumPT - sumP4.Perp2());
}
static double minimum_deltaSumPT(const std::vector<double>& pTs) {
std::vector<double> diff(1 << (pTs.size() - 1), 0.);
for (unsigned i = 0; i < diff.size(); i++)
for (unsigned j = 0; j < pTs.size(); j++)
diff[i] += pTs[j] * (1 - 2 * (int(i >> j) & 1));
return fabs(*min_element(diff.begin(), diff.end(), [](auto x, auto y) { return fabs(x) < fabs(y); }));
}
};
#endif
|