| Geant4 Cross Reference |
1 // -*- C++ -*- 1 // -*- C++ -*-
2 // ------------------------------------------- 2 // ---------------------------------------------------------------------------
3 // 3 //
4 // This file is a part of the CLHEP - a Class 4 // This file is a part of the CLHEP - a Class Library for High Energy Physics.
5 // 5 //
6 // This is the implementation of the HepBoostY 6 // This is the implementation of the HepBoostY class.
7 // 7 //
8 8
>> 9 #ifdef GNUPRAGMA
>> 10 #pragma implementation
>> 11 #endif
>> 12
9 #include "CLHEP/Vector/BoostY.h" 13 #include "CLHEP/Vector/BoostY.h"
10 #include "CLHEP/Vector/Boost.h" 14 #include "CLHEP/Vector/Boost.h"
11 #include "CLHEP/Vector/Rotation.h" 15 #include "CLHEP/Vector/Rotation.h"
12 #include "CLHEP/Vector/LorentzRotation.h" 16 #include "CLHEP/Vector/LorentzRotation.h"
13 <<
14 #include <cmath> <<
15 #include <iostream> <<
16 17
17 namespace CLHEP { 18 namespace CLHEP {
18 19
19 // ---------- Constructors and Assignment: 20 // ---------- Constructors and Assignment:
20 21
21 HepBoostY & HepBoostY::set (double bbeta) { 22 HepBoostY & HepBoostY::set (double bbeta) {
22 double b2 = bbeta*bbeta; 23 double b2 = bbeta*bbeta;
23 if (b2 >= 1) { 24 if (b2 >= 1) {
24 std::cerr << "HepBoostY::set() - " 25 std::cerr << "HepBoostY::set() - "
25 << "Beta supplied to set HepBoostY repre 26 << "Beta supplied to set HepBoostY represents speed >= c." << std::endl;
26 beta_ = 1.0 - 1.0E-8; // NaN 27 beta_ = 1.0 - 1.0E-8; // NaN-proofing
27 gamma_ = 1.0 / std::sqrt(1.0 - b2); 28 gamma_ = 1.0 / std::sqrt(1.0 - b2);
28 return *this; 29 return *this;
29 } 30 }
30 beta_ = bbeta; 31 beta_ = bbeta;
31 gamma_ = 1.0 / std::sqrt(1.0 - b2); 32 gamma_ = 1.0 / std::sqrt(1.0 - b2);
32 return *this; 33 return *this;
33 } 34 }
34 35
35 // ---------- Accessors: 36 // ---------- Accessors:
36 37
37 HepRep4x4 HepBoostY::rep4x4() const { 38 HepRep4x4 HepBoostY::rep4x4() const {
38 double bg = beta_*gamma_; 39 double bg = beta_*gamma_;
39 return HepRep4x4( 1, 0, 0, 0, 40 return HepRep4x4( 1, 0, 0, 0,
40 0, gamma_, 0, bg, 41 0, gamma_, 0, bg,
41 0, 0, 1, 0, 42 0, 0, 1, 0,
42 0, bg, 0, gamma_ ); 43 0, bg, 0, gamma_ );
43 } 44 }
44 45
45 HepRep4x4Symmetric HepBoostY::rep4x4Symmetric( 46 HepRep4x4Symmetric HepBoostY::rep4x4Symmetric() const {
46 double bg = beta_*gamma_; 47 double bg = beta_*gamma_;
47 return HepRep4x4Symmetric ( 1, 0, 48 return HepRep4x4Symmetric ( 1, 0, 0, 0,
48 gamma_, 0, bg, 49 gamma_, 0, bg,
49 1, 0, 50 1, 0,
50 gamma_ ) 51 gamma_ );
51 } 52 }
52 53
53 // ---------- Decomposition: 54 // ---------- Decomposition:
54 55
55 void HepBoostY::decompose (HepRotation & rotat 56 void HepBoostY::decompose (HepRotation & rotation, HepBoost & boost) const {
56 HepAxisAngle vdelta = HepAxisAngle(); 57 HepAxisAngle vdelta = HepAxisAngle();
57 rotation = HepRotation(vdelta); 58 rotation = HepRotation(vdelta);
58 Hep3Vector bbeta = boostVector(); 59 Hep3Vector bbeta = boostVector();
59 boost = HepBoost(bbeta); 60 boost = HepBoost(bbeta);
60 } 61 }
61 62
62 void HepBoostY::decompose (HepAxisAngle & rota 63 void HepBoostY::decompose (HepAxisAngle & rotation, Hep3Vector & boost) const {
63 rotation = HepAxisAngle(); 64 rotation = HepAxisAngle();
64 boost = boostVector(); 65 boost = boostVector();
65 } 66 }
66 67
67 void HepBoostY::decompose (HepBoost & boost, H 68 void HepBoostY::decompose (HepBoost & boost, HepRotation & rotation) const {
68 HepAxisAngle vdelta = HepAxisAngle(); 69 HepAxisAngle vdelta = HepAxisAngle();
69 rotation = HepRotation(vdelta); 70 rotation = HepRotation(vdelta);
70 Hep3Vector bbeta = boostVector(); 71 Hep3Vector bbeta = boostVector();
71 boost = HepBoost(bbeta); 72 boost = HepBoost(bbeta);
72 } 73 }
73 74
74 void HepBoostY::decompose (Hep3Vector & boost, 75 void HepBoostY::decompose (Hep3Vector & boost, HepAxisAngle & rotation) const {
75 rotation = HepAxisAngle(); 76 rotation = HepAxisAngle();
76 boost = boostVector(); 77 boost = boostVector();
77 } 78 }
78 79
79 // ---------- Comparisons: 80 // ---------- Comparisons:
80 81
81 double HepBoostY::distance2( const HepBoost & 82 double HepBoostY::distance2( const HepBoost & b ) const {
82 return b.distance2(*this); 83 return b.distance2(*this);
83 } 84 }
84 85
85 double HepBoostY::distance2( const HepRotation 86 double HepBoostY::distance2( const HepRotation & r ) const {
86 double db2 = norm2(); 87 double db2 = norm2();
87 double dr2 = r.norm2(); 88 double dr2 = r.norm2();
88 return (db2 + dr2); 89 return (db2 + dr2);
89 } 90 }
90 91
91 double HepBoostY::distance2( const HepLorentzR 92 double HepBoostY::distance2( const HepLorentzRotation & lt ) const {
92 HepBoost b1; 93 HepBoost b1;
93 HepRotation r1; 94 HepRotation r1;
94 lt.decompose(b1,r1); 95 lt.decompose(b1,r1);
95 double db2 = distance2(b1); 96 double db2 = distance2(b1);
96 double dr2 = r1.norm2(); 97 double dr2 = r1.norm2();
97 return (db2 + dr2); 98 return (db2 + dr2);
98 } 99 }
99 100
100 bool HepBoostY::isNear (const HepRotation & r, 101 bool HepBoostY::isNear (const HepRotation & r, double epsilon) const {
101 double db2 = norm2(); 102 double db2 = norm2();
102 if (db2 > epsilon*epsilon) return false; 103 if (db2 > epsilon*epsilon) return false;
103 double dr2 = r.norm2(); 104 double dr2 = r.norm2();
104 return (db2+dr2 <= epsilon*epsilon); 105 return (db2+dr2 <= epsilon*epsilon);
105 } 106 }
106 107
107 bool HepBoostY::isNear ( const HepLorentzRotat 108 bool HepBoostY::isNear ( const HepLorentzRotation & lt,
108 double epsilon ) const { 109 double epsilon ) const {
109 HepBoost b1; 110 HepBoost b1;
110 HepRotation r1; 111 HepRotation r1;
111 double db2 = distance2(b1); 112 double db2 = distance2(b1);
112 lt.decompose(b1,r1); 113 lt.decompose(b1,r1);
113 if (db2 > epsilon*epsilon) return false; 114 if (db2 > epsilon*epsilon) return false;
114 double dr2 = r1.norm2(); 115 double dr2 = r1.norm2();
115 return (db2 + dr2); 116 return (db2 + dr2);
116 } 117 }
117 118
118 // ---------- Properties: 119 // ---------- Properties:
119 120
120 void HepBoostY::rectify() { 121 void HepBoostY::rectify() {
121 // Assuming the representation of this is cl 122 // Assuming the representation of this is close to a true pure boost,
122 // but may have drifted due to round-off err 123 // but may have drifted due to round-off error from many operations,
123 // this forms an "exact" pure BoostY matrix 124 // this forms an "exact" pure BoostY matrix for again.
124 125
125 double b2 = beta_*beta_; 126 double b2 = beta_*beta_;
126 if (b2 >= 1) { 127 if (b2 >= 1) {
127 beta_ = 1.0 - 1.0e-8; // Nan-proofing 128 beta_ = 1.0 - 1.0e-8; // Nan-proofing
128 b2 = beta_*beta_; 129 b2 = beta_*beta_;
129 } 130 }
130 gamma_ = 1.0 / std::sqrt(1.0 - b2); 131 gamma_ = 1.0 / std::sqrt(1.0 - b2);
131 } 132 }
132 133
133 // ---------- Application: 134 // ---------- Application:
134 135
135 // ---------- Operations in the group of 4-Rot 136 // ---------- Operations in the group of 4-Rotations
136 137
137 HepBoostY HepBoostY::operator * (const HepBoos 138 HepBoostY HepBoostY::operator * (const HepBoostY & b) const {
138 return HepBoostY ( (beta()+b.beta()) / (1+be 139 return HepBoostY ( (beta()+b.beta()) / (1+beta()*b.beta()) );
139 } 140 }
140 HepLorentzRotation HepBoostY::operator * (cons 141 HepLorentzRotation HepBoostY::operator * (const HepBoost & b) const {
141 HepLorentzRotation me (*this); 142 HepLorentzRotation me (*this);
142 return me*b; 143 return me*b;
143 } 144 }
144 HepLorentzRotation HepBoostY::operator * (cons 145 HepLorentzRotation HepBoostY::operator * (const HepRotation & r) const {
145 HepLorentzRotation me (*this); 146 HepLorentzRotation me (*this);
146 return me*r; 147 return me*r;
147 } 148 }
148 HepLorentzRotation HepBoostY::operator * (cons 149 HepLorentzRotation HepBoostY::operator * (const HepLorentzRotation & lt) const {
149 HepLorentzRotation me (*this); 150 HepLorentzRotation me (*this);
150 return me*lt; 151 return me*lt;
151 } 152 }
152 153
153 // ---------- I/O 154 // ---------- I/O
154 155
155 std::ostream & HepBoostY::print( std::ostream 156 std::ostream & HepBoostY::print( std::ostream & os ) const {
156 os << "Boost in Y direction (beta = " << bet 157 os << "Boost in Y direction (beta = " << beta_
157 << ", gamma = " << gamma_ << ") "; 158 << ", gamma = " << gamma_ << ") ";
158 return os; 159 return os;
159 } 160 }
160 161
161 } // namespace CLHEP 162 } // namespace CLHEP
162 163