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

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