source: trunk/source/geometry/volumes/test/testG4ParameterisedMaterial.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.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-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 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.