source: trunk/source/persistency/ascii/src/G4tgbVolume.cc @ 1348

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

geant4 tag 9.4
File size: 51.1 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: G4tgbVolume.cc,v 1.23 2010/11/05 08:58:36 gcosmo Exp $
28// GEANT4 tag $Name: geant4-09-04-ref-00 $
29//
30//
31// class G4tgbVolume
32
33// History:
34// - Created.                                 P.Arce, CIEMAT (November 2007)
35// -------------------------------------------------------------------------
36
37#include "G4tgbVolume.hh"
38#include "G4tgbVolumeMgr.hh"
39#include "G4tgbMaterialMgr.hh"
40#include "G4tgbRotationMatrixMgr.hh"
41#include "G4tgbPlaceParamLinear.hh"
42#include "G4tgbPlaceParamSquare.hh"
43#include "G4tgbPlaceParamCircle.hh"
44
45#include "G4tgrSolid.hh"
46#include "G4tgrSolidBoolean.hh"
47#include "G4tgrVolume.hh"
48#include "G4tgrVolumeDivision.hh"
49#include "G4tgrVolumeAssembly.hh"
50#include "G4tgrVolumeMgr.hh"
51#include "G4tgrPlace.hh"
52#include "G4tgrPlaceSimple.hh"
53#include "G4tgrPlaceDivRep.hh"
54#include "G4tgrPlaceParameterisation.hh"
55#include "G4tgrUtils.hh"
56
57#include "G4VSolid.hh"
58#include "G4UnionSolid.hh"
59#include "G4SubtractionSolid.hh"
60#include "G4IntersectionSolid.hh"
61#include "G4LogicalVolume.hh"
62#include "G4VPhysicalVolume.hh"
63#include "G4PVPlacement.hh"
64#include "G4PVDivision.hh"
65#include "G4PVReplica.hh"
66#include "G4PVParameterised.hh"
67#include "G4Box.hh"
68#include "G4Tubs.hh"
69#include "G4Cons.hh"
70#include "G4Trap.hh"
71#include "G4Sphere.hh"
72#include "G4Orb.hh"
73#include "G4Trd.hh"
74#include "G4Para.hh"
75#include "G4Torus.hh"
76#include "G4Hype.hh"
77#include "G4Polycone.hh"
78#include "G4Polyhedra.hh"
79#include "G4EllipticalTube.hh"
80#include "G4Ellipsoid.hh"
81#include "G4EllipticalCone.hh"
82#include "G4Hype.hh"
83#include "G4Tet.hh"
84#include "G4TwistedBox.hh"
85#include "G4TwistedTrap.hh"
86#include "G4TwistedTrd.hh"
87#include "G4TwistedTubs.hh"
88#include "G4AssemblyVolume.hh"
89#include "G4BREPSolidBox.hh"
90#include "G4BREPSolidCylinder.hh"
91#include "G4BREPSolidCone.hh"
92#include "G4BREPSolidSphere.hh"
93#include "G4BREPSolidTorus.hh"
94#include "G4BREPSolidPCone.hh"
95#include "G4BREPSolidPolyhedra.hh"
96#include "G4BREPSolidOpenPCone.hh"
97#include "G4TessellatedSolid.hh"
98#include "G4TriangularFacet.hh"
99#include "G4QuadrangularFacet.hh"
100#include "G4ExtrudedSolid.hh"
101
102#include "G4VisExtent.hh"
103#include "G4Material.hh"
104#include "G4RotationMatrix.hh"
105#include "G4ReflectionFactory.hh"
106
107#include "G4VisAttributes.hh"
108#include "G4RegionStore.hh"
109#include "G4tgrMessenger.hh"
110#include "G4UIcommand.hh"
111#include "G4GeometryTolerance.hh"
112
113//-------------------------------------------------------------------
114G4tgbVolume::G4tgbVolume()
115  : theTgrVolume(0), theG4AssemblyVolume(0)
116{
117}
118
119
120//-------------------------------------------------------------------
121G4tgbVolume::~G4tgbVolume()
122{
123}
124
125
126//-------------------------------------------------------------------
127G4tgbVolume::G4tgbVolume( G4tgrVolume* vol)
128{
129  theTgrVolume = vol;
130  theG4AssemblyVolume = 0;
131}
132
133
134//-------------------------------------------------------------------
135void G4tgbVolume::ConstructG4Volumes( const G4tgrPlace* place,
136                                      const G4LogicalVolume* parentLV )
137{
138#ifdef G4VERBOSE
139  if( G4tgrMessenger::GetVerboseLevel() >= 2 )
140  {
141    G4cout << G4endl <<  "@@@ G4tgbVolume::ConstructG4Volumes - " << GetName() << G4endl;
142    if( place && parentLV ) G4cout << "   place in LV " << parentLV->GetName() << G4endl;
143  }
144#endif
145  G4tgbVolumeMgr* g4vmgr = G4tgbVolumeMgr::GetInstance();
146  G4LogicalVolume* logvol = g4vmgr->FindG4LogVol( GetName() );
147  G4bool bFirstCopy = false;
148  if( (logvol == 0) ) 
149  {
150    bFirstCopy = true;
151    if( theTgrVolume->GetType() != "VOLDivision" )
152    {
153      //--- If first time build solid and LogVol
154      G4VSolid* solid = FindOrConstructG4Solid( theTgrVolume->GetSolid() ); 
155      if( solid != 0 )   // for G4AssemblyVolume it is 0
156      {
157        g4vmgr->RegisterMe( solid );
158        logvol = ConstructG4LogVol( solid );
159        g4vmgr->RegisterMe( logvol );
160        g4vmgr->RegisterChildParentLVs( logvol, parentLV ); 
161      }
162      else
163      {
164        return;
165      }
166    }
167    else
168    {
169      return;
170    }
171  } 
172  //--- Construct PhysVol
173  G4VPhysicalVolume* physvol = ConstructG4PhysVol( place, logvol, parentLV );
174  if( physvol != 0 )  // 0 for G4AssemblyVolumes
175  {
176    g4vmgr->RegisterMe( physvol );
177
178    if( logvol == 0 ) // case of divisions
179    {
180      logvol = physvol->GetLogicalVolume();
181    }
182  }
183  else 
184  {
185    return;
186  }
187
188  //--- If first copy build children placements in this LogVol
189  if(bFirstCopy)
190  {
191    std::pair<G4mmapspl::iterator, G4mmapspl::iterator> children
192      = G4tgrVolumeMgr::GetInstance()->GetChildren( GetName() );
193    G4mmapspl::iterator cite; 
194    for( cite = children.first; cite != children.second; cite++ )
195    {
196      //----- Call G4tgrPlace ->constructG4Volumes
197      //---- find G4tgbVolume corresponding to the G4tgrVolume
198      //     pointed by G4tgrPlace
199      G4tgrPlace* pl = const_cast<G4tgrPlace*>((*cite).second);
200      G4tgbVolume* svol = g4vmgr->FindVolume( pl->GetVolume()->GetName() );
201      //--- find copyNo
202#ifdef G4VERBOSE
203      if( G4tgrMessenger::GetVerboseLevel() >= 2 )
204      {
205        G4cout << " G4tgbVolume::ConstructG4Volumes - construct daughter " <<  pl->GetVolume()->GetName() << " # " << pl->GetCopyNo() << G4endl;
206      }
207#endif
208      svol->ConstructG4Volumes( pl, logvol );
209    }
210  }
211
212}
213
214
215
216//-------------------------------------------------------------------
217G4VSolid* G4tgbVolume::FindOrConstructG4Solid( const G4tgrSolid* sol ) 
218{
219  G4double angularTolerance = G4GeometryTolerance::GetInstance()
220                            ->GetAngularTolerance();
221
222  if( sol == 0 ) { return 0; }
223
224#ifdef G4VERBOSE
225  if( G4tgrMessenger::GetVerboseLevel() >= 2 )
226  {
227    G4cout << " G4tgbVolume::FindOrConstructG4Solid():" << G4endl
228           << "   SOLID = " << sol << G4endl
229           << "   " << sol->GetName() << " of type " << sol->GetType()
230           << G4endl;
231  }
232#endif
233
234  //----- Check if solid exists already
235  G4VSolid* solid = G4tgbVolumeMgr::GetInstance()
236                    ->FindG4Solid( sol->GetName() );
237  if( solid ) { return solid; }
238
239  // Give 'sol' as Boolean solids needs to call this method twice
240
241#ifdef G4VERBOSE
242  if( G4tgrMessenger::GetVerboseLevel() >= 2 )
243  {
244    G4cout << " G4tgbVolume::FindOrConstructG4Solid() - "
245           << sol->GetSolidParams().size() << G4endl;
246  }
247#endif
248 
249  std::vector<G4double> solParam;
250
251  // In case of BOOLEAN solids, solidParams are taken from components
252
253  if( sol->GetSolidParams().size() == 1)
254  { 
255    solParam = * sol->GetSolidParams()[ 0 ];
256  }
257
258  //----------- instantiate the appropiate G4VSolid type
259  G4String stype = sol->GetType();
260  G4String sname = sol->GetName();
261
262  if( stype == "BOX" )
263  {
264    CheckNoSolidParams( stype, 3, solParam.size() );
265    solid = new G4Box( sname, solParam[0], solParam[1], solParam[2] ); 
266
267  }
268  else if( stype == "TUBE" )
269  {
270    CheckNoSolidParams( stype, 3, solParam.size() );
271    solid = new G4Tubs( sname, solParam[0], solParam[1], solParam[2],
272                        0.*deg, 360.*deg );
273  }
274  else if( stype == "TUBS" )
275  {
276    CheckNoSolidParams( stype, 5, solParam.size() );
277    G4double phiDelta = solParam[4];
278    if( std::fabs(phiDelta - twopi) < angularTolerance ) { phiDelta = twopi; }
279    solid = new G4Tubs( sname, solParam[0], solParam[1], solParam[2],
280                        solParam[3], phiDelta );
281  }
282  else if( stype == "TRAP" )
283  {
284    if( solParam.size() == 11 )
285    {
286      solid = new G4Trap( sname, solParam[0], solParam[1], solParam[2],
287                          solParam[3], solParam[4], solParam[5], solParam[6],
288                          solParam[7], solParam[8], solParam[9], solParam[10] );
289    }
290    else if( solParam.size() == 4 )
291    {
292      solid = new G4Trap( sname, solParam[0], solParam[1]/deg,
293                          solParam[2]/deg, solParam[3]);
294    }
295    else
296    {
297      G4String ErrMessage1 = "Solid type " + stype;
298      G4String ErrMessage2 = " should have 11 or 4 parameters,\n";
299      G4String ErrMessage3 = "and it has "
300                         + G4UIcommand::ConvertToString(G4int(solParam.size()));
301      G4String ErrMessage = ErrMessage1 + ErrMessage2 + ErrMessage3 + " !";
302      G4Exception("G4tgbVolume::FindOrConstructG4Solid()",
303                  "InvalidSetup", FatalException, ErrMessage);
304      return 0;
305    }
306   
307  }
308  else if( stype == "TRD" )
309  {
310    CheckNoSolidParams( stype, 5, solParam.size() );
311    solid = new G4Trd( sname, solParam[0], solParam[1], solParam[2],
312                       solParam[3], solParam[4] );
313  }
314  else if( stype == "PARA" )
315  {
316    CheckNoSolidParams( stype, 6, solParam.size() );
317    solid = new G4Para( sname, solParam[0], solParam[1], solParam[2],
318                        solParam[3], solParam[4], solParam[5] );
319  }
320  else if( stype == "CONE" )
321  {
322    CheckNoSolidParams( stype, 5, solParam.size() );
323    solid = new G4Cons( sname, solParam[0], solParam[1], solParam[2],
324                        solParam[3], solParam[4], 0., 360.*deg);
325  }
326  else if( stype == "CONS" )
327  {
328    CheckNoSolidParams( stype, 7, solParam.size() );
329    G4double phiDelta = solParam[6];
330    if( std::fabs(phiDelta - twopi) < angularTolerance ) { phiDelta = twopi; }
331    solid = new G4Cons( sname, solParam[0], solParam[1], solParam[2],
332                        solParam[3], solParam[4], solParam[5], phiDelta);
333  }
334  else if( stype == "SPHERE" )
335  {
336    CheckNoSolidParams( stype, 6, solParam.size() );
337    G4double phiDelta = solParam[3];
338    if( std::fabs(phiDelta - twopi) < angularTolerance ) { phiDelta = twopi; }
339    G4double thetaDelta = solParam[5];
340    if( std::fabs(thetaDelta - pi) < angularTolerance ) { thetaDelta = pi; }
341    solid = new G4Sphere( sname, solParam[0], solParam[1], solParam[2],
342                          phiDelta, solParam[4], thetaDelta);
343  }
344  else if( stype == "ORB" )
345  {
346    CheckNoSolidParams( stype, 1, solParam.size() );
347    solid = new G4Orb( sname, solParam[0] );
348  }
349  else if( stype == "TORUS" )
350  {
351    CheckNoSolidParams( stype, 5, solParam.size() );
352    G4double phiDelta = solParam[4];
353    if( std::fabs(phiDelta - twopi) < angularTolerance ) { phiDelta = twopi; }
354    solid = new G4Torus( sname, solParam[0], solParam[1], solParam[2],
355                         solParam[3], phiDelta );
356  }
357  else if( stype == "POLYCONE" )
358  {
359    size_t nplanes = size_t(solParam[2]);
360    G4bool genericPoly = false;
361    if( solParam.size() == 3+nplanes*3 )
362    { 
363      genericPoly = true;
364    }
365    else if( solParam.size() == 3+nplanes*2 )
366    { 
367      genericPoly = false;
368    }
369    else
370    {
371      G4String Err1 = "Solid type " + stype + " should have ";
372      G4String Err2 = G4UIcommand::ConvertToString(G4int(3+nplanes*3))
373                    + " (Z,Rmin,Rmax)\n";
374      G4String Err3 = "or " + G4UIcommand::ConvertToString(G4int(3+nplanes*2));
375      G4String Err4 = " (RZ corners) parameters,\n";
376      G4String Err5 = "and it has "
377                    +  G4UIcommand::ConvertToString(G4int(solParam.size()));
378      G4String ErrMessage = Err1 + Err2 + Err3 + Err4 + Err5 + " !";
379      G4Exception("G4tgbVolume::FindOrConstructG4Solid()",
380                  "InvalidSetup", FatalException, ErrMessage);
381      return 0;
382    }
383
384    if( genericPoly )
385    {
386      std::vector<G4double>* z_p = new std::vector<G4double>;
387      std::vector<G4double>* rmin_p = new std::vector<G4double>;
388      std::vector<G4double>* rmax_p = new std::vector<G4double>;
389      for( size_t ii = 0; ii < nplanes; ii++ )
390      {
391        (*z_p).push_back( solParam[3+3*ii] );
392        (*rmin_p).push_back( solParam[3+3*ii+1] );
393        (*rmax_p).push_back(  solParam[3+3*ii+2] );
394      }
395      G4double phiTotal = solParam[1];
396      if( std::fabs(phiTotal - twopi) < angularTolerance ) { phiTotal = twopi; }
397      solid = new G4Polycone( sname, solParam[0], phiTotal, // start,delta-phi
398                              nplanes, // sections
399                              &((*z_p)[0]), &((*rmin_p)[0]), &((*rmax_p)[0]));
400    }
401    else
402    {
403      std::vector<G4double>* R_c = new std::vector<G4double>;
404      std::vector<G4double>* Z_c = new std::vector<G4double>;
405      for( size_t ii = 0; ii < nplanes; ii++ )
406      {
407        (*R_c).push_back( solParam[3+2*ii] );
408        (*Z_c).push_back( solParam[3+2*ii+1] );
409      }
410      G4double phiTotal = solParam[1];
411      if( std::fabs(phiTotal - twopi) < angularTolerance ) { phiTotal = twopi; }
412      solid = new G4Polycone( sname, solParam[0], phiTotal, // start,delta-phi
413                              nplanes, // sections
414                              &((*R_c)[0]), &((*Z_c)[0]));
415    }
416
417  }
418  else if( stype == "POLYHEDRA" )
419  {
420    size_t nplanes = size_t(solParam[3]);
421    G4bool genericPoly = false;
422    if( solParam.size() == 4+nplanes*3 )
423    { 
424      genericPoly = true;
425    }
426    else if( solParam.size() == 4+nplanes*2 )
427    { 
428      genericPoly = false;
429    }
430    else
431    {
432      G4String Err1 = "Solid type " + stype + " should have ";
433      G4String Err2 = G4UIcommand::ConvertToString(G4int(4+nplanes*3))
434                    + " (Z,Rmin,Rmax)\n";
435      G4String Err3 = "or " + G4UIcommand::ConvertToString(G4int(4+nplanes*2));
436      G4String Err4 = " (RZ corners) parameters,\n";
437      G4String Err5 = "and it has "
438                    + G4UIcommand::ConvertToString(G4int(solParam.size()));
439      G4String ErrMessage = Err1 + Err2 + Err3 + Err4 + Err5 + " !";
440      G4Exception("G4tgbVolume::FindOrConstructG4Solid()",
441                  "InvalidSetup", FatalException, ErrMessage);
442      return 0;
443    }
444   
445    if( genericPoly )
446    {
447      std::vector<G4double>* z_p = new std::vector<G4double>;
448      std::vector<G4double>* rmin_p = new std::vector<G4double>;
449      std::vector<G4double>* rmax_p = new std::vector<G4double>;
450      for( size_t ii = 0; ii < nplanes; ii++ )
451      {
452        (*z_p).push_back( solParam[4+3*ii] );
453        (*rmin_p).push_back( solParam[4+3*ii+1] );
454        (*rmax_p).push_back(  solParam[4+3*ii+2] );
455      }
456      G4double phiTotal = solParam[1];
457      if( std::fabs(phiTotal - twopi) < angularTolerance ) { phiTotal = twopi; }
458      solid = new G4Polyhedra( sname, solParam[0], phiTotal,
459                               G4int(solParam[2]), nplanes,
460                               &((*z_p)[0]), &((*rmin_p)[0]), &((*rmax_p)[0]));
461    }
462    else
463    {
464      std::vector<G4double>* R_c = new std::vector<G4double>;
465      std::vector<G4double>* Z_c = new std::vector<G4double>;
466      for( size_t ii = 0; ii < nplanes; ii++ )
467      {
468        (*R_c).push_back( solParam[4+2*ii] );
469        (*Z_c).push_back( solParam[4+2*ii+1] );
470      }
471      G4double phiTotal = solParam[1];
472      if( std::fabs(phiTotal - twopi) < angularTolerance ) { phiTotal = twopi; }
473      solid = new G4Polyhedra( sname, solParam[0], phiTotal,
474                               G4int(solParam[2]), nplanes,
475                               &((*R_c)[0]), &((*Z_c)[0]));
476    }
477  }
478  else if( stype == "ELLIPTICALTUBE" )
479  {
480    CheckNoSolidParams( stype, 3, solParam.size() );
481    solid = new G4EllipticalTube( sname, solParam[0], solParam[1], solParam[2]);
482  }
483  else if( stype == "ELLIPSOID" )
484  {
485    CheckNoSolidParams( stype, 5, solParam.size() );
486    solid = new G4Ellipsoid( sname, solParam[0], solParam[1], solParam[2],
487                             solParam[3], solParam[4] );
488  }
489  else if( stype == "ELLIPTICALCONE" )
490  {
491    CheckNoSolidParams( stype, 4, solParam.size() );
492    solid = new G4EllipticalCone( sname, solParam[0], solParam[1],
493                                  solParam[2], solParam[3] );
494  }
495  else if( stype == "HYPE" )
496  {
497    CheckNoSolidParams( stype, 5, solParam.size() );
498    solid = new G4Hype( sname, solParam[0], solParam[1], solParam[2],
499                        solParam[3], solParam[4] );
500  }
501  else if( stype == "TET" )
502  {
503    CheckNoSolidParams( stype, 12, solParam.size() );
504    G4ThreeVector anchor(solParam[0], solParam[1], solParam[2]);
505    G4ThreeVector p2(solParam[3], solParam[4], solParam[5]);
506    G4ThreeVector p3(solParam[6], solParam[7], solParam[8]);
507    G4ThreeVector p4(solParam[9], solParam[10], solParam[11]);
508    solid = new G4Tet( sname, anchor, p2, p3, p4 );
509  }
510  else if( stype == "TWISTEDBOX" )
511  {
512    CheckNoSolidParams( stype, 4, solParam.size() );
513    solid = new G4TwistedBox( sname, solParam[0], solParam[1],
514                              solParam[2], solParam[3]);
515  }
516  else if( stype == "TWISTEDTRAP" )
517  {
518    CheckNoSolidParams( stype, 11, solParam.size() );
519    solid = new G4TwistedTrap( sname, solParam[0], solParam[1], solParam[2],
520                        solParam[3], solParam[4], solParam[5], solParam[6],
521                        solParam[7], solParam[8], solParam[9], solParam[10] );
522  }
523  else if( stype == "TWISTEDTRD" )
524  {
525    CheckNoSolidParams( stype, 6, solParam.size() );
526    solid = new G4TwistedTrd( sname, solParam[0], solParam[1], solParam[2],
527                              solParam[3], solParam[4], solParam[5]);
528  }
529  else if( stype == "TWISTEDTUBS" )
530  {
531    CheckNoSolidParams( stype, 5, solParam.size() );
532    G4double phiTotal = solParam[4];
533    if( std::fabs(phiTotal - twopi) < angularTolerance ) { phiTotal = twopi; }
534    solid = new G4TwistedTubs( sname, solParam[0], solParam[1], solParam[2],
535                               solParam[3], phiTotal);
536  }
537  else if( stype == "BREPBOX" )   // EntityType is = "Closed_Shell"
538  {
539    CheckNoSolidParams( stype, 24, solParam.size() );
540    std::vector<G4Point3D> points;
541    for( size_t ii = 0; ii < 8; ii++ )
542    {
543      points.push_back( G4Point3D(solParam[ii*3+0],
544                                  solParam[ii*3+1],
545                                  solParam[ii*3+2]) );
546    }
547    solid = new G4BREPSolidBox( sname, points[0], points[1], points[2],
548                                points[3], points[4], points[5], points[6],
549                                points[7] );
550  }
551  else if( stype == "BREPCYLINDER" )   // EntityType is = "Closed_Shell"
552  {
553    CheckNoSolidParams( stype, 11, solParam.size() );
554    solid = new G4BREPSolidCylinder( sname,
555                  G4ThreeVector( solParam[0], solParam[1], solParam[2] ),
556                  G4ThreeVector( solParam[3], solParam[4], solParam[5] ),
557                  G4ThreeVector( solParam[6], solParam[7], solParam[8] ),
558                  solParam[9], solParam[10] );
559  }
560  else if( stype == "BREPCONE" )   // EntityType is = "Closed_Shell"
561  {
562    CheckNoSolidParams( stype, 12, solParam.size() );
563    solid = new G4BREPSolidCone( sname,
564                  G4ThreeVector( solParam[0], solParam[1], solParam[2] ),
565                  G4ThreeVector( solParam[3], solParam[4], solParam[5] ),
566                  G4ThreeVector( solParam[6], solParam[7], solParam[8] ),
567                  solParam[9], solParam[10], solParam[11] );
568  }
569  else if( stype == "BREPSPHERE" )   // EntityType is = "Closed_Shell"
570  {
571    CheckNoSolidParams( stype, 10, solParam.size() );
572    solid = new G4BREPSolidSphere( sname,
573                  G4ThreeVector( solParam[0], solParam[1], solParam[2] ),
574                  G4ThreeVector( solParam[3], solParam[4], solParam[5] ),
575                  G4ThreeVector( solParam[6], solParam[7], solParam[8] ),
576                  solParam[9] );
577
578  }
579  else if( stype == "BREPTORUS" )   // EntityType is = "Closed_Shell"
580  {
581    CheckNoSolidParams( stype, 11, solParam.size() );
582    solid = new G4BREPSolidTorus( sname,
583                  G4ThreeVector( solParam[0], solParam[1], solParam[2] ),
584                  G4ThreeVector( solParam[3], solParam[4], solParam[5] ),
585                  G4ThreeVector( solParam[6], solParam[7], solParam[8] ),
586                  solParam[9], solParam[10] );
587  }
588  else if( stype == "BREPPCONE" )   // EntityType is = "Closed_Shell"
589  {
590    size_t nplanes = size_t(solParam[2]);
591    CheckNoSolidParams( stype, 4+3*nplanes, solParam.size() );
592    std::vector<G4double>* z_p = new std::vector<G4double>;
593    std::vector<G4double>* rmin_p = new std::vector<G4double>;
594    std::vector<G4double>* rmax_p = new std::vector<G4double>;
595    for( size_t ii = 0; ii < nplanes; ii++ )
596    {
597      (*z_p).push_back( solParam[4+3*ii] );
598      (*rmin_p).push_back( solParam[4+3*ii+1] );
599      (*rmax_p).push_back(  solParam[4+3*ii+2] );
600    }
601    G4double phiTotal = solParam[1];
602    if( std::fabs(phiTotal - twopi) < angularTolerance ) { phiTotal = twopi; }
603    CheckNoSolidParams( stype, 12, solParam.size() );
604    solid = new G4BREPSolidPCone( sname, solParam[0], phiTotal, // start,dph
605                                  nplanes,                      // sections
606                                  solParam[3],                  // z_start
607                                  &((*z_p)[0]), &((*rmin_p)[0]),
608                                  &((*rmax_p)[0]));
609  }
610  else if( stype == "BREPPOLYHEDRA" )   // EntityType is = "Closed_Shell"
611  {
612    size_t nplanes = size_t(solParam[3]);
613    CheckNoSolidParams( stype, 5+3*nplanes, solParam.size() );
614    std::vector<G4double>* z_p = new std::vector<G4double>;
615    std::vector<G4double>* rmin_p = new std::vector<G4double>;
616    std::vector<G4double>* rmax_p = new std::vector<G4double>;
617    for( size_t ii = 0; ii < nplanes; ii++ )
618    {
619      (*z_p).push_back( solParam[5+3*ii] );
620      (*rmin_p).push_back( solParam[5+3*ii+1] );
621      (*rmax_p).push_back(  solParam[5+3*ii+2] );
622    }
623    G4double phiTotal = solParam[1];
624    if( std::fabs(phiTotal - twopi) < angularTolerance ) { phiTotal = twopi; }
625    CheckNoSolidParams( stype, 12, solParam.size() );
626    solid = new G4BREPSolidPolyhedra( sname, solParam[0], phiTotal, // start,dph
627                                      G4int(solParam[2]), // sides
628                                      nplanes,            // sections
629                                      solParam[4],        // z_start
630                                      &((*z_p)[0]), &((*rmin_p)[0]),
631                                      &((*rmax_p)[0]));
632  }
633  else if( stype == "BREPOPENPCONE" )   // EntityType is = "Closed_Shell"
634  {
635    size_t nplanes = size_t(solParam[2]);
636    std::vector<G4double>* z_p = new std::vector<G4double>;
637    std::vector<G4double>* rmin_p = new std::vector<G4double>;
638    std::vector<G4double>* rmax_p = new std::vector<G4double>;
639    for( size_t ii = 0; ii < nplanes; ii++ )
640    {
641      (*z_p).push_back( solParam[4+3*ii] );
642      (*rmin_p).push_back( solParam[4+3*ii+1] );
643      (*rmax_p).push_back(  solParam[4+3*ii+2] );
644    }
645    G4double phiTotal = solParam[1];
646    if( std::fabs(phiTotal - twopi) < angularTolerance ) { phiTotal = twopi; }
647    CheckNoSolidParams( stype, 12, solParam.size() );
648    solid = new G4BREPSolidOpenPCone( sname, solParam[0], phiTotal, // start,dph
649                                      nplanes,                      // sections
650                                      solParam[3],                  // z_start
651                                      &((*z_p)[0]), &((*rmin_p)[0]),
652                                      &((*rmax_p)[0]));
653  }
654  else if( stype == "TESSELLATED" )
655  {
656    G4int nFacets = G4int(solParam[0]);
657    G4int jj = 0;
658    solid = new G4TessellatedSolid(sname);
659    G4TessellatedSolid* solidTS = (G4TessellatedSolid*)(solid);
660    G4VFacet* facet=0;
661   
662    for( G4int ii = 0; ii < nFacets; ii++){
663      G4int nPoints = G4int(solParam[jj+1]);
664      if( G4int(solParam.size()) < jj + nPoints*3 + 2 ) {
665        G4String Err1 = "Too small number of parameters in tesselated solid, it should be at least "  + G4UIcommand::ConvertToString(jj + nPoints*3 + 2 );
666        G4String Err2 = " facet number " + G4UIcommand::ConvertToString(ii);
667        G4String Err3 = " number of parameters is " + G4UIcommand::ConvertToString(G4int(solParam.size()));
668        G4String ErrMessage = Err1 + Err2 + Err3 + " !";
669        G4Exception("G4tgbVolume::FindOrConstructG4Solid()",
670                    "InvalidSetup", FatalException, ErrMessage);
671        return 0;
672      }
673     
674      if( nPoints == 3 ) 
675      {
676        G4ThreeVector pt0(solParam[jj+2],solParam[jj+3],solParam[jj+4]);
677        G4ThreeVector vt1(solParam[jj+5],solParam[jj+6],solParam[jj+7]);
678        G4ThreeVector vt2(solParam[jj+8],solParam[jj+9],solParam[jj+10]);
679        G4FacetVertexType vertexType = ABSOLUTE;
680        if( solParam[jj+11] == 0 ) 
681        {
682          vertexType = ABSOLUTE;
683        } 
684        else if( solParam[jj+11] == 0 ) 
685        {
686          vertexType = RELATIVE;
687        } 
688        else 
689        {
690          G4String Err1 = "Wrong number of vertex type in tesselated solid, it should be 0 =ABSOLUTE) or 1 (=RELATIVE)";
691          G4String Err2 = " facet number " + G4UIcommand::ConvertToString(G4int(ii));
692          G4String Err3 = " vertex type is " + G4UIcommand::ConvertToString(solParam[jj+11]);
693          G4String ErrMessage = Err1 + Err2 + Err3 + " !";
694          G4Exception("G4tgbVolume::FindOrConstructG4Solid()",
695                      "InvalidSetup", FatalException, ErrMessage);
696          return 0;
697        }
698        facet = new G4TriangularFacet( pt0, vt1, vt2, vertexType );
699      } 
700      else if( nPoints == 4 ) 
701      {
702        G4ThreeVector pt0(solParam[jj+2],solParam[jj+3],solParam[jj+4]);
703        G4ThreeVector vt1(solParam[jj+5],solParam[jj+6],solParam[jj+7]);
704        G4ThreeVector vt2(solParam[jj+8],solParam[jj+9],solParam[jj+10]);
705        G4ThreeVector vt3(solParam[jj+11],solParam[jj+12],solParam[jj+13]);
706        G4FacetVertexType vertexType = ABSOLUTE;
707        if( solParam[jj+14] == 0 ) 
708        {
709          vertexType = ABSOLUTE;
710        }
711        else if( solParam[jj+14] == 0 ) 
712        {
713          vertexType = RELATIVE;
714        } 
715        else 
716        {
717          G4String Err1 = "Wrong number of vertex type in tesselated solid, it should be 0 =ABSOLUTE) or 1 (=RELATIVE)";
718          G4String Err2 = " facet number " + G4UIcommand::ConvertToString(G4int(ii));
719          G4String Err3 = " vertex type is " + G4UIcommand::ConvertToString(solParam[jj+14]);
720          G4String ErrMessage = Err1 + Err2 + Err3 + " !";
721          G4Exception("G4tgbVolume::FindOrConstructG4Solid()",
722                      "InvalidSetup", FatalException, ErrMessage);
723          return 0;
724        }
725        facet = new G4QuadrangularFacet( pt0, vt1, vt2, vt3, vertexType );
726      } 
727      else 
728      {
729        G4String Err1 = "Wrong number of points in tesselated solid, it should be 3 or 4";
730        G4String Err2 = " facet number " + G4UIcommand::ConvertToString(G4int(ii));
731        G4String Err3 = " number of points is " + G4UIcommand::ConvertToString(G4int(nPoints));
732        G4String ErrMessage = Err1 + Err2 + Err3 + " !";
733        G4Exception("G4tgbVolume::FindOrConstructG4Solid()",
734                    "InvalidSetup", FatalException, ErrMessage);
735        return 0;
736      }
737     
738      solidTS->AddFacet( facet );
739      jj += nPoints*3 + 2;
740    }
741   
742  } 
743  else if( stype == "EXTRUDED" ) 
744  {
745    std::vector<G4TwoVector> polygonList;
746    std::vector<G4ExtrudedSolid::ZSection> zsectionList;
747    G4int nPolygons = G4int(solParam[0]);
748    G4int ii = 1;
749    G4int nMax = nPolygons*2+1;
750    for( ;ii < nMax; ii+=2 )
751    {
752      polygonList.push_back( G4TwoVector(solParam[ii],solParam[ii+1]) );
753    }
754    G4int nZSections = G4int(solParam[ii]);
755    nMax = nPolygons*2 + nZSections*4 + 2; 
756    ii++;
757    for( ; ii < nMax; ii+=4 )
758    {
759      G4TwoVector offset(solParam[ii+1],solParam[ii+2]);
760      zsectionList.push_back( G4ExtrudedSolid::ZSection(solParam[ii],offset,solParam[ii+3]) );
761    }
762    solid = new G4ExtrudedSolid( sname, polygonList, zsectionList );
763   
764  }
765  else if( stype.substr(0,7) == "Boolean" )
766  {
767    const G4tgrSolidBoolean* solb = dynamic_cast<const G4tgrSolidBoolean*>(sol);
768    if (!solb)
769    {
770      G4Exception("G4tgbVolume::FindOrConstructG4Solid()",
771                  "InvalidSetup", FatalException, "Invalid Solid pointer");
772      return 0;
773    }
774    G4VSolid* sol1 = FindOrConstructG4Solid( solb->GetSolid(0));
775    G4VSolid* sol2 = FindOrConstructG4Solid( solb->GetSolid(1));
776    G4RotationMatrix* relRotMat = G4tgbRotationMatrixMgr::GetInstance()
777      ->FindOrBuildG4RotMatrix( sol->GetRelativeRotMatName() );
778    G4ThreeVector relPlace = solb->GetRelativePlace();
779
780    if( stype == "Boolean_UNION" )
781    {
782      solid = new G4UnionSolid( sname, sol1, sol2, relRotMat, relPlace );
783    }
784    else if( stype == "Boolean_SUBTRACTION" )
785    {
786      solid = new G4SubtractionSolid( sname, sol1, sol2, relRotMat, relPlace );
787    }
788    else if( stype == "Boolean_INTERSECTION" )
789    {
790      solid = new G4IntersectionSolid( sname, sol1, sol2, relRotMat, relPlace );
791    }
792    else
793    {
794      G4String ErrMessage = "Unknown Boolean type " + stype;
795      G4Exception("G4tgbVolume::FindOrConstructG4Solid()",
796                  "InvalidSetup", FatalException, ErrMessage);
797      return 0;
798    }
799  }
800  else
801  {
802    G4String ErrMessage = "Solids of type " + stype
803                        + " not implemented yet, sorry...";
804    G4Exception("G4tgbVolume::FindOrConstructG4Solid()", "NotImplemented",
805                FatalException, ErrMessage);
806    return 0;
807  } 
808 
809#ifdef G4VERBOSE
810  if( G4tgrMessenger::GetVerboseLevel() >= 2 )
811  {
812    G4cout << " G4tgbVolume::FindOrConstructG4Solid()" << G4endl
813           << "   Created solid " << sname
814           << " of type " << solid->GetEntityType() << G4endl; 
815  }
816#endif
817
818#ifdef G4VERBOSE
819    if( G4tgrMessenger::GetVerboseLevel() >= 1 )
820    {
821      G4cout << " Constructing new G4Solid: " 
822             << *solid << G4endl;
823    }
824#endif
825 
826  return solid;
827}
828
829//-------------------------------------------------------------------
830void G4tgbVolume::CheckNoSolidParams( const G4String& solidType,
831                                      const unsigned int NoParamExpected,
832                                      const unsigned int NoParam )
833{
834  if( NoParamExpected != NoParam )
835  {
836    G4String Err1 = "Solid type " + solidType + " should have ";
837    G4String Err2 = G4UIcommand::ConvertToString(G4int(NoParamExpected))
838                  + " parameters,\n";
839    G4String Err3 = "and it has "
840                  + G4UIcommand::ConvertToString(G4int(NoParam));
841    G4String ErrMessage = Err1 + Err2 + Err3 + " !";
842    G4Exception("G4tgbVolume::CheckNoSolidParams()", "InvalidSetup",
843                FatalException, ErrMessage);
844  }
845}
846
847
848//-------------------------------------------------------------------
849G4LogicalVolume* G4tgbVolume::ConstructG4LogVol( const G4VSolid* solid )
850{
851  G4LogicalVolume* logvol;
852
853#ifdef G4VERBOSE
854  if( G4tgrMessenger::GetVerboseLevel() >= 2 )
855  {
856    G4cout << " G4tgbVolume::ConstructG4LogVol() - " << GetName() << G4endl;
857  }
858#endif
859
860  //----------- Get the material first
861  G4Material* mate = G4tgbMaterialMgr::GetInstance()
862              ->FindOrBuildG4Material( theTgrVolume->GetMaterialName() );
863  if( mate == 0 )
864  {
865    G4String ErrMessage = "Material not found "
866                        + theTgrVolume->GetMaterialName()
867                        + " for volume " + GetName() + ".";
868    G4Exception("G4tgbVolume::ConstructG4LogVol()", "InvalidSetup",
869                FatalException, ErrMessage);
870  }
871#ifdef G4VERBOSE
872  if( G4tgrMessenger::GetVerboseLevel() >= 2 )
873  {
874    G4cout << " G4tgbVolume::ConstructG4LogVol() -"
875           << " Material constructed: " << mate->GetName() << G4endl; 
876  }
877#endif
878 
879  //---------- Construct the LV
880  logvol = new G4LogicalVolume( const_cast<G4VSolid*>(solid),
881                                const_cast<G4Material*>(mate), GetName() );
882
883#ifdef G4VERBOSE
884    if( G4tgrMessenger::GetVerboseLevel() >= 1 )
885    {
886      G4cout << " Constructing new G4LogicalVolume: " 
887             << logvol->GetName() << " mate " << mate->GetName() << G4endl;
888    }
889#endif
890 
891  //---------- Set visibility and colour
892  if( !GetVisibility() || GetColour()[0] != -1 )
893  {
894    G4VisAttributes* visAtt = new G4VisAttributes();
895#ifdef G4VERBOSE
896    if( G4tgrMessenger::GetVerboseLevel() >= 1 )
897    {
898      G4cout << " Constructing new G4VisAttributes: " 
899             << *visAtt << G4endl;
900    }
901#endif
902 
903    if( !GetVisibility() )
904    {
905      visAtt->SetVisibility( GetVisibility() );
906    }
907    else if( GetColour()[0] != -1 )
908    {
909      // this else should not be necessary, because if the visibility
910      // is set to off, colour should have no effect. But it does not
911      // work: if you set colour and vis off, it is visualized!?!?!?
912
913      const G4double* col = GetColour();
914      if( col[3] == -1. )
915      {
916        visAtt->SetColour( G4Colour(col[0],col[1],col[2]));
917      }
918      else
919      {
920        visAtt->SetColour( G4Colour(col[0],col[1],col[2],col[3]));
921      }
922    }
923    logvol->SetVisAttributes(visAtt);
924  }
925
926#ifdef G4VERBOSE
927  if( G4tgrMessenger::GetVerboseLevel() >= 2 )
928  {
929    G4cout << " G4tgbVolume::ConstructG4LogVol() -"
930           << " Created logical volume: " << GetName() << G4endl;
931  }
932#endif
933
934  return logvol;
935}
936
937
938//-------------------------------------------------------------------
939G4VPhysicalVolume*
940G4tgbVolume::ConstructG4PhysVol( const G4tgrPlace* place,
941                                 const G4LogicalVolume* currentLV,
942                                 const G4LogicalVolume* parentLV )
943{
944  G4VPhysicalVolume* physvol = 0;
945  G4int copyNo;
946 
947  //----- Case of placement of top volume
948  if( place == 0 )
949  {
950#ifdef G4VERBOSE
951    if( G4tgrMessenger::GetVerboseLevel() >= 2 )
952    {
953      G4cout << " G4tgbVolume::ConstructG4PhysVol() - World: "
954             << GetName() << G4endl;
955    }
956#endif
957    physvol = new G4PVPlacement(0, G4ThreeVector(),
958                                const_cast<G4LogicalVolume*>(currentLV),
959                                GetName(), 0, false, 0, 
960                                theTgrVolume->GetCheckOverlaps());
961#ifdef G4VERBOSE
962    if( G4tgrMessenger::GetVerboseLevel() >= 1 )
963    {
964      G4cout << " Constructing new : G4PVPlacement " 
965             << physvol->GetName() << G4endl;
966    }
967#endif
968  }
969  else
970  { 
971    copyNo = place->GetCopyNo();
972
973#ifdef G4VERBOSE
974    if( G4tgrMessenger::GetVerboseLevel() >= 2 )
975    {
976      G4cout << " G4tgbVolume::ConstructG4PhysVol() - " << GetName() << G4endl
977             << "   inside " << parentLV->GetName() << " copy No: " << copyNo
978             << " of type= " << theTgrVolume->GetType() << G4endl
979             << "   placement type= " << place->GetType() << G4endl;
980    }
981#endif
982   
983    if( theTgrVolume->GetType() == "VOLSimple" )
984    {
985      //----- Get placement
986#ifdef G4VERBOSE
987      if( G4tgrMessenger::GetVerboseLevel() >= 2 )
988      {
989        G4cout << " G4tgbVolume::ConstructG4PhysVol() - Placement type = "
990               << place->GetType() << G4endl;
991      }
992#endif
993     
994      //--------------- If it is  G4tgrPlaceSimple
995      if( place->GetType() == "PlaceSimple" )
996      {
997        //----- Get rotation matrix
998        G4tgrPlaceSimple* placeSimple = (G4tgrPlaceSimple*)place; 
999        G4String rmName = placeSimple->GetRotMatName();
1000
1001        G4RotationMatrix* rotmat = G4tgbRotationMatrixMgr::GetInstance()
1002                                 ->FindOrBuildG4RotMatrix( rmName );
1003        //----- Place volume in mother
1004        G4double check = (rotmat->colX().cross(rotmat->colY()))*rotmat->colZ();
1005        G4double tol = 1.0e-3;
1006        //---- Check that matrix is ortogonal
1007        if (1-std::abs(check)>tol)
1008        {
1009          G4cerr << " Matrix : " << rmName << " " << rotmat->colX()
1010                 << " " << rotmat->colY() << " " << rotmat->colZ() << G4endl
1011                 << " product x X y * z = " << check << " x X y "
1012                 << rotmat->colX().cross(rotmat->colY()) << G4endl;
1013          G4String ErrMessage = "Rotation is not ortogonal " + rmName + " !";
1014          G4Exception("G4tgbVolume::ConstructG4PhysVol()",
1015                      "InvalidSetup", FatalException, ErrMessage);
1016          //---- Check if it is reflection
1017        }
1018        else if (1+check<=tol)
1019        {
1020          G4Translate3D transl = place->GetPlacement();
1021          G4Transform3D trfrm  = transl * G4Rotate3D(*rotmat);
1022          physvol = (G4ReflectionFactory::Instance()->Place(trfrm, GetName(),
1023                     const_cast<G4LogicalVolume*>(currentLV),
1024                     const_cast<G4LogicalVolume*>(parentLV),
1025                     false, copyNo, false )).first;
1026        }
1027        else
1028        {
1029#ifdef G4VERBOSE
1030          if( G4tgrMessenger::GetVerboseLevel() >= 1 )
1031          {
1032            G4cout << "Construction new G4VPhysicalVolume"
1033                   << " through G4ReflectionFactory " << GetName() 
1034                   << " in volume " << parentLV->GetName() 
1035                   << " copyNo " << copyNo
1036                   << " position " << place->GetPlacement() 
1037                   << " ROT " << rotmat->colX() 
1038                   << " " << rotmat->colY() 
1039                   << " " << rotmat->colZ() << G4endl;
1040          }
1041#endif
1042          physvol = new G4PVPlacement( rotmat, place->GetPlacement(),
1043                                       const_cast<G4LogicalVolume*>(currentLV),
1044                                       GetName(),
1045                                       const_cast<G4LogicalVolume*>(parentLV),
1046                                       false, copyNo, 
1047                                       theTgrVolume->GetCheckOverlaps());
1048        }
1049       
1050        //--------------- If it is G4tgrPlaceParam
1051      }
1052      else if( place->GetType() == "PlaceParam" )
1053      {
1054        G4tgrPlaceParameterisation* dp = (G4tgrPlaceParameterisation*)(place);
1055
1056        //----- See what parameterisation type
1057#ifdef G4VERBOSE
1058        if( G4tgrMessenger::GetVerboseLevel() >= 2 )
1059        {
1060          G4cout << " G4tgbVolume::ConstructG4PhysVol() -" << G4endl
1061                 << "   param: " << GetName() << " in " <<  parentLV->GetName()
1062                 << " param type= " << dp->GetParamType() << G4endl;
1063        }
1064#endif
1065       
1066        G4tgbPlaceParameterisation * param=0;
1067       
1068        if( (dp->GetParamType() == "CIRCLE")
1069         || (dp->GetParamType() == "CIRCLE_XY")
1070         || (dp->GetParamType() == "CIRCLE_XZ")
1071         || (dp->GetParamType() == "CIRCLE_YZ") )
1072        { 
1073          param = new G4tgbPlaceParamCircle(dp);
1074         
1075        } 
1076        else if( (dp->GetParamType() == "LINEAR")
1077                || (dp->GetParamType() == "LINEAR_X")
1078                || (dp->GetParamType() == "LINEAR_Y")
1079                || (dp->GetParamType() == "LINEAR_Z") )
1080        {   
1081          param = new G4tgbPlaceParamLinear(dp);
1082         
1083        } 
1084        else if( (dp->GetParamType() == "SQUARE")
1085                || (dp->GetParamType() == "SQUARE_XY")
1086                || (dp->GetParamType() == "SQUARE_XZ")
1087                || (dp->GetParamType() == "SQUARE_YZ") )
1088        {
1089          param = new G4tgbPlaceParamSquare(dp);
1090        }
1091        else
1092        {
1093          G4String ErrMessage = "Parameterisation has wrong type, TYPE: "
1094                              + G4String(dp->GetParamType()) + " !";
1095          G4Exception("G4tgbVolume::ConstructG4PhysVol", "WrongArgument",
1096                      FatalException, ErrMessage);
1097          return 0;
1098        }
1099#ifdef G4VERBOSE
1100        if( G4tgrMessenger::GetVerboseLevel() >= 1 )
1101        {
1102          G4cout << " G4tgbVolume::ConstructG4PhysVol() -" << G4endl
1103                 << "   New G4PVParameterised: " << GetName() << " vol "
1104                 << currentLV->GetName() << " in vol " << parentLV->GetName()
1105                 << " axis " << param->GetAxis() << " nCopies "
1106                 << param->GetNCopies() << G4endl;
1107        }
1108#endif
1109        physvol = new G4PVParameterised(GetName(),
1110                                        const_cast<G4LogicalVolume*>(currentLV),
1111                                        const_cast<G4LogicalVolume*>(parentLV),
1112                                        EAxis(param->GetAxis()),
1113                                        param->GetNCopies(), param);
1114#ifdef G4VERBOSE
1115    if( G4tgrMessenger::GetVerboseLevel() >= 1 )
1116    {
1117      G4cout << " Constructing new G4PVParameterised: " 
1118             << physvol->GetName() << " in volume " << parentLV->GetName() 
1119             << " N copies " << param->GetNCopies() 
1120             << " axis " << param->GetAxis() << G4endl;
1121    }
1122#endif
1123
1124      }
1125      else if( place->GetType() == "PlaceReplica" )
1126      {
1127        //--------------- If it is  PlaceReplica
1128        G4tgrPlaceDivRep* dpr = (G4tgrPlaceDivRep*)place;
1129
1130#ifdef G4VERBOSE
1131        if( G4tgrMessenger::GetVerboseLevel() >= 2 )
1132        {
1133          G4cout << " G4tgbVolume::ConstructG4PhysVol() -" << G4endl
1134                 << "   replica" << " " << currentLV->GetName()
1135                 << " in " <<  parentLV->GetName() 
1136                 << " NDiv " << dpr->GetNDiv() << " Width " << dpr->GetWidth()
1137                 << " offset " << dpr->GetOffset() << G4endl;
1138        }
1139#endif
1140        physvol = new G4PVReplica(GetName(),
1141                                  const_cast<G4LogicalVolume*>(currentLV),
1142                                  const_cast<G4LogicalVolume*>(parentLV),
1143                                  EAxis(dpr->GetAxis()), dpr->GetNDiv(),
1144                                  dpr->GetWidth(), dpr->GetOffset());
1145#ifdef G4VERBOSE
1146    if( G4tgrMessenger::GetVerboseLevel() >= 1 )
1147    {
1148      G4cout << " Constructing new G4PVReplica: " 
1149             << currentLV->GetName()
1150             << " in " <<  parentLV->GetName() 
1151             << " NDiv " << dpr->GetNDiv() << " Width " << dpr->GetWidth()
1152             << " offset " << dpr->GetOffset() << G4endl;
1153    }
1154#endif
1155      }
1156    }
1157    else if( theTgrVolume->GetType() == "VOLDivision" )
1158    {
1159      G4tgrVolumeDivision* volr = (G4tgrVolumeDivision*)theTgrVolume;
1160      G4tgrPlaceDivRep* placeDiv = volr->GetPlaceDivision() ;
1161      G4VSolid* solid = BuildSolidForDivision( parentLV->GetSolid(), placeDiv->GetAxis() );
1162      G4Material* mate = G4tgbMaterialMgr::GetInstance()
1163                  ->FindOrBuildG4Material( theTgrVolume->GetMaterialName() );
1164      G4LogicalVolume* divLV = new G4LogicalVolume(solid,
1165                                                  const_cast<G4Material*>(mate),
1166                                                   GetName() );
1167#ifdef G4VERBOSE
1168      if( G4tgrMessenger::GetVerboseLevel() >= 1 )
1169      {
1170        G4cout << " Constructed new G4LogicalVolume for division: " 
1171               << divLV->GetName() << " mate " << mate->GetName() << G4endl;
1172      }
1173#endif
1174 
1175      G4DivType divType = placeDiv->GetDivType();
1176      switch (divType)
1177      {
1178      case DivByNdiv:
1179        physvol = new G4PVDivision(GetName(), (G4LogicalVolume*)divLV,
1180                                   const_cast<G4LogicalVolume*>(parentLV),
1181                                   placeDiv->GetAxis(), placeDiv->GetNDiv(),
1182                                   placeDiv->GetOffset());
1183#ifdef G4VERBOSE
1184        if( G4tgrMessenger::GetVerboseLevel() >= 1 )
1185        {
1186          G4cout << " Constructing new G4PVDivision by number of divisions: " 
1187                 << GetName() << " in " <<  parentLV->GetName() 
1188                 << " axis " << placeDiv->GetAxis() 
1189                 << " Ndiv " << placeDiv->GetNDiv()
1190                 << " offset " << placeDiv->GetOffset() << G4endl;
1191        }
1192#endif
1193        break;
1194      case DivByWidth:
1195        physvol = new G4PVDivision(GetName(), (G4LogicalVolume*)divLV,
1196                                   const_cast<G4LogicalVolume*>(parentLV),
1197                                   placeDiv->GetAxis(), placeDiv->GetWidth(),
1198                                   placeDiv->GetOffset());
1199#ifdef G4VERBOSE
1200        if( G4tgrMessenger::GetVerboseLevel() >= 1 )
1201        {
1202          G4cout << " Constructing new G4PVDivision by width: " 
1203                 << GetName() << " in " <<  parentLV->GetName() 
1204                 << " axis " << placeDiv->GetAxis() 
1205                 << " width " << placeDiv->GetWidth()
1206                 << " offset " << placeDiv->GetOffset() << G4endl;
1207        }
1208#endif
1209        break;
1210      case DivByNdivAndWidth:
1211        physvol = new G4PVDivision(GetName(), (G4LogicalVolume*)divLV,
1212                                   const_cast<G4LogicalVolume*>(parentLV),
1213                                   placeDiv->GetAxis(), placeDiv->GetNDiv(),
1214                                   placeDiv->GetWidth(),
1215                                   placeDiv->GetOffset());
1216#ifdef G4VERBOSE
1217        if( G4tgrMessenger::GetVerboseLevel() >= 1 )
1218        {
1219          G4cout << " Constructing new G4PVDivision by width"
1220                 << " and number of divisions: " 
1221                 << GetName() << " in " <<  parentLV->GetName() 
1222                 << " axis " << placeDiv->GetAxis() 
1223                 << " Ndiv " << placeDiv->GetNDiv()
1224                 << " width " << placeDiv->GetWidth()
1225                 << " offset " << placeDiv->GetOffset() << G4endl;
1226        }
1227#endif
1228        break;
1229      }
1230    }
1231    else if( theTgrVolume->GetType() == "VOLAssembly" )
1232    {
1233      // Define one layer as one assembly volume
1234      G4tgrVolumeAssembly * tgrAssembly = (G4tgrVolumeAssembly *)theTgrVolume;
1235
1236      if( !theG4AssemblyVolume )
1237      {
1238        theG4AssemblyVolume = new G4AssemblyVolume;
1239#ifdef G4VERBOSE
1240        if( G4tgrMessenger::GetVerboseLevel() >= 1 )
1241          {
1242            G4cout << " Constructing new G4AssemblyVolume: " 
1243                   << " number of assembly components " 
1244                   <<  tgrAssembly->GetNoComponents() << G4endl;
1245          }
1246#endif       
1247        G4tgbVolumeMgr* g4vmgr = G4tgbVolumeMgr::GetInstance();
1248        for( G4int ii = 0; ii < tgrAssembly->GetNoComponents(); ii++ )
1249        {
1250          // Rotation and translation of a plate inside the assembly
1251
1252          G4ThreeVector transl =  tgrAssembly->GetComponentPos(ii);
1253          G4String rmName = tgrAssembly->GetComponentRM(ii);
1254          G4RotationMatrix* rotmat = G4tgbRotationMatrixMgr::GetInstance()
1255                                   ->FindOrBuildG4RotMatrix( rmName );
1256         
1257          //----- Get G4LogicalVolume of component
1258          G4String lvname = tgrAssembly->GetComponentName(ii);
1259          G4LogicalVolume* logvol = g4vmgr->FindG4LogVol( lvname);
1260          if( logvol == 0 )
1261          {
1262            g4vmgr->FindVolume( lvname )->ConstructG4Volumes( 0, 0);
1263            logvol = g4vmgr->FindG4LogVol( lvname, true );
1264          }
1265          // Fill the assembly by the plates
1266          theG4AssemblyVolume->AddPlacedVolume( logvol, transl, rotmat );
1267#ifdef G4VERBOSE
1268          if( G4tgrMessenger::GetVerboseLevel() >= 1 )
1269          {
1270            G4cout << " G4AssemblyVolume->AddPlacedVolume " << ii 
1271                   << " " << logvol->GetName()
1272                   << " translation " << transl 
1273                   << " rotmat " << rotmat->colX() 
1274                   << " " << rotmat->colY() 
1275                   << " " << rotmat->colZ() << G4endl;
1276          }
1277#endif
1278        }
1279      }
1280
1281      // Rotation and Translation of the assembly inside the world
1282
1283      G4tgrPlaceSimple* placeSimple = (G4tgrPlaceSimple*)place; 
1284      G4String rmName = placeSimple->GetRotMatName();
1285      G4RotationMatrix* rotmat = G4tgbRotationMatrixMgr::GetInstance()
1286                               ->FindOrBuildG4RotMatrix( rmName );
1287      G4ThreeVector transl = place->GetPlacement();
1288
1289      G4LogicalVolume* parentLV_nonconst =
1290                       const_cast<G4LogicalVolume*>(parentLV);
1291      theG4AssemblyVolume->MakeImprint( parentLV_nonconst, transl, rotmat );
1292 
1293    }
1294    else   // If it is G4tgrVolumeAssembly
1295    {
1296      G4String ErrMessage = "Volume type not supported: "
1297                          + theTgrVolume->GetType() + ", sorry...";
1298      G4Exception("G4tgbVolume::ConstructG4PhysVol()", "NotImplemented",
1299                  FatalException, ErrMessage);
1300    }   
1301  } 
1302
1303  return physvol;
1304}
1305
1306
1307//-------------------------------------------------------------------
1308G4VSolid* G4tgbVolume::BuildSolidForDivision( G4VSolid* parentSolid, EAxis axis  )
1309{
1310  G4VSolid* solid=0;
1311  G4double redf = (parentSolid->GetExtent().GetXmax()-parentSolid->GetExtent().GetXmin());
1312  redf = std::min(redf,parentSolid->GetExtent().GetYmax()-parentSolid->GetExtent().GetYmin());
1313  redf = std::min(redf,parentSolid->GetExtent().GetZmax()-parentSolid->GetExtent().GetZmin());
1314  redf *= 0.001; //make daugther much smaller, to fit in parent
1315
1316  if( parentSolid->GetEntityType() == "G4Box" )
1317  {
1318    G4Box* psolid = (G4Box*)(parentSolid);
1319    solid = new G4Box(GetName(), psolid->GetXHalfLength()*redf,
1320                                 psolid->GetZHalfLength()*redf,
1321                                 psolid->GetZHalfLength()*redf);
1322  } 
1323  else if ( parentSolid->GetEntityType() == "G4Tubs" )
1324  {
1325    G4Tubs* psolid = (G4Tubs*)(parentSolid);
1326    solid = new G4Tubs( GetName(), psolid->GetInnerRadius()*redf,
1327                                   psolid->GetOuterRadius()*redf,
1328                                   psolid->GetZHalfLength()*redf,
1329                                   psolid->GetSPhi(), psolid->GetDPhi());
1330  } 
1331  else if ( parentSolid->GetEntityType() == "G4Cons" )
1332  {
1333    G4Cons* psolid = (G4Cons*)(parentSolid);
1334    solid = new G4Cons( GetName(), psolid->GetInnerRadiusMinusZ()*redf,
1335                                   psolid->GetOuterRadiusMinusZ()*redf,
1336                                   psolid->GetInnerRadiusPlusZ()*redf,
1337                                   psolid->GetOuterRadiusPlusZ()*redf,
1338                                   psolid->GetZHalfLength()*redf,
1339                                   psolid->GetSPhi(), psolid->GetDPhi());
1340  } 
1341  else if ( parentSolid->GetEntityType() == "G4Trd" )
1342  {
1343    G4Trd* psolid = (G4Trd*)(parentSolid);
1344    G4double mpDx1 = psolid->GetXHalfLength1();
1345    G4double mpDx2 = psolid->GetXHalfLength2();
1346
1347    if( axis == kXAxis && std::fabs(mpDx1 - mpDx2) > G4GeometryTolerance::GetInstance()->GetSurfaceTolerance() )
1348    {
1349      solid = new G4Trap( GetName(), psolid->GetZHalfLength()*redf, 
1350                          psolid->GetYHalfLength1()*redf, 
1351                          psolid->GetXHalfLength2()*redf, 
1352                          psolid->GetXHalfLength1()*redf );
1353    } 
1354    else
1355    {
1356      solid = new G4Trd( GetName(), psolid->GetXHalfLength1()*redf,
1357                         psolid->GetXHalfLength2()*redf,
1358                         psolid->GetYHalfLength1()*redf,
1359                         psolid->GetYHalfLength2()*redf,
1360                         psolid->GetZHalfLength()*redf);
1361    }
1362   
1363  } 
1364  else if ( parentSolid->GetEntityType() == "G4Para" )
1365  {
1366    G4Para* psolid = (G4Para*)(parentSolid);
1367    solid = new G4Para( GetName(), psolid->GetXHalfLength()*redf,
1368                                   psolid->GetYHalfLength()*redf,
1369                                   psolid->GetZHalfLength()*redf,
1370                                   std::atan(psolid->GetTanAlpha()),
1371                                   psolid->GetSymAxis().theta(),
1372                                   psolid->GetSymAxis().phi() ); 
1373  } 
1374  else if ( parentSolid->GetEntityType() == "G4Polycone" )
1375  {
1376    G4Polycone* psolid = (G4Polycone*)(parentSolid);
1377    G4PolyconeHistorical origParam = *(psolid->GetOriginalParameters());
1378    for( G4int ii = 0; ii < origParam.Num_z_planes; ii++ )
1379    {
1380      origParam.Rmin[ii] = origParam.Rmin[ii]*redf;
1381      origParam.Rmax[ii] = origParam.Rmax[ii]*redf;
1382    }
1383    solid = new G4Polycone( GetName(), psolid->GetStartPhi(),
1384                                       psolid->GetEndPhi(),
1385                            origParam.Num_z_planes, origParam.Z_values,
1386                            origParam.Rmin, origParam.Rmax);
1387
1388  } 
1389  else if ( parentSolid->GetEntityType() == "G4Polyhedra" )
1390  {
1391    G4Polyhedra* psolid = (G4Polyhedra*)(parentSolid);
1392    G4PolyhedraHistorical origParam = *(psolid->GetOriginalParameters());
1393    for( G4int ii = 0; ii < origParam.Num_z_planes; ii++ )
1394    {
1395      origParam.Rmin[ii] = origParam.Rmin[ii]*redf;
1396      origParam.Rmax[ii] = origParam.Rmax[ii]*redf;
1397    }
1398    solid = new G4Polyhedra( GetName(), psolid->GetStartPhi(),
1399                                        psolid->GetEndPhi(),
1400                                        psolid->GetNumSide(),
1401                             origParam.Num_z_planes, origParam.Z_values,
1402                             origParam.Rmin, origParam.Rmax);
1403  }
1404  else
1405  { 
1406    G4String ErrMessage = "Solid type not supported. VOLUME= " + GetName()
1407                        + " Solid type= " + parentSolid->GetEntityType() + "\n"
1408                        + "Only supported types are: G4Box, G4Tubs, G4Cons,"
1409                        + " G4Trd, G4Para, G4Polycone, G4Polyhedra.";
1410    G4Exception("G4tgbVolume::BuildSolidForDivision()", "NotImplemented",
1411                FatalException, ErrMessage);
1412    return 0;
1413  }
1414
1415#ifdef G4VERBOSE
1416    if( G4tgrMessenger::GetVerboseLevel() >= 1 )
1417    {
1418      G4cout << " Constructing new G4Solid for division: " 
1419             << *solid << G4endl;
1420    }
1421#endif
1422  return solid;
1423}
1424 
Note: See TracBrowser for help on using the repository browser.