diff --git a/geometry/detector/ThinQuartz/DetectorArray/DetectorArray.gdml b/geometry/detector/ThinQuartz/DetectorArray/DetectorArray.gdml index ad25af5d4..92bcad72e 100755 --- a/geometry/detector/ThinQuartz/DetectorArray/DetectorArray.gdml +++ b/geometry/detector/ThinQuartz/DetectorArray/DetectorArray.gdml @@ -1,1142 +1,1142 @@ - - - - - + + + + + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/include/remollDetectorConstruction.hh b/include/remollDetectorConstruction.hh index 047a2dcba..c39aea00e 100644 --- a/include/remollDetectorConstruction.hh +++ b/include/remollDetectorConstruction.hh @@ -7,6 +7,7 @@ #include "G4GenericMessenger.hh" #include "G4Types.hh" #include "G4Version.hh" +#include "remollGDMLReadStructure.hh" #include #include @@ -25,7 +26,7 @@ class remollDetectorConstruction : public G4VUserDetectorConstruction { public: - remollDetectorConstruction(const G4String& name, const G4String& gdmlfile); + remollDetectorConstruction(const G4String& name, const G4String& gdmlfile, remollGDMLReadStructure* rs=nullptr); virtual ~remollDetectorConstruction(); private: diff --git a/include/remollGDMLReadStructure.hh b/include/remollGDMLReadStructure.hh new file mode 100644 index 000000000..6d461bc59 --- /dev/null +++ b/include/remollGDMLReadStructure.hh @@ -0,0 +1,25 @@ +// -*- coding: utf-8 -*- +// vim: ai ts=2 sts=2 et sw=2 ft=cpp +// author : Prakash +// date : 2024-02-18 + +#pragma once + +#include +#include + +#include "G4GDMLReadStructure.hh" + +class remollGDMLReadStructure : public G4GDMLReadStructure { + public: + remollGDMLReadStructure(); + + void RotationRead(const xercesc::DOMElement* const vectorElement, G4RotationMatrix& rot); + G4LogicalVolume* FileRead(const xercesc::DOMElement* const fileElement); + void PhysvolRead( const xercesc::DOMElement* const physvolElement, G4AssemblyVolume* pAssembly=0); + void ParametersRead(const xercesc::DOMElement* const element); + void Volume_contentRead( const xercesc::DOMElement* const volumeElement) override; + + private: + std::map>> rotations; +}; diff --git a/remoll.cc b/remoll.cc index b4def9fe3..34f81c52f 100644 --- a/remoll.cc +++ b/remoll.cc @@ -28,6 +28,7 @@ typedef G4RunManager RunManager; #include "remollActionInitialization.hh" #include "remollDetectorConstruction.hh" #include "remollParallelConstruction.hh" +#include "remollGDMLReadStructure.hh" #include "remollSearchPath.hh" @@ -155,8 +156,9 @@ int main(int argc, char** argv) { remollIO::GetInstance(); // Detector geometry + remollGDMLReadStructure* rs = new remollGDMLReadStructure(); G4String material_name = "material"; - remollDetectorConstruction* detector = new remollDetectorConstruction(material_name, geometry_gdmlfile); + remollDetectorConstruction* detector = new remollDetectorConstruction(material_name, geometry_gdmlfile,rs); // Parallel world geometry G4String parallel_name = "parallel"; // Note: name must correspond with name of G4ParallelWorldPhysics remollParallelConstruction* parallel = new remollParallelConstruction(parallel_name, parallel_gdmlfile); diff --git a/src/remollDetectorConstruction.cc b/src/remollDetectorConstruction.cc index 09caaf0cb..5228c01d8 100644 --- a/src/remollDetectorConstruction.cc +++ b/src/remollDetectorConstruction.cc @@ -28,6 +28,7 @@ // GDML export #include "G4GDMLParser.hh" +#include "remollGDMLReadStructure.hh" //CADMesh #ifdef __USE_CADMESH @@ -53,8 +54,9 @@ G4ThreadLocal remollGlobalField* remollDetectorConstruction::fGlobalField = 0; G4UserLimits* remollDetectorConstruction::fKryptoniteUserLimits = new G4UserLimits(0,0,0,DBL_MAX,DBL_MAX); -remollDetectorConstruction::remollDetectorConstruction(const G4String& name, const G4String& gdmlfile) +remollDetectorConstruction::remollDetectorConstruction(const G4String& name, const G4String& gdmlfile, remollGDMLReadStructure* read) : fVerboseLevel(0), + fGDMLParser(read), fGDMLValidate(false), fGDMLOverlapCheck(true), fGDMLPath(""), diff --git a/src/remollGDMLReadStructure.cc b/src/remollGDMLReadStructure.cc new file mode 100644 index 000000000..3a6c6a180 --- /dev/null +++ b/src/remollGDMLReadStructure.cc @@ -0,0 +1,502 @@ +// -*- coding: utf-8 -*- +// vim: ai ts=2 sts=2 et sw=2 ft=cpp +// author : Prakash [pgautam@jlab.org] +// date : 2024-02-18 + +#include "G4UnitsTable.hh" +#include "G4LogicalVolume.hh" +#include "G4VPhysicalVolume.hh" +#include "G4AssemblyVolume.hh" +#include "G4ReflectionFactory.hh" + +#include "remollGDMLReadStructure.hh" + +remollGDMLReadStructure::remollGDMLReadStructure():G4GDMLReadStructure() +{ + auto rotatex=[](G4RotationMatrix& rot,G4ThreeVector angle) { rot.rotateX(angle.x()); }; + auto rotatey=[](G4RotationMatrix& rot,G4ThreeVector angle) { rot.rotateY(angle.y()); }; + auto rotatez=[](G4RotationMatrix& rot,G4ThreeVector angle) { rot.rotateZ(angle.z()); }; + rotations = { + {"xyz", {rotatex,rotatey,rotatez}}, + {"xzy", {rotatex,rotatez,rotatey}}, + {"yxz", {rotatey,rotatex,rotatez}}, + {"yzx", {rotatey,rotatez,rotatex}}, + {"zxy", {rotatez,rotatex,rotatey}}, + {"zyx", {rotatez,rotatey,rotatex}} + }; + +} + +void remollGDMLReadStructure::RotationRead(const xercesc::DOMElement* const vectorElement, G4RotationMatrix& rot) +{ + G4double unit = 1.0; + G4ThreeVector vec; + std::string order = "xyz"; + + const xercesc::DOMNamedNodeMap* const attributes = vectorElement->getAttributes(); + XMLSize_t attributeCount = attributes->getLength(); + + for(XMLSize_t attribute_index = 0; attribute_index < attributeCount; ++attribute_index) + { + xercesc::DOMNode* attribute_node = attributes->item(attribute_index); + + if(attribute_node->getNodeType() != xercesc::DOMNode::ATTRIBUTE_NODE) + { + continue; + } + + const xercesc::DOMAttr* const attribute = dynamic_cast(attribute_node); + if(attribute == nullptr) + { + G4Exception("G4GDMLRead::VectorRead()", "InvalidRead", FatalException, "No attribute found!"); + return; + } + const G4String attName = Transcode(attribute->getName()); + const G4String attValue = Transcode(attribute->getValue()); + + if(attName == "unit") + { + unit = G4UnitDefinition::GetValueOf(attValue); + } + else if(attName == "x") + { + vec.setX(eval.Evaluate(attValue)); + } + else if(attName == "y") + { + vec.setY(eval.Evaluate(attValue)); + } + else if(attName == "z") + { + vec.setZ(eval.Evaluate(attValue)); + } + else if (attName == "order" ){ + order = attValue; + } + } + + vec *= unit; + + for(auto rotate : rotations.at(order)) rotate(rot,vec); + +} + +G4LogicalVolume* remollGDMLReadStructure::FileRead(const xercesc::DOMElement* const fileElement) +{ + G4String name; + G4String volname; + + const xercesc::DOMNamedNodeMap* const attributes = fileElement->getAttributes(); + XMLSize_t attributeCount = attributes->getLength(); + + for(XMLSize_t attribute_index = 0; attribute_index < attributeCount; ++attribute_index) + { + xercesc::DOMNode* attribute_node = attributes->item(attribute_index); + + if(attribute_node->getNodeType() != xercesc::DOMNode::ATTRIBUTE_NODE) + { + continue; + } + + const xercesc::DOMAttr* const attribute = + dynamic_cast(attribute_node); + if(attribute == nullptr) + { + G4Exception("G4GDMLReadStructure::FileRead()", "InvalidRead", + FatalException, "No attribute found!"); + return nullptr; + } + const G4String attName = Transcode(attribute->getName()); + const G4String attValue = Transcode(attribute->getValue()); + + if(attName == "name") + { + name = attValue; + } + else if(attName == "volname") + { + volname = attValue; + } + } + + const G4bool isModule = true; + remollGDMLReadStructure structure; + structure.Read(name, validate, isModule); + + // Register existing auxiliar information defined in child module + // + const G4GDMLAuxMapType* aux = structure.GetAuxMap(); + if(!aux->empty()) + { + for(auto pos = aux->cbegin(); pos != aux->cend(); ++pos) + { + auxMap.insert(std::make_pair(pos->first, pos->second)); + } + } + + if(volname.empty()) + { + return structure.GetVolume(structure.GetSetup("Default")); + } + else + { + return structure.GetVolume(structure.GenerateName(volname)); + } +} + + +void remollGDMLReadStructure::PhysvolRead( const xercesc::DOMElement* const physvolElement, G4AssemblyVolume* pAssembly) +{ + G4String name; + G4LogicalVolume* logvol = nullptr; + G4AssemblyVolume* assembly = nullptr; + G4ThreeVector position(0.0, 0.0, 0.0); + G4ThreeVector rotation(0.0, 0.0, 0.0); + G4RotationMatrix orotm; + bool orotation_requested = false; + G4ThreeVector scale(1.0, 1.0, 1.0); + G4int copynumber = 0; + + const xercesc::DOMNamedNodeMap* const attributes = physvolElement->getAttributes(); + XMLSize_t attributeCount = attributes->getLength(); + + for(XMLSize_t attribute_index = 0; attribute_index < attributeCount; ++attribute_index) + { + xercesc::DOMNode* attribute_node = attributes->item(attribute_index); + + if(attribute_node->getNodeType() != xercesc::DOMNode::ATTRIBUTE_NODE) + { + continue; + } + + const xercesc::DOMAttr* const attribute = + dynamic_cast(attribute_node); + if(attribute == nullptr) + { + G4Exception("G4GDMLReadStructure::PhysvolRead()", "InvalidRead", + FatalException, "No attribute found!"); + return; + } + const G4String attName = Transcode(attribute->getName()); + const G4String attValue = Transcode(attribute->getValue()); + + if(attName == "name") + { + name = attValue; + } + if(attName == "copynumber") + { + copynumber = eval.EvaluateInteger(attValue); + } + } + + for(xercesc::DOMNode* iter = physvolElement->getFirstChild(); iter != nullptr; iter = iter->getNextSibling()) + { + if(iter->getNodeType() != xercesc::DOMNode::ELEMENT_NODE) + { + continue; + } + + const xercesc::DOMElement* const child = dynamic_cast(iter); + if(child == nullptr) + { + G4Exception("G4GDMLReadStructure::PhysvolRead()", "InvalidRead", FatalException, "No child found!"); + return; + } + const G4String tag = Transcode(child->getTagName()); + + if(tag == "volumeref") + { + const G4String& child_name = GenerateName(RefRead(child)); assembly = GetAssembly(child_name); + if(assembly == nullptr) + { + logvol = GetVolume(child_name); + } + } + else if(tag == "file") + { + logvol = FileRead(child); + } + else if(tag == "position") + { + VectorRead(child, position); + } + else if(tag == "rotation") + { + VectorRead(child, rotation); + } + else if(tag == "orotation") + { + RotationRead(child, orotm); + orotation_requested = true; + } + else if(tag == "scale") + { + VectorRead(child, scale); + } + else if(tag == "positionref") + { + position = GetPosition(GenerateName(RefRead(child))); + } + else if(tag == "rotationref") + { + rotation = GetRotation(GenerateName(RefRead(child))); + } + else if(tag == "scaleref") + { + scale = GetScale(GenerateName(RefRead(child))); + } + else + { + G4String error_msg = "Unknown tag in physvol: " + tag; + G4Exception("G4GDMLReadStructure::PhysvolRead()", "ReadError", + FatalException, error_msg); + return; + } + } + + G4Transform3D transform; + if(!orotation_requested){ + transform = G4Transform3D(GetRotationMatrix(rotation).inverse(), position); + } else { + transform = G4Transform3D(orotm.inverse(),position); + } + + transform = transform * G4Scale3D(scale.x(), scale.y(), scale.z()); + + if(pAssembly != nullptr) // Fill assembly structure + { + if(assembly != nullptr) // Case of recursive assemblies + { + pAssembly->AddPlacedAssembly(assembly, transform); + } + if(logvol == nullptr) + { + return; + } + pAssembly->AddPlacedVolume(logvol, transform); + } + else // Generate physical volume tree or do assembly imprint + { + if(assembly != nullptr) + { + assembly->MakeImprint(pMotherLogical, transform, 0, check); + } + else + { + if(logvol == nullptr) + { + return; + } + G4String pv_name = logvol->GetName() + "_PV"; + G4PhysicalVolumesPair pair = G4ReflectionFactory::Instance()->Place(transform, pv_name, logvol, pMotherLogical, false, copynumber, check); + + if(pair.first != nullptr) + { + GeneratePhysvolName(name, pair.first); + } + if(pair.second != nullptr) + { + GeneratePhysvolName(name, pair.second); + } + } + } +} + +void remollGDMLReadStructure::ParametersRead(const xercesc::DOMElement* const element) +{ + G4ThreeVector rotation(0.0, 0.0, 0.0); + G4ThreeVector position(0.0, 0.0, 0.0); + bool orotation_requested = false; + G4RotationMatrix orotm; + + G4GDMLParameterisation::PARAMETER parameter; + parameter.pRot = new G4RotationMatrix(); + + for(xercesc::DOMNode* iter = element->getFirstChild(); iter != nullptr; iter = iter->getNextSibling()) + { + if(iter->getNodeType() != xercesc::DOMNode::ELEMENT_NODE) + { + continue; + } + + const xercesc::DOMElement* const child = + dynamic_cast(iter); + if(child == nullptr) + { + G4Exception("G4GDMLReadParamvol::ParametersRead()", "InvalidRead", + FatalException, "No child found!"); + return; + } + const G4String tag = Transcode(child->getTagName()); + if(tag == "rotation") + { + VectorRead(child, rotation); + } + else if(tag == "orotation") + { + RotationRead(child,orotm); + orotation_requested = true; // makes orotation get priority over normal rotation. + } + else if(tag == "position") + { + VectorRead(child, position); + } + else if(tag == "positionref") + { + position = GetPosition(GenerateName(RefRead(child))); + } + else if(tag == "rotationref") + { + rotation = GetRotation(GenerateName(RefRead(child))); + } + else if(tag == "box_dimensions") + { + Box_dimensionsRead(child, parameter); + } + else if(tag == "trd_dimensions") + { + Trd_dimensionsRead(child, parameter); + } + else if(tag == "trap_dimensions") + { + Trap_dimensionsRead(child, parameter); + } + else if(tag == "tube_dimensions") + { + Tube_dimensionsRead(child, parameter); + } + else if(tag == "cone_dimensions") + { + Cone_dimensionsRead(child, parameter); + } + else if(tag == "sphere_dimensions") + { + Sphere_dimensionsRead(child, parameter); + } + else if(tag == "orb_dimensions") + { + Orb_dimensionsRead(child, parameter); + } + else if(tag == "torus_dimensions") + { + Torus_dimensionsRead(child, parameter); + } + else if(tag == "ellipsoid_dimensions") + { + Ellipsoid_dimensionsRead(child, parameter); + } + else if(tag == "para_dimensions") + { + Para_dimensionsRead(child, parameter); + } + else if(tag == "polycone_dimensions") + { + Polycone_dimensionsRead(child, parameter); + } + else if(tag == "polyhedra_dimensions") + { + Polyhedra_dimensionsRead(child, parameter); + } + else if(tag == "hype_dimensions") + { + Hype_dimensionsRead(child, parameter); + } + else + { + G4String error_msg = "Unknown tag in parameters: " + tag; + G4Exception("G4GDMLReadParamvol::ParametersRead()", "ReadError", + FatalException, error_msg); + } + } + + + if(orotation_requested){ + parameter.pRot = std::move(&orotm); + } else { + parameter.pRot->rotateX(rotation.x()); + parameter.pRot->rotateY(rotation.y()); + parameter.pRot->rotateZ(rotation.z()); + } + + + parameter.position = position; + + parameterisation->AddParameter(parameter); + +} + +void remollGDMLReadStructure::Volume_contentRead( const xercesc::DOMElement* const volumeElement) +{ + for(xercesc::DOMNode* iter = volumeElement->getFirstChild(); iter != nullptr; iter = iter->getNextSibling()) + { + if(iter->getNodeType() != xercesc::DOMNode::ELEMENT_NODE) + { + continue; + } + + const xercesc::DOMElement* const child = + dynamic_cast(iter); + if(child == nullptr) + { + G4Exception("G4GDMLReadStructure::Volume_contentRead()", "InvalidRead", + FatalException, "No child found!"); + return; + } + const G4String tag = Transcode(child->getTagName()); + + if((tag == "auxiliary") || (tag == "materialref") || (tag == "solidref")) + { + // These are already processed in VolumeRead() + } + else if(tag == "paramvol") + { + ParamvolRead(child, pMotherLogical); + } + else if(tag == "physvol") + { + PhysvolRead(child); + } + else if(tag == "replicavol") + { + G4int number = 1; + const xercesc::DOMNamedNodeMap* const attributes = child->getAttributes(); + XMLSize_t attributeCount = attributes->getLength(); + for(XMLSize_t attribute_index = 0; attribute_index < attributeCount; ++attribute_index) + { + xercesc::DOMNode* attribute_node = attributes->item(attribute_index); + if(attribute_node->getNodeType() != xercesc::DOMNode::ATTRIBUTE_NODE) + { + continue; + } + const xercesc::DOMAttr* const attribute = + dynamic_cast(attribute_node); + if(attribute == nullptr) + { + G4Exception("G4GDMLReadStructure::Volume_contentRead()", + "InvalidRead", FatalException, "No attribute found!"); + return; + } + const G4String attName = Transcode(attribute->getName()); + const G4String attValue = Transcode(attribute->getValue()); + if(attName == "number") + { + number = eval.EvaluateInteger(attValue); + } + } + ReplicavolRead(child, number); + } + else if(tag == "divisionvol") + { + DivisionvolRead(child); + } + else if(tag == "loop") + { + LoopRead(child, &G4GDMLRead::Volume_contentRead); + } + else + { + G4cout << "Treating unknown GDML tag in volume '" << tag << "' as GDML extension..." << G4endl; + } + } +}