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