ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
QAG4Util.h
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file QAG4Util.h
1 #ifndef QA_QAG4UTIL_H
2 #define QA_QAG4UTIL_H
3 
10 #include <g4main/PHG4Hit.h>
11 
12 #include <cmath>
13 
14 // utility functions
15 namespace QAG4Util
16 {
18  template <class T>
19  inline constexpr T square(T x)
20  {
21  return x * x;
22  }
23 
25  template <class T>
26  inline constexpr T get_r(T x, T y)
27  {
28  return std::sqrt(square(x) + square(y));
29  }
30 
32  template <class T>
33  inline const T delta_phi(T phi1, T phi2)
34  {
35  auto out = phi1 - phi2;
36  while (out >= M_PI) out -= 2 * M_PI;
37  while (out < -M_PI) out += 2 * M_PI;
38  return out;
39  }
40 
42 
43  template <float (PHG4Hit::*accessor)(int) const>
44  float interpolate(std::set<PHG4Hit*> hits, float rextrap)
45  {
46  // calculate all terms needed for the interpolation
47  // need to use double everywhere here due to numerical divergences
48  double sw = 0;
49  double swr = 0;
50  double swr2 = 0;
51  double swx = 0;
52  double swrx = 0;
53 
54  bool valid(false);
55  for (const auto& hit : hits)
56  {
57  const double x0 = (hit->*accessor)(0);
58  const double x1 = (hit->*accessor)(1);
59  if (std::isnan(x0) || std::isnan(x1)) continue;
60 
61  const double w = hit->get_edep();
62  if (w < 0) continue;
63 
64  valid = true;
65  const double r0 = get_r(hit->get_x(0), hit->get_y(0));
66  const double r1 = get_r(hit->get_x(1), hit->get_y(1));
67 
68  sw += w * 2;
69  swr += w * (r0 + r1);
70  swr2 += w * (square(r0) + square(r1));
71  swx += w * (x0 + x1);
72  swrx += w * (r0 * x0 + r1 * x1);
73  }
74 
75  if (!valid) return NAN;
76 
77  const auto alpha = (sw * swrx - swr * swx);
78  const auto beta = (swr2 * swx - swr * swrx);
79  const auto denom = (sw * swr2 - square(swr));
80 
81  return (alpha * rextrap + beta) / denom;
82  }
83 
84 } // namespace QAG4Util
85 
86 #endif