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

Last change on this file since 1199 was 1035, checked in by garnier, 17 years ago

dossiers oublies

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