source: trunk/source/geometry/solids/CSG/test/testG4Orb.cc@ 1350

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

geant4 tag 9.4
File size: 12.8 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: testG4Orb.cc,v 1.7 2009/01/29 16:54:51 grichine Exp $
28// GEANT4 tag $Name: geant4-09-04-ref-00 $
29//
30// G4Orb 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 G4Orb.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
49#include "G4ThreeVector.hh"
50#include "G4RotationMatrix.hh"
51#include "G4AffineTransform.hh"
52#include "G4VoxelLimits.hh"
53#include "G4Orb.hh"
54
55//const G4double kApproxEqualTolerance = kCarTolerance;
56
57// Return true if the double check is approximately equal to target
58//
59// Process:
60//
61// Return true is difference < kApproxEqualTolerance
62
63//G4bool ApproxEqual(const G4double check,const G4double target)
64//{
65// return (std::fabs(check-target)<kApproxEqualTolerance) ? true : false ;
66//}
67
68// Return true if the 3vector check is approximately equal to target
69//G4bool ApproxEqual(const G4ThreeVector& check, const G4ThreeVector& target)
70//{
71// return (ApproxEqual(check.x(),target.x())&&
72// ApproxEqual(check.y(),target.y())&&
73// ApproxEqual(check.z(),target.z()))? true : false;
74//}
75
76///////////////////////////////////////////////////////////////////
77//
78// Dave's auxiliary function
79
80const G4String OutputInside(const EInside a)
81{
82 switch(a)
83 {
84 case kInside: return "Inside";
85 case kOutside: return "Outside";
86 case kSurface: return "Surface";
87 }
88 return "????";
89}
90
91
92
93int main(void)
94{
95 G4double Dist, vol, volCheck;
96 G4ThreeVector pzero(0,0,0),px(30,0,0),py(0,30,0),pz(0,0,30);
97 G4ThreeVector Pmx(-30,0,0),pmy(0,-30,0),pmz(0,0,-30);
98 G4ThreeVector pbigx(100,0,0),pbigy(0,100,0),pbigz(0,0,100);
99 G4ThreeVector pbigmx(-100,0,0),pbigmy(0,-100,0),pbigmz(0,0,-100);
100
101 G4ThreeVector ponrmin1(45,0,0),ponrmax1(50,0,0),ponzmax(0,0,50),
102 ponrmin2(45/std::sqrt(2.),45/std::sqrt(2.),0),
103 ponrmin3(0,0,-45),ponrminJ(0,0,-300),ponrmaxJ(0,0,-500),
104 ponrmax2(50/std::sqrt(2.),50/std::sqrt(2.),0);
105 G4ThreeVector ponphi1(48/std::sqrt(2.),-48/std::sqrt(2.),0),
106 ponphi2(48/std::sqrt(2.),48/std::sqrt(2.),0),
107 pInPhi(48*0.866,-24,0),
108 pOverPhi(-48/std::sqrt(2.),48/std::sqrt(2.),0);
109 G4ThreeVector pontheta1(0,48*std::sin(pi/4),48*std::cos(pi/4)),
110 pontheta2(0,48*std::sin(pi/4),-48*std::cos(pi/4));
111
112 G4ThreeVector ptestphi1(-100,-45/std::sqrt(2.),0),
113 ptestphi2(-100,45/std::sqrt(2.),0);
114
115 G4ThreeVector ptesttheta1(0,48/std::sqrt(2.),100),
116 ptesttheta2(0,48/std::sqrt(2.),-100);
117
118 G4ThreeVector vx(1,0,0),vy(0,1,0),vz(0,0,1);
119 G4ThreeVector vmx(-1,0,0),vmy(0,-1,0),vmz(0,0,-1);
120 G4ThreeVector vxy(1/std::sqrt(2.),1/std::sqrt(2.),0),vmxmy(-1/std::sqrt(2.),-1/std::sqrt(2.),0);
121 G4ThreeVector vxmy(1/std::sqrt(2.),-1/std::sqrt(2.),0),vmxy(-1/std::sqrt(2.),1/std::sqrt(2.),0);
122 G4ThreeVector v345exit1(-0.8,0.6,0),v345exit2(0.8,0.6,0),
123 v345exit3(0.6,0.8,0);
124 G4ThreeVector norm,*pNorm;
125 G4bool *pgoodNorm,goodNorm,calcNorm=true;
126
127 pNorm=&norm;
128 pgoodNorm=&goodNorm;
129
130 G4Orb s1("Solid G4Orb",50);
131
132 G4Orb s10("s10",0.018*mm);
133
134 G4Orb* solidD1= new G4Orb("D1", 2.7*cm);
135
136 G4ThreeVector s10p(0.01160957408065766*mm,
137 0.01308205826682229*mm,
138 0.004293345210644617*mm);
139
140 G4ThreeVector positionKip(-5.1427112977805294,
141 25.37671986392073,
142 -7.6534050889634502);
143 G4ThreeVector directionKip(-0.90020960513809678,
144 -0.052042885743481759,
145 0.43233575477931813);
146
147 G4Orb G650("G650",5.0e-05*mm);
148
149 G4Orb b1046("b1046",4800*km);
150
151 G4ThreeVector positionG650(-2.61756492392351e-06*mm,
152 -4.992317579421979e-05*mm,
153 -9.474848146062698e-07*mm);
154 G4ThreeVector directionG650(0.6967354558785861,
155 0.716168068935222,
156 0.0407799161261308);
157
158
159
160
161
162 G4cout.precision(16);
163
164
165#ifdef NDEBUG
166 G4Exception("FAIL: *** Assertions must be compiled in! ***");
167#endif
168
169 //////////////// Check name /////////////////////////
170
171 assert(s1.GetName()=="Solid G4Orb");
172
173
174 // check cubic volume cases
175
176 vol = s1.GetCubicVolume();
177 volCheck = 4*pi*50*50*50/3;
178 assert(ApproxEqual(vol,volCheck));
179
180// Check G4Orb::Inside
181
182
183 EInside inside = s10.Inside(s10p);
184 G4cout<<"s10.Inside(s10p) = "
185 <<OutputInside(inside)<<G4endl ;
186 G4cout<<"p radius = "
187 <<s10p.mag()<<G4endl ;
188
189 inside = G650.Inside(positionG650);
190 G4cout<<"G650.Inside(positionG650) = "
191 <<OutputInside(inside)<<G4endl ;
192
193 assert(s1.Inside(pzero)==kInside);
194 // assert(s1.Inside(pz)==kInside);
195
196// Checking G4Orb::SurfaceNormal
197 norm=s1.SurfaceNormal(ponrmax1);
198 assert(ApproxEqual(norm,vx));
199
200// Checking G4Orb::DistanceToOut(P)
201 Dist=s1.DistanceToOut(pzero);
202 assert(ApproxEqual(Dist,50));
203 Dist=s1.DistanceToOut(ponrmax1);
204 assert(ApproxEqual(Dist,0));
205
206// Checking G4Orb::DistanceToOut(p,v)
207
208
209 Dist=s1.DistanceToOut(pz,vz,calcNorm,pgoodNorm,pNorm);
210 G4cout<<"Dist=s1.DistanceToOut(pz,vz) = "<<Dist<<G4endl;
211
212 Dist=s1.DistanceToOut(ponrmax1,vx,calcNorm,pgoodNorm,pNorm);
213 *pNorm=pNorm->unit();
214 assert(ApproxEqual(Dist,0)&&*pgoodNorm&&ApproxEqual(*pNorm,vx));
215
216
217 Dist=s1.DistanceToOut(pzero,vx,calcNorm,pgoodNorm,pNorm);
218 assert(ApproxEqual(Dist,50)&&ApproxEqual(pNorm->unit(),vx)&&*pgoodNorm);
219 Dist=s1.DistanceToOut(pzero,vmx,calcNorm,pgoodNorm,pNorm);
220 assert(ApproxEqual(Dist,50)&&ApproxEqual(pNorm->unit(),vmx)&&*pgoodNorm);
221 Dist=s1.DistanceToOut(pzero,vy,calcNorm,pgoodNorm,pNorm);
222 assert(ApproxEqual(Dist,50)&&ApproxEqual(pNorm->unit(),vy)&&*pgoodNorm);
223 Dist=s1.DistanceToOut(pzero,vmy,calcNorm,pgoodNorm,pNorm);
224 assert(ApproxEqual(Dist,50)&&ApproxEqual(pNorm->unit(),vmy)&&*pgoodNorm);
225 Dist=s1.DistanceToOut(pzero,vz,calcNorm,pgoodNorm,pNorm);
226 assert(ApproxEqual(Dist,50)&&ApproxEqual(pNorm->unit(),vz)&&*pgoodNorm);
227 Dist=s1.DistanceToOut(pzero,vmz,calcNorm,pgoodNorm,pNorm);
228 assert(ApproxEqual(Dist,50)&&ApproxEqual(pNorm->unit(),vmz)&&*pgoodNorm);
229 Dist=s1.DistanceToOut(pzero,vxy,calcNorm,pgoodNorm,pNorm);
230 assert(ApproxEqual(Dist,50)&&ApproxEqual(pNorm->unit(),vxy)&&*pgoodNorm);
231
232 Dist=s1.DistanceToOut(pzero,vmz,calcNorm,pgoodNorm,pNorm);
233 assert(ApproxEqual(Dist,50)&&ApproxEqual(pNorm->unit(),vmz)&&*pgoodNorm);
234 Dist=s1.DistanceToOut(pzero,vxy,calcNorm,pgoodNorm,pNorm);
235 assert(ApproxEqual(Dist,50)&&ApproxEqual(pNorm->unit(),vxy)&&*pgoodNorm);
236
237
238
239
240// Checking G4Orb::DistanceToIn(P)
241
242 Dist=s1.DistanceToIn(ponrmax1);
243 assert(ApproxEqual(Dist,0));
244
245// Checking G4Orb::DistanceToIn(P,V)
246
247 Dist=s1.DistanceToIn(ponzmax,vz);
248 G4cout<<"s1.DistanceToIn(ponzmax,vz) = "<<Dist<<G4endl;
249 Dist=s1.DistanceToIn(pbigy,vy);
250 assert(Dist==kInfinity);
251 Dist=s1.DistanceToIn(pbigy,vmy);
252 assert(ApproxEqual(Dist,50));
253 Dist=b1046.DistanceToIn(G4ThreeVector(0.,0.,4800*km),vmz);
254 G4cout<<"b1046.DistanceToIn(G4ThreeVector(0.,0.,4800*km),vmz) = "<<Dist<<G4endl;
255 assert(ApproxEqual(Dist,0));
256
257 // Checking In/Out both return zeros
258
259 Dist = solidD1->DistanceToIn(positionKip, directionKip);
260 G4cout << " Dist In (pos, uKipp) = " << Dist << G4endl;
261
262 Dist = solidD1->DistanceToOut(positionKip, directionKip);
263 G4cout << " Dist Out (pos, uKipp) = " << Dist << G4endl;
264
265 G4cout << " Dot( pos, uKipp) = " << positionKip.dot(directionKip)<< G4endl;
266
267 inside = solidD1->Inside(positionKip);
268 G4cout<<"solidD1->Inside(positionKip) = " << OutputInside(inside)<<G4endl ;
269
270
271
272 ///////////////////////////////////////////////////////////////////////////
273
274
275// CalculateExtent
276 G4VoxelLimits limit; // Unlimited
277 G4RotationMatrix noRot;
278 G4AffineTransform origin;
279 G4double min,max;
280 assert(s1.CalculateExtent(kXAxis,limit,origin,min,max));
281 assert(min<=-50&&max>=50);
282 assert(s1.CalculateExtent(kYAxis,limit,origin,min,max));
283 assert(min<=-50&&max>=50);
284 assert(s1.CalculateExtent(kZAxis,limit,origin,min,max));
285 assert(min<=-50&&max>=50);
286
287 G4ThreeVector pmxmymz(-100,-110,-120);
288 G4AffineTransform tPosOnly(pmxmymz);
289 assert(s1.CalculateExtent(kXAxis,limit,tPosOnly,min,max));
290 assert(min<=-150&&max>=-50);
291 assert(s1.CalculateExtent(kYAxis,limit,tPosOnly,min,max));
292 assert(min<=-160&&max>=-60);
293 assert(s1.CalculateExtent(kZAxis,limit,tPosOnly,min,max));
294 assert(min<=-170&&max>=-70);
295
296
297 G4RotationMatrix r90Z;
298 r90Z.rotateZ(halfpi);
299 G4AffineTransform tRotZ(r90Z,pzero);
300 assert(s1.CalculateExtent(kXAxis,limit,tRotZ,min,max));
301 assert(min<=-50&&max>=50);
302 assert(s1.CalculateExtent(kYAxis,limit,tRotZ,min,max));
303 assert(min<=-50&&max>=50);
304 assert(s1.CalculateExtent(kZAxis,limit,tRotZ,min,max));
305 assert(min<=-50&&max>=50);
306
307
308// Check that clipped away
309 G4VoxelLimits xClip;
310 xClip.AddLimit(kXAxis,-100,-60);
311 assert(!s1.CalculateExtent(kXAxis,xClip,origin,min,max));
312
313// Assert clipped to volume
314 G4VoxelLimits allClip;
315 allClip.AddLimit(kXAxis,-5,+5);
316 allClip.AddLimit(kYAxis,-5,+5);
317 allClip.AddLimit(kZAxis,-5,+5);
318 G4RotationMatrix genRot;
319 genRot.rotateX(pi/6);
320 genRot.rotateY(pi/6);
321 genRot.rotateZ(pi/6);
322 G4AffineTransform tGen(genRot,vx);
323 assert(s1.CalculateExtent(kXAxis,allClip,tGen,min,max));
324 assert(min<=-5&&max>=5);
325 assert(s1.CalculateExtent(kYAxis,allClip,tGen,min,max));
326 assert(min<=-5&&max>=5);
327 assert(s1.CalculateExtent(kZAxis,allClip,tGen,min,max));
328 assert(min<=-5&&max>=5);
329
330
331// Test z clipping ok
332 for (G4double zTest=-100;zTest<100;zTest+=9)
333 {
334 G4VoxelLimits zTestClip;
335 zTestClip.AddLimit(kZAxis,-kInfinity,zTest);
336 if (zTest<-50)
337 {
338 assert(!s1.CalculateExtent(kZAxis,zTestClip,origin,min,max));
339 }
340 else
341 {
342 assert(s1.CalculateExtent(kZAxis,zTestClip,origin,min,max));
343
344 G4double testMin=-50;
345 G4double testMax=(zTest<50) ? zTest : 50;
346 assert (ApproxEqual(min,testMin) && ApproxEqual(max,testMax));
347 }
348 }
349
350// Test y clipping ok
351 for (G4double xTest=-100;xTest<100;xTest+=9)
352 {
353 G4VoxelLimits xTestClip;
354 xTestClip.AddLimit(kXAxis,-kInfinity,xTest);
355 if (xTest<-50)
356 {
357 assert(!s1.CalculateExtent(kYAxis,xTestClip,origin,min,max));
358 }
359 else
360 {
361 assert(s1.CalculateExtent(kYAxis,xTestClip,origin,min,max));
362// Calc max y coordinate
363 G4double testMax=(xTest<0) ? std::sqrt(50*50-xTest*xTest) : 50;
364 assert (ApproxEqual(min,-testMax) && ApproxEqual(max,testMax));
365 }
366 }
367
368// Test x clipping ok
369 for (G4double yTest=-100;yTest<100;yTest+=9)
370 {
371 G4VoxelLimits yTestClip;
372 yTestClip.AddLimit(kYAxis,-kInfinity,yTest);
373 if (yTest<-50)
374 {
375 assert(!s1.CalculateExtent(kXAxis,yTestClip,origin,min,max));
376 }
377 else
378 {
379 assert(s1.CalculateExtent(kXAxis,yTestClip,origin,min,max));
380// Calc max y coordinate
381 G4double testMax=(yTest<0) ? std::sqrt(50*50-yTest*yTest) : 50;
382 assert (ApproxEqual(min,-testMax) && ApproxEqual(max,testMax));
383 }
384 }
385 return 0;
386}
387
388
389
390
391
392
393
394
395
396
397
Note: See TracBrowser for help on using the repository browser.