Converting TEvtGen to native cmake
[u/mrichter/AliRoot.git] / TEvtGen / EvtGen / EvtGenBase / EvtVector4R.cpp
1 //--------------------------------------------------------------------------
2 //
3 // Environment:
4 //      This software is part of the EvtGen package developed jointly
5 //      for the BaBar and CLEO collaborations.  If you use all or part
6 //      of it, please give an appropriate acknowledgement.
7 //
8 // Copyright Information: See EvtGen/COPYRIGHT
9 //      Copyright (C) 1998      Caltech, UCSB
10 //
11 // Module: EvtVector4R.cc
12 //
13 // Description: Real implementation of 4-vectors
14 //
15 // Modification history:
16 //
17 //    DJL/RYD  September 25, 1996            Module created
18 //
19 //------------------------------------------------------------------------
20 // 
21 #include "EvtGenBase/EvtPatches.hh"
22 #include <iostream>
23 #include <cmath>
24 #include "EvtGenBase/EvtVector4R.hh"
25 #include "EvtGenBase/EvtVector3R.hh"
26 #include "EvtGenBase/EvtVector4C.hh"
27 #include "EvtGenBase/EvtTensor4C.hh"
28
29 using std::ostream;
30
31 EvtVector4R::EvtVector4R() {
32   v[0] = 0.0; v[1] = 0.0; v[2] = 0.0; v[3] = 0.0;
33 }
34
35 EvtVector4R::EvtVector4R(double e,double p1,double p2, double p3){
36   
37   v[0]=e; v[1]=p1; v[2]=p2; v[3]=p3;
38 }
39
40 double EvtVector4R::mass() const{
41
42   double m2=v[0]*v[0]-v[1]*v[1]-v[2]*v[2]-v[3]*v[3];
43
44   if (m2>0.0) {
45     return sqrt(m2);
46   }
47   else{
48     return 0.0;
49   }
50 }
51
52
53 EvtVector4R rotateEuler(const EvtVector4R& rs,
54                         double alpha,double beta,double gamma){
55
56   EvtVector4R tmp(rs);
57   tmp.applyRotateEuler(alpha,beta,gamma);
58   return tmp;
59
60 }
61
62 EvtVector4R boostTo(const EvtVector4R& rs,
63                     const EvtVector4R& p4, bool inverse){
64
65   EvtVector4R tmp(rs);
66   tmp.applyBoostTo(p4, inverse);
67   return tmp;
68
69 }
70
71 EvtVector4R boostTo(const EvtVector4R& rs,
72                     const EvtVector3R& boost, bool inverse){
73
74   EvtVector4R tmp(rs);
75   tmp.applyBoostTo(boost, inverse);
76   return tmp;
77
78 }
79
80
81
82 void EvtVector4R::applyRotateEuler(double phi,double theta,double ksi){
83
84   double sp=sin(phi);
85   double st=sin(theta);
86   double sk=sin(ksi);
87   double cp=cos(phi);
88   double ct=cos(theta);
89   double ck=cos(ksi);
90
91   double x=( ck*ct*cp-sk*sp)*v[1]+( -sk*ct*cp-ck*sp)*v[2]+st*cp*v[3];
92   double y=( ck*ct*sp+sk*cp)*v[1]+(-sk*ct*sp+ck*cp)*v[2]+st*sp*v[3];
93   double z=-ck*st*v[1]+sk*st*v[2]+ct*v[3];
94
95   v[1]=x;
96   v[2]=y;
97   v[3]=z;
98   
99 }
100
101 ostream& operator<<(ostream& s, const EvtVector4R& v){
102
103   s<<"("<<v.v[0]<<","<<v.v[1]<<","<<v.v[2]<<","<<v.v[3]<<")";
104
105   return s;
106
107 }
108
109 void EvtVector4R::applyBoostTo(const EvtVector4R& p4, bool inverse){
110
111   double e=p4.get(0);
112
113   EvtVector3R boost(p4.get(1)/e,p4.get(2)/e,p4.get(3)/e);
114
115   applyBoostTo(boost, inverse);
116
117   return;
118
119 }
120
121 void EvtVector4R::applyBoostTo(const EvtVector3R& boost, bool inverse){
122
123   double bx,by,bz,gamma,b2;
124
125   bx=boost.get(0);
126   by=boost.get(1);
127   bz=boost.get(2);
128
129   double bxx=bx*bx;
130   double byy=by*by;
131   double bzz=bz*bz;
132
133   b2=bxx+byy+bzz;
134
135   if (b2 > 0.0 && b2 < 1.0) {
136
137     gamma=1.0/sqrt(1.0-b2);
138
139     double gb2=(gamma-1.0)/b2;
140
141     double gb2xy=gb2*bx*by;
142     double gb2xz=gb2*bx*bz;
143     double gb2yz=gb2*by*bz;
144
145     double gbx=gamma*bx;
146     double gby=gamma*by;
147     double gbz=gamma*bz;
148
149     double e2=v[0];
150     double px2=v[1];
151     double py2=v[2];
152     double pz2=v[3];
153
154     if ( inverse ) {
155       v[0]=gamma*e2-gbx*px2-gby*py2-gbz*pz2;
156       
157       v[1]=-gbx*e2+gb2*bxx*px2+px2+gb2xy*py2+gb2xz*pz2;
158  
159       v[2]=-gby*e2+gb2*byy*py2+py2+gb2xy*px2+gb2yz*pz2;
160  
161       v[3]=-gbz*e2+gb2*bzz*pz2+pz2+gb2yz*py2+gb2xz*px2;
162     }
163     else {
164       v[0]=gamma*e2+gbx*px2+gby*py2+gbz*pz2;
165       
166       v[1]=gbx*e2+gb2*bxx*px2+px2+gb2xy*py2+gb2xz*pz2;
167  
168       v[2]=gby*e2+gb2*byy*py2+py2+gb2xy*px2+gb2yz*pz2;
169  
170       v[3]=gbz*e2+gb2*bzz*pz2+pz2+gb2yz*py2+gb2xz*px2;
171     }
172   }
173
174 }
175
176 EvtVector4R EvtVector4R::cross( const EvtVector4R& p2 ){
177
178   //Calcs the cross product.  Added by djl on July 27, 1995.
179   //Modified for real vectros by ryd Aug 28-96
180
181   EvtVector4R temp;
182   
183   temp.v[0] = 0.0; 
184   temp.v[1] = v[2]*p2.v[3] - v[3]*p2.v[2];
185   temp.v[2] = v[3]*p2.v[1] - v[1]*p2.v[3];
186   temp.v[3] = v[1]*p2.v[2] - v[2]*p2.v[1];
187
188   return temp;
189 }
190
191 double EvtVector4R::d3mag() const
192
193 // returns the 3 momentum mag.
194 {
195   double temp;
196
197   temp = v[1]*v[1]+v[2]*v[2]+v[3]*v[3];
198
199   temp = sqrt( temp );
200
201   return temp;
202 } // r3mag
203
204 double EvtVector4R::dot ( const EvtVector4R& p2 )const{
205
206   //Returns the dot product of the 3 momentum.  Added by
207   //djl on July 27, 1995.  for real!!!
208
209   double temp;
210
211   temp = v[1]*p2.v[1];
212   temp += v[2]*p2.v[2];
213   temp += v[3]*p2.v[3];
214  
215   return temp;
216
217 } //dot
218
219 // Functions below added by AJB
220
221 // Calculate ( \vec{p1} cross \vec{p2} ) \cdot \vec{p3} in rest frame of object
222 double EvtVector4R::scalartripler3( const EvtVector4R& p1,
223         const EvtVector4R& p2, const EvtVector4R& p3 ) const
224 {
225   EvtVector4C lc=dual(EvtGenFunctions::directProd(*this, p1)).cont2(p2);
226   EvtVector4R  l(real(lc.get(0)), real(lc.get(1)), real(lc.get(2)),
227                  real(lc.get(3)));
228   
229   return -1.0/mass() * (l * p3);
230 }
231
232 // Calculate the 3-d dot product of 4-vectors p1 and p2 in the rest frame of
233 // 4-vector p0
234 double EvtVector4R::dotr3( const EvtVector4R& p1, const EvtVector4R& p2 ) const
235 {
236     return 1/mass2() * ((*this) * p1) * ((*this) * p2) - p1 * p2;
237 }
238
239 // Calculate the 3-d magnitude squared of 4-vector p1 in the rest frame of
240 // 4-vector p0
241 double EvtVector4R::mag2r3( const EvtVector4R& p1 ) const
242 {
243     return Square((*this) * p1)/mass2() - p1.mass2();
244 }
245
246 // Calculate the 3-d magnitude 4-vector p1 in the rest frame of 4-vector p0.
247 double EvtVector4R::magr3( const EvtVector4R& p1 ) const
248 {
249     return sqrt(mag2r3(p1));
250 }
251
252
253
254
255
256
257