source: trunk/source/geometry/volumes/test/testG4NestedParameterised.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: 14.3 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: testG4NestedParameterised.cc,v 1.5 2006/06/29 18:58:36 gunter Exp $
28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
29//
30//
31//   Test the Navigation in geometry with parameterised volumes (which
32//    include rotations as well as translations).
33//   Locate & Step within simple boxlike geometry, both
34//   with and without voxels. Parameterised volumes are included.
35//   Started from testG4Navigator1.cc
36
37#include <assert.h>
38#include "G4ios.hh"
39#include "ApproxEqual.hh"
40
41// Global defs
42#include "globals.hh"
43
44#include "G4LogicalVolume.hh"
45#include "G4VPhysicalVolume.hh"
46#include "G4PVPlacement.hh"
47#include "G4PVParameterised.hh"
48#include "G4VPVParameterisation.hh"
49#include "G4VNestedParameterisation.hh"
50#include "G4Box.hh"
51
52#include "G4GeometryManager.hh"
53
54#include "G4RotationMatrix.hh"
55#include "G4ThreeVector.hh"
56
57// Sample Parameterisation
58class MoveNRotate : public G4VNestedParameterisation
59{
60 public:
61  MoveNRotate(G4double twistAngle)
62  { 
63    fTwistAngle= twistAngle;
64    fRotationVec= new G4RotationMatrix();
65  }
66
67  virtual ~MoveNRotate() { delete fRotationVec; }
68
69  G4double GetTwistAngle() { return fTwistAngle; }
70  void     SetTwistAngle(G4double newAngle ) { fTwistAngle= newAngle; }
71
72private:
73  void ComputeTransformation(const G4int n,
74                             G4VPhysicalVolume* pRep) const
75  {
76    pRep->SetTranslation(G4ThreeVector(0,n*100,0));
77    *fRotationVec = G4RotationMatrix();             // Unit matrix
78    fRotationVec->rotateZ( n * fTwistAngle );
79    pRep->SetRotation( fRotationVec );
80  }
81 
82  virtual void ComputeDimensions(G4Box &pBox,
83                                 const G4int,
84                                 const G4VPhysicalVolume*, 
85                                 const G4VTouchable* parentTouch
86                                 ) const
87  {
88    G4int no_parent= parentTouch->GetCopyNumber();
89    G4double half_len=0.0; 
90    if( no_parent == 0 ) { 
91      half_len= 10.0*mm; 
92    } else { 
93      half_len=  7.5*mm;
94    }
95
96    pBox.SetXHalfLength(half_len);
97    pBox.SetYHalfLength(half_len);
98    pBox.SetZHalfLength(half_len);
99  }
100 
101  G4int       GetNumberOfMaterials() const { return 0; } 
102  G4Material* GetMaterial(G4int)     const { return 0; } 
103
104  virtual void ComputeDimensions(G4Tubs &,
105                                 const G4int ,
106                                 const G4VPhysicalVolume*) const {}
107  virtual void ComputeDimensions(G4Trd &, 
108                                 const G4int,
109                                 const G4VPhysicalVolume*) const {}
110  virtual void ComputeDimensions(G4Cons &,
111                                 const G4int ,
112                                 const G4VPhysicalVolume*) const {}
113  virtual void ComputeDimensions(G4Trap &,
114                                 const G4int ,
115                                 const G4VPhysicalVolume*) const {}
116  virtual void ComputeDimensions(G4Hype &,
117                                 const G4int ,
118                                 const G4VPhysicalVolume*) const {}
119  virtual void ComputeDimensions(G4Orb &,
120                                 const G4int ,
121                                 const G4VPhysicalVolume*) const {}
122  virtual void ComputeDimensions(G4Sphere &,
123                                 const G4int ,
124                                 const G4VPhysicalVolume*) const {}
125  virtual void ComputeDimensions(G4Torus &,
126                                 const G4int ,
127                                 const G4VPhysicalVolume*) const {}
128  virtual void ComputeDimensions(G4Para &,
129                                 const G4int ,
130                                 const G4VPhysicalVolume*) const {}
131  virtual void ComputeDimensions(G4Polycone &,
132                                 const G4int ,
133                                 const G4VPhysicalVolume*) const {}
134  virtual void ComputeDimensions(G4Polyhedra &,
135                                 const G4int ,
136                                 const G4VPhysicalVolume*) const {}
137
138  // Mandatory method, required as reason for this class
139  virtual G4Material* ComputeMaterial(G4VPhysicalVolume *currentVol,
140                                      const G4int no_lev, 
141                                      const G4VTouchable *parentTouch);
142 private:
143    G4RotationMatrix *fRotationVec;
144    G4double fTwistAngle;
145};
146 
147G4Material* MoveNRotate::ComputeMaterial(G4VPhysicalVolume *currentVol,
148                                         const G4int no_lev, 
149                                         const G4VTouchable *parentTouch)
150{
151    // Get the information about the parent volume
152    G4int no_parent= parentTouch->GetReplicaNumber(); 
153
154    G4int no_total= no_parent + no_lev; 
155    if (no_total == 0) no_total += 32; 
156
157    G4Material *material= 0;
158    // Can add a material here ... that depends on no_parent & no_lev
159
160    G4LogicalVolume* currentLogVol= currentVol->GetLogicalVolume(); 
161
162    currentLogVol->SetMaterial( material ); 
163
164    return material;
165}
166
167
168G4double    angle1= 15.0*pi/180.;
169MoveNRotate myParam(angle1);
170
171
172// Build simple geometry:
173// 4 small cubes (G4Boxes) are positioned inside a larger cuboid
174G4VPhysicalVolume* BuildGeometry()
175{
176
177    // The world volume
178    //
179    G4Box *myBigBox= new G4Box ("Big Cube", 500, 500, 500);
180
181    G4LogicalVolume *worldLog=new G4LogicalVolume(myBigBox,0,
182                                                  "World",0,0,0);
183                                // Logical with no material,field,
184                                // sensitive detector or user limits
185   
186    G4PVPlacement *worldPhys=new G4PVPlacement(0,G4ThreeVector(0,0,0),
187                                               "World",worldLog,
188                                               0,false,0);
189                                // Note: no mother pointer set
190
191
192    // A set of boxes
193    G4Box *myBox=new G4Box("cube",10,10,10);  // Dimensions can change
194    G4LogicalVolume *boxLog=new G4LogicalVolume(myBox,0,
195                                                "Rotating Box",0,0,0);
196
197//  G4PVParameterised *paramP=
198                              new G4PVParameterised("Rotating Blocks",
199                                                    boxLog,
200                                                    worldPhys, //OR worldLog,
201                                                    kYAxis,
202                                                    3,
203                                                    &myParam);
204    // Copies 0, 1 & 2 will exist   
205
206    return worldPhys;
207}
208
209//
210// Test LocateGlobalPointAndSetup
211//
212G4bool testG4Navigator1(G4VPhysicalVolume *pTopNode)
213{
214    MyNavigator myNav;
215    G4VPhysicalVolume *located;
216    myNav.SetWorldVolume(pTopNode);
217
218    assert(!myNav.LocateGlobalPointAndSetup(G4ThreeVector(kInfinity,0,0),0, false));
219    located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(100,100,100),0,false);
220    assert(located->GetName()=="World");
221
222    assert(!myNav.LocateGlobalPointAndSetup(G4ThreeVector(kInfinity,0,0)));
223
224//
225    located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(0,-5,-5),0,false);
226    assert(located->GetName()=="Rotating Blocks");
227    assert(located->GetCopyNo()== 0);
228    assert(ApproxEqual(myNav.CurrentLocalCoordinate(),G4ThreeVector(0,-5,-5)));
229    G4cout << " Local coords = " << myNav.CurrentLocalCoordinate() << G4endl;
230
231    located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(0,100,5));
232    assert(located->GetName()=="Rotating Blocks");
233    assert(located->GetCopyNo()== 1);
234    G4cout << " Local coords = " << myNav.CurrentLocalCoordinate() << G4endl;
235//    assert(ApproxEqual(myNav.CurrentLocalCoordinate(),
236//                       G4ThreeVector(0,0,10)));
237   
238// Check that outside point causes stack to unwind
239    assert(!myNav.LocateGlobalPointAndSetup(G4ThreeVector(kInfinity,0,0)));
240
241// Check parameterised volumes
242
243// Replication 0
244    located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(0,5,5));
245    assert(located->GetName()=="Rotating Blocks");
246    assert(located->GetCopyNo()== 0);
247    located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(0,15,15));
248    assert(located->GetName()=="World");
249
250// Replication 1
251    located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(0,105,5));
252    assert(located->GetName()=="Rotating Blocks");
253    assert(located->GetCopyNo()== 1);
254    located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(0,0,-17));
255    assert(located->GetName()=="World");
256
257// Replication 2
258    located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(0,205,5));
259    assert(located->GetName()=="Rotating Blocks");
260    assert(located->GetCopyNo()== 2);
261    located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(15,15,-18));
262    assert(located->GetName()=="World");
263
264    return true;
265}
266
267
268//
269// Test Stepping
270//
271G4bool testG4Navigator2(G4VPhysicalVolume *pTopNode)
272{
273    MyNavigator myNav;
274    G4VPhysicalVolume *located;
275    G4double Step,physStep,safety;
276    G4ThreeVector xHat(1,0,0),yHat(0,1,0),zHat(0,0,1);
277    G4ThreeVector mxHat(-1,0,0),myHat(0,-1,0),mzHat(0,0,-1);
278   
279    myNav.SetWorldVolume(pTopNode);
280 
281//
282// Test location & Step computation
283// 
284    G4ThreeVector  StartPoint(-50,0,-5);
285    located=myNav.LocateGlobalPointAndSetup( StartPoint ); 
286    assert(located->GetName()=="World");
287    physStep=kInfinity;
288    Step=myNav.ComputeStep( StartPoint, mxHat,physStep,safety);  // -x dir
289    assert(ApproxEqual(Step,450));
290    // assert(ApproxEqual(safety,40));
291    // assert(safety>=0);
292
293    StartPoint= G4ThreeVector(-15,0,-5);
294    located=myNav.LocateGlobalPointAndSetup( StartPoint ); 
295    assert(located->GetName()=="World");
296    physStep=kInfinity;
297    Step=myNav.ComputeStep( StartPoint,xHat,physStep,safety); // +x dir
298    assert(ApproxEqual(Step,5));
299//    assert(ApproxEqual(safety,5));
300    assert(safety>=0);
301    myNav.SetGeometricallyLimitedStep();
302    G4ThreeVector EndPoint = StartPoint + Step * xHat;
303    located=myNav.LocateGlobalPointAndSetup(EndPoint,0,true);
304    assert(located->GetName()=="Rotating Blocks");
305
306    located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(0,0,-40));
307    assert(located->GetName()=="World");
308    physStep=kInfinity;
309    Step=myNav.ComputeStep(G4ThreeVector(0,0,-40),zHat,physStep,safety);
310    assert(ApproxEqual(Step,30));
311//    assert(ApproxEqual(safety,5));
312    assert(safety>=0);
313
314    located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(0,0, 40));
315    assert(located->GetName()=="World");
316    physStep=kInfinity;
317    Step=myNav.ComputeStep(G4ThreeVector(0,0,40),mzHat,physStep,safety);
318    assert(ApproxEqual(Step,30));
319//    assert(ApproxEqual(safety,5));
320    assert(safety>=0);
321
322
323//
324// Test moving through series of volumes
325//
326    StartPoint= G4ThreeVector(0,-20,0);
327    located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(0,-20,0));
328    assert(located->GetName()=="World");
329   
330    // Replication 0 block
331    //
332    physStep=kInfinity;
333    Step=myNav.ComputeStep(G4ThreeVector(0,-20,0),yHat,physStep,safety);
334    assert(ApproxEqual(Step,10));
335    EndPoint= StartPoint + Step * yHat;   //  Should be  0, -10, 0
336    assert(ApproxEqual( 0, (EndPoint-G4ThreeVector(0,-10,0)).mag()) );
337    // assert(ApproxEqual(safety,0));
338   
339    myNav.SetGeometricallyLimitedStep();
340    located=myNav.LocateGlobalPointAndSetup(EndPoint) ;
341    assert(located->GetName()=="Rotating Blocks");
342    Step=myNav.ComputeStep(EndPoint,yHat,physStep,safety);
343    assert(ApproxEqual(Step,20));
344    assert(ApproxEqual(safety,0));
345    myNav.SetGeometricallyLimitedStep();
346    EndPoint += Step * yHat;   //  Should be  0, +10, 0
347    located=myNav.LocateGlobalPointAndSetup( EndPoint );
348    assert(located->GetName()=="World");
349
350    // Replication 1 block
351    //
352    StartPoint= EndPoint;
353    physStep=kInfinity;
354    Step=myNav.ComputeStep(StartPoint,yHat,physStep,safety);
355    assert(ApproxEqual(Step,90.-10./std::cos(angle1)));
356    EndPoint= StartPoint + Step * yHat;   //  Should near  0, 90, 0
357    assert(safety<=Step);
358    myNav.SetGeometricallyLimitedStep();
359    located=myNav.LocateGlobalPointAndSetup(EndPoint) ;
360    assert(located->GetName()=="Rotating Blocks");
361   
362    StartPoint= EndPoint;
363    physStep=kInfinity;
364    Step=myNav.ComputeStep(StartPoint,yHat,physStep,safety);
365    assert(ApproxEqual(Step,20./std::cos(angle1)));
366    assert(ApproxEqual(safety,0));
367    myNav.SetGeometricallyLimitedStep();
368    EndPoint += Step * yHat;   //  Should be near 0, 110, 0
369    located=myNav.LocateGlobalPointAndSetup( EndPoint );
370    assert(located->GetName()=="World");
371
372    // Replication 2 block
373    //
374    StartPoint= EndPoint;
375    physStep=kInfinity;
376    Step=myNav.ComputeStep(StartPoint,yHat,physStep,safety);
377    assert(ApproxEqual(Step,100.-10.*(1./std::cos(angle1)+1./std::cos(2.*angle1))));
378    EndPoint= StartPoint + Step * yHat;   //  Should near  0, 190, 0
379    assert(safety<=Step);
380    myNav.SetGeometricallyLimitedStep();
381    located=myNav.LocateGlobalPointAndSetup(EndPoint);
382    assert(located->GetName()=="Rotating Blocks");
383   
384    StartPoint= EndPoint;
385    physStep=kInfinity;
386    Step=myNav.ComputeStep(StartPoint,yHat,physStep,safety);
387    assert(ApproxEqual(Step,20./std::cos(2.*angle1)));
388    assert(ApproxEqual(safety,0));
389    myNav.SetGeometricallyLimitedStep();
390    EndPoint += Step * yHat;   //  Should be near 0, 110, 0
391    located=myNav.LocateGlobalPointAndSetup( EndPoint );
392    assert(located->GetName()=="World");
393
394    // Edge of the world
395    //
396    StartPoint= EndPoint;
397    physStep=kInfinity;
398    Step=myNav.ComputeStep(StartPoint,yHat,physStep,safety);
399    assert(ApproxEqual(Step, 300. - 10./std::cos(2.*angle1) ));
400    assert(ApproxEqual(safety,0));
401    myNav.SetGeometricallyLimitedStep();
402    EndPoint += Step * yHat;   //  Should be near 0, 110, 0
403    located=myNav.LocateGlobalPointAndSetup( EndPoint );
404    assert(!located);
405
406
407    return true;
408}
409
410int main()
411{
412    G4VPhysicalVolume *myTopNode;
413    myTopNode=BuildGeometry();  // Build the geometry
414    G4GeometryManager::GetInstance()->CloseGeometry(false);
415    testG4Navigator1(myTopNode);
416    testG4Navigator2(myTopNode);
417// Repeat tests but with full voxels
418    G4GeometryManager::GetInstance()->OpenGeometry();
419    G4GeometryManager::GetInstance()->CloseGeometry(true);
420    testG4Navigator1(myTopNode);
421    testG4Navigator2(myTopNode);
422
423    G4GeometryManager::GetInstance()->OpenGeometry();
424    return 0;
425}
426
427
428
429
Note: See TracBrowser for help on using the repository browser.