Context Navigation

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

Last change on this file since 1230 was 1228, checked in by garnier, 16 years ago
update geant4.9.3 tag
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.12 2009/11/19 13:29:04 arce Exp $
28	// GEANT4 tag $Name: geant4-09-03 $
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	//------------------------------------------------------------------------
79	G4tgbGeometryDumper* G4tgbGeometryDumper::theInstance = 0;
80
81	//------------------------------------------------------------------------
82	G4tgbGeometryDumper::G4tgbGeometryDumper()
83	{
84
85	theRotationNumber = 0;
86	}
87
88	//------------------------------------------------------------------------
89	G4tgbGeometryDumper* G4tgbGeometryDumper::GetInstance()
90	{
91	if( theInstance == 0 ){
92	theInstance = new G4tgbGeometryDumper;
93	}
94
95	return theInstance;
96
97	}
98
99	//------------------------------------------------------------------------
100	void 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	//---------------------------------------------------------------------
109	G4VPhysicalVolume* 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	//------------------------------------------------------------------------
135	G4tgbGeometryDumper::~G4tgbGeometryDumper()
136	{
137	}
138
139	//------------------------------------------------------------------------
140	void 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	//------------------------------------------------------------------------
205	void
206	G4tgbGeometryDumper::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	//------------------------------------------------------------------------
268	void 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	//------------------------------------------------------------------------
394	void 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	//------------------------------------------------------------------------
460	G4String
461	G4tgbGeometryDumper::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	//------------------------------------------------------------------------
501	G4String 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	//------------------------------------------------------------------------
553	void 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	//------------------------------------------------------------------------
601	void 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	//------------------------------------------------------------------------
618	G4String 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( "SUBTRACTION", solid );
645
646	} else if (solidType == "INTERSECTIONSOLID") {
647	DumpBooleanVolume( "INTERSECTION", 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	//------------------------------------------------------------------------
667	void 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	//------------------------------------------------------------------------
717	void 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	//------------------------------------------------------------------------
729	std::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 > 180deg ) { angphi -= 360deg; }
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 > 180deg ) angphi -= 360deg;
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	//------------------------------------------------------------------------
945	G4String 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	//------------------------------------------------------------------------
994	std::vector<G4VPhysicalVolume*>
995	G4tgbGeometryDumper::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	//------------------------------------------------------------------------
1020	G4String 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	//------------------------------------------------------------------------
1032	G4double 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	//-----------------------------------------------------------------------
1042	G4double 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	//-----------------------------------------------------------------------
1053	G4String 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	//------------------------------------------------------------------------
1077	G4String 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	//------------------------------------------------------------------------
1088	G4String 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	//------------------------------------------------------------------------
1100	G4String 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	//------------------------------------------------------------------------
1153	template< class TYP > G4String G4tgbGeometryDumper::
1154	GetObjectName( 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	//------------------------------------------------------------------------
1196	G4bool 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	//-----------------------------------------------------------------------
1226	G4bool 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	//-----------------------------------------------------------------------
1257	G4String
1258	G4tgbGeometryDumper::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	//------------------------------------------------------------------------
1276	G4bool
1277	G4tgbGeometryDumper::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	//------------------------------------------------------------------------
1293	const 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.

Download in other formats: