Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/persistency/ascii/src/G4tgbVolume.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

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