|
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
|
#ifndef DataFormats_Math_deltaR_h
#define DataFormats_Math_deltaR_h
/* functions to compute deltaR
*
* Ported from original code in RecoJets
* by Fedor Ratnikov, FNAL
*/
#include "DataFormats/Math/interface/deltaPhi.h"
#include <cmath>
namespace reco {
// assumption is that eta and phi are cached AND phi is computed using std::atan2
// type is the type of T1::eta();
template <typename T1, typename T2>
constexpr auto deltaR2(const T1& t1, const T2& t2) -> decltype(t1.eta()) {
typedef decltype(t1.eta()) Float;
Float p1 = t1.phi();
Float p2 = t2.phi();
Float e1 = t1.eta();
Float e2 = t2.eta();
auto dp = std::abs(p1 - p2);
if (dp > Float(M_PI))
dp -= Float(2 * M_PI);
return (e1 - e2) * (e1 - e2) + dp * dp;
}
// do not use it: always cut in deltaR2!
template <typename T1, typename T2>
constexpr auto deltaR(const T1& t1, const T2& t2) -> decltype(t1.eta()) {
return std::sqrt(deltaR2(t1, t2));
}
//prefer the above...
template <class T1, class T2, class T3, class T4>
constexpr T1 deltaR2(T1 eta1, T2 phi1, T3 eta2, T4 phi2) {
T1 deta = eta1 - eta2;
T1 dphi = std::abs(phi1 - phi2);
if (dphi > T1(M_PI))
dphi -= T1(2 * M_PI);
return deta * deta + dphi * dphi;
}
// to be avoided
template <class T1, class T2, class T3, class T4>
constexpr T1 deltaR(T1 eta1, T2 phi1, T3 eta2, T4 phi2) {
return std::sqrt(deltaR2(eta1, phi1, eta2, phi2));
}
} // namespace reco
// woderful! VI
using reco::deltaR;
using reco::deltaR2;
// obsolete use lambdas (and cut in deltaR2!)
template <typename T1, typename T2 = T1>
struct DeltaR {
constexpr double operator()(const T1& t1, const T2& t2) const { return reco::deltaR(t1, t2); }
};
#endif
|