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

Last change on this file since 1342 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: 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.