source: trunk/source/geometry/navigation/test/testG4NavigatorSafety.cc

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

update to last version 4.9.4
File size: 9.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: testG4NavigatorSafety.cc,v 1.6 2010/11/09 10:48:27 gcosmo Exp $
28// GEANT4 tag $Name: geant4-09-04-ref-00 $
29//
30//
31// Create a tubular "calorimeter". Generate random points along x, y & z
32// axes, printing location, steps & safeties. Compare results of standard
33// voxel safety calculations with "exact safety" computed values.
34//
35// Optional arguments:
36// - Number of points to generate [Default: 10000]
37// - Initial random seed modulo 256 [Default: CLHEP default]
38
39#include <assert.h>
40#include "G4ios.hh"
41#include <stdlib.h>
42#include <vector>
43
44#include "globals.hh"
45
46#include "G4Timer.hh"
47#include "ApproxEqual.hh"
48
49#include "G4Navigator.hh"
50
51#include "G4LogicalVolume.hh"
52#include "G4VPhysicalVolume.hh"
53#include "G4PVPlacement.hh"
54#include "G4Box.hh"
55#include "G4Tubs.hh"
56
57#include "G4GeometryManager.hh"
58
59#include "G4RotationMatrix.hh"
60#include "G4ThreeVector.hh"
61#include "Randomize.hh"
62
63// Parameters for building a tubular calorimeter:
64// an array of interlocking complete tubes inside a box
65//
66static const G4double kTubeHalfHeight = 10*mm;
67static const G4double kTubeRadius = 5*mm;
68static const G4double kTubeNoRow = 10;
69static const G4double kTubeNoColumn = 11; // Odd for symmetrical array
70
71static const G4double kBoxDx=kTubeNoRow*kTubeRadius;
72static const G4double yDelta=2.0*kTubeRadius*std::sin(pi/3.0);
73static const G4double kBoxDy=(kTubeNoColumn-1)*yDelta*0.5+kTubeRadius;
74static const G4double kBoxDz=kTubeHalfHeight;
75static const G4double kWorldhxsize = kBoxDx+0.1*mm;
76static const G4double kWorldhysize = kBoxDy+0.1*mm;
77static const G4double kWorldhzsize = kBoxDz+0.1*mm;
78
79G4bool compare = false;
80std::vector<G4ThreeVector> kPoints;
81std::vector<std::pair<G4double, G4double> > kSafeties;
82
83G4VPhysicalVolume* BuildGeometry()
84{
85 G4double bigXStart=-(kTubeNoRow-1)*kTubeRadius;
86 G4double smallXStart=bigXStart+kTubeRadius;
87
88 G4double bigYStart=-(kTubeNoColumn-1)*yDelta*0.5;
89 G4double smallYStart=bigYStart+yDelta;
90
91 G4int row,column;
92
93 // Solids ==============================
94 G4Box *worldBox = new G4Box ("World Box",kWorldhxsize,kWorldhysize,kWorldhzsize);
95 // World box
96
97 G4Box *calBox = new G4Box ("Cal Box",kBoxDx,kBoxDy,kBoxDz);
98 G4Tubs *calTube = new G4Tubs("Cal Tube",0,kTubeRadius,kTubeHalfHeight,0,360);
99
100 // Logical Volumes ------------------------------
101 G4LogicalVolume *myWorldLog=
102 new G4LogicalVolume(worldBox,0,"World",0,0,0);
103 // Logical with no material,field, sensitive detector or user limits
104
105 G4PVPlacement *myWorldPhys=
106 new G4PVPlacement(0,G4ThreeVector(0,0,0),"World",myWorldLog,0,false,0);
107 // World physical volume
108
109 G4LogicalVolume *myDetectorLog=
110 new G4LogicalVolume(calBox,0,"DetectorLog",0,0,0);
111 // Logical with no material,field, sensitive detector or user limits
112
113 G4PVPlacement *myDetectorPhys=
114 new G4PVPlacement(0,G4ThreeVector(0,0,0),"DetectorPhys",
115 myDetectorLog,myWorldPhys,false,0);
116 // Detector physical volume placed in the world
117
118 G4LogicalVolume *calTubLog=
119 new G4LogicalVolume(calTube,0,"Cal Crystal",0,0,0);
120
121 G4String tname("Target");
122 G4int copyNo=0;
123 for (column=0;column<kTubeNoColumn;column+=2)
124 {
125 for (row=0;row<kTubeNoRow;row++)
126 {
127// G4PVPlacement *calPhys=
128 new G4PVPlacement(0,G4ThreeVector(bigXStart+row*kTubeRadius*2.0,
129 bigYStart+column*yDelta,0),
130 tname,calTubLog,myDetectorPhys,false,copyNo++);
131 }
132 }
133 for (column=0;column<kTubeNoColumn-1;column+=2)
134 {
135 for (row=0;row<kTubeNoRow-1;row++)
136 {
137// G4PVPlacement *calPhys=
138 new G4PVPlacement(0,G4ThreeVector(smallXStart+row*kTubeRadius*2.0,
139 smallYStart+column*yDelta),
140 tname,calTubLog,myDetectorPhys,false,copyNo++);
141 }
142 }
143 return myWorldPhys;
144}
145
146void generatePoints(G4int n)
147{
148 for (int i=0; i<n; i++)
149 {
150 G4ThreeVector p(CLHEP::RandFlat::shoot(-kWorldhxsize,kWorldhxsize),
151 CLHEP::RandFlat::shoot(-kWorldhysize,kWorldhysize),
152 CLHEP::RandFlat::shoot(-kWorldhzsize,kWorldhzsize));
153 kPoints.push_back(p);
154 }
155}
156
157void computeApproxSafeties(G4VPhysicalVolume *pTopNode)
158{
159 MyNavigator myNav;
160 myNav.SetWorldVolume(pTopNode);
161 std::vector<G4ThreeVector>::const_iterator pos;
162 for (pos=kPoints.begin(); pos!=kPoints.end(); pos++)
163 {
164 myNav.LocateGlobalPointAndSetup(*pos);
165 G4double safety = myNav.ComputeSafety(*pos);
166 std::pair<G4double,G4double> s(safety,0.);
167 kSafeties.push_back(s);
168 }
169}
170
171void computeExactSafeties(G4VPhysicalVolume *pTopNode)
172{
173 G4int i=0;
174 G4VPhysicalVolume *pPhysVol=0;
175
176 MyNavigator myNav;
177 myNav.SetWorldVolume(pTopNode);
178
179 // Enable use of Best Safety Estimate -- ie as exact as solids allow
180 myNav.UseBestSafety( true );
181
182 myNav.SetVerboseLevel( 1 );
183 myNav.CheckMode( true );
184
185 G4ThreeVector center( 0., 0., 0. );
186 G4cout << " Trial point= " << center << G4endl;
187 pPhysVol= myNav.LocateGlobalPointAndSetup( center );
188 // G4double saf= myNav.ComputeSafety( center );
189
190 std::vector<G4ThreeVector>::const_iterator pos;
191 for (pos=kPoints.begin(); pos!=kPoints.end(); pos++)
192 {
193 G4cout << " Trial point= " << *pos << G4endl;
194 // Relocate point
195 pPhysVol= myNav.LocateGlobalPointAndSetup(*pos);
196
197 G4cout << G4endl;
198 G4cout << " ============================================== " << G4endl;
199 G4cout << " Calculating 'exact' safety for " << *pos << G4endl;
200 G4double safety = myNav.ComputeSafety(*pos);
201 kSafeties[i].second = safety;
202
203 G4cout << " Old safety = " << kSafeties[i].first << G4endl;
204 G4cout << " Exact safety = " << safety << G4endl;
205
206 i++;
207 }
208 compare = true;
209}
210
211void compareSafeties()
212{
213 G4int n = kPoints.size();
214 if (!compare)
215 {
216 G4cout << "Printing out non-zero safety values computed ..." << G4endl;
217 for (G4int i=0; i<n; i++)
218 {
219 G4double safety = kSafeties[i].first;
220 if (safety)
221 {
222 G4cout << i << ". - Point: " << kPoints[i]
223 << " - Safety: " << safety << G4endl;
224 }
225 }
226 }
227 else
228 {
229 G4int diffs=0;
230 G4cout << "Printing out cases of different safety values ..." << G4endl;
231 for (G4int i=0; i<n; i++)
232 {
233 if (kSafeties[i].first != kSafeties[i].second)
234 {
235 G4cout << i << ". - Point: " << kPoints[i]
236 << " - Approx safety: " << kSafeties[i].first
237 << " - Exact safety: " << kSafeties[i].second << G4endl;
238 diffs++;
239 }
240 }
241 G4cout << "Total number of differences: " << diffs << G4endl;
242 }
243}
244
245int main(int argc, char* argv[])
246{
247 G4Timer timer;
248 G4int iter=10000;
249
250 if (argc==2)
251 {
252 G4int num = atoi(argv[1]);
253 if (num>=0) { iter = num; }
254 else { G4cout << ">>> Invalid number of iterations in input!" << G4endl
255 << " Sticking to: " << iter << " ..." << G4endl; }
256 }
257 else if (argc==3)
258 {
259 G4int num = atoi(argv[1]);
260 if (num>=0) { iter = num; }
261 else { G4cout << ">>> Invalid number of iterations in input!" << G4endl
262 << " Sticking to: " << iter << " ..." << G4endl; }
263 G4long seed = atoi(argv[2]);
264 if (seed>0) { CLHEP::HepRandom::setTheSeed(seed); }
265 else { G4cout << ">>> Invalid negative random seed in input!" << G4endl
266 << " Sticking to default: "
267 << CLHEP::HepRandom::getTheSeed() << " ..." << G4endl; }
268 }
269
270 G4VPhysicalVolume *myTopNode;
271 myTopNode=BuildGeometry(); // Build the geometry
272 // Do not close the geometry --> the voxels will limit safety
273
274 timer.Start();
275
276 generatePoints(iter);
277
278 computeApproxSafeties(myTopNode);
279
280 G4GeometryManager::GetInstance()->CloseGeometry(); // Voxelise the geometry
281 computeExactSafeties(myTopNode);
282
283 // Check
284 compareSafeties();
285
286 timer.Stop();
287
288 G4cout << timer << G4endl;
289
290 G4GeometryManager::GetInstance()->OpenGeometry();
291
292 return 0;
293}
Note: See TracBrowser for help on using the repository browser.