1 // -*- C++ -*-
2 // ---------------------------------------------------------------------------
3 //
4 // This file is a part of the CLHEP - a Class Library for High Energy Physics.
5 //
6 // This is the implementation of those parts of the HepLorentzRotation class
7 // which involve decomposition into Boost*Rotation.
8
9 #ifdef GNUPRAGMA
10 #pragma implementation
11 #endif
12
13 #include "CLHEP/Vector/defs.h"
14 #include "CLHEP/Vector/LorentzRotation.h"
15
16 namespace CLHEP {
17
18 // ---------- Decomposition:
19
decompose(HepBoost & bboost,HepRotation & rotation) const20 void HepLorentzRotation::decompose
21 (HepBoost & bboost, HepRotation & rotation) const {
22
23 // The boost will be the pure boost based on column 4 of the transformation
24 // matrix. Since the constructor takes the beta vector, and not beta*gamma,
25 // we first divide through by gamma = the tt element. This of course can
26 // never be zero since the last row has t**2 - v**2 = +1.
27
28 Hep3Vector betaVec ( xt(), yt(), zt() );
29 betaVec *= 1.0 / tt();
30 bboost.set( betaVec );
31
32 // The rotation will be inverse of B times T.
33
34 HepBoost B( -betaVec );
35 HepLorentzRotation R( B * *this );
36
37 HepRep3x3 m1 ( R.xx(), R.xy(), R.xz(),
38 R.yx(), R.yy(), R.yz(),
39 R.zx(), R.zy(), R.zz() );
40 rotation.set( m1 );
41 rotation.rectify();
42
43 return;
44
45 }
46
decompose(Hep3Vector & bboost,HepAxisAngle & rotation) const47 void HepLorentzRotation::decompose
48 (Hep3Vector & bboost, HepAxisAngle & rotation) const {
49 HepRotation r;
50 HepBoost b;
51 decompose(b,r);
52 bboost = b.boostVector();
53 rotation = r.axisAngle();
54 return;
55 }
56
decompose(HepRotation & rotation,HepBoost & bboost) const57 void HepLorentzRotation::decompose
58 (HepRotation & rotation, HepBoost & bboost) const {
59
60 // In this case the pure boost is based on row 4 of the matrix.
61
62 Hep3Vector betaVec( tx(), ty(), tz() );
63 betaVec *= 1.0 / tt();
64 bboost.set( betaVec );
65
66 // The rotation will be T times the inverse of B.
67
68 HepBoost B( -betaVec );
69 HepLorentzRotation R( *this * B );
70
71 HepRep3x3 m1 ( R.xx(), R.xy(), R.xz(),
72 R.yx(), R.yy(), R.yz(),
73 R.zx(), R.zy(), R.zz() );
74 rotation.set( m1 );
75 rotation.rectify();
76 return;
77
78 }
79
decompose(HepAxisAngle & rotation,Hep3Vector & bboost) const80 void HepLorentzRotation::decompose
81 (HepAxisAngle & rotation, Hep3Vector & bboost) const {
82 HepRotation r;
83 HepBoost b;
84 decompose(r,b);
85 rotation = r.axisAngle();
86 bboost = b.boostVector();
87 return;
88 }
89
distance2(const HepBoost & b) const90 double HepLorentzRotation::distance2( const HepBoost & b ) const {
91 HepBoost b1;
92 HepRotation r1;
93 decompose( b1, r1 );
94 double db2 = b1.distance2( b );
95 double dr2 = r1.norm2();
96 return ( db2 + dr2 );
97 }
98
distance2(const HepRotation & r) const99 double HepLorentzRotation::distance2( const HepRotation & r ) const {
100 HepBoost b1;
101 HepRotation r1;
102 decompose( b1, r1 );
103 double db2 = b1.norm2( );
104 double dr2 = r1.distance2( r );
105 return ( db2 + dr2 );
106 }
107
distance2(const HepLorentzRotation & lt) const108 double HepLorentzRotation::distance2(
109 const HepLorentzRotation & lt ) const {
110 HepBoost b1;
111 HepRotation r1;
112 decompose( b1, r1 );
113 HepBoost b2;
114 HepRotation r2;
115 lt.decompose (b2, r2);
116 double db2 = b1.distance2( b2 );
117 double dr2 = r1.distance2( r2 );
118 return ( db2 + dr2 );
119 }
120
howNear(const HepBoost & b) const121 double HepLorentzRotation::howNear( const HepBoost & b ) const {
122 return std::sqrt( distance2( b ) );
123 }
howNear(const HepRotation & r) const124 double HepLorentzRotation::howNear( const HepRotation & r ) const {
125 return std::sqrt( distance2( r ) );
126 }
howNear(const HepLorentzRotation & lt) const127 double HepLorentzRotation::howNear( const HepLorentzRotation & lt )const {
128 return std::sqrt( distance2( lt ) );
129 }
130
isNear(const HepBoost & b,double epsilon) const131 bool HepLorentzRotation::isNear(
132 const HepBoost & b, double epsilon ) const {
133 HepBoost b1;
134 HepRotation r1;
135 decompose( b1, r1 );
136 double db2 = b1.distance2(b);
137 if ( db2 > epsilon*epsilon ) {
138 return false; // Saves the time-consuming Rotation::norm2
139 }
140 double dr2 = r1.norm2();
141 return ( (db2 + dr2) <= epsilon*epsilon );
142 }
143
isNear(const HepRotation & r,double epsilon) const144 bool HepLorentzRotation::isNear(
145 const HepRotation & r, double epsilon ) const {
146 HepBoost b1;
147 HepRotation r1;
148 decompose( b1, r1 );
149 double db2 = b1.norm2();
150 if ( db2 > epsilon*epsilon ) {
151 return false; // Saves the time-consuming Rotation::distance2
152 }
153 double dr2 = r1.distance2(r);
154 return ( (db2 + dr2) <= epsilon*epsilon );
155 }
156
isNear(const HepLorentzRotation & lt,double epsilon) const157 bool HepLorentzRotation::isNear(
158 const HepLorentzRotation & lt, double epsilon ) const {
159 HepBoost b1;
160 HepRotation r1;
161 decompose( b1, r1 );
162 HepBoost b2;
163 HepRotation r2;
164 lt.decompose (b2, r2);
165 double db2 = b1.distance2(b2);
166 if ( db2 > epsilon*epsilon ) {
167 return false; // Saves the time-consuming Rotation::distance2
168 }
169 double dr2 = r1.distance2(r2);
170 return ( (db2 + dr2) <= epsilon*epsilon );
171 }
172
norm2() const173 double HepLorentzRotation::norm2() const {
174 HepBoost b;
175 HepRotation r;
176 decompose( b, r );
177 return b.norm2() + r.norm2();
178 }
179
rectify()180 void HepLorentzRotation::rectify() {
181
182 // Assuming the representation of this is close to a true LT,
183 // but may have drifted due to round-off error from many operations,
184 // this forms an "exact" orthosymplectic matrix for the LT again.
185
186 // There are several ways to do this, all equivalent to lowest order in
187 // the corrected error. We choose to form an LT based on the inverse boost
188 // extracted from row 4, and left-multiply by LT to form what would be
189 // a rotation if the LT were kosher. We drop the possibly non-zero t
190 // components of that, rectify that rotation and multiply back by the boost.
191
192 Hep3Vector beta (tx(), ty(), tz());
193 double gam = tt(); // NaN-proofing
194 if ( gam <= 0 ) {
195 ZMthrowA ( ZMxpvImproperTransformation (
196 "rectify() on a transformation with tt() <= 0 - will not help!" ));
197 gam = 1;
198 }
199 beta *= 1.0/gam;
200 HepLorentzRotation R = (*this) * HepBoost(-beta);
201
202 HepRep3x3 m1 ( R.xx(), R.xy(), R.xz(),
203 R.yx(), R.yy(), R.yz(),
204 R.zx(), R.zy(), R.zz() );
205
206 HepRotation Rgood (m1);
207 Rgood.rectify();
208
209 set ( Rgood, HepBoost(beta) );
210 }
211
212 } // namespace CLHEP
213