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