source: trunk/source/geometry/volumes/test/testExitNormalNav.cc@ 1337

Last change on this file since 1337 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: 13.2 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: testExitNormalNav.cc,v 1.7 2006/06/29 18:58:25 gunter Exp $
28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
29//
30// Testing the product of Exit Normal of the Navigator for
31// simple hierarchial geometry.
32// ( replicas, parameterised volumes currently not included )
33//
34// First version: J. Apostolakis, 18th June 2002
35
36#include <assert.h>
37#include "ApproxEqual.hh"
38
39// Global defs
40#include "globals.hh"
41
42#include "G4LogicalVolume.hh"
43#include "G4VPhysicalVolume.hh"
44#include "G4PVPlacement.hh"
45// #include "G4PVParameterised.hh"
46// #include "G4VPVParameterisation.hh"
47#include "G4Box.hh"
48
49#include "G4GeometryManager.hh"
50
51#include "G4RotationMatrix.hh"
52#include "G4ThreeVector.hh"
53
54// Build simple geometry:
55// 6 small cubes inside a slab (all G4Boxes)
56// 3 slabs are positioned inside the world cuboid
57
58G4VPhysicalVolume* BuildGeometry()
59{
60
61 // Rotations in X
62 G4RotationMatrix *prot90d_X, *prot180d_X, *prot270d_X;
63 prot90d_X = new G4RotationMatrix();
64 prot180d_X = new G4RotationMatrix();
65 prot270d_X = new G4RotationMatrix();
66 prot90d_X->rotateX(pi*0.5);
67 prot180d_X->rotateX(pi);
68 prot270d_X->rotateX(pi*1.5);
69
70 // Rotations in Y
71 G4RotationMatrix *prot90d_Y, *prot180d_Y, *prot270d_Y;
72 prot90d_Y = new G4RotationMatrix();
73 prot180d_Y = new G4RotationMatrix();
74 prot270d_Y = new G4RotationMatrix();
75 prot90d_Y->rotateY(pi*0.5);
76 prot180d_Y->rotateY(pi);
77 prot270d_Y->rotateY(pi*1.5);
78
79 // Rotations in Z
80 G4RotationMatrix *prot90d_Z, *prot180d_Z, *prot270d_Z;
81 prot90d_Z = new G4RotationMatrix();
82 prot180d_Z = new G4RotationMatrix();
83 prot270d_Z = new G4RotationMatrix();
84 prot90d_Z->rotateZ(pi*0.5);
85 prot180d_Z->rotateZ(pi);
86 prot270d_Z->rotateZ(-pi*0.5);
87
88 // Solids
89 G4Box *myBigBox= new G4Box("BigBox-World",200.*cm,200.*cm,200.*cm);
90 G4Box *Slab= new G4Box("slab",17.5*cm,10.*cm,7.5*cm);
91 G4Box *inCube10= new G4Box("Cube ten",5.*cm,5.*cm,5.*cm);
92 G4Box *smallCube= new G4Box("Small cube", 0.5*cm, 0.5*cm, 0.5*cm);
93
94 // World
95 G4LogicalVolume *worldLog=new G4LogicalVolume(myBigBox,0,
96 "WorldLV",0,0,0);
97 // Logical with no material,field,
98 // sensitive detector or user limits
99
100 G4PVPlacement *worldPhys=new G4PVPlacement(0,G4ThreeVector(0,0,0),
101 "WorldPV",worldLog,
102 0,false,0);
103 // Note: no mother pointer set
104
105 // Slab volume
106 G4LogicalVolume *slabLog=new G4LogicalVolume(Slab, 0, "Slab-logV");
107
108 // Inner volume
109 G4LogicalVolume *boxLog=
110 new G4LogicalVolume(inCube10, 0, "Cube10-logV");
111
112 // Smallest cube volume
113 G4LogicalVolume *smallLog=
114 new G4LogicalVolume(smallCube, 0, "smallCube1-lV");
115
116 // Place small cubes inside Cube ten
117 new G4PVPlacement(prot90d_Y,
118 G4ThreeVector( -4.5*cm, 0.0, 0.0),
119 smallLog,
120 "smallBackX",
121 boxLog,
122 false,
123 0);
124 new G4PVPlacement(prot180d_Y,
125 G4ThreeVector( 4.5*cm, 0.0, 0.0),
126 smallLog,
127 "smallFrontX",
128 boxLog,
129 false,
130 1);
131 new G4PVPlacement(prot90d_Z,
132 G4ThreeVector( 0.0, -4.5*cm, 0.0),
133 smallLog,
134 "smallBackY",
135 boxLog,
136 false,
137 2);
138 new G4PVPlacement(prot90d_X,
139 G4ThreeVector( 0.0, 4.5*cm, 0.0),
140 smallLog,
141 "smallFrontY",
142 boxLog,
143 false,
144 3);
145
146 // Fill the slab with inner volumes
147 //
148
149 G4ThreeVector centerPositionFirst(12.5*cm,-5*cm,0.0);
150 new G4PVPlacement(prot90d_Y,
151 centerPositionFirst,
152 boxLog,
153 "Lower Front",
154 slabLog,
155 false,
156 0);
157
158 G4ThreeVector centerPositionSecond(12.5*cm, 5*cm,0.0);
159 new G4PVPlacement(prot180d_Z,
160 centerPositionSecond,
161 boxLog,
162 "Upper Front",
163 slabLog,
164 false,
165 1);
166
167 G4ThreeVector centerPositionThird(-12.5*cm, 5*cm,0.0);
168 new G4PVPlacement(prot270d_X,
169 centerPositionThird,
170 boxLog,
171 "Upper Back",
172 slabLog,
173 false,
174 2);
175
176
177 G4ThreeVector centerPositionFourth(-12.0*cm, -5*cm,0.0);
178 new G4PVPlacement(prot180d_Z,
179 centerPositionFourth,
180 boxLog,
181 "Lower Back",
182 slabLog,
183 false,
184 3);
185
186
187 G4ThreeVector centerPositionFifth(-2.5*cm, 5*cm,0.0);
188 new G4PVPlacement(prot90d_Y,
189 centerPositionFifth,
190 boxLog,
191 "Upper Mid-Back",
192 slabLog,
193 false,
194 4);
195
196
197 G4ThreeVector centerPositionSixth( 2.5*cm, -5*cm,0.0);
198 new G4PVPlacement(prot90d_X,
199 centerPositionSixth,
200 boxLog,
201 "Lower Mid-Front",
202 slabLog,
203 false,
204 5);
205
206 // Placement of Slabs in World Volume
207 //
208 G4ThreeVector slabPositionOne( -27.5*cm, 0.0,0.0);
209 new G4PVPlacement(prot180d_Y,
210 slabPositionOne,
211 "Back-Slab1",
212 slabLog,
213 worldPhys,
214 false,
215 1);
216
217 G4ThreeVector slabPositionTwo( 0.0, 0.0, 0.0);
218 new G4PVPlacement(prot90d_Z,
219 slabPositionTwo,
220 "Upright-Middle-Slab2",
221 slabLog,
222 worldPhys,
223 false,
224 2);
225
226
227 G4ThreeVector slabPositionThree( 27.5*cm, 0.0, 0.0);
228 new G4PVPlacement(prot180d_Z,
229 slabPositionThree,
230 "Front-Slab3",
231 slabLog,
232 worldPhys,
233 false,
234 3);
235
236 return worldPhys;
237}
238
239//
240// Test LocateGlobalPointAndSetup
241//
242G4bool testG4Navigator1(G4VPhysicalVolume *pTopNode)
243{
244 MyNavigator myNav;
245 G4VPhysicalVolume *located;
246 myNav.SetWorldVolume(pTopNode);
247
248 assert(!myNav.LocateGlobalPointAndSetup(G4ThreeVector(1000*cm,0,0),0,false));
249
250 located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(1800.0*mm,0,0),0,false);
251 assert(located->GetName()=="WorldPV");
252
253 assert(!myNav.LocateGlobalPointAndSetup(G4ThreeVector(1000.*cm,0,0)));
254
255 return true;
256
257 // Can add more location checks here, like the old ones below.
258
259// Check relative search that causes backup one level and then search down:
260// Nonrel' finds Target 3, then rel' with point in Target 5 finds Target 5
261 located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(15.*cm,0,-5.*cm),0,false);
262 assert(located->GetName()=="Upper Front");
263
264 located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(0,-15.*cm,20.*cm));
265 assert(located->GetName()=="Target 5");
266 assert(ApproxEqual(myNav.CurrentLocalCoordinate(),G4ThreeVector(0,0,10)));
267// Check that outside point causes stack to unwind
268 assert(!myNav.LocateGlobalPointAndSetup(G4ThreeVector(kInfinity,0,0)));
269
270 return true;
271}
272
273int verbose= 1;
274//
275// Test Stepping
276//
277G4bool testExitNormal(G4VPhysicalVolume *pTopNode,
278 G4ThreeVector initialPoint,
279 G4ThreeVector direction,
280 G4ThreeVector expectedExitNorm)
281{
282 MyNavigator myNav;
283 G4VPhysicalVolume *located;
284 G4double Step,physStep,safety;
285
286 myNav.SetWorldVolume(pTopNode);
287//
288// Test location & Step computation
289//
290 G4ThreeVector initPoint(initialPoint), newPoint(0,0,0);
291 // G4ThreeVector direction= xHat;
292 G4bool valid;
293
294 if( verbose ){
295 G4cout << "Initial step " << G4endl;
296 G4cout << "-Initial Point = " << initPoint << G4endl;
297 }
298
299 located=myNav.LocateGlobalPointAndSetup(initPoint);
300 assert(located->GetName()=="WorldPV");
301 if( verbose )
302 G4cout << "-Located: Location before is " << located->GetName() << G4endl;
303
304 physStep=kInfinity;
305 Step=myNav.ComputeStep(initPoint, direction, physStep, safety);
306 if( verbose ){
307 G4cout << "-Moved: Step was = " << Step << " expected " << 5.0 * cm << G4endl;
308 G4cout << " safety= " << safety << G4endl;
309 }
310 assert(ApproxEqual(Step,5.0*cm));
311 // assert(ApproxEqual(safety,50.0));
312 assert(safety>=0.0);
313 assert(safety<=50.0);
314
315 newPoint= initPoint + Step * direction;
316 G4ThreeVector localNormal = myNav.GetLocalExitNormal(&valid);
317 assert(valid);
318 G4ThreeVector globalNormal = myNav.GetLocalToGlobalTransform().TransformAxis(localNormal);
319 assert( globalNormal == expectedExitNorm );
320
321 myNav.SetGeometricallyLimitedStep();
322 located=myNav.LocateGlobalPointAndSetup(initPoint);
323 assert(located->GetName()!="WorldPV");
324 if( verbose )
325 G4cout << "-Located: Location is " << located->GetName() << G4endl;
326
327 // Next Steps
328 G4int istep;
329 for ( istep=0; istep < 15; istep++ ){
330
331 initPoint= newPoint;
332
333 if( verbose ){
334 G4cout << "Sub step " << istep << G4endl;
335 G4cout << "-Initial Point = " << initPoint << G4endl;
336 G4cout << "-Location before is " << located->GetName() << G4endl;
337 }
338
339 physStep=kInfinity;
340 Step=myNav.ComputeStep(initPoint, direction, physStep, safety);
341 if( verbose )
342 G4cout << "-Moved: Step was = " << Step << G4endl;
343 assert( Step <= 10.0*cm);
344 assert(ApproxEqual(safety,0.0));
345 assert(safety>=0);
346
347 newPoint= initPoint + Step * direction;
348 G4ThreeVector localNormal = myNav.GetLocalExitNormal(&valid);
349 assert(valid);
350 G4ThreeVector globalNormal = myNav.GetLocalToGlobalTransform().TransformAxis(localNormal);
351
352 if( 0 ) { // globalNormal != G4ThreeVector(1.0,0.0,0.0) ){
353 G4cout << " **Problem** with pre-relocation normals: " << G4endl;
354 G4cout << " *Point = " << newPoint << G4endl;
355 G4cout << " *localNorm = " << localNormal << G4endl;
356 G4cout << " *globalNorm = " << globalNormal << G4endl;
357 }
358 // assert( globalNormal == G4ThreeVector(1.0,0.0,0.0) );
359
360 myNav.SetGeometricallyLimitedStep();
361 located=myNav.LocateGlobalPointAndSetup(newPoint);
362 // assert(located->GetName()!="WorldPV");
363 if( verbose )
364 G4cout << "-Located: Location after is " << located->GetName() << G4endl;
365
366 localNormal = myNav.GetLocalExitNormal(&valid);
367 assert(valid);
368 globalNormal = myNav.GetLocalToGlobalTransform().TransformAxis(localNormal);
369 if( verbose ) {
370 G4cout << "Post-relocation normals: " << G4endl;
371 G4cout << " Point = " << newPoint << G4endl;
372 G4cout << " Location after is " << located->GetName() << G4endl;
373 G4cout << " localNorm = " << localNormal << G4endl;
374 G4cout << " globalNorm = " << globalNormal << G4endl;
375 G4cout << " expectedExitNorm = " << expectedExitNorm << G4endl;
376 }
377 // assert( ApproxEqual( globalNormal, expectedExitNorm ) );
378
379 }
380
381 return true;
382}
383
384G4bool testExitNormalNav()
385{
386 G4VPhysicalVolume *myTopNode;
387 const G4ThreeVector xHat(1,0,0),yHat(0,1,0),zHat(0,0,1);
388 const G4ThreeVector mxHat(-1,0,0),myHat(0,-1,0),mzHat(0,0,-1);
389
390 G4ThreeVector initPointMinusX(-50.0*cm,0.01*cm,0.);
391 G4ThreeVector initPointPluxX(50.0*cm, -0.01*cm,0.);
392
393 myTopNode=BuildGeometry(); // Build the geometry
394 G4GeometryManager::GetInstance()->CloseGeometry(false);
395 testG4Navigator1(myTopNode);
396 testExitNormal(myTopNode, initPointMinusX, xHat, xHat);
397 testExitNormal(myTopNode, initPointPluxX, mxHat, mxHat);
398 testExitNormal(myTopNode, G4ThreeVector(-50.0*cm,2.0*cm,0.0), xHat, xHat);
399 testExitNormal(myTopNode, G4ThreeVector(-50.0*cm,-2.0*cm,0.0), xHat, xHat);
400
401// Repeat tests but with full voxels
402 G4GeometryManager::GetInstance()->OpenGeometry();
403 G4GeometryManager::GetInstance()->CloseGeometry(true);
404 testG4Navigator1(myTopNode);
405 testExitNormal(myTopNode, initPointMinusX, xHat, xHat);
406 testExitNormal(myTopNode, initPointPluxX, mxHat, mxHat);
407 testExitNormal(myTopNode, G4ThreeVector(-50.0*cm,2.0*cm,0.0), xHat, xHat);
408 testExitNormal(myTopNode, G4ThreeVector(-50.0*cm,-2.0*cm,0.0), xHat, xHat);
409 G4GeometryManager::GetInstance()->OpenGeometry();
410 return true;
411}
412
413int main()
414{
415 assert(testExitNormalNav());
416 return 0;
417}
418
Note: See TracBrowser for help on using the repository browser.