source: trunk/source/geometry/solids/CSG/test/testG4Sphere2.cc@ 1350

Last change on this file since 1350 was 1347, checked in by garnier, 15 years ago

geant4 tag 9.4
File size: 14.2 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26
27// $Id: testG4Sphere2.cc,v 1.6 2007/05/18 10:24:32 gcosmo Exp $
28// GEANT4 tag $Name: geant4-09-04-ref-00 $
29//
30// G4Sphere Test File
31//
32// o Basic asserts on each function +
33// awkward cases for tracking / geom algorithms
34//
35// o Add tests on dicovering bugs in G4Sphere.cc...
36//
37// History:
38// 28.03.95 P.Kent Initial version
39// 20.10.96 V.Grichine Final modifications to commit
40
41#include "G4ios.hh"
42#include <assert.h>
43#include <cmath>
44#include "globals.hh"
45#include "geomdefs.hh"
46#include "G4GeometryTolerance.hh"
47
48#include "ApproxEqual.hh"
49
50#include "G4ThreeVector.hh"
51#include "G4RotationMatrix.hh"
52#include "G4AffineTransform.hh"
53#include "G4VoxelLimits.hh"
54#include "G4Sphere.hh"
55
56//const G4double kApproxEqualTolerance = kCarTolerance;
57
58// Return true if the double check is approximately equal to target
59//
60// Process:
61//
62// Return true is difference < kApproxEqualTolerance
63
64//G4bool ApproxEqual(const G4double check,const G4double target)
65//{
66// return (std::fabs(check-target)<kApproxEqualTolerance) ? true : false ;
67//}
68
69// Return true if the 3vector check is approximately equal to target
70//G4bool ApproxEqual(const G4ThreeVector& check, const G4ThreeVector& target)
71//{
72// return (ApproxEqual(check.x(),target.x())&&
73// ApproxEqual(check.y(),target.y())&&
74// ApproxEqual(check.z(),target.z()))? true : false;
75//}
76
77///////////////////////////////////////////////////////////////////
78//
79// Dave's auxiliary function
80
81const G4String OutputInside(const EInside a)
82{
83 switch(a)
84 {
85 case kInside: return "Inside";
86 case kOutside: return "Outside";
87 case kSurface: return "Surface";
88 }
89 return "????";
90}
91
92
93
94int main(void)
95{
96 G4double kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
97
98 G4ThreeVector pzero(0,0,0),px(30,0,0),py(0,30,0),pz(0,0,30);
99 G4ThreeVector Pmx(-30,0,0),pmy(0,-30,0),pmz(0,0,-30);
100 G4ThreeVector pbigx(100,0,0),pbigy(0,100,0),pbigz(0,0,100);
101 G4ThreeVector pbigmx(-100,0,0),pbigmy(0,-100,0),pbigmz(0,0,-100);
102
103 G4ThreeVector ponrmin1(45,0,0),ponrmax1(50,0,0),ponzmax(0,0,50),
104 ponrmin2(45/std::sqrt(2.),45/std::sqrt(2.),0),
105 ponrmin3(0,0,-45),ponrminJ(0,0,-300),ponrmaxJ(0,0,-500),
106 ponrmax2(50/std::sqrt(2.),50/std::sqrt(2.),0);
107 G4ThreeVector ponphi1(48/std::sqrt(2.),-48/std::sqrt(2.),0),
108 ponphi2(48/std::sqrt(2.),48/std::sqrt(2.),0),
109 pInPhi(48*0.866,-24,0),
110 pOverPhi(-48/std::sqrt(2.),48/std::sqrt(2.),0);
111 G4ThreeVector pontheta1(0,48*std::sin(pi/4),48*std::cos(pi/4)),
112 pontheta2(0,48*std::sin(pi/4),-48*std::cos(pi/4));
113
114 G4ThreeVector ptestphi1(-100,-45/std::sqrt(2.),0),
115 ptestphi2(-100,45/std::sqrt(2.),0);
116 G4ThreeVector ptesttheta1(0,48/std::sqrt(2.),100),
117 ptesttheta2(0,48/std::sqrt(2.),-100);
118
119 G4ThreeVector vx(1,0,0),vy(0,1,0),vz(0,0,1);
120 G4ThreeVector vmx(-1,0,0),vmy(0,-1,0),vmz(0,0,-1);
121 G4ThreeVector vxy(1/std::sqrt(2.),1/std::sqrt(2.),0),
122 vmxmy(-1/std::sqrt(2.),-1/std::sqrt(2.),0);
123 G4ThreeVector vxmy(1/std::sqrt(2.),-1/std::sqrt(2.),0),
124 vmxy(-1/std::sqrt(2.),1/std::sqrt(2.),0);
125 G4ThreeVector v345exit1(-0.8,0.6,0),v345exit2(0.8,0.6,0),
126 v345exit3(0.6,0.8,0);
127 G4ThreeVector norm,*pNorm;
128 G4bool *pgoodNorm,goodNorm;
129
130 pNorm=&norm;
131 pgoodNorm=&goodNorm;
132
133
134 G4bool checkPoint( const G4Sphere& pSph, G4ThreeVector origin,
135 G4double d, G4ThreeVector dir, EInside exp);
136
137 G4Sphere SpAroundX("SpAroundX", 10.*mm, 1000.*mm,
138 -1.0*degree, 2.0*degree, 0.*degree, 180.0*degree );
139
140 G4double sinOneDeg = std::sin( 1.0 * degree );
141 G4double radOne = 100.0 * mm;
142
143 G4ThreeVector ptPhiSurfExct= G4ThreeVector( radOne * std::cos( -1.0 * degree ) ,
144 - radOne * sinOneDeg,
145 0.0 );
146 G4cout << " Starting from point " << ptPhiSurfExct << G4endl;
147 G4cout << " using direction " << vy << G4endl;
148 checkPoint( SpAroundX, ptPhiSurfExct, -radOne * kAngTolerance * 1.5,
149 vy, kOutside);
150 checkPoint( SpAroundX, ptPhiSurfExct, -radOne * kAngTolerance * 0.49,
151 vy, kSurface);
152 checkPoint( SpAroundX, ptPhiSurfExct, -radOne * kAngTolerance * 0.25,
153 vy, kSurface);
154 checkPoint( SpAroundX, ptPhiSurfExct, 0.0,
155 vy, kSurface);
156 checkPoint( SpAroundX, ptPhiSurfExct, radOne * kAngTolerance * 0.25,
157 vy, kSurface);
158 checkPoint( SpAroundX, ptPhiSurfExct, radOne * kAngTolerance * 0.49,
159 vy, kSurface);
160 checkPoint( SpAroundX, ptPhiSurfExct, radOne * kAngTolerance * 1.5,
161 vy, kInside);
162
163 // Try one that has a 'deep' negative phi section
164 // --> Vlad. Grichine test case, 30 Oct 2003
165 //
166 G4cout << G4endl << G4endl << "" << G4endl;
167 G4cout << "========================================================= " << G4endl;
168
169 G4Sphere SphDeepNeg("DeepNegPhiSphere", 10.*mm, 1000.*mm,
170 -270.0*degree, 280.0*degree, // start Phi, delta Phi
171 0.*degree, 180.0*degree ); // start Theta, delta Theta
172 G4double phiPoint = 160.0 * degree;
173 G4ThreeVector StartPt( radOne * std::cos(phiPoint), radOne * std::sin(phiPoint), 0.0);
174 G4cout << "For sphere " << SphDeepNeg.GetName() << G4endl;
175 G4cout << " Starting from point " << ptPhiSurfExct << G4endl;
176
177 checkPoint( SphDeepNeg, StartPt, 0.0, vy, kInside);
178
179 // Try the edges
180 G4ThreeVector NegEdgePt( radOne * std::cos(-270.0*degree),
181 radOne * std::sin(-270.0*degree), 0.0);
182 G4ThreeVector PosEdgePt( radOne * std::cos(10.0*degree),
183 radOne * std::sin(10.0*degree), 0.0);
184
185 G4cout << "--------------------------------------------------------" << G4endl;
186 G4cout << " New point " << NegEdgePt << " should be at Neg edge of -270.0 degrees " << G4endl;
187 checkPoint( SphDeepNeg, NegEdgePt, 0.0, -vx, kSurface);
188 checkPoint( SphDeepNeg, NegEdgePt, radOne*kAngTolerance * 0.25, -vx, kSurface);
189 checkPoint( SphDeepNeg, NegEdgePt, -radOne*kAngTolerance * 0.25, -vx, kSurface);
190 checkPoint( SphDeepNeg, NegEdgePt, radOne*kAngTolerance * 1.25, -vx, kInside);
191 checkPoint( SphDeepNeg, NegEdgePt, -radOne*kAngTolerance * 1.25, -vx, kOutside);
192
193 G4cout << "--------------------------------------------------------" << G4endl;
194 G4cout << " New point " << PosEdgePt << " should be at Pos edge of +10.0 degrees " << G4endl;
195 checkPoint( SphDeepNeg, PosEdgePt, 0.0, -vy, kSurface);
196 checkPoint( SphDeepNeg, PosEdgePt, radOne*kAngTolerance * 0.25, -vy, kSurface);
197 checkPoint( SphDeepNeg, PosEdgePt, -radOne*kAngTolerance * 0.25, -vy, kSurface);
198 checkPoint( SphDeepNeg, PosEdgePt, -radOne*kAngTolerance * 1.25, -vy, kOutside);
199 checkPoint( SphDeepNeg, PosEdgePt, radOne*kAngTolerance * 1.25, -vy, kInside);
200
201 G4double radMax= 1000.0 * mm;
202 NegEdgePt = G4ThreeVector( radMax * std::cos(-270.0*degree),
203 radMax * std::sin(-270.0*degree), 0.0);
204 G4cout << "--------------------------------------------------------" << G4endl;
205 G4cout << " New point " << NegEdgePt << " should be at RadMax / Neg edge of -270.0 degrees " << G4endl;
206 checkPoint( SphDeepNeg, NegEdgePt, 0.0, -vx, kSurface);
207 checkPoint( SphDeepNeg, NegEdgePt, radMax*kAngTolerance * 0.25, -vx, kSurface);
208 checkPoint( SphDeepNeg, NegEdgePt, -radMax*kAngTolerance * 0.25, -vx, kSurface);
209 checkPoint( SphDeepNeg, NegEdgePt, radMax*kAngTolerance * 1.25, -vx, kSurface);
210 checkPoint( SphDeepNeg, NegEdgePt, -radMax*kAngTolerance * 1.25, -vx, kOutside);
211
212 PosEdgePt= G4ThreeVector( radMax * std::cos(10.0*degree),
213 radMax * std::sin(10.0*degree), 0.0);
214
215 G4cout << "--------------------------------------------------------" << G4endl;
216 G4cout << " New point " << PosEdgePt << " should be at RadMax Pos edge of +10.0 degrees " << G4endl;
217 checkPoint( SphDeepNeg, PosEdgePt, 0.0, -vy, kSurface);
218 checkPoint( SphDeepNeg, PosEdgePt, radMax*kAngTolerance * 0.25, -vy, kSurface);
219 checkPoint( SphDeepNeg, PosEdgePt, -radMax*kAngTolerance * 0.25, -vy, kSurface);
220 checkPoint( SphDeepNeg, PosEdgePt, -radMax*kAngTolerance * 1.25, -vy, kOutside);
221 checkPoint( SphDeepNeg, PosEdgePt, radMax*kAngTolerance * 1.25, -vy, kSurface);
222
223 G4ThreeVector NormInDir = - std::cos(10.0*degree) * vy + std::sin(10.0*degree) * vx;
224
225 checkPoint( SphDeepNeg, PosEdgePt, 0.0, -vy, kSurface);
226 checkPoint( SphDeepNeg, PosEdgePt, radMax*kAngTolerance * 0.25, NormInDir, kSurface);
227 checkPoint( SphDeepNeg, PosEdgePt, -radMax*kAngTolerance * 0.25, NormInDir, kSurface);
228 checkPoint( SphDeepNeg, PosEdgePt, -radMax*kAngTolerance * 1.25, NormInDir, kOutside);
229 checkPoint( SphDeepNeg, PosEdgePt, radMax*kAngTolerance * 1.25, NormInDir, kSurface);
230
231 return 0;
232}
233
234// Given the sphere 'rSphere', a point 'origin'
235// --> --> -->
236// the function below checks that the point pnew=( origin + dist * direction )
237// - the point (p+dist*dir) is located in the part of the solid given by 'expectedInResult'
238// - and that from there, the DistanceToIn along 'dir' is not negative or Infinite
239//
240// Use cases expected:
241// - 'origin' is on/near a surface and 'direction' is pointing towards the inside of the solid
242
243G4bool
244checkPoint( const G4Sphere &rSphere,
245 G4ThreeVector origin,
246 G4double dist,
247 G4ThreeVector direction,
248 EInside expectedInResult)
249{
250 G4int verbose = 0, verboseErr= 2;
251 G4bool testAll = false ; // if false, do not call DistToIn if Outside etc.
252
253 G4double kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
254
255 G4ThreeVector newPoint;
256 G4double distIn = -1.0, distOut = -1.0;
257
258 newPoint = origin + dist * direction;
259
260 G4int oldPrecision= G4cout.precision(10);
261 // G4cout << " --- Sphere " << rSphere.GetName() << "" << G4endl;
262
263 if( verbose > 0 )
264 {
265 G4cout << G4endl;
266 if (verbose > 2 ) G4cout << " Sphere " << rSphere.GetName();
267 G4cout.precision(10);
268 if (verbose > 1 ) G4cout << " dir= " << direction;
269 G4cout << " dist= " << dist;
270 }
271 EInside inSphere= rSphere.Inside( newPoint ) ;
272 /*======*/
273 G4cout.precision(15);
274 // G4cout << " NewPoint " << newPoint << " is "
275 G4bool goodIn= (inSphere == expectedInResult);
276
277 if ( !goodIn )
278 {
279 G4cout << " ************ Unexpected Result for Inside *************** " << G4endl;
280 }
281 if ( verbose || !goodIn )
282 {
283 G4cout << " New point "
284 << " is " << OutputInside( inSphere )
285 << " vs " << OutputInside( expectedInResult )
286 << " expected." << G4endl ;
287 }
288 G4bool goodDistIn = true;
289
290 distIn = rSphere.DistanceToIn( newPoint, direction );
291 /*===========*/
292 if ( verbose ) G4cout << " DistToIn (p, dir) = " << distIn << G4endl;
293
294 if( (inSphere == kOutside) && (distIn < 0.0 ) // Cannot use 0.5*kCarTolerance for Angular tolerance!!
295 )
296 {
297 G4cout << " ********** Unexpected Result for DistanceToIn from outside ********* " << G4endl;
298 // G4cout << " It should be " << G4endl;
299 goodDistIn = false;
300 }
301 if( (inSphere == kSurface ) &&
302 ( (distIn < 0.0) || (distIn >= kInfinity )) )
303 {
304 G4cout << " ********** Unexpected Result for DistanceToIn on surface ********* " << G4endl;
305 // if ( (distIn != 0.0) )
306 // - Can check that the return value must be 0.0
307 // But in general case the direction can be away from the solid,
308 // and then a finite or kInfinity answer is correct
309 // --> must check the direction against the normal
310 // in order to perform this check in general case.
311
312 goodDistIn = false;
313 }
314 if ( verbose || !goodDistIn )
315 {
316 G4cout << " DistToIn (p, dir) = " << distIn << G4endl;
317 }
318 G4bool good= (goodIn && goodDistIn);
319
320 if ( !good )
321 {
322 // There was an error -- document the use case!
323 G4cout << " --- Sphere " << rSphere.GetName() << "" << G4endl;
324
325 if ( verboseErr > 1 )
326 {
327 G4cout << " dist= " << dist ; // << G4endl;
328 G4cout << " Direction= " << direction ; // << G4endl;
329 G4cout << " dist/kAngTolerance = " << dist / kAngTolerance << G4endl;
330 G4cout << " Original pt= " << origin << G4endl;
331 }
332 G4cout << " Actual-point= " << newPoint << G4endl;
333 G4cout << " Rho= " << G4ThreeVector(newPoint.x(), newPoint.y(), 0.).mag() << G4endl;
334 }
335 if( testAll || (inSphere!=kOutside) )
336 {
337 distOut = rSphere.DistanceToOut( newPoint, direction );
338 /*=============*/
339 if ( verbose ) G4cout << " DistToOut (p, dir) = " << distOut << G4endl;
340 }
341 G4cout.precision(oldPrecision);
342
343 return good;
344}
Note: See TracBrowser for help on using the repository browser.