| Geant4 Cross Reference |
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 // G4GDMLWriteStructure implementation
27 //
28 // Author: Zoltan Torzsok, November 2007
29 // --------------------------------------------------------------------
30
31 #include "G4GDMLWriteStructure.hh"
32 #include "G4GDMLEvaluator.hh"
33
34 #include "G4Material.hh"
35 #include "G4ReflectedSolid.hh"
36 #include "G4DisplacedSolid.hh"
37 #include "G4LogicalVolumeStore.hh"
38 #include "G4PhysicalVolumeStore.hh"
39 #include "G4ReflectionFactory.hh"
40 #include "G4PVDivision.hh"
41 #include "G4PVReplica.hh"
42 #include "G4Region.hh"
43 #include "G4OpticalSurface.hh"
44 #include "G4LogicalSkinSurface.hh"
45 #include "G4LogicalBorderSurface.hh"
46
47 #include "G4ProductionCuts.hh"
48 #include "G4ProductionCutsTable.hh"
49 #include "G4Gamma.hh"
50 #include "G4Electron.hh"
51 #include "G4Positron.hh"
52 #include "G4Proton.hh"
53
54 #include "G4VSensitiveDetector.hh"
55 #include "G4AssemblyStore.hh"
56 #include "G4AssemblyVolume.hh"
57
58 G4int G4GDMLWriteStructure::levelNo = 0; // Counter for level being exported
59
60 // --------------------------------------------------------------------
61 G4GDMLWriteStructure::G4GDMLWriteStructure()
62 : G4GDMLWriteParamvol()
63 , maxLevel(INT_MAX)
64 {
65 reflFactory = G4ReflectionFactory::Instance();
66 }
67
68 // --------------------------------------------------------------------
69 G4GDMLWriteStructure::~G4GDMLWriteStructure()
70 {
71 }
72
73 // --------------------------------------------------------------------
74 void G4GDMLWriteStructure::DivisionvolWrite(
75 xercesc::DOMElement* volumeElement, const G4PVDivision* const divisionvol)
76 {
77 EAxis axis = kUndefined;
78 G4int number = 0;
79 G4double width = 0.0;
80 G4double offset = 0.0;
81 G4bool consuming = false;
82
83 divisionvol->GetReplicationData(axis, number, width, offset, consuming);
84 axis = divisionvol->GetDivisionAxis();
85
86 G4String unitString("mm");
87 G4String axisString("kUndefined");
88 if(axis == kXAxis)
89 {
90 axisString = "kXAxis";
91 }
92 else if(axis == kYAxis)
93 {
94 axisString = "kYAxis";
95 }
96 else if(axis == kZAxis)
97 {
98 axisString = "kZAxis";
99 }
100 else if(axis == kRho)
101 {
102 axisString = "kRho";
103 }
104 else if(axis == kPhi)
105 {
106 axisString = "kPhi";
107 unitString = "rad";
108 }
109
110 const G4String name = GenerateName(divisionvol->GetName(), divisionvol);
111 const G4String volumeref =
112 GenerateName(divisionvol->GetLogicalVolume()->GetName(),
113 divisionvol->GetLogicalVolume());
114
115 xercesc::DOMElement* divisionvolElement = NewElement("divisionvol");
116 divisionvolElement->setAttributeNode(NewAttribute("axis", axisString));
117 divisionvolElement->setAttributeNode(NewAttribute("number", number));
118 divisionvolElement->setAttributeNode(NewAttribute("width", width));
119 divisionvolElement->setAttributeNode(NewAttribute("offset", offset));
120 divisionvolElement->setAttributeNode(NewAttribute("unit", unitString));
121 xercesc::DOMElement* volumerefElement = NewElement("volumeref");
122 volumerefElement->setAttributeNode(NewAttribute("ref", volumeref));
123 divisionvolElement->appendChild(volumerefElement);
124 volumeElement->appendChild(divisionvolElement);
125 }
126
127 // --------------------------------------------------------------------
128 void G4GDMLWriteStructure::PhysvolWrite(xercesc::DOMElement* volumeElement,
129 const G4VPhysicalVolume* const physvol,
130 const G4Transform3D& T,
131 const G4String& ModuleName)
132 {
133 HepGeom::Scale3D scale;
134 HepGeom::Rotate3D rotate;
135 HepGeom::Translate3D translate;
136
137 T.getDecomposition(scale, rotate, translate);
138
139 const G4ThreeVector scl(scale(0, 0), scale(1, 1), scale(2, 2));
140 const G4ThreeVector rot = GetAngles(rotate.getRotation());
141 const G4ThreeVector pos = T.getTranslation();
142
143 const G4String name = GenerateName(physvol->GetName(), physvol);
144 const G4int copynumber = physvol->GetCopyNo();
145
146 xercesc::DOMElement* physvolElement = NewElement("physvol");
147 physvolElement->setAttributeNode(NewAttribute("name", name));
148 if(copynumber)
149 {
150 physvolElement->setAttributeNode(NewAttribute("copynumber", copynumber));
151 }
152
153 volumeElement->appendChild(physvolElement);
154
155 G4LogicalVolume* lv;
156 // Is it reflected?
157 if(reflFactory->IsReflected(physvol->GetLogicalVolume()))
158 {
159 lv = reflFactory->GetConstituentLV(physvol->GetLogicalVolume());
160 }
161 else
162 {
163 lv = physvol->GetLogicalVolume();
164 }
165
166 const G4String volumeref = GenerateName(lv->GetName(), lv);
167
168 if(ModuleName.empty())
169 {
170 xercesc::DOMElement* volumerefElement = NewElement("volumeref");
171 volumerefElement->setAttributeNode(NewAttribute("ref", volumeref));
172 physvolElement->appendChild(volumerefElement);
173 }
174 else
175 {
176 xercesc::DOMElement* fileElement = NewElement("file");
177 fileElement->setAttributeNode(NewAttribute("name", ModuleName));
178 fileElement->setAttributeNode(NewAttribute("volname", volumeref));
179 physvolElement->appendChild(fileElement);
180 }
181
182 if(std::fabs(pos.x()) > kLinearPrecision ||
183 std::fabs(pos.y()) > kLinearPrecision ||
184 std::fabs(pos.z()) > kLinearPrecision)
185 {
186 PositionWrite(physvolElement, name + "_pos", pos);
187 }
188 if(std::fabs(rot.x()) > kAngularPrecision ||
189 std::fabs(rot.y()) > kAngularPrecision ||
190 std::fabs(rot.z()) > kAngularPrecision)
191 {
192 RotationWrite(physvolElement, name + "_rot", rot);
193 }
194 if(std::fabs(scl.x() - 1.0) > kRelativePrecision ||
195 std::fabs(scl.y() - 1.0) > kRelativePrecision ||
196 std::fabs(scl.z() - 1.0) > kRelativePrecision)
197 {
198 ScaleWrite(physvolElement, name + "_scl", scl);
199 }
200 }
201
202 // --------------------------------------------------------------------
203 void G4GDMLWriteStructure::ReplicavolWrite(
204 xercesc::DOMElement* volumeElement, const G4VPhysicalVolume* const replicavol)
205 {
206 EAxis axis = kUndefined;
207 G4int number = 0;
208 G4double width = 0.0;
209 G4double offset = 0.0;
210 G4bool consuming = false;
211 G4String unitString("mm");
212
213 replicavol->GetReplicationData(axis, number, width, offset, consuming);
214
215 const G4String volumeref = GenerateName(
216 replicavol->GetLogicalVolume()->GetName(), replicavol->GetLogicalVolume());
217
218 xercesc::DOMElement* replicavolElement = NewElement("replicavol");
219 replicavolElement->setAttributeNode(NewAttribute("number", number));
220 xercesc::DOMElement* volumerefElement = NewElement("volumeref");
221 volumerefElement->setAttributeNode(NewAttribute("ref", volumeref));
222 replicavolElement->appendChild(volumerefElement);
223 xercesc::DOMElement* replicateElement = NewElement("replicate_along_axis");
224 replicavolElement->appendChild(replicateElement);
225
226 xercesc::DOMElement* dirElement = NewElement("direction");
227 if(axis == kXAxis)
228 {
229 dirElement->setAttributeNode(NewAttribute("x", "1"));
230 }
231 else if(axis == kYAxis)
232 {
233 dirElement->setAttributeNode(NewAttribute("y", "1"));
234 }
235 else if(axis == kZAxis)
236 {
237 dirElement->setAttributeNode(NewAttribute("z", "1"));
238 }
239 else if(axis == kRho)
240 {
241 dirElement->setAttributeNode(NewAttribute("rho", "1"));
242 }
243 else if(axis == kPhi)
244 {
245 dirElement->setAttributeNode(NewAttribute("phi", "1"));
246 unitString = "rad";
247 }
248 replicateElement->appendChild(dirElement);
249
250 xercesc::DOMElement* widthElement = NewElement("width");
251 widthElement->setAttributeNode(NewAttribute("value", width));
252 widthElement->setAttributeNode(NewAttribute("unit", unitString));
253 replicateElement->appendChild(widthElement);
254
255 xercesc::DOMElement* offsetElement = NewElement("offset");
256 offsetElement->setAttributeNode(NewAttribute("value", offset));
257 offsetElement->setAttributeNode(NewAttribute("unit", unitString));
258 replicateElement->appendChild(offsetElement);
259
260 volumeElement->appendChild(replicavolElement);
261 }
262
263 // --------------------------------------------------------------------
264 void G4GDMLWriteStructure::AssemblyWrite(xercesc::DOMElement* volumeElement,
265 const G4int assemblyID)
266 {
267 G4AssemblyStore* assemblies = G4AssemblyStore::GetInstance();
268 G4AssemblyVolume* myassembly = assemblies->GetAssembly(assemblyID);
269
270 xercesc::DOMElement* assemblyElement = NewElement("assembly");
271 G4String name = "Assembly_" + std::to_string(assemblyID);
272
273 assemblyElement->setAttributeNode(NewAttribute("name", name));
274
275 auto vit = myassembly->GetTripletsIterator();
276
277 G4int depth = 0;
278 const G4String ModuleName;
279
280 for(std::size_t i5 = 0; i5 < myassembly->TotalTriplets(); ++i5)
281 {
282 G4LogicalVolume* lvol = (*vit).GetVolume();
283 if (lvol == nullptr)
284 {
285 G4String message = "Nested assemblies not yet supported for exporting. Sorry!";
286 G4Exception("G4GDMLWriteStructure::AssemblyWrite()", "InvalidSetup",
287 FatalException, message);
288 return;
289 }
290 TraverseVolumeTree(lvol, depth + 1);
291
292 const G4ThreeVector rot = GetAngles((*vit).GetRotation()->inverse());
293 const G4ThreeVector pos = (*vit).GetTranslation();
294
295 const G4String pname =
296 GenerateName((*vit).GetVolume()->GetName() + "_pv", &(*vit));
297
298 xercesc::DOMElement* physvolElement = NewElement("physvol");
299 physvolElement->setAttributeNode(NewAttribute("name", pname));
300
301 assemblyElement->appendChild(physvolElement);
302
303 const G4String volumeref =
304 GenerateName((*vit).GetVolume()->GetName(), (*vit).GetVolume());
305
306 xercesc::DOMElement* volumerefElement = NewElement("volumeref");
307 volumerefElement->setAttributeNode(NewAttribute("ref", volumeref));
308 physvolElement->appendChild(volumerefElement);
309
310 if(std::fabs(pos.x()) > kLinearPrecision ||
311 std::fabs(pos.y()) > kLinearPrecision ||
312 std::fabs(pos.z()) > kLinearPrecision)
313 {
314 PositionWrite(physvolElement,name+"_position_" + std::to_string(i5), pos);
315 }
316
317 if(std::fabs(rot.x()) > kAngularPrecision ||
318 std::fabs(rot.y()) > kAngularPrecision ||
319 std::fabs(rot.z()) > kAngularPrecision)
320 {
321 RotationWrite(physvolElement,name+"_rotation_" + std::to_string(i5), rot);
322 }
323 ++vit;
324 }
325
326 volumeElement->appendChild(assemblyElement);
327 }
328
329 // --------------------------------------------------------------------
330 void G4GDMLWriteStructure::BorderSurfaceCache(
331 const G4LogicalBorderSurface* const bsurf)
332 {
333 if(bsurf == nullptr)
334 {
335 return;
336 }
337
338 const G4SurfaceProperty* psurf = bsurf->GetSurfaceProperty();
339
340 // Generate the new element for border-surface
341 //
342 const G4String& bsname = GenerateName(bsurf->GetName(), bsurf);
343 const G4String& psname = GenerateName(psurf->GetName(), psurf);
344 xercesc::DOMElement* borderElement = NewElement("bordersurface");
345 borderElement->setAttributeNode(NewAttribute("name", bsname));
346 borderElement->setAttributeNode(NewAttribute("surfaceproperty", psname));
347
348 const G4String volumeref1 =
349 GenerateName(bsurf->GetVolume1()->GetName(), bsurf->GetVolume1());
350 const G4String volumeref2 =
351 GenerateName(bsurf->GetVolume2()->GetName(), bsurf->GetVolume2());
352 xercesc::DOMElement* volumerefElement1 = NewElement("physvolref");
353 xercesc::DOMElement* volumerefElement2 = NewElement("physvolref");
354 volumerefElement1->setAttributeNode(NewAttribute("ref", volumeref1));
355 volumerefElement2->setAttributeNode(NewAttribute("ref", volumeref2));
356 borderElement->appendChild(volumerefElement1);
357 borderElement->appendChild(volumerefElement2);
358
359 if(FindOpticalSurface(psurf))
360 {
361 const G4OpticalSurface* opsurf =
362 dynamic_cast<const G4OpticalSurface*>(psurf);
363 if(opsurf == nullptr)
364 {
365 G4Exception("G4GDMLWriteStructure::BorderSurfaceCache()", "InvalidSetup",
366 FatalException, "No optical surface found!");
367 return;
368 }
369 OpticalSurfaceWrite(solidsElement, opsurf);
370 }
371
372 borderElementVec.push_back(borderElement);
373 }
374
375 // --------------------------------------------------------------------
376 void G4GDMLWriteStructure::SkinSurfaceCache(
377 const G4LogicalSkinSurface* const ssurf)
378 {
379 if(ssurf == nullptr)
380 {
381 return;
382 }
383
384 const G4SurfaceProperty* psurf = ssurf->GetSurfaceProperty();
385
386 // Generate the new element for border-surface
387 //
388 const G4String& ssname = GenerateName(ssurf->GetName(), ssurf);
389 const G4String& psname = GenerateName(psurf->GetName(), psurf);
390 xercesc::DOMElement* skinElement = NewElement("skinsurface");
391 skinElement->setAttributeNode(NewAttribute("name", ssname));
392 skinElement->setAttributeNode(NewAttribute("surfaceproperty", psname));
393
394 const G4String volumeref = GenerateName(ssurf->GetLogicalVolume()->GetName(),
395 ssurf->GetLogicalVolume());
396 xercesc::DOMElement* volumerefElement = NewElement("volumeref");
397 volumerefElement->setAttributeNode(NewAttribute("ref", volumeref));
398 skinElement->appendChild(volumerefElement);
399
400 if(FindOpticalSurface(psurf))
401 {
402 const G4OpticalSurface* opsurf =
403 dynamic_cast<const G4OpticalSurface*>(psurf);
404 if(opsurf == nullptr)
405 {
406 G4Exception("G4GDMLWriteStructure::SkinSurfaceCache()", "InvalidSetup",
407 FatalException, "No optical surface found!");
408 return;
409 }
410 OpticalSurfaceWrite(solidsElement, opsurf);
411 }
412
413 skinElementVec.push_back(skinElement);
414 }
415
416 // --------------------------------------------------------------------
417 G4bool G4GDMLWriteStructure::FindOpticalSurface(const G4SurfaceProperty* psurf)
418 {
419 const G4OpticalSurface* osurf = dynamic_cast<const G4OpticalSurface*>(psurf);
420 auto pos = std::find(opt_vec.cbegin(), opt_vec.cend(), osurf);
421 if(pos != opt_vec.cend())
422 {
423 return false;
424 } // item already created!
425
426 opt_vec.push_back(osurf); // cache it for future reference
427 return true;
428 }
429
430 // --------------------------------------------------------------------
431 const G4LogicalSkinSurface* G4GDMLWriteStructure::GetSkinSurface(
432 const G4LogicalVolume* const lvol)
433 {
434 G4LogicalSkinSurface* surf = nullptr;
435 std::size_t nsurf = G4LogicalSkinSurface::GetNumberOfSkinSurfaces();
436 if(nsurf)
437 {
438 const G4LogicalSkinSurfaceTable* stable =
439 G4LogicalSkinSurface::GetSurfaceTable();
440 auto pos = stable->find(lvol);
441 if(pos != stable->cend())
442 {
443 surf = pos->second;
444 }
445 }
446 return surf;
447 }
448
449 // --------------------------------------------------------------------
450 const G4LogicalBorderSurface* G4GDMLWriteStructure::GetBorderSurface(
451 const G4VPhysicalVolume* const pvol)
452 {
453 G4LogicalBorderSurface* surf = nullptr;
454 std::size_t nsurf = G4LogicalBorderSurface::GetNumberOfBorderSurfaces();
455 if(nsurf)
456 {
457 const G4LogicalBorderSurfaceTable* btable =
458 G4LogicalBorderSurface::GetSurfaceTable();
459 for(auto pos = btable->cbegin(); pos != btable->cend(); ++pos)
460 {
461 if(pvol == pos->first.first) // just the first in the couple
462 { // could be enough?
463 surf = pos->second; // break;
464 BorderSurfaceCache(surf);
465 }
466 }
467 }
468 return surf;
469 }
470
471 // --------------------------------------------------------------------
472 void G4GDMLWriteStructure::SurfacesWrite()
473 {
474 #ifdef G4VERBOSE
475 G4cout << "G4GDML: Writing surfaces..." << G4endl;
476 #endif
477 for(auto pos = skinElementVec.cbegin();
478 pos != skinElementVec.cend(); ++pos)
479 {
480 structureElement->appendChild(*pos);
481 }
482 for(auto pos = borderElementVec.cbegin();
483 pos != borderElementVec.cend(); ++pos)
484 {
485 structureElement->appendChild(*pos);
486 }
487 }
488
489 // --------------------------------------------------------------------
490 void G4GDMLWriteStructure::StructureWrite(xercesc::DOMElement* gdmlElement)
491 {
492 #ifdef G4VERBOSE
493 G4cout << "G4GDML: Writing structure..." << G4endl;
494 #endif
495
496 // filling the list of phys volumes that are parts of assemblies
497
498 G4AssemblyStore* assemblies = G4AssemblyStore::GetInstance();
499
500 for(auto it = assemblies->cbegin(); it != assemblies->cend(); ++it)
501 {
502 auto vit = (*it)->GetVolumesIterator();
503
504 for(std::size_t i5 = 0; i5 < (*it)->TotalImprintedVolumes(); ++i5)
505 {
506 G4String pvname = (*vit)->GetName();
507 std::size_t pos = pvname.find("_impr_") + 6;
508 G4String impID = pvname.substr(pos);
509
510 pos = impID.find("_");
511 impID = impID.substr(0, pos);
512
513 assemblyVolMap[*vit] = (*it)->GetAssemblyID();
514 imprintsMap[*vit] = std::atoi(impID.c_str());
515 ++vit;
516 }
517 }
518
519 structureElement = NewElement("structure");
520 gdmlElement->appendChild(structureElement);
521 }
522
523 // --------------------------------------------------------------------
524 G4Transform3D G4GDMLWriteStructure::TraverseVolumeTree(
525 const G4LogicalVolume* const volumePtr, const G4int depth)
526 {
527 if(VolumeMap().find(volumePtr) != VolumeMap().cend())
528 {
529 return VolumeMap()[volumePtr]; // Volume is already processed
530 }
531
532 G4VSolid* solidPtr = volumePtr->GetSolid();
533 G4Transform3D R, invR;
534 G4int trans = 0;
535
536 std::map<const G4LogicalVolume*, G4GDMLAuxListType>::iterator auxiter;
537
538 ++levelNo;
539
540 while(true) // Solve possible displacement/reflection
541 { // of the referenced solid!
542 if(trans > maxTransforms)
543 {
544 G4String ErrorMessage = "Referenced solid in volume '" +
545 volumePtr->GetName() +
546 "' was displaced/reflected too many times!";
547 G4Exception("G4GDMLWriteStructure::TraverseVolumeTree()", "InvalidSetup",
548 FatalException, ErrorMessage);
549 }
550
551 if(G4ReflectedSolid* refl = dynamic_cast<G4ReflectedSolid*>(solidPtr))
552 {
553 R = R * refl->GetTransform3D();
554 solidPtr = refl->GetConstituentMovedSolid();
555 ++trans;
556 continue;
557 }
558
559 if(G4DisplacedSolid* disp = dynamic_cast<G4DisplacedSolid*>(solidPtr))
560 {
561 R = R * G4Transform3D(disp->GetObjectRotation(),
562 disp->GetObjectTranslation());
563 solidPtr = disp->GetConstituentMovedSolid();
564 ++trans;
565 continue;
566 }
567
568 break;
569 }
570
571 // check if it is a reflected volume
572 G4LogicalVolume* tmplv = const_cast<G4LogicalVolume*>(volumePtr);
573
574 if(reflFactory->IsReflected(tmplv))
575 {
576 tmplv = reflFactory->GetConstituentLV(tmplv);
577 if(VolumeMap().find(tmplv) != VolumeMap().cend())
578 {
579 return R; // Volume is already processed
580 }
581 }
582
583 // Only compute the inverse when necessary!
584 //
585 if(trans > 0)
586 {
587 invR = R.inverse();
588 }
589
590 const G4String name = GenerateName(tmplv->GetName(), tmplv);
591
592 G4String materialref = "NULL";
593
594 if(volumePtr->GetMaterial())
595 {
596 materialref = GenerateName(volumePtr->GetMaterial()->GetName(),
597 volumePtr->GetMaterial());
598 }
599
600 const G4String solidref = GenerateName(solidPtr->GetName(), solidPtr);
601
602 xercesc::DOMElement* volumeElement = NewElement("volume");
603 volumeElement->setAttributeNode(NewAttribute("name", name));
604 xercesc::DOMElement* materialrefElement = NewElement("materialref");
605 materialrefElement->setAttributeNode(NewAttribute("ref", materialref));
606 volumeElement->appendChild(materialrefElement);
607 xercesc::DOMElement* solidrefElement = NewElement("solidref");
608 solidrefElement->setAttributeNode(NewAttribute("ref", solidref));
609 volumeElement->appendChild(solidrefElement);
610
611 std::size_t daughterCount = volumePtr->GetNoDaughters();
612
613 if(levelNo == maxLevel) // Stop exporting if reached levels limit
614 {
615 daughterCount = 0;
616 }
617
618 std::map<G4int, std::vector<G4int> > assemblyIDToAddedImprints;
619
620 for(std::size_t i = 0; i < daughterCount; ++i) // Traverse all the children!
621 {
622 const G4VPhysicalVolume* const physvol = volumePtr->GetDaughter(i);
623 const G4String ModuleName = Modularize(physvol, depth);
624
625 G4Transform3D daughterR;
626
627 if(ModuleName.empty()) // Check if subtree requested to be
628 { // a separate module!
629 daughterR = TraverseVolumeTree(physvol->GetLogicalVolume(), depth + 1);
630 }
631 else
632 {
633 G4GDMLWriteStructure writer;
634 daughterR = writer.Write(ModuleName, physvol->GetLogicalVolume(),
635 SchemaLocation, depth + 1);
636 }
637
638 if(const G4PVDivision* const divisionvol =
639 dynamic_cast<const G4PVDivision*>(physvol)) // Is it division?
640 {
641 if(!G4Transform3D::Identity.isNear(invR * daughterR, kRelativePrecision))
642 {
643 G4String ErrorMessage = "Division volume in '" + name +
644 "' can not be related to reflected solid!";
645 G4Exception("G4GDMLWriteStructure::TraverseVolumeTree()",
646 "InvalidSetup", FatalException, ErrorMessage);
647 }
648 DivisionvolWrite(volumeElement, divisionvol);
649 }
650 else
651 {
652 if(physvol->IsParameterised()) // Is it a paramvol?
653 {
654 if(!G4Transform3D::Identity.isNear(invR * daughterR,
655 kRelativePrecision))
656 {
657 G4String ErrorMessage = "Parameterised volume in '" + name +
658 "' can not be related to reflected solid!";
659 G4Exception("G4GDMLWriteStructure::TraverseVolumeTree()",
660 "InvalidSetup", FatalException, ErrorMessage);
661 }
662 ParamvolWrite(volumeElement, physvol);
663 }
664 else
665 {
666 if(physvol->IsReplicated()) // Is it a replicavol?
667 {
668 if(!G4Transform3D::Identity.isNear(invR * daughterR,
669 kRelativePrecision))
670 {
671 G4String ErrorMessage = "Replica volume in '" + name +
672 "' can not be related to reflected solid!";
673 G4Exception("G4GDMLWriteStructure::TraverseVolumeTree()",
674 "InvalidSetup", FatalException, ErrorMessage);
675 }
676 ReplicavolWrite(volumeElement, physvol);
677 }
678 else // Is it a physvol or an assembly?
679 {
680 if(assemblyVolMap.find(physvol) != assemblyVolMap.cend())
681 {
682 G4int assemblyID = assemblyVolMap[physvol];
683
684 G4String assemblyref = "Assembly_" + std::to_string(assemblyID);
685
686 // here I need to retrieve the imprint ID
687
688 G4int imprintID = imprintsMap[physvol];
689
690 // there are 2 steps:
691 //
692 // 1) add assembly to the structure if that has not yet been done
693 // (but after the constituents volumes have been added)
694 //
695
696 if(std::find(addedAssemblies.cbegin(), addedAssemblies.cend(),
697 assemblyID) == addedAssemblies.cend())
698 {
699 AssemblyWrite(structureElement, assemblyID);
700 addedAssemblies.push_back(assemblyID);
701 assemblyIDToAddedImprints[assemblyID] = std::vector<G4int>();
702 }
703
704 // 2) add the assembly (as physical volume) to the mother volume
705 // (but only once), using it's original position and rotation.
706 //
707
708 // here I need a check if assembly has been already added to the
709 // mother volume
710 std::vector<G4int>& addedImprints = assemblyIDToAddedImprints[assemblyID];
711 if(std::find(addedImprints.cbegin(), addedImprints.cend(),
712 imprintID) == addedImprints.cend())
713 {
714 G4String imprintname = "Imprint_" + std::to_string(imprintID) + "_";
715 imprintname = GenerateName(imprintname, physvol);
716
717 // I need to get those two from the container of imprints from
718 // the assembly I have the imprint ID, I need to get pos and rot
719 //
720 G4Transform3D& transf = G4AssemblyStore::GetInstance()
721 ->GetAssembly(assemblyID)
722 ->GetImprintTransformation(imprintID);
723
724 HepGeom::Scale3D scale;
725 HepGeom::Rotate3D rotate;
726 HepGeom::Translate3D translate;
727
728 transf.getDecomposition(scale, rotate, translate);
729
730 const G4ThreeVector scl(scale(0, 0), scale(1, 1), scale(2, 2));
731 const G4ThreeVector rot = GetAngles(rotate.getRotation().inverse());
732 const G4ThreeVector pos = transf.getTranslation();
733
734 // here I need a normal physvol referencing to my assemblyref
735
736 xercesc::DOMElement* physvolElement = NewElement("physvol");
737 physvolElement->setAttributeNode(NewAttribute("name", imprintname));
738
739 xercesc::DOMElement* volumerefElement = NewElement("volumeref");
740 volumerefElement->setAttributeNode(NewAttribute("ref", assemblyref));
741 physvolElement->appendChild(volumerefElement);
742
743 if(std::fabs(pos.x()) > kLinearPrecision ||
744 std::fabs(pos.y()) > kLinearPrecision ||
745 std::fabs(pos.z()) > kLinearPrecision)
746 {
747 PositionWrite(physvolElement, imprintname + "_pos", pos);
748 }
749 if(std::fabs(rot.x()) > kAngularPrecision ||
750 std::fabs(rot.y()) > kAngularPrecision ||
751 std::fabs(rot.z()) > kAngularPrecision)
752 {
753 RotationWrite(physvolElement, imprintname + "_rot", rot);
754 }
755 if(std::fabs(scl.x() - 1.0) > kRelativePrecision ||
756 std::fabs(scl.y() - 1.0) > kRelativePrecision ||
757 std::fabs(scl.z() - 1.0) > kRelativePrecision)
758 {
759 ScaleWrite(physvolElement, name + "_scl", scl);
760 }
761
762 volumeElement->appendChild(physvolElement);
763 //
764 addedImprints.push_back(imprintID);
765 }
766 }
767 else // not part of assembly, so a normal physical volume
768 {
769 G4RotationMatrix rot;
770 if(physvol->GetFrameRotation() != nullptr)
771 {
772 rot = *(physvol->GetFrameRotation());
773 }
774 G4Transform3D P(rot, physvol->GetObjectTranslation());
775
776 PhysvolWrite(volumeElement, physvol, invR * P * daughterR,
777 ModuleName);
778 }
779 }
780 }
781 }
782 // BorderSurfaceCache(GetBorderSurface(physvol));
783 GetBorderSurface(physvol);
784 }
785
786 if(cexport)
787 {
788 ExportEnergyCuts(volumePtr);
789 }
790 // Add optional energy cuts
791
792 if(sdexport)
793 {
794 ExportSD(volumePtr);
795 }
796 // Add optional SDs
797
798 // Here write the auxiliary info
799 //
800 auxiter = auxmap.find(volumePtr);
801 if(auxiter != auxmap.cend())
802 {
803 AddAuxInfo(&(auxiter->second), volumeElement);
804 }
805
806 structureElement->appendChild(volumeElement);
807 // Append the volume AFTER traversing the children so that
808 // the order of volumes will be correct!
809
810 VolumeMap()[tmplv] = R;
811
812 AddExtension(volumeElement, volumePtr);
813 // Add any possible user defined extension attached to a volume
814
815 AddMaterial(volumePtr->GetMaterial());
816 // Add the involved materials and solids!
817
818 AddSolid(solidPtr);
819
820 SkinSurfaceCache(GetSkinSurface(volumePtr));
821
822 return R;
823 }
824
825 // --------------------------------------------------------------------
826 void G4GDMLWriteStructure::AddVolumeAuxiliary(G4GDMLAuxStructType myaux,
827 const G4LogicalVolume* const lvol)
828 {
829 auto pos = auxmap.find(lvol);
830
831 if(pos == auxmap.cend())
832 {
833 auxmap[lvol] = G4GDMLAuxListType();
834 }
835
836 auxmap[lvol].push_back(myaux);
837 }
838
839 // --------------------------------------------------------------------
840 void G4GDMLWriteStructure::SetEnergyCutsExport(G4bool fcuts)
841 {
842 cexport = fcuts;
843 }
844
845 // --------------------------------------------------------------------
846 void G4GDMLWriteStructure::ExportEnergyCuts(const G4LogicalVolume* const lvol)
847 {
848 G4GDMLEvaluator eval;
849 G4ProductionCuts* pcuts = lvol->GetRegion()->GetProductionCuts();
850 G4ProductionCutsTable* ctab = G4ProductionCutsTable::GetProductionCutsTable();
851 G4Gamma* gamma = G4Gamma::Gamma();
852 G4Electron* eminus = G4Electron::Electron();
853 G4Positron* eplus = G4Positron::Positron();
854 G4Proton* proton = G4Proton::Proton();
855
856 G4double gamma_cut = ctab->ConvertRangeToEnergy(
857 gamma, lvol->GetMaterial(), pcuts->GetProductionCut("gamma"));
858 G4double eminus_cut = ctab->ConvertRangeToEnergy(
859 eminus, lvol->GetMaterial(), pcuts->GetProductionCut("e-"));
860 G4double eplus_cut = ctab->ConvertRangeToEnergy(
861 eplus, lvol->GetMaterial(), pcuts->GetProductionCut("e+"));
862 G4double proton_cut = ctab->ConvertRangeToEnergy(
863 proton, lvol->GetMaterial(), pcuts->GetProductionCut("proton"));
864
865 G4GDMLAuxStructType gammainfo = { "gammaECut",
866 eval.ConvertToString(gamma_cut), "MeV", 0 };
867 G4GDMLAuxStructType eminusinfo = { "electronECut",
868 eval.ConvertToString(eminus_cut), "MeV",
869 0 };
870 G4GDMLAuxStructType eplusinfo = { "positronECut",
871 eval.ConvertToString(eplus_cut), "MeV", 0 };
872 G4GDMLAuxStructType protinfo = { "protonECut",
873 eval.ConvertToString(proton_cut), "MeV", 0 };
874
875 AddVolumeAuxiliary(gammainfo, lvol);
876 AddVolumeAuxiliary(eminusinfo, lvol);
877 AddVolumeAuxiliary(eplusinfo, lvol);
878 AddVolumeAuxiliary(protinfo, lvol);
879 }
880
881 // --------------------------------------------------------------------
882 void G4GDMLWriteStructure::SetSDExport(G4bool fsd)
883 {
884 sdexport = fsd;
885 }
886
887 // --------------------------------------------------------------------
888 void G4GDMLWriteStructure::ExportSD(const G4LogicalVolume* const lvol)
889 {
890 G4VSensitiveDetector* sd = lvol->GetMasterSensitiveDetector();
891
892 if(sd != nullptr)
893 {
894 G4String SDname = sd->GetName();
895
896 G4GDMLAuxStructType SDinfo = { "SensDet", SDname, "", 0 };
897 AddVolumeAuxiliary(SDinfo, lvol);
898 }
899 }
900
901 // --------------------------------------------------------------------
902 G4int G4GDMLWriteStructure::GetMaxExportLevel() const
903 {
904 return maxLevel;
905 }
906
907 // --------------------------------------------------------------------
908 void G4GDMLWriteStructure::SetMaxExportLevel(G4int level)
909 {
910 if(level <= 0)
911 {
912 G4Exception("G4GDMLWriteStructure::TraverseVolumeTree()", "InvalidSetup",
913 FatalException, "Levels to export must be greater than zero!");
914 return;
915 }
916 maxLevel = level;
917 levelNo = 0;
918 }
919