ParticleFlow/src/FWPFMaths.cc

0001 #include <cmath>
0002 #include "TMath.h"
0003 #include "TEveVector.h"
0004
0005 namespace FWPFMaths {
0006   //______________________________________________________________________________
0007   float sgn(float val) { return (val < 0) ? -1 : 1; }
0008
0009   //______________________________________________________________________________
0010   TEveVector cross(const TEveVector &v1, const TEveVector &v2) {
0011     TEveVector vec;
0012
0013     vec.fX = (v1.fY * v2.fZ) - (v1.fZ * v2.fY);
0014     vec.fY = (v1.fZ * v2.fX) - (v1.fX * v2.fZ);
0015     vec.fZ = (v1.fX * v2.fY) - (v1.fY * v2.fX);
0016
0017     return vec;
0018   }
0019
0020   //______________________________________________________________________________
0021   float dot(const TEveVector &v1, const TEveVector &v2) {
0022     float result = (v1.fX * v2.fX) + (v1.fY * v2.fY) + (v1.fZ * v1.fZ);
0023
0024     return result;
0025   }
0026
0027   //______________________________________________________________________________
0028   TEveVector lineCircleIntersect(const TEveVector &v1, const TEveVector &v2, float r) {
0029     // Definitions
0030     float x, y;
0031     float dx = v1.fX - v2.fX;
0032     float dy = v1.fY - v2.fY;
0033     float dr = sqrt((dx * dx) + (dy * dy));
0034     float D = ((v2.fX * v1.fY) - (v1.fX * v2.fY));
0035
0036     float rtDescrim = sqrt(((r * r) * (dr * dr)) - (D * D));
0037
0038     if (dy < 0)  // Going down
0039     {
0040       x = (D * dy) - ((sgn(dy) * dx) * rtDescrim);
0041       x /= (dr * dr);
0042
0043       y = (-D * dx) - (fabs(dy) * rtDescrim);
0044       y /= (dr * dr);
0045     } else {
0046       x = (D * dy) + ((sgn(dy) * dx) * rtDescrim);
0047       x /= (dr * dr);
0048
0049       y = (-D * dx) + (fabs(dy) * rtDescrim);
0050       y /= (dr * dr);
0051     }
0052
0053     TEveVector result = TEveVector(x, y, 0.001);
0054     return result;
0055   }
0056
0057   //______________________________________________________________________________
0058   TEveVector lineLineIntersect(const TEveVector &p1, const TEveVector &p2, const TEveVector &p3, const TEveVector &p4) {
0059     TEveVector a = p2 - p1;
0060     TEveVector b = p4 - p3;
0061     TEveVector c = p3 - p1;
0062     TEveVector result;
0063     float s, val;
0064
0065     s = dot(cross(c, b), cross(a, b));
0066     val = dot(cross(a, b), cross(a, b));
0067     s /= val;
0068
0069     result = p1 + (a * s);
0070     return result;
0071   }
0072
0073   //______________________________________________________________________________
0074   float linearInterpolation(const TEveVector &p1, const TEveVector &p2, float z) {
0075     float y;
0076
0077     y = ((z - fabs(p1.fZ)) * p2.fY) + ((fabs(p2.fZ) - z) * p1.fY);
0078     y /= (fabs(p2.fZ) - fabs(p1.fZ));
0079
0080     return y;
0081   }
0082
0083   //______________________________________________________________________________
0084   bool checkIntersect(const TEveVector &p, float r) {
0085     float h = sqrt((p.fX * p.fX) + (p.fY * p.fY));
0086
0087     if (h >= r)
0088       return true;
0089
0090     return false;
0091   }
0092
0093   //______________________________________________________________________________
0094   float calculateEt(const TEveVector &centre, float e) {
0095     TEveVector vec = centre;
0096     float et;
0097
0098     vec.Normalize();
0099     vec *= e;
0100     et = vec.Perp();
0101
0102     return et;
0103   }
0104 }  // namespace FWPFMaths