diff --git a/inc/TRestGeant4GeometryInfo.h b/inc/TRestGeant4GeometryInfo.h index 56b617a..8f2b7ed 100644 --- a/inc/TRestGeant4GeometryInfo.h +++ b/inc/TRestGeant4GeometryInfo.h @@ -2,6 +2,7 @@ #ifndef REST_TRESTGEANT4GEOMETRYINFO_H #define REST_TRESTGEANT4GEOMETRYINFO_H +#include #include #include @@ -12,10 +13,11 @@ class G4VPhysicalVolume; class TRestGeant4GeometryInfo { - ClassDef(TRestGeant4GeometryInfo, 3); + ClassDef(TRestGeant4GeometryInfo, 4); private: bool fIsAssembly = false; + TString fPathSeparator = "_"; std::map fVolumeNameMap = {}; std::map fVolumeNameReverseMap = {}; @@ -34,12 +36,14 @@ class TRestGeant4GeometryInfo { std::vector fGdmlLogicalNames; std::map - fGeant4PhysicalNameToNewPhysicalNameMap; /* - * only makes sense when using assembly - */ + fNewPhysicalToGeant4PhysicalNameMap; // reverse map of fGeant4PhysicalNameToNewPhysicalNameMap + + std::map fGeant4AssemblyImprintToGdmlNameMap; + std::map> fGdmlAssemblyToChildrenGeant4ToGdmlPhysicalNameMap; + std::map fGeant4AssemblyImprintToAssemblyLogicalNameMap; std::map fPhysicalToLogicalVolumeMap; - std::map > fLogicalToPhysicalMap; + std::map> fLogicalToPhysicalMap; // many physical volumes can point to one single logical std::map fLogicalToMaterialMap; std::map fPhysicalToPositionInWorldMap; @@ -49,7 +53,9 @@ class TRestGeant4GeometryInfo { void PopulateFromGdml(const TString&); + TString GetAlternativePathFromGeant4Path(const TString&) const; TString GetAlternativeNameFromGeant4PhysicalName(const TString&) const; + std::set GetAlternativeNamesFromGeant4PhysicalName(const TString&) const; TString GetGeant4PhysicalNameFromAlternativeName(const TString&) const; std::vector GetAllLogicalVolumes() const; @@ -86,7 +92,8 @@ class TRestGeant4GeometryInfo { } inline bool IsAssembly() const { return fIsAssembly; } - + inline TString GetPathSeparator() const { return fPathSeparator; } + void SetPathSeparator(const TString& separator) { fPathSeparator = separator; } void InsertVolumeName(Int_t id, const TString& volumeName); TString GetVolumeFromID(Int_t id) const; diff --git a/src/TRestGeant4GeometryInfo.cxx b/src/TRestGeant4GeometryInfo.cxx index 48dfca9..1b72569 100644 --- a/src/TRestGeant4GeometryInfo.cxx +++ b/src/TRestGeant4GeometryInfo.cxx @@ -9,6 +9,8 @@ #include +#include "TRestStringHelper.h" + using namespace std; namespace myXml { @@ -34,9 +36,22 @@ TString GetNodeAttribute(TXMLEngine xml, XMLNodePointer_t node, const TString& a } return {}; } +XMLNodePointer_t GetChildByAttributeValue(TXMLEngine xml, XMLNodePointer_t node, const TString& attributeName, + const TString& attributeValue) { + XMLNodePointer_t child = xml.GetChild(node); + while (child) { + TString childAttributeValue = GetNodeAttribute(xml, child, attributeName); + if (childAttributeValue.EqualTo(attributeValue)) { + return child; + } + child = xml.GetNext(child); + } + return nullptr; +} void AddVolumesRecursively(vector* physicalNames, vector* logicalNames, const vector& children, map& nameTable, - map>& childrenTable, const TString& name = "") { + map>& childrenTable, const TString& name = "", + const TString& separator = "_") { if (children.empty()) { physicalNames->push_back(name); const auto logicalVolumeName = nameTable[name]; @@ -44,10 +59,10 @@ void AddVolumesRecursively(vector* physicalNames, vector* logi return; } for (const auto& childName : children) { - const auto newName = name + "_" + childName; + const auto newName = name + separator + childName; nameTable[newName] = nameTable[childName]; // needed to retrieve logical volume name AddVolumesRecursively(physicalNames, logicalNames, childrenTable[nameTable[childName]], nameTable, - childrenTable, newName); + childrenTable, newName, separator); } } } // namespace myXml @@ -65,25 +80,172 @@ void TRestGeant4GeometryInfo::PopulateFromGdml(const TString& gdmlFilename) { << " not found" << endl; exit(1); } - map nameTable; - map> childrenTable; XMLNodePointer_t mainNode = xml.DocGetRootElement(xmldoc); XMLNodePointer_t structure = myXml::FindChildByName(xml, mainNode, "structure"); XMLNodePointer_t child = xml.GetChild(structure); - while (child) { - TString name = xml.GetNodeName(child); + + /* When a PV is placed from an assembly, its daughter physical volumes are imprinted into the mother + volume where you are placing the assembly. This daughter PV are named following the format: + "av_WWW_impr_XXX_YYY_ZZZ". This first loop over the gdml structure is to get the map "assemblyName" -> + "av_WWW" + */ + size_t assemblyCounter = 0; // track the WWW + map gdmlToGeant4AssemblyNameMap; // e.g. "assemblyName" -> "av_1" + map gdmlAssemblyNameToImprintCounterMap; // to track the number of imprints per assembly + RESTDebug << "Searching for assemblies..." << RESTendl; + while (child) { // loop over the direct children of structure (logical volumes and assemblies) + TString name = xml.GetNodeName(child); // name of the type of children: "volume" or "assembly" + if (!name.EqualTo("assembly")) { + child = xml.GetNext(child); + continue; + } + TString assemblyName = myXml::GetNodeAttribute(xml, child, "name"); + gdmlToGeant4AssemblyNameMap[assemblyName] = + "av_" + to_string(++assemblyCounter); // first assembly is av_1 + gdmlAssemblyNameToImprintCounterMap[assemblyName] = 0; // initialize with the assembly found + RESTDebug << "Assembly found: " << assemblyName << " → " << gdmlToGeant4AssemblyNameMap[assemblyName] + << RESTendl; + + /* The daughter physical volumes defined after a nested assembly gets the imprint number XXX+1 + instead of the XXX of its mother assembly imprint (which I think is a bug in Geant4). + To avoid having several "av_WWW_impr_XXX+n" pointing to different GDML names, + you must define all the daughter PV from normal LV before all the nested assembly PV.*/ + bool hasNestedAssemblies = false; + std::map childrenGeant4toGdmlMap; + size_t childrenPVCounter = 0; + auto physicalVolumeNode = + xml.GetChild(child); // children of a volume or assembly are physical volumes + while (physicalVolumeNode) { + auto physicalVolumeName = myXml::GetNodeAttribute(xml, physicalVolumeNode, "name"); + auto volumeRefNode = + xml.GetChild(physicalVolumeNode); // this are volumeref, position and rotation + while (volumeRefNode) { + TString volumeRefNodeName = + xml.GetNodeName(volumeRefNode); // "volumeref", "position" or "rotation" + if (volumeRefNodeName.EqualTo("volumeref")) { + TString refName = + myXml::GetNodeAttribute(xml, volumeRefNode, "ref"); // the logical volume name + TString geant4Name = + refName + "_pv_" + + to_string(childrenPVCounter++); // first pv is logicalVolumeName_pv_0 + childrenGeant4toGdmlMap[geant4Name] = physicalVolumeName; + if (gdmlToGeant4AssemblyNameMap.count(refName) > 0) { + hasNestedAssemblies = true; + } else { + if (hasNestedAssemblies) { + RESTError << "TRestGeant4GeometryInfo::PopulateFromGdml - Assembly '" + << assemblyName << "' contains physical volumes from normal " + << "(i.e. non-assembly) logical volumes defined after physical volumes " + << "from assemblies. Due to a Geant4 bug you cannot do this. " + << "Please define the physical volumes from the assemblies last." + << RESTendl; + } + } + } + volumeRefNode = xml.GetNext(volumeRefNode); + } + physicalVolumeNode = xml.GetNext(physicalVolumeNode); + } + fGdmlAssemblyToChildrenGeant4ToGdmlPhysicalNameMap[assemblyName] = childrenGeant4toGdmlMap; + child = xml.GetNext(child); + } + + /*Recursive function to obtain the prefix 'av_WWW_impr_XXX' of the daughters of the assemblies imprints. + When a PV is placed from an assembly, its daughter physical volumes are imprinted into the mother + volume where you are placing the assembly. This daughter PV are named following the format: + "av_WWW_impr_XXX_YYY_ZZZ". But, if one of those daughter PV is itself an assembly, its own daughter PV are + named (wrongly in my opinion) "av_WWW_impr_XXX+1_yyy_zzz". This behavior is propagated down the chain + to the final child assembly which does not have any daughter assembly. So, the godFatherAssembly is the + highest assembly which begins this chain and its "av_WWW" is used for all its consecutive assembly + children imprints and each of this assembly children adds +1 to the imprint number of the + godFatherAssembly.*/ + std::function ProcessNestedAssembliesRecursively = + [&](const XMLNodePointer_t parentNode, const TString godFatherAssemblyName, const TString pathSoFar) { + auto physicalVolumeNode = xml.GetChild(parentNode); + // godFatherAssemblyName = parentNode logical name // the highest assembly volume in the nested + // chain + while (physicalVolumeNode) { + auto physicalVolumeName = myXml::GetNodeAttribute(xml, physicalVolumeNode, "name"); + auto volumeRefNode = + xml.GetChild(physicalVolumeNode); // this are volumeref, position and rotation + while (volumeRefNode) { + TString volumeRefNodeName = + xml.GetNodeName(volumeRefNode); // "volumeref", "position" or "rotation" + if (volumeRefNodeName.EqualTo("volumeref")) { + TString refName = + myXml::GetNodeAttribute(xml, volumeRefNode, "ref"); // the logical volume name + if (gdmlToGeant4AssemblyNameMap.count(refName) > 0) { + // it's an assembly + TString newGodFatherAssemblyName = godFatherAssemblyName; + if (newGodFatherAssemblyName.IsNull()) { + // start assembly children chain with this assembly as godFather + newGodFatherAssemblyName = refName; + } + size_t imprintCounter = + ++gdmlAssemblyNameToImprintCounterMap[newGodFatherAssemblyName]; + TString imprint = gdmlToGeant4AssemblyNameMap[newGodFatherAssemblyName] + + "_impr_" + to_string(imprintCounter); + TString path = + pathSoFar + (pathSoFar.IsNull() ? "" : fPathSeparator) + physicalVolumeName; + fGeant4AssemblyImprintToGdmlNameMap[imprint] = path; + // Continue the assembly children chain with its correspondant godFatherAssembly + // and path + auto assemblyNode = + myXml::GetChildByAttributeValue(xml, structure, "name", refName); + ProcessNestedAssembliesRecursively(assemblyNode, newGodFatherAssemblyName, path); + } else { + // its a regular logical volume + // Regular children resets the godFatherAssembly and path + auto assemblyNode = + myXml::GetChildByAttributeValue(xml, structure, "name", refName); + ProcessNestedAssembliesRecursively(assemblyNode, "", ""); + } + } + volumeRefNode = xml.GetNext(volumeRefNode); + } + physicalVolumeNode = xml.GetNext(physicalVolumeNode); + } + }; + + // We loop a second time over the gdml structure to get the imprint of each assembly + // into fGeant4AssemblyImprintToGdmlNameMap: e.g. "av_2_impr_5" -> "shielding/vessel" + auto worldNode = myXml::GetChildByAttributeValue(xml, structure, "name", "world"); + if (!worldNode) { + worldNode = myXml::GetChildByAttributeValue(xml, structure, "name", "World"); + if (!worldNode) { + cout << "Could not find world volume in GDML, please name it either 'World' or 'world'" << endl; + exit(1); + } + } + TString godFatherAssemblyName = ""; // the highest assembly volume in the nested chain + TString pathSoFar = ""; // the path for the nested assemblies + ProcessNestedAssembliesRecursively(worldNode, godFatherAssemblyName, pathSoFar); + + // We loop a third time over gdml structure to get the physical and logical volumes names + map nameTable; + map> childrenTable; + child = xml.GetChild(structure); + while (child) { // loop over the direct children of structure (logical volumes and assemblies) + TString name = xml.GetNodeName(child); // name of the type of children: "volume" or "assembly" TString volumeName = myXml::GetNodeAttribute(xml, child, "name"); - auto physicalVolumeNode = xml.GetChild(child); + auto physicalVolumeNode = + xml.GetChild(child); // children of a volume or assembly are physical volumes childrenTable[volumeName] = {}; while (physicalVolumeNode) { auto physicalVolumeName = myXml::GetNodeAttribute(xml, physicalVolumeNode, "name"); - auto volumeRefNode = xml.GetChild(physicalVolumeNode); + auto volumeRefNode = + xml.GetChild(physicalVolumeNode); // this are volumeref, position and rotation while (volumeRefNode) { - TString volumeRefNodeName = xml.GetNodeName(volumeRefNode); + TString volumeRefNodeName = + xml.GetNodeName(volumeRefNode); // "volumeref", "position" or "rotation" if (volumeRefNodeName.EqualTo("volumeref")) { - TString refName = myXml::GetNodeAttribute(xml, volumeRefNode, "ref"); + TString refName = + myXml::GetNodeAttribute(xml, volumeRefNode, "ref"); // the logical volume name nameTable[physicalVolumeName] = refName; childrenTable[volumeName].push_back(physicalVolumeName); + if (gdmlToGeant4AssemblyNameMap.count(refName) > 0) { + fGeant4AssemblyImprintToAssemblyLogicalNameMap[physicalVolumeName] = refName; + } } volumeRefNode = xml.GetNext(volumeRefNode); } @@ -106,7 +268,7 @@ void TRestGeant4GeometryInfo::PopulateFromGdml(const TString& gdmlFilename) { for (const auto& topName : childrenTable[worldVolumeName]) { auto children = childrenTable[nameTable[topName]]; myXml::AddVolumesRecursively(&fGdmlNewPhysicalNames, &fGdmlLogicalNames, children, nameTable, - childrenTable, topName); + childrenTable, topName, fPathSeparator); } xml.FreeDoc(xmldoc); @@ -128,22 +290,96 @@ void TRestGeant4GeometryInfo::PopulateFromGdml(const TString& gdmlFilename) { << endl; exit(1); } + + /* + // print gdmlToGeant4AssemblyNameMap + cout << "GDML to Geant4 Assembly Name Map:" << endl; + for (const auto& kv : gdmlToGeant4AssemblyNameMap) { + cout << "\t" << kv.first << " → " << kv.second << endl; + } + + //print gdmlAssemblyNameToImprintCounterMap + cout << "GDML Assembly Name to Imprint Counter Map:" << endl; + for (const auto& kv : gdmlAssemblyNameToImprintCounterMap) { + cout << "\t" << kv.first << " → " << kv.second << endl; + } + + //print fGeant4AssemblyImprintToGdmlNameMap + cout << "GDML Assembly Imprint to GDML Name Map:" << endl; + for (const auto& kv : fGeant4AssemblyImprintToGdmlNameMap) { + cout << "\t" << kv.first << " → " << kv.second << endl; + } + + cout << "GDML assembly children Geant4 to GDML Physical Name Map:" << endl; + for (const auto& kv : fGdmlAssemblyToChildrenGeant4ToGdmlPhysicalNameMap) { + cout << "\tAssembly: " << kv.first << endl; + for (const auto& childKv : kv.second) { + cout << "\t\t" << childKv.first << " → " << childKv.second << endl; + } + } + + cout << "Assembly imprint to Assembly Logical Name Map:" << endl; + for (const auto& kv : fGeant4AssemblyImprintToAssemblyLogicalNameMap) { + cout << "\t" << kv.first << " → " << kv.second << endl; + } + */ +} + +////////////////////////////////////////////////////////////////////////// +/// \brief Converts a Geant4 volume path to a GDML volume path. The path is +/// the concatenation of the names of the nested volumes from the world to the desired volume. +/// This method is meant mainly meant to handle the conversion of assembly imprints names to the +/// GDML assembly names used in the GDML file. +/// \param geant4Path The Geant4 volume path. +/// \return The corresponding GDML volume path using the PV names defined in the GDML file. +TString TRestGeant4GeometryInfo::GetAlternativePathFromGeant4Path(const TString& geant4Path) const { + std::string convertedPath = geant4Path.Data(); + // convert the Geant4 assembly imprint convention to GDML name + // e.g. av_1_impr_2_childLV_pv_1 → assemblyName/childLV_pv_1 + for (const auto& [gvName, gdmlName] : fGeant4AssemblyImprintToGdmlNameMap) { + convertedPath = Replace(convertedPath, (gvName + "_").Data(), (gdmlName + fPathSeparator).Data()); + } + // convert the children names inside assemblies to physical volume names in GDML + // e.g. assemblyName/childLV_pv_1 → assemblyName/childPVname + for (const auto& [gvImprint, assemblyLogicalName] : fGeant4AssemblyImprintToAssemblyLogicalNameMap) { + auto gvImprintPosition = convertedPath.find(gvImprint.Data()); + if (gvImprintPosition != std::string::npos) { + for (const auto& [childGeant4Name, childGdmlName] : + fGdmlAssemblyToChildrenGeant4ToGdmlPhysicalNameMap.at(assemblyLogicalName)) { + TString toReplace = gvImprint + fPathSeparator + childGeant4Name; + TString replaceBy = gvImprint + fPathSeparator + childGdmlName; + convertedPath = Replace(convertedPath, toReplace.Data(), replaceBy.Data()); + } + } + } + return TString(convertedPath); } TString TRestGeant4GeometryInfo::GetAlternativeNameFromGeant4PhysicalName( const TString& geant4PhysicalName) const { - if (fGeant4PhysicalNameToNewPhysicalNameMap.count(geant4PhysicalName) > 0) { - return fGeant4PhysicalNameToNewPhysicalNameMap.at(geant4PhysicalName); + for (const auto& kv : fNewPhysicalToGeant4PhysicalNameMap) { + if (kv.second.EqualTo(geant4PhysicalName)) { + return kv.first; + } } return geant4PhysicalName; } +set TRestGeant4GeometryInfo::GetAlternativeNamesFromGeant4PhysicalName( + const TString& geant4PhysicalName) const { + set alternativeNames; + for (const auto& kv : fNewPhysicalToGeant4PhysicalNameMap) { + if (kv.second.EqualTo(geant4PhysicalName)) { + alternativeNames.insert(kv.first); + } + } + return alternativeNames; +} + TString TRestGeant4GeometryInfo::GetGeant4PhysicalNameFromAlternativeName( const TString& alternativeName) const { - for (const auto& kv : fGeant4PhysicalNameToNewPhysicalNameMap) { - if (kv.second == alternativeName) { - return kv.first; - } + if (fNewPhysicalToGeant4PhysicalNameMap.count(alternativeName) > 0) { + return fNewPhysicalToGeant4PhysicalNameMap.at(alternativeName); } return ""; } @@ -222,8 +458,8 @@ std::vector TRestGeant4GeometryInfo::GetAllPhysicalVolumes() const { std::vector TRestGeant4GeometryInfo::GetAllAlternativePhysicalVolumes() const { auto volumes = std::vector(); - for (const auto& kv : fPhysicalToLogicalVolumeMap) { - volumes.emplace_back(GetAlternativeNameFromGeant4PhysicalName(kv.first)); + for (const auto& kv : fNewPhysicalToGeant4PhysicalNameMap) { + volumes.emplace_back(kv.first); } return volumes; diff --git a/src/TRestGeant4Metadata.cxx b/src/TRestGeant4Metadata.cxx index 55ba277..04eb90a 100644 --- a/src/TRestGeant4Metadata.cxx +++ b/src/TRestGeant4Metadata.cxx @@ -829,6 +829,9 @@ void TRestGeant4Metadata::InitFromConfigFile() { } fGeometryPath = GetParameter("geometryPath", ""); + auto volPathSeparator = + GetParameter("volPathSeparator", "_"); // for TRestGeant4GeometryInfo::fPathSeparator + fGeant4GeometryInfo.SetPathSeparator(volPathSeparator); fStoreHadronicTargetInfo = StringToBool(GetParameter("storeHadronicTargetInfo", "false")); @@ -1387,6 +1390,8 @@ void TRestGeant4Metadata::ReadDetector() { vector physicalVolumes; if (!fGeant4GeometryInfo.IsValidGdmlName(name)) { if (fGeant4GeometryInfo.IsValidLogicalVolume(name)) { + RESTDebug << "Volume name '" << name << "' is a valid logical volume. " + << "Inserting all physical volumes from it." << RESTendl; for (const auto& physical : fGeant4GeometryInfo.GetAllPhysicalVolumesFromLogical(name)) { physicalVolumes.emplace_back( fGeant4GeometryInfo.GetAlternativeNameFromGeant4PhysicalName(physical)); @@ -1394,22 +1399,32 @@ void TRestGeant4Metadata::ReadDetector() { } // does not match a explicit (gdml) physical name or a logical name, perhaps its a regular exp if (physicalVolumes.empty()) { + RESTDebug << "Volume name '" << name << "' is not a valid logical volume. " + << "Trying to match as regular expression for physical volumes." << RESTendl; for (const auto& physical : fGeant4GeometryInfo.GetAllPhysicalVolumesMatchingExpression(name)) { + RESTExtreme << "Volume name '" << name << "' matches physical volume '" << physical << "'" + << RESTendl; physicalVolumes.emplace_back( fGeant4GeometryInfo.GetAlternativeNameFromGeant4PhysicalName(physical)); } } if (physicalVolumes.empty()) { + RESTDebug << "Volume name '" << name + << "' is not a valid logical volume neither physical volume regex. " + << "Trying to match as regular expression for logical volumes." << RESTendl; for (const auto& logical : fGeant4GeometryInfo.GetAllLogicalVolumesMatchingExpression(name)) { for (const auto& physical : fGeant4GeometryInfo.GetAllPhysicalVolumesFromLogical(logical)) { + RESTExtreme << "Volume name '" << name << "' matches logical volume '" << logical + << "' with physical volume '" << physical << "'" << RESTendl; physicalVolumes.emplace_back( fGeant4GeometryInfo.GetAlternativeNameFromGeant4PhysicalName(physical)); } } } } else { + RESTDebug << "Volume name '" << name << "' is a valid physical volume." << RESTendl; physicalVolumes.push_back(name); } @@ -1450,7 +1465,8 @@ void TRestGeant4Metadata::ReadDetector() { for (const auto& physical : physicalVolumes) { RESTInfo << "Adding " << (isSensitive ? "sensitive" : "active") << " volume from RML: '" - << physical << (chance != 1 ? " with change: " + to_string(chance) : " ") << RESTendl; + << physical << "'" << (chance != 1 ? " with chance: " + to_string(chance) : " ") + << RESTendl; SetActiveVolume(physical, chance, maxStep); if (isSensitive) { InsertSensitiveVolume(physical);