source: trunk/source/geometry/divisions/test/testG4PVDivision.cc @ 1347

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

geant4 tag 9.4
File size: 27.5 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: 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
68G4String theSolidTypeStr;
69G4String thePVTypeStr;
70G4double theWorldDim = 1*m;
71G4int numberOfPoints = 1000;
72G4int theNReplicas;
73G4double theWidth;
74G4double theOffset;
75G4VSolid* theParentSolid;
76std::vector<G4LogicalVolume*> theParentLogs;
77std::vector<G4VPhysicalVolume*> theParentPhyss;
78std::vector<G4VSolid*> theChildSolids;
79std::vector<G4LogicalVolume*> theChildLogs;
80std::vector<EAxis> theAxis;
81std::vector<G4double> theWidths;
82std::vector<G4double> theExtraPars;
83G4int theDivType;
84
85enum SolidType{g4box, g4trd, g4tube, g4tubs, g4cone, g4cons, g4polycone, g4polyhedra };
86enum PVType{pvDivision, pvReplica, pvPlacement };
87
88//--------------------------------------------------------------------------
89void initialize();
90void calculateParentSolid( SolidType solType );
91void calculateChildSolids( SolidType solType, G4VSolid* pSolid );
92void calculateAxis( SolidType solType );
93void buildOutputName( SolidType& soltype, PVType pvtype );
94void starttest( const std::vector<G4String>& vsarg );
95G4int checkNumberOfArguments( const G4String& st, G4int narg );
96PVType getPVType( const G4String& pvt );
97SolidType getSolidType( const G4String& st );
98void generateRandomPoints();
99void generateScanPointsForBox();
100void generateScanPointsForTube();
101void generateScanPointsForTrd();
102void generateScanPointsForPolycone();
103G4VPhysicalVolume* BuildGeometry( SolidType solType, PVType pvType );
104G4bool testG4Navigator1(G4VPhysicalVolume *pTopNode);
105G4bool 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//--------------------------------------------------------------------------
114int 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//--------------------------------------------------------------------------
126void 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//--------------------------------------------------------------------------
196void 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//--------------------------------------------------------------------------
210void 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//--------------------------------------------------------------------------
252void 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//--------------------------------------------------------------------------
294void 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//--------------------------------------------------------------------------
340void generateScanPointsForPolycone()
341{
342  generateScanPointsForTube();
343}
344
345//--------------------------------------------------------------------------
346// Build simple geometry:
347//  world is
348G4VPhysicalVolume* 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//
428G4bool 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//
459G4bool 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//--------------------------------------------------------------------------
491void initialize()
492{
493 theNReplicas = 10;
494 theOffset = 0.;
495}
496
497//--------------------------------------------------------------------------
498void 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//--------------------------------------------------------------------------
536void 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//--------------------------------------------------------------------------
627void 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//--------------------------------------------------------------------------
649SolidType 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//--------------------------------------------------------------------------
676G4int 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//--------------------------------------------------------------------------
710PVType 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}
Note: See TracBrowser for help on using the repository browser.