Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DVector2.h
Go to the documentation of this file.
1 // This is a replacement for TVector2 that uses sse2 instructions
2 
3 #ifndef _DVector2_
4 #define _DVector2_
5 
6 #ifndef USE_SSE2
7 
8 #include <TVector2.h>
9 typedef TVector2 DVector2;
10 
11 #else
12 
13 #include <math.h>
14 #include <emmintrin.h>
15 #include <iostream>
16 #include <iomanip>
17 #include <align_16.h>
18 #define RAD2DEG 180./M_PI
19 using namespace std;
20 
21 
22 class DVector2{
23  public:
24  DVector2()
25  : vec( ALIGNED_16_BLOCK_PTR(union dvec, 1, vec) )
26  {
27  vec->v=_mm_setzero_pd();
28  }
29  DVector2(const double a, const double b)
30  : vec( ALIGNED_16_BLOCK_PTR(union dvec, 1, vec) )
31  {
32  vec->v=_mm_setr_pd(a,b);
33  }
34  DVector2(const __m128d v)
35  : vec( ALIGNED_16_BLOCK_PTR(union dvec, 1, vec) )
36  {
37  vec->v=v;
38  }
39  DVector2(const DVector2& dv)
40  : vec( ALIGNED_16_BLOCK_PTR(union dvec, 1, vec) )
41  {
42  vec->v=dv.vec->v;
43  }
44  ~DVector2(){};
45 
46  __m128d GetV() const {return vec->v;}
47 
48  // Assignment
49  DVector2& operator=(const DVector2& dv){
50  vec->v=dv.vec->v;
51  return *this;
52  }
53 
54  // Access by indices
55  double &operator() (int row){
56  return vec->d[row];
57  }
58  double operator() (int row) const{
59  return vec->d[row];
60  }
61 
62  void Set(double a, double b){
63  vec->v=_mm_setr_pd(a,b);
64  }
65 
66  double X() const {return vec->d[0];}
67  double Y() const {return vec->d[1];}
68 
69  // Vector subtraction
70  DVector2 operator-(const DVector2 &v1) const{
71  return DVector2(_mm_sub_pd(GetV(),v1.GetV()));
72  }
73  DVector2 &operator-=(const DVector2 &v1){
74  vec->v=_mm_sub_pd(GetV(),v1.GetV());
75  return *this;
76  }
77  // Vector addition
78  DVector2 operator+(const DVector2 &v1) const{
79  return DVector2(_mm_add_pd(GetV(),v1.GetV()));
80  }
81  DVector2 &operator+=(const DVector2 &v1){
82  vec->v=_mm_add_pd(GetV(),v1.GetV());
83  return *this;
84  }
85  // division by a double
86  DVector2 &operator/=(const double c) {
87  uint32_t padded_space1[8]; // enough space for a block of 16 bytes
88  // that is force-aligned on a 16-byte boundary
89  struct dvec1{
90  __m128d val;
91  }* const scale = reinterpret_cast<struct dvec1 *>
92  ((reinterpret_cast<uintptr_t>(padded_space1)+15) & ~ 15);
93  scale->val=_mm_set1_pd(1./c);
94  vec->v=_mm_mul_pd(vec->v,scale->val);
95  return *this;
96  }
97  // Dot product
98  double operator*(const DVector2 &v1) const{
99  return (X()*v1.X()+Y()*v1.Y());
100  }
101  // multiplication by a double
102  DVector2 &operator*=(const double c){
103  uint32_t padded_space1[8]; // enough space for a block of 16 bytes
104  // that is force-aligned on a 16-byte boundary
105  struct dvec1{
106  __m128d val;
107  }* const scale = reinterpret_cast<struct dvec1 *>
108  ((reinterpret_cast<uintptr_t>(padded_space1)+15) & ~ 15);
109  scale->val=_mm_set1_pd(c);
110  vec->v=_mm_mul_pd(scale->val,GetV());
111  return *this;
112  }
113  double Mod2() const {return (X()*X()+Y()*Y());}
114  double Mod() const {return sqrt(X()*X()+Y()*Y());}
115  double Phi() const {return atan2(Y(),X());}
116 
117  // Angular difference between two vectors
118  double DeltaPhi(const DVector2 &v1) const {
119  double twopi=2.*M_PI;
120  double dphi=Phi()-v1.Phi();
121  while (dphi>=M_PI) dphi-=twopi;
122  while (dphi<-M_PI) dphi+=twopi;
123  return dphi;
124  }
125 
126  // return phi angle between 0 and 2pi
127  double Phi_0_2pi(double angle){
128  double twopi=2.*M_PI;
129  while (angle>=twopi) angle-=twopi;
130  while (angle<0) angle+=twopi;
131  return angle;
132  }
133 
134 
135  void Print(){
136  cout << "DVector2 (x,y)=("<<X()<<","<<Y()<<") (rho,phi)=("<< Mod()
137  <<","<<RAD2DEG*Phi()<<")"<<endl;
138  }
139 
140 
141  private:
142  union dvec{
143  __m128d v;
144  double d[2];
145  };
146  ALIGNED_16_BLOCK(union dvec, 1, vec)
147 };
148 
149 //Scaling by a double
150 inline DVector2 operator*(const double c,const DVector2 &v1){
151  ALIGNED_16_BLOCK_WITH_PTR(__m128d, 1, p)
152  __m128d &scale=p[0];
153  scale=_mm_set1_pd(c);
154  return DVector2(_mm_mul_pd(scale,v1.GetV()));
155 }
156 inline DVector2 operator*(const DVector2 &v1,const double c){
157  ALIGNED_16_BLOCK_WITH_PTR(__m128d, 1, p)
158  __m128d &scale=p[0];
159  scale=_mm_set1_pd(c);
160  return DVector2(_mm_mul_pd(scale,v1.GetV()));
161 }
162 // Division by a double
163 inline DVector2 operator/(const DVector2 &v1,const double c){
164  ALIGNED_16_BLOCK_WITH_PTR(__m128d, 1, p)
165  __m128d &scale=p[0];
166  scale=_mm_set1_pd(1./c);
167  return DVector2(_mm_mul_pd(v1.GetV(),scale));
168 }
169 
170 #endif // USE_SSE2
171 #endif // _DVector2_
DMatrix2x1 operator*(const double c, const DMatrix2x1 &M)
Definition: DMatrix2x1.h:58
#define ALIGNED_16_BLOCK_PTR(TYPE, NUM, PTR)
Definition: align_16.h:24
TVector2 DVector2
Definition: DVector2.h:9
DVector2S operator+(const DVector2S &vec1, const DVector2S &vec2)
Definition: DVector2S.h:54
locHist_ADCmulti Print()
#define c
#define ALIGNED_16_BLOCK_WITH_PTR(TYPE, NUM, PTR)
Definition: align_16.h:27
#define X(str)
Definition: hddm-c.cpp:83
const double RAD2DEG
double sqrt(double)
DVector2S operator-(const DVector2S &vec1, const DVector2S &vec2)
Definition: DVector2S.h:59
#define ALIGNED_16_BLOCK(TYPE, NUM, PTR)
Definition: align_16.h:19