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