packages/QAUtils/QAUtil.h

0001 #ifndef QA_QAG4UTIL_H
0002 #define QA_QAG4UTIL_H
0003
0004 /*!
0005  * \file QAG4Util.h
0006  * \brief some common utility functions used in detector specific evaluation modules
0007  * \author Hugo Pereira Da Costa <hugo.pereira-da-costa@cea.fr>
0008  */
0009
0010 #include <g4main/PHG4Hit.h>
0011 #include <set>
0012 #include <cmath>
0013
0014 // utility functions
0015 namespace QAG4Util
0016 {
0017   /// square
0018   template <class T>
0019   inline constexpr T square(T x)
0020   {
0021     return x * x;
0022   }
0023
0024   /// radius
0025   template <class T>
0026   inline constexpr T get_r(T x, T y)
0027   {
0028     return std::sqrt(square(x) + square(y));
0029   }
0030
0031   /// angle difference between [-pi, pi[
0032   template <class T>
0033   inline const T delta_phi(T phi1, T phi2)
0034   {
0035     auto out = phi1 - phi2;
0036     while (out >= M_PI) out -= 2 * M_PI;
0037     while (out < -M_PI) out += 2 * M_PI;
0038     return out;
0039   }
0040
0041   /// calculate the best average of member function called on all members in collection, extrapolated at a given radius
0042   /** each hit is weighted by its deposited energy */
0043   template <float (PHG4Hit::*accessor)(int) const>
0044   float interpolate(std::set<PHG4Hit*> hits, float rextrap)
0045   {
0046     // calculate all terms needed for the interpolation
0047     // need to use double everywhere here due to numerical divergences
0048     double sw = 0;
0049     double swr = 0;
0050     double swr2 = 0;
0051     double swx = 0;
0052     double swrx = 0;
0053
0054     bool valid(false);
0055     for (const auto& hit : hits)
0056     {
0057       const double x0 = (hit->*accessor)(0);
0058       const double x1 = (hit->*accessor)(1);
0059       if (std::isnan(x0) || std::isnan(x1)) continue;
0060
0061       const double w = hit->get_edep();
0062       if (w < 0) continue;
0063
0064       valid = true;
0065       const double r0 = get_r(hit->get_x(0), hit->get_y(0));
0066       const double r1 = get_r(hit->get_x(1), hit->get_y(1));
0067
0068       sw += w * 2;
0069       swr += w * (r0 + r1);
0070       swr2 += w * (square(r0) + square(r1));
0071       swx += w * (x0 + x1);
0072       swrx += w * (r0 * x0 + r1 * x1);
0073     }
0074
0075     if (!valid) return NAN;
0076
0077     const auto alpha = (sw * swrx - swr * swx);
0078     const auto beta = (swr2 * swx - swr * swrx);
0079     const auto denom = (sw * swr2 - square(swr));
0080
0081     return (alpha * rextrap + beta) / denom;
0082   }
0083
0084 }  // namespace QAG4Util
0085
0086 #endif