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