source: trunk/source/persistency/ascii/src/G4tgbGeometryDumper.cc@ 1357

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

geant4 tag 9.4

File size: 40.6 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27// $Id: G4tgbGeometryDumper.cc,v 1.15 2010/11/02 11:13:05 gcosmo Exp $
28// GEANT4 tag $Name: geant4-09-04-ref-00 $
29//
30//
31// class G4tgbGeometryDumper
32
33// History:
34// - Created. P.Arce, CIEMAT (November 2007)
35// -------------------------------------------------------------------------
36
37#include "G4tgbGeometryDumper.hh"
38
39#include "G4tgrMessenger.hh"
40
41#include "G4UIcommand.hh"
42#include "G4Material.hh"
43#include "G4Element.hh"
44#include "G4VSolid.hh"
45#include "G4Box.hh"
46#include "G4Tubs.hh"
47#include "G4Cons.hh"
48#include "G4Trap.hh"
49#include "G4Sphere.hh"
50#include "G4Orb.hh"
51#include "G4Trd.hh"
52#include "G4Para.hh"
53#include "G4Torus.hh"
54#include "G4Hype.hh"
55#include "G4Polycone.hh"
56#include "G4Polyhedra.hh"
57#include "G4EllipticalTube.hh"
58#include "G4Ellipsoid.hh"
59#include "G4EllipticalCone.hh"
60#include "G4Hype.hh"
61#include "G4Tet.hh"
62#include "G4TwistedBox.hh"
63#include "G4TwistedTrap.hh"
64#include "G4TwistedTrd.hh"
65#include "G4TwistedTubs.hh"
66#include "G4PVPlacement.hh"
67#include "G4PVParameterised.hh"
68#include "G4PVReplica.hh"
69#include "G4BooleanSolid.hh"
70#include "G4ReflectionFactory.hh"
71#include "G4ReflectedSolid.hh"
72#include "G4LogicalVolumeStore.hh"
73#include "G4PhysicalVolumeStore.hh"
74#include "G4GeometryTolerance.hh"
75#include "G4VPVParameterisation.hh"
76#include <iomanip>
77
78//------------------------------------------------------------------------
79G4tgbGeometryDumper* G4tgbGeometryDumper::theInstance = 0;
80
81//------------------------------------------------------------------------
82G4tgbGeometryDumper::G4tgbGeometryDumper()
83 : theFile(0), theRotationNumber(0)
84{
85}
86
87//------------------------------------------------------------------------
88G4tgbGeometryDumper* G4tgbGeometryDumper::GetInstance()
89{
90 if( theInstance == 0 ){
91 theInstance = new G4tgbGeometryDumper;
92 }
93
94 return theInstance;
95
96}
97
98//------------------------------------------------------------------------
99void G4tgbGeometryDumper::DumpGeometry( const G4String& fname )
100{
101 theFile = new std::ofstream(fname);
102
103 G4VPhysicalVolume* pv = GetTopPhysVol();
104 DumpPhysVol( pv ); // dump volume and recursively it will dump all hierarchy
105}
106
107//---------------------------------------------------------------------
108G4VPhysicalVolume* G4tgbGeometryDumper::GetTopPhysVol()
109{
110 G4PhysicalVolumeStore* pvstore = G4PhysicalVolumeStore::GetInstance();
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//------------------------------------------------------------------------
134G4tgbGeometryDumper::~G4tgbGeometryDumper()
135{
136}
137
138//------------------------------------------------------------------------
139void G4tgbGeometryDumper::DumpPhysVol( G4VPhysicalVolume* pv )
140{
141
142 //--- Dump logical volume first
143 G4LogicalVolume* lv = pv->GetLogicalVolume();
144
145 G4ReflectionFactory* reffact = G4ReflectionFactory::Instance();
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//------------------------------------------------------------------------
204void
205G4tgbGeometryDumper::DumpPVPlacement( G4VPhysicalVolume* pv,
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
214 G4ReflectionFactory* reffact = G4ReflectionFactory::Instance();
215 G4LogicalVolume* lv = pv->GetLogicalVolume();
216 if( reffact->IsReflected( lv ) )
217 {
218#ifdef G4VERBOSE
219 if( G4tgrMessenger::GetVerboseLevel() >= 1 )
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 CLHEP::HepRep3x3 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 );
240 G4ThreeVector pos = pv->GetTranslation();
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 << " "
257 << SubstituteRefl(AddQuotes(pv->GetMotherLogical()->GetName()))
258 << " " << AddQuotes(rotName) << " "
259 << pos.x() << " " << pos.y() << " " << pos.z() << G4endl;
260
261 thePhysVols[fullname] = pv;
262 }
263}
264
265
266//------------------------------------------------------------------------
267void G4tgbGeometryDumper::DumpPVParameterised( G4PVParameterised* pv )
268{
269 G4String pvName = pv->GetName();
270
271 EAxis axis;
272 G4int nReplicas;
273 G4double width;
274 G4double offset;
275 G4bool consuming;
276 pv->GetReplicationData(axis, nReplicas, width, offset, consuming);
277
278 G4VPVParameterisation* param = pv->GetParameterisation();
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//------------------------------------------------------------------------
393void G4tgbGeometryDumper::DumpPVReplica( G4PVReplica* pv,
394 const G4String& lvName )
395{
396 EAxis axis;
397 G4int nReplicas;
398 G4double width;
399 G4double offset;
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//------------------------------------------------------------------------
459G4String
460G4tgbGeometryDumper::DumpLogVol( G4LogicalVolume* lv, G4String extraName,
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//------------------------------------------------------------------------
500G4String G4tgbGeometryDumper::DumpMaterial( G4Material* mat )
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) << " "
542 << mat->GetIonisation()->GetMeanExcitationEnergy()/eV
543 << "*eV" << G4endl;
544
545 theMaterials[mateName] = mat;
546
547 return mateName;
548}
549
550
551//------------------------------------------------------------------------
552void G4tgbGeometryDumper::DumpElement( G4Element* ele)
553{
554 G4String elemName = GetObjectName(ele,theElements);
555
556 if( theElements.find( elemName ) != theElements.end() ) // alredy dumped
557 {
558 return;
559 }
560
561 //--- Add symbol name: Material mixtures store the components as elements
562 // (even if the input are materials), but without symbol
563 //
564 G4String symbol = ele->GetSymbol();
565 if( symbol == "" || symbol == " " )
566 {
567 symbol = elemName;
568 }
569
570 if( ele->GetNumberOfIsotopes() == 0 )
571 {
572 (*theFile) << ":ELEM " << AddQuotes(elemName) << " "
573 << AddQuotes(symbol) << " " << ele->GetZ() << " "
574 << ele->GetA()/(g/mole) << " " << G4endl;
575 }
576 else
577 {
578 const G4IsotopeVector* isots = ele->GetIsotopeVector();
579 for (size_t ii = 0; ii < ele->GetNumberOfIsotopes(); ii++)
580 {
581 DumpIsotope( (*isots)[ii] );
582 }
583
584 (*theFile) << ":ELEM_FROM_ISOT " << AddQuotes(elemName) << " "
585 << AddQuotes(symbol) << " " << ele->GetNumberOfIsotopes()
586 << G4endl;
587 const G4double* fractions = ele->GetRelativeAbundanceVector();
588 for (size_t ii = 0; ii < ele->GetNumberOfIsotopes(); ii++)
589 {
590 (*theFile) << " "
591 << AddQuotes(GetObjectName((*isots)[ii],theIsotopes)) << " "
592 << fractions[ii] << G4endl;
593 }
594 }
595 theElements[elemName] = ele;
596}
597
598
599//------------------------------------------------------------------------
600void G4tgbGeometryDumper::DumpIsotope( G4Isotope* isot)
601{
602 G4String isotName = GetObjectName(isot,theIsotopes);
603 if( theIsotopes.find( isotName ) != theIsotopes.end() ) // alredy dumped
604 {
605 return;
606 }
607
608 (*theFile) << ":ISOT " << AddQuotes(isotName) << " "
609 << isot->GetZ() << " " << isot->GetN() << " "
610 << isot->GetA()/(g/mole) << " " << G4endl;
611
612 theIsotopes[isotName] = isot;
613}
614
615
616//------------------------------------------------------------------------
617G4String G4tgbGeometryDumper::DumpSolid( G4VSolid* solid,
618 const G4String& extraName )
619{
620 G4String solidName;
621 if( extraName == "" )
622 {
623 solidName = GetObjectName(solid,theSolids);
624 }
625 else
626 {
627 solidName = solid->GetName()+extraName;
628 }
629
630 if( theSolids.find( solidName ) != theSolids.end() ) // alredy dumped
631 {
632 return solidName;
633 }
634
635 G4String solidType = solid->GetEntityType();
636 solidType = GetTGSolidType( solidType );
637
638 if (solidType == "UNIONSOLID")
639 {
640 DumpBooleanVolume( "UNION", solid );
641
642 } else if (solidType == "SUBTRACTIONSOLID") {
643 DumpBooleanVolume( "SUBTRACTION", solid );
644
645 } else if (solidType == "INTERSECTIONSOLID") {
646 DumpBooleanVolume( "INTERSECTION", solid );
647
648 } else if (solidType == "REFLECTEDSOLID") {
649 G4ReflectedSolid* solidrefl = dynamic_cast<G4ReflectedSolid*>(solid);
650 if (!solidrefl)
651 {
652 G4Exception("G4tgbGeometryDumper::DumpSolid()",
653 "InvalidType", FatalException, "Invalid reflected solid!");
654 return solidName;
655 }
656 G4VSolid* solidori = solidrefl->GetConstituentMovedSolid();
657 DumpSolid( solidori );
658 }
659 else
660 {
661 (*theFile) << ":SOLID " << AddQuotes(solidName) << " ";
662 (*theFile) << AddQuotes(solidType) << " ";
663 DumpSolidParams( solid );
664 theSolids[solidName] = solid;
665 }
666
667 return solidName;
668}
669
670
671//------------------------------------------------------------------------
672void G4tgbGeometryDumper::DumpBooleanVolume( const G4String& solidType,
673 G4VSolid* so )
674{
675 G4BooleanSolid * bso = dynamic_cast < G4BooleanSolid * > (so);
676 if (!bso) { return; }
677 G4VSolid* solid0 = bso->GetConstituentSolid( 0 );
678 G4VSolid* solid1 = bso->GetConstituentSolid( 1 );
679 G4DisplacedSolid* solid1Disp = 0;
680 G4bool displaced = dynamic_cast<G4DisplacedSolid*>(solid1);
681 if( displaced )
682 {
683 solid1Disp = dynamic_cast<G4DisplacedSolid*>(solid1);
684 if (solid1Disp) { solid1 = solid1Disp->GetConstituentMovedSolid(); }
685 }
686 DumpSolid( solid0 );
687 DumpSolid( solid1 );
688
689 G4String rotName;
690 G4ThreeVector pos;
691 if( displaced )
692 {
693 pos = solid1Disp->GetObjectTranslation(); // translation is of mother frame
694 rotName = DumpRotationMatrix( new G4RotationMatrix( (solid1Disp->
695 GetTransform().NetRotation()).inverse() ) );
696 }
697 else // no displacement
698 {
699 rotName = DumpRotationMatrix( new G4RotationMatrix );
700 pos = G4ThreeVector();
701 }
702
703 G4String bsoName = GetObjectName(so,theSolids);
704 if( theSolids.find( bsoName ) != theSolids.end() ) return; // alredy dumped
705 G4String solid0Name = FindSolidName( solid0 );
706 G4String solid1Name = FindSolidName( solid1 );
707
708 (*theFile) << ":SOLID "
709 << AddQuotes(bsoName) << " "
710 << AddQuotes(solidType) << " "
711 << AddQuotes(solid0Name) << " "
712 << AddQuotes(solid1Name) << " "
713 << AddQuotes(rotName) << " "
714 << approxTo0(pos.x()) << " "
715 << approxTo0(pos.y()) << " "
716 << approxTo0(pos.z()) << " " << G4endl;
717
718 theSolids[bsoName] = bso;
719}
720
721
722//------------------------------------------------------------------------
723void G4tgbGeometryDumper::DumpSolidParams( G4VSolid * so)
724{
725 std::vector<G4double> params = GetSolidParams( so );
726 for( size_t ii = 0 ; ii < params.size(); ii++ )
727 {
728 (*theFile) << params[ii] << " " ;
729 }
730 (*theFile) << G4endl;
731}
732
733
734//------------------------------------------------------------------------
735std::vector<G4double> G4tgbGeometryDumper::GetSolidParams( const G4VSolid * so)
736{
737 std::vector<G4double> params;
738
739 G4String solidType = so->GetEntityType();
740 solidType = GetTGSolidType( solidType );
741
742 if (solidType == "BOX") {
743 const G4Box * sb = dynamic_cast < const G4Box*>(so);
744 if (sb) {
745 params.push_back( sb->GetXHalfLength() );
746 params.push_back( sb->GetYHalfLength() );
747 params.push_back( sb->GetZHalfLength() );
748 }
749 } else if (solidType == "TUBS") {
750 const G4Tubs * tu = dynamic_cast < const G4Tubs * > (so);
751 if (tu) {
752 params.push_back( tu->GetInnerRadius() );
753 params.push_back( tu->GetOuterRadius() );
754 params.push_back( tu->GetZHalfLength() );
755 params.push_back( tu->GetStartPhiAngle()/deg );
756 params.push_back( tu->GetDeltaPhiAngle()/deg );
757 }
758 } else if (solidType == "TRAP") {
759 const G4Trap * trp = dynamic_cast < const G4Trap * > (so);
760 if (trp) {
761 G4ThreeVector symAxis(trp->GetSymAxis());
762 G4double theta = symAxis.theta()/deg;
763 G4double phi = symAxis.phi()/deg;
764 params.push_back( trp->GetZHalfLength() );
765 params.push_back( theta );
766 params.push_back( phi);
767 params.push_back( trp->GetYHalfLength1() );
768 params.push_back( trp->GetXHalfLength1() );
769 params.push_back( trp->GetXHalfLength2() );
770 params.push_back( std::atan(trp->GetTanAlpha1())/deg );
771 params.push_back( trp->GetYHalfLength2() );
772 params.push_back( trp->GetXHalfLength3() );
773 params.push_back( trp->GetXHalfLength4() );
774 params.push_back( std::atan(trp->GetTanAlpha2())/deg );
775 }
776 } else if (solidType == "TRD") {
777 const G4Trd * tr = dynamic_cast < const G4Trd * > (so);
778 if (tr) {
779 params.push_back( tr->GetXHalfLength1() );
780 params.push_back( tr->GetXHalfLength2() );
781 params.push_back( tr->GetYHalfLength1() );
782 params.push_back( tr->GetYHalfLength2() );
783 params.push_back( tr->GetZHalfLength());
784 }
785 } else if (solidType == "PARA") {
786 const G4Para * para = dynamic_cast < const G4Para * > (so);
787 if (para) {
788 G4double phi = 0.;
789 if(para->GetSymAxis().z()!=1.0)
790 { phi = std::atan(para->GetSymAxis().y()/para->GetSymAxis().x()); }
791 params.push_back( para->GetXHalfLength());
792 params.push_back( para->GetYHalfLength());
793 params.push_back( para->GetZHalfLength());
794 params.push_back( std::atan(para->GetTanAlpha())/deg);
795 params.push_back( std::acos(para->GetSymAxis().z())/deg);
796 params.push_back( phi/deg);
797 }
798 } else if (solidType == "CONS") {
799 const G4Cons * cn = dynamic_cast < const G4Cons * > (so);
800 if (cn) {
801 params.push_back( cn->GetInnerRadiusMinusZ() );
802 params.push_back( cn->GetOuterRadiusMinusZ() );
803 params.push_back( cn->GetInnerRadiusPlusZ() );
804 params.push_back( cn->GetOuterRadiusPlusZ() );
805 params.push_back( cn->GetZHalfLength() );
806 params.push_back( cn->GetStartPhiAngle()/deg );
807 params.push_back( cn->GetDeltaPhiAngle()/deg );
808 }
809 } else if (solidType == "SPHERE") {
810 const G4Sphere * sphere = dynamic_cast < const G4Sphere * > (so);
811 if (sphere) {
812 params.push_back( sphere->GetInnerRadius());
813 params.push_back( sphere->GetOuterRadius());
814 params.push_back( sphere->GetStartPhiAngle()/deg);
815 params.push_back( sphere->GetDeltaPhiAngle()/deg);
816 params.push_back( sphere->GetStartThetaAngle()/deg);
817 params.push_back( sphere->GetDeltaThetaAngle()/deg);
818 }
819 } else if (solidType == "ORB") {
820 const G4Orb * orb = dynamic_cast < const G4Orb * > (so);
821 if (orb) {
822 params.push_back( orb->GetRadius());
823 }
824 } else if (solidType == "TORUS") {
825 const G4Torus * torus = dynamic_cast < const G4Torus * > (so);
826 if (torus) {
827 params.push_back( torus->GetRmin());
828 params.push_back( torus->GetRmax());
829 params.push_back( torus->GetRtor());
830 params.push_back( torus->GetSPhi()/deg);
831 params.push_back( torus->GetDPhi()/deg);
832 }
833 } else if (solidType == "POLYCONE") {
834 //--- Dump RZ corners, as original parameters will not be present
835 // if it was build from RZ corners
836 const G4Polycone * pc = dynamic_cast < const G4Polycone * > (so);
837 if (pc) {
838 G4double angphi = pc->GetStartPhi()/deg;
839 if( angphi > 180*deg ) { angphi -= 360*deg; }
840 G4int ncor = pc->GetNumRZCorner();
841 params.push_back( angphi );
842 params.push_back( pc->GetOriginalParameters()->Opening_angle/deg );
843 params.push_back( ncor );
844
845 for( G4int ii = 0; ii < ncor; ii++ )
846 {
847 params.push_back( pc->GetCorner(ii).r );
848 params.push_back( pc->GetCorner(ii).z );
849 }
850 }
851 } else if (solidType == "POLYHEDRA") {
852 //--- Dump RZ corners, as original parameters will not be present
853 // if it was build from RZ corners
854 const G4Polyhedra * ph = (dynamic_cast < const G4Polyhedra * > (so));
855 if (ph) {
856 G4double angphi = ph->GetStartPhi()/deg;
857 if( angphi > 180*deg ) angphi -= 360*deg;
858
859 G4int ncor = ph->GetNumRZCorner();
860
861 params.push_back( angphi );
862 params.push_back( ph->GetOriginalParameters()->Opening_angle/deg );
863 params.push_back( ph->GetNumSide() );
864 params.push_back( ncor );
865
866 for( G4int ii = 0; ii < ncor; ii++ )
867 {
868 params.push_back( ph->GetCorner(ii).r );
869 params.push_back( ph->GetCorner(ii).z );
870 }
871 }
872 } else if (solidType == "ELLIPTICALTUBE") {
873 const G4EllipticalTube * eltu =
874 dynamic_cast < const G4EllipticalTube * > (so);
875 if (eltu) {
876 params.push_back( eltu->GetDx());
877 params.push_back( eltu->GetDy());
878 params.push_back( eltu->GetDz());
879 }
880 } else if (solidType == "ELLIPSOID" ){
881 const G4Ellipsoid* dso = dynamic_cast < const G4Ellipsoid * > (so);
882 if (dso) {
883 params.push_back( dso->GetSemiAxisMax(0) );
884 params.push_back( dso->GetSemiAxisMax(1) );
885 params.push_back( dso->GetSemiAxisMax(2) );
886 params.push_back( dso->GetZBottomCut() );
887 params.push_back( dso->GetZTopCut() );
888 }
889 } else if (solidType == "ELLIPTICAL_CONE") {
890 const G4EllipticalCone * elco =
891 dynamic_cast < const G4EllipticalCone * > (so);
892 if (elco) {
893 params.push_back( elco-> GetSemiAxisX() );
894 params.push_back( elco-> GetSemiAxisY() );
895 params.push_back( elco-> GetZMax() );
896 params.push_back( elco-> GetZTopCut() );
897 }
898 } else if (solidType == "HYPE") {
899 const G4Hype* hype = dynamic_cast < const G4Hype * > (so);
900 if (hype) {
901 params.push_back( hype->GetInnerRadius());
902 params.push_back( hype->GetOuterRadius());
903 params.push_back( hype->GetInnerStereo()/deg);
904 params.push_back( hype->GetOuterStereo()/deg);
905 params.push_back( 2*hype->GetZHalfLength());
906 }
907// } else if( solidType == "TET" ) {
908
909 } else if( solidType == "TWISTEDBOX" ) {
910 const G4TwistedBox* tbox = dynamic_cast < const G4TwistedBox * > (so);
911 if (tbox) {
912 params.push_back( tbox->GetPhiTwist()/deg );
913 params.push_back( tbox->GetXHalfLength() );
914 params.push_back( tbox->GetYHalfLength() );
915 params.push_back( tbox->GetZHalfLength() );
916 }
917 } else if( solidType == "TWISTEDTRAP" ) {
918 const G4TwistedTrap * ttrap = dynamic_cast < const G4TwistedTrap * > (so);
919 if (ttrap) {
920 params.push_back( ttrap->GetPhiTwist()/deg );
921 params.push_back( ttrap->GetZHalfLength() );
922 params.push_back( ttrap->GetPolarAngleTheta()/deg );
923 params.push_back( ttrap->GetAzimuthalAnglePhi()/deg );
924 params.push_back( ttrap->GetY1HalfLength() );
925 params.push_back( ttrap->GetX1HalfLength() );
926 params.push_back( ttrap->GetX2HalfLength() );
927 params.push_back( ttrap->GetY2HalfLength() );
928 params.push_back( ttrap->GetX3HalfLength() );
929 params.push_back( ttrap->GetX4HalfLength() );
930 params.push_back( ttrap->GetTiltAngleAlpha()/deg );
931 }
932 } else if( solidType == "TWISTEDTRD" ) {
933 const G4TwistedTrd * ttrd = dynamic_cast < const G4TwistedTrd * > (so);
934 if (ttrd) {
935 params.push_back( ttrd->GetX1HalfLength());
936 params.push_back( ttrd->GetX2HalfLength() );
937 params.push_back( ttrd->GetY1HalfLength() );
938 params.push_back( ttrd->GetY2HalfLength() );
939 params.push_back( ttrd->GetZHalfLength() );
940 params.push_back( ttrd->GetPhiTwist()/deg );
941 }
942 } else if( solidType == "TWISTEDTUBS" ) {
943 const G4TwistedTubs * ttub = dynamic_cast < const G4TwistedTubs * > (so);
944 if (ttub) {
945 params.push_back( ttub->GetInnerRadius() );
946 params.push_back( ttub->GetOuterRadius() );
947 params.push_back( ttub->GetZHalfLength() );
948 params.push_back( ttub->GetDPhi()/deg );
949 params.push_back( ttub->GetPhiTwist()/deg );
950 }
951 }
952 else
953 {
954 G4String ErrMessage = "Solid type not supported, sorry... " + solidType;
955 G4Exception("G4tgbGeometryDumpe::DumpSolidParams()",
956 "NotImplemented", FatalException, ErrMessage);
957 }
958
959 return params;
960}
961
962
963//------------------------------------------------------------------------
964G4String G4tgbGeometryDumper::DumpRotationMatrix( G4RotationMatrix* rotm )
965{
966 if (!rotm) { rotm = new G4RotationMatrix(); }
967
968 G4double de = MatDeterminant(rotm);
969 G4String rotName = LookForExistingRotation( rotm );
970 if( rotName != "" ) { return rotName; }
971
972 G4ThreeVector v(1.,1.,1.);
973 if (de < -0.9 ) // a reflection ....
974 {
975 (*theFile) << ":ROTM ";
976 rotName = "RRM";
977 rotName += G4UIcommand::ConvertToString(theRotationNumber++);
978
979 (*theFile) << AddQuotes(rotName) << std::setprecision(9) << " "
980 << approxTo0(rotm->xx()) << " "
981 << approxTo0(rotm->yx()) << " "
982 << approxTo0(rotm->zx()) << " "
983 << approxTo0(rotm->xy()) << " "
984 << approxTo0(rotm->yy()) << " "
985 << approxTo0(rotm->zy()) << " "
986 << approxTo0(rotm->xz()) << " "
987 << approxTo0(rotm->yz()) << " "
988 << approxTo0(rotm->zz()) << G4endl;
989 }
990 else if(de > 0.9 ) // a rotation ....
991 {
992 (*theFile) << ":ROTM ";
993 rotName = "RM";
994 rotName += G4UIcommand::ConvertToString(theRotationNumber++);
995
996 (*theFile) << AddQuotes(rotName) << " "
997 << approxTo0(rotm->thetaX()/deg) << " "
998 << approxTo0(rotm->phiX()/deg) << " "
999 << approxTo0(rotm->thetaY()/deg) << " "
1000 << approxTo0(rotm->phiY()/deg) << " "
1001 << approxTo0(rotm->thetaZ()/deg) << " "
1002 << approxTo0(rotm->phiZ()/deg) << G4endl;
1003 }
1004
1005 theRotMats[rotName] = rotm;
1006
1007 return rotName;
1008}
1009
1010
1011//------------------------------------------------------------------------
1012std::vector<G4VPhysicalVolume*>
1013G4tgbGeometryDumper::GetPVChildren( G4LogicalVolume* lv )
1014{
1015 G4PhysicalVolumeStore* pvstore = G4PhysicalVolumeStore::GetInstance();
1016 G4PhysicalVolumeStore::const_iterator ite;
1017 std::vector<G4VPhysicalVolume*> children;
1018 for( ite = pvstore->begin(); ite != pvstore->end(); ite++ )
1019 {
1020 if( (*ite)->GetMotherLogical() == lv )
1021 {
1022 children.push_back( *ite );
1023#ifdef G4VERBOSE
1024 if( G4tgrMessenger::GetVerboseLevel() >= 1 )
1025 {
1026 G4cout << " G4tgbGeometryDumper::GetPVChildren() - adding children: "
1027 << (*ite)->GetName() << " of " << lv->GetName() << G4endl;
1028 }
1029#endif
1030 }
1031 }
1032
1033 return children;
1034}
1035
1036
1037//------------------------------------------------------------------------
1038G4String G4tgbGeometryDumper::GetTGSolidType( const G4String& solidType )
1039{
1040 G4String newsolidType = solidType.substr(2,solidType.length() );
1041 for( size_t ii = 0; ii < newsolidType.length(); ii++ )
1042 {
1043 newsolidType[ii] = toupper(newsolidType[ii] );
1044 }
1045 return newsolidType;
1046}
1047
1048
1049//------------------------------------------------------------------------
1050G4double G4tgbGeometryDumper::MatDeterminant(G4RotationMatrix * ro)
1051{
1052 CLHEP::HepRep3x3 r = ro->rep3x3();
1053 return r.xx_*(r.yy_*r.zz_ - r.zy_*r.yz_)
1054 - r.yx_*(r.xy_*r.zz_ - r.zy_*r.xz_)
1055 + r.zx_*(r.xy_*r.yz_ - r.yy_*r.xz_);
1056}
1057
1058
1059//-----------------------------------------------------------------------
1060G4double G4tgbGeometryDumper::approxTo0( G4double val )
1061{
1062 G4double precision = G4GeometryTolerance::GetInstance()
1063 ->GetSurfaceTolerance();
1064
1065 if( std::fabs(val) < precision ) { val = 0; }
1066 return val;
1067}
1068
1069
1070//-----------------------------------------------------------------------
1071G4String G4tgbGeometryDumper::AddQuotes( const G4String& str )
1072{
1073 //--- look if there is a separating blank
1074
1075 G4bool bBlank = FALSE;
1076 size_t siz = str.length();
1077 for( size_t ii = 0; ii < siz; ii++ )
1078 {
1079 if( str.substr(ii,1) == " " )
1080 {
1081 bBlank = TRUE;
1082 break;
1083 }
1084 }
1085 G4String str2 = str;
1086 if( bBlank )
1087 {
1088 str2 = G4String("\"") + str2 + G4String("\"");
1089 }
1090 return str2;
1091}
1092
1093
1094//------------------------------------------------------------------------
1095G4String G4tgbGeometryDumper::SupressRefl( G4String name )
1096{
1097 G4int irefl = name.rfind("_refl");
1098 if( irefl != -1 )
1099 {
1100 name = name.substr( 0, irefl );
1101 }
1102 return name;
1103}
1104
1105//------------------------------------------------------------------------
1106G4String G4tgbGeometryDumper::SubstituteRefl( G4String name )
1107{
1108 G4int irefl = name.rfind("_refl");
1109 if( irefl != -1 )
1110 {
1111 name = name.substr( 0, irefl ) + "_REFL";
1112 }
1113 return name;
1114}
1115
1116
1117//------------------------------------------------------------------------
1118G4String G4tgbGeometryDumper::GetIsotopeName( G4Isotope* isot )
1119{
1120 G4String isotName = isot->GetName();
1121 // first look if this is isotope is already dumped,
1122 // with original isotope name or new one
1123 //
1124 std::map<G4String,G4Isotope*>::const_iterator ite;
1125 for( ite = theIsotopes.begin(); ite != theIsotopes.end(); ite++ )
1126 {
1127 if( isot == (*ite).second ) { return (*ite).first; }
1128 }
1129
1130 // Now look if there is another isotope dumped with same name,
1131 // and if found add _N to the name
1132 //
1133 ite = theIsotopes.find( isotName );
1134 if( ite != theIsotopes.end() ) // Isotope found with same name
1135 {
1136 G4Isotope* isotold = (*ite).second;
1137 if( isot != isotold ) // new isotope it is not the really
1138 { // the same one as isotope found
1139 if( !Same2G4Isotopes(isot, isotold))
1140 { // if the two have same data, use the old one
1141 G4int ii = 2; // G4Nist does names isotopes of same element
1142 // with same name
1143 for(;;ii++)
1144 {
1145 G4String newIsotName = isotName + "_"
1146 + G4UIcommand::ConvertToString(ii);
1147 std::map<G4String,G4Isotope*>::const_iterator ite2 =
1148 theIsotopes.find( newIsotName );
1149 if( ite2 == theIsotopes.end() )
1150 {
1151 isotName = newIsotName;
1152 break;
1153 }
1154 else
1155 {
1156 if( Same2G4Isotopes( isot, (*ite2).second ) )
1157 {
1158 isotName = newIsotName;
1159 break;
1160 }
1161 }
1162 }
1163 }
1164 }
1165 }
1166 return isotName;
1167}
1168
1169
1170//------------------------------------------------------------------------
1171template< class TYP > G4String G4tgbGeometryDumper::
1172GetObjectName( TYP* obj, std::map<G4String,TYP*> objectsDumped )
1173{
1174 G4String objName = obj->GetName();
1175
1176 // first look if this is objecy is already dumped,
1177 // with original object name or new one
1178 //
1179 typename std::map<G4String,TYP*>::const_iterator ite;
1180 for( ite = objectsDumped.begin(); ite != objectsDumped.end(); ite++ )
1181 {
1182 if( obj == (*ite).second ) { return (*ite).first; }
1183 }
1184
1185 // Now look if there is another object dumped with same name,
1186 // and if found add _N to the name
1187 //
1188 ite = objectsDumped.find( objName );
1189
1190 if( ite != objectsDumped.end() ) // Object found with same name
1191 {
1192 TYP* objold = (*ite).second;
1193 if( obj != objold ) // new object it is not the really
1194 { // the same one as object found
1195 G4int ii = 2;
1196 for(;;ii++)
1197 {
1198 G4String newObjName = objName + "_" + G4UIcommand::ConvertToString(ii);
1199 typename std::map<G4String,TYP*>::const_iterator ite2 =
1200 objectsDumped.find( newObjName );
1201 if( ite2 == objectsDumped.end() )
1202 {
1203 objName = newObjName;
1204 break;
1205 }
1206 }
1207 }
1208 }
1209 return objName;
1210}
1211
1212
1213//------------------------------------------------------------------------
1214G4bool G4tgbGeometryDumper::CheckIfLogVolExists( const G4String& name,
1215 G4LogicalVolume* pt )
1216{
1217 if( theLogVols.find( name ) != theLogVols.end() )
1218 {
1219 G4LogicalVolume* lvnew = (*(theLogVols.find(name))).second;
1220 if( lvnew != pt )
1221 {
1222 /*
1223 //---- Reflected volumes are repeated
1224
1225 G4ReflectionFactory* reffact = G4ReflectionFactory::Instance();
1226 if( !reffact->IsReflected( pt ) && !reffact->IsReflected( lvnew ) )
1227 {
1228 G4String ErrMessage = "LogVol found but not same as before: " + name;
1229 G4Exception("G4tgbGeometryDumper::CheckIfLogVolExists()",
1230 "InvalidSetup", FatalException, ErrMessage);
1231 }
1232 */
1233 }
1234 return 1;
1235 }
1236 else
1237 {
1238 return 0;
1239 }
1240}
1241
1242
1243//-----------------------------------------------------------------------
1244G4bool G4tgbGeometryDumper::CheckIfPhysVolExists( const G4String& name,
1245 G4VPhysicalVolume* pt )
1246{
1247#ifdef G4VERBOSE
1248 if( G4tgrMessenger::GetVerboseLevel() >= 1 )
1249 {
1250 G4cout << " G4tgbGeometryDumper::CheckIfPhysVolExists() - "
1251 << name << G4endl;
1252 }
1253#endif
1254 if( thePhysVols.find( name ) != thePhysVols.end() )
1255 {
1256 if( (*(thePhysVols.find(name))).second != pt )
1257 {
1258 // G4String ErrMessage = "Placement found but not same as before: "
1259 // + name;
1260 // G4Exception("G4tgbGeometryDumper::CheckIfPhysVolExists()",
1261 // "InvalidSetup", FatalException, ErrMessage);
1262 G4cerr << " G4tgbGeometryDumper::CheckIfPhysVolExists () -"
1263 << " Placement found but not same as before : " << name << G4endl;
1264 }
1265 return 1;
1266 }
1267 else
1268 {
1269 return 0;
1270 }
1271}
1272
1273
1274//-----------------------------------------------------------------------
1275G4String
1276G4tgbGeometryDumper::LookForExistingRotation( const G4RotationMatrix* rotm )
1277{
1278 G4String rmName = "";
1279
1280 std::map<G4String,G4RotationMatrix*>::const_iterator ite;
1281 for( ite = theRotMats.begin(); ite != theRotMats.end(); ite++ )
1282 {
1283 if( (*ite).second->isNear( *rotm ) )
1284 {
1285 rmName = (*ite).first;
1286 break;
1287 }
1288 }
1289 return rmName;
1290}
1291
1292
1293//------------------------------------------------------------------------
1294G4bool
1295G4tgbGeometryDumper::Same2G4Isotopes( G4Isotope* isot1, G4Isotope* isot2 )
1296{
1297 if ( (isot1->GetZ() != isot2->GetZ())
1298 || (isot1->GetN() != isot2->GetN())
1299 || (isot1->GetA() != isot2->GetA()) )
1300 {
1301 return 0;
1302 }
1303 else
1304 {
1305 return 1;
1306 }
1307}
1308
1309
1310//------------------------------------------------------------------------
1311const G4String& G4tgbGeometryDumper::FindSolidName( G4VSolid* solid )
1312{
1313 std::map<G4String,G4VSolid*>::const_iterator ite;
1314 for( ite = theSolids.begin(); ite != theSolids.end(); ite++ )
1315 {
1316 if( solid == (*ite).second ) { return (*ite).first; }
1317 }
1318
1319 if( ite == theSolids.end() )
1320 {
1321 G4Exception("G4tgbGeometryDumper::FindSolidName()", "ReadError",
1322 FatalException, "Programming error.");
1323 }
1324 return (*ite).first;
1325}
Note: See TracBrowser for help on using the repository browser.