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