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: testG4PVDivision.cc,v 1.5 2009/05/14 14:19:32 ivana Exp $ |
---|
28 | // GEANT4 tag $Name: geant4-09-04-ref-00 $ |
---|
29 | // |
---|
30 | // test for G4PVDivision classes |
---|
31 | // |
---|
32 | // 26.05.03 - P.Arce Initial version |
---|
33 | // ******************************************************************** |
---|
34 | |
---|
35 | #include <assert.h> |
---|
36 | #include <fstream> |
---|
37 | #include <vector> |
---|
38 | |
---|
39 | #include "G4ios.hh" |
---|
40 | #include "globals.hh" |
---|
41 | |
---|
42 | #include "G4Navigator.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 | #include "G4Trd.hh" |
---|
51 | #include "G4Trap.hh" |
---|
52 | #include "G4Tubs.hh" |
---|
53 | #include "G4Cons.hh" |
---|
54 | #include "G4Para.hh" |
---|
55 | #include "G4Polycone.hh" |
---|
56 | #include "G4Polyhedra.hh" |
---|
57 | |
---|
58 | #include "G4GeometryManager.hh" |
---|
59 | |
---|
60 | #include "G4RotationMatrix.hh" |
---|
61 | #include "G4ThreeVector.hh" |
---|
62 | |
---|
63 | #include "Randomize.hh" |
---|
64 | #include "G4PVReplica.hh" |
---|
65 | |
---|
66 | #include "G4PVDivision.hh" |
---|
67 | |
---|
68 | G4String theSolidTypeStr; |
---|
69 | G4String thePVTypeStr; |
---|
70 | G4double theWorldDim = 1*m; |
---|
71 | G4int numberOfPoints = 1000; |
---|
72 | G4int theNReplicas; |
---|
73 | G4double theWidth; |
---|
74 | G4double theOffset; |
---|
75 | G4VSolid* theParentSolid; |
---|
76 | std::vector<G4LogicalVolume*> theParentLogs; |
---|
77 | std::vector<G4VPhysicalVolume*> theParentPhyss; |
---|
78 | std::vector<G4VSolid*> theChildSolids; |
---|
79 | std::vector<G4LogicalVolume*> theChildLogs; |
---|
80 | std::vector<EAxis> theAxis; |
---|
81 | std::vector<G4double> theWidths; |
---|
82 | std::vector<G4double> theExtraPars; |
---|
83 | G4int theDivType; |
---|
84 | |
---|
85 | enum SolidType{g4box, g4trd, g4tube, g4tubs, g4cone, g4cons, g4polycone, g4polyhedra }; |
---|
86 | enum PVType{pvDivision, pvReplica, pvPlacement }; |
---|
87 | |
---|
88 | //-------------------------------------------------------------------------- |
---|
89 | void initialize(); |
---|
90 | void calculateParentSolid( SolidType solType ); |
---|
91 | void calculateChildSolids( SolidType solType, G4VSolid* pSolid ); |
---|
92 | void calculateAxis( SolidType solType ); |
---|
93 | void buildOutputName( SolidType& soltype, PVType pvtype ); |
---|
94 | void starttest( const std::vector<G4String>& vsarg ); |
---|
95 | G4int checkNumberOfArguments( const G4String& st, G4int narg ); |
---|
96 | PVType getPVType( const G4String& pvt ); |
---|
97 | SolidType getSolidType( const G4String& st ); |
---|
98 | void generateRandomPoints(); |
---|
99 | void generateScanPointsForBox(); |
---|
100 | void generateScanPointsForTube(); |
---|
101 | void generateScanPointsForTrd(); |
---|
102 | void generateScanPointsForPolycone(); |
---|
103 | G4VPhysicalVolume* BuildGeometry( SolidType solType, PVType pvType ); |
---|
104 | G4bool testG4Navigator1(G4VPhysicalVolume *pTopNode); |
---|
105 | G4bool testG4Navigator2(G4VPhysicalVolume *pTopNode); |
---|
106 | // World geometry is a box 1 X 1 X 3. |
---|
107 | // User select the type of solid |
---|
108 | // Inside it three solids of the chosen type are placed along Z |
---|
109 | // Each of these solids is divided along a different axis |
---|
110 | // |
---|
111 | // Then a set of points is generated and it is tested in which division copy they are and which is the DistanceToOut in the three directions (X,Y,Z) |
---|
112 | |
---|
113 | //-------------------------------------------------------------------------- |
---|
114 | int main( G4int argc, char** argv ) |
---|
115 | { |
---|
116 | // first argument is type of divisioning (repli/divis), second is type of solid |
---|
117 | std::vector<G4String> vsarg; |
---|
118 | for( G4int jj = 0; jj < argc; jj ++ ) { |
---|
119 | vsarg.push_back( G4String(argv[jj] ) ); |
---|
120 | } |
---|
121 | |
---|
122 | starttest( vsarg ); |
---|
123 | } |
---|
124 | |
---|
125 | //-------------------------------------------------------------------------- |
---|
126 | void starttest( const std::vector<G4String>& vsarg ) |
---|
127 | { |
---|
128 | G4int narg = vsarg.size(); |
---|
129 | if( narg == 1 ) { |
---|
130 | thePVTypeStr = "divis"; |
---|
131 | theSolidTypeStr = "box"; |
---|
132 | } else if( narg == 2 ){ |
---|
133 | // wrong number |
---|
134 | checkNumberOfArguments( " ", narg ); |
---|
135 | } else { |
---|
136 | G4int divTypeSet = checkNumberOfArguments( vsarg[2], narg ); |
---|
137 | if( divTypeSet ) { |
---|
138 | theDivType = atoi( vsarg[divTypeSet] ); |
---|
139 | } else { |
---|
140 | theDivType = 0; |
---|
141 | } |
---|
142 | thePVTypeStr = G4String(vsarg[1]); |
---|
143 | theSolidTypeStr = G4String(vsarg[2]); |
---|
144 | } |
---|
145 | if( narg > 3 ){ |
---|
146 | for( G4int ii = 3; ii < narg; ii++ ){ |
---|
147 | theExtraPars.push_back( atof( vsarg[ii] ) ); |
---|
148 | } |
---|
149 | } |
---|
150 | PVType pvtype = getPVType( thePVTypeStr ); |
---|
151 | SolidType soltype = getSolidType( theSolidTypeStr ); |
---|
152 | |
---|
153 | // PVType pvtype = pvDivision; |
---|
154 | initialize(); |
---|
155 | if( soltype == g4tube || soltype == g4tubs ) { |
---|
156 | // generateRandomPoints(); |
---|
157 | generateScanPointsForTube(); |
---|
158 | }else if( soltype == g4cone || soltype == g4cons ) { |
---|
159 | // generateRandomPoints(); |
---|
160 | generateScanPointsForTube(); |
---|
161 | } else if( soltype == g4box ) { |
---|
162 | // generateRandomPoints(); |
---|
163 | generateScanPointsForBox(); |
---|
164 | } else if( soltype == g4trd ) { |
---|
165 | // generateRandomPoints(); |
---|
166 | generateScanPointsForTrd(); |
---|
167 | } else if( soltype == g4polycone ) { |
---|
168 | // generateRandomPoints(); |
---|
169 | generateScanPointsForPolycone(); |
---|
170 | } else { |
---|
171 | generateRandomPoints(); |
---|
172 | } |
---|
173 | |
---|
174 | G4VPhysicalVolume *myTopNode; |
---|
175 | myTopNode=BuildGeometry( soltype, pvtype ); // Build the geometry |
---|
176 | G4GeometryManager::GetInstance()->CloseGeometry(false); |
---|
177 | testG4Navigator1(myTopNode); |
---|
178 | testG4Navigator2(myTopNode); |
---|
179 | // Repeat tests but with full voxels |
---|
180 | G4GeometryManager::GetInstance()->OpenGeometry(); |
---|
181 | G4GeometryManager::GetInstance()->CloseGeometry(true); |
---|
182 | testG4Navigator1(myTopNode); |
---|
183 | testG4Navigator2(myTopNode); |
---|
184 | |
---|
185 | G4cout << " theParentSolid " << G4endl; |
---|
186 | G4cout << *theParentSolid << G4endl; |
---|
187 | for( size_t ii = 0; ii < theChildSolids.size(); ii++) { |
---|
188 | G4cout << " theChildSolid after tracking " << " "<< ii << G4endl; |
---|
189 | G4cout << *theChildSolids[ii] << G4endl; |
---|
190 | } |
---|
191 | |
---|
192 | G4GeometryManager::GetInstance()->OpenGeometry(); |
---|
193 | } |
---|
194 | |
---|
195 | //-------------------------------------------------------------------------- |
---|
196 | void generateRandomPoints() |
---|
197 | { |
---|
198 | std::ofstream fout("points.lis"); |
---|
199 | G4double posX, posY, posZ; |
---|
200 | |
---|
201 | for( G4int ii = 0; ii < numberOfPoints; ii++ ) { |
---|
202 | posX = CLHEP::RandFlat::shoot(-theWorldDim, theWorldDim); |
---|
203 | posY = CLHEP::RandFlat::shoot(-theWorldDim, theWorldDim); |
---|
204 | posZ = CLHEP::RandFlat::shoot(-3*theWorldDim, 3*theWorldDim); |
---|
205 | fout << posX << " " << posY << " " << posZ << G4endl; |
---|
206 | } |
---|
207 | } |
---|
208 | |
---|
209 | //-------------------------------------------------------------------------- |
---|
210 | void generateScanPointsForBox() |
---|
211 | { |
---|
212 | std::ofstream fout("points.lis"); |
---|
213 | G4int ii; |
---|
214 | |
---|
215 | G4int nPointsPerReplica = 2; |
---|
216 | G4int numberOfPoints = theNReplicas*nPointsPerReplica; |
---|
217 | // For division along X |
---|
218 | G4ThreeVector centre(0.,0.,-2*theWorldDim); |
---|
219 | for( ii = 0; ii < numberOfPoints; ii++ ) { |
---|
220 | // any Z, any Y |
---|
221 | G4ThreeVector pR( 0., theWorldDim/100., theWorldDim/100. ); |
---|
222 | G4double X = -theWorldDim + (ii+0.001) * 2*theWorldDim/numberOfPoints; |
---|
223 | pR.setX( X ); |
---|
224 | pR += centre; |
---|
225 | fout << pR.x() << " " << pR.y() << " " << pR.z() << G4endl; |
---|
226 | } |
---|
227 | |
---|
228 | // For division along Y |
---|
229 | centre = G4ThreeVector(0.,0.,0.); |
---|
230 | for( ii = 0; ii < numberOfPoints; ii++ ) { |
---|
231 | // any X, any Z |
---|
232 | G4ThreeVector pR( theWorldDim/100., 0., theWorldDim/100. ); |
---|
233 | G4double Y = -theWorldDim + (ii+0.001) * 2*theWorldDim/numberOfPoints; |
---|
234 | pR.setY( Y ); |
---|
235 | pR += centre; |
---|
236 | fout << pR.x() << " " << pR.y() << " " << pR.z() << G4endl; |
---|
237 | } |
---|
238 | |
---|
239 | // For division along Z |
---|
240 | centre = G4ThreeVector(0.,0.,2*theWorldDim); |
---|
241 | for( ii = 0; ii < numberOfPoints; ii++ ) { |
---|
242 | // any X, any Y |
---|
243 | G4ThreeVector pR( theWorldDim/100., 0., theWorldDim/100. ); |
---|
244 | G4double Z = -theWorldDim + (ii+0.001) * 2*theWorldDim/numberOfPoints; |
---|
245 | pR.setZ( Z ); |
---|
246 | pR += centre; |
---|
247 | fout << pR.x() << " " << pR.y() << " " << pR.z() << G4endl; |
---|
248 | } |
---|
249 | } |
---|
250 | |
---|
251 | //-------------------------------------------------------------------------- |
---|
252 | void generateScanPointsForTrd() |
---|
253 | { |
---|
254 | std::ofstream fout("points.lis"); |
---|
255 | G4int ii; |
---|
256 | |
---|
257 | G4int nPointsPerReplica = 2; |
---|
258 | G4int numberOfPoints = theNReplicas*nPointsPerReplica; |
---|
259 | // For division along X |
---|
260 | G4ThreeVector centre(0.,0.,-2*theWorldDim); |
---|
261 | for( ii = 0; ii < numberOfPoints; ii++ ) { |
---|
262 | // any Z, any Y |
---|
263 | G4ThreeVector pR( 0., theWorldDim/100., theWorldDim/100. ); |
---|
264 | G4double X = -theWorldDim*0.75 + (ii+0.001) * 1.5*theWorldDim/numberOfPoints; |
---|
265 | pR.setX( X ); |
---|
266 | pR += centre; |
---|
267 | fout << pR.x() << " " << pR.y() << " " << pR.z() << G4endl; |
---|
268 | } |
---|
269 | |
---|
270 | // For division along Y |
---|
271 | centre = G4ThreeVector(0.,0.,0.); |
---|
272 | for( ii = 0; ii < numberOfPoints; ii++ ) { |
---|
273 | // any X, any Z |
---|
274 | G4ThreeVector pR( theWorldDim/100., 0., theWorldDim/100. ); |
---|
275 | G4double Y = -theWorldDim*0.75 + (ii+0.001) * 1.5*theWorldDim/numberOfPoints; |
---|
276 | pR.setY( Y ); |
---|
277 | pR += centre; |
---|
278 | fout << pR.x() << " " << pR.y() << " " << pR.z() << G4endl; |
---|
279 | } |
---|
280 | |
---|
281 | // For division along Z |
---|
282 | centre = G4ThreeVector(0.,0.,2*theWorldDim); |
---|
283 | for( ii = 0; ii < numberOfPoints; ii++ ) { |
---|
284 | // any X, any Y |
---|
285 | G4ThreeVector pR( theWorldDim/100., 0., theWorldDim/100. ); |
---|
286 | G4double Z = -theWorldDim + (ii+0.001) * 2*theWorldDim/numberOfPoints; |
---|
287 | pR.setZ( Z ); |
---|
288 | pR += centre; |
---|
289 | fout << pR.x() << " " << pR.y() << " " << pR.z() << G4endl; |
---|
290 | } |
---|
291 | } |
---|
292 | |
---|
293 | //-------------------------------------------------------------------------- |
---|
294 | void generateScanPointsForTube() |
---|
295 | { |
---|
296 | std::ofstream fout("points.lis"); |
---|
297 | G4int ii; |
---|
298 | |
---|
299 | G4int nPointsPerReplica = 2; |
---|
300 | G4int numberOfPoints = theNReplicas*nPointsPerReplica; |
---|
301 | G4cout << " numberOfPoints " << numberOfPoints << G4endl; |
---|
302 | // For division along R |
---|
303 | G4ThreeVector centre(0.,0.,-2*theWorldDim); |
---|
304 | for( ii = 0; ii < numberOfPoints; ii++ ) { |
---|
305 | // any Z, phi close initial phi |
---|
306 | G4double phi = 1.*deg; |
---|
307 | //t? if( theExtraPars.size() > 0 ) phi = (theExtraPars[0]+1.)*deg; |
---|
308 | G4ThreeVector pR( std::cos(phi), std::sin(phi), theWorldDim/100. ); |
---|
309 | G4double rho = (ii+0.001) * theWorldDim/numberOfPoints; |
---|
310 | pR.setRho( rho ); |
---|
311 | pR += centre; |
---|
312 | fout << pR.x() << " " << pR.y() << " " << pR.z() << G4endl; |
---|
313 | } |
---|
314 | |
---|
315 | // For division along phi |
---|
316 | centre = G4ThreeVector(0.,0.,0.); |
---|
317 | for( ii = 0; ii < numberOfPoints; ii++ ) { |
---|
318 | G4double phi = (ii+0.001) * 360*deg/numberOfPoints; |
---|
319 | // any rho (=1), any Z |
---|
320 | G4ThreeVector pR( std::cos(phi), std::sin(phi), theWorldDim/100. ); |
---|
321 | pR += centre; |
---|
322 | fout << pR.x() << " " << pR.y() << " " << pR.z() << G4endl; |
---|
323 | } |
---|
324 | |
---|
325 | // For division along Z |
---|
326 | centre = G4ThreeVector(0.,0.,2*theWorldDim); |
---|
327 | for( ii = 0; ii < numberOfPoints; ii++ ) { |
---|
328 | //any rho (=1), phi close initial phi |
---|
329 | G4double phi = 1.*deg; |
---|
330 | //t? if( theExtraPars.size() > 0 ) phi = (theExtraPars[0]+1.)*deg; |
---|
331 | G4ThreeVector pR( std::cos(phi), std::sin(phi),0.); |
---|
332 | G4double Z = -theWorldDim + (ii+0.001) * 2*theWorldDim / numberOfPoints; |
---|
333 | pR.setZ( Z ); |
---|
334 | pR += centre; |
---|
335 | fout << pR.x() << " " << pR.y() << " " << pR.z() << G4endl; |
---|
336 | } |
---|
337 | } |
---|
338 | |
---|
339 | //-------------------------------------------------------------------------- |
---|
340 | void generateScanPointsForPolycone() |
---|
341 | { |
---|
342 | generateScanPointsForTube(); |
---|
343 | } |
---|
344 | |
---|
345 | //-------------------------------------------------------------------------- |
---|
346 | // Build simple geometry: |
---|
347 | // world is |
---|
348 | G4VPhysicalVolume* BuildGeometry( SolidType solType, PVType pvType ) |
---|
349 | { |
---|
350 | G4int ii; |
---|
351 | // The world volume |
---|
352 | // |
---|
353 | G4Box *worldSolid= new G4Box ("Big Cube", theWorldDim, theWorldDim, 3*theWorldDim); |
---|
354 | |
---|
355 | G4LogicalVolume *worldLog=new G4LogicalVolume(worldSolid,0, |
---|
356 | "World",0,0,0); |
---|
357 | // Logical with no material,field, |
---|
358 | // sensitive detector or user limits |
---|
359 | |
---|
360 | G4PVPlacement *worldPhys=new G4PVPlacement(0,G4ThreeVector(0,0,0), |
---|
361 | "World",worldLog, |
---|
362 | 0,false,0); |
---|
363 | // Note: no parent pointer set |
---|
364 | |
---|
365 | // build three boxes |
---|
366 | // A set of boxes |
---|
367 | calculateParentSolid( solType ); |
---|
368 | |
---|
369 | //parent logical volumes do not depend on SolidType, PVType |
---|
370 | theParentLogs.push_back( new G4LogicalVolume(theParentSolid,0, "Parent1",0,0,0) ); |
---|
371 | theParentLogs.push_back( new G4LogicalVolume(theParentSolid,0, "Parent2",0,0,0) ); |
---|
372 | theParentLogs.push_back( new G4LogicalVolume(theParentSolid,0, "Parent3",0,0,0) ); |
---|
373 | |
---|
374 | //parent physical volumes positions do not depend on SolidType, PVType |
---|
375 | for( ii = 0; ii < 3; ii++ ) { |
---|
376 | theParentPhyss.push_back( new G4PVPlacement( 0, G4ThreeVector(0.,0.,(ii-1)*2*theWorldDim), theParentLogs[ii] , "parent", worldLog, 0, 0 ) ); |
---|
377 | } |
---|
378 | |
---|
379 | // children |
---|
380 | calculateChildSolids( solType, theParentSolid ); |
---|
381 | calculateAxis( solType ); |
---|
382 | for( ii = 0; ii < 3; ii++ ) { |
---|
383 | theChildLogs.push_back( new G4LogicalVolume(theChildSolids[ii], 0, "child",0,0,0) ); |
---|
384 | } |
---|
385 | |
---|
386 | if( pvType == pvDivision ) { |
---|
387 | for( ii = 0; ii < 3; ii++ ) { |
---|
388 | if( theDivType == 0 ) { |
---|
389 | new G4PVDivision("child", theChildLogs[ii], theParentLogs[ii], |
---|
390 | theAxis[ii], |
---|
391 | theNReplicas, |
---|
392 | theWidths[ii], |
---|
393 | theOffset ); |
---|
394 | G4cout << "division NDIVandWIDTH axis " |
---|
395 | << theNReplicas << " " << theWidths[ii]<< " " << theOffset << " " << theAxis[ii] << G4endl; |
---|
396 | } else if( theDivType == 1 ) { |
---|
397 | new G4PVDivision("child", theChildLogs[ii], theParentLogs[ii], |
---|
398 | theAxis[ii], |
---|
399 | theWidths[ii], |
---|
400 | theOffset ); |
---|
401 | G4cout << "division WIDTH " << theWidths[ii]<< " " << theOffset << G4endl; |
---|
402 | }else if( theDivType == 2 ) { |
---|
403 | new G4PVDivision("child", theChildLogs[ii], theParentLogs[ii], |
---|
404 | theAxis[ii], |
---|
405 | theNReplicas, |
---|
406 | theOffset ); |
---|
407 | G4cout << "division NDIV " << theNReplicas << " " << theOffset << G4endl; |
---|
408 | } |
---|
409 | } |
---|
410 | } else if( pvType == pvReplica ) { |
---|
411 | |
---|
412 | for( ii = 0; ii < 3; ii++ ) { |
---|
413 | new G4PVReplica("child", theChildLogs[ii], theParentLogs[ii], |
---|
414 | theAxis[ii], |
---|
415 | theNReplicas, |
---|
416 | theWidths[ii], |
---|
417 | theOffset ); |
---|
418 | G4cout << "replica " << theNReplicas << " " << theWidths[ii]<< " " << theOffset << G4endl; |
---|
419 | } |
---|
420 | } |
---|
421 | return worldPhys; |
---|
422 | } |
---|
423 | |
---|
424 | //-------------------------------------------------------------------------- |
---|
425 | // |
---|
426 | // Test LocateGlobalPointAndSetup |
---|
427 | // |
---|
428 | G4bool testG4Navigator1(G4VPhysicalVolume *pTopNode) |
---|
429 | { |
---|
430 | G4Navigator myNav; |
---|
431 | G4VPhysicalVolume *located; |
---|
432 | myNav.SetWorldVolume(pTopNode); |
---|
433 | |
---|
434 | G4String foutname = "points." + thePVTypeStr + "." + theSolidTypeStr + ".out"; |
---|
435 | std::ofstream fout(foutname); |
---|
436 | std::ifstream fin("points.lis"); |
---|
437 | G4double posX, posY, posZ; |
---|
438 | for( G4int ii = 0; ii < numberOfPoints; ii++ ) { |
---|
439 | fin >> posX >> posY >> posZ; |
---|
440 | if( fin.eof() ) break; |
---|
441 | located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(posX,posY,posZ),0,false); |
---|
442 | G4ThreeVector pos(posX,posY,posZ); |
---|
443 | fout << ii << pos << " " << located->GetName() << " " << located->GetCopyNo(); |
---|
444 | if( theSolidTypeStr == "tubs" || theSolidTypeStr == "tube" ){ |
---|
445 | // fout << " " << pos.phi()/deg << G4endl; |
---|
446 | fout << G4endl; |
---|
447 | } else { |
---|
448 | fout << G4endl; |
---|
449 | } |
---|
450 | } |
---|
451 | |
---|
452 | return true; |
---|
453 | } |
---|
454 | |
---|
455 | //-------------------------------------------------------------------------- |
---|
456 | // |
---|
457 | // Test Stepping |
---|
458 | // |
---|
459 | G4bool testG4Navigator2(G4VPhysicalVolume *pTopNode) |
---|
460 | { |
---|
461 | G4Navigator myNav; |
---|
462 | G4VPhysicalVolume *located; |
---|
463 | G4double Step,physStep,safety; |
---|
464 | G4ThreeVector* Hat = new G4ThreeVector[3]; |
---|
465 | Hat[0] = G4ThreeVector(1,0,0); |
---|
466 | Hat[1] = G4ThreeVector(0,1,0); |
---|
467 | Hat[2] = G4ThreeVector(0,0,1); |
---|
468 | |
---|
469 | myNav.SetWorldVolume(pTopNode); |
---|
470 | |
---|
471 | G4String foutname = "step." + thePVTypeStr + "." + theSolidTypeStr + ".out"; |
---|
472 | std::ofstream fout(foutname); |
---|
473 | std::ifstream fin("points.lis"); |
---|
474 | G4double posX, posY, posZ; |
---|
475 | for( G4int ii = 0; ii < numberOfPoints; ii++ ) { |
---|
476 | if( fin.eof() ) break; |
---|
477 | fin >> posX >> posY >> posZ; |
---|
478 | located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(posX,posY,posZ)); |
---|
479 | physStep=kInfinity; |
---|
480 | for( G4int jj = 0; jj < 3; jj++ ) { |
---|
481 | Step=myNav.ComputeStep(G4ThreeVector(posX,posY,posZ),Hat[jj],physStep,safety); |
---|
482 | fout << Step << " "; |
---|
483 | } |
---|
484 | fout << G4endl; |
---|
485 | } |
---|
486 | |
---|
487 | return true; |
---|
488 | } |
---|
489 | |
---|
490 | //-------------------------------------------------------------------------- |
---|
491 | void initialize() |
---|
492 | { |
---|
493 | theNReplicas = 10; |
---|
494 | theOffset = 0.; |
---|
495 | } |
---|
496 | |
---|
497 | //-------------------------------------------------------------------------- |
---|
498 | void calculateParentSolid( SolidType solType ) |
---|
499 | { |
---|
500 | if( solType == g4box ) { |
---|
501 | theParentSolid = new G4Box("parent", theWorldDim, theWorldDim, theWorldDim); |
---|
502 | }else if( solType == g4trd ) { |
---|
503 | theParentSolid = new G4Trd("parent", theWorldDim/2., theWorldDim, theWorldDim/2., theWorldDim, theWorldDim); |
---|
504 | }else if( solType == g4tube ) { |
---|
505 | theParentSolid = new G4Tubs("parent", 0., theWorldDim, theWorldDim, 0., 360.*deg); |
---|
506 | }else if( solType == g4tubs ) { |
---|
507 | theParentSolid = new G4Tubs("parent", 0., theWorldDim, theWorldDim, theExtraPars[0]*deg, theExtraPars[1]*deg); |
---|
508 | }else if( solType == g4cone ) { |
---|
509 | theParentSolid = new G4Cons("parent", 1.E-6*mm, theWorldDim, 1.E-6*mm, theWorldDim/2., theWorldDim, 0., 360.*deg); |
---|
510 | }else if( solType == g4cons ) { |
---|
511 | theParentSolid = new G4Cons("parent", 1.E-6*mm, theWorldDim/2., 1.E-6*mm, theWorldDim, theWorldDim, theExtraPars[0]*deg, theExtraPars[1]*deg); |
---|
512 | }else if( solType == g4polycone ) { |
---|
513 | |
---|
514 | G4double* zPlane = new G4double[4]; |
---|
515 | zPlane[0] = -theWorldDim; zPlane[1] = -theWorldDim/2.; zPlane[2] = theWorldDim/2; zPlane[3] = theWorldDim; |
---|
516 | G4double* rInner = new G4double[4]; |
---|
517 | rInner[0] = theWorldDim/10.; rInner[1] = 0.; rInner[2] = 0.; rInner[3] = theWorldDim/10.; |
---|
518 | G4double* rOuter = new G4double[4]; |
---|
519 | rOuter[0] = theWorldDim; rOuter[1] = theWorldDim*0.9; rOuter[2] = theWorldDim*0.9; rOuter[3] = theWorldDim; |
---|
520 | |
---|
521 | theParentSolid = new G4Polycone("parent", 0., 360.*deg, 4, zPlane, rInner, rOuter); |
---|
522 | }else if( solType == g4polyhedra ) { |
---|
523 | |
---|
524 | G4double* zPlane = new G4double[4]; |
---|
525 | zPlane[0] = -theWorldDim; zPlane[1] = -theWorldDim/2.; zPlane[2] = theWorldDim/2; zPlane[3] = theWorldDim; |
---|
526 | G4double* rInner = new G4double[4]; |
---|
527 | rInner[0] = theWorldDim/10.; rInner[1] = 0.; rInner[2] = 0.; rInner[3] = theWorldDim/10.; |
---|
528 | G4double* rOuter = new G4double[4]; |
---|
529 | rOuter[0] = theWorldDim; rOuter[1] = theWorldDim*0.9; rOuter[2] = theWorldDim*0.9; rOuter[3] = theWorldDim; |
---|
530 | |
---|
531 | theParentSolid = new G4Polyhedra("parent", 0., 360.*deg, theNReplicas, 4, zPlane, rInner, rOuter); |
---|
532 | } |
---|
533 | } |
---|
534 | |
---|
535 | //-------------------------------------------------------------------------- |
---|
536 | void calculateChildSolids( SolidType solType, G4VSolid* pSolid ) |
---|
537 | { |
---|
538 | if( solType == g4box ) { |
---|
539 | theWidths.push_back( 2*theWorldDim / theNReplicas ); |
---|
540 | theWidths.push_back( 2*theWorldDim / theNReplicas ); |
---|
541 | theWidths.push_back( 2*theWorldDim / theNReplicas ); |
---|
542 | |
---|
543 | theChildSolids.push_back( new G4Box("child", theWorldDim/theNReplicas, theWorldDim, theWorldDim) ); |
---|
544 | theChildSolids.push_back( new G4Box("child", theWorldDim, theWorldDim/theNReplicas, theWorldDim) ); |
---|
545 | theChildSolids.push_back( new G4Box("child", theWorldDim, theWorldDim, theWorldDim/theNReplicas) ); |
---|
546 | |
---|
547 | }else if( solType == g4trd ) { |
---|
548 | theWidths.push_back( 1.5*theWorldDim / theNReplicas ); |
---|
549 | theWidths.push_back( 1.5*theWorldDim / theNReplicas ); |
---|
550 | theWidths.push_back( 2*theWorldDim / theNReplicas ); |
---|
551 | |
---|
552 | theChildSolids.push_back( new G4Trap("child", theWorldDim/2./theNReplicas, theWorldDim/theNReplicas, theWorldDim, theWorldDim, theWorldDim) ); // solid dimensions will be recalculated |
---|
553 | theChildSolids.push_back( new G4Trap("child", theWorldDim, theWorldDim, theWorldDim/2./theNReplicas, theWorldDim/theNReplicas, theWorldDim) ); // solid dimensions will be recalculated |
---|
554 | theChildSolids.push_back( new G4Trd("child", theWorldDim, theWorldDim, theWorldDim, theWorldDim, theWorldDim/theNReplicas) ); |
---|
555 | |
---|
556 | }else if( solType == g4tube || solType == g4tubs ) { |
---|
557 | G4Tubs* pTubs = reinterpret_cast<G4Tubs*>( pSolid ); |
---|
558 | |
---|
559 | theWidths.push_back( (pTubs->GetInnerRadius()+pTubs->GetOuterRadius()) / theNReplicas ); |
---|
560 | theWidths.push_back( pTubs->GetDeltaPhiAngle() / theNReplicas ); |
---|
561 | theWidths.push_back( 2*pTubs->GetZHalfLength() / theNReplicas ); |
---|
562 | |
---|
563 | theChildSolids.push_back( new G4Tubs("child", pTubs->GetInnerRadius(), pTubs->GetInnerRadius()+theWidths[0], pTubs->GetZHalfLength(), pTubs->GetStartPhiAngle(), pTubs->GetDeltaPhiAngle() )); |
---|
564 | theChildSolids.push_back( new G4Tubs("child", pTubs->GetInnerRadius(), pTubs->GetOuterRadius(), pTubs->GetZHalfLength(), pTubs->GetStartPhiAngle(), theWidths[1] )); |
---|
565 | theChildSolids.push_back( new G4Tubs("child", pTubs->GetInnerRadius(), pTubs->GetOuterRadius(), theWidths[2]/2., pTubs->GetStartPhiAngle(), pTubs->GetDeltaPhiAngle() )); |
---|
566 | }else if( solType == g4cone || solType == g4cons ) { |
---|
567 | G4Cons* pCons = reinterpret_cast<G4Cons*>( pSolid ); |
---|
568 | |
---|
569 | theWidths.push_back( (pCons->GetInnerRadiusMinusZ()+pCons->GetOuterRadiusMinusZ()) / theNReplicas ); |
---|
570 | theWidths.push_back( pCons->GetDeltaPhiAngle() / theNReplicas ); |
---|
571 | theWidths.push_back( 2*pCons->GetZHalfLength() / theNReplicas ); |
---|
572 | |
---|
573 | G4double fwidthPlus = 0.; |
---|
574 | if( pCons->GetInnerRadiusMinusZ() != 0. ) fwidthPlus = theWidths[0] * pCons->GetInnerRadiusPlusZ() / pCons->GetInnerRadiusMinusZ(); |
---|
575 | theChildSolids.push_back( new G4Cons("child", pCons->GetInnerRadiusMinusZ(), pCons->GetInnerRadiusMinusZ()+theWidths[0], pCons->GetInnerRadiusPlusZ(), pCons->GetInnerRadiusPlusZ()+fwidthPlus, pCons->GetZHalfLength(), pCons->GetStartPhiAngle(), pCons->GetDeltaPhiAngle() )); |
---|
576 | theChildSolids.push_back( new G4Cons("child", pCons->GetInnerRadiusMinusZ(), pCons->GetOuterRadiusMinusZ(), pCons->GetInnerRadiusPlusZ(), pCons->GetOuterRadiusPlusZ(), pCons->GetZHalfLength(), pCons->GetStartPhiAngle(), theWidths[1] )); |
---|
577 | theChildSolids.push_back( new G4Cons("child", pCons->GetInnerRadiusMinusZ(), pCons->GetOuterRadiusMinusZ(), pCons->GetInnerRadiusPlusZ(), pCons->GetOuterRadiusPlusZ(), theWidths[2]/2., pCons->GetStartPhiAngle(), pCons->GetDeltaPhiAngle() )); |
---|
578 | |
---|
579 | }else if( solType == g4polycone ) { |
---|
580 | G4Polycone* pPCone = reinterpret_cast<G4Polycone*>( pSolid ); |
---|
581 | G4PolyconeHistorical* origparam = pPCone->GetOriginalParameters() ; |
---|
582 | G4int nz = origparam->Num_z_planes; |
---|
583 | theWidths.push_back( (origparam->Rmax[0] - origparam->Rmin[0] ) / theNReplicas ); |
---|
584 | theWidths.push_back( (pPCone->GetEndPhi() - pPCone->GetStartPhi() ) / theNReplicas ); |
---|
585 | theWidths.push_back( (origparam->Z_values[1] - origparam->Z_values[0] ) / theNReplicas ); |
---|
586 | |
---|
587 | G4double* zPlane = new G4double[4]; |
---|
588 | zPlane[0] = -theWorldDim; zPlane[1] = -theWorldDim/2.; zPlane[2] = theWorldDim/2; zPlane[3] = theWorldDim; |
---|
589 | G4double* rInner = new G4double[4]; |
---|
590 | rInner[0] = theWorldDim/10.; rInner[1] = 0.; rInner[2] = 0.; rInner[3] = theWorldDim/10.; |
---|
591 | G4double* rOuter = new G4double[4]; |
---|
592 | rOuter[0] = theWorldDim; rOuter[1] = theWorldDim*0.9; rOuter[2] = theWorldDim*0.9; rOuter[3] = theWorldDim; |
---|
593 | |
---|
594 | theChildSolids.push_back( new G4Polycone("child", 0., 360.*deg, nz, zPlane, rInner, rOuter) ); |
---|
595 | theChildSolids.push_back( new G4Polycone("child", 0., 360.*deg, nz, zPlane, rInner, rOuter) ); |
---|
596 | theChildSolids.push_back( new G4Polycone("child", 0., 360.*deg, nz, zPlane, rInner, rOuter) ); |
---|
597 | |
---|
598 | }else if( solType == g4polyhedra ) { |
---|
599 | G4Polyhedra* pPhedra = reinterpret_cast<G4Polyhedra*>( pSolid ); |
---|
600 | G4PolyhedraHistorical* origparam = pPhedra->GetOriginalParameters() ; |
---|
601 | G4int ns = origparam->numSide; |
---|
602 | G4int nz = origparam->Num_z_planes; |
---|
603 | theWidths.push_back( (origparam->Rmax[0] - origparam->Rmin[0] ) / theNReplicas ); |
---|
604 | theWidths.push_back( (pPhedra->GetEndPhi() - pPhedra->GetStartPhi() ) / theNReplicas ); |
---|
605 | theWidths.push_back( (origparam->Z_values[1] - origparam->Z_values[0] ) / theNReplicas ); |
---|
606 | |
---|
607 | G4double* zPlane = new G4double[4]; |
---|
608 | zPlane[0] = -theWorldDim; zPlane[1] = -theWorldDim/2.; zPlane[2] = theWorldDim/2; zPlane[3] = theWorldDim; |
---|
609 | G4double* rInner = new G4double[4]; |
---|
610 | rInner[0] = theWorldDim/10.; rInner[1] = 0.; rInner[2] = 0.; rInner[3] = theWorldDim/10.; |
---|
611 | G4double* rOuter = new G4double[4]; |
---|
612 | rOuter[0] = theWorldDim; rOuter[1] = theWorldDim*0.9; rOuter[2] = theWorldDim*0.9; rOuter[3] = theWorldDim; |
---|
613 | |
---|
614 | theChildSolids.push_back( new G4Polyhedra("child", 0., 360.*deg, ns, nz, zPlane, rInner, rOuter) ); |
---|
615 | theChildSolids.push_back( new G4Polyhedra("child", 0., 360.*deg, ns, nz, zPlane, rInner, rOuter) ); |
---|
616 | theChildSolids.push_back( new G4Polyhedra("child", 0., 360.*deg, ns, nz, zPlane, rInner, rOuter) ); |
---|
617 | |
---|
618 | } |
---|
619 | for( size_t ii = 0; ii < theChildSolids.size(); ii++) { |
---|
620 | G4cout << theChildSolids[0] << "child solid built 0 " << *theChildSolids[0] << G4endl; |
---|
621 | G4cout << " theChildSolid built " << ii << G4endl; |
---|
622 | G4cout << *theChildSolids[ii] << G4endl; |
---|
623 | } |
---|
624 | } |
---|
625 | |
---|
626 | //-------------------------------------------------------------------------- |
---|
627 | void calculateAxis( SolidType solType ) |
---|
628 | { |
---|
629 | if( solType == g4box ) { |
---|
630 | theAxis.push_back( kXAxis ); |
---|
631 | theAxis.push_back( kYAxis ); |
---|
632 | theAxis.push_back( kZAxis ); |
---|
633 | }else if( solType == g4trd ) { |
---|
634 | theAxis.push_back( kXAxis ); |
---|
635 | theAxis.push_back( kYAxis ); |
---|
636 | theAxis.push_back( kZAxis ); |
---|
637 | }else if( solType == g4tube || solType == g4tubs ) { |
---|
638 | theAxis.push_back( kRho ); |
---|
639 | theAxis.push_back( kPhi ); |
---|
640 | theAxis.push_back( kZAxis ); |
---|
641 | }else if( solType == g4cone || solType == g4cons || solType == g4polycone || solType == g4polyhedra ) { |
---|
642 | theAxis.push_back( kRho ); |
---|
643 | theAxis.push_back( kPhi ); |
---|
644 | theAxis.push_back( kZAxis ); |
---|
645 | } |
---|
646 | } |
---|
647 | |
---|
648 | //-------------------------------------------------------------------------- |
---|
649 | SolidType getSolidType( const G4String& st ) |
---|
650 | { |
---|
651 | SolidType stype = g4box; |
---|
652 | |
---|
653 | if( st == "box" ) { |
---|
654 | stype = g4box; |
---|
655 | }else if( st == "trd" ){ |
---|
656 | stype = g4trd; |
---|
657 | }else if( st == "tube" ){ |
---|
658 | stype = g4tube; |
---|
659 | }else if( st == "tubs" ){ |
---|
660 | stype = g4tubs; |
---|
661 | }else if( st == "cone" ){ |
---|
662 | stype = g4cone; |
---|
663 | }else if( st == "cons" ){ |
---|
664 | stype = g4cons; |
---|
665 | }else if( st == "polycone" ){ |
---|
666 | stype = g4polycone; |
---|
667 | }else if( st == "polyhedra" ){ |
---|
668 | stype = g4polyhedra; |
---|
669 | } else { |
---|
670 | G4Exception(" Input solid type can only be 'box', 'trd', 'tube', 'tubs', 'cons', 'polycone', 'polyhedra' " ); |
---|
671 | } |
---|
672 | return stype; |
---|
673 | } |
---|
674 | |
---|
675 | //-------------------------------------------------------------------------- |
---|
676 | G4int checkNumberOfArguments( const G4String& st, G4int narg ) |
---|
677 | { |
---|
678 | G4bool nok = false; |
---|
679 | if( narg == 1 ) return 0; |
---|
680 | G4int nArgExcess; |
---|
681 | if( st == "box" || st == "trd" || st == "tube" || st == "cone" ){ |
---|
682 | nArgExcess = narg - 3; |
---|
683 | }else if( st == "tubs" || st == "cons" ){ |
---|
684 | nArgExcess = narg - 5; |
---|
685 | }else if( st == "polycone" ){ |
---|
686 | nArgExcess = narg - 3; |
---|
687 | }else if( st == "polyhedra" ){ |
---|
688 | nArgExcess = narg - 3; |
---|
689 | } else { |
---|
690 | nArgExcess = 10; |
---|
691 | } |
---|
692 | |
---|
693 | G4int divTypeSet = 0; |
---|
694 | if( nArgExcess == 0 ) { |
---|
695 | nok = true; |
---|
696 | } else if( nArgExcess == 1 ) { |
---|
697 | divTypeSet = narg - 1; //last argument sets the division type |
---|
698 | nok = true; |
---|
699 | } |
---|
700 | |
---|
701 | if( !nok ){ |
---|
702 | G4cerr << " Number Of arguments " << narg << " for solid type " << st << G4endl; |
---|
703 | G4Exception(" Number of arguments incorrect for input solid type, please check method checkNumberOfArguments " ); |
---|
704 | } |
---|
705 | |
---|
706 | return divTypeSet; |
---|
707 | } |
---|
708 | |
---|
709 | //-------------------------------------------------------------------------- |
---|
710 | PVType getPVType( const G4String& pvt ) |
---|
711 | { |
---|
712 | G4cout << pvt << G4endl; |
---|
713 | PVType vtype = pvPlacement; |
---|
714 | |
---|
715 | if( pvt == "divis"){ |
---|
716 | vtype = pvDivision; |
---|
717 | }else if( pvt == "repli") { |
---|
718 | vtype = pvReplica; |
---|
719 | }else if( pvt == "place") { |
---|
720 | vtype = pvPlacement; |
---|
721 | } else { |
---|
722 | G4Exception(" Input PV type can only be 'divis', 'repli' or 'place' " ); |
---|
723 | } |
---|
724 | return vtype; |
---|
725 | } |
---|