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

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

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