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

Last change on this file since 1358 was 1347, checked in by garnier, 14 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.