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

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

geant4 tag 9.4
File size: 47.4 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.30 2010/03/24 09:50:03 grichine 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
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 EInside inside;
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
295 G4Sphere s18_03_10("s18_03_10", 47*mm, 48*mm, 0*deg, 200*deg, 80*deg, 100*deg);
296 G4ThreeVector p18_03_10(43.37539710867407*mm, 16.12036885157033*mm, -8.224548565698871*mm);
297 G4ThreeVector v18_03_10(-0.4161958548132512, 0.6603942936714858, -0.6250283092167956);
298
299 inside = s18_03_10.Inside(p18_03_10) ;
300 G4cout<<"s18_03_10.Inside(p18_03_10 ... = "<<OutputInside(inside)<<G4endl ;
301
302
303
304#ifdef NDEBUG
305 G4Exception("FAIL: *** Assertions must be compiled in! ***");
306#endif
307
308 G4cout.precision(20);
309
310 //////////////// Check name /////////////////////////
311
312 assert(s1.GetName()=="Solid G4Sphere");
313
314 // check cubic volume
315
316 vol = s1.GetCubicVolume();
317 volCheck = 4*pi*125000/3;
318
319 assert( ApproxEqual(vol,volCheck));
320
321 assert(ApproxEqual(s1.GetCubicVolume(),4*pi*125000/3));
322
323
324 // Some user application cases
325
326 Dist = s4.DistanceToOut(ponrmin3,vmz,calcNorm,pgoodNorm,pNorm) ;
327 // G4cout<<"s4.DistanceToOut(ponrmin3,vmz,calcNorm,pgoodNorm,pNorm) = "<<Dist
328 // <<G4endl ;
329 Dist = s6.DistanceToOut(ponrminJ,vmz,calcNorm,pgoodNorm,pNorm) ;
330 // G4cout<<"s6.DistanceToOut(ponrminJ,vmz,calcNorm,pgoodNorm,pNorm) = "<<Dist
331 // <<G4endl ;
332 Dist = s6.DistanceToOut(ponrmaxJ,vmz,calcNorm,pgoodNorm,pNorm) ;
333 // G4cout<<"s6.DistanceToOut(ponrmaxJ,vmz,calcNorm,pgoodNorm,pNorm) = "<<Dist
334 // <<G4endl ;
335
336 Dist = s7.DistanceToOut(G4ThreeVector(1399.984667238032,
337 5.9396696802500299,
338 -2.7661927818688308),
339 G4ThreeVector(0.50965504781062942,
340 -0.80958849145715217,
341 0.29123565499656401),
342 calcNorm,pgoodNorm,pNorm) ;
343 // G4cout<<"s7.DistanceToOut(shereP,sphereV,calcNorm,pgoodNorm,pNorm) = "<<Dist
344 // <<G4endl ;
345
346 Dist = s7.DistanceToIn(G4ThreeVector(1399.984667238032,
347 5.9396696802500299,
348 -2.7661927818688308),
349 G4ThreeVector(0.50965504781062942,
350 -0.80958849145715217,
351 0.29123565499656401)) ;
352 // G4cout<<"s7.DistanceToIn(shereP,sphereV,calcNorm,pgoodNorm,pNorm) = "<<Dist
353 // <<G4endl ;
354
355
356
357// Check G4Sphere::Inside
358
359 inside = s7.Inside(G4ThreeVector(1399.984667238032,
360 5.9396696802500299,
361 -2.7661927818688308) ) ;
362 // G4cout<<"s7.Inside(G4ThreeVector(1399.98466 ... = "
363 // <<OutputInside(inside)<<G4endl ;
364
365 inside = s8.Inside(G4ThreeVector(-249.5020724528353*mm,
366 26.81253142743162*mm,
367 114.8988524453591*mm ) ) ;
368 // G4cout<<"s8.Inside(G4ThreeVector(-249.5020 ... = "<<OutputInside(inside)<<G4endl ;
369 inside = b216.Inside(p216);
370 // G4cout<<"b216.Inside(p216) = "<<OutputInside(inside)<<G4endl ;
371 inside = s1.Inside(pz);
372 // G4cout<<"s1.Inside(pz) = "<<OutputInside(inside)<<G4endl ;
373
374 inside = s9.Inside(s9p);
375 // G4cout<<"s9.Inside(s9p) = "<<OutputInside(inside)<<G4endl ;
376
377 inside = b402.Inside(p402);
378 // G4cout<<"p402.Inside(p402) = "<<OutputInside(inside)<<G4endl ;
379
380 inside = s10.Inside(s10p);
381 // G4cout<<"s10.Inside(s10p) = "<<OutputInside(inside)<<G4endl ;
382 // G4cout<<"p radius = "<<s10p.mag()<<G4endl ;
383
384 inside = s11.Inside(ps11);
385 // G4cout<<"s11.Inside(ps11) = "<<OutputInside(inside)<<G4endl ;
386 // G4cout<<"ps11.mag() = "<<ps11.mag()<<G4endl ;
387
388 inside = sLHCB.Inside(pLHCB);
389 // G4cout<<"sLHCB.Inside(pLHCB) = "<<OutputInside(inside)<<G4endl ;
390 // G4cout<<"pLHCB.mag() = "<<pLHCB.mag()<<G4endl ;
391
392 inside = spAroundX.Inside(ptPhiMinus);
393 // G4cout<<"spAroundX.Inside(ptPhiMinus) = "<<OutputInside(inside)<<G4endl ;
394 inside = b658.Inside(p658);
395 G4cout<<"b658.Inside(p658) = "<<OutputInside(inside)<<G4endl ;
396
397 assert(s1.Inside(pzero)==kInside);
398 // assert(s1.Inside(pz)==kInside);
399 assert(s2.Inside(pzero)==kOutside);
400 assert(s2.Inside(ponrmin2)==kSurface);
401 assert(s2.Inside(ponrmax2)==kSurface);
402 assert(s3.Inside(pontheta1)==kSurface);
403 assert(s3.Inside(pontheta2)==kSurface);
404 assert(s4.Inside(ponphi1)==kSurface);
405 assert(s4.Inside(ponphi1)==kSurface);
406 assert(s4.Inside(pOverPhi)==kOutside);
407 assert(s4.Inside(pInPhi)==kInside);
408 assert(s5.Inside(pbigz)==kOutside);
409
410 assert(s41.Inside(pmx)==kSurface);
411 assert(s42.Inside(pmx)==kSurface);
412
413// Checking G4Sphere::SurfaceNormal
414 G4double p2=1./std::sqrt(2.),p3=1./std::sqrt(3.);
415
416
417 norm=sn1.SurfaceNormal(G4ThreeVector(0.,0.,50.));
418 assert(ApproxEqual(norm,G4ThreeVector(p3,p3,p3)));
419 norm=sn1.SurfaceNormal(G4ThreeVector(0.,0.,0.));
420 assert(ApproxEqual(norm,G4ThreeVector(p2,p2,0.)));
421
422 norm=sn11.SurfaceNormal(G4ThreeVector(0.,0.,0.));
423 assert(ApproxEqual(norm,G4ThreeVector(0.,0.,-1.)));
424 norm=sn12.SurfaceNormal(G4ThreeVector(0.,0.,0.));
425 assert(ApproxEqual(norm,G4ThreeVector(0.,0.,-1.)));
426 norm=sn13.SurfaceNormal(G4ThreeVector(0.,0.,0.));
427 assert(ApproxEqual(norm,G4ThreeVector(0.,0.,1.)));
428
429 norm=sn2.SurfaceNormal(G4ThreeVector(-45.,0.,0.));
430 assert(ApproxEqual(norm,G4ThreeVector(p2,-p2,0.)));
431
432
433
434 norm=s1.SurfaceNormal(ponrmax1);
435 assert(ApproxEqual(norm,vx));
436
437// Checking G4Sphere::DistanceToOut(P)
438 Dist=s1.DistanceToOut(pzero);
439 assert(ApproxEqual(Dist,50));
440 Dist=s1.DistanceToOut(ponrmax1);
441 assert(ApproxEqual(Dist,0));
442
443// Checking G4Sphere::DistanceToOut(p,v)
444
445
446 Dist=s1.DistanceToOut(pz,vz,calcNorm,pgoodNorm,pNorm);
447 // G4cout<<"Dist=s1.DistanceToOut(pz,vz) = "<<Dist<<G4endl;
448 assert(ApproxEqual(Dist,20.));
449
450 Dist=s1.DistanceToOut(ponrmax1,vx,calcNorm,pgoodNorm,pNorm);
451 *pNorm=pNorm->unit();
452 assert(ApproxEqual(Dist,0)&&*pgoodNorm&&ApproxEqual(*pNorm,vx));
453
454 Dist=s2.DistanceToOut(ponrmin1,vx,calcNorm,pgoodNorm,pNorm);
455 *pNorm=pNorm->unit();
456 assert(ApproxEqual(Dist,5)&&*pgoodNorm&&ApproxEqual(*pNorm,vx));
457
458 Dist=s2.DistanceToOut(ponrmax2,vx,calcNorm,pgoodNorm,pNorm);
459 assert(ApproxEqual(Dist,0)&&*pgoodNorm&&ApproxEqual(*pNorm,vxy));
460
461 Dist=s1.DistanceToOut(pzero,vx,calcNorm,pgoodNorm,pNorm);
462 assert(ApproxEqual(Dist,50)&&ApproxEqual(pNorm->unit(),vx)&&*pgoodNorm);
463 Dist=s1.DistanceToOut(pzero,vmx,calcNorm,pgoodNorm,pNorm);
464 assert(ApproxEqual(Dist,50)&&ApproxEqual(pNorm->unit(),vmx)&&*pgoodNorm);
465 Dist=s1.DistanceToOut(pzero,vy,calcNorm,pgoodNorm,pNorm);
466 assert(ApproxEqual(Dist,50)&&ApproxEqual(pNorm->unit(),vy)&&*pgoodNorm);
467 Dist=s1.DistanceToOut(pzero,vmy,calcNorm,pgoodNorm,pNorm);
468 assert(ApproxEqual(Dist,50)&&ApproxEqual(pNorm->unit(),vmy)&&*pgoodNorm);
469 Dist=s1.DistanceToOut(pzero,vz,calcNorm,pgoodNorm,pNorm);
470 assert(ApproxEqual(Dist,50)&&ApproxEqual(pNorm->unit(),vz)&&*pgoodNorm);
471 Dist=s1.DistanceToOut(pzero,vmz,calcNorm,pgoodNorm,pNorm);
472 assert(ApproxEqual(Dist,50)&&ApproxEqual(pNorm->unit(),vmz)&&*pgoodNorm);
473 Dist=s1.DistanceToOut(pzero,vxy,calcNorm,pgoodNorm,pNorm);
474 assert(ApproxEqual(Dist,50)&&ApproxEqual(pNorm->unit(),vxy)&&*pgoodNorm);
475
476 Dist=s4.DistanceToOut(ponphi1,vx,calcNorm,pgoodNorm,pNorm);
477 //assert(ApproxEqual(Dist,0)&&ApproxEqual(pNorm->unit(),vmxmy)&&*pgoodNorm);
478 Dist=s4.DistanceToOut(ponphi2,vx,calcNorm,pgoodNorm,pNorm);
479 // assert(ApproxEqual(Dist,0)&&ApproxEqual(pNorm->unit(),vmxy)&&*pgoodNorm);
480 Dist=s3.DistanceToOut(pontheta1,vz,calcNorm,pgoodNorm,pNorm);
481 // assert(ApproxEqual(Dist,0)&&ApproxEqual(pNorm->unit(),vy)&&*pgoodNorm);
482 Dist=s32.DistanceToOut(pontheta1,vmz,calcNorm,pgoodNorm,pNorm);
483 //assert(ApproxEqual(Dist,0)&&ApproxEqual(pNorm->unit(),vmy)&&*pgoodNorm);
484 Dist=s32.DistanceToOut(pontheta1,vz,calcNorm,pgoodNorm,pNorm);
485 //assert(ApproxEqual(Dist,50)&&ApproxEqual(pNorm->unit(),vz)&&*pgoodNorm);
486 Dist=s1.DistanceToOut(pzero,vmz,calcNorm,pgoodNorm,pNorm);
487 assert(ApproxEqual(Dist,50)&&ApproxEqual(pNorm->unit(),vmz)&&*pgoodNorm);
488 Dist=s1.DistanceToOut(pzero,vxy,calcNorm,pgoodNorm,pNorm);
489 assert(ApproxEqual(Dist,50)&&ApproxEqual(pNorm->unit(),vxy)&&*pgoodNorm);
490
491
492 Dist=s2.DistanceToOut(ponrmin1,vxy,calcNorm,pgoodNorm,pNorm);
493 // G4cout<<"Dist=s2.DistanceToOut(pormin1,vxy) = "<<Dist<<G4endl;
494
495 Dist=s2.DistanceToOut(ponrmax1,vmx,calcNorm,pgoodNorm,pNorm);
496 // G4cout<<"Dist=s2.DistanceToOut(ponxside,vmx) = "<<Dist<<G4endl;
497 Dist=s2.DistanceToOut(ponrmax1,vmxmy,calcNorm,pgoodNorm,pNorm);
498 // G4cout<<"Dist=s2.DistanceToOut(ponxside,vmxmy) = "<<Dist<<G4endl;
499 Dist=s2.DistanceToOut(ponrmax1,vz,calcNorm,pgoodNorm,pNorm);
500 // G4cout<<"Dist=s2.DistanceToOut(ponxside,vz) = "<<Dist<<G4endl;
501
502 // Dist=s2.DistanceToOut(pbigx,vx,calcNorm,pgoodNorm,pNorm);
503 // G4cout<<"Dist=s2.DistanceToOut(pbigx,vx) = "<<Dist<<G4endl;
504 // Dist=s2.DistanceToOut(pbigx,vxy,calcNorm,pgoodNorm,pNorm);
505 // G4cout<<"Dist=s2.DistanceToOut(pbigx,vxy) = "<<Dist<<G4endl;
506 // Dist=s2.DistanceToOut(pbigx,vz,calcNorm,pgoodNorm,pNorm);
507 // G4cout<<"Dist=s2.DistanceToOut(pbigx,vz) = "<<Dist<<G4endl;
508 //Test Distance for phi section
509 Dist=sn22.DistanceToOut(G4ThreeVector(0.,49.,0.),vmy,calcNorm,pgoodNorm,pNorm);
510 assert(ApproxEqual(Dist,49.));
511 Dist=sn22.DistanceToOut(G4ThreeVector(-45.,0.,0.),vx,calcNorm,pgoodNorm,pNorm);
512 assert(ApproxEqual(Dist,45.));
513G4cout<<"Dist from Center ="<<sn22.DistanceToOut(G4ThreeVector(0.,49.,0),G4ThreeVector(0,-1,0))<<G4endl;
514G4cout<<"Dist from Center ="<<sn22.DistanceToOut(G4ThreeVector(-45.,0.,0),G4ThreeVector(1,0,0))<<G4endl;
515
516 //
517 Dist=b216.DistanceToOut(p216,v216,calcNorm,pgoodNorm,pNorm);
518
519 // G4cout<<"b216.DistanceToOut(p216,v216,... = "<<Dist<<G4endl;
520
521 // call from outside
522 // Dist=sAlex.DistanceToOut(psAlex,vsAlex,calcNorm,pgoodNorm,pNorm);
523 // G4cout<<"sAlex.DistanceToOut(psAlex,vsAlex,... = "<<Dist<<G4endl;
524
525
526 // Dist=s9.DistanceToIn(s9p,s9v);
527 // G4cout<<"s9.DistanceToIn(s9p,s9v,... = "<<Dist<<G4endl;
528 Dist=s13.DistanceToOut(G4ThreeVector(20.,0.,0.),vz,calcNorm,pgoodNorm,pNorm);
529 // G4cout<<"s13.DistanceToOut(G4ThreeVector(20.,0.,0.),vz... = "<<Dist<<G4endl;
530 assert(ApproxEqual(Dist,34.641016151377549));
531
532 Dist=s13.DistanceToOut(G4ThreeVector(20.,0.,0.),vmz,calcNorm,pgoodNorm,pNorm);
533 // G4cout<<"s13.DistanceToOut(G4ThreeVector(20.,0.,0.),vmz... = "<<Dist<<G4endl;
534 assert(ApproxEqual(Dist,11.547005383792508));
535
536 Dist=s14.DistanceToOut(G4ThreeVector(20.,0.,0.),vz,calcNorm,pgoodNorm,pNorm);
537 // G4cout<<"s14.DistanceToOut(G4ThreeVector(20.,0.,0.),vz... = "<<Dist<<G4endl;
538 assert(ApproxEqual(Dist,11.547005383792508));
539
540 Dist=s14.DistanceToOut(G4ThreeVector(20.,0.,0.),vmz,calcNorm,pgoodNorm,pNorm);
541 // G4cout<<"s14.DistanceToOut(G4ThreeVector(20.,0.,0.),vmz... = "<<Dist<<G4endl;
542 assert(ApproxEqual(Dist,34.641016151377549));
543
544 Dist=sb830.DistanceToOut(pb830,vb830,calcNorm,pgoodNorm,pNorm);
545 G4cout<<"sb830.DistanceToOut(pb830,vb830... = "<<Dist<<G4endl;
546 inside = sb830.Inside(pb830+Dist*vb830);
547 G4cout<<"sb830.Inside(pb830+Dist*vb830) = "<<OutputInside(inside)<<G4endl ;
548
549 Dist=sn11.DistanceToOut( G4ThreeVector(0.,0.,20.), vmz, calcNorm, pgoodNorm, pNorm);
550 // G4cout<<"sn11.DistanceToOut(G4ThreeVector(0.,0.,20.),vmz... = "<<Dist<<G4endl;
551 assert(ApproxEqual(Dist,20.));
552
553 Dist=sn11.DistanceToOut(G4ThreeVector(0.,0.,0.),vmz,calcNorm,pgoodNorm,pNorm);
554 // G4cout<<"sn11.DistanceToOut(G4ThreeVector(0.,0.,0.),vmz... = "<<Dist<<G4endl;
555 assert(ApproxEqual(Dist,0.));
556
557 Dist=sn11.DistanceToOut(G4ThreeVector(0.,0.,0.),vz,calcNorm,pgoodNorm,pNorm);
558 // G4cout<<"sn11.DistanceToOut(G4ThreeVector(0.,0.,0.),vmz... = "<<Dist<<G4endl;
559 assert(ApproxEqual(Dist,50.));
560
561 Dist=sn11.DistanceToOut(G4ThreeVector(10.,0.,0.),vmz,calcNorm,pgoodNorm,pNorm);
562 // G4cout<<"sn11.DistanceToOut(G4ThreeVector(10.,0.,0.),vmz... = "<<Dist<<G4endl;
563 assert(ApproxEqual(Dist,0.));
564
565 Dist=sn12.DistanceToOut(G4ThreeVector(0.,0.,20.),vxmz,calcNorm,pgoodNorm,pNorm);
566 // G4cout<<"sn12.DistanceToOut(G4ThreeVector(0.,0.,20.),vxmz... = "<<Dist<<G4endl;
567 assert(ApproxEqual(Dist,20./std::sqrt(2.)));
568
569 Dist=sn12.DistanceToOut(G4ThreeVector(0.,0.,0.),vmz,calcNorm,pgoodNorm,pNorm);
570 // G4cout<<"sn12.DistanceToOut(G4ThreeVector(0.,0.,0.),vmz... = "<<Dist<<G4endl;
571 assert(ApproxEqual(Dist,0.));
572
573 Dist=sn12.DistanceToOut(G4ThreeVector(10.,0.,10.),vz,calcNorm,pgoodNorm,pNorm);
574 // G4cout<<"sn12.DistanceToOut(G4ThreeVector(10.,0.,10.),vz... = "<<Dist<<G4endl;
575 assert(ApproxEqual(Dist,38.989794855663561179));
576
577
578 Dist=sn12.DistanceToOut(G4ThreeVector(10.,0.,10.),vmz,calcNorm,pgoodNorm,pNorm);
579 // G4cout<<"sn12.DistanceToOut(G4ThreeVector(10.,0.,10.),vmz... = "<<Dist<<G4endl;
580 assert(ApproxEqual(Dist,0.));
581
582
583 Dist=sn14.DistanceToOut(G4ThreeVector(10.,0.,10.),vz,calcNorm,pgoodNorm,pNorm);
584 // G4cout<<"sn14.DistanceToOut(G4ThreeVector(10.,0.,10.),vz... = "<<Dist<<G4endl;
585 assert(ApproxEqual(Dist,0.));
586
587
588 Dist=sn14.DistanceToOut(G4ThreeVector(10.,0.,5.),vz,calcNorm,pgoodNorm,pNorm);
589 // G4cout<<"sn14.DistanceToOut(G4ThreeVector(10.,0.,5.),vz... = "<<Dist<<G4endl;
590 assert(ApproxEqual(Dist,5.));
591
592
593 Dist=sn14.DistanceToOut(G4ThreeVector(10.,0.,5.),vmxz,calcNorm,pgoodNorm,pNorm);
594 // G4cout<<"sn14.DistanceToOut(G4ThreeVector(10.,0.,5.),vmxz... = "<<Dist<<G4endl;
595 assert(ApproxEqual(Dist,3.5355339059327381968));
596
597 Dist=sn14.DistanceToOut(G4ThreeVector(10.,0.,10.),vmxz,calcNorm,pgoodNorm,pNorm);
598 // G4cout<<"sn14.DistanceToOut(G4ThreeVector(10.,0.,10.),vmxz... = "<<Dist<<G4endl;
599 assert(ApproxEqual(Dist,0.));
600
601
602
603 iMax = 10000;
604 rMax = 49.;
605 zP = -1.;
606
607 for( i = 0; i < iMax; i++)
608 {
609 rRand = rMax*G4UniformRand();
610 phi = twopi*G4UniformRand();
611 pRand = G4ThreeVector( rRand*std::cos(phi), rRand*std::sin(phi), zP);
612
613 cosTheta = 0.99*G4UniformRand();
614
615 sinTheta = std::sqrt((1.-cosTheta)*(1.+cosTheta));
616 phi = twopi*G4UniformRand();
617 vRand = G4ThreeVector( sinTheta*std::cos(phi), sinTheta*std::sin(phi), cosTheta);
618
619 Dist = sn15.DistanceToOut(pRand,vRand,calcNorm,pgoodNorm,pNorm);
620
621 if( Dist*cosTheta > -2.*zP )
622 {
623 G4cout<<"sn15.DistanceToOut(pRand,vRand,..., i = "<<i<<G4endl;
624 G4cout<<pRand.x()<<", "<<pRand.y()<<", "<<pRand.z()<<G4endl;
625 G4cout<<vRand.x()<<", "<<vRand.y()<<", "<<vRand.z()<<G4endl;
626 break;
627 }
628 }
629 pRand = G4ThreeVector( -9.3576454272980740257, 10.436745988128407703, -1.);
630 vRand = G4ThreeVector( -0.064906521070262901407, 0.35249758429187183495, 0.93355910181999202102);
631 Dist = sn15.DistanceToOut(pRand,vRand,calcNorm,pgoodNorm,pNorm);
632 // G4cout<<"sn15.DistanceToOut(pRand,vRand... = "<<Dist<<G4endl;
633 assert(ApproxEqual(Dist,1.3409673565670394701));
634
635
636 zP = 1.;
637
638 for( i = 0; i < iMax; i++)
639 {
640 rRand = rMax*G4UniformRand();
641 phi = twopi*G4UniformRand();
642 pRand = G4ThreeVector( rRand*std::cos(phi), rRand*std::sin(phi), zP);
643
644 cosTheta = -0.99*G4UniformRand();
645
646 sinTheta = std::sqrt((1.-cosTheta)*(1.+cosTheta));
647 phi = twopi*G4UniformRand();
648 vRand = G4ThreeVector( sinTheta*std::cos(phi), sinTheta*std::sin(phi), cosTheta);
649
650 Dist = sn16.DistanceToOut(pRand,vRand,calcNorm,pgoodNorm,pNorm);
651
652 if( -Dist*cosTheta > 2.*zP )
653 {
654 G4cout<<"sn16.DistanceToOut(pRand,vRand,..., i = "<<i<<G4endl;
655 G4cout<<pRand.x()<<", "<<pRand.y()<<", "<<pRand.z()<<G4endl;
656 G4cout<<vRand.x()<<", "<<vRand.y()<<", "<<vRand.z()<<G4endl;
657 break;
658 }
659 }
660
661
662
663 pRand = G4ThreeVector( -18.872396210750576273, 11.573660000258405134, 1);
664 vRand = G4ThreeVector( -0.85674495247509219187, -0.51247617288646407641, -0.057933225631713866632);
665 Dist = sn16.DistanceToOut(pRand,vRand,calcNorm,pgoodNorm,pNorm);
666 // G4cout<<"sn16.DistanceToOut(pRand,vRand... = "<<Dist<<G4endl;
667 assert(ApproxEqual(Dist,8.9849541128705734394));
668
669
670
671 zP = -2.;
672
673 for( i = 0; i < iMax; i++)
674 {
675 rRand = rMax*G4UniformRand();
676 phi = twopi*G4UniformRand();
677
678 pRand = G4ThreeVector( rRand*std::cos(phi)*std::sin(91.*deg),
679 rRand*std::sin(phi)*std::sin(91.*deg),
680 rRand*std::cos(91.*deg));
681
682 cosTheta = std::cos( 180.*deg-89.*deg*G4UniformRand() );
683
684 sinTheta = std::sqrt((1.-cosTheta)*(1.+cosTheta));
685 phi = twopi*G4UniformRand();
686 vRand = G4ThreeVector( sinTheta*std::cos(phi), sinTheta*std::sin(phi), cosTheta);
687
688 Dist = sn17.DistanceToOut(pRand,vRand,calcNorm,pgoodNorm,pNorm);
689
690 if( -Dist*cosTheta > 50. )
691 {
692 G4cout<<"sn17.DistanceToOut(pRand,vRand,..., i = "<<i<<G4endl;
693 G4cout<<pRand.x()<<", "<<pRand.y()<<", "<<pRand.z()<<G4endl;
694 G4cout<<vRand.x()<<", "<<vRand.y()<<", "<<vRand.z()<<G4endl;
695 break;
696 }
697 }
698
699 pRand = G4ThreeVector(16.20075504802145, -22.42917454903122, -41.6469406430184);
700 vRand = G4ThreeVector( -0.280469198715188, 0.4463870649961534, 0.8497503261406731 );
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 // assert(ApproxEqual(Dist,8.9849541128705734394));
708
709 pRand = G4ThreeVector(2.469342694407535, -0.5746362009323557, -0.04425421860130281);
710 vRand = G4ThreeVector( -0.2513630044708909, 0.4396138160169732, -0.8622971255607673);
711 norm = sn17.SurfaceNormal(pRand);
712 G4cout<<"norm = sn17.SurfaceNormal(pRand) = "<<norm.x()<<", "<<norm.y()<<", "<<norm.z()<<G4endl;
713 inside = sn17.Inside(pRand);
714 G4cout<<"sn17.Inside(pRand) = "<<OutputInside(inside)<<G4endl ;
715 Dist = sn17.DistanceToOut(pRand,vRand,calcNorm,pgoodNorm,pNorm);
716 G4cout<<"sn17.DistanceToOut(pRand,vRand,...) = "<<Dist<<G4endl;
717 Dist = sn17.DistanceToIn(pRand,vRand);
718 G4cout<<"sn17.DistanceToIn(pRand,vRand) = "<<Dist<<G4endl;
719 // assert(ApproxEqual(Dist,8.9849541128705734394));
720
721
722
723
724// Checking G4Sphere::DistanceToIn(P)
725
726 Dist=s2.DistanceToIn(pzero);
727 assert(ApproxEqual(Dist,45));
728 Dist=s1.DistanceToIn(ponrmax1);
729 assert(ApproxEqual(Dist,0));
730
731// Checking G4Sphere::DistanceToIn(P,V)
732
733
734 zP = 3.;
735
736 for( i = 0; i < iMax; i++)
737 {
738 rRand = 0.5*rMax*G4UniformRand();
739 phi = twopi*G4UniformRand();
740 pRand = G4ThreeVector( rRand*std::cos(phi), rRand*std::sin(phi), zP);
741
742 cosTheta = std::cos( 180.*deg - 78.*deg*G4UniformRand() );
743 cosTheta = -0.99; // std::cos( 180.*deg - 78.*deg*G4UniformRand() );
744
745 sinTheta = std::sqrt((1.-cosTheta)*(1.+cosTheta));
746 phi = twopi*G4UniformRand();
747 vRand = G4ThreeVector( sinTheta*std::cos(phi), sinTheta*std::sin(phi), cosTheta);
748
749 Dist = sn17.DistanceToIn(pRand,vRand);
750
751 if( -Dist*cosTheta > 2.*zP )
752 {
753 G4cout<<"sn17.DistanceToIn(pRand,vRand), i = "<<i<<G4endl;
754 G4cout<<pRand.x()<<", "<<pRand.y()<<", "<<pRand.z()<<G4endl;
755 G4cout<<vRand.x()<<", "<<vRand.y()<<", "<<vRand.z()<<G4endl;
756 break;
757 }
758 }
759
760
761
762 pRand = G4ThreeVector( -10.276604989981144911, -4.7072500741022338389, 1);
763 vRand = G4ThreeVector( -0.28873162920537848164, 0.51647890054938394577, -0.80615357816219324061);
764 Dist = sn17.DistanceToIn(pRand,vRand);
765 G4cout<<"sn17.DistanceToIn(pRand,vRand) = "<<Dist<<G4endl;
766 // assert(ApproxEqual(Dist,8.9849541128705734394));
767
768
769
770
771 Dist=s1.DistanceToIn(ponzmax,vz);
772 G4cout<<"s1.DistanceToIn(ponzmax,vz) = "<<Dist<<G4endl;
773 Dist=s1.DistanceToIn(pbigy,vy);
774 assert(Dist==kInfinity);
775 Dist=s1.DistanceToIn(pbigy,vmy);
776 assert(ApproxEqual(Dist,50));
777
778 Dist=s2.DistanceToIn(pzero,vy);
779 G4cout<<"s2.DistanceToIn(pzero,vx) = "<<Dist<<G4endl;
780 assert(ApproxEqual(Dist,45));
781 Dist=s2.DistanceToIn(pzero,vmy);
782 assert(ApproxEqual(Dist,45));
783 Dist=s2.DistanceToIn(ponrmin1,vx);
784 // G4cout<<"s2.DistanceToIn(ponmin1,vx) = "<<Dist<<G4endl;
785 assert(Dist==0);
786 Dist=s2.DistanceToIn(ponrmin1,vmx);
787 assert(ApproxEqual(Dist,90));
788
789 Dist=s2.DistanceToIn(ponrmin2,vx);
790 assert(Dist==0);
791 Dist=s2.DistanceToIn(ponrmin2,vmx);
792 assert(ApproxEqual(Dist,90/std::sqrt(2.)));
793
794
795 Dist=s3.DistanceToIn(ptesttheta1,vmz);
796 // G4cout<<"s3.DistanceToIn(ptesttheta1,vmz) = "<<Dist<<G4endl;
797 assert(ApproxEqual(Dist,100-48/std::sqrt(2.)));
798 Dist=s3.DistanceToIn(pontheta1,vz);
799 // G4cout<<"s3.DistanceToIn(pontheta1,vz) = "<<Dist<<G4endl;
800 assert(Dist==kInfinity);
801 Dist=s3.DistanceToIn(pontheta1,vmz);
802 assert(Dist==0);
803 Dist=s3.DistanceToIn(pontheta2,vz);
804 // G4cout<<"s3.DistanceToIn(pontheta2,vz) = "<<Dist<<G4endl;
805 assert(Dist==0);
806 Dist=s3.DistanceToIn(pontheta2,vmz);
807 assert(Dist==kInfinity);
808 Dist=s32.DistanceToIn(pontheta1,vz);
809 // G4cout<<"s32.DistanceToIn(pontheta1,vz) = "<<Dist<<G4endl;
810 assert(Dist==0);
811 Dist=s32.DistanceToIn(pontheta1,vmz);
812 // G4cout<<"s32.DistanceToIn(pontheta1,vmz) = "<<Dist<<G4endl;
813 assert(Dist==kInfinity);
814 Dist=s33.DistanceToIn(pontheta2,vz);
815 // G4cout<<"s33.DistanceToIn(pontheta2,vz) = "<<Dist<<G4endl;
816 assert(Dist==kInfinity);
817 Dist=s33.DistanceToIn(pontheta2,vmz);
818 // G4cout<<"s33.DistanceToIn(pontheta2,vmz) = "<<Dist<<G4endl;
819 assert(Dist==0);
820
821 Dist=s4.DistanceToIn(pbigy,vmy);
822 assert(Dist==kInfinity);
823 Dist=s4.DistanceToIn(pbigz,vmz);
824 assert(ApproxEqual(Dist,50));
825 Dist=s4.DistanceToIn(pzero,vy);
826 assert(Dist==kInfinity);
827 Dist=s4.DistanceToIn(pzero,vx);
828 assert(ApproxEqual(Dist,45));
829
830 Dist=s4.DistanceToIn(ptestphi1,vx);
831 assert(ApproxEqual(Dist,100+45/std::sqrt(2.)));
832 Dist=s4.DistanceToIn(ponphi1,vmxmy);
833 assert(Dist==kInfinity);
834 Dist=s4.DistanceToIn(ponphi1,vxy);
835 // G4cout<<"s4.DistanceToIn(ponphi1,vxy) = "<<Dist<<G4endl;
836 assert(ApproxEqual(Dist,0));
837
838 Dist=s4.DistanceToIn(ptestphi2,vx);
839 assert(ApproxEqual(Dist,100+45/std::sqrt(2.)));
840 Dist=s4.DistanceToIn(ponphi2,vmxy);
841 assert(Dist==kInfinity);
842 Dist=s4.DistanceToIn(ponphi2,vxmy);
843 // G4cout<<"s4.DistanceToIn(ponphi2,vxmy) = "<<Dist<<G4endl;
844 assert(ApproxEqual(Dist,0));
845
846 Dist=s3.DistanceToIn(pzero,vx);
847 assert(ApproxEqual(Dist,45));
848 Dist=s3.DistanceToIn(ptesttheta1,vmz);
849 assert(ApproxEqual(Dist,100-48/std::sqrt(2.)));
850 Dist=b216.DistanceToIn(p216,v216);
851 // G4cout<<"b216.DistanceToIn(p216,v216) = "<<Dist<<G4endl;
852
853 Dist=b658.DistanceToIn(pzero,vz);
854 G4cout<<"b658.DistanceToIn(pzero,vz) = "<<Dist<<G4endl;
855
856
857
858 Dist=s13.DistanceToIn(G4ThreeVector(20.,0.,-70.),vz);
859 // G4cout<<"s13.DistanceToIn(G4ThreeVector(20.,0.,-70.),vz... = "<<Dist<<G4endl;
860 assert(ApproxEqual(Dist,58.452994616207498));
861
862 Dist=s13.DistanceToIn(G4ThreeVector(20.,0.,70.),vmz);
863 // G4cout<<"s13.DistanceToIn(G4ThreeVector(20.,0.,70.),vmz... = "<<Dist<<G4endl;
864 assert(ApproxEqual(Dist,35.358983848622451));
865
866
867 Dist=s14.DistanceToIn(G4ThreeVector(20.,0.,-70.),vz);
868 // G4cout<<"s14.DistanceToIn(G4ThreeVector(20.,0.,-70.),vz... = "<<Dist<<G4endl;
869 assert(ApproxEqual(Dist,35.358983848622451));
870
871 Dist=s14.DistanceToIn(G4ThreeVector(20.,0.,70.),vmz);
872 // G4cout<<"s14.DistanceToIn(G4ThreeVector(20.,0.,70.),vmz... = "<<Dist<<G4endl;
873 assert(ApproxEqual(Dist,58.452994616207498));
874
875 Dist=b1046.DistanceToIn(G4ThreeVector(0.,0.,4800*km),vmz);
876 G4cout<<"b1046.DistanceToIn(G4ThreeVector(0.,0.,4800*km),vmz... = "<<Dist<<G4endl;
877 // assert(ApproxEqual(Dist,0.));
878
879
880
881
882
883 ///////////////////////////////////////////////////////////////////////////
884
885
886// CalculateExtent
887 G4VoxelLimits limit; // Unlimited
888 G4RotationMatrix noRot;
889 G4AffineTransform origin;
890 G4double min,max;
891 assert(s1.CalculateExtent(kXAxis,limit,origin,min,max));
892 assert(min<=-50&&max>=50);
893 assert(s1.CalculateExtent(kYAxis,limit,origin,min,max));
894 assert(min<=-50&&max>=50);
895 assert(s1.CalculateExtent(kZAxis,limit,origin,min,max));
896 assert(min<=-50&&max>=50);
897
898 assert(s2.CalculateExtent(kXAxis,limit,origin,min,max));
899 assert(min<=-50&&max>=50);
900 assert(s2.CalculateExtent(kYAxis,limit,origin,min,max));
901 assert(min<=-50&&max>=50);
902 assert(s2.CalculateExtent(kZAxis,limit,origin,min,max));
903 assert(min<=-50&&max>=50);
904
905 assert(s3.CalculateExtent(kXAxis,limit,origin,min,max));
906 assert(min<=-50&&max>=50);
907 assert(s3.CalculateExtent(kYAxis,limit,origin,min,max));
908 assert(min<=-50&&max>=50);
909 assert(s3.CalculateExtent(kZAxis,limit,origin,min,max));
910 assert(min<=-50/std::sqrt(2.)&&max>=50/std::sqrt(2.));
911
912 G4ThreeVector pmxmymz(-100,-110,-120);
913 G4AffineTransform tPosOnly(pmxmymz);
914 assert(s1.CalculateExtent(kXAxis,limit,tPosOnly,min,max));
915 assert(min<=-150&&max>=-50);
916 assert(s1.CalculateExtent(kYAxis,limit,tPosOnly,min,max));
917 assert(min<=-160&&max>=-60);
918 assert(s1.CalculateExtent(kZAxis,limit,tPosOnly,min,max));
919 assert(min<=-170&&max>=-70);
920
921 assert(s3.CalculateExtent(kXAxis,limit,tPosOnly,min,max));
922 assert(min<=-150&&max>=-50);
923 assert(s3.CalculateExtent(kYAxis,limit,tPosOnly,min,max));
924 assert(min<=-160&&max>=-60);
925 assert(s3.CalculateExtent(kZAxis,limit,tPosOnly,min,max));
926 assert(min<=-170+50/std::sqrt(2.)&&max>=-70-50/std::sqrt(2.));
927
928 G4RotationMatrix r90Z;
929 r90Z.rotateZ(halfpi);
930 G4AffineTransform tRotZ(r90Z,pzero);
931 assert(s1.CalculateExtent(kXAxis,limit,tRotZ,min,max));
932 assert(min<=-50&&max>=50);
933 assert(s1.CalculateExtent(kYAxis,limit,tRotZ,min,max));
934 assert(min<=-50&&max>=50);
935 assert(s1.CalculateExtent(kZAxis,limit,tRotZ,min,max));
936 assert(min<=-50&&max>=50);
937
938 assert(s2.CalculateExtent(kXAxis,limit,tRotZ,min,max));
939 assert(min<=-50&&max>=50);
940 assert(s2.CalculateExtent(kYAxis,limit,tRotZ,min,max));
941 assert(min<=-50&&max>=50);
942 assert(s2.CalculateExtent(kZAxis,limit,tRotZ,min,max));
943 assert(min<=-50&&max>=50);
944
945 assert(s3.CalculateExtent(kXAxis,limit,tRotZ,min,max));
946 assert(min<=-50&&max>=50);
947 assert(s3.CalculateExtent(kYAxis,limit,tRotZ,min,max));
948 assert(min<=-50&&max>=50);
949 assert(s3.CalculateExtent(kZAxis,limit,tRotZ,min,max));
950 assert(min<=-50/std::sqrt(2.)&&max>=50/std::sqrt(2.));
951
952// Check that clipped away
953 G4VoxelLimits xClip;
954 xClip.AddLimit(kXAxis,-100,-60);
955 assert(!s1.CalculateExtent(kXAxis,xClip,origin,min,max));
956
957// Assert clipped to volume
958 G4VoxelLimits allClip;
959 allClip.AddLimit(kXAxis,-5,+5);
960 allClip.AddLimit(kYAxis,-5,+5);
961 allClip.AddLimit(kZAxis,-5,+5);
962 G4RotationMatrix genRot;
963 genRot.rotateX(pi/6);
964 genRot.rotateY(pi/6);
965 genRot.rotateZ(pi/6);
966 G4AffineTransform tGen(genRot,vx);
967 assert(s1.CalculateExtent(kXAxis,allClip,tGen,min,max));
968 assert(min<=-5&&max>=5);
969 assert(s1.CalculateExtent(kYAxis,allClip,tGen,min,max));
970 assert(min<=-5&&max>=5);
971 assert(s1.CalculateExtent(kZAxis,allClip,tGen,min,max));
972 assert(min<=-5&&max>=5);
973
974 assert(s2.CalculateExtent(kXAxis,allClip,tGen,min,max));
975 assert(min<=-5&&max>=5);
976 assert(s2.CalculateExtent(kYAxis,allClip,tGen,min,max));
977 assert(min<=-5&&max>=5);
978 assert(s2.CalculateExtent(kZAxis,allClip,tGen,min,max));
979 assert(min<=-5&&max>=5);
980
981 s34.CalculateExtent(kXAxis,allClip,tGen,min,max);
982 // G4cout<<"s34.CalculateExtent(kXAxis,allClip,tGen,min,max)"<<G4endl ;
983 // G4cout<<"min = "<<min<<" max = "<<max<<G4endl ;
984
985 s34.CalculateExtent(kYAxis,allClip,tGen,min,max);
986 // G4cout<<"s3.CalculateExtent(kYAxis,allClip,tGen,min,max)"<<G4endl ;
987 // G4cout<<"min = "<<min<<" max = "<<max<<G4endl ;
988
989 s34.CalculateExtent(kZAxis,allClip,tGen,min,max);
990 // G4cout<<"s34.CalculateExtent(kZAxis,allClip,tGen,min,max)"<<G4endl ;
991 // G4cout<<"min = "<<min<<" max = "<<max<<G4endl ;
992 assert(s34.CalculateExtent(kXAxis,allClip,tGen,min,max));
993 assert(min<=-5&&max>=5);
994 assert(s34.CalculateExtent(kYAxis,allClip,tGen,min,max));
995 assert(min<=-5&&max>=5);
996 assert(s34.CalculateExtent(kZAxis,allClip,tGen,min,max));
997 assert(min<=-5&&max>=5);
998
999// Test z clipping ok
1000 for (G4double zTest=-100;zTest<100;zTest+=9)
1001 {
1002 G4VoxelLimits zTestClip;
1003 zTestClip.AddLimit(kZAxis,-kInfinity,zTest);
1004 if (zTest<-50)
1005 {
1006 assert(!s1.CalculateExtent(kZAxis,zTestClip,origin,min,max));
1007 assert(!s2.CalculateExtent(kZAxis,zTestClip,origin,min,max));
1008 }
1009 else
1010 {
1011 assert(s1.CalculateExtent(kZAxis,zTestClip,origin,min,max));
1012 assert(s2.CalculateExtent(kZAxis,zTestClip,origin,min,max));
1013 G4double testMin=-50;
1014 G4double testMax=(zTest<50) ? zTest : 50;
1015 assert (ApproxEqual(min,testMin)
1016 &&ApproxEqual(max,testMax));
1017 }
1018 }
1019
1020// Test y clipping ok
1021 for (G4double xTest=-100;xTest<100;xTest+=9)
1022 {
1023 G4VoxelLimits xTestClip;
1024 xTestClip.AddLimit(kXAxis,-kInfinity,xTest);
1025 if (xTest<-50)
1026 {
1027 assert(!s1.CalculateExtent(kYAxis,xTestClip,origin,min,max));
1028 }
1029 else
1030 {
1031 assert(s1.CalculateExtent(kYAxis,xTestClip,origin,min,max));
1032// Calc max y coordinate
1033 G4double testMax=(xTest<0) ? std::sqrt(50*50-xTest*xTest) : 50;
1034 assert (ApproxEqual(min,-testMax)
1035 &&ApproxEqual(max,testMax));
1036 }
1037 }
1038
1039// Test x clipping ok
1040 for (G4double yTest=-100;yTest<100;yTest+=9)
1041 {
1042 G4VoxelLimits yTestClip;
1043 yTestClip.AddLimit(kYAxis,-kInfinity,yTest);
1044 if (yTest<-50)
1045 {
1046 assert(!s1.CalculateExtent(kXAxis,yTestClip,origin,min,max));
1047 }
1048 else
1049 {
1050 assert(s1.CalculateExtent(kXAxis,yTestClip,origin,min,max));
1051// Calc max y coordinate
1052 G4double testMax=(yTest<0) ? std::sqrt(50*50-yTest*yTest) : 50;
1053 assert (ApproxEqual(min,-testMax)
1054 &&ApproxEqual(max,testMax));
1055 }
1056 }
1057
1058 G4bool checkPoint( const G4Sphere& pSph, G4ThreeVector origin,
1059 G4double d, G4ThreeVector dir, EInside exp);
1060
1061 G4Sphere SpAroundX("SpAroundX", 10.*mm, 1000.*mm, -1.0*degree, 2.0*degree, 0.*degree, 180.0*degree );
1062
1063 G4double sinOneDeg = std::sin( 1.0 * degree );
1064 radOne = 100.0 * mm;
1065
1066 G4ThreeVector ptPhiSurfExct= G4ThreeVector( radOne * std::cos( -1.0 * degree ) ,
1067 - radOne * sinOneDeg,
1068 0.0 );
1069 G4cout << " Starting from point " << ptPhiSurfExct << G4endl;
1070 G4cout << " using direction " << vy << G4endl;
1071 checkPoint( SpAroundX, ptPhiSurfExct, -radOne * kAngTolerance * 1.5,
1072 vy, kOutside);
1073 checkPoint( SpAroundX, ptPhiSurfExct, -radOne * kAngTolerance * 0.49,
1074 vy, kSurface);
1075 checkPoint( SpAroundX, ptPhiSurfExct, -radOne * kAngTolerance * 0.25,
1076 vy, kSurface);
1077 checkPoint( SpAroundX, ptPhiSurfExct, 0.0,
1078 vy, kSurface);
1079 checkPoint( SpAroundX, ptPhiSurfExct, radOne * kAngTolerance * 0.25,
1080 vy, kSurface);
1081 checkPoint( SpAroundX, ptPhiSurfExct, radOne * kAngTolerance * 0.49,
1082 vy, kSurface);
1083 checkPoint( SpAroundX, ptPhiSurfExct, radOne * kAngTolerance * 1.5,
1084 vy, kInside);
1085
1086 // Try one that has a 'deep' negative phi section
1087 // --> Vlad. Grichine test case, 30 Oct 2003
1088 //
1089 G4cout << G4endl << G4endl << "" << G4endl;
1090 G4cout << "========================================================= " << G4endl;
1091
1092 G4Sphere SphDeepNeg("DeepNegPhiSphere", 10.*mm, 1000.*mm,
1093 -270.0*degree, 280.0*degree, // start Phi, delta Phi
1094 0.*degree, 180.0*degree ); // start Theta, delta Theta
1095 G4double phiPoint = 160.0 * degree;
1096 G4ThreeVector StartPt( radOne * std::cos(phiPoint), radOne * std::sin(phiPoint), 0.0);
1097 G4cout << "For sphere " << SphDeepNeg.GetName() << G4endl;
1098 G4cout << " Starting from point " << ptPhiSurfExct << G4endl;
1099
1100 checkPoint( SphDeepNeg, StartPt, 0.0, vy, kInside);
1101
1102 // Try the edges
1103 G4ThreeVector NegEdgePt( radOne * std::cos(-270.0*degree), radOne * std::sin(-270.0*degree), 0.0);
1104 G4ThreeVector PosEdgePt( radOne * std::cos(10.0*degree), radOne * std::sin(10.0*degree), 0.0);
1105
1106 G4cout << "--------------------------------------------------------" << G4endl;
1107 G4cout << " New point " << NegEdgePt << " should be at Neg edge of -270.0 degrees " << G4endl;
1108 checkPoint( SphDeepNeg, NegEdgePt, 0.0, -vx, kSurface);
1109 checkPoint( SphDeepNeg, NegEdgePt, radOne*kAngTolerance * 0.25, -vx, kSurface);
1110 checkPoint( SphDeepNeg, NegEdgePt, -radOne*kAngTolerance * 0.25, -vx, kSurface);
1111 checkPoint( SphDeepNeg, NegEdgePt, radOne*kAngTolerance * 1.25, -vx, kInside);
1112 checkPoint( SphDeepNeg, NegEdgePt, -radOne*kAngTolerance * 1.25, -vx, kOutside);
1113
1114 G4cout << "--------------------------------------------------------" << G4endl;
1115 G4cout << " New point " << PosEdgePt << " should be at Pos edge of +10.0 degrees " << G4endl;
1116 checkPoint( SphDeepNeg, PosEdgePt, 0.0, -vy, kSurface);
1117 checkPoint( SphDeepNeg, PosEdgePt, radOne*kAngTolerance * 0.25, -vy, kSurface);
1118 checkPoint( SphDeepNeg, PosEdgePt, -radOne*kAngTolerance * 0.25, -vy, kSurface);
1119 checkPoint( SphDeepNeg, PosEdgePt, -radOne*kAngTolerance * 1.25, -vy, kOutside);
1120 checkPoint( SphDeepNeg, PosEdgePt, radOne*kAngTolerance * 1.25, -vy, kInside);
1121
1122
1123
1124 // G4double sTheta = 135*degree, pxt = -10., pyt = 10.;
1125
1126 // G4cout<<"sTheta = "<<sTheta/degree<<" degree"<<G4endl;
1127 // G4cout<<"std::tan(sTheta) = "<<std::tan(sTheta)<<G4endl;
1128 // G4cout<<"std::sin(sTheta) = "<<std::sin(sTheta)<<G4endl;
1129 // G4cout<<"std::atan(sTheta) = "<<std::atan(sTheta)/degree<<G4endl;
1130 // G4cout<<"std::atan2(pyt,pxt) = pyt/pxt = "<<std::atan2(pyt,pxt)/degree<<G4endl;
1131
1132 return 0;
1133}
1134
1135// Given the sphere 'rSphere', a point 'origin'
1136// --> --> -->
1137// the function below checks that the point pnew=( origin + dist * direction )
1138// - the point (p+dist*dir) is located in the part of the solid given by 'expectedInResult'
1139// - and that from there, the DistanceToIn along 'dir' is not negative or Infinite
1140//
1141// Use cases expected:
1142// - 'origin' is on/near a surface and 'direction' is pointing towards the inside of the solid
1143
1144G4bool
1145checkPoint( const G4Sphere &rSphere,
1146 G4ThreeVector origin,
1147 G4double dist,
1148 G4ThreeVector direction,
1149 EInside expectedInResult)
1150{
1151 G4int verbose = 0;
1152
1153 G4ThreeVector newPoint;
1154 G4double distIn=-1.0, distOut=-1.0;
1155
1156 newPoint = origin + dist * direction;
1157
1158 G4int oldPrecision= G4cout.precision(10);
1159 // G4cout << " --- Sphere " << rSphere.GetName() << "" << G4endl;
1160 if( verbose > 0 ) {
1161 G4cout << G4endl;
1162 if (verbose > 2 ) G4cout << " Sphere " << rSphere.GetName();
1163 G4cout.precision(10);
1164 if (verbose > 1 ) G4cout << " dir= " << direction;
1165 G4cout << " dist= " << dist;
1166 }
1167
1168 EInside inSphere= rSphere.Inside( newPoint ) ;
1169 /*======*/
1170 G4cout.precision(15);
1171 // G4cout << " NewPoint " << newPoint << " is "
1172 G4bool goodIn= (inSphere == expectedInResult) ;
1173 if ( !goodIn ) {
1174 G4cout << " ************ Unexpected Result for Inside *************** " << G4endl;
1175 }
1176 if ( verbose || !goodIn ) {
1177 G4cout << " New point "
1178 << " is " << OutputInside( inSphere )
1179 << " vs " << OutputInside( expectedInResult )
1180 << " expected." << G4endl ;
1181 }
1182
1183 G4bool goodDistIn = true;
1184
1185 distIn = rSphere.DistanceToIn( newPoint, direction );
1186 /*===========*/
1187 if ( verbose ) G4cout << " DistToIn (p, dir) = " << distIn << G4endl;
1188 if( (inSphere == kOutside)
1189 && (distIn < 0.0 ) // Cannot use 0.5*kCarTolerance for Angular tolerance!!
1190 ){
1191 G4cout << " ********** Unexpected Result for DistanceToIn from outside ********* " << G4endl;
1192 // G4cout << " It should be " << G4endl;
1193 goodDistIn = false;
1194 }
1195 if( (inSphere == kSurface )
1196 && ( (distIn < 0.0) || (distIn >= kInfinity ))
1197 ){
1198 G4cout << " ********** Unexpected Result for DistanceToIn on surface ********* " << G4endl;
1199 // if ( (distIn != 0.0) )
1200 // - Can check that the return value must be 0.0
1201 // But in general case the direction can be away from the solid,
1202 // and then a finite or kInfinity answer is correct
1203 // --> must check the direction against the normal
1204 // in order to perform this check in general case.
1205
1206 goodDistIn = false;
1207 }
1208 if ( verbose || !goodDistIn ) {
1209 G4cout << " DistToIn (p, dir) = " << distIn << G4endl;
1210 }
1211
1212 G4bool good= (goodIn && goodDistIn);
1213 if ( !good ){
1214 // There was an error -- document the use case!
1215 G4cout << " --- Sphere " << rSphere.GetName() << "" << G4endl;
1216 G4cout << " Origin= " << origin << G4endl;
1217 G4cout << " Direction= " << direction << G4endl;
1218 G4cout << " dist= " << dist;
1219 G4cout << " Actual-point= " << newPoint << G4endl;
1220 }
1221
1222 distOut = rSphere.DistanceToOut( newPoint, direction );
1223 /*=============*/
1224 if ( verbose ) G4cout << " DistToOut (p, dir) = " << distOut << G4endl;
1225
1226 G4cout.precision(oldPrecision);
1227
1228 return good;
1229}
1230
Note: See TracBrowser for help on using the repository browser.