source: trunk/source/geometry/solids/Boolean/test/testG4UnionSolid.cc@ 1330

Last change on this file since 1330 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: 26.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
29#include <assert.h>
30#include <cmath>
31
32#include "globals.hh"
33#include "geomdefs.hh"
34
35#include "ApproxEqual.hh"
36
37#include "G4ThreeVector.hh"
38#include "G4RotationMatrix.hh"
39#include "G4AffineTransform.hh"
40#include "G4VoxelLimits.hh"
41
42#include "G4Box.hh"
43#include "G4Cons.hh"
44#include "G4Para.hh"
45#include "G4Sphere.hh"
46#include "G4Torus.hh"
47#include "G4Trap.hh"
48#include "G4Trd.hh"
49#include "G4Tubs.hh"
50
51#include "G4UnionSolid.hh"
52// #include "G4DisplacedSolid.hh"
53
54///////////////////////////////////////////////////////////////////
55//
56// Dave's auxiliary function
57
58const G4String OutputInside(const EInside a)
59{
60 switch(a)
61 {
62 case kInside: return "Inside";
63 case kOutside: return "Outside";
64 case kSurface: return "Surface";
65 }
66 return "????";
67}
68
69
70int main()
71{
72 G4ThreeVector pzero(0,0,0);
73 G4ThreeVector ponxside(20,0,0),ponyside(0,30,0),ponzside(0,0,40),
74 ponb2x(10,0,0),ponb2y(0,10,0),ponb2z(0,0,10),
75 ponb2mx(-10,0,0),ponb2my(0,-10,0),ponb2mz(0,0,-10);
76 G4ThreeVector ponmxside(-20,0,0),ponmyside(0,-30,0),ponmzside(0,0,-40);
77 G4ThreeVector ponzsidey(0,25,40),ponmzsidey(0,25,-40),
78 ponb2zy(0,5,10),ponb2mzy(0,5,-10) ;
79
80 G4ThreeVector pbigx(100,0,0),pbigy(0,100,0),pbigz(0,0,100);
81 G4ThreeVector pbigmx(-100,0,0),pbigmy(0,-100,0),pbigmz(0,0,-100);
82
83 G4ThreeVector vx(1,0,0),vy(0,1,0),vz(0,0,1);
84 G4ThreeVector vmx(-1,0,0),vmy(0,-1,0),vmz(0,0,-1);
85 G4ThreeVector vxy(1/std::sqrt(2.0),1/std::sqrt(2.0),0);
86 G4ThreeVector vmxy(-1/std::sqrt(2.0),1/std::sqrt(2.0),0);
87 G4ThreeVector vmxmy(-1/std::sqrt(2.0),-1/std::sqrt(2.0),0);
88 G4ThreeVector vxmy(1/std::sqrt(2.0),-1/std::sqrt(2.0),0);
89
90 G4ThreeVector vxz(1/std::sqrt(2.0),0,1/std::sqrt(2.0));
91 G4ThreeVector vxmz(1/std::sqrt(2.0),0,-1/std::sqrt(2.0));
92 G4ThreeVector vmxz(-1/std::sqrt(2.0),0,1/std::sqrt(2.0));
93 G4ThreeVector vmxmz(-1/std::sqrt(2.0),0,-1/std::sqrt(2.0));
94
95 G4double dist, vol, volCheck, diff;
96 G4ThreeVector *pNorm,norm;
97 G4bool *pgoodNorm,goodNorm,calcNorm=true;
98
99 pNorm=&norm;
100 pgoodNorm=&goodNorm;
101
102 G4RotationMatrix identity, xRot ;
103
104// NOTE: xRot = rotation such that x axis->-x axis & y axis->-y axis
105
106 xRot.rotateZ(-pi) ;
107
108 G4Transform3D transform(xRot,G4ThreeVector(0,40,0)) ;
109
110 G4Box b1("Test Box #1",20,30,40);
111 G4Box b2("Test Box #2",10,10,10);
112 G4Box b3("Test Box #3",10,20,50);
113 G4Box* b4 = new G4Box("b4",50,50,50) ;
114 G4Box* b5 = new G4Box("b5",10,10,60) ;
115
116 G4Tubs t1("Solid Tube #1",0,50,50,0,360);
117 G4Tubs t2("Hole Tube #2",45,50,50,0,360);
118 G4Tubs t3("disk",0,50,5,0,360);
119 G4Tubs t4("t4",0,45,50,0,360);
120
121 G4Cons c1("Hollow Full Tube",50,100,50,100,50,0,2*pi),
122 c2("Full Cone",0,50,0,100,50,0,2*pi) ;
123
124 G4UnionSolid b1Ub2("b1Unionb2",&b1,&b2),
125 t1Ub2("t1Unionb2",&t1,&b2),
126 c2Ub2("c2Unionb2",&c2,&b2);
127
128 G4UnionSolid b4Ub5("b4Ub5",b4,b5);
129
130 G4UnionSolid t2Ut4("t2Ut4",&t2,&t4);
131
132 // With placement
133
134 G4UnionSolid b1Ub3("b1Unionb3",&b1,&b3,transform);
135
136 G4UnionSolid t1Ub3("t1Unionb3",&t1,&b3,transform);
137
138 G4UnionSolid t2Ut1Trans("t2Ut1Trans",&t2,&t1,
139 &identity,G4ThreeVector(0,0,110));
140
141 G4UnionSolid t2Ut1Ut1("t2Ut1Ut1",&t2Ut1Trans,&t1,
142 &identity,G4ThreeVector(0,0,-110));
143
144 G4UnionSolid t2Ut3("t2Ut3",&t2,&t3,
145 &identity,G4ThreeVector(0,0,50)) ;
146
147 G4UnionSolid t2Ut3Ut3("t2Ut3Ut3",&t2Ut3,&t3,
148 &identity,G4ThreeVector(0,0,-50)) ;
149
150 G4UnionSolid t3Ut3("t3Ut3",&t3,&t3,
151 &identity,G4ThreeVector(0,0,10)) ;
152
153 G4UnionSolid t3Ut3Ut3("t3Ut3Ut3",&t3Ut3,&t3,
154 &identity,G4ThreeVector(0,0,-10)) ;
155
156
157 G4UnionSolid b4Ub4xtouch("b4Ub4xtouch",b4,b4,&identity,G4ThreeVector(100.,0,0));
158
159
160 G4Box * box1 = new G4Box("Box1",1092.500000,240.103374,92.000000);
161 G4Box * box2 = new G4Box("Box2",540.103374,792.500000,92.000000);
162
163 G4double L1 = 1104;
164
165 G4UnionSolid envelope("ECShapeBoxes",
166 box1,
167 box2,
168 0,
169 G4ThreeVector(-L1/2.,
170 L1/2.,
171 0.) );
172
173 G4ThreeVector pPrePre(-301.4252442474112,903.7985455093324,3013.122572566068-2924);
174 G4ThreeVector pPre(-301.5459060972926,904.1596390942879,3014.325585765819-2924);
175 G4ThreeVector pPos(-301.7139066520877,904.6621760944698,3016-2924);
176
177 // Vacuum Cross
178
179 G4double startPhi = 0.*deg;
180 G4double deltaPhi = 360.*deg;
181 G4double minRadius = 0.*cm;
182
183
184 //Sample changer (carbon-fiber cross)
185
186 //G4double minRadiusTUV = 3.9*cm; // Vertical carbon fiber tube
187 //G4double maxRadiusTUV = 4.1*cm;
188 //G4double halfLengthTUV = 47.5*cm;
189 G4double maxRadiusVACTUV = 3.9*cm;
190 G4double halfLengthVACTUV = 47.5*cm;
191
192 //G4double minRadiusTUH = 2.5*cm; // ??? Horizontal carbon fiber tube
193 //G4double maxRadiusTUH = 2.7*cm;
194 //G4double halfLengthTUH = 44.0*cm;
195 G4double maxRadiusVACTUH = 2.5*cm;
196 G4double halfLengthVACTUH = 44.0*cm;
197
198
199 G4RotationMatrix rmVacCross;
200 rmVacCross.rotateX(90.*deg);
201 G4RotationMatrix rmVACTUV;
202 rmVACTUV.rotateX(90.*deg);
203
204 G4Tubs* solidVACTUH = new G4Tubs("VACTUH",minRadius,maxRadiusVACTUH,
205 halfLengthVACTUH,startPhi,deltaPhi);
206
207 G4Tubs* solidVACTUV = new G4Tubs("VACTUV",minRadius,maxRadiusVACTUV,
208 halfLengthVACTUV,startPhi,deltaPhi);
209
210 G4UnionSolid* solidVacCross = new G4UnionSolid("VacCross",solidVACTUH,solidVACTUV,
211 &rmVacCross,G4ThreeVector());
212
213
214 // check cubic volume
215
216 vol = b1Ub2.GetCubicVolume();
217 volCheck = 8.*20.*30.*40.;
218 // G4cout<<"vol = "<<vol<<"; volCheck = "<<volCheck<<G4endl;
219 assert(ApproxEqual(vol,volCheck));
220
221
222 // t2Ut4.SetCubVolStatistics(10000000);
223 vol = t2Ut4.GetCubicVolume();
224 volCheck = 2.*pi*50.*50.*50.;
225 diff = std::abs(vol-volCheck)/volCheck;
226 G4cout<<"vol = "<<vol<<"; volCheck = "<<volCheck<<"; rel. diff = "<<diff<<G4endl;
227 // assert(ApproxEqual(vol,volCheck));
228
229
230
231// Check Inside
232 EInside side = solidVacCross->Inside(G4ThreeVector( 2.3391096733692156,
233 1.1325145357709736,
234 -0.85000000000000009));
235 G4cout<<"solidVacCross->Inside(G4ThreeVector( 2.3391... = "
236 <<OutputInside(side)<<G4endl;
237
238 side = envelope.Inside(pPrePre);
239 G4cout<<"envelope.Inside(pPrePre) = "<<OutputInside(side)<<G4endl;
240
241 side = envelope.Inside(pPre);
242 G4cout<<"envelope.Inside(pPre) = "<<OutputInside(side)<<G4endl;
243
244 side = envelope.Inside(pPos);
245 G4cout<<"envelope.Inside(pPos) = "<<OutputInside(side)<<G4endl;
246
247 side = b4Ub4xtouch.Inside(G4ThreeVector(50.,0,0));
248 G4cout<<"b4Ub4xtouch.Inside(G4ThreeVector(50.,0,0)) = "<<OutputInside(side)<<G4endl;
249
250 side = b4Ub4xtouch.Inside(G4ThreeVector(50.,50.,0));
251 G4cout<<"b4Ub4xtouch.Inside(G4ThreeVector(50.,50.,0)) = "<<OutputInside(side)<<G4endl;
252
253
254
255
256
257 assert(b1Ub2.Inside(pzero)==kInside);
258 assert(b1Ub2.Inside(pbigz)==kOutside);
259 assert(b1Ub2.Inside(ponxside)==kSurface);
260 assert(b1Ub2.Inside(ponyside)==kSurface);
261 assert(b1Ub2.Inside(ponzside)==kSurface);
262
263 assert(t1Ub2.Inside(pzero)==kInside);
264 assert(t1Ub2.Inside(pbigz)==kOutside);
265 assert(t1Ub2.Inside(ponb2x)==kInside);
266 assert(t1Ub2.Inside(ponb2y)==kInside);
267 assert(t1Ub2.Inside(ponb2z)==kInside);
268
269 assert(c2Ub2.Inside(pzero)==kInside);
270 assert(c2Ub2.Inside(pbigz)==kOutside);
271 assert(c2Ub2.Inside(ponb2x)==kInside);
272 assert(c2Ub2.Inside(ponb2y)==kInside);
273 assert(c2Ub2.Inside(ponb2z)==kInside);
274
275 // With placement
276
277 assert(t1Ub3.Inside(pzero)==kInside);
278 assert(t1Ub3.Inside(G4ThreeVector(0,60,0))==kSurface);
279 assert(t1Ub3.Inside(G4ThreeVector(50,0,0))==kSurface);
280 assert(t1Ub3.Inside(G4ThreeVector(0,65,0))==kOutside);
281 assert(t1Ub3.Inside(G4ThreeVector(0,50,0))==kInside);
282
283
284// Check Surface Normal
285
286 G4ThreeVector normal;
287
288 normal=b1Ub2.SurfaceNormal(ponxside);
289 assert(ApproxEqual(normal,G4ThreeVector(1,0,0)));
290
291 normal=b1Ub2.SurfaceNormal(ponmxside);
292 assert(ApproxEqual(normal,G4ThreeVector(-1,0,0)));
293
294 normal=b1Ub2.SurfaceNormal(ponyside);
295 assert(ApproxEqual(normal,G4ThreeVector(0,1,0)));
296
297 normal=b1Ub2.SurfaceNormal(ponmyside);
298 assert(ApproxEqual(normal,G4ThreeVector(0,-1,0)));
299
300 normal=b1Ub2.SurfaceNormal(ponzside);
301 assert(ApproxEqual(normal,G4ThreeVector(0,0,1)));
302
303 normal=b1Ub2.SurfaceNormal(ponmzside);
304 assert(ApproxEqual(normal,G4ThreeVector(0,0,-1)));
305
306 // Wirh placemenet
307
308 normal=t1Ub3.SurfaceNormal(G4ThreeVector(0,60,0));
309 assert(ApproxEqual(normal,G4ThreeVector(0,1,0)));
310
311 normal=t1Ub3.SurfaceNormal(G4ThreeVector(10,55,0));
312 assert(ApproxEqual(normal,G4ThreeVector(1,0,0)));
313
314 normal=t1Ub3.SurfaceNormal(G4ThreeVector(-10,55,0));
315 assert(ApproxEqual(normal,G4ThreeVector(-1,0,0)));
316
317 normal=t1Ub3.SurfaceNormal(G4ThreeVector(50.0/std::sqrt(2.0),50.0/std::sqrt(2.0),0));
318 assert(ApproxEqual(normal,vxy));
319
320
321// DistanceToOut(P)
322
323 dist=b1Ub2.DistanceToOut(pzero);
324 assert(ApproxEqual(dist,20));
325
326 dist=b1Ub2.DistanceToOut(vx);
327 assert(ApproxEqual(dist,19));
328
329 // With placement
330
331 dist=t1Ub3.DistanceToOut(pzero);
332 assert(ApproxEqual(dist,50));
333
334 dist=t1Ub3.DistanceToOut(G4ThreeVector(0,55,0));
335 assert(ApproxEqual(dist,5));
336
337 dist=t1Ub3.DistanceToOut(G4ThreeVector(45,0,0));
338 assert(ApproxEqual(dist,5));
339
340 dist=t1Ub3.DistanceToOut(G4ThreeVector(0,45,0));
341 assert(ApproxEqual(dist,10));
342
343 dist=t1Ub3.DistanceToOut(G4ThreeVector(0,-45,0));
344 assert(ApproxEqual(dist,5));
345
346
347// DistanceToOut(P,V)
348
349 dist=b4Ub5.DistanceToOut(G4ThreeVector(50,0,0),vmx,
350 calcNorm,pgoodNorm,pNorm);
351 assert(ApproxEqual(dist,100));
352 // G4cout<<"(100) b4Ub5.DistanceToOut(G4ThreeVector(50,0,0),vmx) = "
353 // <<dist<<G4endl ;
354
355 dist=b4Ub5.DistanceToOut(G4ThreeVector(-10,0,-60),vxz,
356 calcNorm,pgoodNorm,pNorm);
357 assert(ApproxEqual(dist,84.8528));
358 // G4cout<<"(84.85) b4Ub5.DistanceToOut(G4ThreeVector(-10,0,-60),vxz) = "
359 // <<dist<<G4endl ;
360
361
362 dist=b1Ub2.DistanceToOut(pzero,vx,calcNorm,pgoodNorm,pNorm);
363 assert(ApproxEqual(dist,20)&&ApproxEqual(*pNorm,vx)); // &&*pgoodNorm);
364
365 dist=b1Ub2.DistanceToOut(pzero,vmx,calcNorm,pgoodNorm,pNorm);
366 assert(ApproxEqual(dist,20)&&ApproxEqual(norm,vmx)); // &&*pgoodNorm);
367
368 dist=b1Ub2.DistanceToOut(pzero,vy,calcNorm,pgoodNorm,pNorm);
369 assert(ApproxEqual(dist,30)&&ApproxEqual(norm,vy)); // &&*pgoodNorm);
370
371 dist=b1Ub2.DistanceToOut(pzero,vmy,calcNorm,pgoodNorm,pNorm);
372 assert(ApproxEqual(dist,30)&&ApproxEqual(norm,vmy)); // &&*pgoodNorm);
373
374 dist=b1Ub2.DistanceToOut(pzero,vz,calcNorm,pgoodNorm,pNorm);
375 assert(ApproxEqual(dist,40)&&ApproxEqual(norm,vz)); // &&*pgoodNorm);
376
377 dist=b1Ub2.DistanceToOut(pzero,vmz,calcNorm,pgoodNorm,pNorm);
378 assert(ApproxEqual(dist,40)&&ApproxEqual(norm,vmz)); // &&*pgoodNorm);
379
380 dist=b1Ub2.DistanceToOut(pzero,vxy,calcNorm,pgoodNorm,pNorm);
381 assert(ApproxEqual(dist,std::sqrt(800.))); // &&*pgoodNorm);
382
383 dist=b1Ub2.DistanceToOut(ponb2x,vx,calcNorm,pgoodNorm,pNorm);
384 assert(ApproxEqual(dist,10)&&ApproxEqual(*pNorm,vx)); // &&*pgoodNorm);
385
386 dist=b1Ub2.DistanceToOut(ponb2x,vmx,calcNorm,pgoodNorm,pNorm);
387 assert(ApproxEqual(dist,30)&&ApproxEqual(*pNorm,vmx)); // &&*pgoodNorm);
388
389 dist=b1.DistanceToOut(ponxside,vy,calcNorm,pgoodNorm,pNorm);
390// G4cout<<"b1.DistanceToOut(ponxside,vy) = "<<dist<<G4endl;
391// assert(ApproxEqual(dist,0)&&ApproxEqual(*pNorm,vy)); // &&*pgoodNorm);
392
393 dist=b1Ub2.DistanceToOut(ponb2mx,vmx,calcNorm,pgoodNorm,pNorm);
394 assert(ApproxEqual(dist,10)&&ApproxEqual(norm,vmx)); //&&*pgoodNorm);
395
396 dist=b1Ub2.DistanceToOut(ponb2y,vy,calcNorm,pgoodNorm,pNorm);
397 assert(ApproxEqual(dist,20)&&ApproxEqual(norm,vy)); // &&*pgoodNorm);
398
399 dist=b1Ub2.DistanceToOut(ponb2my,vmy,calcNorm,pgoodNorm,pNorm);
400 assert(ApproxEqual(dist,20)&&ApproxEqual(norm,vmy)) ; // &&*pgoodNorm);
401
402 dist=b1Ub2.DistanceToOut(ponb2z,vz,calcNorm,pgoodNorm,pNorm);
403 assert(ApproxEqual(dist,30)&&ApproxEqual(norm,vz)) ; // &&*pgoodNorm);
404
405 dist=b1Ub2.DistanceToOut(ponb2mz,vmz,calcNorm,pgoodNorm,pNorm);
406 assert(ApproxEqual(dist,30)&&ApproxEqual(norm,vmz)) ; // &&*pgoodNorm);
407
408 // With placement
409
410 dist=t1Ub3.DistanceToOut(G4ThreeVector(0,55,0),vmy,
411 calcNorm,pgoodNorm,pNorm);
412 assert(ApproxEqual(dist,105)&&ApproxEqual(norm,vmy)); // &&*pgoodNorm);
413
414 dist=t1Ub3.DistanceToOut(G4ThreeVector(0,55,0),vy,
415 calcNorm,pgoodNorm,pNorm);
416 assert(ApproxEqual(dist,5)&&ApproxEqual(norm,vy)); // &&*pgoodNorm);
417
418 dist=t1Ub3.DistanceToOut(G4ThreeVector(0,55,0),vx,
419 calcNorm,pgoodNorm,pNorm);
420 assert(ApproxEqual(dist,10)&&ApproxEqual(norm,vx)); // &&*pgoodNorm);
421
422 dist=t1Ub3.DistanceToOut(G4ThreeVector(0,55,0),vmx,
423 calcNorm,pgoodNorm,pNorm);
424 assert(ApproxEqual(dist,10)&&ApproxEqual(norm,vmx)); // &&*pgoodNorm);
425
426 dist=t1Ub3.DistanceToOut(G4ThreeVector(0,55,0),vz,
427 calcNorm,pgoodNorm,pNorm);
428 assert(ApproxEqual(dist,50)&&ApproxEqual(norm,vz)); // &&*pgoodNorm);
429
430 dist=t2Ut1Ut1.DistanceToOut(G4ThreeVector(50,0,0),vmx,
431 calcNorm,pgoodNorm,pNorm);
432 assert(ApproxEqual(dist,5));
433 // G4cout<<"(5) t2Ut1Ut1.DistanceToOut(G4ThreeVector(50,0,0),vmx) = "
434 // <<dist<<G4endl ;
435
436 dist=t2Ut1Ut1.DistanceToOut(G4ThreeVector(50,0,0),vx,
437 calcNorm,pgoodNorm,pNorm);
438 assert(ApproxEqual(dist,0));
439 // G4cout<<"(0) t2Ut1Ut1.DistanceToOut(G4ThreeVector(50,0,0),vx) = "
440 // <<dist<<G4endl ;
441
442 dist=t2Ut1Ut1.DistanceToOut(G4ThreeVector(50,0,70),vmx,
443 calcNorm,pgoodNorm,pNorm);
444 assert(ApproxEqual(dist,100));
445 // G4cout<<"(100) t2Ut1Ut1.DistanceToOut(G4ThreeVector(50,0,70),vmx) = "
446 // <<dist<<G4endl ;
447
448 dist=t2Ut1Ut1.DistanceToOut(G4ThreeVector(50,0,70),vx,
449 calcNorm,pgoodNorm,pNorm);
450 assert(ApproxEqual(dist,0));
451 // G4cout<<"(0) t2Ut1Ut1.DistanceToOut(G4ThreeVector(50,0,70),vx) = "
452 // <<dist<<G4endl ;
453
454 dist=t2Ut1Ut1.DistanceToOut(G4ThreeVector(0,0,60),vmz,
455 calcNorm,pgoodNorm,pNorm);
456 assert(ApproxEqual(dist,0));
457 // G4cout<<"(0) t2Ut1Ut1.DistanceToOut(G4ThreeVector(0,0,60),vmz) = "
458 // <<dist<<G4endl ;
459
460 dist=t2Ut1Ut1.DistanceToOut(G4ThreeVector(0,0,60),vz,
461 calcNorm,pgoodNorm,pNorm);
462 assert(ApproxEqual(dist,100));
463 // G4cout<<"(100) t2Ut1Ut1.DistanceToOut(G4ThreeVector(0,0,60),vz) = "
464 // <<dist<<G4endl ;
465
466 dist=t2Ut3Ut3.DistanceToOut(G4ThreeVector(-50,0,-48),vxmz,
467 calcNorm,pgoodNorm,pNorm);
468 assert(ApproxEqual(dist,9.89949));
469 // G4cout<<"(9.89) t2Ut3Ut3.DistanceToOut(G4ThreeVector(-50,0,-48),vxmz) = "
470 // <<dist<<G4endl ;
471
472 dist=t2Ut3Ut3.DistanceToOut(G4ThreeVector(-50,0,-50),vxz,
473 calcNorm,pgoodNorm,pNorm);
474 assert(ApproxEqual(dist,7.07107));
475 // G4cout<<"(7.07) t2Ut3Ut3.DistanceToOut(G4ThreeVector(-50,0,-50),vxz) = "
476 // <<dist<<G4endl ;
477
478 dist=t2Ut3Ut3.DistanceToOut(G4ThreeVector(-50,0,-50),vxmz,
479 calcNorm,pgoodNorm,pNorm);
480 assert(ApproxEqual(dist,7.07107));
481 // G4cout<<"(7.07) t2Ut3Ut3.DistanceToOut(G4ThreeVector(-50,0,-50),vxmz) = "
482 // <<dist<<G4endl ;
483
484 dist=t2Ut3Ut3.DistanceToOut(G4ThreeVector(-50,0,-52),vxmz,
485 calcNorm,pgoodNorm,pNorm);
486 assert(ApproxEqual(dist,4.24264));
487 // G4cout<<"(4.24) t2Ut3Ut3.DistanceToOut(G4ThreeVector(-50,0,-52),vxmz) = "
488 // <<dist<<G4endl ;
489
490 dist=t3Ut3Ut3.DistanceToOut(G4ThreeVector(0,0,-15),vz,
491 calcNorm,pgoodNorm,pNorm);
492 assert(ApproxEqual(dist,30));
493 // G4cout<<"(30) t3Ut3Ut3.DistanceToOut(G4ThreeVector(0,0,-15),vz) = "
494 // <<dist<<G4endl ;
495
496 dist=t2Ut4.DistanceToOut(G4ThreeVector(-50,0,0),vx,
497 calcNorm,pgoodNorm,pNorm);
498 assert(ApproxEqual(dist,100));
499 // G4cout<<"(100) t2Ut4.DistanceToOut(G4ThreeVector(-50,0,0),vx) = "
500 // <<dist<<G4endl ;
501
502 dist=t2Ut4.DistanceToOut(G4ThreeVector(-15,-45,0),vx,
503 calcNorm,pgoodNorm,pNorm);
504 assert(ApproxEqual(dist,36.7945));
505 // G4cout<<"(36.8) t2Ut4.DistanceToOut(G4ThreeVector(-15,-45,0),vx) = "
506 // <<dist<<G4endl ;
507
508 dist=t2Ut4.DistanceToOut(G4ThreeVector(0,-45,0),vx,
509 calcNorm,pgoodNorm,pNorm);
510 assert(ApproxEqual(dist,21.7945));
511
512 // G4cout<<"(21.8) t2Ut4.DistanceToOut(G4ThreeVector(0,-45,0),vx) = "
513 // <<dist<<G4endl ;
514
515 const G4double dx=15.3125;
516 const G4double dy=100;
517 const G4double dz=3.75;
518 const G4double theta=0.283794;
519
520 const G4double radius=2;
521 const G4double length=200; // cm
522
523 // Building solid
524
525 G4Para* para=new G4Para("Para",dx*cm,dy*cm,dz*cm,0,theta,0);
526 G4Tubs* tube = new G4Tubs("Tube",0.0,radius*cm*0.99, length*cm/2.0, 0.0, 2*M_PI);
527
528 G4ThreeVector transv(-164.062,0,-17.5); // mm
529 G4Transform3D trans1=G4Translate3D(transv*mm)*G4RotateX3D(M_PI/2);
530
531 G4VSolid* solid=new G4UnionSolid("solid",para,tube,trans1);
532
533 // Checking points
534 G4ThreeVector point1(-142.188*mm,0,-0.5*mm); // mm
535 G4ThreeVector point2(-142.188*mm,0,-1*mm); // mm
536 G4ThreeVector direction(-1,0,0);
537
538 if (solid->Inside(point1)==kInside) {
539 // G4cout<<"Test point "<<point1<<" is inside. That's right!"<<G4endl;
540 }
541 else {
542 // G4cout<<"Test point "<<point1<<" is not inside. That's wrong!"<<G4endl;
543 }
544
545 dist=solid->DistanceToOut(point1,direction);
546
547 // G4cout<<"Distance is "<<dist<<G4endl;
548
549 if (solid->Inside(point2)==kInside) {
550 // G4cout<<"Test point "<<point2<<" is inside. That's right!"<<G4endl;
551 }
552 else {
553 // G4cout<<"Test point "<<point2<<" is not inside. That's wrong!"<<G4endl;
554 }
555
556 // dist=solid->DistanceToOut(point2,direction);
557
558 // G4cout<<"Distance is "<<dist<<G4endl;
559
560
561
562
563// DistanceToIn(P)
564
565 dist=b1Ub2.DistanceToIn(pbigx);
566 assert(ApproxEqual(dist,80));
567 // G4cout<<"b1Ub2.DistanceToIn(pbigx) = "<<dist<<G4endl ;
568 // G4cout<<"b1.DistanceToIn(pbigx) = "<<b1.DistanceToIn(pbigx)<<G4endl ;
569 // G4cout<<"b2.DistanceToIn(pbigx) = "<<b2.DistanceToIn(pbigx)<<G4endl ;
570
571 dist=b1Ub2.DistanceToIn(ponxside);
572 assert(ApproxEqual(dist,0));
573
574 dist=b1Ub2.DistanceToIn(ponxside);
575 assert(ApproxEqual(dist,0));
576
577 dist=b1Ub2.DistanceToIn(pbigmy);
578 assert(ApproxEqual(dist,70));
579
580 dist=b1Ub2.DistanceToIn(pbigz);
581 assert(ApproxEqual(dist,60));
582
583 dist=b1Ub2.DistanceToIn(pbigmz);
584 assert(ApproxEqual(dist,60));
585
586 // With placement
587
588 dist=t1Ub3.DistanceToIn(pbigmz);
589 assert(ApproxEqual(dist,50));
590
591 dist=t1Ub3.DistanceToIn(pbigy);
592 assert(ApproxEqual(dist,40));
593
594 dist=t1Ub3.DistanceToIn(G4ThreeVector(0,65,0));
595 assert(ApproxEqual(dist,5));
596
597 dist=t1Ub3.DistanceToIn(G4ThreeVector(0,60,0));
598 assert(ApproxEqual(dist,0));
599
600
601// DistanceToIn(P,V)
602
603 dist=b1Ub2.DistanceToIn(pbigx,vmx);
604 assert(ApproxEqual(dist,80));
605
606 dist=b1Ub2.DistanceToIn(pbigmx,vx);
607 assert(ApproxEqual(dist,80));
608
609 dist=b1Ub2.DistanceToIn(pbigy,vmy);
610 assert(ApproxEqual(dist,70));
611
612 dist=b1Ub2.DistanceToIn(pbigmy,vy);
613 assert(ApproxEqual(dist,70));
614
615 dist=b1Ub2.DistanceToIn(pbigz,vmz);
616 assert(ApproxEqual(dist,60));
617
618 dist=b1Ub2.DistanceToIn(pbigmz,vz);
619 assert(ApproxEqual(dist,60));
620
621 dist=b1Ub2.DistanceToIn(pbigx,vxy);
622 assert(ApproxEqual(dist,kInfinity));
623
624 dist=b1Ub2.DistanceToIn(pbigmx,vxy);
625 assert(ApproxEqual(dist,kInfinity));
626
627
628 // With placement
629
630 dist=t1Ub3.DistanceToIn(pbigx,vmx);
631 assert(ApproxEqual(dist,50));
632
633 dist=t1Ub3.DistanceToIn(pbigy,vmy);
634 assert(ApproxEqual(dist,40));
635
636 dist=t1Ub3.DistanceToIn(pbigmy,vy);
637 assert(ApproxEqual(dist,50));
638
639 dist=t1Ub3.DistanceToIn(pbigmx,vx);
640 assert(ApproxEqual(dist,50));
641
642 dist=t2Ut1Trans.DistanceToIn(pzero,vx);
643 // G4cout<<"t2Ut1Trans.DistanceToIn(pzero,vx) = "<<dist<<G4endl ;
644 assert(ApproxEqual(dist,45));
645
646 dist=t2Ut1Trans.DistanceToIn(pzero,vz);
647 // G4cout<<"t2Ut1Trans.DistanceToIn(pzero,vz) = "<<dist<<G4endl ;
648 assert(ApproxEqual(dist,60));
649
650 dist=t2Ut1Trans.DistanceToIn(pzero,vmz);
651 // G4cout<<"t2Ut1Trans.DistanceToIn(pzero,vmz) = "<<dist<<G4endl ;
652 assert(ApproxEqual(dist,kInfinity));
653
654 dist=t2Ut1Trans.DistanceToIn(G4ThreeVector(0,0,55),vx);
655 // G4cout<<"t2Ut1Trans.DistanceToIn(G4ThreeVector(0,0,55),vx) = "
656 // <<dist<<G4endl ;
657 assert(ApproxEqual(dist,kInfinity));
658
659 dist=t2Ut1Ut1.DistanceToIn(pzero,vmz);
660 // G4cout<<"t2Ut1Ut1.DistanceToIn(pzero,vmz) = "<<dist<<G4endl ;
661 assert(ApproxEqual(dist,60));
662
663
664 G4cout<<"Tracking functions are OK"<<G4endl ;
665
666
667// CalculateExtent
668
669 G4VoxelLimits limit; // Unlimited
670 G4RotationMatrix noRot;
671 G4AffineTransform origin;
672 G4double min,max;
673
674 assert(b1.CalculateExtent(kXAxis,limit,origin,min,max));
675 assert(ApproxEqual(min,-20)&&ApproxEqual(max,20));
676
677 assert(b1.CalculateExtent(kYAxis,limit,origin,min,max));
678 assert(ApproxEqual(min,-30)&&ApproxEqual(max,30));
679
680 assert(b1.CalculateExtent(kZAxis,limit,origin,min,max));
681 assert(ApproxEqual(min,-40)&&ApproxEqual(max,40));
682
683 assert(b1Ub3.CalculateExtent(kXAxis,limit,origin,min,max));
684 assert(ApproxEqual(min,-20)&&ApproxEqual(max,20));
685
686 assert(b1Ub3.CalculateExtent(kYAxis,limit,origin,min,max));
687 assert(ApproxEqual(min,-30)&&ApproxEqual(max,60));
688
689 assert(b1Ub3.CalculateExtent(kZAxis,limit,origin,min,max));
690 assert(ApproxEqual(min,-50)&&ApproxEqual(max,50));
691
692 G4ThreeVector pmxmymz(-100,-110,-120);
693 G4AffineTransform tPosOnly(pmxmymz);
694
695 assert(b1.CalculateExtent(kXAxis,limit,tPosOnly,min,max));
696 assert(ApproxEqual(min,-120)&&ApproxEqual(max,-80));
697
698 assert(b1.CalculateExtent(kYAxis,limit,tPosOnly,min,max));
699 assert(ApproxEqual(min,-140)&&ApproxEqual(max,-80));
700
701 assert(b1.CalculateExtent(kZAxis,limit,tPosOnly,min,max));
702 assert(ApproxEqual(min,-160)&&ApproxEqual(max,-80));
703
704 G4RotationMatrix r90Z;
705 r90Z.rotateZ(pi/2);
706 G4AffineTransform tRotZ(r90Z,pzero);
707
708 assert(b1.CalculateExtent(kXAxis,limit,tRotZ,min,max));
709 assert(ApproxEqual(min,-30)&&ApproxEqual(max,30));
710
711 assert(b1.CalculateExtent(kYAxis,limit,tRotZ,min,max));
712 assert(ApproxEqual(min,-20)&&ApproxEqual(max,20));
713
714 assert(b1.CalculateExtent(kZAxis,limit,tRotZ,min,max));
715 assert(ApproxEqual(min,-40)&&ApproxEqual(max,40));
716
717// Check that clipped away
718
719 G4VoxelLimits xClip;
720 xClip.AddLimit(kXAxis,-100,-50);
721
722 assert(!b1.CalculateExtent(kXAxis,xClip,origin,min,max));
723
724// Assert clipped to volume
725
726 G4VoxelLimits allClip;
727 allClip.AddLimit(kXAxis,-5,+5);
728 allClip.AddLimit(kYAxis,-5,+5);
729 allClip.AddLimit(kZAxis,-5,+5);
730
731 G4RotationMatrix genRot;
732 genRot.rotateX(pi/6);
733 genRot.rotateY(pi/6);
734 genRot.rotateZ(pi/6);
735 G4AffineTransform tGen(genRot,vx);
736
737 assert(b1.CalculateExtent(kXAxis,allClip,tGen,min,max));
738 assert(ApproxEqual(min,-5)&&ApproxEqual(max,5));
739
740 assert(b1.CalculateExtent(kYAxis,allClip,tGen,min,max));
741 assert(ApproxEqual(min,-5)&&ApproxEqual(max,5));
742
743 assert(b1.CalculateExtent(kZAxis,allClip,tGen,min,max));
744 assert(ApproxEqual(min,-5)&&ApproxEqual(max,5));
745
746 G4VoxelLimits buggyClip2;
747 buggyClip2.AddLimit(kXAxis,5,15);
748
749 assert(b1.CalculateExtent(kXAxis,buggyClip2,origin,min,max));
750 assert(ApproxEqual(min,5)&&ApproxEqual(max,15));
751
752 assert(b1.CalculateExtent(kYAxis,buggyClip2,origin,min,max));
753 assert(ApproxEqual(min,-30)&&ApproxEqual(max,30));
754
755 assert(b1.CalculateExtent(kZAxis,buggyClip2,origin,min,max));
756 assert(ApproxEqual(min,-40)&&ApproxEqual(max,40));
757
758 buggyClip2.AddLimit(kYAxis,5,15);
759
760 assert(b1.CalculateExtent(kXAxis,buggyClip2,origin,min,max));
761 assert(ApproxEqual(min,5)&&ApproxEqual(max,15));
762
763 assert(b1.CalculateExtent(kYAxis,buggyClip2,origin,min,max));
764 assert(ApproxEqual(min,5)&&ApproxEqual(max,15));
765
766 assert(b1.CalculateExtent(kZAxis,buggyClip2,origin,min,max));
767 assert(ApproxEqual(min,-40)&&ApproxEqual(max,40));
768
769 G4VoxelLimits buggyClip1;
770 buggyClip1.AddLimit(kXAxis,-5,+5);
771
772 assert(b1.CalculateExtent(kXAxis,buggyClip1,origin,min,max));
773 assert(ApproxEqual(min,-5)&&ApproxEqual(max,5));
774
775 assert(b1.CalculateExtent(kYAxis,buggyClip1,origin,min,max));
776 assert(ApproxEqual(min,-30)&&ApproxEqual(max,30));
777
778 assert(b1.CalculateExtent(kZAxis,buggyClip1,origin,min,max));
779 assert(ApproxEqual(min,-40)&&ApproxEqual(max,40));
780
781 buggyClip1.AddLimit(kYAxis,-5,+5);
782
783 assert(b1.CalculateExtent(kXAxis,buggyClip1,origin,min,max));
784 assert(ApproxEqual(min,-5)&&ApproxEqual(max,5));
785
786 assert(b1.CalculateExtent(kYAxis,buggyClip1,origin,min,max));
787 assert(ApproxEqual(min,-5)&&ApproxEqual(max,5));
788
789 assert(b1.CalculateExtent(kZAxis,buggyClip1,origin,min,max));
790 assert(ApproxEqual(min,-40)&&ApproxEqual(max,40));
791
792 G4cout<<"CalculateExtent is OK"<<G4endl ;
793
794
795 return 0 ;
796}
Note: See TracBrowser for help on using the repository browser.