source: trunk/source/geometry/solids/Boolean/test/testG4SubtractionSolid.cc@ 1350

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

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
File size: 23.6 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
28#include <assert.h>
29#include <cmath>
30
31#include "globals.hh"
32#include "geomdefs.hh"
33
34#include "Randomize.hh"
35
36#include "ApproxEqual.hh"
37
38#include "G4ThreeVector.hh"
39#include "G4RotationMatrix.hh"
40#include "G4AffineTransform.hh"
41#include "G4VoxelLimits.hh"
42
43#include "G4Box.hh"
44#include "G4Cons.hh"
45#include "G4Para.hh"
46#include "G4Sphere.hh"
47#include "G4Orb.hh"
48#include "G4Ellipsoid.hh"
49#include "G4Torus.hh"
50#include "G4Trap.hh"
51#include "G4Trd.hh"
52#include "G4Tubs.hh"
53
54#include "G4SubtractionSolid.hh"
55#include "G4IntersectionSolid.hh"
56
57
58// #include "G4DisplacedSolid.hh"
59
60///////////////////////////////////////////////////////////////////
61//
62// Dave's auxiliary function
63
64const G4String OutputInside(const EInside a)
65{
66 switch(a)
67 {
68 case kInside: return "Inside";
69 case kOutside: return "Outside";
70 case kSurface: return "Surface";
71 }
72 return "????";
73}
74
75
76int main()
77{
78 G4int i;
79 G4ThreeVector pzero(0,0,0), p;
80
81 G4ThreeVector ponxside(20,0,0),ponyside(0,30,0),ponzside(0,0,40),
82
83 ponb2x(10,0,0), ponb2y(0,10,0), ponb2z(0,0,10),
84
85 ponb2mx(-10,0,0),ponb2my(0,-10,0),ponb2mz(0,0,-10);
86
87 G4ThreeVector ponmxside(-20,0,0),ponmyside(0,-30,0),ponmzside(0,0,-40);
88
89 G4ThreeVector ponzsidey(0,25,40),ponmzsidey(0,25,-40),
90
91 ponb2zy(0,5,10),ponb2mzy(0,5,-10) ;
92
93 G4ThreeVector pbigx(100,0,0),pbigy(0,100,0),pbigz(0,0,100);
94
95 G4ThreeVector pbigmx(-100,0,0),pbigmy(0,-100,0),pbigmz(0,0,-100);
96
97 G4ThreeVector vx(1,0,0),vy(0,1,0),vz(0,0,1);
98
99 G4ThreeVector vmx(-1,0,0),vmy(0,-1,0),vmz(0,0,-1);
100
101 G4ThreeVector vxy(1/std::sqrt(2.0),1/std::sqrt(2.0),0);
102
103 G4ThreeVector vmxy(-1/std::sqrt(2.0),1/std::sqrt(2.0),0);
104
105 G4ThreeVector vmxmy(-1/std::sqrt(2.0),-1/std::sqrt(2.0),0);
106
107 G4ThreeVector vxmy(1/std::sqrt(2.0),-1/std::sqrt(2.0),0);
108
109 G4ThreeVector vxmz(1/std::sqrt(2.0),0,-1/std::sqrt(2.0));
110
111 G4double dist;
112 G4ThreeVector *pNorm,norm;
113 G4bool *pgoodNorm,goodNorm,calcNorm=true;
114
115 pNorm=&norm;
116 pgoodNorm=&goodNorm;
117
118 G4RotationMatrix identity, xRot ;
119
120// NOTE: xRot = rotation such that x axis->-x axis & y axis->-y axis
121
122 xRot.rotateZ(-pi);
123
124 G4Transform3D transform(xRot,G4ThreeVector(0,30,0)) ;
125 G4Transform3D transform2(xRot,G4ThreeVector(50,0,0)) ;
126
127 G4Box b1("Test Box #1",20.,30.,40.);
128 G4Box b2("Test Box #2",10.,10.,10.);
129 G4Box b3("Test Box #3",10.,20.,50.);
130 G4Box b4("Test Box #4",20.,20.,40.);
131 G4Box b5("Test Box #4",50.,50.,50.);
132
133 G4Tubs t1("Solid Tube #1",0,50.,50.,0.,360.*degree);
134
135 // t2\t3 for DistanceToIn
136
137 G4Tubs t2("Hole Tube #2",50.,60.,50.,0.,2.*pi);
138
139 G4Tubs t3("Hole Tube #3",45.,55.,50.,pi/4.,pi*3./2.);
140
141 G4Cons c1("Hollow Full Tube",50.,100.,50.,100.,50.,0,2*pi),
142
143 c2("Full Cone",0,50.,0,100.,50.,0,2*pi) ;
144
145 G4SubtractionSolid b1Sb2("b1Sb2",&b1,&b2),
146
147 t1Sb2("t1Sb2",&t1,&b2),
148
149 t2St3("t2St3",&t2,&t3),
150
151 c2Sb2("c2Sb2",&c2,&b2) ;
152
153 // Tube t1 \ box t3, which was rotated by -pi and translated +30y
154
155 G4SubtractionSolid t1Sb3("t1Subtractionb3",&t1,&b3,transform) ;
156 G4SubtractionSolid b1Sb4("t1Subtractionb3",&b1,&b4,transform) ;
157 G4SubtractionSolid b5St1("b4St1",&b5,&t1,transform2) ;
158
159 G4SubtractionSolid b1Sb2touch("b1Sb4touch",&b1,&b2,&identity,G4ThreeVector(10.,0,0)) ;
160
161
162
163
164
165 G4Tubs* tube4 = new G4Tubs( "OuterFrame",
166 1.0*m, // inner radius
167 1.1*m, // outer radius
168 0.01*m, // half-thickness in z
169 -15*deg, // start angle
170 30*deg ); // total angle
171
172 G4Box* tube5 = new G4Box( "Cutout", // name
173 0.02*m, // half-width (x)
174 0.25*m, // half-height (y)
175 0.01001*m ); // half-thickness (z)
176
177 G4Transform3D tran2 = G4Translate3D( 1.03*m, 0.0, 0.0 );
178
179 G4VSolid* solid = new G4SubtractionSolid( "drcExample",tube4,tube5, tran2 );
180
181 G4Cons* cone3 = new G4Cons( "OuterFrame",
182 0.6*m, // pRmin1
183 1.0*m, // pRmax1
184 0.2*m, // pRmin2
185 0.8*m, // pRmax2
186 0.2*m,
187 0*deg,
188 180*deg );
189
190 G4Cons* cone4 = new G4Cons( "OuterFrame",
191 0.6*m, // pRmin1
192 1.0*m, // pRmax1
193 0.2*m, // pRmin2
194 0.8*m, // pRmax2
195 0.2*m,
196 0*deg,
197 180*deg );
198 G4RotationMatrix rotmat3;
199 rotmat3.rotateY(pi/4.0);
200 G4Transform3D tran3 = G4Transform3D(rotmat3,G4ThreeVector(0.0,0.0,0.0));
201
202 G4VSolid* c3Ic4 = new G4SubtractionSolid( "Example", cone3, cone4, tran3 );
203
204
205 //G4Torus* insp = new G4Torus("Isp", 7.5*cm, 8.1*cm, 15.6*cm, 0, 2*pi);
206 //G4Tubs* tmptube = new G4Tubs("TMPT", 7.5*cm, 15.6*cm, 8.1*cm, 0, 2*pi);
207
208 G4RotationMatrix* rotDz5 = new G4RotationMatrix();
209 rotDz5->rotateZ(pi/2.);
210
211 // G4SubtractionSolid* inn = new G4SubtractionSolid("Innerpart", insp,
212 // tmptube, rotDz5, G4ThreeVector(0,0,0));
213
214 // LHCB RICH-detector problem:
215
216 G4ThreeVector lhcbTransl(8606.,173396.33,100110.42);
217 G4RotationMatrix lhcbMat;
218 lhcbMat.rotateX(-60*degree);
219 G4Transform3D lhcbTran = G4Transform3D(lhcbMat,lhcbTransl);
220 G4Sphere* lhcbSphere = new G4Sphere("lhcbSphere",8600*mm,
221 8606*mm,
222 -1.6991355*degree,
223 3.3982711*degree,
224 88.528559*degree,
225 2.9428812*degree );
226
227 G4Box* lhcbBox = new G4Box("lhcbBox",100000*mm,
228 100000*mm,
229 200000*mm );
230
231
232 G4VSolid* lhcbSub = new G4SubtractionSolid( "lhcbSub", lhcbSphere,
233 lhcbBox, lhcbTran );
234
235 G4ThreeVector lhcbP(7298.73956059975,-394.9290932077643,1497.43923688271);
236 G4ThreeVector lhcbV(0.6062807093945133,0.06537967965081047,-0.7925586406726276);
237
238
239 G4Orb* largeOrb = new G4Orb("largeOrb",11.*mm);
240 G4Orb* smallOrb = new G4Orb("smallOrb",10.*mm);
241
242 G4VSolid* orbShell = new G4SubtractionSolid( "orbShell", largeOrb,
243 smallOrb, 0,0 );
244
245 G4Ellipsoid* largeEll = new G4Ellipsoid("largeEll",11.*mm, 11*mm, 11*mm); //, -11*mm, 11*mm);
246 G4Ellipsoid* smallEll = new G4Ellipsoid("smallEll",10.*mm, 10*mm, 10*mm); //, -10*mm, 10*mm);
247
248 G4VSolid* ellShell = new G4SubtractionSolid( "orbShell", largeEll,
249 smallEll, 0,0 );
250
251 G4Box* orbBox = new G4Box("orbBox",20*mm, 20*mm, 20*mm );
252
253 G4RotationMatrix orbMat;
254 G4ThreeVector orbTransl(0.,0.,-20.);
255 G4Transform3D orbTran = G4Transform3D(orbMat, orbTransl);
256
257 G4VSolid* orbMirror = new G4IntersectionSolid( "orbMirror", orbShell,
258 orbBox, orbTran);
259
260 G4VSolid* ellMirror = new G4IntersectionSolid( "orbMirror", ellShell,
261 orbBox, orbTran);
262
263 G4cout.precision(16) ;
264
265// Check Inside
266
267 assert(b1Sb2.Inside(pzero)==kOutside);
268 assert(b1Sb2.Inside(pbigz)==kOutside);
269 assert(b1Sb2.Inside(ponb2x)==kSurface);
270 assert(b1Sb2.Inside(ponb2y)==kSurface);
271 assert(b1Sb2.Inside(ponb2z)==kSurface);
272
273 assert(t1Sb3.Inside(pzero)==kInside);
274 assert(t1Sb3.Inside(pbigz)==kOutside);
275 assert(t1Sb3.Inside(ponb2x)==kInside);
276 assert(t1Sb3.Inside(ponb2y)==kSurface);
277 assert(t1Sb3.Inside(ponb2z)==kInside);
278
279 assert(c2Sb2.Inside(pzero)==kOutside);
280 assert(c2Sb2.Inside(pbigz)==kOutside);
281 assert(c2Sb2.Inside(ponb2x)==kSurface);
282 assert(c2Sb2.Inside(ponb2y)==kSurface);
283 assert(c2Sb2.Inside(ponb2z)==kSurface);
284
285 assert(b1Sb2touch.Inside(G4ThreeVector(20.,0,0))==kOutside);
286 assert(b1Sb2touch.Inside(G4ThreeVector(20.,10.,0))==kSurface);
287
288
289
290// Check Surface Normal
291
292 G4ThreeVector normal;
293
294 normal=b1Sb2.SurfaceNormal(ponb2x);
295 assert(ApproxEqual(normal,G4ThreeVector(-1,0,0)));
296
297 normal=b1Sb2.SurfaceNormal(ponb2mx);
298 assert(ApproxEqual(normal,G4ThreeVector(1,0,0)));
299
300 normal=b1Sb2.SurfaceNormal(ponb2y);
301 assert(ApproxEqual(normal,G4ThreeVector(0,-1,0)));
302
303 normal=b1Sb2.SurfaceNormal(ponb2my);
304 assert(ApproxEqual(normal,G4ThreeVector(0,1,0)));
305
306 normal=b1Sb2.SurfaceNormal(ponb2z);
307 assert(ApproxEqual(normal,G4ThreeVector(0,0,-1)));
308
309 normal=b1Sb2.SurfaceNormal(ponb2mz);
310 assert(ApproxEqual(normal,G4ThreeVector(0,0,1)));
311
312 normal=b1Sb2.SurfaceNormal(ponb2zy);
313 assert(ApproxEqual(normal,G4ThreeVector(0,0,-1)));
314
315 normal=b1Sb2.SurfaceNormal(ponb2mzy);
316 assert(ApproxEqual(normal,G4ThreeVector(0,0,1)));
317
318
319 normal=t1Sb3.SurfaceNormal(ponb2y);
320 assert(ApproxEqual(normal,G4ThreeVector(0,1,0)));
321
322
323// DistanceToOut(P)
324
325 dist=b1Sb2.DistanceToOut(ponb2x);
326 assert(ApproxEqual(dist,0));
327
328 dist=b1Sb2.DistanceToOut(ponb2y);
329 assert(ApproxEqual(dist,0));
330
331 dist=b1Sb2.DistanceToOut(ponb2z);
332 assert(ApproxEqual(dist,0));
333
334 dist=b1Sb2.DistanceToOut(ponxside);
335 assert(ApproxEqual(dist,0));
336
337 dist=t1Sb3.DistanceToOut(ponb2y);
338 assert(ApproxEqual(dist,0));
339
340 dist=t1Sb3.DistanceToOut(pzero);
341 assert(ApproxEqual(dist,10));
342
343 dist=t1Sb3.DistanceToOut(ponb2x);
344 assert(ApproxEqual(dist,10));
345
346 dist=t1Sb3.DistanceToOut(ponb2z);
347 assert(ApproxEqual(dist,10));
348
349
350// DistanceToOut(P,V)
351
352 dist=b1Sb2.DistanceToOut(ponb2x,vx,calcNorm,pgoodNorm,pNorm);
353 assert(ApproxEqual(dist,10)&&ApproxEqual(*pNorm,vx)&&*pgoodNorm);
354
355 dist=b1Sb2.DistanceToOut(ponxside,vmx,calcNorm,pgoodNorm,pNorm);
356 assert(ApproxEqual(dist,10)&&ApproxEqual(norm,vmx)); // &&*pgoodNorm);
357
358 dist=b1Sb2.DistanceToOut(ponb2y,vy,calcNorm,pgoodNorm,pNorm);
359 assert(ApproxEqual(dist,20)&&ApproxEqual(norm,vy)&&*pgoodNorm);
360
361 dist=b1Sb2.DistanceToOut(ponyside,vmy,calcNorm,pgoodNorm,pNorm);
362 assert(ApproxEqual(dist,20)&&ApproxEqual(norm,vmy)); // &&*pgoodNorm);
363
364 dist=b1Sb2.DistanceToOut(ponb2z,vz,calcNorm,pgoodNorm,pNorm);
365 assert(ApproxEqual(dist,30)&&ApproxEqual(norm,vz)&&*pgoodNorm);
366
367 dist=b1Sb2.DistanceToOut(ponzside,vmz,calcNorm,pgoodNorm,pNorm);
368 assert(ApproxEqual(dist,30)&&ApproxEqual(norm,vmz)); // &&*pgoodNorm);
369
370 dist=b1Sb2.DistanceToOut(ponb2x,vxy,calcNorm,pgoodNorm,pNorm);
371 assert(ApproxEqual(dist,std::sqrt(200.))&&*pgoodNorm);
372
373 dist=b1Sb2.DistanceToOut(ponb2x,vx,calcNorm,pgoodNorm,pNorm);
374 assert(ApproxEqual(dist,10)&&ApproxEqual(*pNorm,vx)&&*pgoodNorm);
375
376 dist=b1Sb2.DistanceToOut(ponb2x,vmx,calcNorm,pgoodNorm,pNorm);
377 assert(ApproxEqual(dist,0)&&ApproxEqual(*pNorm,vmx)); // &&*pgoodNorm);
378
379 dist=b1.DistanceToOut(ponxside,vy,calcNorm,pgoodNorm,pNorm);
380
381// G4cout<<"b1.DistanceToOut(ponxside,vy) = "<<dist<<G4endl;
382// assert(ApproxEqual(dist,0)&&ApproxEqual(*pNorm,vy)&&*pgoodNorm);
383
384 dist=b1Sb2.DistanceToOut(ponb2mx,vmx,calcNorm,pgoodNorm,pNorm);
385 assert(ApproxEqual(dist,10)&&ApproxEqual(norm,vmx)&&*pgoodNorm);
386
387 dist=b1Sb2.DistanceToOut(ponb2y,vy,calcNorm,pgoodNorm,pNorm);
388 assert(ApproxEqual(dist,20)&&ApproxEqual(norm,vy)&&*pgoodNorm);
389
390 dist=b1Sb2.DistanceToOut(ponb2my,vmy,calcNorm,pgoodNorm,pNorm);
391 assert(ApproxEqual(dist,20)&&ApproxEqual(norm,vmy)&&*pgoodNorm);
392
393 dist=b1Sb2.DistanceToOut(ponb2z,vz,calcNorm,pgoodNorm,pNorm);
394 assert(ApproxEqual(dist,30)&&ApproxEqual(norm,vz)&&*pgoodNorm);
395
396 dist=b1Sb2.DistanceToOut(ponb2mz,vmz,calcNorm,pgoodNorm,pNorm);
397 assert(ApproxEqual(dist,30)&&ApproxEqual(norm,vmz)&&*pgoodNorm);
398
399 // With placement
400
401 dist=t1Sb3.DistanceToOut(pzero,vy,calcNorm,pgoodNorm,pNorm);
402 assert(ApproxEqual(dist,10)&&ApproxEqual(norm,vy)); // &&*pgoodNorm);
403
404 dist=t1Sb3.DistanceToOut(pzero,vx,calcNorm,pgoodNorm,pNorm);
405 assert(ApproxEqual(dist,50)&&ApproxEqual(norm,vx)&&*pgoodNorm);
406
407 dist=t1Sb3.DistanceToOut(pzero,vmy,calcNorm,pgoodNorm,pNorm);
408 assert(ApproxEqual(dist,50)&&ApproxEqual(norm,vmy)&&*pgoodNorm);
409
410 dist=t1Sb3.DistanceToOut(pzero,vz,calcNorm,pgoodNorm,pNorm);
411 assert(ApproxEqual(dist,50)&&ApproxEqual(norm,vz)&&*pgoodNorm);
412
413
414//DistanceToIn(P)
415
416 dist=b1Sb2.DistanceToIn(pbigx);
417 assert(ApproxEqual(dist,80));
418
419 dist=b1Sb2.DistanceToIn(ponb2x);
420 assert(ApproxEqual(dist,0));
421
422 dist=b1Sb2.DistanceToIn(ponxside);
423 assert(ApproxEqual(dist,0));
424
425 dist=b1Sb2.DistanceToIn(pbigmy);
426 assert(ApproxEqual(dist,70));
427
428 dist=b1Sb2.DistanceToIn(pbigz);
429 assert(ApproxEqual(dist,60));
430
431 dist=b1Sb2.DistanceToIn(pbigmz);
432 assert(ApproxEqual(dist,60));
433
434 // With placement
435
436 dist=t1Sb3.DistanceToIn(pbigy);
437 assert(ApproxEqual(dist,50));
438
439 dist=t1Sb3.DistanceToIn(pbigmy);
440 assert(ApproxEqual(dist,50));
441
442 dist=t1Sb3.DistanceToIn(G4ThreeVector(0,20,0));
443 assert(ApproxEqual(dist,10));
444
445 dist=t1Sb3.DistanceToIn(G4ThreeVector(0,15,0));
446 assert(ApproxEqual(dist,5));
447
448// DistanceToIn(P,V)
449
450 dist=b1Sb2.DistanceToIn(pbigx,vmx);
451 assert(ApproxEqual(dist,80));
452
453 dist=b1Sb2.DistanceToIn(pbigmx,vx);
454 assert(ApproxEqual(dist,80));
455
456 dist=b1Sb2.DistanceToIn(pbigy,vmy);
457 assert(ApproxEqual(dist,70));
458
459 dist=b1Sb2.DistanceToIn(pbigmy,vy);
460 assert(ApproxEqual(dist,70));
461
462 dist=b1Sb2.DistanceToIn(pbigz,vmz);
463 assert(ApproxEqual(dist,60));
464
465 dist=b1Sb2.DistanceToIn(pbigmz,vz);
466 assert(ApproxEqual(dist,60));
467
468 dist=b1Sb2.DistanceToIn(pbigx,vxy);
469 assert(ApproxEqual(dist,kInfinity));
470
471 dist=b1Sb2.DistanceToIn(pbigmx,vxy);
472 assert(ApproxEqual(dist,kInfinity));
473
474 // Cases with a number of iterations inside 'while'
475
476 dist=t2St3.DistanceToIn(G4ThreeVector(2.5,-52.5,0),vy) ;
477 assert(ApproxEqual(dist,107.443));
478// G4cout<<"t2St3.DistanceToIn(G4ThreeVector(2.5,-52.5,0),vy) = "<<dist<<G4endl ;
479
480 dist=t2St3.DistanceToIn(G4ThreeVector(2.5,-62.5,0),vy) ;
481 assert(ApproxEqual(dist,2.55211));
482// G4cout<<"t2St3.DistanceToIn(G4ThreeVector(2.5,-62.5,0),vy) = "<<dist<<G4endl ;
483
484 // With placement
485
486 dist=t1Sb3.DistanceToIn(G4ThreeVector(0,100,0),vmy);
487 G4cout<<"t1Sb3.DistanceToIn(G4ThreeVector(0,100,0),vmy) = "<<dist<<G4endl; // 90?
488 // assert(ApproxEqual(dist,90)); // was 50
489
490 dist=t1Sb3.DistanceToIn(G4ThreeVector(0,36,0),vmy);
491 assert(ApproxEqual(dist,26));
492
493 dist=t1Sb3.DistanceToIn(G4ThreeVector(0,36,0),vy);
494 assert(ApproxEqual(dist,kInfinity));
495
496 dist=t1Sb3.DistanceToIn(G4ThreeVector(0,36,0),vx);
497 assert(ApproxEqual(dist,10));
498
499 dist=solid->DistanceToIn(
500 G4ThreeVector(1462.639634076807,-1183.908767887379,1926.359643354072),
501 G4ThreeVector(-0.1985206962900923,0.5639538694097884,-0.8015894000810041));
502 // G4cout<<"solid.DistanceToIn(p,v) = "<<dist<<G4endl ;
503 assert(ApproxEqual(dist,2390.7));
504
505 dist=solid->DistanceToIn(
506 G4ThreeVector(1572.774744032383,-1704.268426706159,1922.635294938558),
507 G4ThreeVector(-0.2113262591952278,0.6627445069650045,-0.7184086098191366));
508 // G4cout<<"solid.DistanceToIn(p,v) = "<<dist<<G4endl ;
509 assert(ApproxEqual(dist,2663.0611));
510
511 dist=c3Ic4->DistanceToIn(
512 G4ThreeVector(967.602560486646,-1211.626221342084,610.1637191085789),
513 G4ThreeVector(-0.1174412716483462,0.8263033300403027,-0.5508451274885945));
514 // assert(ApproxEqual(dist,kInfinity));
515 G4cout<<"c3Ic4->DistanceToIn = "<<dist<<G4endl ;
516
517 dist=lhcbSub->DistanceToIn(lhcbP,lhcbV);
518 // G4cout<<"lhcbSub->DistanceToIn(lhcbP,lhcbV) = "<<dist<<G4endl ;
519
520 G4ThreeVector orbP = G4ThreeVector(0.,0.,-9.999);
521
522
523 G4double orbTheta = 30*degree;
524 G4double orbPhi = 30*degree;
525
526 G4ThreeVector orbV = G4ThreeVector(std::sin(orbTheta)*std::cos(orbPhi),
527 std::sin(orbTheta)*std::sin(orbPhi),
528 std::cos(orbTheta));
529
530
531 EInside insideP = smallOrb->Inside(orbP);
532 G4cout<<"smallOrb->Inside(orbP) = "<<OutputInside(insideP)<<G4endl ;
533 dist = smallOrb->DistanceToOut(orbP,orbV);
534 G4cout<<"smallOrb->DistanceToOut(orbP,orbV) = "<<dist<<G4endl ;
535
536 insideP = orbShell->Inside(orbP);
537 G4cout<<"orbShell->Inside(orbP) = "<<OutputInside(insideP)<<G4endl ;
538 dist = orbShell->DistanceToIn(orbP,orbV);
539 G4cout<<"orbShell->DistanceToIn(orbP,orbV) = "<<dist<<G4endl ;
540
541 insideP = orbMirror->Inside(orbP);
542 G4cout<<"orbMirror->Inside(orbP) = "<<OutputInside(insideP)<<G4endl ;
543 dist = orbMirror->DistanceToIn(orbP,orbV);
544 G4cout<<"orbMirror->DistanceToIn(orbP,orbV) = "<<dist<<G4endl ;
545
546 insideP = ellShell->Inside(orbP);
547 G4cout<<"ellShell->Inside(orbP) = "<<OutputInside(insideP)<<G4endl ;
548 dist = ellShell->DistanceToIn(orbP,orbV);
549 G4cout<<"ellShell->DistanceToIn(orbP,orbV) = "<<dist<<G4endl ;
550
551 insideP = ellMirror->Inside(orbP);
552 G4cout<<"ellMirror->Inside(orbP) = "<<OutputInside(insideP)<<G4endl ;
553 dist = ellMirror->DistanceToIn(orbP,orbV);
554 G4cout<<"ellMirror->DistanceToIn(orbP,orbV) = "<<dist<<G4endl ;
555
556 for(i=0;i<3000;i++)
557 {
558 orbTheta = 45*degree*G4UniformRand();
559 orbPhi = 360*degree*G4UniformRand();
560 orbV = G4ThreeVector(std::sin(orbTheta)*std::cos(orbPhi),
561 std::sin(orbTheta)*std::sin(orbPhi),
562 std::cos(orbTheta));
563 dist = ellMirror->DistanceToIn(orbP,orbV);
564
565 if(dist != kInfinity)
566 {
567 G4cout<<i<<"\t"<<"ellMirror->DistanceToIn(orbP,orbV) = "<<dist<<G4endl ;
568 }
569 }
570 G4cout<<G4endl;
571
572 G4cout<<"Tracking functions are OK"<<G4endl ;
573
574
575// Point on surface
576
577 G4cout<<G4endl;
578 // G4cout<<"Point on surface of t1Sb3"<<G4endl<<G4endl;
579 G4cout<<"Point on surface of b5St1"<<G4endl<<G4endl;
580 for(i=0;i<30;i++)
581 {
582 // p = t1Sb3.GetPointOnSurface();
583 p = b5St1.GetPointOnSurface();
584 G4cout<<p.x()<<"\t\t"<<p.y()<<"\t\t"<<p.z()<<G4endl;
585 }
586 G4cout<<G4endl;
587
588
589
590// CalculateExtent
591
592 G4VoxelLimits limit; // Unlimited
593 G4RotationMatrix noRot;
594 G4AffineTransform origin;
595 G4double min,max;
596
597 assert(b1.CalculateExtent(kXAxis,limit,origin,min,max));
598 assert(ApproxEqual(min,-20)&&ApproxEqual(max,20));
599
600 assert(b1.CalculateExtent(kYAxis,limit,origin,min,max));
601 assert(ApproxEqual(min,-30)&&ApproxEqual(max,30));
602
603 assert(b1.CalculateExtent(kZAxis,limit,origin,min,max));
604 assert(ApproxEqual(min,-40)&&ApproxEqual(max,40));
605
606 assert(b1Sb4.CalculateExtent(kXAxis,limit,origin,min,max));
607 assert(ApproxEqual(min,-20)&&ApproxEqual(max,20));
608
609 assert(b1Sb4.CalculateExtent(kYAxis,limit,origin,min,max));
610 // G4cout<<"min of b1Sb4.CalculateExtent(kYAxis,limit,origin,min,max) = "
611 // <<min<<G4endl ;
612 // G4cout<<"max of b1Sb4.CalculateExtent(kYAxis,limit,origin,min,max) = "
613 // <<max<<G4endl ;
614 assert(ApproxEqual(min,-30)&&ApproxEqual(max,30));
615
616 assert(b1Sb4.CalculateExtent(kZAxis,limit,origin,min,max));
617 assert(ApproxEqual(min,-40)&&ApproxEqual(max,40));
618
619
620
621 G4ThreeVector pmxmymz(-100,-110,-120);
622 G4AffineTransform tPosOnly(pmxmymz);
623
624 assert(b1.CalculateExtent(kXAxis,limit,tPosOnly,min,max));
625 assert(ApproxEqual(min,-120)&&ApproxEqual(max,-80));
626
627 assert(b1.CalculateExtent(kYAxis,limit,tPosOnly,min,max));
628 assert(ApproxEqual(min,-140)&&ApproxEqual(max,-80));
629
630 assert(b1.CalculateExtent(kZAxis,limit,tPosOnly,min,max));
631 assert(ApproxEqual(min,-160)&&ApproxEqual(max,-80));
632
633 G4RotationMatrix r90Z;
634 r90Z.rotateZ(pi/2);
635 G4AffineTransform tRotZ(r90Z,pzero);
636
637 assert(b1.CalculateExtent(kXAxis,limit,tRotZ,min,max));
638 assert(ApproxEqual(min,-30)&&ApproxEqual(max,30));
639
640 assert(b1.CalculateExtent(kYAxis,limit,tRotZ,min,max));
641 assert(ApproxEqual(min,-20)&&ApproxEqual(max,20));
642
643 assert(b1.CalculateExtent(kZAxis,limit,tRotZ,min,max));
644 assert(ApproxEqual(min,-40)&&ApproxEqual(max,40));
645
646// Check that clipped away
647
648 G4VoxelLimits xClip;
649 xClip.AddLimit(kXAxis,-100,-50);
650
651 assert(!b1.CalculateExtent(kXAxis,xClip,origin,min,max));
652
653// Assert clipped to volume
654
655 G4VoxelLimits allClip;
656 allClip.AddLimit(kXAxis,-5,+5);
657 allClip.AddLimit(kYAxis,-5,+5);
658 allClip.AddLimit(kZAxis,-5,+5);
659
660 G4RotationMatrix genRot;
661 genRot.rotateX(pi/6);
662 genRot.rotateY(pi/6);
663 genRot.rotateZ(pi/6);
664 G4AffineTransform tGen(genRot,vx);
665
666 assert(b1.CalculateExtent(kXAxis,allClip,tGen,min,max));
667 assert(ApproxEqual(min,-5)&&ApproxEqual(max,5));
668
669 assert(b1.CalculateExtent(kYAxis,allClip,tGen,min,max));
670 assert(ApproxEqual(min,-5)&&ApproxEqual(max,5));
671
672 assert(b1.CalculateExtent(kZAxis,allClip,tGen,min,max));
673 assert(ApproxEqual(min,-5)&&ApproxEqual(max,5));
674
675 G4VoxelLimits buggyClip2;
676 buggyClip2.AddLimit(kXAxis,5,15);
677
678 assert(b1.CalculateExtent(kXAxis,buggyClip2,origin,min,max));
679 assert(ApproxEqual(min,5)&&ApproxEqual(max,15));
680
681 assert(b1.CalculateExtent(kYAxis,buggyClip2,origin,min,max));
682 assert(ApproxEqual(min,-30)&&ApproxEqual(max,30));
683
684 assert(b1.CalculateExtent(kZAxis,buggyClip2,origin,min,max));
685 assert(ApproxEqual(min,-40)&&ApproxEqual(max,40));
686
687 buggyClip2.AddLimit(kYAxis,5,15);
688
689 assert(b1.CalculateExtent(kXAxis,buggyClip2,origin,min,max));
690 assert(ApproxEqual(min,5)&&ApproxEqual(max,15));
691
692 assert(b1.CalculateExtent(kYAxis,buggyClip2,origin,min,max));
693 assert(ApproxEqual(min,5)&&ApproxEqual(max,15));
694
695 assert(b1.CalculateExtent(kZAxis,buggyClip2,origin,min,max));
696 assert(ApproxEqual(min,-40)&&ApproxEqual(max,40));
697
698 G4VoxelLimits buggyClip1;
699 buggyClip1.AddLimit(kXAxis,-5,+5);
700
701 assert(b1.CalculateExtent(kXAxis,buggyClip1,origin,min,max));
702 assert(ApproxEqual(min,-5)&&ApproxEqual(max,5));
703
704 assert(b1.CalculateExtent(kYAxis,buggyClip1,origin,min,max));
705 assert(ApproxEqual(min,-30)&&ApproxEqual(max,30));
706
707 assert(b1.CalculateExtent(kZAxis,buggyClip1,origin,min,max));
708 assert(ApproxEqual(min,-40)&&ApproxEqual(max,40));
709
710 buggyClip1.AddLimit(kYAxis,-5,+5);
711
712 assert(b1.CalculateExtent(kXAxis,buggyClip1,origin,min,max));
713 assert(ApproxEqual(min,-5)&&ApproxEqual(max,5));
714
715 assert(b1.CalculateExtent(kYAxis,buggyClip1,origin,min,max));
716 assert(ApproxEqual(min,-5)&&ApproxEqual(max,5));
717
718 assert(b1.CalculateExtent(kZAxis,buggyClip1,origin,min,max));
719 assert(ApproxEqual(min,-40)&&ApproxEqual(max,40));
720
721 G4cout<<"CalculateExtent is OK "<<G4endl ;
722
723 delete rotDz5;
724 G4int out =0 ;
725 return out ;
726}
727
728
729
730
731
732
733
734
Note: See TracBrowser for help on using the repository browser.