getGeoInformation.py

Go to the documentation of this file.

#!/usr/bin/env python
# SPDX-License-Identifier: LGPL-3.0-or-later
# SPDX-FileCopyrightText: Copyright CERN for the benefit of the SHiP Collaboration
 
# prints z-coordinators of SHiP detector volumes
# WARNING: printing the entire geometry takes a lot of time
# 24-02-2015 comments to EvH
 
import operator
from argparse import ArgumentParser
from array import array
 
import ROOT
 
 
def local2Global(n):
    Info = {}
    nav = ROOT.gGeoManager.GetCurrentNavigator()
    nav.cd(n)
    Info["node"] = nav.GetCurrentNode()
    Info["path"] = n
    tmp = Info["node"].GetVolume().GetShape()
    Info["material"] = Info["node"].GetVolume().GetMaterial().GetName()
    if options.moreInfo:
        x = ROOT.gGeoManager.GetVerboseLevel()
        ROOT.gGeoManager.SetVerboseLevel(0)
        Info["weight"] = Info["node"].GetVolume().Weight()  # kg
        Info["cubicmeter"] = Info["node"].GetVolume().Capacity() / 1000000.0  #
        ROOT.gGeoManager.SetVerboseLevel(x)
    o = [tmp.GetOrigin()[0], tmp.GetOrigin()[1], tmp.GetOrigin()[2]]
    Info["locorign"] = o
    local = array("d", o)
    globOrigin = array("d", [0, 0, 0])
    nav.LocalToMaster(local, globOrigin)
    Info["origin"] = globOrigin
    shifts = [
        [-tmp.GetDX() + o[0], o[1], o[2]],
        [tmp.GetDX() + o[0], o[1], o[2]],
        [o[0], -tmp.GetDY() + o[1], o[2]],
        [o[0], tmp.GetDY() + o[1], o[2]],
        [o[0], o[1], -tmp.GetDZ() + o[2]],
        [o[0], o[1], tmp.GetDZ() + o[2]],
    ]
    shifted = []
    for s in shifts:
        local = array("d", s)
        glob = array("d", [0, 0, 0])
        nav.LocalToMaster(local, glob)
        shifted.append([glob[0], glob[1], glob[2]])
    Info["boundingbox"] = {}
    for j in range(3):
        jmin = 1e30
        jmax = -1e30
        for s in shifted:
            if s[j] < jmin:
                jmin = s[j]
            if s[j] > jmax:
                jmax = s[j]
        Info["boundingbox"][j] = [jmin, jmax]
    return Info
 
 
def print_info(path: str, node, level: int, currentlevel: int, print_sub_det_info=False) -> None:
    sub_nodes = {}
    fullInfo = {}
    for subnode in node.GetNodes():
        name = subnode.GetName()
        fullInfo[name] = local2Global(path + "/" + name)
        sub_nodes[name] = fullInfo[name]["origin"][2]
 
    for name, _ in sorted(sub_nodes.items(), key=operator.itemgetter(1)):
        boundingbox = fullInfo[name]["boundingbox"]
 
        format_string = (
            "{:<28s}: z={:10.4F}cm  dZ={:10.4F}cm  [{:10.4F}   {:10.4F}]"
            + " dx={:10.4F}cm [{:10.4F}   {:10.4F}] dy={:10.4F}cm [{:10.4F}   {:10.4F}] {:>20s}"
        )
 
        format_variable = [
            "   " * int(currentlevel) + name,
            fullInfo[name]["origin"][2],
            abs(boundingbox[2][0] - boundingbox[2][1]) / 2.0,
            boundingbox[2][0],
            boundingbox[2][1],
            abs(boundingbox[0][0] - boundingbox[0][1]) / 2.0,
            boundingbox[0][0],
            boundingbox[0][1],
            abs(boundingbox[1][0] - boundingbox[1][1]) / 2.0,
            boundingbox[1][0],
            boundingbox[1][1],
            fullInfo[name]["material"],
        ]
 
        if options.moreInfo:
            cubicmeter = fullInfo[name]["cubicmeter"]
            weight = fullInfo[name]["weight"]
            format_string += " {:10.1F}kg {:10.1F}m3"
            format_variable.extend([weight, cubicmeter])
 
        print(format_string.format(*format_variable))
 
        if options.volume in ["", name]:
            print_sub_det_info = True
 
        if print_sub_det_info and currentlevel < level and fullInfo[name]["node"].GetNodes():
            print_info(fullInfo[name]["path"], fullInfo[name]["node"], level, currentlevel + 1, print_sub_det_info)
 
        if currentlevel == 0:
            print_sub_det_info = False
 
 
parser = ArgumentParser()
parser.add_argument("-g", "--geometry", help="Input geometry file", required=True)
parser.add_argument("-l", "--level", help="Max subnode level", default=0)
parser.add_argument("-v", "--volume", help="Name of node to expand", default="")
parser.add_argument("-X", "--moreInfo", help="Print weight and capacity", action="store_true")
 
options = parser.parse_args()
fname = options.geometry
fgeom = ROOT.TFile.Open(fname)
fGeo = fgeom.Get("FAIRGeom")
top = fGeo.GetTopVolume()
 
 
if options.moreInfo:
    print(
        "   Detector element             z(midpoint)     halflength       volume-start volume-end   dx"
        "                x-start       x-end       dy                y-start       y-end         material          weight  capacity"
    )
else:
    print(
        "   Detector element             z(midpoint)     halflength       volume-start volume-end   dx"
        "                x-start       x-end       dy                y-start       y-end         material"
    )
 
currentlevel = 0
print_info("", top, int(options.level), currentlevel)