JetProducers/interface/VirtualJetProducerHelper.h

0001 #ifndef RecoJets_JetProducers_interface_VirtualJetProducerHelper_h
0002 #define RecoJets_JetProducers_interface_VirtualJetProducerHelper_h
0003
0004 #include <cmath>
0005 #include <algorithm>
0006
0007 namespace reco {
0008
0009   namespace helper {
0010
0011     namespace VirtualJetProducerHelper {
0012
0013       // Area of intersection of two unit-radius disks with centers separated by r12.
0014       inline double intersection(double r12) {
0015         if (r12 == 0)
0016           return M_PI;
0017         if (r12 >= 2)
0018           return 0;
0019         return 2 * std::acos(0.5 * r12) - 0.5 * r12 * sqrt(std::max(0.0, 4 - r12 * r12));
0020       }
0021
0022       // Area of intersection of three unit-radius disks with centers separated by r12, r23, r13.
0023       inline double intersection(double r12, double r23, double r13) {
0024         if (r12 >= 2 || r23 >= 2 || r13 >= 2)
0025           return 0;
0026         const double r12_2 = r12 * r12;
0027         const double r13_2 = r13 * r13;
0028         const double temp = (r12_2 + r13_2 - r23 * r23);
0029         const double T2 = std::max(0.0, 4 * r12_2 * r13_2 - temp * temp);
0030         const double common = 0.5 * (intersection(r12) + intersection(r13) + intersection(r23) - M_PI + 0.5 * sqrt(T2));
0031         return common;
0032       }
0033       inline double intersection(double r12, double r23, double r13, double a12, double a23, double a13) {
0034         if (r12 >= 2 || r23 >= 2 || r13 >= 2)
0035           return 0;
0036         const double r12_2 = r12 * r12;
0037         const double r13_2 = r13 * r13;
0038         const double temp = (r12_2 + r13_2 - r23 * r23);
0039         const double T2 = std::max(0.0, 4 * r12_2 * r13_2 - temp * temp);
0040         const double common = 0.5 * (a12 + a13 + a23 - M_PI + 0.5 * sqrt(T2));
0041         return common;
0042       }
0043
0044     }  // namespace VirtualJetProducerHelper
0045   }  // namespace helper
0046 }  // namespace reco
0047 #endif