source: trunk/source/geometry/volumes/test/testG4ParameterisedMaterial.cc@ 1349

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

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