| 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 //
27 // << 27 // $Id: G4ErrorCylSurfaceTarget.cc 66356 2012-12-18 09:02:32Z gcosmo $
28 // Created: P.Arce, September 2004 << 28 //
>> 29 //
>> 30 // --------------------------------------------------------------------
>> 31 // GEANT 4 class implementation file
29 // ------------------------------------------- 32 // --------------------------------------------------------------------
30 33
31 #include "G4ErrorCylSurfaceTarget.hh" 34 #include "G4ErrorCylSurfaceTarget.hh"
32 #include "G4GeometryTolerance.hh" 35 #include "G4GeometryTolerance.hh"
33 36
34 #ifdef G4VERBOSE 37 #ifdef G4VERBOSE
35 #include "G4ErrorPropagatorData.hh" // for ver << 38 #include "G4ErrorPropagatorData.hh" //for verbosity checking
36 #endif 39 #endif
37 40
38 #include "geomdefs.hh" 41 #include "geomdefs.hh"
39 #include "G4Normal3D.hh" <<
40 #include "G4Plane3D.hh" 42 #include "G4Plane3D.hh"
41 43
42 //-------------------------------------------- 44 //---------------------------------------------------------------------
43 45
44 G4ErrorCylSurfaceTarget:: 46 G4ErrorCylSurfaceTarget::
45 G4ErrorCylSurfaceTarget( const G4double& radiu 47 G4ErrorCylSurfaceTarget( const G4double& radius,
46 const G4ThreeVector& 48 const G4ThreeVector& trans,
47 const G4RotationMatri 49 const G4RotationMatrix& rotm )
48 : fradius(radius) 50 : fradius(radius)
49 { 51 {
50 theType = G4ErrorTarget_CylindricalSurface; 52 theType = G4ErrorTarget_CylindricalSurface;
51 53
52 ftransform = G4AffineTransform( rotm.inverse 54 ftransform = G4AffineTransform( rotm.inverse(), -trans );
53 #ifdef G4VERBOSE 55 #ifdef G4VERBOSE
54 if(G4ErrorPropagatorData::verbose() >= 2 ) << 56 if(G4ErrorPropagatorData::verbose() >= 2 ) {
55 { <<
56 Dump( " $$$ creating G4ErrorCylSurfaceTarg 57 Dump( " $$$ creating G4ErrorCylSurfaceTarget ");
57 } 58 }
58 #endif 59 #endif
59 } 60 }
60 61
61 //-------------------------------------------- 62 //---------------------------------------------------------------------
62 63
63 G4ErrorCylSurfaceTarget:: 64 G4ErrorCylSurfaceTarget::
64 G4ErrorCylSurfaceTarget( const G4double& radiu 65 G4ErrorCylSurfaceTarget( const G4double& radius,
65 const G4AffineTransfo 66 const G4AffineTransform& trans )
66 : fradius(radius), ftransform(trans.Inverse( 67 : fradius(radius), ftransform(trans.Inverse())
67 { 68 {
68 theType = G4ErrorTarget_CylindricalSurface; 69 theType = G4ErrorTarget_CylindricalSurface;
69 70
70 #ifdef G4VERBOSE 71 #ifdef G4VERBOSE
71 if(G4ErrorPropagatorData::verbose() >= 2 ) 72 if(G4ErrorPropagatorData::verbose() >= 2 )
72 { 73 {
73 Dump( " $$$ creating G4ErrorCylSurfaceTarg 74 Dump( " $$$ creating G4ErrorCylSurfaceTarget ");
74 } 75 }
75 #endif 76 #endif
>> 77 }
>> 78
>> 79 //---------------------------------------------------------------------
>> 80
>> 81 G4ErrorCylSurfaceTarget::~G4ErrorCylSurfaceTarget()
>> 82 {
76 } 83 }
77 84
78 //-------------------------------------------- 85 //---------------------------------------------------------------------
79 86
80 G4double G4ErrorCylSurfaceTarget:: 87 G4double G4ErrorCylSurfaceTarget::
81 GetDistanceFromPoint( const G4ThreeVector& poi 88 GetDistanceFromPoint( const G4ThreeVector& point,
82 const G4ThreeVector& dir 89 const G4ThreeVector& dir ) const
83 { 90 {
84 if( dir.mag() == 0. ) 91 if( dir.mag() == 0. )
85 { 92 {
86 G4Exception("G4ErrorCylSurfaceTarget::GetD 93 G4Exception("G4ErrorCylSurfaceTarget::GetDistanceFromPoint()",
87 "GeomMgt0003", FatalException, 94 "GeomMgt0003", FatalException, "Direction is zero !");
88 } 95 }
89 96
90 //----- Get intersection point 97 //----- Get intersection point
91 G4ThreeVector localPoint = ftransform.Transf 98 G4ThreeVector localPoint = ftransform.TransformPoint( point );
92 G4ThreeVector localDir = ftransform.Transfor 99 G4ThreeVector localDir = ftransform.TransformAxis( dir );
93 G4ThreeVector inters = IntersectLocal(localP 100 G4ThreeVector inters = IntersectLocal(localPoint, localDir);
94 101
95 G4double dist = (localPoint-inters).mag(); 102 G4double dist = (localPoint-inters).mag();
96 103
97 #ifdef G4VERBOSE 104 #ifdef G4VERBOSE
98 if(G4ErrorPropagatorData::verbose() >= 3 ) 105 if(G4ErrorPropagatorData::verbose() >= 3 )
99 { 106 {
100 G4cout << " G4ErrorCylSurfaceTarget::GetDi 107 G4cout << " G4ErrorCylSurfaceTarget::GetDistanceFromPoint():" << G4endl
101 << " Global point " << point << " d 108 << " Global point " << point << " dir " << dir << G4endl
102 << " Intersection " << inters << G4 109 << " Intersection " << inters << G4endl
103 << " Distance " << dist << G4endl; 110 << " Distance " << dist << G4endl;
104 Dump( " CylSurface: " ); 111 Dump( " CylSurface: " );
105 } 112 }
106 #endif 113 #endif
107 114
108 return dist; 115 return dist;
109 } 116 }
110 117
111 118
112 //-------------------------------------------- 119 //---------------------------------------------------------------------
113 120
114 G4double G4ErrorCylSurfaceTarget:: 121 G4double G4ErrorCylSurfaceTarget::
115 GetDistanceFromPoint( const G4ThreeVector& poi 122 GetDistanceFromPoint( const G4ThreeVector& point ) const
116 { 123 {
117 G4ThreeVector localPoint = ftransform.Transf 124 G4ThreeVector localPoint = ftransform.TransformPoint( point );
118 125
119 #ifdef G4VERBOSE 126 #ifdef G4VERBOSE
120 if(G4ErrorPropagatorData::verbose() >= 3 ) 127 if(G4ErrorPropagatorData::verbose() >= 3 )
121 { 128 {
122 G4cout << " G4ErrorCylSurfaceTarget::GetDi 129 G4cout << " G4ErrorCylSurfaceTarget::GetDistanceFromPoint:" << G4endl
123 << " Global point " << point << G4e 130 << " Global point " << point << G4endl
124 << " Distance " << fradius - localP 131 << " Distance " << fradius - localPoint.perp() << G4endl;
125 Dump( " CylSurface: " ); 132 Dump( " CylSurface: " );
126 } 133 }
127 #endif 134 #endif
128 135
129 return fradius - localPoint.perp(); 136 return fradius - localPoint.perp();
130 } 137 }
131 138
132 //-------------------------------------------- 139 //---------------------------------------------------------------------
133 140
134 G4ThreeVector G4ErrorCylSurfaceTarget:: 141 G4ThreeVector G4ErrorCylSurfaceTarget::
135 IntersectLocal( const G4ThreeVector& localPoin 142 IntersectLocal( const G4ThreeVector& localPoint,
136 const G4ThreeVector& localDir 143 const G4ThreeVector& localDir ) const
137 { 144 {
138 G4double eqa = localDir.x()*localDir.x()+loc 145 G4double eqa = localDir.x()*localDir.x()+localDir.y()*localDir.y();
139 G4double eqb = 2*(localPoint.x()*localDir.x( 146 G4double eqb = 2*(localPoint.x()*localDir.x()+localPoint.y()*localDir.y());
140 G4double eqc = -fradius*fradius+localPoint.x 147 G4double eqc = -fradius*fradius+localPoint.x()*localPoint.x()
141 +localPoint.y()*localPoint.y( 148 +localPoint.y()*localPoint.y();
142 G4int inside = (localPoint.perp() > fradius) 149 G4int inside = (localPoint.perp() > fradius) ? -1 : 1;
143 G4double lambda; 150 G4double lambda;
144 151
145 if( eqa*inside > 0. ) 152 if( eqa*inside > 0. )
146 { 153 {
147 lambda = (-eqb + std::sqrt(eqb*eqb-4*eqa*e 154 lambda = (-eqb + std::sqrt(eqb*eqb-4*eqa*eqc) ) / (2.*eqa);
148 } 155 }
149 else if( eqa*inside < 0. ) 156 else if( eqa*inside < 0. )
150 { 157 {
151 lambda = (-eqb - std::sqrt(eqb*eqb-4*eqa*e 158 lambda = (-eqb - std::sqrt(eqb*eqb-4*eqa*eqc) ) / (2.*eqa);
152 } 159 }
153 else 160 else
154 { 161 {
155 if( eqb != 0. ) 162 if( eqb != 0. )
156 { 163 {
157 lambda = -eqc/eqb; 164 lambda = -eqc/eqb;
158 } 165 }
159 else 166 else
160 { 167 {
161 std::ostringstream message; 168 std::ostringstream message;
162 message << "Intersection not possible !" 169 message << "Intersection not possible !" << G4endl
163 << " Point: " << localP 170 << " Point: " << localPoint << ", direction: "
164 << localDir; 171 << localDir;
165 Dump( " CylSurface: " ); 172 Dump( " CylSurface: " );
166 G4Exception("G4ErrorCylSurfaceTarget::In 173 G4Exception("G4ErrorCylSurfaceTarget::IntersectLocal()",
167 "GeomMgt1002", JustWarning, 174 "GeomMgt1002", JustWarning, message);
168 lambda = kInfinity; 175 lambda = kInfinity;
169 } 176 }
170 } 177 }
171 178
172 G4ThreeVector inters = localPoint + lambda*l 179 G4ThreeVector inters = localPoint + lambda*localDir/localDir.mag();
173 180
174 #ifdef G4VERBOSE 181 #ifdef G4VERBOSE
175 if(G4ErrorPropagatorData::verbose() >= 4 ) { 182 if(G4ErrorPropagatorData::verbose() >= 4 ) {
176 G4cout << " G4ErrorCylSurfaceTarget::Inter 183 G4cout << " G4ErrorCylSurfaceTarget::IntersectLocal " << inters << " "
177 << inters.perp() << " localPoint " 184 << inters.perp() << " localPoint " << localPoint << " localDir "
178 << localDir << G4endl; 185 << localDir << G4endl;
179 } 186 }
180 #endif 187 #endif
181 188
182 return inters; 189 return inters;
183 } 190 }
184 191
185 //-------------------------------------------- 192 //---------------------------------------------------------------------
186 193
187 G4Plane3D G4ErrorCylSurfaceTarget:: 194 G4Plane3D G4ErrorCylSurfaceTarget::
188 GetTangentPlane( const G4ThreeVector& point ) 195 GetTangentPlane( const G4ThreeVector& point ) const
189 { 196 {
190 G4ThreeVector localPoint = ftransform.Transf 197 G4ThreeVector localPoint = ftransform.TransformPoint( point );
191 198
192 // check that point is at cylinder surface 199 // check that point is at cylinder surface
193 // 200 //
194 if( std::fabs( localPoint.perp() - fradius ) 201 if( std::fabs( localPoint.perp() - fradius )
195 > 1000.*G4GeometryTolerance::GetInstance 202 > 1000.*G4GeometryTolerance::GetInstance()->GetSurfaceTolerance() )
196 { 203 {
197 std::ostringstream message; 204 std::ostringstream message;
198 message << "Local point not at surface !" 205 message << "Local point not at surface !" << G4endl
199 << " Point: " << point << 206 << " Point: " << point << ", local: " << localPoint
200 << G4endl 207 << G4endl
201 << " is not at surface, b 208 << " is not at surface, but far away by: "
202 << localPoint.perp() - fradius << 209 << localPoint.perp() - fradius << " !";
203 G4Exception("G4ErrorCylSurfaceTarget::GetT 210 G4Exception("G4ErrorCylSurfaceTarget::GetTangentPlane()",
204 "GeomMgt1002", JustWarning, me 211 "GeomMgt1002", JustWarning, message);
205 } 212 }
206 213
207 G4Normal3D normal = localPoint - ftransform. 214 G4Normal3D normal = localPoint - ftransform.NetTranslation();
208 215
209 return G4Plane3D( normal, point ); 216 return G4Plane3D( normal, point );
210 } 217 }
211 218
212 219
213 //-------------------------------------------- 220 //---------------------------------------------------------------------
214 221
215 void G4ErrorCylSurfaceTarget::Dump( const G4St 222 void G4ErrorCylSurfaceTarget::Dump( const G4String& msg ) const
216 { 223 {
217 G4cout << msg << " radius " << fradius 224 G4cout << msg << " radius " << fradius
218 << " centre " << ftransform.Ne 225 << " centre " << ftransform.NetTranslation()
219 << " rotation " << ftransform. 226 << " rotation " << ftransform.NetRotation() << G4endl;
220 } 227 }
221 228