| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // G4ErrorCylSurfaceTarget class implementatio 26 // G4ErrorCylSurfaceTarget class implementation
27 // 27 //
28 // Created: P.Arce, September 2004 28 // Created: P.Arce, September 2004
29 // ------------------------------------------- 29 // --------------------------------------------------------------------
30 30
31 #include "G4ErrorCylSurfaceTarget.hh" 31 #include "G4ErrorCylSurfaceTarget.hh"
32 #include "G4GeometryTolerance.hh" 32 #include "G4GeometryTolerance.hh"
33 33
34 #ifdef G4VERBOSE 34 #ifdef G4VERBOSE
35 #include "G4ErrorPropagatorData.hh" // for ver 35 #include "G4ErrorPropagatorData.hh" // for verbosity checking
36 #endif 36 #endif
37 37
38 #include "geomdefs.hh" 38 #include "geomdefs.hh"
39 #include "G4Normal3D.hh" <<
40 #include "G4Plane3D.hh" 39 #include "G4Plane3D.hh"
41 40
42 //-------------------------------------------- 41 //---------------------------------------------------------------------
43 42
44 G4ErrorCylSurfaceTarget:: 43 G4ErrorCylSurfaceTarget::
45 G4ErrorCylSurfaceTarget( const G4double& radiu 44 G4ErrorCylSurfaceTarget( const G4double& radius,
46 const G4ThreeVector& 45 const G4ThreeVector& trans,
47 const G4RotationMatri 46 const G4RotationMatrix& rotm )
48 : fradius(radius) 47 : fradius(radius)
49 { 48 {
50 theType = G4ErrorTarget_CylindricalSurface; 49 theType = G4ErrorTarget_CylindricalSurface;
51 50
52 ftransform = G4AffineTransform( rotm.inverse 51 ftransform = G4AffineTransform( rotm.inverse(), -trans );
53 #ifdef G4VERBOSE 52 #ifdef G4VERBOSE
54 if(G4ErrorPropagatorData::verbose() >= 2 ) 53 if(G4ErrorPropagatorData::verbose() >= 2 )
55 { 54 {
56 Dump( " $$$ creating G4ErrorCylSurfaceTarg 55 Dump( " $$$ creating G4ErrorCylSurfaceTarget ");
57 } 56 }
58 #endif 57 #endif
59 } 58 }
60 59
61 //-------------------------------------------- 60 //---------------------------------------------------------------------
62 61
63 G4ErrorCylSurfaceTarget:: 62 G4ErrorCylSurfaceTarget::
64 G4ErrorCylSurfaceTarget( const G4double& radiu 63 G4ErrorCylSurfaceTarget( const G4double& radius,
65 const G4AffineTransfo 64 const G4AffineTransform& trans )
66 : fradius(radius), ftransform(trans.Inverse( 65 : fradius(radius), ftransform(trans.Inverse())
67 { 66 {
68 theType = G4ErrorTarget_CylindricalSurface; 67 theType = G4ErrorTarget_CylindricalSurface;
69 68
70 #ifdef G4VERBOSE 69 #ifdef G4VERBOSE
71 if(G4ErrorPropagatorData::verbose() >= 2 ) 70 if(G4ErrorPropagatorData::verbose() >= 2 )
72 { 71 {
73 Dump( " $$$ creating G4ErrorCylSurfaceTarg 72 Dump( " $$$ creating G4ErrorCylSurfaceTarget ");
74 } 73 }
75 #endif 74 #endif
>> 75 }
>> 76
>> 77 //---------------------------------------------------------------------
>> 78
>> 79 G4ErrorCylSurfaceTarget::~G4ErrorCylSurfaceTarget()
>> 80 {
76 } 81 }
77 82
78 //-------------------------------------------- 83 //---------------------------------------------------------------------
79 84
80 G4double G4ErrorCylSurfaceTarget:: 85 G4double G4ErrorCylSurfaceTarget::
81 GetDistanceFromPoint( const G4ThreeVector& poi 86 GetDistanceFromPoint( const G4ThreeVector& point,
82 const G4ThreeVector& dir 87 const G4ThreeVector& dir ) const
83 { 88 {
84 if( dir.mag() == 0. ) 89 if( dir.mag() == 0. )
85 { 90 {
86 G4Exception("G4ErrorCylSurfaceTarget::GetD 91 G4Exception("G4ErrorCylSurfaceTarget::GetDistanceFromPoint()",
87 "GeomMgt0003", FatalException, 92 "GeomMgt0003", FatalException, "Direction is zero !");
88 } 93 }
89 94
90 //----- Get intersection point 95 //----- Get intersection point
91 G4ThreeVector localPoint = ftransform.Transf 96 G4ThreeVector localPoint = ftransform.TransformPoint( point );
92 G4ThreeVector localDir = ftransform.Transfor 97 G4ThreeVector localDir = ftransform.TransformAxis( dir );
93 G4ThreeVector inters = IntersectLocal(localP 98 G4ThreeVector inters = IntersectLocal(localPoint, localDir);
94 99
95 G4double dist = (localPoint-inters).mag(); 100 G4double dist = (localPoint-inters).mag();
96 101
97 #ifdef G4VERBOSE 102 #ifdef G4VERBOSE
98 if(G4ErrorPropagatorData::verbose() >= 3 ) 103 if(G4ErrorPropagatorData::verbose() >= 3 )
99 { 104 {
100 G4cout << " G4ErrorCylSurfaceTarget::GetDi 105 G4cout << " G4ErrorCylSurfaceTarget::GetDistanceFromPoint():" << G4endl
101 << " Global point " << point << " d 106 << " Global point " << point << " dir " << dir << G4endl
102 << " Intersection " << inters << G4 107 << " Intersection " << inters << G4endl
103 << " Distance " << dist << G4endl; 108 << " Distance " << dist << G4endl;
104 Dump( " CylSurface: " ); 109 Dump( " CylSurface: " );
105 } 110 }
106 #endif 111 #endif
107 112
108 return dist; 113 return dist;
109 } 114 }
110 115
111 116
112 //-------------------------------------------- 117 //---------------------------------------------------------------------
113 118
114 G4double G4ErrorCylSurfaceTarget:: 119 G4double G4ErrorCylSurfaceTarget::
115 GetDistanceFromPoint( const G4ThreeVector& poi 120 GetDistanceFromPoint( const G4ThreeVector& point ) const
116 { 121 {
117 G4ThreeVector localPoint = ftransform.Transf 122 G4ThreeVector localPoint = ftransform.TransformPoint( point );
118 123
119 #ifdef G4VERBOSE 124 #ifdef G4VERBOSE
120 if(G4ErrorPropagatorData::verbose() >= 3 ) 125 if(G4ErrorPropagatorData::verbose() >= 3 )
121 { 126 {
122 G4cout << " G4ErrorCylSurfaceTarget::GetDi 127 G4cout << " G4ErrorCylSurfaceTarget::GetDistanceFromPoint:" << G4endl
123 << " Global point " << point << G4e 128 << " Global point " << point << G4endl
124 << " Distance " << fradius - localP 129 << " Distance " << fradius - localPoint.perp() << G4endl;
125 Dump( " CylSurface: " ); 130 Dump( " CylSurface: " );
126 } 131 }
127 #endif 132 #endif
128 133
129 return fradius - localPoint.perp(); 134 return fradius - localPoint.perp();
130 } 135 }
131 136
132 //-------------------------------------------- 137 //---------------------------------------------------------------------
133 138
134 G4ThreeVector G4ErrorCylSurfaceTarget:: 139 G4ThreeVector G4ErrorCylSurfaceTarget::
135 IntersectLocal( const G4ThreeVector& localPoin 140 IntersectLocal( const G4ThreeVector& localPoint,
136 const G4ThreeVector& localDir 141 const G4ThreeVector& localDir ) const
137 { 142 {
138 G4double eqa = localDir.x()*localDir.x()+loc 143 G4double eqa = localDir.x()*localDir.x()+localDir.y()*localDir.y();
139 G4double eqb = 2*(localPoint.x()*localDir.x( 144 G4double eqb = 2*(localPoint.x()*localDir.x()+localPoint.y()*localDir.y());
140 G4double eqc = -fradius*fradius+localPoint.x 145 G4double eqc = -fradius*fradius+localPoint.x()*localPoint.x()
141 +localPoint.y()*localPoint.y( 146 +localPoint.y()*localPoint.y();
142 G4int inside = (localPoint.perp() > fradius) 147 G4int inside = (localPoint.perp() > fradius) ? -1 : 1;
143 G4double lambda; 148 G4double lambda;
144 149
145 if( eqa*inside > 0. ) 150 if( eqa*inside > 0. )
146 { 151 {
147 lambda = (-eqb + std::sqrt(eqb*eqb-4*eqa*e 152 lambda = (-eqb + std::sqrt(eqb*eqb-4*eqa*eqc) ) / (2.*eqa);
148 } 153 }
149 else if( eqa*inside < 0. ) 154 else if( eqa*inside < 0. )
150 { 155 {
151 lambda = (-eqb - std::sqrt(eqb*eqb-4*eqa*e 156 lambda = (-eqb - std::sqrt(eqb*eqb-4*eqa*eqc) ) / (2.*eqa);
152 } 157 }
153 else 158 else
154 { 159 {
155 if( eqb != 0. ) 160 if( eqb != 0. )
156 { 161 {
157 lambda = -eqc/eqb; 162 lambda = -eqc/eqb;
158 } 163 }
159 else 164 else
160 { 165 {
161 std::ostringstream message; 166 std::ostringstream message;
162 message << "Intersection not possible !" 167 message << "Intersection not possible !" << G4endl
163 << " Point: " << localP 168 << " Point: " << localPoint << ", direction: "
164 << localDir; 169 << localDir;
165 Dump( " CylSurface: " ); 170 Dump( " CylSurface: " );
166 G4Exception("G4ErrorCylSurfaceTarget::In 171 G4Exception("G4ErrorCylSurfaceTarget::IntersectLocal()",
167 "GeomMgt1002", JustWarning, 172 "GeomMgt1002", JustWarning, message);
168 lambda = kInfinity; 173 lambda = kInfinity;
169 } 174 }
170 } 175 }
171 176
172 G4ThreeVector inters = localPoint + lambda*l 177 G4ThreeVector inters = localPoint + lambda*localDir/localDir.mag();
173 178
174 #ifdef G4VERBOSE 179 #ifdef G4VERBOSE
175 if(G4ErrorPropagatorData::verbose() >= 4 ) { 180 if(G4ErrorPropagatorData::verbose() >= 4 ) {
176 G4cout << " G4ErrorCylSurfaceTarget::Inter 181 G4cout << " G4ErrorCylSurfaceTarget::IntersectLocal " << inters << " "
177 << inters.perp() << " localPoint " 182 << inters.perp() << " localPoint " << localPoint << " localDir "
178 << localDir << G4endl; 183 << localDir << G4endl;
179 } 184 }
180 #endif 185 #endif
181 186
182 return inters; 187 return inters;
183 } 188 }
184 189
185 //-------------------------------------------- 190 //---------------------------------------------------------------------
186 191
187 G4Plane3D G4ErrorCylSurfaceTarget:: 192 G4Plane3D G4ErrorCylSurfaceTarget::
188 GetTangentPlane( const G4ThreeVector& point ) 193 GetTangentPlane( const G4ThreeVector& point ) const
189 { 194 {
190 G4ThreeVector localPoint = ftransform.Transf 195 G4ThreeVector localPoint = ftransform.TransformPoint( point );
191 196
192 // check that point is at cylinder surface 197 // check that point is at cylinder surface
193 // 198 //
194 if( std::fabs( localPoint.perp() - fradius ) 199 if( std::fabs( localPoint.perp() - fradius )
195 > 1000.*G4GeometryTolerance::GetInstance 200 > 1000.*G4GeometryTolerance::GetInstance()->GetSurfaceTolerance() )
196 { 201 {
197 std::ostringstream message; 202 std::ostringstream message;
198 message << "Local point not at surface !" 203 message << "Local point not at surface !" << G4endl
199 << " Point: " << point << 204 << " Point: " << point << ", local: " << localPoint
200 << G4endl 205 << G4endl
201 << " is not at surface, b 206 << " is not at surface, but far away by: "
202 << localPoint.perp() - fradius << 207 << localPoint.perp() - fradius << " !";
203 G4Exception("G4ErrorCylSurfaceTarget::GetT 208 G4Exception("G4ErrorCylSurfaceTarget::GetTangentPlane()",
204 "GeomMgt1002", JustWarning, me 209 "GeomMgt1002", JustWarning, message);
205 } 210 }
206 211
207 G4Normal3D normal = localPoint - ftransform. 212 G4Normal3D normal = localPoint - ftransform.NetTranslation();
208 213
209 return G4Plane3D( normal, point ); 214 return G4Plane3D( normal, point );
210 } 215 }
211 216
212 217
213 //-------------------------------------------- 218 //---------------------------------------------------------------------
214 219
215 void G4ErrorCylSurfaceTarget::Dump( const G4St 220 void G4ErrorCylSurfaceTarget::Dump( const G4String& msg ) const
216 { 221 {
217 G4cout << msg << " radius " << fradius 222 G4cout << msg << " radius " << fradius
218 << " centre " << ftransform.Ne 223 << " centre " << ftransform.NetTranslation()
219 << " rotation " << ftransform. 224 << " rotation " << ftransform.NetRotation() << G4endl;
220 } 225 }
221 226