source: trunk/source/geometry/navigation/test/testG4PathFinder.cc@ 1339

Last change on this file since 1339 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: 18.0 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: testG4PathFinder.cc,v 1.8 2007/02/13 16:15:34 japost Exp $
28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
29//
30//
31// Moving with PathFinder in simple boxlike geometry,
32
33// TO DO // both with and without voxels.
34
35#include <assert.h>
36#include "ApproxEqual.hh"
37
38// Global defs
39#include "globals.hh"
40
41// To create the geometry
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// For the unit under test
55// #include "G4Navigator.hh"
56#include "G4PathFinder.hh"
57
58#include <iomanip>
59
60G4bool printing= true; // false for minimum output, true to debug/check
61
62// Build simple geometry:
63// 4 small cubes + 1 slab (all G4Boxes) are positioned inside a larger cuboid
64G4VPhysicalVolume* BuildGeometry()
65{
66
67 G4Box *myBigBox= new G4Box ("BixBox",250*cm,250*cm,200*cm);
68
69 G4LogicalVolume *worldLog=
70 new G4LogicalVolume(myBigBox,0, "World",0,0,0);
71 // No material, field, sensitive detector or user limits
72 G4PVPlacement *worldPhys=
73 new G4PVPlacement(0,G4ThreeVector(0,0,0),
74 "World",worldLog, 0,false,0);
75 // No displacement, no mother pointer - world!
76
77 G4Box *myBox=new G4Box("cube-10",10*cm,10*cm,10*cm);
78 G4LogicalVolume *boxLog=
79 new G4LogicalVolume(myBox,0, "Crystal Box",0,0,0);
80
81 new G4PVPlacement(0,G4ThreeVector(15*cm,0,0), "Target 01",boxLog,
82 worldPhys,false,1);
83 new G4PVPlacement(0,G4ThreeVector(25*cm,0,0), "Target 02",boxLog,
84 worldPhys,false,2);
85 new G4PVPlacement(0,G4ThreeVector(-15*cm,0,0), "Target 03",boxLog,
86 worldPhys,false,3);
87 new G4PVPlacement(0,G4ThreeVector(-25*cm,0,0), "Target 04",boxLog,
88 worldPhys,false,4);
89 new G4PVPlacement(0,G4ThreeVector(0,15*cm,0), "Target 11",boxLog,
90 worldPhys,false,11);
91 new G4PVPlacement(0,G4ThreeVector(0,25*cm,0), "Target 12",boxLog,
92 worldPhys,false,12);
93 new G4PVPlacement(0,G4ThreeVector(0,-15*cm,0), "Target 13",boxLog,
94 worldPhys,false,13);
95 new G4PVPlacement(0,G4ThreeVector(0,-25*cm,0), "Target 14",boxLog,
96 worldPhys,false,14);
97 new G4PVPlacement(0,G4ThreeVector(0,0,15*cm), "Target 21",boxLog,
98 worldPhys,false,21);
99 new G4PVPlacement(0,G4ThreeVector(0,0,25*cm), "Target 22",boxLog,
100 worldPhys,false,22);
101 new G4PVPlacement(0,G4ThreeVector(0,0,-15*cm), "Target 23",boxLog,
102 worldPhys,false,23);
103 new G4PVPlacement(0,G4ThreeVector(0,0,-25*cm), "Target 24",boxLog,
104 worldPhys,false,24);
105
106
107 new G4PVPlacement(0,G4ThreeVector(-15*cm,15*cm,-10*cm), "Target 101",boxLog,
108 worldPhys,false,1);
109 new G4PVPlacement(0,G4ThreeVector(-15*cm,-15*cm,-10*cm), "Target 102",boxLog,
110 worldPhys,false,2);
111
112 G4Box *mySlab= new G4Box("slab",10*cm,25*cm,10*cm);
113 G4LogicalVolume *slabLog=new G4LogicalVolume(mySlab,0,
114 "Crystal Slab",0,0,0);
115 // G4PVPlacement *offYPhys=
116 new G4PVPlacement(0,G4ThreeVector(15*cm,0,-10*cm), "Target 3",slabLog,
117 worldPhys,false,3);
118
119 new G4PVPlacement(0,G4ThreeVector(0,15*cm,10*cm), "Target 4",boxLog,
120 worldPhys,false,4);
121
122 new G4PVPlacement(0,G4ThreeVector(0,-15*cm,10*cm), "Target 5",boxLog,
123 worldPhys,false,5);
124
125
126 return worldPhys;
127}
128
129#include "G4TransportationManager.hh"
130
131G4int navIdExpected= 0; // New convention in G4TransportationManager
132
133//
134// Test LocateGlobalPointAndSetup
135//
136G4PathFinder* setupPathFinder(G4VPhysicalVolume *pTopNode)
137{
138 MyNavigator myNav;
139 G4Navigator* pNav;
140
141 G4int navId;
142 G4bool useMyNav= false;
143 // G4bool overwriteNav= false; // Default
144 G4bool overwriteNav= true; // Second option 'new'
145
146 pNav= new G4Navigator(); // (&myNav);
147 // pNav= (&myNav);
148
149 if( useMyNav ){
150 pNav= (&myNav);
151 myNav.SetWorldVolume(pTopNode);
152 }else{
153 pNav->SetWorldVolume(pTopNode);
154 }
155
156 G4PathFinder* pathFinder= G4PathFinder::GetInstance();
157 //=========== --------------------------
158 static G4TransportationManager* transportMgr=
159 G4TransportationManager::GetTransportationManager();
160
161 G4Navigator*
162 origNav= transportMgr->GetNavigatorForTracking();
163 origNav->SetWorldVolume(pTopNode);
164
165 navId= transportMgr->ActivateNavigator( origNav );
166 G4cout << " navId for original Navigator for tracking is " << navId << G4endl;
167
168 if( overwriteNav ){
169 transportMgr->DeActivateNavigator( origNav );
170 G4cout << " DeActivated Navigator " << origNav << G4endl;
171
172 // G4bool registered= transportMgr->RegisterNavigator( pNav );
173 // assert( registered );
174 transportMgr->SetNavigatorForTracking(pNav);
175 G4cout << " Setting new Navigator for Tracking " << pNav << G4endl;
176 // transportMgr->SetWorldVolume( pTopNode );
177
178 navId= transportMgr->ActivateNavigator( pNav );
179 G4cout << " navId for new Navigator for tracking is " << navId << G4endl;
180 assert ( navId == navIdExpected );
181 }
182
183 return pathFinder;
184}
185
186G4bool testG4PathFinder1(G4VPhysicalVolume *) // pTopNode)
187{
188 G4VPhysicalVolume *located;
189 static G4TransportationManager* transportMgr=
190 G4TransportationManager::GetTransportationManager();
191 static G4PathFinder* pathFinder= G4PathFinder::GetInstance();
192
193 G4Navigator *pNav= transportMgr->GetNavigatorForTracking();
194
195 G4int navId= navIdExpected; // As checked/ asserted above !
196
197 G4ThreeVector position( 0., 0., 0.), dirUx(1.,0.,0.);
198 pathFinder->PrepareNewTrack( position, dirUx );
199 //========== ---------------
200 // Also locates !!
201
202 // Need the volume information
203 // pathFinder->Locate( endPoint, endDirection );
204 //========== ------
205 located= pathFinder->GetLocatedVolume( navId );
206
207 G4double t0=0.0, Ekin=100.00*MeV, restMass= 0.511*MeV, charge= -1, magDipole=0.0, s=0.0;
208 G4ThreeVector Spin(0.,0.,0.);
209 G4FieldTrack startFT( position, t0, dirUx, Ekin, restMass, charge, Spin, magDipole, s),
210 endFT('a'); // Default values for return
211
212 G4int stepNo=0;
213 G4double steplen= 123.00 * mm, safetyRet=0.0;
214 ELimited limited;
215 G4cout << " test: input FT = " << startFT << G4endl
216 << " "
217 << " len = " << steplen
218 << " navId = " << navId << " stepNo = " << stepNo << G4endl;
219 G4double stepdone=
220 pathFinder->ComputeStep( startFT, steplen, navId, stepNo, safetyRet, limited, endFT, located );
221 //========== -----------
222
223
224 G4ThreeVector endPoint = endFT.GetPosition();
225 G4ThreeVector endDirection= endFT.GetMomentumDirection();
226 pathFinder->Locate( endPoint, endDirection );
227 //========== ------
228 located= pathFinder->GetLocatedVolume( navId );
229 if (located && printing ) {
230 G4cout << " Located in volume " << located->GetName()
231 << " id= " << located->GetCopyNo() << G4endl;
232 }
233
234 G4double stepAgain=
235 pathFinder->ComputeStep( startFT, steplen, navId, stepNo, safetyRet, limited, endFT, located );
236 // Should not move, since the 'stepNo' is the same !!
237
238 G4double endSafety= pathFinder->ComputeSafety( endFT.GetPosition() );
239
240 G4cout.precision(5);
241 do{
242 startFT= endFT;
243
244 stepdone=
245 pathFinder->ComputeStep( startFT, steplen, navId, ++stepNo, safetyRet, limited, endFT, located );
246 pathFinder->Locate( endFT.GetPosition(), endFT.GetMomentumDirection() );
247 located= pathFinder->GetLocatedVolume( navId );
248
249 if( std::fabs(safetyRet-endSafety) > 1.0e-9 * safetyRet ) {
250 G4cerr << " Problem in safety at new point "
251 << " Last endpoint returned " << endSafety
252 << " while new computeStep gives "
253 << G4endl;
254 }
255
256 endSafety= pathFinder->ComputeSafety( endFT.GetPosition() );
257
258 G4double truestep= (stepdone < steplen ) ? stepdone : steplen;
259 if ( printing ) {
260 G4cout << " Step " << std::setw(3) << stepNo << " "
261 << " start-safety= " << std::setw(7) << safetyRet << " "
262 << " step-length= " << std::setw(7) << truestep / mm << " mm "
263 << " to " << endFT.GetPosition() << " " ;
264 if( located ) {
265 G4cout
266 << " new volume= " << std::setw(10) << located->GetName()
267 << " copyNo= " << located->GetCopyNo();
268 }
269 G4cout << " end-safety= " << endSafety;
270 if( located ) G4cout << G4endl;
271 }
272 } while ( located );
273
274 G4cout << " Last Step it exits the World. Located = " << located << G4endl;
275 G4cout << G4endl;
276 // G4cout << " Step " << stepNo << " is the last one, it exits the World." << G4endl;
277 return true;
278
279 //// OLD Checks ......
280 assert(!pNav->LocateGlobalPointAndSetup(G4ThreeVector(kInfinity,0,0),0,false));
281 located=pNav->LocateGlobalPointAndSetup(G4ThreeVector(0,0,0),0,false);
282 assert(located->GetName()=="World");
283
284
285 return true;
286
287 // transportMgr->DeActivateNavigator( pNav );
288 // ----- // transportMgr->DeRegisterNavigator( pNav ); // Cannot de-register Mass Navigator
289 // exit(1);
290
291 //-------------------- END --------------------------------
292 // G4VPhysicalVolume *located;
293 MyNavigator& myNav= dynamic_cast<MyNavigator&>(*pNav);
294
295 assert(!myNav.LocateGlobalPointAndSetup(G4ThreeVector(kInfinity,0,0),0,false));
296 located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(0,0,0),0,false);
297 assert(located->GetName()=="World");
298
299 assert(!myNav.LocateGlobalPointAndSetup(G4ThreeVector(kInfinity,0,0)));
300
301// Check relative search that causes backup one level and then search down:
302// Nonrel' finds Target 3, then rel' with point in Target 5 finds Target 5
303 located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(15,0,-10),0,false);
304 assert(located->GetName()=="Vari' Blocks");
305
306 located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(0,-15,20));
307 assert(located->GetName()=="Target 5");
308 assert(ApproxEqual(myNav.CurrentLocalCoordinate(),G4ThreeVector(0,0,10)));
309// Check that outside point causes stack to unwind
310 assert(!myNav.LocateGlobalPointAndSetup(G4ThreeVector(kInfinity,0,0)));
311
312// Check parameterised volumes
313
314// Replication 0
315 located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(15,-15,-10));
316 assert(located->GetName()=="Vari' Blocks");
317 located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(15,-15,-16));
318 assert(located->GetName()=="Target 3");
319
320// Replication 1
321 located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(15,0,-10));
322 assert(located->GetName()=="Vari' Blocks");
323 located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(15,0,-17));
324 assert(located->GetName()=="Target 3");
325
326// Replication 2
327 located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(15,15,-10));
328 assert(located->GetName()=="Vari' Blocks");
329 located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(15,15,-18));
330 assert(located->GetName()=="Target 3");
331
332 return true;
333}
334
335
336//
337// Test Stepping
338//
339G4bool testG4Navigator2(G4VPhysicalVolume *pTopNode)
340{
341 MyNavigator myNav;
342 G4VPhysicalVolume *located;
343 G4double Step,physStep,safety;
344 G4ThreeVector xHat(1,0,0),yHat(0,1,0),zHat(0,0,1);
345 G4ThreeVector mxHat(-1,0,0),myHat(0,-1,0),mzHat(0,0,-1);
346
347 myNav.SetWorldVolume(pTopNode);
348
349//
350// Test location & Step computation
351//
352 located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(0,0,-10));
353 assert(located->GetName()=="World");
354 physStep=kInfinity;
355 Step=myNav.ComputeStep(G4ThreeVector(0,0,-10),mxHat,physStep,safety);
356 assert(ApproxEqual(Step,25));
357// assert(ApproxEqual(safety,5));
358 assert(safety>=0);
359
360 located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(0,0,-10));
361 assert(located->GetName()=="World");
362 physStep=kInfinity;
363 Step=myNav.ComputeStep(G4ThreeVector(0,0,-10),xHat,physStep,safety);
364 assert(ApproxEqual(Step,5));
365// assert(ApproxEqual(safety,5));
366 assert(safety>=0);
367 myNav.SetGeometricallyLimitedStep();
368 located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(5,0,-10),0,true);
369 assert(located->GetName()=="Vari' Blocks");
370
371 located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(0,0,-10));
372 assert(located->GetName()=="World");
373 physStep=kInfinity;
374 Step=myNav.ComputeStep(G4ThreeVector(0,0,-10),zHat,physStep,safety);
375 assert(ApproxEqual(Step,30));
376// assert(ApproxEqual(safety,5));
377 assert(safety>=0);
378
379 located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(0,0,-10));
380 assert(located->GetName()=="World");
381 physStep=kInfinity;
382 Step=myNav.ComputeStep(G4ThreeVector(0,0,-10),mzHat,physStep,safety);
383 assert(ApproxEqual(Step,10));
384// assert(ApproxEqual(safety,5));
385 assert(safety>=0);
386
387
388//
389// Test stepping through common boundaries
390//
391 located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(-7,7,-20));
392 assert(located->GetName()=="Target 1");
393 physStep=kInfinity;
394 Step=myNav.ComputeStep(G4ThreeVector(-7,7,-20),zHat,physStep,safety);
395 assert(ApproxEqual(Step,20));
396 assert(ApproxEqual(safety,0));
397 myNav.SetGeometricallyLimitedStep();
398 located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(-7,7,0));
399 assert(located->GetName()=="Target 4");
400 Step=myNav.ComputeStep(G4ThreeVector(-7,7,0),zHat,physStep,safety);
401 assert(ApproxEqual(Step,20));
402 assert(ApproxEqual(safety,0));
403 myNav.SetGeometricallyLimitedStep();
404 located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(-7,7,20));
405 assert(!located);
406
407//
408// Test mother limited Step
409//
410 located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(-25,0,10));
411 assert(located->GetName()=="World");
412 physStep=kInfinity;
413 Step=myNav.ComputeStep(G4ThreeVector(-25,0,10),xHat,physStep,safety);
414 assert(ApproxEqual(Step,50));
415 assert(ApproxEqual(safety,0));
416
417//
418// Test stepping through parameterised volumes
419//
420 located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(15,-25,-10),0,false);
421 assert(located->GetName()=="Target 3");
422 physStep=kInfinity;
423 Step=myNav.ComputeStep(G4ThreeVector(15,-25,-10),yHat,physStep,safety);
424 assert(ApproxEqual(Step,5));
425 assert(ApproxEqual(safety,0));
426 myNav.SetGeometricallyLimitedStep();
427 located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(15,-20,-10));
428 assert(located->GetName()=="Vari' Blocks");
429 Step=myNav.ComputeStep(G4ThreeVector(15,-20,-10),yHat,physStep,safety);
430 assert(ApproxEqual(Step,10));
431 assert(ApproxEqual(safety,0));
432 myNav.SetGeometricallyLimitedStep();
433 located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(15,-10,-10));
434 assert(located->GetName()=="Target 3");
435 Step=myNav.ComputeStep(G4ThreeVector(15,-10,-10),yHat,physStep,safety);
436 assert(ApproxEqual(Step,4));
437 assert(ApproxEqual(safety,0));
438 myNav.SetGeometricallyLimitedStep();
439 located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(15,-6,-10));
440 assert(located->GetName()=="Vari' Blocks");
441 Step=myNav.ComputeStep(G4ThreeVector(15,-6,-10),yHat,physStep,safety);
442 assert(ApproxEqual(Step,12));
443 assert(ApproxEqual(safety,0));
444 myNav.SetGeometricallyLimitedStep();
445 located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(15,6,-10));
446 assert(located->GetName()=="Target 3");
447 Step=myNav.ComputeStep(G4ThreeVector(15,6,-10),yHat,physStep,safety);
448 assert(ApproxEqual(Step,2));
449 assert(ApproxEqual(safety,0));
450 myNav.SetGeometricallyLimitedStep();
451 located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(15,8,-10));
452 assert(located->GetName()=="Vari' Blocks");
453 Step=myNav.ComputeStep(G4ThreeVector(15,8,-10),yHat,physStep,safety);
454 assert(ApproxEqual(Step,14));
455 assert(ApproxEqual(safety,0));
456 myNav.SetGeometricallyLimitedStep();
457 located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(15,22,-10));
458 assert(located->GetName()=="Target 3");
459 Step=myNav.ComputeStep(G4ThreeVector(15,22,-10),yHat,physStep,safety);
460 assert(ApproxEqual(Step,3));
461 assert(ApproxEqual(safety,0));
462 myNav.SetGeometricallyLimitedStep();
463 located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(15,25,-10));
464 assert(!located);
465
466 return true;
467}
468
469int main()
470{
471 G4VPhysicalVolume *myTopNode;
472 myTopNode=BuildGeometry(); // Build the geometry
473 G4GeometryManager::GetInstance()->CloseGeometry(false);
474 G4PathFinder* pPathFinder;
475
476 pPathFinder= setupPathFinder(myTopNode);
477
478 testG4PathFinder1(myTopNode);
479 // testG4Navigator2(myTopNode);
480// Repeat tests but with full voxels
481 G4GeometryManager::GetInstance()->OpenGeometry();
482 return 0;
483 // exit(1);
484
485 // --- Future 2nd test
486
487 G4GeometryManager::GetInstance()->CloseGeometry(true);
488 testG4PathFinder1(myTopNode);
489 // testG4Navigator2(myTopNode);
490
491 G4GeometryManager::GetInstance()->OpenGeometry();
492 return 0;
493}
Note: See TracBrowser for help on using the repository browser.