File indexing completed on 2024-04-06 12:09:25
0001 #ifndef ALPHA_T_H
0002 #define ALPHA_T_H
0003
0004 #include <cmath>
0005 #include <vector>
0006 #include <numeric>
0007 #include <functional>
0008 #include <algorithm>
0009
0010 struct alpha_T {
0011 template <class LVs>
0012 double operator()(const LVs& p4s) const {
0013 typedef typename LVs::value_type LorentzV;
0014 if (p4s.size() < 2)
0015 return 0;
0016
0017 std::vector<double> pTs;
0018 transform(p4s.begin(), p4s.end(), back_inserter(pTs), std::mem_fn(&LorentzV::Pt));
0019
0020 const double DsumPT = minimum_deltaSumPT(pTs);
0021 const double sumPT = accumulate(pTs.begin(), pTs.end(), double(0));
0022 const LorentzV sumP4 = accumulate(p4s.begin(), p4s.end(), LorentzV());
0023
0024 return 0.5 * (sumPT - DsumPT) / sqrt(sumPT * sumPT - sumP4.Perp2());
0025 }
0026
0027 static double minimum_deltaSumPT(const std::vector<double>& pTs) {
0028 std::vector<double> diff(1 << (pTs.size() - 1), 0.);
0029 for (unsigned i = 0; i < diff.size(); i++)
0030 for (unsigned j = 0; j < pTs.size(); j++)
0031 diff[i] += pTs[j] * (1 - 2 * (int(i >> j) & 1));
0032
0033 return fabs(*min_element(diff.begin(), diff.end(), [](auto x, auto y) { return fabs(x) < fabs(y); }));
0034 }
0035 };
0036
0037 #endif
