EcalClusters/interface/AnglesUtil.h

0001 #ifndef INC_ANGLESUTIL
0002 #define INC_ANGLESUTIL
0003 ///////////////////////////////////////////////////////////////////////////////
0004 // File: AnglesUtil.hpp
0005 //
0006 // Purpose:  Provide useful functions for calculating angles and eta
0007 //
0008 // Created:   4-NOV-1998  Serban Protopopescu
0009 // History:   replaced old KinemUtil
0010 // Modified:  23-July-2000 Add rapidity calculation
0011 //            14-Aug-2003 converted all functions to double precision
0012 ///////////////////////////////////////////////////////////////////////////////
0013 // Dependencies (#includes)
0014
0015 #include <cmath>
0016 #include <cstdlib>
0017
0018 namespace kinem {
0019   const double PI = 2.0 * acos(0.);
0020   const double TWOPI = 2.0 * PI;
0021   const float ETA_LIMIT = 15.0;
0022   const float EPSILON = 1.E-10;
0023
0024   //  calculate phi from x, y
0025   inline double phi(double x, double y);
0026   //  calculate phi for a line defined by xy1 and xy2 (xy2-xy1)
0027   inline double phi(double xy1[2], double xy2[2]);
0028   inline double phi(float xy1[2], float xy2[2]);
0029
0030   //  calculate theta from x, y, z
0031   inline double theta(double x, double y, double z);
0032   //  calculate theta for a line defined by xyz1 and xyz2 (xyz2-xyz1)
0033   inline double theta(double xyz1[3], double xyz2[3]);
0034   inline double theta(float xyz1[3], float xyz2[3]);
0035   //  calculate theta from eta
0036   inline double theta(double etap);
0037
0038   //  calculate eta from x, y, z (return also theta)
0039   inline double eta(double x, double y, double z);
0040   //  calculate eta for a line defined by xyz1 and xyz2 (xyz2-xyz1)
0041   inline double eta(double xyz1[3], double xyz2[3]);
0042   inline double eta(float xyz1[3], float xyz2[3]);
0043   //  calculate eta from theta
0044   inline double eta(double th);
0045
0046   // calculate rapidity from E, pz
0047   inline double y(double E, double pz);
0048
0049   //  calculate phi1-phi2 keeping value between 0 and pi
0050   inline double delta_phi(double ph11, double phi2);
0051   // calculate phi1-phi2 keeping value between -pi and pi
0052   inline double signed_delta_phi(double ph11, double phi2);
0053   // calculate eta1 - eta2
0054   inline double delta_eta(double eta1, double eta2);
0055
0056   //  calculate deltaR
0057   inline double delta_R(double eta1, double phi1, double eta2, double phi2);
0058
0059   //  calculate unit vectors given two points
0060   inline void uvectors(double u[3], double xyz1[3], double xyz2[3]);
0061   inline void uvectors(float u[3], float xyz1[3], float xyz2[3]);
0062
0063   inline double tanl_from_theta(double theta);
0064   inline double theta_from_tanl(double tanl);
0065 }  // namespace kinem
0066
0067 inline double kinem::phi(double x, double y) {
0068   double PHI = atan2(y, x);
0069   return (PHI >= 0) ? PHI : kinem::TWOPI + PHI;
0070 }
0071
0072 inline double kinem::phi(float xy1[2], float xy2[2]) {
0073   double dxy1[2] = {xy1[0], xy1[1]};
0074   double dxy2[2] = {xy2[0], xy2[1]};
0075   return phi(dxy1, dxy2);
0076 }
0077
0078 inline double kinem::phi(double xy1[2], double xy2[2]) {
0079   double x = xy2[0] - xy1[0];
0080   double y = xy2[1] - xy1[1];
0081   return phi(x, y);
0082 }
0083
0084 inline double kinem::delta_phi(double phi1, double phi2) {
0085   double PHI = fabs(phi1 - phi2);
0086   return (PHI <= PI) ? PHI : kinem::TWOPI - PHI;
0087 }
0088
0089 inline double kinem::delta_eta(double eta1, double eta2) { return eta1 - eta2; }
0090
0091 inline double kinem::signed_delta_phi(double phi1, double phi2) {
0092   double phia = phi1;
0093   if (phi1 > PI)
0094     phia = phi1 - kinem::TWOPI;
0095   double phib = phi2;
0096   if (phi2 > PI)
0097     phib = phi2 - kinem::TWOPI;
0098   double dphi = phia - phib;
0099   if (dphi > PI)
0100     dphi -= kinem::TWOPI;
0101   if (dphi < -PI)
0102     dphi += kinem::TWOPI;
0103   return dphi;
0104 }
0105
0106 inline double kinem::delta_R(double eta1, double phi1, double eta2, double phi2) {
0107   double deta = eta1 - eta2;
0108   double dphi = kinem::delta_phi(phi1, phi2);
0109   return sqrt(deta * deta + dphi * dphi);
0110 }
0111
0112 inline double kinem::theta(double xyz1[3], double xyz2[3]) {
0113   double x = xyz2[0] - xyz1[0];
0114   double y = xyz2[1] - xyz1[1];
0115   double z = xyz2[2] - xyz1[2];
0116   return theta(x, y, z);
0117 }
0118
0119 inline double kinem::theta(float xyz1[3], float xyz2[3]) {
0120   double dxyz1[3] = {xyz1[0], xyz1[1], xyz1[2]};
0121   double dxyz2[3] = {xyz2[0], xyz2[1], xyz2[2]};
0122   return theta(dxyz1, dxyz2);
0123 }
0124
0125 inline double kinem::theta(double x, double y, double z) { return atan2(sqrt(x * x + y * y), z); }
0126
0127 inline double kinem::theta(double etap) { return 2.0 * atan(exp(-etap)); }
0128
0129 inline double kinem::eta(double xyz1[3], double xyz2[3]) {
0130   double x = xyz2[0] - xyz1[0];
0131   double y = xyz2[1] - xyz1[1];
0132   double z = xyz2[2] - xyz1[2];
0133   return eta(x, y, z);
0134 }
0135
0136 inline double kinem::eta(float xyz1[3], float xyz2[3]) {
0137   double dxyz1[3] = {xyz1[0], xyz1[1], xyz1[2]};
0138   double dxyz2[3] = {xyz2[0], xyz2[1], xyz2[2]};
0139   return eta(dxyz1, dxyz2);
0140 }
0141
0142 inline double kinem::eta(double x, double y, double z) {
0143   return 0.5 * log((sqrt(x * x + y * y + z * z) + z + EPSILON) / (sqrt(x * x + y * y + z * z) - z + EPSILON));
0144 }
0145
0146 inline double kinem::eta(double th) {
0147   if (th == 0)
0148     return ETA_LIMIT;
0149   if (th >= PI - 0.0001)
0150     return -ETA_LIMIT;
0151   return -log(tan(th / 2.0));
0152 }
0153
0154 inline double kinem::y(double E, double pz) { return 0.5 * log((E + pz + EPSILON) / (E - pz + EPSILON)); }
0155
0156 inline void kinem::uvectors(double u[3], double xyz1[3], double xyz2[3]) {
0157   double xdiff = xyz2[0] - xyz1[0];
0158   double ydiff = xyz2[1] - xyz1[1];
0159   double zdiff = xyz2[2] - xyz1[2];
0160   double s = sqrt(xdiff * xdiff + ydiff * ydiff + zdiff * zdiff);
0161   if (s > 0) {
0162     u[0] = xdiff / s;
0163     u[1] = ydiff / s;
0164     u[2] = zdiff / s;
0165   } else {
0166     u[0] = 0;
0167     u[1] = 0;
0168     u[2] = 0;
0169   }
0170 }
0171
0172 inline void kinem::uvectors(float u[3], float xyz1[3], float xyz2[3]) {
0173   double du[3];
0174   double dxyz1[3] = {xyz1[0], xyz1[1], xyz1[2]};
0175   double dxyz2[3] = {xyz2[0], xyz2[1], xyz2[2]};
0176   uvectors(du, dxyz1, dxyz2);
0177   u[0] = du[0];
0178   u[1] = du[1];
0179   u[2] = du[2];
0180 }
0181
0182 inline double kinem::tanl_from_theta(double theta) { return tan(PI / 2.0 - theta); }
0183
0184 inline double kinem::theta_from_tanl(double tanl) { return PI / 2.0 - atan(tanl); }
0185
0186 #endif  // INC_ANGLESUTIL