source: trunk/source/geometry/solids/CSG/test/testG4Sphere.cc@ 1323

Last change on this file since 1323 was 1316, checked in by garnier, 15 years ago

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
File size: 47.1 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: testG4Sphere.cc,v 1.29 2009/05/14 13:22:44 tnikitin Exp $
28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
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
47#include "ApproxEqual.hh"
48#include "G4GeometryTolerance.hh"
49
50#include "G4ThreeVector.hh"
51#include "G4RotationMatrix.hh"
52#include "G4AffineTransform.hh"
53#include "G4VoxelLimits.hh"
54#include "G4Sphere.hh"
55#include "Randomize.hh"
56
57//const G4double kApproxEqualTolerance = kCarTolerance;
58
59// Return true if the double check is approximately equal to target
60//
61// Process:
62//
63// Return true is difference < kApproxEqualTolerance
64
65//G4bool ApproxEqual(const G4double check,const G4double target)
66//{
67// return (std::fabs(check-target)<kApproxEqualTolerance) ? true : false ;
68//}
69
70// Return true if the 3vector check is approximately equal to target
71//G4bool ApproxEqual(const G4ThreeVector& check, const G4ThreeVector& target)
72//{
73// return (ApproxEqual(check.x(),target.x())&&
74// ApproxEqual(check.y(),target.y())&&
75// ApproxEqual(check.z(),target.z()))? true : false;
76//}
77
78///////////////////////////////////////////////////////////////////
79//
80// Dave's auxiliary function
81
82const G4String OutputInside(const EInside a)
83{
84 switch(a)
85 {
86 case kInside: return "Inside";
87 case kOutside: return "Outside";
88 case kSurface: return "Surface";
89 }
90 return "????";
91}
92
93
94
95int main(void)
96{
97 G4double kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
98
99 G4int i, iMax;
100 G4double Dist, vol, volCheck;
101 G4double phi, cosTheta, sinTheta, rMax, rRand, zP;
102
103 G4ThreeVector pzero(0,0,0),px(30,0,0),py(0,30,0),pz(0,0,30);
104 G4ThreeVector pmx(-30,0,0),pmy(0,-30,0),pmz(0,0,-30);
105 G4ThreeVector pbigx(100,0,0),pbigy(0,100,0),pbigz(0,0,100);
106 G4ThreeVector pbigmx(-100,0,0),pbigmy(0,-100,0),pbigmz(0,0,-100);
107
108 G4ThreeVector ponrmin1(45,0,0),ponrmax1(50,0,0),ponzmax(0,0,50),
109 ponrmin2(45/std::sqrt(2.),45/std::sqrt(2.),0),
110 ponrmin3(0,0,-45),ponrminJ(0,0,-300),ponrmaxJ(0,0,-500),
111 ponrmax2(50/std::sqrt(2.),50/std::sqrt(2.),0);
112 G4ThreeVector ponphi1(48/std::sqrt(2.),-48/std::sqrt(2.),0),
113 ponphi2(48/std::sqrt(2.),48/std::sqrt(2.),0),
114 pInPhi(48*0.866,-24,0),
115 pOverPhi(-48/std::sqrt(2.),48/std::sqrt(2.),0);
116 G4ThreeVector pontheta1(0,48*std::sin(pi/4),48*std::cos(pi/4)),
117 pontheta2(0,48*std::sin(pi/4),-48*std::cos(pi/4));
118
119 G4ThreeVector ptestphi1(-100,-45/std::sqrt(2.),0),
120 ptestphi2(-100,45/std::sqrt(2.),0);
121
122 G4ThreeVector ptesttheta1(0,48/std::sqrt(2.),100),
123 ptesttheta2(0,48/std::sqrt(2.),-100);
124
125 // Directions
126
127 G4ThreeVector vx(1,0,0),vy(0,1,0),vz(0,0,1);
128 G4ThreeVector vmx(-1,0,0),vmy(0,-1,0),vmz(0,0,-1);
129 G4ThreeVector vxy(1/std::sqrt(2.),1/std::sqrt(2.),0),vmxmy(-1/std::sqrt(2.),-1/std::sqrt(2.),0);
130 G4ThreeVector vxmy(1/std::sqrt(2.),-1/std::sqrt(2.),0),vmxy(-1/std::sqrt(2.),1/std::sqrt(2.),0);
131 G4ThreeVector vxmz(1/std::sqrt(2.),0.,-1/std::sqrt(2.)),vymz(0.,1/std::sqrt(2.),-1/std::sqrt(2.));
132 G4ThreeVector vmxmz(-1/std::sqrt(2.),0.,-1/std::sqrt(2.)),vmxz(-1/std::sqrt(2.),0.,1/std::sqrt(2.));
133
134 G4ThreeVector v345exit1(-0.8,0.6,0),v345exit2(0.8,0.6,0),
135 v345exit3(0.6,0.8,0);
136
137
138
139
140
141 G4ThreeVector pRand, vRand, norm, *pNorm;
142 G4bool *pgoodNorm,goodNorm,calcNorm=true;
143
144 pNorm=&norm;
145 pgoodNorm=&goodNorm;
146
147 G4Sphere s1("Solid G4Sphere",0,50,0,twopi,0,pi);
148 G4Sphere sn1("sn1",0,50,halfpi,3.*halfpi,0,pi);
149
150 // Theta sections
151
152 G4Sphere sn11("sn11",0,50,0,twopi,0.,halfpi);
153 G4Sphere sn12("sn12",0,50,0,twopi,0.,0.25*pi);
154 G4Sphere sn13("sn13",0,50,0,twopi,0.75*pi,0.25*pi);
155 G4Sphere sn14("sn14",0,50,0,twopi,0.25*pi,0.75*pi);
156 G4Sphere sn15("sn15",0,50,0,twopi,89.*deg,91.*deg);
157 G4Sphere sn16("sn16",0,50,0,twopi,0.*deg,89.*deg);
158 G4Sphere sn17("sn17",0,50,0,twopi,91.*deg,89.*deg);
159
160 G4Sphere s2("Spherical Shell",45,50,0,twopi,0,pi);
161 G4Sphere sn2("sn2",45,50,halfpi,halfpi,0,pi);
162 G4Sphere sn22("sn22",0,50,halfpi,halfpi,0,pi);
163
164
165
166 G4Sphere s3("Band (theta segment)",45,50,0,twopi,pi/4,halfpi);
167 G4Sphere s32("Band (theta segment2)",45,50,0,twopi,0,pi/4);
168 G4Sphere s33("Band (theta segment1)",45,50,0,twopi,pi*3/4,pi/4);
169 G4Sphere s34("Band (theta segment)",4,50,0,twopi,pi/4,halfpi);
170 G4Sphere s4("Band (phi segment)",45,50,-pi/4,halfpi,0,twopi);
171 // G4cout<<"s4.fSPhi = "<<s4.GetSPhi()<<G4endl;
172 G4Sphere s41("Band (phi segment)",5,50,-pi,3.*pi/2.,0,twopi);
173 G4Sphere s42("Band (phi segment)",5,50,-pi/2,3.*pi/2.,0,twopi);
174 G4Sphere s5("Patch (phi/theta seg)",45,50,-pi/4,halfpi,pi/4,halfpi);
175
176 G4Sphere s6("John example",300,500,0,5.76,0,pi) ;
177 G4Sphere s7("sphere7",1400.,1550.,0.022321428571428572,0.014642857142857141,
178 1.5631177553663251,0.014642857142857141 );
179 G4Sphere s8("sphere",278.746*mm, 280.0*mm, 0.0*degree, 360.0*degree,
180 0.0*degree, 90.0*degree);
181
182 G4Sphere b216("b216", 1400.0, 1550.0,
183 0.022321428571428572,
184 0.014642857142857141,
185 1.578117755366325,
186 0.014642857142867141);
187
188
189
190 G4Sphere s9("s9",0*mm,410*mm,0*degree,360*degree,90*degree,90*degree);
191
192 G4Sphere b402("b402", 475*mm, 480*mm,
193 0*degree,360*degree,17.8*degree,144.4*degree);
194
195 G4Sphere b1046("b1046", 4750*km, 4800*km,
196 0*degree,360*degree,0*degree,90*degree);
197
198
199
200G4ThreeVector p402(471.7356120367253*mm, 51.95081450791341*mm, 5.938043020529463*mm);
201G4ThreeVector v402(-0.519985502840818, 0.2521089719986221, 0.8161226274728446);
202
203
204
205G4ThreeVector p216(1549.9518578505142,1.2195415370970153,-12.155289555510985),
206 v216(-0.61254821852534425,-0.51164551429243466,-0.60249775741147549);
207
208
209G4ThreeVector s9p(384.8213314370455*mm,
210 134.264386151667*mm,
211 -44.56026800002064*mm);
212
213G4ThreeVector s9v(-0.6542770611918751,
214 -0.0695116921641141,
215 -0.7530535517814154);
216
217 G4Sphere s10("s10",0*mm,0.018*mm,0*degree,360*degree,0*degree,180*degree);
218
219 G4ThreeVector s10p(0.01160957408065766*mm,
220 0.01308205826682229*mm,0.004293345210644617*mm);
221
222 G4Sphere s11("s11",5000000.*mm,
223 3700000000.*mm,
224 0*degree,360*degree,0*degree,180*degree);
225
226 // G4ThreeVector ps11(-1184000000.*mm,-3477212676.843337059*mm,-444000000.*mm);
227
228 // G4ThreeVector ps11( -3339032195.112830162*mm, -1480000000*mm, -592000000*mm );
229
230 G4ThreeVector ps11( -3072559844.81995153427124*mm, -1924000000*mm, -740000000*mm );
231
232 G4Sphere sAlex("sAlex",500.*mm,
233 501.*mm,
234 0*degree,360*degree,0*degree,180*degree);
235
236 G4ThreeVector psAlex(-360.4617031263808*mm,
237 -158.1198807105035*mm,308.326878333183*mm);
238
239 G4ThreeVector vsAlex(-0.7360912456240805,-0.4955800202572754,0.4610532741813497 );
240
241
242 G4Sphere sLHCB("sLHCB",8600*mm, 8606*mm,
243 -1.699135525184141*degree,
244 3.398271050368263*degree,88.52855940538514*degree,2.942881189229715*degree );
245
246 G4ThreeVector pLHCB(8600.242072535835,-255.1193517702246,-69.0010277128286);
247
248
249
250 G4Sphere b658("b658", 209.6*mm, 211.2658*mm,
251 0.0*degree, 360*degree, 0.0*degree, 90*degree);
252
253
254 G4ThreeVector p658(-35.69953348982516*mm, 198.3183279249958, 56.30959457033987 );
255 G4ThreeVector v658(-.2346058124516908,-0.9450502890785083,0.2276841318065671);
256
257
258 G4Sphere spAroundX("SpAroundX", 10.*mm, 1000.*mm, -1.0*degree,
259 2.0*degree,
260 0.*degree, 180.0*degree );
261
262 G4double radOne = 100.0*mm;
263 G4double angle = -1.0*degree - 0.25*kAngTolerance;
264 G4ThreeVector ptPhiMinus= G4ThreeVector( radOne*std::cos(angle) ,
265 radOne*std::sin(angle),
266 0.0 );
267
268 // spheres for theta cone intersections
269
270 G4Sphere s13("Band (theta segment)",5,50,0,twopi,pi/6.,halfpi);
271 G4Sphere s14("Band (theta segment)",5,50,0,twopi,pi/3.,halfpi);
272
273
274 // b. 830
275
276 G4double mainInnerRadius = 21.45 * cm ;
277 G4double mainOuterRadius = 85.0 * cm ;
278 G4double minTheta = 18.0 * degree ;
279
280
281
282
283 G4Sphere sb830( "mainSp",
284 mainInnerRadius,
285 mainOuterRadius,
286 0.0, M_PI*2,
287 minTheta,
288 M_PI - 2*minTheta);
289
290
291G4ThreeVector pb830(81.61117212,-27.77179755,196.4143423);
292 G4ThreeVector vb830(0.1644697995,0.18507236,0.9688642354);
293
294#ifdef NDEBUG
295 G4Exception("FAIL: *** Assertions must be compiled in! ***");
296#endif
297
298 G4cout.precision(20);
299
300 //////////////// Check name /////////////////////////
301
302 assert(s1.GetName()=="Solid G4Sphere");
303
304 // check cubic volume
305
306 vol = s1.GetCubicVolume();
307 volCheck = 4*pi*125000/3;
308
309 assert( ApproxEqual(vol,volCheck));
310
311 assert(ApproxEqual(s1.GetCubicVolume(),4*pi*125000/3));
312
313
314 // Some user application cases
315
316 Dist = s4.DistanceToOut(ponrmin3,vmz,calcNorm,pgoodNorm,pNorm) ;
317 // G4cout<<"s4.DistanceToOut(ponrmin3,vmz,calcNorm,pgoodNorm,pNorm) = "<<Dist
318 // <<G4endl ;
319 Dist = s6.DistanceToOut(ponrminJ,vmz,calcNorm,pgoodNorm,pNorm) ;
320 // G4cout<<"s6.DistanceToOut(ponrminJ,vmz,calcNorm,pgoodNorm,pNorm) = "<<Dist
321 // <<G4endl ;
322 Dist = s6.DistanceToOut(ponrmaxJ,vmz,calcNorm,pgoodNorm,pNorm) ;
323 // G4cout<<"s6.DistanceToOut(ponrmaxJ,vmz,calcNorm,pgoodNorm,pNorm) = "<<Dist
324 // <<G4endl ;
325
326 Dist = s7.DistanceToOut(G4ThreeVector(1399.984667238032,
327 5.9396696802500299,
328 -2.7661927818688308),
329 G4ThreeVector(0.50965504781062942,
330 -0.80958849145715217,
331 0.29123565499656401),
332 calcNorm,pgoodNorm,pNorm) ;
333 // G4cout<<"s7.DistanceToOut(shereP,sphereV,calcNorm,pgoodNorm,pNorm) = "<<Dist
334 // <<G4endl ;
335
336 Dist = s7.DistanceToIn(G4ThreeVector(1399.984667238032,
337 5.9396696802500299,
338 -2.7661927818688308),
339 G4ThreeVector(0.50965504781062942,
340 -0.80958849145715217,
341 0.29123565499656401)) ;
342 // G4cout<<"s7.DistanceToIn(shereP,sphereV,calcNorm,pgoodNorm,pNorm) = "<<Dist
343 // <<G4endl ;
344
345
346
347// Check G4Sphere::Inside
348
349 EInside inside = s7.Inside(G4ThreeVector(1399.984667238032,
350 5.9396696802500299,
351 -2.7661927818688308) ) ;
352 // G4cout<<"s7.Inside(G4ThreeVector(1399.98466 ... = "
353 // <<OutputInside(inside)<<G4endl ;
354
355 inside = s8.Inside(G4ThreeVector(-249.5020724528353*mm,
356 26.81253142743162*mm,
357 114.8988524453591*mm ) ) ;
358 // G4cout<<"s8.Inside(G4ThreeVector(-249.5020 ... = "<<OutputInside(inside)<<G4endl ;
359 inside = b216.Inside(p216);
360 // G4cout<<"b216.Inside(p216) = "<<OutputInside(inside)<<G4endl ;
361 inside = s1.Inside(pz);
362 // G4cout<<"s1.Inside(pz) = "<<OutputInside(inside)<<G4endl ;
363
364 inside = s9.Inside(s9p);
365 // G4cout<<"s9.Inside(s9p) = "<<OutputInside(inside)<<G4endl ;
366
367 inside = b402.Inside(p402);
368 // G4cout<<"p402.Inside(p402) = "<<OutputInside(inside)<<G4endl ;
369
370 inside = s10.Inside(s10p);
371 // G4cout<<"s10.Inside(s10p) = "<<OutputInside(inside)<<G4endl ;
372 // G4cout<<"p radius = "<<s10p.mag()<<G4endl ;
373
374 inside = s11.Inside(ps11);
375 // G4cout<<"s11.Inside(ps11) = "<<OutputInside(inside)<<G4endl ;
376 // G4cout<<"ps11.mag() = "<<ps11.mag()<<G4endl ;
377
378 inside = sLHCB.Inside(pLHCB);
379 // G4cout<<"sLHCB.Inside(pLHCB) = "<<OutputInside(inside)<<G4endl ;
380 // G4cout<<"pLHCB.mag() = "<<pLHCB.mag()<<G4endl ;
381
382 inside = spAroundX.Inside(ptPhiMinus);
383 // G4cout<<"spAroundX.Inside(ptPhiMinus) = "<<OutputInside(inside)<<G4endl ;
384 inside = b658.Inside(p658);
385 G4cout<<"b658.Inside(p658) = "<<OutputInside(inside)<<G4endl ;
386
387 assert(s1.Inside(pzero)==kInside);
388 // assert(s1.Inside(pz)==kInside);
389 assert(s2.Inside(pzero)==kOutside);
390 assert(s2.Inside(ponrmin2)==kSurface);
391 assert(s2.Inside(ponrmax2)==kSurface);
392 assert(s3.Inside(pontheta1)==kSurface);
393 assert(s3.Inside(pontheta2)==kSurface);
394 assert(s4.Inside(ponphi1)==kSurface);
395 assert(s4.Inside(ponphi1)==kSurface);
396 assert(s4.Inside(pOverPhi)==kOutside);
397 assert(s4.Inside(pInPhi)==kInside);
398 assert(s5.Inside(pbigz)==kOutside);
399
400 assert(s41.Inside(pmx)==kSurface);
401 assert(s42.Inside(pmx)==kSurface);
402
403// Checking G4Sphere::SurfaceNormal
404 G4double p2=1./std::sqrt(2.),p3=1./std::sqrt(3.);
405
406
407 norm=sn1.SurfaceNormal(G4ThreeVector(0.,0.,50.));
408 assert(ApproxEqual(norm,G4ThreeVector(p3,p3,p3)));
409 norm=sn1.SurfaceNormal(G4ThreeVector(0.,0.,0.));
410 assert(ApproxEqual(norm,G4ThreeVector(p2,p2,0.)));
411
412 norm=sn11.SurfaceNormal(G4ThreeVector(0.,0.,0.));
413 assert(ApproxEqual(norm,G4ThreeVector(0.,0.,-1.)));
414 norm=sn12.SurfaceNormal(G4ThreeVector(0.,0.,0.));
415 assert(ApproxEqual(norm,G4ThreeVector(0.,0.,-1.)));
416 norm=sn13.SurfaceNormal(G4ThreeVector(0.,0.,0.));
417 assert(ApproxEqual(norm,G4ThreeVector(0.,0.,1.)));
418
419 norm=sn2.SurfaceNormal(G4ThreeVector(-45.,0.,0.));
420 assert(ApproxEqual(norm,G4ThreeVector(p2,-p2,0.)));
421
422
423
424 norm=s1.SurfaceNormal(ponrmax1);
425 assert(ApproxEqual(norm,vx));
426
427// Checking G4Sphere::DistanceToOut(P)
428 Dist=s1.DistanceToOut(pzero);
429 assert(ApproxEqual(Dist,50));
430 Dist=s1.DistanceToOut(ponrmax1);
431 assert(ApproxEqual(Dist,0));
432
433// Checking G4Sphere::DistanceToOut(p,v)
434
435
436 Dist=s1.DistanceToOut(pz,vz,calcNorm,pgoodNorm,pNorm);
437 // G4cout<<"Dist=s1.DistanceToOut(pz,vz) = "<<Dist<<G4endl;
438 assert(ApproxEqual(Dist,20.));
439
440 Dist=s1.DistanceToOut(ponrmax1,vx,calcNorm,pgoodNorm,pNorm);
441 *pNorm=pNorm->unit();
442 assert(ApproxEqual(Dist,0)&&*pgoodNorm&&ApproxEqual(*pNorm,vx));
443
444 Dist=s2.DistanceToOut(ponrmin1,vx,calcNorm,pgoodNorm,pNorm);
445 *pNorm=pNorm->unit();
446 assert(ApproxEqual(Dist,5)&&*pgoodNorm&&ApproxEqual(*pNorm,vx));
447
448 Dist=s2.DistanceToOut(ponrmax2,vx,calcNorm,pgoodNorm,pNorm);
449 assert(ApproxEqual(Dist,0)&&*pgoodNorm&&ApproxEqual(*pNorm,vxy));
450
451 Dist=s1.DistanceToOut(pzero,vx,calcNorm,pgoodNorm,pNorm);
452 assert(ApproxEqual(Dist,50)&&ApproxEqual(pNorm->unit(),vx)&&*pgoodNorm);
453 Dist=s1.DistanceToOut(pzero,vmx,calcNorm,pgoodNorm,pNorm);
454 assert(ApproxEqual(Dist,50)&&ApproxEqual(pNorm->unit(),vmx)&&*pgoodNorm);
455 Dist=s1.DistanceToOut(pzero,vy,calcNorm,pgoodNorm,pNorm);
456 assert(ApproxEqual(Dist,50)&&ApproxEqual(pNorm->unit(),vy)&&*pgoodNorm);
457 Dist=s1.DistanceToOut(pzero,vmy,calcNorm,pgoodNorm,pNorm);
458 assert(ApproxEqual(Dist,50)&&ApproxEqual(pNorm->unit(),vmy)&&*pgoodNorm);
459 Dist=s1.DistanceToOut(pzero,vz,calcNorm,pgoodNorm,pNorm);
460 assert(ApproxEqual(Dist,50)&&ApproxEqual(pNorm->unit(),vz)&&*pgoodNorm);
461 Dist=s1.DistanceToOut(pzero,vmz,calcNorm,pgoodNorm,pNorm);
462 assert(ApproxEqual(Dist,50)&&ApproxEqual(pNorm->unit(),vmz)&&*pgoodNorm);
463 Dist=s1.DistanceToOut(pzero,vxy,calcNorm,pgoodNorm,pNorm);
464 assert(ApproxEqual(Dist,50)&&ApproxEqual(pNorm->unit(),vxy)&&*pgoodNorm);
465
466 Dist=s4.DistanceToOut(ponphi1,vx,calcNorm,pgoodNorm,pNorm);
467 //assert(ApproxEqual(Dist,0)&&ApproxEqual(pNorm->unit(),vmxmy)&&*pgoodNorm);
468 Dist=s4.DistanceToOut(ponphi2,vx,calcNorm,pgoodNorm,pNorm);
469 // assert(ApproxEqual(Dist,0)&&ApproxEqual(pNorm->unit(),vmxy)&&*pgoodNorm);
470 Dist=s3.DistanceToOut(pontheta1,vz,calcNorm,pgoodNorm,pNorm);
471 // assert(ApproxEqual(Dist,0)&&ApproxEqual(pNorm->unit(),vy)&&*pgoodNorm);
472 Dist=s32.DistanceToOut(pontheta1,vmz,calcNorm,pgoodNorm,pNorm);
473 //assert(ApproxEqual(Dist,0)&&ApproxEqual(pNorm->unit(),vmy)&&*pgoodNorm);
474 Dist=s32.DistanceToOut(pontheta1,vz,calcNorm,pgoodNorm,pNorm);
475 //assert(ApproxEqual(Dist,50)&&ApproxEqual(pNorm->unit(),vz)&&*pgoodNorm);
476 Dist=s1.DistanceToOut(pzero,vmz,calcNorm,pgoodNorm,pNorm);
477 assert(ApproxEqual(Dist,50)&&ApproxEqual(pNorm->unit(),vmz)&&*pgoodNorm);
478 Dist=s1.DistanceToOut(pzero,vxy,calcNorm,pgoodNorm,pNorm);
479 assert(ApproxEqual(Dist,50)&&ApproxEqual(pNorm->unit(),vxy)&&*pgoodNorm);
480
481
482 Dist=s2.DistanceToOut(ponrmin1,vxy,calcNorm,pgoodNorm,pNorm);
483 // G4cout<<"Dist=s2.DistanceToOut(pormin1,vxy) = "<<Dist<<G4endl;
484
485 Dist=s2.DistanceToOut(ponrmax1,vmx,calcNorm,pgoodNorm,pNorm);
486 // G4cout<<"Dist=s2.DistanceToOut(ponxside,vmx) = "<<Dist<<G4endl;
487 Dist=s2.DistanceToOut(ponrmax1,vmxmy,calcNorm,pgoodNorm,pNorm);
488 // G4cout<<"Dist=s2.DistanceToOut(ponxside,vmxmy) = "<<Dist<<G4endl;
489 Dist=s2.DistanceToOut(ponrmax1,vz,calcNorm,pgoodNorm,pNorm);
490 // G4cout<<"Dist=s2.DistanceToOut(ponxside,vz) = "<<Dist<<G4endl;
491
492 // Dist=s2.DistanceToOut(pbigx,vx,calcNorm,pgoodNorm,pNorm);
493 // G4cout<<"Dist=s2.DistanceToOut(pbigx,vx) = "<<Dist<<G4endl;
494 // Dist=s2.DistanceToOut(pbigx,vxy,calcNorm,pgoodNorm,pNorm);
495 // G4cout<<"Dist=s2.DistanceToOut(pbigx,vxy) = "<<Dist<<G4endl;
496 // Dist=s2.DistanceToOut(pbigx,vz,calcNorm,pgoodNorm,pNorm);
497 // G4cout<<"Dist=s2.DistanceToOut(pbigx,vz) = "<<Dist<<G4endl;
498 //Test Distance for phi section
499 Dist=sn22.DistanceToOut(G4ThreeVector(0.,49.,0.),vmy,calcNorm,pgoodNorm,pNorm);
500 assert(ApproxEqual(Dist,49.));
501 Dist=sn22.DistanceToOut(G4ThreeVector(-45.,0.,0.),vx,calcNorm,pgoodNorm,pNorm);
502 assert(ApproxEqual(Dist,45.));
503 G4cout<<"Dist from Center ="<<sn22.DistanceToOut(G4ThreeVector(0.,49.,0),G4ThreeVector(0,-1,0))<<G4endl;
504 G4cout<<"Dist from Center ="<<sn22.DistanceToOut(G4ThreeVector(-45.,0.,0),G4ThreeVector(1,0,0))<<G4endl;
505
506 //
507 Dist=b216.DistanceToOut(p216,v216,calcNorm,pgoodNorm,pNorm);
508
509 // G4cout<<"b216.DistanceToOut(p216,v216,... = "<<Dist<<G4endl;
510
511 // call from outside
512 // Dist=sAlex.DistanceToOut(psAlex,vsAlex,calcNorm,pgoodNorm,pNorm);
513 // G4cout<<"sAlex.DistanceToOut(psAlex,vsAlex,... = "<<Dist<<G4endl;
514
515
516 // Dist=s9.DistanceToIn(s9p,s9v);
517 // G4cout<<"s9.DistanceToIn(s9p,s9v,... = "<<Dist<<G4endl;
518 Dist=s13.DistanceToOut(G4ThreeVector(20.,0.,0.),vz,calcNorm,pgoodNorm,pNorm);
519 // G4cout<<"s13.DistanceToOut(G4ThreeVector(20.,0.,0.),vz... = "<<Dist<<G4endl;
520 assert(ApproxEqual(Dist,34.641016151377549));
521
522 Dist=s13.DistanceToOut(G4ThreeVector(20.,0.,0.),vmz,calcNorm,pgoodNorm,pNorm);
523 // G4cout<<"s13.DistanceToOut(G4ThreeVector(20.,0.,0.),vmz... = "<<Dist<<G4endl;
524 assert(ApproxEqual(Dist,11.547005383792508));
525
526 Dist=s14.DistanceToOut(G4ThreeVector(20.,0.,0.),vz,calcNorm,pgoodNorm,pNorm);
527 // G4cout<<"s14.DistanceToOut(G4ThreeVector(20.,0.,0.),vz... = "<<Dist<<G4endl;
528 assert(ApproxEqual(Dist,11.547005383792508));
529
530 Dist=s14.DistanceToOut(G4ThreeVector(20.,0.,0.),vmz,calcNorm,pgoodNorm,pNorm);
531 // G4cout<<"s14.DistanceToOut(G4ThreeVector(20.,0.,0.),vmz... = "<<Dist<<G4endl;
532 assert(ApproxEqual(Dist,34.641016151377549));
533
534 Dist=sb830.DistanceToOut(pb830,vb830,calcNorm,pgoodNorm,pNorm);
535 G4cout<<"sb830.DistanceToOut(pb830,vb830... = "<<Dist<<G4endl;
536 inside = sb830.Inside(pb830+Dist*vb830);
537 G4cout<<"sb830.Inside(pb830+Dist*vb830) = "<<OutputInside(inside)<<G4endl ;
538
539 Dist=sn11.DistanceToOut( G4ThreeVector(0.,0.,20.), vmz, calcNorm, pgoodNorm, pNorm);
540 // G4cout<<"sn11.DistanceToOut(G4ThreeVector(0.,0.,20.),vmz... = "<<Dist<<G4endl;
541 assert(ApproxEqual(Dist,20.));
542
543 Dist=sn11.DistanceToOut(G4ThreeVector(0.,0.,0.),vmz,calcNorm,pgoodNorm,pNorm);
544 // G4cout<<"sn11.DistanceToOut(G4ThreeVector(0.,0.,0.),vmz... = "<<Dist<<G4endl;
545 assert(ApproxEqual(Dist,0.));
546
547 Dist=sn11.DistanceToOut(G4ThreeVector(0.,0.,0.),vz,calcNorm,pgoodNorm,pNorm);
548 // G4cout<<"sn11.DistanceToOut(G4ThreeVector(0.,0.,0.),vmz... = "<<Dist<<G4endl;
549 assert(ApproxEqual(Dist,50.));
550
551 Dist=sn11.DistanceToOut(G4ThreeVector(10.,0.,0.),vmz,calcNorm,pgoodNorm,pNorm);
552 // G4cout<<"sn11.DistanceToOut(G4ThreeVector(10.,0.,0.),vmz... = "<<Dist<<G4endl;
553 assert(ApproxEqual(Dist,0.));
554
555 Dist=sn12.DistanceToOut(G4ThreeVector(0.,0.,20.),vxmz,calcNorm,pgoodNorm,pNorm);
556 // G4cout<<"sn12.DistanceToOut(G4ThreeVector(0.,0.,20.),vxmz... = "<<Dist<<G4endl;
557 assert(ApproxEqual(Dist,20./std::sqrt(2.)));
558
559 Dist=sn12.DistanceToOut(G4ThreeVector(0.,0.,0.),vmz,calcNorm,pgoodNorm,pNorm);
560 // G4cout<<"sn12.DistanceToOut(G4ThreeVector(0.,0.,0.),vmz... = "<<Dist<<G4endl;
561 assert(ApproxEqual(Dist,0.));
562
563 Dist=sn12.DistanceToOut(G4ThreeVector(10.,0.,10.),vz,calcNorm,pgoodNorm,pNorm);
564 // G4cout<<"sn12.DistanceToOut(G4ThreeVector(10.,0.,10.),vz... = "<<Dist<<G4endl;
565 assert(ApproxEqual(Dist,38.989794855663561179));
566
567
568 Dist=sn12.DistanceToOut(G4ThreeVector(10.,0.,10.),vmz,calcNorm,pgoodNorm,pNorm);
569 // G4cout<<"sn12.DistanceToOut(G4ThreeVector(10.,0.,10.),vmz... = "<<Dist<<G4endl;
570 assert(ApproxEqual(Dist,0.));
571
572
573 Dist=sn14.DistanceToOut(G4ThreeVector(10.,0.,10.),vz,calcNorm,pgoodNorm,pNorm);
574 // G4cout<<"sn14.DistanceToOut(G4ThreeVector(10.,0.,10.),vz... = "<<Dist<<G4endl;
575 assert(ApproxEqual(Dist,0.));
576
577
578 Dist=sn14.DistanceToOut(G4ThreeVector(10.,0.,5.),vz,calcNorm,pgoodNorm,pNorm);
579 // G4cout<<"sn14.DistanceToOut(G4ThreeVector(10.,0.,5.),vz... = "<<Dist<<G4endl;
580 assert(ApproxEqual(Dist,5.));
581
582
583 Dist=sn14.DistanceToOut(G4ThreeVector(10.,0.,5.),vmxz,calcNorm,pgoodNorm,pNorm);
584 // G4cout<<"sn14.DistanceToOut(G4ThreeVector(10.,0.,5.),vmxz... = "<<Dist<<G4endl;
585 assert(ApproxEqual(Dist,3.5355339059327381968));
586
587 Dist=sn14.DistanceToOut(G4ThreeVector(10.,0.,10.),vmxz,calcNorm,pgoodNorm,pNorm);
588 // G4cout<<"sn14.DistanceToOut(G4ThreeVector(10.,0.,10.),vmxz... = "<<Dist<<G4endl;
589 assert(ApproxEqual(Dist,0.));
590
591
592
593 iMax = 10000;
594 rMax = 49.;
595 zP = -1.;
596
597 for( i = 0; i < iMax; i++)
598 {
599 rRand = rMax*G4UniformRand();
600 phi = twopi*G4UniformRand();
601 pRand = G4ThreeVector( rRand*std::cos(phi), rRand*std::sin(phi), zP);
602
603 cosTheta = 0.99*G4UniformRand();
604
605 sinTheta = std::sqrt((1.-cosTheta)*(1.+cosTheta));
606 phi = twopi*G4UniformRand();
607 vRand = G4ThreeVector( sinTheta*std::cos(phi), sinTheta*std::sin(phi), cosTheta);
608
609 Dist = sn15.DistanceToOut(pRand,vRand,calcNorm,pgoodNorm,pNorm);
610
611 if( Dist*cosTheta > -2.*zP )
612 {
613 G4cout<<"sn15.DistanceToOut(pRand,vRand,..., i = "<<i<<G4endl;
614 G4cout<<pRand.x()<<", "<<pRand.y()<<", "<<pRand.z()<<G4endl;
615 G4cout<<vRand.x()<<", "<<vRand.y()<<", "<<vRand.z()<<G4endl;
616 break;
617 }
618 }
619 pRand = G4ThreeVector( -9.3576454272980740257, 10.436745988128407703, -1.);
620 vRand = G4ThreeVector( -0.064906521070262901407, 0.35249758429187183495, 0.93355910181999202102);
621 Dist = sn15.DistanceToOut(pRand,vRand,calcNorm,pgoodNorm,pNorm);
622 // G4cout<<"sn15.DistanceToOut(pRand,vRand... = "<<Dist<<G4endl;
623 assert(ApproxEqual(Dist,1.3409673565670394701));
624
625
626 zP = 1.;
627
628 for( i = 0; i < iMax; i++)
629 {
630 rRand = rMax*G4UniformRand();
631 phi = twopi*G4UniformRand();
632 pRand = G4ThreeVector( rRand*std::cos(phi), rRand*std::sin(phi), zP);
633
634 cosTheta = -0.99*G4UniformRand();
635
636 sinTheta = std::sqrt((1.-cosTheta)*(1.+cosTheta));
637 phi = twopi*G4UniformRand();
638 vRand = G4ThreeVector( sinTheta*std::cos(phi), sinTheta*std::sin(phi), cosTheta);
639
640 Dist = sn16.DistanceToOut(pRand,vRand,calcNorm,pgoodNorm,pNorm);
641
642 if( -Dist*cosTheta > 2.*zP )
643 {
644 G4cout<<"sn16.DistanceToOut(pRand,vRand,..., i = "<<i<<G4endl;
645 G4cout<<pRand.x()<<", "<<pRand.y()<<", "<<pRand.z()<<G4endl;
646 G4cout<<vRand.x()<<", "<<vRand.y()<<", "<<vRand.z()<<G4endl;
647 break;
648 }
649 }
650
651
652
653 pRand = G4ThreeVector( -18.872396210750576273, 11.573660000258405134, 1);
654 vRand = G4ThreeVector( -0.85674495247509219187, -0.51247617288646407641, -0.057933225631713866632);
655 Dist = sn16.DistanceToOut(pRand,vRand,calcNorm,pgoodNorm,pNorm);
656 // G4cout<<"sn16.DistanceToOut(pRand,vRand... = "<<Dist<<G4endl;
657 assert(ApproxEqual(Dist,8.9849541128705734394));
658
659
660
661 zP = -2.;
662
663 for( i = 0; i < iMax; i++)
664 {
665 rRand = rMax*G4UniformRand();
666 phi = twopi*G4UniformRand();
667
668 pRand = G4ThreeVector( rRand*std::cos(phi)*std::sin(91.*deg),
669 rRand*std::sin(phi)*std::sin(91.*deg),
670 rRand*std::cos(91.*deg));
671
672 cosTheta = std::cos( 180.*deg-89.*deg*G4UniformRand() );
673
674 sinTheta = std::sqrt((1.-cosTheta)*(1.+cosTheta));
675 phi = twopi*G4UniformRand();
676 vRand = G4ThreeVector( sinTheta*std::cos(phi), sinTheta*std::sin(phi), cosTheta);
677
678 Dist = sn17.DistanceToOut(pRand,vRand,calcNorm,pgoodNorm,pNorm);
679
680 if( -Dist*cosTheta > 50. )
681 {
682 G4cout<<"sn17.DistanceToOut(pRand,vRand,..., i = "<<i<<G4endl;
683 G4cout<<pRand.x()<<", "<<pRand.y()<<", "<<pRand.z()<<G4endl;
684 G4cout<<vRand.x()<<", "<<vRand.y()<<", "<<vRand.z()<<G4endl;
685 break;
686 }
687 }
688
689 pRand = G4ThreeVector(16.20075504802145, -22.42917454903122, -41.6469406430184);
690 vRand = G4ThreeVector( -0.280469198715188, 0.4463870649961534, 0.8497503261406731 );
691 norm = sn17.SurfaceNormal(pRand);
692 G4cout<<"norm = sn17.SurfaceNormal(pRand) = "<<norm.x()<<", "<<norm.y()<<", "<<norm.z()<<G4endl;
693 inside = sn17.Inside(pRand);
694 G4cout<<"sn17.Inside(pRand) = "<<OutputInside(inside)<<G4endl ;
695 Dist = sn17.DistanceToOut(pRand,vRand,calcNorm,pgoodNorm,pNorm);
696 G4cout<<"sn17.DistanceToOut(pRand,vRand,...) = "<<Dist<<G4endl;
697 // assert(ApproxEqual(Dist,8.9849541128705734394));
698
699 pRand = G4ThreeVector(2.469342694407535, -0.5746362009323557, -0.04425421860130281);
700 vRand = G4ThreeVector( -0.2513630044708909, 0.4396138160169732, -0.8622971255607673);
701 norm = sn17.SurfaceNormal(pRand);
702 G4cout<<"norm = sn17.SurfaceNormal(pRand) = "<<norm.x()<<", "<<norm.y()<<", "<<norm.z()<<G4endl;
703 inside = sn17.Inside(pRand);
704 G4cout<<"sn17.Inside(pRand) = "<<OutputInside(inside)<<G4endl ;
705 Dist = sn17.DistanceToOut(pRand,vRand,calcNorm,pgoodNorm,pNorm);
706 G4cout<<"sn17.DistanceToOut(pRand,vRand,...) = "<<Dist<<G4endl;
707 Dist = sn17.DistanceToIn(pRand,vRand);
708 G4cout<<"sn17.DistanceToIn(pRand,vRand) = "<<Dist<<G4endl;
709 // assert(ApproxEqual(Dist,8.9849541128705734394));
710
711
712
713
714// Checking G4Sphere::DistanceToIn(P)
715
716 Dist=s2.DistanceToIn(pzero);
717 assert(ApproxEqual(Dist,45));
718 Dist=s1.DistanceToIn(ponrmax1);
719 assert(ApproxEqual(Dist,0));
720
721// Checking G4Sphere::DistanceToIn(P,V)
722
723
724 zP = 3.;
725
726 for( i = 0; i < iMax; i++)
727 {
728 rRand = 0.5*rMax*G4UniformRand();
729 phi = twopi*G4UniformRand();
730 pRand = G4ThreeVector( rRand*std::cos(phi), rRand*std::sin(phi), zP);
731
732 cosTheta = std::cos( 180.*deg - 78.*deg*G4UniformRand() );
733 cosTheta = -0.99; // std::cos( 180.*deg - 78.*deg*G4UniformRand() );
734
735 sinTheta = std::sqrt((1.-cosTheta)*(1.+cosTheta));
736 phi = twopi*G4UniformRand();
737 vRand = G4ThreeVector( sinTheta*std::cos(phi), sinTheta*std::sin(phi), cosTheta);
738
739 Dist = sn17.DistanceToIn(pRand,vRand);
740
741 if( -Dist*cosTheta > 2.*zP )
742 {
743 G4cout<<"sn17.DistanceToIn(pRand,vRand), i = "<<i<<G4endl;
744 G4cout<<pRand.x()<<", "<<pRand.y()<<", "<<pRand.z()<<G4endl;
745 G4cout<<vRand.x()<<", "<<vRand.y()<<", "<<vRand.z()<<G4endl;
746 break;
747 }
748 }
749
750
751
752 pRand = G4ThreeVector( -10.276604989981144911, -4.7072500741022338389, 1);
753 vRand = G4ThreeVector( -0.28873162920537848164, 0.51647890054938394577, -0.80615357816219324061);
754 Dist = sn17.DistanceToIn(pRand,vRand);
755 G4cout<<"sn17.DistanceToIn(pRand,vRand) = "<<Dist<<G4endl;
756 // assert(ApproxEqual(Dist,8.9849541128705734394));
757
758
759
760
761 Dist=s1.DistanceToIn(ponzmax,vz);
762 G4cout<<"s1.DistanceToIn(ponzmax,vz) = "<<Dist<<G4endl;
763 Dist=s1.DistanceToIn(pbigy,vy);
764 assert(Dist==kInfinity);
765 Dist=s1.DistanceToIn(pbigy,vmy);
766 assert(ApproxEqual(Dist,50));
767
768 Dist=s2.DistanceToIn(pzero,vy);
769 G4cout<<"s2.DistanceToIn(pzero,vx) = "<<Dist<<G4endl;
770 assert(ApproxEqual(Dist,45));
771 Dist=s2.DistanceToIn(pzero,vmy);
772 assert(ApproxEqual(Dist,45));
773 Dist=s2.DistanceToIn(ponrmin1,vx);
774 // G4cout<<"s2.DistanceToIn(ponmin1,vx) = "<<Dist<<G4endl;
775 assert(Dist==0);
776 Dist=s2.DistanceToIn(ponrmin1,vmx);
777 assert(ApproxEqual(Dist,90));
778
779 Dist=s2.DistanceToIn(ponrmin2,vx);
780 assert(Dist==0);
781 Dist=s2.DistanceToIn(ponrmin2,vmx);
782 assert(ApproxEqual(Dist,90/std::sqrt(2.)));
783
784
785 Dist=s3.DistanceToIn(ptesttheta1,vmz);
786 // G4cout<<"s3.DistanceToIn(ptesttheta1,vmz) = "<<Dist<<G4endl;
787 assert(ApproxEqual(Dist,100-48/std::sqrt(2.)));
788 Dist=s3.DistanceToIn(pontheta1,vz);
789 // G4cout<<"s3.DistanceToIn(pontheta1,vz) = "<<Dist<<G4endl;
790 assert(Dist==kInfinity);
791 Dist=s3.DistanceToIn(pontheta1,vmz);
792 assert(Dist==0);
793 Dist=s3.DistanceToIn(pontheta2,vz);
794 // G4cout<<"s3.DistanceToIn(pontheta2,vz) = "<<Dist<<G4endl;
795 assert(Dist==0);
796 Dist=s3.DistanceToIn(pontheta2,vmz);
797 assert(Dist==kInfinity);
798 Dist=s32.DistanceToIn(pontheta1,vz);
799 // G4cout<<"s32.DistanceToIn(pontheta1,vz) = "<<Dist<<G4endl;
800 assert(Dist==0);
801 Dist=s32.DistanceToIn(pontheta1,vmz);
802 // G4cout<<"s32.DistanceToIn(pontheta1,vmz) = "<<Dist<<G4endl;
803 assert(Dist==kInfinity);
804 Dist=s33.DistanceToIn(pontheta2,vz);
805 // G4cout<<"s33.DistanceToIn(pontheta2,vz) = "<<Dist<<G4endl;
806 assert(Dist==kInfinity);
807 Dist=s33.DistanceToIn(pontheta2,vmz);
808 // G4cout<<"s33.DistanceToIn(pontheta2,vmz) = "<<Dist<<G4endl;
809 assert(Dist==0);
810
811 Dist=s4.DistanceToIn(pbigy,vmy);
812 assert(Dist==kInfinity);
813 Dist=s4.DistanceToIn(pbigz,vmz);
814 assert(ApproxEqual(Dist,50));
815 Dist=s4.DistanceToIn(pzero,vy);
816 assert(Dist==kInfinity);
817 Dist=s4.DistanceToIn(pzero,vx);
818 assert(ApproxEqual(Dist,45));
819
820 Dist=s4.DistanceToIn(ptestphi1,vx);
821 assert(ApproxEqual(Dist,100+45/std::sqrt(2.)));
822 Dist=s4.DistanceToIn(ponphi1,vmxmy);
823 assert(Dist==kInfinity);
824 Dist=s4.DistanceToIn(ponphi1,vxy);
825 // G4cout<<"s4.DistanceToIn(ponphi1,vxy) = "<<Dist<<G4endl;
826 assert(ApproxEqual(Dist,0));
827
828 Dist=s4.DistanceToIn(ptestphi2,vx);
829 assert(ApproxEqual(Dist,100+45/std::sqrt(2.)));
830 Dist=s4.DistanceToIn(ponphi2,vmxy);
831 assert(Dist==kInfinity);
832 Dist=s4.DistanceToIn(ponphi2,vxmy);
833 // G4cout<<"s4.DistanceToIn(ponphi2,vxmy) = "<<Dist<<G4endl;
834 assert(ApproxEqual(Dist,0));
835
836 Dist=s3.DistanceToIn(pzero,vx);
837 assert(ApproxEqual(Dist,45));
838 Dist=s3.DistanceToIn(ptesttheta1,vmz);
839 assert(ApproxEqual(Dist,100-48/std::sqrt(2.)));
840 Dist=b216.DistanceToIn(p216,v216);
841 // G4cout<<"b216.DistanceToIn(p216,v216) = "<<Dist<<G4endl;
842
843 Dist=b658.DistanceToIn(pzero,vz);
844 G4cout<<"b658.DistanceToIn(pzero,vz) = "<<Dist<<G4endl;
845
846
847
848 Dist=s13.DistanceToIn(G4ThreeVector(20.,0.,-70.),vz);
849 // G4cout<<"s13.DistanceToIn(G4ThreeVector(20.,0.,-70.),vz... = "<<Dist<<G4endl;
850 assert(ApproxEqual(Dist,58.452994616207498));
851
852 Dist=s13.DistanceToIn(G4ThreeVector(20.,0.,70.),vmz);
853 // G4cout<<"s13.DistanceToIn(G4ThreeVector(20.,0.,70.),vmz... = "<<Dist<<G4endl;
854 assert(ApproxEqual(Dist,35.358983848622451));
855
856
857 Dist=s14.DistanceToIn(G4ThreeVector(20.,0.,-70.),vz);
858 // G4cout<<"s14.DistanceToIn(G4ThreeVector(20.,0.,-70.),vz... = "<<Dist<<G4endl;
859 assert(ApproxEqual(Dist,35.358983848622451));
860
861 Dist=s14.DistanceToIn(G4ThreeVector(20.,0.,70.),vmz);
862 // G4cout<<"s14.DistanceToIn(G4ThreeVector(20.,0.,70.),vmz... = "<<Dist<<G4endl;
863 assert(ApproxEqual(Dist,58.452994616207498));
864
865 Dist=b1046.DistanceToIn(G4ThreeVector(0.,0.,4800*km),vmz);
866 G4cout<<"b1046.DistanceToIn(G4ThreeVector(0.,0.,4800*km),vmz... = "<<Dist<<G4endl;
867 // assert(ApproxEqual(Dist,0.));
868
869
870
871
872
873 ///////////////////////////////////////////////////////////////////////////
874
875
876// CalculateExtent
877 G4VoxelLimits limit; // Unlimited
878 G4RotationMatrix noRot;
879 G4AffineTransform origin;
880 G4double min,max;
881 assert(s1.CalculateExtent(kXAxis,limit,origin,min,max));
882 assert(min<=-50&&max>=50);
883 assert(s1.CalculateExtent(kYAxis,limit,origin,min,max));
884 assert(min<=-50&&max>=50);
885 assert(s1.CalculateExtent(kZAxis,limit,origin,min,max));
886 assert(min<=-50&&max>=50);
887
888 assert(s2.CalculateExtent(kXAxis,limit,origin,min,max));
889 assert(min<=-50&&max>=50);
890 assert(s2.CalculateExtent(kYAxis,limit,origin,min,max));
891 assert(min<=-50&&max>=50);
892 assert(s2.CalculateExtent(kZAxis,limit,origin,min,max));
893 assert(min<=-50&&max>=50);
894
895 assert(s3.CalculateExtent(kXAxis,limit,origin,min,max));
896 assert(min<=-50&&max>=50);
897 assert(s3.CalculateExtent(kYAxis,limit,origin,min,max));
898 assert(min<=-50&&max>=50);
899 assert(s3.CalculateExtent(kZAxis,limit,origin,min,max));
900 assert(min<=-50/std::sqrt(2.)&&max>=50/std::sqrt(2.));
901
902 G4ThreeVector pmxmymz(-100,-110,-120);
903 G4AffineTransform tPosOnly(pmxmymz);
904 assert(s1.CalculateExtent(kXAxis,limit,tPosOnly,min,max));
905 assert(min<=-150&&max>=-50);
906 assert(s1.CalculateExtent(kYAxis,limit,tPosOnly,min,max));
907 assert(min<=-160&&max>=-60);
908 assert(s1.CalculateExtent(kZAxis,limit,tPosOnly,min,max));
909 assert(min<=-170&&max>=-70);
910
911 assert(s3.CalculateExtent(kXAxis,limit,tPosOnly,min,max));
912 assert(min<=-150&&max>=-50);
913 assert(s3.CalculateExtent(kYAxis,limit,tPosOnly,min,max));
914 assert(min<=-160&&max>=-60);
915 assert(s3.CalculateExtent(kZAxis,limit,tPosOnly,min,max));
916 assert(min<=-170+50/std::sqrt(2.)&&max>=-70-50/std::sqrt(2.));
917
918 G4RotationMatrix r90Z;
919 r90Z.rotateZ(halfpi);
920 G4AffineTransform tRotZ(r90Z,pzero);
921 assert(s1.CalculateExtent(kXAxis,limit,tRotZ,min,max));
922 assert(min<=-50&&max>=50);
923 assert(s1.CalculateExtent(kYAxis,limit,tRotZ,min,max));
924 assert(min<=-50&&max>=50);
925 assert(s1.CalculateExtent(kZAxis,limit,tRotZ,min,max));
926 assert(min<=-50&&max>=50);
927
928 assert(s2.CalculateExtent(kXAxis,limit,tRotZ,min,max));
929 assert(min<=-50&&max>=50);
930 assert(s2.CalculateExtent(kYAxis,limit,tRotZ,min,max));
931 assert(min<=-50&&max>=50);
932 assert(s2.CalculateExtent(kZAxis,limit,tRotZ,min,max));
933 assert(min<=-50&&max>=50);
934
935 assert(s3.CalculateExtent(kXAxis,limit,tRotZ,min,max));
936 assert(min<=-50&&max>=50);
937 assert(s3.CalculateExtent(kYAxis,limit,tRotZ,min,max));
938 assert(min<=-50&&max>=50);
939 assert(s3.CalculateExtent(kZAxis,limit,tRotZ,min,max));
940 assert(min<=-50/std::sqrt(2.)&&max>=50/std::sqrt(2.));
941
942// Check that clipped away
943 G4VoxelLimits xClip;
944 xClip.AddLimit(kXAxis,-100,-60);
945 assert(!s1.CalculateExtent(kXAxis,xClip,origin,min,max));
946
947// Assert clipped to volume
948 G4VoxelLimits allClip;
949 allClip.AddLimit(kXAxis,-5,+5);
950 allClip.AddLimit(kYAxis,-5,+5);
951 allClip.AddLimit(kZAxis,-5,+5);
952 G4RotationMatrix genRot;
953 genRot.rotateX(pi/6);
954 genRot.rotateY(pi/6);
955 genRot.rotateZ(pi/6);
956 G4AffineTransform tGen(genRot,vx);
957 assert(s1.CalculateExtent(kXAxis,allClip,tGen,min,max));
958 assert(min<=-5&&max>=5);
959 assert(s1.CalculateExtent(kYAxis,allClip,tGen,min,max));
960 assert(min<=-5&&max>=5);
961 assert(s1.CalculateExtent(kZAxis,allClip,tGen,min,max));
962 assert(min<=-5&&max>=5);
963
964 assert(s2.CalculateExtent(kXAxis,allClip,tGen,min,max));
965 assert(min<=-5&&max>=5);
966 assert(s2.CalculateExtent(kYAxis,allClip,tGen,min,max));
967 assert(min<=-5&&max>=5);
968 assert(s2.CalculateExtent(kZAxis,allClip,tGen,min,max));
969 assert(min<=-5&&max>=5);
970
971 s34.CalculateExtent(kXAxis,allClip,tGen,min,max);
972 // G4cout<<"s34.CalculateExtent(kXAxis,allClip,tGen,min,max)"<<G4endl ;
973 // G4cout<<"min = "<<min<<" max = "<<max<<G4endl ;
974
975 s34.CalculateExtent(kYAxis,allClip,tGen,min,max);
976 // G4cout<<"s3.CalculateExtent(kYAxis,allClip,tGen,min,max)"<<G4endl ;
977 // G4cout<<"min = "<<min<<" max = "<<max<<G4endl ;
978
979 s34.CalculateExtent(kZAxis,allClip,tGen,min,max);
980 // G4cout<<"s34.CalculateExtent(kZAxis,allClip,tGen,min,max)"<<G4endl ;
981 // G4cout<<"min = "<<min<<" max = "<<max<<G4endl ;
982 assert(s34.CalculateExtent(kXAxis,allClip,tGen,min,max));
983 assert(min<=-5&&max>=5);
984 assert(s34.CalculateExtent(kYAxis,allClip,tGen,min,max));
985 assert(min<=-5&&max>=5);
986 assert(s34.CalculateExtent(kZAxis,allClip,tGen,min,max));
987 assert(min<=-5&&max>=5);
988
989// Test z clipping ok
990 for (G4double zTest=-100;zTest<100;zTest+=9)
991 {
992 G4VoxelLimits zTestClip;
993 zTestClip.AddLimit(kZAxis,-kInfinity,zTest);
994 if (zTest<-50)
995 {
996 assert(!s1.CalculateExtent(kZAxis,zTestClip,origin,min,max));
997 assert(!s2.CalculateExtent(kZAxis,zTestClip,origin,min,max));
998 }
999 else
1000 {
1001 assert(s1.CalculateExtent(kZAxis,zTestClip,origin,min,max));
1002 assert(s2.CalculateExtent(kZAxis,zTestClip,origin,min,max));
1003 G4double testMin=-50;
1004 G4double testMax=(zTest<50) ? zTest : 50;
1005 assert (ApproxEqual(min,testMin)
1006 &&ApproxEqual(max,testMax));
1007 }
1008 }
1009
1010// Test y clipping ok
1011 for (G4double xTest=-100;xTest<100;xTest+=9)
1012 {
1013 G4VoxelLimits xTestClip;
1014 xTestClip.AddLimit(kXAxis,-kInfinity,xTest);
1015 if (xTest<-50)
1016 {
1017 assert(!s1.CalculateExtent(kYAxis,xTestClip,origin,min,max));
1018 }
1019 else
1020 {
1021 assert(s1.CalculateExtent(kYAxis,xTestClip,origin,min,max));
1022// Calc max y coordinate
1023 G4double testMax=(xTest<0) ? std::sqrt(50*50-xTest*xTest) : 50;
1024 assert (ApproxEqual(min,-testMax)
1025 &&ApproxEqual(max,testMax));
1026 }
1027 }
1028
1029// Test x clipping ok
1030 for (G4double yTest=-100;yTest<100;yTest+=9)
1031 {
1032 G4VoxelLimits yTestClip;
1033 yTestClip.AddLimit(kYAxis,-kInfinity,yTest);
1034 if (yTest<-50)
1035 {
1036 assert(!s1.CalculateExtent(kXAxis,yTestClip,origin,min,max));
1037 }
1038 else
1039 {
1040 assert(s1.CalculateExtent(kXAxis,yTestClip,origin,min,max));
1041// Calc max y coordinate
1042 G4double testMax=(yTest<0) ? std::sqrt(50*50-yTest*yTest) : 50;
1043 assert (ApproxEqual(min,-testMax)
1044 &&ApproxEqual(max,testMax));
1045 }
1046 }
1047
1048 G4bool checkPoint( const G4Sphere& pSph, G4ThreeVector origin,
1049 G4double d, G4ThreeVector dir, EInside exp);
1050
1051 G4Sphere SpAroundX("SpAroundX", 10.*mm, 1000.*mm, -1.0*degree, 2.0*degree, 0.*degree, 180.0*degree );
1052
1053 G4double sinOneDeg = std::sin( 1.0 * degree );
1054 radOne = 100.0 * mm;
1055
1056 G4ThreeVector ptPhiSurfExct= G4ThreeVector( radOne * std::cos( -1.0 * degree ) ,
1057 - radOne * sinOneDeg,
1058 0.0 );
1059 G4cout << " Starting from point " << ptPhiSurfExct << G4endl;
1060 G4cout << " using direction " << vy << G4endl;
1061 checkPoint( SpAroundX, ptPhiSurfExct, -radOne * kAngTolerance * 1.5,
1062 vy, kOutside);
1063 checkPoint( SpAroundX, ptPhiSurfExct, -radOne * kAngTolerance * 0.49,
1064 vy, kSurface);
1065 checkPoint( SpAroundX, ptPhiSurfExct, -radOne * kAngTolerance * 0.25,
1066 vy, kSurface);
1067 checkPoint( SpAroundX, ptPhiSurfExct, 0.0,
1068 vy, kSurface);
1069 checkPoint( SpAroundX, ptPhiSurfExct, radOne * kAngTolerance * 0.25,
1070 vy, kSurface);
1071 checkPoint( SpAroundX, ptPhiSurfExct, radOne * kAngTolerance * 0.49,
1072 vy, kSurface);
1073 checkPoint( SpAroundX, ptPhiSurfExct, radOne * kAngTolerance * 1.5,
1074 vy, kInside);
1075
1076 // Try one that has a 'deep' negative phi section
1077 // --> Vlad. Grichine test case, 30 Oct 2003
1078 //
1079 G4cout << G4endl << G4endl << "" << G4endl;
1080 G4cout << "========================================================= " << G4endl;
1081
1082 G4Sphere SphDeepNeg("DeepNegPhiSphere", 10.*mm, 1000.*mm,
1083 -270.0*degree, 280.0*degree, // start Phi, delta Phi
1084 0.*degree, 180.0*degree ); // start Theta, delta Theta
1085 G4double phiPoint = 160.0 * degree;
1086 G4ThreeVector StartPt( radOne * std::cos(phiPoint), radOne * std::sin(phiPoint), 0.0);
1087 G4cout << "For sphere " << SphDeepNeg.GetName() << G4endl;
1088 G4cout << " Starting from point " << ptPhiSurfExct << G4endl;
1089
1090 checkPoint( SphDeepNeg, StartPt, 0.0, vy, kInside);
1091
1092 // Try the edges
1093 G4ThreeVector NegEdgePt( radOne * std::cos(-270.0*degree), radOne * std::sin(-270.0*degree), 0.0);
1094 G4ThreeVector PosEdgePt( radOne * std::cos(10.0*degree), radOne * std::sin(10.0*degree), 0.0);
1095
1096 G4cout << "--------------------------------------------------------" << G4endl;
1097 G4cout << " New point " << NegEdgePt << " should be at Neg edge of -270.0 degrees " << G4endl;
1098 checkPoint( SphDeepNeg, NegEdgePt, 0.0, -vx, kSurface);
1099 checkPoint( SphDeepNeg, NegEdgePt, radOne*kAngTolerance * 0.25, -vx, kSurface);
1100 checkPoint( SphDeepNeg, NegEdgePt, -radOne*kAngTolerance * 0.25, -vx, kSurface);
1101 checkPoint( SphDeepNeg, NegEdgePt, radOne*kAngTolerance * 1.25, -vx, kInside);
1102 checkPoint( SphDeepNeg, NegEdgePt, -radOne*kAngTolerance * 1.25, -vx, kOutside);
1103
1104 G4cout << "--------------------------------------------------------" << G4endl;
1105 G4cout << " New point " << PosEdgePt << " should be at Pos edge of +10.0 degrees " << G4endl;
1106 checkPoint( SphDeepNeg, PosEdgePt, 0.0, -vy, kSurface);
1107 checkPoint( SphDeepNeg, PosEdgePt, radOne*kAngTolerance * 0.25, -vy, kSurface);
1108 checkPoint( SphDeepNeg, PosEdgePt, -radOne*kAngTolerance * 0.25, -vy, kSurface);
1109 checkPoint( SphDeepNeg, PosEdgePt, -radOne*kAngTolerance * 1.25, -vy, kOutside);
1110 checkPoint( SphDeepNeg, PosEdgePt, radOne*kAngTolerance * 1.25, -vy, kInside);
1111
1112
1113
1114 // G4double sTheta = 135*degree, pxt = -10., pyt = 10.;
1115
1116 // G4cout<<"sTheta = "<<sTheta/degree<<" degree"<<G4endl;
1117 // G4cout<<"std::tan(sTheta) = "<<std::tan(sTheta)<<G4endl;
1118 // G4cout<<"std::sin(sTheta) = "<<std::sin(sTheta)<<G4endl;
1119 // G4cout<<"std::atan(sTheta) = "<<std::atan(sTheta)/degree<<G4endl;
1120 // G4cout<<"std::atan2(pyt,pxt) = pyt/pxt = "<<std::atan2(pyt,pxt)/degree<<G4endl;
1121
1122 return 0;
1123}
1124
1125// Given the sphere 'rSphere', a point 'origin'
1126// --> --> -->
1127// the function below checks that the point pnew=( origin + dist * direction )
1128// - the point (p+dist*dir) is located in the part of the solid given by 'expectedInResult'
1129// - and that from there, the DistanceToIn along 'dir' is not negative or Infinite
1130//
1131// Use cases expected:
1132// - 'origin' is on/near a surface and 'direction' is pointing towards the inside of the solid
1133
1134G4bool
1135checkPoint( const G4Sphere &rSphere,
1136 G4ThreeVector origin,
1137 G4double dist,
1138 G4ThreeVector direction,
1139 EInside expectedInResult)
1140{
1141 G4int verbose = 0;
1142
1143 G4ThreeVector newPoint;
1144 G4double distIn=-1.0, distOut=-1.0;
1145
1146 newPoint = origin + dist * direction;
1147
1148 G4int oldPrecision= G4cout.precision(10);
1149 // G4cout << " --- Sphere " << rSphere.GetName() << "" << G4endl;
1150 if( verbose > 0 ) {
1151 G4cout << G4endl;
1152 if (verbose > 2 ) G4cout << " Sphere " << rSphere.GetName();
1153 G4cout.precision(10);
1154 if (verbose > 1 ) G4cout << " dir= " << direction;
1155 G4cout << " dist= " << dist;
1156 }
1157
1158 EInside inSphere= rSphere.Inside( newPoint ) ;
1159 /*======*/
1160 G4cout.precision(15);
1161 // G4cout << " NewPoint " << newPoint << " is "
1162 G4bool goodIn= (inSphere == expectedInResult) ;
1163 if ( !goodIn ) {
1164 G4cout << " ************ Unexpected Result for Inside *************** " << G4endl;
1165 }
1166 if ( verbose || !goodIn ) {
1167 G4cout << " New point "
1168 << " is " << OutputInside( inSphere )
1169 << " vs " << OutputInside( expectedInResult )
1170 << " expected." << G4endl ;
1171 }
1172
1173 G4bool goodDistIn = true;
1174
1175 distIn = rSphere.DistanceToIn( newPoint, direction );
1176 /*===========*/
1177 if ( verbose ) G4cout << " DistToIn (p, dir) = " << distIn << G4endl;
1178 if( (inSphere == kOutside)
1179 && (distIn < 0.0 ) // Cannot use 0.5*kCarTolerance for Angular tolerance!!
1180 ){
1181 G4cout << " ********** Unexpected Result for DistanceToIn from outside ********* " << G4endl;
1182 // G4cout << " It should be " << G4endl;
1183 goodDistIn = false;
1184 }
1185 if( (inSphere == kSurface )
1186 && ( (distIn < 0.0) || (distIn >= kInfinity ))
1187 ){
1188 G4cout << " ********** Unexpected Result for DistanceToIn on surface ********* " << G4endl;
1189 // if ( (distIn != 0.0) )
1190 // - Can check that the return value must be 0.0
1191 // But in general case the direction can be away from the solid,
1192 // and then a finite or kInfinity answer is correct
1193 // --> must check the direction against the normal
1194 // in order to perform this check in general case.
1195
1196 goodDistIn = false;
1197 }
1198 if ( verbose || !goodDistIn ) {
1199 G4cout << " DistToIn (p, dir) = " << distIn << G4endl;
1200 }
1201
1202 G4bool good= (goodIn && goodDistIn);
1203 if ( !good ){
1204 // There was an error -- document the use case!
1205 G4cout << " --- Sphere " << rSphere.GetName() << "" << G4endl;
1206 G4cout << " Origin= " << origin << G4endl;
1207 G4cout << " Direction= " << direction << G4endl;
1208 G4cout << " dist= " << dist;
1209 G4cout << " Actual-point= " << newPoint << G4endl;
1210 }
1211
1212 distOut = rSphere.DistanceToOut( newPoint, direction );
1213 /*=============*/
1214 if ( verbose ) G4cout << " DistToOut (p, dir) = " << distOut << G4endl;
1215
1216 G4cout.precision(oldPrecision);
1217
1218 return good;
1219}
1220
Note: See TracBrowser for help on using the repository browser.