ParticleFlow/src/FWPFMaths.cc

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