Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/persistency/ascii/src/G4tgbGeometryDumper.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 // G4tgbGeometryDumper implementation
 27 //
 28 // Author: P.Arce, CIEMAT (November 2007)
 29 // --------------------------------------------------------------------
 30 
 31 #include "G4tgbGeometryDumper.hh"
 32 
 33 #include "G4tgrMessenger.hh"
 34 
 35 #include "G4UIcommand.hh"
 36 #include "G4Material.hh"
 37 #include "G4Element.hh"
 38 #include "G4VSolid.hh"
 39 #include "G4Box.hh"
 40 #include "G4Tubs.hh"
 41 #include "G4Cons.hh"
 42 #include "G4Trap.hh"
 43 #include "G4Sphere.hh"
 44 #include "G4Orb.hh"
 45 #include "G4Trd.hh"
 46 #include "G4Para.hh"
 47 #include "G4Torus.hh"
 48 #include "G4Hype.hh"
 49 #include "G4Polycone.hh"
 50 #include "G4GenericPolycone.hh"
 51 #include "G4Polyhedra.hh"
 52 #include "G4EllipticalTube.hh"
 53 #include "G4Ellipsoid.hh"
 54 #include "G4EllipticalCone.hh"
 55 #include "G4Hype.hh"
 56 #include "G4Tet.hh"
 57 #include "G4TwistedBox.hh"
 58 #include "G4TwistedTrap.hh"
 59 #include "G4TwistedTrd.hh"
 60 #include "G4TwistedTubs.hh"
 61 #include "G4MultiUnion.hh"
 62 #include "G4ScaledSolid.hh"
 63 #include "G4PVPlacement.hh"
 64 #include "G4PVParameterised.hh"
 65 #include "G4PVReplica.hh"
 66 #include "G4BooleanSolid.hh"
 67 #include "G4ReflectionFactory.hh"
 68 #include "G4ReflectedSolid.hh"
 69 #include "G4LogicalVolumeStore.hh"
 70 #include "G4PhysicalVolumeStore.hh"
 71 #include "G4GeometryTolerance.hh"
 72 #include "G4VPVParameterisation.hh"
 73 #include "G4SystemOfUnits.hh"
 74 #include <iomanip>
 75 
 76 // --------------------------------------------------------------------
 77 G4ThreadLocal G4tgbGeometryDumper* G4tgbGeometryDumper::theInstance = nullptr;
 78 
 79 // --------------------------------------------------------------------
 80 G4tgbGeometryDumper::G4tgbGeometryDumper()
 81 {
 82 }
 83 
 84 // --------------------------------------------------------------------
 85 G4tgbGeometryDumper* G4tgbGeometryDumper::GetInstance()
 86 {
 87   if(theInstance == nullptr)
 88   {
 89     theInstance = new G4tgbGeometryDumper;
 90   }
 91   return theInstance;
 92 }
 93 
 94 // --------------------------------------------------------------------
 95 void G4tgbGeometryDumper::DumpGeometry(const G4String& fname)
 96 {
 97   theFile = new std::ofstream(fname);
 98 
 99   G4VPhysicalVolume* pv = GetTopPhysVol();
100   DumpPhysVol(pv);  // dump volume and recursively it will dump all hierarchy
101 }
102 
103 // --------------------------------------------------------------------
104 G4VPhysicalVolume* G4tgbGeometryDumper::GetTopPhysVol()
105 {
106   G4PhysicalVolumeStore* pvstore = G4PhysicalVolumeStore::GetInstance();
107   G4VPhysicalVolume* pv = *(pvstore->cbegin());
108   for(;;)
109   {
110     G4LogicalVolume* lv = pv->GetMotherLogical();
111     if(lv == 0)
112     {
113       break;
114     }
115 
116     //----- look for one PV of this LV
117     for(auto ite = pvstore->cbegin(); ite != pvstore->cend(); ++ite)
118     {
119       pv = (*ite);
120       if(pv->GetLogicalVolume() == lv)
121       {
122         break;
123       }
124     }
125   }
126 
127   return pv;
128 }
129 
130 // --------------------------------------------------------------------
131 G4tgbGeometryDumper::~G4tgbGeometryDumper()
132 {
133 }
134 
135 // --------------------------------------------------------------------
136 void G4tgbGeometryDumper::DumpPhysVol(G4VPhysicalVolume* pv)
137 {
138   //--- Dump logical volume first
139   G4LogicalVolume* lv = pv->GetLogicalVolume();
140 
141   G4ReflectionFactory* reffact = G4ReflectionFactory::Instance();
142 
143   //--- It is not needed to dump _refl volumes created when parent is reflected
144   // !!WARNING : it must be avoided to reflect a volume hierarchy if children
145   //             has also been reflected, as both will have same name
146 
147   if(reffact->IsReflected(lv) && reffact->IsReflected(pv->GetMotherLogical()))
148   {
149     return;
150   }
151 
152   G4bool bVolExists = CheckIfLogVolExists(lv->GetName(), lv);
153 
154   //---- Construct this PV
155   if(pv->GetMotherLogical() != nullptr)  // not WORLD volume
156   {
157     if(!pv->IsReplicated())
158     {
159       G4String lvName = lv->GetName();
160       if(!bVolExists)
161       {
162         lvName = DumpLogVol(lv);
163       }
164       DumpPVPlacement(pv, lvName);
165     }
166     else if(pv->IsParameterised())
167     {
168       G4PVParameterised* pvparam = (G4PVParameterised*) (pv);
169       DumpPVParameterised(pvparam);
170     }
171     else
172     {
173       G4String lvName = lv->GetName();
174       if(!bVolExists)
175       {
176         lvName = DumpLogVol(lv);
177       }
178       G4PVReplica* pvrepl = (G4PVReplica*) (pv);
179       DumpPVReplica(pvrepl, lvName);
180     }
181   }
182   else
183   {
184     DumpLogVol(lv);
185   }
186 
187   if(!bVolExists)
188   {
189     //---- Construct PV's who has this LV as mother
190     std::vector<G4VPhysicalVolume*> pvChildren = GetPVChildren(lv);
191     for(auto ite = pvChildren.cbegin(); ite != pvChildren.cend(); ++ite)
192     {
193       DumpPhysVol(*ite);
194     }
195   }
196 }
197 
198 // --------------------------------------------------------------------
199 void G4tgbGeometryDumper::DumpPVPlacement(G4VPhysicalVolume* pv,
200                                           const G4String& lvName, G4int copyNo)
201 {
202   G4String pvName = pv->GetName();
203 
204   G4RotationMatrix* rotMat = pv->GetRotation();
205   if(rotMat == nullptr)
206     rotMat = new G4RotationMatrix();
207 
208   //---- Check if it is reflected
209   G4ReflectionFactory* reffact = G4ReflectionFactory::Instance();
210   G4LogicalVolume* lv          = pv->GetLogicalVolume();
211   if(reffact->IsReflected(lv))
212   {
213 #ifdef G4VERBOSE
214     if(G4tgrMessenger::GetVerboseLevel() >= 1)
215     {
216       G4cout << " G4tgbGeometryDumper::DumpPVPlacement() - Reflected volume: "
217              << pv->GetName() << G4endl;
218     }
219 #endif
220     G4ThreeVector colx = rotMat->colX();
221     G4ThreeVector coly = rotMat->colY();
222     G4ThreeVector colz = rotMat->colZ();
223     // apply a Z reflection (reflection matrix is decomposed in new
224     // reflection-free rotation + z-reflection)
225     colz *= -1.;
226     G4Rep3x3 rottemp(colx.x(), coly.x(), colz.x(), colx.y(), coly.y(), colz.y(),
227                      colx.z(), coly.z(), colz.z());
228     // matrix representation (inverted)
229     *rotMat = G4RotationMatrix(rottemp);
230     *rotMat = (*rotMat).inverse();
231     pvName += "_refl";
232   }
233   const G4String& rotName  = DumpRotationMatrix(rotMat);
234   G4ThreeVector pos = pv->GetTranslation();
235 
236   if(copyNo == -999)  // for parameterisations copy number is provided
237   {
238     copyNo = pv->GetCopyNo();
239   }
240 
241   const G4String& fullname = pvName + "#" + G4UIcommand::ConvertToString(copyNo)
242                            + "/" + pv->GetMotherLogical()->GetName();
243 
244   if(!CheckIfPhysVolExists(fullname, pv))
245   {
246     (*theFile) << ":PLACE " << SubstituteRefl(AddQuotes(lvName)) << " "
247                << copyNo << " "
248                << SubstituteRefl(AddQuotes(pv->GetMotherLogical()->GetName()))
249                << " " << AddQuotes(rotName) << " " << pos.x() << " " << pos.y()
250                << " " << pos.z() << G4endl;
251 
252     thePhysVols[fullname] = pv;
253   }
254 }
255 
256 // --------------------------------------------------------------------
257 void G4tgbGeometryDumper::DumpPVParameterised(G4PVParameterised* pv)
258 {
259   G4String pvName = pv->GetName();
260 
261   EAxis axis;
262   G4int nReplicas;
263   G4double width;
264   G4double offset;
265   G4bool consuming;
266   pv->GetReplicationData(axis, nReplicas, width, offset, consuming);
267 
268   G4VPVParameterisation* param = pv->GetParameterisation();
269 
270   G4LogicalVolume* lv             = pv->GetLogicalVolume();
271   G4VSolid* solid1st              = param->ComputeSolid(0, pv);
272   G4Material* mate1st             = param->ComputeMaterial(0, pv);
273   std::vector<G4double> params1st = GetSolidParams(solid1st);
274   std::vector<G4double> newParams;
275   G4VSolid* newSolid = solid1st;
276   G4String lvName;
277 
278   for(G4int ii = 0; ii < nReplicas; ++ii)
279   {
280     G4Material* newMate = param->ComputeMaterial(ii, pv);
281     if(solid1st->GetEntityType() == "G4Box")
282     {
283       G4Box* box = (G4Box*) (solid1st);
284       param->ComputeDimensions(*box, ii, pv);
285       newParams = GetSolidParams(box);
286       newSolid  = (G4VSolid*) box;
287     }
288     else if(solid1st->GetEntityType() == "G4Tubs")
289     {
290       G4Tubs* tubs = (G4Tubs*) (solid1st);
291       param->ComputeDimensions(*tubs, ii, pv);
292       newParams = GetSolidParams(tubs);
293       newSolid  = (G4VSolid*) tubs;
294     }
295     else if(solid1st->GetEntityType() == "G4Trd")
296     {
297       G4Trd* trd = (G4Trd*) (solid1st);
298       param->ComputeDimensions(*trd, ii, pv);
299       newParams = GetSolidParams(trd);
300       newSolid  = (G4VSolid*) trd;
301     }
302     else if(solid1st->GetEntityType() == "G4Trap")
303     {
304       G4Trap* trap = (G4Trap*) (solid1st);
305       param->ComputeDimensions(*trap, ii, pv);
306       newParams = GetSolidParams(trap);
307       newSolid  = (G4VSolid*) trap;
308     }
309     else if(solid1st->GetEntityType() == "G4Cons")
310     {
311       G4Cons* cons = (G4Cons*) (solid1st);
312       param->ComputeDimensions(*cons, ii, pv);
313       newParams = GetSolidParams(cons);
314       newSolid  = (G4VSolid*) cons;
315     }
316     else if(solid1st->GetEntityType() == "G4Sphere")
317     {
318       G4Sphere* sphere = (G4Sphere*) (solid1st);
319       param->ComputeDimensions(*sphere, ii, pv);
320       newParams = GetSolidParams(sphere);
321       newSolid  = (G4VSolid*) sphere;
322     }
323     else if(solid1st->GetEntityType() == "G4Orb")
324     {
325       G4Orb* orb = (G4Orb*) (solid1st);
326       param->ComputeDimensions(*orb, ii, pv);
327       newParams = GetSolidParams(orb);
328       newSolid  = (G4VSolid*) orb;
329     }
330     else if(solid1st->GetEntityType() == "G4Torus")
331     {
332       G4Torus* torus = (G4Torus*) (solid1st);
333       param->ComputeDimensions(*torus, ii, pv);
334       newParams = GetSolidParams(torus);
335       newSolid  = (G4VSolid*) torus;
336     }
337     else if(solid1st->GetEntityType() == "G4Para")
338     {
339       G4Para* para = (G4Para*) (solid1st);
340       param->ComputeDimensions(*para, ii, pv);
341       newParams = GetSolidParams(para);
342       newSolid  = (G4VSolid*) para;
343     }
344     else if(solid1st->GetEntityType() == "G4Polycone")
345     {
346       G4Polycone* polycone = (G4Polycone*) (solid1st);
347       param->ComputeDimensions(*polycone, ii, pv);
348       newParams = GetSolidParams(polycone);
349       newSolid  = (G4VSolid*) polycone;
350     }
351     else if(solid1st->GetEntityType() == "G4Polyhedra")
352     {
353       G4Polyhedra* polyhedra = (G4Polyhedra*) (solid1st);
354       param->ComputeDimensions(*polyhedra, ii, pv);
355       newParams = GetSolidParams(polyhedra);
356       newSolid  = (G4VSolid*) polyhedra;
357     }
358     else if(solid1st->GetEntityType() == "G4Hype")
359     {
360       G4Hype* hype = (G4Hype*) (solid1st);
361       param->ComputeDimensions(*hype, ii, pv);
362       newParams = GetSolidParams(hype);
363       newSolid  = (G4VSolid*) hype;
364     }
365     if(ii == 0 || mate1st != newMate || params1st[0] != newParams[0])
366     {
367       G4String extraName = "";
368       if(ii != 0)
369       {
370         extraName = "#" + G4UIcommand::ConvertToString(ii) + "/" +
371                     pv->GetMotherLogical()->GetName();
372       }
373       lvName = DumpLogVol(lv, extraName, newSolid, newMate);
374     }
375 
376     param->ComputeTransformation(ii, pv);
377     DumpPVPlacement(pv, lvName, ii);
378   }
379 }
380 
381 // --------------------------------------------------------------------
382 void G4tgbGeometryDumper::DumpPVReplica(G4PVReplica* pv, const G4String& lvName)
383 {
384   EAxis axis;
385   G4int nReplicas;
386   G4double width;
387   G4double offset;
388   G4bool consuming;
389   pv->GetReplicationData(axis, nReplicas, width, offset, consuming);
390   G4String axisName;
391   switch(axis)
392   {
393     case kXAxis:
394       axisName = "X";
395       break;
396     case kYAxis:
397       axisName = "Y";
398       break;
399     case kZAxis:
400       axisName = "Z";
401       break;
402     case kRho:
403       axisName = "R";
404       break;
405     case kPhi:
406       axisName = "PHI";
407       break;
408     case kRadial3D:
409     case kUndefined:
410       G4String ErrMessage =
411         "Unknown axis of replication for volume" + pv->GetName();
412       G4Exception("G4tgbGeometryDumper::DumpPVReplica", "Wrong axis ",
413                   FatalException, ErrMessage);
414       break;
415   }
416 
417   const G4String& fullname = lvName + "/" + pv->GetMotherLogical()->GetName();
418 
419   if(!CheckIfPhysVolExists(fullname, pv))
420   {
421     (*theFile) << ":REPL " << SubstituteRefl(AddQuotes(lvName)) << " "
422                << SubstituteRefl(AddQuotes(pv->GetMotherLogical()->GetName()))
423                << " " << axisName << " " << nReplicas;
424     if(axis != kPhi)
425     {
426       (*theFile) << " " << width << " " << offset << G4endl;
427     }
428     else
429     {
430       (*theFile) << " " << width / deg << "*deg"
431                  << " " << offset / deg << "*deg" << G4endl;
432     }
433 
434     thePhysVols[fullname] = pv;
435   }
436 }
437 
438 // --------------------------------------------------------------------
439 G4String G4tgbGeometryDumper::DumpLogVol(G4LogicalVolume* lv,
440                                          const G4String& extraName,
441                                          G4VSolid* solid,
442                                          G4Material* mate)
443 {
444   G4String lvName;
445 
446   if(extraName == "")  //--- take out the '_refl' in the name
447   {
448     lvName = GetObjectName(lv, theLogVols);
449   }
450   else
451   {
452     lvName = lv->GetName() + extraName;
453   }
454 
455   if(theLogVols.find(lvName) != theLogVols.cend())  // alredy dumped
456   {
457     return lvName;
458   }
459 
460   if(solid == nullptr)
461   {
462     solid = lv->GetSolid();
463   }
464 
465   //---- Dump solid
466   const G4String& solidName = DumpSolid(solid, extraName);
467 
468   //---- Dump material
469   if(mate == nullptr)
470   {
471     mate = lv->GetMaterial();
472   }
473   const G4String& mateName = DumpMaterial(mate);
474 
475   //---- Dump logical volume (solid + material)
476   (*theFile) << ":VOLU " << SubstituteRefl(AddQuotes(lvName)) << " "
477              << SupressRefl(AddQuotes(solidName)) << " " << AddQuotes(mateName)
478              << G4endl;
479 
480   theLogVols[lvName] = lv;
481 
482   return lvName;
483 }
484 
485 // --------------------------------------------------------------------
486 G4String G4tgbGeometryDumper::DumpMaterial(G4Material* mat)
487 {
488   const G4String& mateName = GetObjectName(mat, theMaterials);
489   if(theMaterials.find(mateName) != theMaterials.cend())  // alredy dumped
490   {
491     return mateName;
492   }
493 
494   std::size_t numElements = mat->GetNumberOfElements();
495   G4double density   = mat->GetDensity() / g * cm3;
496 
497   // start tag
498   //
499   if(numElements == 1)
500   {
501     (*theFile) << ":MATE " << AddQuotes(mateName) << " " << mat->GetZ() << " "
502                << mat->GetA() / (g / mole) << " " << density << G4endl;
503   }
504   else
505   {
506     const G4ElementVector* elems = mat->GetElementVector();
507     const G4double* fractions    = mat->GetFractionVector();
508     for(std::size_t ii = 0; ii < numElements; ++ii)
509     {
510       DumpElement(const_cast<G4Element*>((*elems)[ii]));
511     }
512 
513     (*theFile) << ":MIXT " << AddQuotes(mateName) << " " << density << " "
514                << numElements << G4endl;
515     // close start element tag and get ready to do composit "parts"
516     for(std::size_t ii = 0; ii < numElements; ++ii)
517     {
518       (*theFile) << "   " << AddQuotes(GetObjectName(const_cast<G4Element*>((*elems)[ii]), theElements))
519                  << " " << fractions[ii] << G4endl;
520     }
521   }
522 
523   (*theFile) << ":MATE_MEE " << AddQuotes(mateName) << " "
524              << mat->GetIonisation()->GetMeanExcitationEnergy() / eV << "*eV"
525              << G4endl;
526 
527   (*theFile) << ":MATE_TEMPERATURE " << AddQuotes(mateName) << " "
528              << mat->GetTemperature() / kelvin << "*kelvin" << G4endl;
529 
530   (*theFile) << ":MATE_PRESSURE " << AddQuotes(mateName) << " "
531              << mat->GetPressure() / atmosphere << "*atmosphere" << G4endl;
532 
533   G4State state = mat->GetState();
534   G4String stateStr;
535   switch(state)
536   {
537     case kStateUndefined:
538       stateStr = "Undefined";
539       break;
540     case kStateSolid:
541       stateStr = "Solid";
542       break;
543     case kStateLiquid:
544       stateStr = "Liquid";
545       break;
546     case kStateGas:
547       stateStr = "Gas";
548       break;
549   }
550 
551   (*theFile) << ":MATE_STATE " << AddQuotes(mateName) << " " << stateStr
552              << G4endl;
553 
554   theMaterials[mateName] = mat;
555 
556   return mateName;
557 }
558 
559 // --------------------------------------------------------------------
560 void G4tgbGeometryDumper::DumpElement(G4Element* ele)
561 {
562   const G4String& elemName = GetObjectName(ele, theElements);
563 
564   if(theElements.find(elemName) != theElements.cend())  // alredy dumped
565   {
566     return;
567   }
568 
569   //--- Add symbol name: Material mixtures store the components as elements
570   //    (even if the input are materials), but without symbol
571   //
572   G4String symbol = ele->GetSymbol();
573   if(symbol == "" || symbol == " ")
574   {
575     symbol = elemName;
576   }
577 
578   if(ele->GetNumberOfIsotopes() == 0)
579   {
580     (*theFile) << ":ELEM " << AddQuotes(elemName) << " " << AddQuotes(symbol)
581                << " " << ele->GetZ() << " " << ele->GetA() / (g / mole) << " "
582                << G4endl;
583   }
584   else
585   {
586     const G4IsotopeVector* isots = ele->GetIsotopeVector();
587     for(std::size_t ii = 0; ii < ele->GetNumberOfIsotopes(); ++ii)
588     {
589       DumpIsotope((*isots)[ii]);
590     }
591 
592     (*theFile) << ":ELEM_FROM_ISOT " << AddQuotes(elemName) << " "
593                << AddQuotes(symbol) << " " << ele->GetNumberOfIsotopes()
594                << G4endl;
595     const G4double* fractions = ele->GetRelativeAbundanceVector();
596     for(std::size_t ii = 0; ii < ele->GetNumberOfIsotopes(); ++ii)
597     {
598       (*theFile) << "   " << AddQuotes(GetObjectName((*isots)[ii], theIsotopes))
599                  << " " << fractions[ii] << G4endl;
600     }
601   }
602   theElements[elemName] = ele;
603 }
604 
605 // --------------------------------------------------------------------
606 void G4tgbGeometryDumper::DumpIsotope(G4Isotope* isot)
607 {
608   const G4String& isotName = GetObjectName(isot, theIsotopes);
609   if(theIsotopes.find(isotName) != theIsotopes.cend())  // alredy dumped
610   {
611     return;
612   }
613 
614   (*theFile) << ":ISOT " << AddQuotes(isotName) << " " << isot->GetZ() << " "
615              << isot->GetN() << " " << isot->GetA() / (g / mole) << " "
616              << G4endl;
617 
618   theIsotopes[isotName] = isot;
619 }
620 
621 // --------------------------------------------------------------------
622 G4String G4tgbGeometryDumper::DumpSolid(G4VSolid* solid,
623                                         const G4String& extraName)
624 {
625   G4String solidName;
626   if(extraName == "")
627   {
628     solidName = GetObjectName(solid, theSolids);
629   }
630   else
631   {
632     solidName = solid->GetName() + extraName;
633   }
634 
635   if(theSolids.find(solidName) != theSolids.cend())  // alredy dumped
636   {
637     return solidName;
638   }
639 
640   G4String solidType = solid->GetEntityType();
641   solidType          = GetTGSolidType(solidType);
642   
643   if(solidType == "UNIONSOLID")
644   {
645     DumpBooleanVolume("UNION", solid);
646   }
647   else if(solidType == "SUBTRACTIONSOLID")
648   {
649     DumpBooleanVolume("SUBTRACTION", solid);
650   }
651   else if(solidType == "INTERSECTIONSOLID")
652   {
653     DumpBooleanVolume("INTERSECTION", solid);
654   }
655   else if(solidType == "REFLECTEDSOLID")
656   {
657     G4ReflectedSolid* solidrefl = dynamic_cast<G4ReflectedSolid*>(solid);
658     if(solidrefl == nullptr)
659     {
660       G4Exception("G4tgbGeometryDumper::DumpSolid()", "InvalidType",
661                   FatalException, "Invalid reflected solid!");
662       return solidName;
663     }
664     G4VSolid* solidori = solidrefl->GetConstituentMovedSolid();
665     DumpSolid(solidori);
666   }
667   else if(solidType == "MULTIUNION")
668   {
669     DumpMultiUnionVolume(solid);
670   }
671   else if(solidType == "SCALEDSOLID")
672   {
673     DumpScaledVolume(solid);
674   }
675   else
676   {
677     (*theFile) << ":SOLID " << AddQuotes(solidName) << " ";
678     (*theFile) << AddQuotes(solidType) << " ";
679     DumpSolidParams( solid );
680     theSolids[solidName] = solid;
681   }
682 
683   return solidName;
684 }
685 
686 // --------------------------------------------------------------------
687 void G4tgbGeometryDumper::DumpBooleanVolume(const G4String& solidType,
688                                             G4VSolid* so)
689 {
690   G4BooleanSolid* bso = dynamic_cast<G4BooleanSolid*>(so);
691   if(bso == nullptr)
692   {
693     return;
694   }
695   G4VSolid* solid0             = bso->GetConstituentSolid(0);
696   G4VSolid* solid1             = bso->GetConstituentSolid(1);
697   G4DisplacedSolid* solid1Disp = nullptr;
698   G4bool displaced             = dynamic_cast<G4DisplacedSolid*>(solid1);
699   if(displaced)
700   {
701     solid1Disp = dynamic_cast<G4DisplacedSolid*>(solid1);
702     if(solid1Disp != nullptr)
703     {
704       solid1 = solid1Disp->GetConstituentMovedSolid();
705     }
706     else
707     {
708       return;
709     }
710   }
711   DumpSolid(solid0);
712   DumpSolid(solid1);
713 
714   G4String rotName;
715   G4ThreeVector pos;
716   if(displaced)
717   {
718     pos = solid1Disp->GetObjectTranslation();  // translation is of mother frame
719     rotName = DumpRotationMatrix(new G4RotationMatrix(
720       (solid1Disp->GetTransform().NetRotation()).inverse()));
721   }
722   else  // no displacement
723   {
724     rotName = DumpRotationMatrix(new G4RotationMatrix);
725     pos     = G4ThreeVector();
726   }
727 
728   const G4String& bsoName = GetObjectName(so, theSolids);
729   if(theSolids.find(bsoName) != theSolids.cend()) return;  // alredy dumped
730   const G4String& solid0Name = FindSolidName(solid0);
731   const G4String& solid1Name = FindSolidName(solid1);
732 
733   (*theFile) << ":SOLID " << AddQuotes(bsoName) << " " << AddQuotes(solidType)
734              << " " << AddQuotes(solid0Name) << " " << AddQuotes(solid1Name)
735              << " " << AddQuotes(rotName) << " " << approxTo0(pos.x()) << " "
736              << approxTo0(pos.y()) << " " << approxTo0(pos.z()) << " "
737              << G4endl;
738 
739   theSolids[bsoName] = bso;
740 }
741 
742 // --------------------------------------------------------------------
743 void G4tgbGeometryDumper::DumpMultiUnionVolume( G4VSolid* so)
744 {
745   const G4MultiUnion* muun = dynamic_cast<const G4MultiUnion*>(so);
746   if(muun != nullptr)
747     {
748       G4int nSolids = muun->GetNumberOfSolids();
749       std::vector<G4String> rotList;
750       for( G4int iso = 0; iso < nSolids; iso++ ) {
751   G4Transform3D trans = muun->GetTransformation(iso);
752   const G4String& rotName = DumpRotationMatrix( new G4RotationMatrix(trans.getRotation()));
753   rotList.push_back(rotName);
754   G4VSolid* solN = muun->GetSolid(iso);
755   DumpSolid(solN);
756       }
757       const G4String& bsoName = GetObjectName(const_cast<G4VSolid*>(so), theSolids);
758       (*theFile) << ":SOLID " << AddQuotes(bsoName) << " MULTIUNION "
759      << nSolids;
760       
761       for( G4int iso = 0; iso < nSolids; ++iso ) {
762   G4VSolid* solN = muun->GetSolid(iso);
763   G4Transform3D trans = muun->GetTransformation(iso);
764   G4ThreeVector pos = trans.getTranslation();  // translation is of mother frame
765   (*theFile) << " " <<  solN->GetName()
766        << " " << " " << rotList[iso]
767        << " " << approxTo0(pos.x())
768        << " " << approxTo0(pos.y())
769        << " " << approxTo0(pos.z());
770       }
771       (*theFile) << G4endl;
772       
773     }
774 }
775 
776 // --------------------------------------------------------------------
777 void G4tgbGeometryDumper::DumpScaledVolume( G4VSolid* so)
778 {
779   const G4ScaledSolid* ssol = dynamic_cast<const G4ScaledSolid*>(so);
780   if(ssol != nullptr)
781     {
782       G4VSolid* unscaledSolid = ssol->GetUnscaledSolid();
783       G4Scale3D scaleTransf = ssol->GetScaleTransform();
784       G4String bsoName = GetObjectName(const_cast<G4VSolid*>(so), theSolids);
785       (*theFile) << ":SOLID " << AddQuotes(bsoName) << " SCALED "
786      << unscaledSolid->GetName() << " "
787      << scaleTransf.xx() << " "
788      << scaleTransf.yy() << " "
789      << scaleTransf.zz() << G4endl;
790     }
791 }
792 
793 // --------------------------------------------------------------------
794 void G4tgbGeometryDumper::DumpSolidParams(G4VSolid* so)
795 {
796   std::vector<G4double> params = GetSolidParams(so);
797   for(std::size_t ii = 0; ii < params.size(); ++ii)
798   {
799     (*theFile) << params[ii] << " ";
800   }
801   (*theFile) << G4endl;
802 }
803 
804 // --------------------------------------------------------------------
805 std::vector<G4double> G4tgbGeometryDumper::GetSolidParams(const G4VSolid* so)
806 {
807   std::vector<G4double> params;
808 
809   G4String solidType = so->GetEntityType();
810   solidType          = GetTGSolidType(solidType);
811 
812   if(solidType == "BOX")
813   {
814     const G4Box* sb = dynamic_cast<const G4Box*>(so);
815     if(sb != nullptr)
816     {
817       params.push_back(sb->GetXHalfLength());
818       params.push_back(sb->GetYHalfLength());
819       params.push_back(sb->GetZHalfLength());
820     }
821   }
822   else if(solidType == "TUBS")
823   {
824     const G4Tubs* tu = dynamic_cast<const G4Tubs*>(so);
825     if(tu != nullptr)
826     {
827       params.push_back(tu->GetInnerRadius());
828       params.push_back(tu->GetOuterRadius());
829       params.push_back(tu->GetZHalfLength());
830       params.push_back(tu->GetStartPhiAngle() / deg);
831       params.push_back(tu->GetDeltaPhiAngle() / deg);
832     }
833   }
834   else if(solidType == "TRAP")
835   {
836     const G4Trap* trp = dynamic_cast<const G4Trap*>(so);
837     if(trp != nullptr)
838     {
839       G4ThreeVector symAxis(trp->GetSymAxis());
840       params.push_back(trp->GetZHalfLength());
841       params.push_back(symAxis.theta() / deg);
842       params.push_back(symAxis.phi() / deg);
843       params.push_back(trp->GetYHalfLength1());
844       params.push_back(trp->GetXHalfLength1());
845       params.push_back(trp->GetXHalfLength2());
846       params.push_back(std::atan(trp->GetTanAlpha1()) / deg);
847       params.push_back(trp->GetYHalfLength2());
848       params.push_back(trp->GetXHalfLength3());
849       params.push_back(trp->GetXHalfLength4());
850       params.push_back(std::atan(trp->GetTanAlpha2()) / deg);
851     }
852   }
853   else if(solidType == "TRD")
854   {
855     const G4Trd* tr = dynamic_cast<const G4Trd*>(so);
856     if(tr != nullptr)
857     {
858       params.push_back(tr->GetXHalfLength1());
859       params.push_back(tr->GetXHalfLength2());
860       params.push_back(tr->GetYHalfLength1());
861       params.push_back(tr->GetYHalfLength2());
862       params.push_back(tr->GetZHalfLength());
863     }
864   }
865   else if(solidType == "PARA")
866   {
867     const G4Para* para = dynamic_cast<const G4Para*>(so);
868     if(para != nullptr)
869     {
870       G4ThreeVector symAxis(para->GetSymAxis());
871       params.push_back(para->GetXHalfLength());
872       params.push_back(para->GetYHalfLength());
873       params.push_back(para->GetZHalfLength());
874       params.push_back(std::atan(para->GetTanAlpha()) / deg);
875       params.push_back(symAxis.theta() / deg);
876       params.push_back(symAxis.phi() / deg);
877     }
878   }
879   else if(solidType == "CONS")
880   {
881     const G4Cons* cn = dynamic_cast<const G4Cons*>(so);
882     if(cn != nullptr)
883     {
884       params.push_back(cn->GetInnerRadiusMinusZ());
885       params.push_back(cn->GetOuterRadiusMinusZ());
886       params.push_back(cn->GetInnerRadiusPlusZ());
887       params.push_back(cn->GetOuterRadiusPlusZ());
888       params.push_back(cn->GetZHalfLength());
889       params.push_back(cn->GetStartPhiAngle() / deg);
890       params.push_back(cn->GetDeltaPhiAngle() / deg);
891     }
892   }
893   else if(solidType == "SPHERE")
894   {
895     const G4Sphere* sphere = dynamic_cast<const G4Sphere*>(so);
896     if(sphere != nullptr)
897     {
898       params.push_back(sphere->GetInnerRadius());
899       params.push_back(sphere->GetOuterRadius());
900       params.push_back(sphere->GetStartPhiAngle() / deg);
901       params.push_back(sphere->GetDeltaPhiAngle() / deg);
902       params.push_back(sphere->GetStartThetaAngle() / deg);
903       params.push_back(sphere->GetDeltaThetaAngle() / deg);
904     }
905   }
906   else if(solidType == "ORB")
907   {
908     const G4Orb* orb = dynamic_cast<const G4Orb*>(so);
909     if(orb != nullptr)
910     {
911       params.push_back(orb->GetRadius());
912     }
913   }
914   else if(solidType == "TORUS")
915   {
916     const G4Torus* torus = dynamic_cast<const G4Torus*>(so);
917     if(torus != nullptr)
918     {
919       params.push_back(torus->GetRmin());
920       params.push_back(torus->GetRmax());
921       params.push_back(torus->GetRtor());
922       params.push_back(torus->GetSPhi() / deg);
923       params.push_back(torus->GetDPhi() / deg);
924     }
925   }
926   else if(solidType == "POLYCONE")
927   {
928     //--- Dump RZ corners, as original parameters will not be present
929     //    if it was build from RZ corners
930     const G4Polycone* plc = dynamic_cast<const G4Polycone*>(so);
931     if(plc != nullptr)
932     {
933       G4double angphi = plc->GetStartPhi() / deg;
934       if(angphi > 180 * deg)
935       {
936         angphi -= 360 * deg;
937       }
938       G4double endphi = plc->GetEndPhi() / deg;
939       if(endphi > 180 * deg)
940       {
941         endphi -= 360 * deg;
942       }
943       params.push_back(angphi);
944       params.push_back(endphi - angphi);
945       //      params.push_back(plc->GetOriginalParameters()->Opening_angle / deg);
946       G4int ncor = plc->GetNumRZCorner();
947       params.push_back(ncor);
948 
949       for(G4int ii = 0; ii < ncor; ++ii)
950       {
951         params.push_back(plc->GetCorner(ii).r);
952         params.push_back(plc->GetCorner(ii).z);
953       }
954     }
955   }
956   else if(solidType == "GENERICPOLYCONE")
957   {
958     //--- Dump RZ corners
959     const G4GenericPolycone* plc = dynamic_cast<const G4GenericPolycone*>(so);
960     if(plc != nullptr)
961     {
962       G4double angphi = plc->GetStartPhi() / deg;
963       if(angphi > 180 * deg)
964       {
965         angphi -= 360 * deg;
966       }
967       G4double endphi = plc->GetEndPhi() / deg;
968       if(endphi > 180 * deg)
969       {
970         endphi -= 360 * deg;
971       }
972       params.push_back(angphi);
973       params.push_back(endphi - angphi);
974       G4int ncor = plc->GetNumRZCorner();
975       params.push_back(ncor);
976 
977       for(G4int ii = 0; ii < ncor; ++ii)
978       {
979         params.push_back(plc->GetCorner(ii).r);
980         params.push_back(plc->GetCorner(ii).z);
981       }
982     }
983   }
984   else if(solidType == "POLYHEDRA")
985   {
986     //--- Dump RZ corners, as original parameters will not be present
987     //    if it was build from RZ corners
988     const G4Polyhedra* ph = (dynamic_cast<const G4Polyhedra*>(so));
989     if(ph != nullptr)
990     {
991       G4double angphi = ph->GetStartPhi() / deg;
992       if(angphi > 180 * deg)
993         angphi -= 360 * deg;
994 
995       G4int ncor = ph->GetNumRZCorner();
996 
997       params.push_back(angphi);
998       params.push_back(ph->GetOriginalParameters()->Opening_angle / deg);
999       params.push_back(ph->GetNumSide());
1000       params.push_back(ncor);
1001 
1002       for(G4int ii = 0; ii < ncor; ++ii)
1003       {
1004         params.push_back(ph->GetCorner(ii).r);
1005         params.push_back(ph->GetCorner(ii).z);
1006       }
1007     }
1008   }
1009   else if(solidType == "ELLIPTICALTUBE")
1010   {
1011     const G4EllipticalTube* eltu = dynamic_cast<const G4EllipticalTube*>(so);
1012     if(eltu != nullptr)
1013     {
1014       params.push_back(eltu->GetDx());
1015       params.push_back(eltu->GetDy());
1016       params.push_back(eltu->GetDz());
1017     }
1018   }
1019   else if(solidType == "ELLIPSOID")
1020   {
1021     const G4Ellipsoid* dso = dynamic_cast<const G4Ellipsoid*>(so);
1022     if(dso != nullptr)
1023     {
1024       params.push_back(dso->GetSemiAxisMax(0));
1025       params.push_back(dso->GetSemiAxisMax(1));
1026       params.push_back(dso->GetSemiAxisMax(2));
1027       params.push_back(dso->GetZBottomCut());
1028       params.push_back(dso->GetZTopCut());
1029     }
1030   }
1031   else if(solidType == "ELLIPTICAL_CONE")
1032   {
1033     const G4EllipticalCone* elco = dynamic_cast<const G4EllipticalCone*>(so);
1034     if(elco != nullptr)
1035     {
1036       params.push_back(elco->GetSemiAxisX());
1037       params.push_back(elco->GetSemiAxisY());
1038       params.push_back(elco->GetZMax());
1039       params.push_back(elco->GetZTopCut());
1040     }
1041   }
1042   else if(solidType == "HYPE")
1043   {
1044     const G4Hype* hype = dynamic_cast<const G4Hype*>(so);
1045     if(hype != nullptr)
1046     {
1047       params.push_back(hype->GetInnerRadius());
1048       params.push_back(hype->GetOuterRadius());
1049       params.push_back(hype->GetInnerStereo() / deg);
1050       params.push_back(hype->GetOuterStereo() / deg);
1051       params.push_back(2 * hype->GetZHalfLength());
1052     }
1053     //  } else if( solidType == "TET" ) {
1054   }
1055   else if(solidType == "TWISTEDBOX")
1056   {
1057     const G4TwistedBox* tbox = dynamic_cast<const G4TwistedBox*>(so);
1058     if(tbox != nullptr)
1059     {
1060       params.push_back(tbox->GetPhiTwist() / deg);
1061       params.push_back(tbox->GetXHalfLength());
1062       params.push_back(tbox->GetYHalfLength());
1063       params.push_back(tbox->GetZHalfLength());
1064     }
1065   }
1066   else if(solidType == "TWISTEDTRAP")
1067   {
1068     const G4TwistedTrap* ttrap = dynamic_cast<const G4TwistedTrap*>(so);
1069     if(ttrap != nullptr)
1070     {
1071       params.push_back(ttrap->GetPhiTwist() / deg);
1072       params.push_back(ttrap->GetZHalfLength());
1073       params.push_back(ttrap->GetPolarAngleTheta() / deg);
1074       params.push_back(ttrap->GetAzimuthalAnglePhi() / deg);
1075       params.push_back(ttrap->GetY1HalfLength());
1076       params.push_back(ttrap->GetX1HalfLength());
1077       params.push_back(ttrap->GetX2HalfLength());
1078       params.push_back(ttrap->GetY2HalfLength());
1079       params.push_back(ttrap->GetX3HalfLength());
1080       params.push_back(ttrap->GetX4HalfLength());
1081       params.push_back(ttrap->GetTiltAngleAlpha() / deg);
1082     }
1083   }
1084   else if(solidType == "TWISTEDTRD")
1085   {
1086     const G4TwistedTrd* ttrd = dynamic_cast<const G4TwistedTrd*>(so);
1087     if(ttrd != nullptr)
1088     {
1089       params.push_back(ttrd->GetX1HalfLength());
1090       params.push_back(ttrd->GetX2HalfLength());
1091       params.push_back(ttrd->GetY1HalfLength());
1092       params.push_back(ttrd->GetY2HalfLength());
1093       params.push_back(ttrd->GetZHalfLength());
1094       params.push_back(ttrd->GetPhiTwist() / deg);
1095     }
1096   }
1097   else if(solidType == "TWISTEDTUBS")
1098   {
1099     const G4TwistedTubs* ttub = dynamic_cast<const G4TwistedTubs*>(so);
1100     if(ttub != nullptr)
1101     {
1102       params.push_back(ttub->GetInnerRadius());
1103       params.push_back(ttub->GetOuterRadius());
1104       params.push_back(ttub->GetZHalfLength());
1105       params.push_back(ttub->GetDPhi() / deg);
1106       params.push_back(ttub->GetPhiTwist() / deg);
1107     }
1108   }
1109   else
1110   {
1111     const G4String& ErrMessage = "Solid type not supported, sorry... " + solidType;
1112     G4Exception("G4tgbGeometryDumper::DumpSolidParams()", "NotImplemented",
1113                 FatalException, ErrMessage);
1114   }
1115 
1116   return params;
1117 }
1118 
1119 // --------------------------------------------------------------------
1120 G4String G4tgbGeometryDumper::DumpRotationMatrix(G4RotationMatrix* rotm)
1121 {
1122   if(rotm == nullptr)
1123   {
1124     rotm = new G4RotationMatrix();
1125   }
1126 
1127   G4double de      = MatDeterminant(rotm);
1128   G4String rotName = LookForExistingRotation(rotm);
1129   if(rotName != "")
1130   {
1131     return rotName;
1132   }
1133 
1134   G4ThreeVector v(1., 1., 1.);
1135   if(de < -0.9)  // a reflection ....
1136   {
1137     (*theFile) << ":ROTM ";
1138     rotName = "RRM";
1139     rotName += G4UIcommand::ConvertToString(theRotationNumber++);
1140 
1141     (*theFile) << AddQuotes(rotName) << std::setprecision(9) << " "
1142                << approxTo0(rotm->xx()) << " " << approxTo0(rotm->yx()) << " "
1143                << approxTo0(rotm->zx()) << " " << approxTo0(rotm->xy()) << " "
1144                << approxTo0(rotm->yy()) << " " << approxTo0(rotm->zy()) << " "
1145                << approxTo0(rotm->xz()) << " " << approxTo0(rotm->yz()) << " "
1146                << approxTo0(rotm->zz()) << G4endl;
1147   }
1148   else if(de > 0.9)  // a rotation ....
1149   {
1150     (*theFile) << ":ROTM ";
1151     rotName = "RM";
1152     rotName += G4UIcommand::ConvertToString(theRotationNumber++);
1153 
1154     (*theFile) << AddQuotes(rotName) << " " << approxTo0(rotm->thetaX() / deg)
1155                << " " << approxTo0(rotm->phiX() / deg) << " "
1156                << approxTo0(rotm->thetaY() / deg) << " "
1157                << approxTo0(rotm->phiY() / deg) << " "
1158                << approxTo0(rotm->thetaZ() / deg) << " "
1159                << approxTo0(rotm->phiZ() / deg) << G4endl;
1160   }
1161 
1162   theRotMats[rotName] = rotm;
1163 
1164   return rotName;
1165 }
1166 
1167 // --------------------------------------------------------------------
1168 std::vector<G4VPhysicalVolume*>
1169 G4tgbGeometryDumper::GetPVChildren(G4LogicalVolume* lv)
1170 {
1171   G4PhysicalVolumeStore* pvstore = G4PhysicalVolumeStore::GetInstance();
1172   std::vector<G4VPhysicalVolume*> children;
1173   for(auto ite = pvstore->cbegin(); ite != pvstore->cend(); ++ite)
1174   {
1175     if((*ite)->GetMotherLogical() == lv)
1176     {
1177       children.push_back(*ite);
1178 #ifdef G4VERBOSE
1179       if(G4tgrMessenger::GetVerboseLevel() >= 1)
1180       {
1181         G4cout << " G4tgbGeometryDumper::GetPVChildren() - adding children: "
1182                << (*ite)->GetName() << " of " << lv->GetName() << G4endl;
1183       }
1184 #endif
1185     }
1186   }
1187 
1188   return children;
1189 }
1190 
1191 // --------------------------------------------------------------------
1192 G4String G4tgbGeometryDumper::GetTGSolidType(const G4String& solidType)
1193 {
1194   G4String newsolidType = solidType.substr(2, solidType.length());
1195   for(G4int ii = 0; ii < (G4int)newsolidType.length(); ++ii)
1196   {
1197     newsolidType[ii] = (char)std::toupper(newsolidType[ii]);
1198   }
1199   return newsolidType;
1200 }
1201 
1202 // --------------------------------------------------------------------
1203 G4double G4tgbGeometryDumper::MatDeterminant(G4RotationMatrix* ro)
1204 {
1205   G4Rep3x3 r = ro->rep3x3();
1206   return r.xx_ * (r.yy_ * r.zz_ - r.zy_ * r.yz_) -
1207          r.yx_ * (r.xy_ * r.zz_ - r.zy_ * r.xz_) +
1208          r.zx_ * (r.xy_ * r.yz_ - r.yy_ * r.xz_);
1209 }
1210 
1211 // --------------------------------------------------------------------
1212 G4double G4tgbGeometryDumper::approxTo0(G4double val)
1213 {
1214   G4double precision =
1215     G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
1216 
1217   if(std::fabs(val) < precision)
1218   {
1219     val = 0.0;
1220   }
1221   return val;
1222 }
1223 
1224 // --------------------------------------------------------------------
1225 G4String G4tgbGeometryDumper::AddQuotes(const G4String& str)
1226 {
1227   //--- look if there is a separating blank
1228 
1229   G4bool bBlank = FALSE;
1230   std::size_t siz = str.length();
1231   for(std::size_t ii = 0; ii < siz; ++ii)
1232   {
1233     if(str.substr(ii, 1) == " ")
1234     {
1235       bBlank = TRUE;
1236       break;
1237     }
1238   }
1239   G4String str2 = str;
1240   if(bBlank)
1241   {
1242     str2 = G4String("\"") + str2 + G4String("\"");
1243   }
1244   return str2;
1245 }
1246 
1247 // --------------------------------------------------------------------
1248 G4String G4tgbGeometryDumper::SupressRefl(G4String name)
1249 {
1250   std::size_t irefl = name.rfind("_refl");
1251   if(irefl != G4String::npos)
1252   {
1253     name = name.substr(0, irefl);
1254   }
1255   return name;
1256 }
1257 
1258 // --------------------------------------------------------------------
1259 G4String G4tgbGeometryDumper::SubstituteRefl(G4String name)
1260 {
1261   std::size_t irefl = name.rfind("_refl");
1262   if(irefl != G4String::npos)
1263   {
1264     name = name.substr(0, irefl) + "_REFL";
1265   }
1266   return name;
1267 }
1268 
1269 // --------------------------------------------------------------------
1270 G4String G4tgbGeometryDumper::GetIsotopeName(G4Isotope* isot)
1271 {
1272   G4String isotName = isot->GetName();
1273   // first look if this is isotope is already dumped,
1274   // with original isotope name or new one
1275   //
1276   std::map<G4String, G4Isotope*>::const_iterator ite;
1277   for(ite = theIsotopes.cbegin(); ite != theIsotopes.cend(); ++ite)
1278   {
1279     if(isot == (*ite).second)
1280     {
1281       return (*ite).first;
1282     }
1283   }
1284 
1285   // Now look if there is another isotope dumped with same name,
1286   // and if found add _N to the name
1287   //
1288   ite = theIsotopes.find(isotName);
1289   if(ite != theIsotopes.cend())  // Isotope found with same name
1290   {
1291     G4Isotope* isotold = (*ite).second;
1292     if(isot != isotold)  // new isotope it is not the really
1293     {                    // the same one as isotope found
1294       if(!Same2G4Isotopes(isot, isotold))
1295       {                // if the two have same data, use the old one
1296         G4int ii = 2;  // G4Nist does names isotopes of same element
1297                        // with same name
1298         for(;; ++ii)
1299         {
1300           const G4String& newIsotName =
1301             isotName + "_" + G4UIcommand::ConvertToString(ii);
1302           std::map<G4String, G4Isotope*>::const_iterator ite2 =
1303             theIsotopes.find(newIsotName);
1304           if(ite2 == theIsotopes.cend())
1305           {
1306             isotName = newIsotName;
1307             break;
1308           }
1309           else
1310           {
1311             if(Same2G4Isotopes(isot, (*ite2).second))
1312             {
1313               isotName = newIsotName;
1314               break;
1315             }
1316           }
1317         }
1318       }
1319     }
1320   }
1321   return isotName;
1322 }
1323 
1324 // --------------------------------------------------------------------
1325 template <class TYP>
1326 G4String G4tgbGeometryDumper::GetObjectName(
1327   TYP* obj, std::map<G4String, TYP*> objectsDumped)
1328 {
1329   G4String objName = obj->GetName();
1330 
1331   // first look if this is objecy is already dumped,
1332   // with original object name or new one
1333   //
1334   typename std::map<G4String, TYP*>::const_iterator ite;
1335   for(ite = objectsDumped.cbegin(); ite != objectsDumped.cend(); ++ite)
1336   {
1337     if(obj == (*ite).second)
1338     {
1339       return (*ite).first;
1340     }
1341   }
1342 
1343   // Now look if there is another object dumped with same name,
1344   // and if found add _N to the name
1345   //
1346   ite = objectsDumped.find(objName);
1347 
1348   if(ite != objectsDumped.cend())  // Object found with same name
1349   {
1350     TYP* objold = (*ite).second;
1351     if(obj != objold)  // new object it is not the really
1352     {                  // the same one as object found
1353       G4int ii = 2;
1354       for(;; ++ii)
1355       {
1356         const G4String& newObjName = objName + "_" + G4UIcommand::ConvertToString(ii);
1357         typename std::map<G4String, TYP*>::const_iterator ite2 =
1358           objectsDumped.find(newObjName);
1359         if(ite2 == objectsDumped.cend())
1360         {
1361           objName = newObjName;
1362           break;
1363         }
1364       }
1365     }
1366   }
1367   return objName;
1368 }
1369 
1370 // --------------------------------------------------------------------
1371 G4bool G4tgbGeometryDumper::CheckIfLogVolExists(const G4String& name,
1372                                                 G4LogicalVolume* pt)
1373 {
1374   if(theLogVols.find(name) != theLogVols.cend())
1375   {
1376     G4LogicalVolume* lvnew = (*(theLogVols.find(name))).second;
1377     if(lvnew != pt)
1378     {
1379       /*
1380       //---- Reflected volumes are repeated
1381 
1382       G4ReflectionFactory* reffact = G4ReflectionFactory::Instance();
1383       if( !reffact->IsReflected( pt ) && !reffact->IsReflected( lvnew ) )
1384       {
1385         G4String ErrMessage = "LogVol found but not same as before: " + name;
1386         G4Exception("G4tgbGeometryDumper::CheckIfLogVolExists()",
1387                     "InvalidSetup", FatalException, ErrMessage);
1388       }
1389       */
1390     }
1391     return true;
1392   }
1393   else
1394   {
1395     return false;
1396   }
1397 }
1398 
1399 // --------------------------------------------------------------------
1400 G4bool G4tgbGeometryDumper::CheckIfPhysVolExists(const G4String& name,
1401                                                  G4VPhysicalVolume* pt)
1402 {
1403 #ifdef G4VERBOSE
1404   if(G4tgrMessenger::GetVerboseLevel() >= 1)
1405   {
1406     G4cout << " G4tgbGeometryDumper::CheckIfPhysVolExists() - " << name
1407            << G4endl;
1408   }
1409 #endif
1410   if(thePhysVols.find(name) != thePhysVols.cend())
1411   {
1412     if((*(thePhysVols.find(name))).second != pt)
1413     {
1414       // G4String ErrMessage = "Placement found but not same as before: "
1415       //                     + name;
1416       // G4Exception("G4tgbGeometryDumper::CheckIfPhysVolExists()",
1417       //             "InvalidSetup", FatalException, ErrMessage);
1418       G4cerr << " G4tgbGeometryDumper::CheckIfPhysVolExists () -"
1419              << " Placement found but not same as before : " << name << G4endl;
1420     }
1421     return true;
1422   }
1423   else
1424   {
1425     return false;
1426   }
1427 }
1428 
1429 // --------------------------------------------------------------------
1430 G4String
1431 G4tgbGeometryDumper::LookForExistingRotation(const G4RotationMatrix* rotm)
1432 {
1433   G4String rmName = "";
1434 
1435   for(auto ite = theRotMats.cbegin(); ite != theRotMats.cend(); ++ite)
1436   {
1437     if((*ite).second->isNear(*rotm))
1438     {
1439       rmName = (*ite).first;
1440       break;
1441     }
1442   }
1443   return rmName;
1444 }
1445 
1446 // --------------------------------------------------------------------
1447 G4bool G4tgbGeometryDumper::Same2G4Isotopes(G4Isotope* isot1, G4Isotope* isot2)
1448 {
1449   if((isot1->GetZ() != isot2->GetZ()) || (isot1->GetN() != isot2->GetN()) ||
1450      (isot1->GetA() != isot2->GetA()))
1451   {
1452     return false;
1453   }
1454   else
1455   {
1456     return true;
1457   }
1458 }
1459 
1460 // --------------------------------------------------------------------
1461 const G4String& G4tgbGeometryDumper::FindSolidName(G4VSolid* solid)
1462 {
1463   std::map<G4String, G4VSolid*>::const_iterator ite;
1464   for(ite = theSolids.cbegin(); ite != theSolids.cend(); ++ite)
1465   {
1466     if(solid == (*ite).second)
1467     {
1468       return (*ite).first;
1469     }
1470   }
1471 
1472   if(ite == theSolids.cend())
1473   {
1474     G4Exception("G4tgbGeometryDumper::FindSolidName()", "ReadError",
1475                 FatalException, "Programming error.");
1476   }
1477   return (*ite).first;
1478 }
1479