source: trunk/source/geometry/navigation/test/testG4PathFinder.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: 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.