source: trunk/source/geometry/volumes/test/testG4NestedParameterised.cc@ 1350

Last change on this file since 1350 was 1347, checked in by garnier, 15 years ago

geant4 tag 9.4
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-ref-00 $
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.