Line data Source code
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 0 : EvtVector4R::EvtVector4R() {
32 0 : v[0] = 0.0; v[1] = 0.0; v[2] = 0.0; v[3] = 0.0;
33 0 : }
34 :
35 0 : EvtVector4R::EvtVector4R(double e,double p1,double p2, double p3){
36 :
37 0 : v[0]=e; v[1]=p1; v[2]=p2; v[3]=p3;
38 0 : }
39 :
40 : double EvtVector4R::mass() const{
41 :
42 0 : double m2=v[0]*v[0]-v[1]*v[1]-v[2]*v[2]-v[3]*v[3];
43 :
44 0 : if (m2>0.0) {
45 0 : return sqrt(m2);
46 : }
47 : else{
48 0 : return 0.0;
49 : }
50 0 : }
51 :
52 :
53 : EvtVector4R rotateEuler(const EvtVector4R& rs,
54 : double alpha,double beta,double gamma){
55 :
56 0 : EvtVector4R tmp(rs);
57 0 : tmp.applyRotateEuler(alpha,beta,gamma);
58 0 : return tmp;
59 :
60 : }
61 :
62 : EvtVector4R boostTo(const EvtVector4R& rs,
63 : const EvtVector4R& p4, bool inverse){
64 :
65 0 : EvtVector4R tmp(rs);
66 0 : tmp.applyBoostTo(p4, inverse);
67 0 : return tmp;
68 :
69 : }
70 :
71 : EvtVector4R boostTo(const EvtVector4R& rs,
72 : const EvtVector3R& boost, bool inverse){
73 :
74 0 : EvtVector4R tmp(rs);
75 0 : tmp.applyBoostTo(boost, inverse);
76 0 : return tmp;
77 :
78 : }
79 :
80 :
81 :
82 : void EvtVector4R::applyRotateEuler(double phi,double theta,double ksi){
83 :
84 0 : double sp=sin(phi);
85 0 : double st=sin(theta);
86 0 : double sk=sin(ksi);
87 0 : double cp=cos(phi);
88 0 : double ct=cos(theta);
89 0 : double ck=cos(ksi);
90 :
91 0 : double x=( ck*ct*cp-sk*sp)*v[1]+( -sk*ct*cp-ck*sp)*v[2]+st*cp*v[3];
92 0 : double y=( ck*ct*sp+sk*cp)*v[1]+(-sk*ct*sp+ck*cp)*v[2]+st*sp*v[3];
93 0 : double z=-ck*st*v[1]+sk*st*v[2]+ct*v[3];
94 :
95 0 : v[1]=x;
96 0 : v[2]=y;
97 0 : v[3]=z;
98 :
99 0 : }
100 :
101 : ostream& operator<<(ostream& s, const EvtVector4R& v){
102 :
103 0 : s<<"("<<v.v[0]<<","<<v.v[1]<<","<<v.v[2]<<","<<v.v[3]<<")";
104 :
105 0 : return s;
106 :
107 : }
108 :
109 : void EvtVector4R::applyBoostTo(const EvtVector4R& p4, bool inverse){
110 :
111 0 : double e=p4.get(0);
112 :
113 0 : EvtVector3R boost(p4.get(1)/e,p4.get(2)/e,p4.get(3)/e);
114 :
115 0 : applyBoostTo(boost, inverse);
116 :
117 : return;
118 :
119 0 : }
120 :
121 : void EvtVector4R::applyBoostTo(const EvtVector3R& boost, bool inverse){
122 :
123 : double bx,by,bz,gamma,b2;
124 :
125 0 : bx=boost.get(0);
126 0 : by=boost.get(1);
127 0 : bz=boost.get(2);
128 :
129 0 : double bxx=bx*bx;
130 0 : double byy=by*by;
131 0 : double bzz=bz*bz;
132 :
133 0 : b2=bxx+byy+bzz;
134 :
135 0 : if (b2 > 0.0 && b2 < 1.0) {
136 :
137 0 : gamma=1.0/sqrt(1.0-b2);
138 :
139 0 : double gb2=(gamma-1.0)/b2;
140 :
141 0 : double gb2xy=gb2*bx*by;
142 0 : double gb2xz=gb2*bx*bz;
143 0 : double gb2yz=gb2*by*bz;
144 :
145 0 : double gbx=gamma*bx;
146 0 : double gby=gamma*by;
147 0 : double gbz=gamma*bz;
148 :
149 0 : double e2=v[0];
150 0 : double px2=v[1];
151 0 : double py2=v[2];
152 0 : double pz2=v[3];
153 :
154 0 : if ( inverse ) {
155 0 : v[0]=gamma*e2-gbx*px2-gby*py2-gbz*pz2;
156 :
157 0 : v[1]=-gbx*e2+gb2*bxx*px2+px2+gb2xy*py2+gb2xz*pz2;
158 :
159 0 : v[2]=-gby*e2+gb2*byy*py2+py2+gb2xy*px2+gb2yz*pz2;
160 :
161 0 : v[3]=-gbz*e2+gb2*bzz*pz2+pz2+gb2yz*py2+gb2xz*px2;
162 0 : }
163 : else {
164 0 : v[0]=gamma*e2+gbx*px2+gby*py2+gbz*pz2;
165 :
166 0 : v[1]=gbx*e2+gb2*bxx*px2+px2+gb2xy*py2+gb2xz*pz2;
167 :
168 0 : v[2]=gby*e2+gb2*byy*py2+py2+gb2xy*px2+gb2yz*pz2;
169 :
170 0 : v[3]=gbz*e2+gb2*bzz*pz2+pz2+gb2yz*py2+gb2xz*px2;
171 : }
172 0 : }
173 :
174 0 : }
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 0 : EvtVector4R temp;
182 :
183 0 : temp.v[0] = 0.0;
184 0 : temp.v[1] = v[2]*p2.v[3] - v[3]*p2.v[2];
185 0 : temp.v[2] = v[3]*p2.v[1] - v[1]*p2.v[3];
186 0 : temp.v[3] = v[1]*p2.v[2] - v[2]*p2.v[1];
187 :
188 0 : return temp;
189 : }
190 :
191 : double EvtVector4R::d3mag() const
192 :
193 : // returns the 3 momentum mag.
194 : {
195 : double temp;
196 :
197 0 : temp = v[1]*v[1]+v[2]*v[2]+v[3]*v[3];
198 :
199 0 : temp = sqrt( temp );
200 :
201 0 : 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 0 : temp = v[1]*p2.v[1];
212 0 : temp += v[2]*p2.v[2];
213 0 : temp += v[3]*p2.v[3];
214 :
215 0 : 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 0 : EvtVector4C lc=dual(EvtGenFunctions::directProd(*this, p1)).cont2(p2);
226 0 : EvtVector4R l(real(lc.get(0)), real(lc.get(1)), real(lc.get(2)),
227 0 : real(lc.get(3)));
228 :
229 0 : return -1.0/mass() * (l * p3);
230 0 : }
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 0 : 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 0 : 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 0 : return sqrt(mag2r3(p1));
250 : }
251 :
252 :
253 :
254 :
255 :
256 :
257 :
|
