eic-opticks/analytic/gdml.py

0001 #!/usr/bin/env python
0002 #
0003 # Copyright (c) 2019 Opticks Team. All Rights Reserved.
0004 #
0005 # This file is part of Opticks
0006 # (see https://bitbucket.org/simoncblyth/opticks).
0007 #
0008 # Licensed under the Apache License, Version 2.0 (the "License");
0009 # you may not use this file except in compliance with the License.
0010 # You may obtain a copy of the License at
0011 #
0012 #   http://www.apache.org/licenses/LICENSE-2.0
0013 #
0014 # Unless required by applicable law or agreed to in writing, software
0015 # distributed under the License is distributed on an "AS IS" BASIS,
0016 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
0017 # See the License for the specific language governing permissions and
0018 # limitations under the License.
0019 #
0020
0021 """
0022 gdml.py : parsing GDML
0023 =========================
0024
0025
0026 """
0027 import os, re, logging, math, collections, argparse
0028
0029 log = logging.getLogger(__name__)
0030
0031 #from opticks.ana.main import opticks_main
0032
0033 from opticks.ana.key import key_
0034 from opticks.ana.base import u_, b_, d_
0035 from opticks.ana.nbase import find_ranges, np_fromstring, np_tostring, np_digest
0036 from opticks.ana.shape import X, SEllipsoid, STubs, STorus, SCons, SPolycone, SSubtractionSolid, SUnionSolid, SIntersectionSolid
0037
0038 from opticks.analytic.csg import CSG
0039 from opticks.analytic.treebuilder import TreeBuilder
0040 from opticks.analytic.glm import make_trs, make_transform
0041 from opticks.analytic.prism import make_trapezoid
0042
0043
0044 import numpy as np
0045 import lxml.etree as ET
0046 import lxml.html as HT
0047
0048
0049 tostring_ = lambda _:ET.tostring(_)
0050 exists_ = lambda _:os.path.exists(os.path.expandvars(_))
0051 parse_ = lambda _:ET.parse(os.path.expandvars(_)).getroot()
0052 fparse_ = lambda _:HT.fragments_fromstring(file(os.path.expandvars(_)).read())
0053 pp_ = lambda d:"\n".join([" %30s : %f " % (k,d[k]) for k in sorted(d.keys())])
0054
0055 unref_ = lambda n:n[:-9] if n[-9:-7] == '0x' else n   # HamamatsuR12860_PMT_20inch_body_solid_1_90x32aa10 trim the 0x32aa10
0056
0057
0058 def construct_transform(obj):
0059     tla = obj.position.xyz if obj.position is not None else None
0060     rot = obj.rotation.xyz if obj.rotation is not None else None
0061     sca = obj.scale.xyz if obj.scale is not None else None
0062     order = "trs"
0063
0064     #elem = filter(None, [tla,rot,sca])
0065     #if len(elem) > 1:
0066     #    log.warning("construct_transform multi %s " % repr(obj))
0067     #pass
0068
0069     return make_transform( order, tla, rot, sca , three_axis_rotate=True, transpose_rotation=True, suppress_identity=False, dtype=np.float32 )
0070
0071
0072 class G(object):
0073     """
0074     Base wrapper type for lxml parsed GDML elements
0075     """
0076     pvtype = 'Physvol'
0077     lvtype = 'Volume'
0078     postype = 'position'
0079
0080     typ = property(lambda self:self.__class__.__name__)
0081     name  = property(lambda self:self.elem.attrib.get('name',None))
0082     xml = property(lambda self:tostring_(self.elem))
0083     gidx = property(lambda self:"[%d]" % self.idx if hasattr(self, 'idx') else '[]' )  # materials, volumes and top level solids have idx attribute
0084
0085
0086
0087     is_primitive = property(lambda self:issubclass(self.__class__ , Primitive ))
0088     is_boolean   = property(lambda self:issubclass(self.__class__ , Boolean ))
0089
0090     def _get_shortname(self):
0091         """/dd/Geometry/PMT/lvPmtHemi0xc133740 -> lvPmtHemi"""
0092         base = self.name.split("/")[-1]
0093         return unref_(base)
0094     shortname = property(_get_shortname)
0095
0096     def att(self, name, default=None, typ=None):
0097         assert typ is not None
0098         v = self.elem.attrib.get(name, default)
0099         return typ(v)
0100
0101     def __init__(self, elem, g=None):
0102         """
0103         :param elem: lxml element
0104         :param g:
0105         """
0106         self.elem = elem
0107         self.g = g
0108
0109     def findall_(self, expr):
0110         """
0111         lxml findall result elements are wrapped in the class appropriate to their tags,
0112         note the global g gets passed into all Elem
0113
0114         g.kls is a dict associating tag names to classes, Elem
0115         is fallback, all the classes have signature of  elem-instance, g
0116
0117         """
0118         wrap_ = lambda e:self.g.kls.get(e.tag,G)(e,self.g)
0119         fa = list(map(wrap_, self.elem.xpath(expr) ))    # formerly used findall, moved to xpath for selective expressions
0120         kln = self.__class__.__name__
0121         name = self.name
0122         log.debug("findall_ from %s:%s expr:%s returned %s " % (kln, name, expr, len(fa)))
0123         return fa
0124
0125     def findone_(self, expr):
0126         """simple dumb lxml find"""
0127         all_ = self.findall_(expr)
0128         assert len(all_) == 1
0129         return all_[0]
0130
0131     def find1_(self, expr):
0132         all_ = self.findall_(expr)
0133         assert len(all_) in [0,1]
0134         return all_[0] if len(all_) == 1 else None
0135
0136     def find_(self, expr):
0137         e = self.elem.find(expr)
0138         wrap_ = lambda e:self.g.kls.get(e.tag,G)(e,self.g)
0139         return wrap_(e) if e is not None else None
0140
0141     def __repr__(self):
0142         return "path %s path_label %s " % (self.path, self.path_label)
0143
0144     def __repr__0(self):
0145         return "%15s : %s " % ( self.elem.tag, repr(self.elem.attrib) )
0146
0147     def __str__(self):
0148         return self.xml
0149
0150
0151
0152 class Material(G):
0153     state = property(lambda self:self.att('state', '', typ=str ))
0154     def __repr__(self):
0155         return "%s %s %s %s" % (self.gidx, self.typ, self.name, self.state )
0156
0157
0158 class Transform(G):
0159     unit = property(lambda self:self.att('unit', "", typ=str ))
0160     x = property(lambda self:self.att('x', 0, typ=float))
0161     y = property(lambda self:self.att('y', 0, typ=float))
0162     z = property(lambda self:self.att('z', 0, typ=float))
0163     xyz = property(lambda self:[self.x, self.y, self.z] )
0164
0165     def __repr__(self):
0166         return "%s %s %s %s %s " % (self.typ, self.unit, self.x, self.y, self.z )
0167
0168 class Position(Transform):
0169     pass
0170 class Rotation(Transform):
0171     pass
0172 class Scale(Transform):
0173     pass
0174
0175
0176
0177 class Geometry(G):
0178     """
0179     Boolean and Primitive are Geometry subclasses
0180     """
0181     def as_ncsg(self):
0182         assert 0, "Geometry.as_ncsg needs to be overridden in the subclass: %s " % self.__class__
0183
0184     def as_shape(self, **kwa):
0185         assert 0, "Geometry.as_shape needs to be overridden in the subclass: %s " % self.__class__
0186
0187     def _get_subsolids(self):
0188         ss = []
0189         def subsolids_r(solid, top=False):
0190             if not top:
0191                 ss.append(solid.idx)
0192             pass
0193             if solid.is_boolean:
0194                 subsolids_r(solid.first)
0195                 subsolids_r(solid.second)
0196             pass
0197         pass
0198         subsolids_r(self, top=True)
0199         return ss
0200     subsolids = property(_get_subsolids)
0201     subsolidranges = property(lambda self:list(find_ranges(sorted(self.subsolids))))
0202     subsolidcount = property(lambda self:len(self.subsolids))
0203
0204     def find_solids(self, klsname):
0205         """
0206         GDML boolean aware find
0207
0208         :param klsname: eg Union, Ellipsoid
0209         :return ss: list of solids
0210         """
0211         ss = []
0212         def find_solids_r(solid):
0213             if solid.__class__.__name__ == klsname:
0214                 ss.append(solid)
0215             pass
0216             if solid.is_boolean:
0217                 find_solids_r(solid.first)
0218                 find_solids_r(solid.second)
0219             pass
0220         pass
0221         find_solids_r(self)
0222         return ss
0223
0224     def __repr__(self):
0225         """for indents need to sub_traverse from root solid of interest"""
0226         sub_depth = getattr(self, 'sub_depth', 0)
0227         sub_indent = "   " * sub_depth
0228         sub_name = getattr(self, 'sub_name', self.name)
0229         sub_prefix = getattr(self, 'sub_prefix', None)
0230
0231         sub_root = getattr(self, 'sub_root', None)
0232         is_sub_root = sub_root == self
0233         sub_root_marker = "[%s]" % sub_prefix if sub_root == self else ""
0234
0235         right_xyz = self.right_xyz if self.is_boolean else None
0236         right_xyz = " : right_xyz:%s/%s/%6.3f" % tuple(right_xyz) if not right_xyz is None else ""
0237
0238         line = "%d%s %s %s %s%s " % (sub_depth, sub_indent, self.gidx, self.typ, sub_root_marker, sub_name )
0239         line = "%-45s %s" % (line, right_xyz )
0240
0241         if is_sub_root:
0242             line = "*:%s" % line
0243         pass
0244
0245         lrep_ = lambda label,obj:"%s:%r"%(label,obj)
0246         lines = [line]
0247         if self.is_boolean:
0248             lines.append( lrep_("l",self.first) )
0249             lines.append( lrep_("r",self.second) )
0250         pass
0251         return "\n".join(lines)
0252
0253
0254
0255 class Boolean(Geometry):
0256     all_transforms = []
0257
0258     @classmethod
0259     def SaveBuffer(cls):
0260         num_boolean_transforms = len(cls.all_transforms)
0261         path = os.path.expandvars("$TMP/Boolean_all_transforms.npy")
0262         print("Boolean.all_transforms %d save to %s  " % ( num_boolean_transforms, path ))
0263         if num_boolean_transforms > 0:
0264              tbuf = np.vstack(cls.all_transforms).reshape(-1,4,4)
0265              np.save(path, tbuf)
0266         pass
0267     pass
0268
0269     firstref = property(lambda self:self.elem.find("first").attrib["ref"])
0270     secondref = property(lambda self:self.elem.find("second").attrib["ref"])
0271
0272     position = property(lambda self:self.find1_("position"))
0273     rotation = property(lambda self:self.find1_("rotation"))
0274     scale = None
0275     secondtransform = property(lambda self:construct_transform(self))
0276
0277     first = property(lambda self:self.g.solids[self.firstref])
0278     second = property(lambda self:self.g.solids[self.secondref])
0279
0280     @classmethod
0281     def SubTraverse(cls, root):
0282         """
0283         :param root: Geometry instance, Boolean or Primitive
0284
0285         Traverses from root setting "temporary" properties, sub_root, sub_depth, sub_name
0286
0287         As solids are reused within others cannot define a fixed
0288         depth or even a fixed parent link.
0289         Instead can only define temporary depths and parents within
0290         the context of a traversal from some root.
0291
0292         """
0293         names = []
0294         def collect_names_r(node, depth):
0295             names.append(node.name)
0296             if node.is_boolean:
0297                 collect_names_r(node.first, depth+1)
0298                 collect_names_r(node.second, depth+1)
0299             pass
0300         pass
0301         collect_names_r(root, 0)
0302         prefix = os.path.commonprefix(names)
0303         log.debug("\n".join(["names", "prefix:%s" % prefix] + names))
0304
0305         subnames = []
0306         def sub_traverse_r(node, depth):
0307
0308             sub_name = unref_(node.name[len(prefix):])
0309             subnames.append(sub_name)
0310
0311             node.sub_root = root
0312             node.sub_depth = depth
0313             node.sub_prefix = prefix
0314             node.sub_name = sub_name
0315
0316             if node.is_boolean:
0317                 node.first.sub_parent = node
0318                 node.second.sub_parent = node
0319                 sub_traverse_r(node.first, depth+1)
0320                 sub_traverse_r(node.second, depth+1)
0321             pass
0322         pass
0323         sub_traverse_r(root, 0)  # 2nd traverse doing sub_ labelling
0324         log.debug("\n".join(["subnames"] + subnames))
0325
0326     def sub_traverse(self):
0327         self.SubTraverse(self)
0328
0329     def as_ncsg(self):
0330         if not hasattr(self.first, 'as_ncsg'):
0331             print(self.first)
0332         pass
0333         left = self.first.as_ncsg()
0334         right = self.second.as_ncsg()
0335
0336         assert left, " left fail as_ncsg for first : %r self: %r " % (self.first, self)
0337         assert right, "right fail as_ncsg for second : %r self: %r " % (self.second, self)
0338
0339         transform = self.secondtransform
0340
0341         right.transform = transform
0342         self.__class__.all_transforms.append(transform)
0343
0344         cn = CSG(self.operation, name=self.name)
0345         cn.left = left
0346         cn.right = right
0347         return cn
0348
0349     def _get_right_xyz(self):
0350         transform = self.secondtransform
0351         no_rotation = np.all(np.eye(3) == transform[:3,:3])
0352         assert no_rotation, transform
0353         xyz = transform[3,:3]
0354         return xyz
0355     right_xyz = property(_get_right_xyz)
0356
0357     def as_shape(self, **kwa):
0358         """
0359         Only looking at single level of transform
0360         """
0361         left = self.first.as_shape(**kwa)
0362         right = self.second.as_shape(**kwa)
0363         transform = self.secondtransform
0364
0365         no_rotation = np.all(np.eye(3) == transform[:3,:3])
0366         assert no_rotation, transform
0367         xyz = transform[3,:3]
0368
0369         log.debug(xyz)
0370
0371         no_xy_shifts = xyz[0] == xyz[1] == 0
0372         if not no_xy_shifts:
0373             log.fatal("xy_shifts not handled %s " % repr(xyz))
0374         pass
0375         #assert no_xy_shifts, ("only z shifts handled ", xyz)
0376
0377
0378         assert left, " left fail as_shape for first : %r self: %r " % (self.first, self)
0379         assert right, "right fail as_shape for second : %r self: %r " % (self.second, self)
0380
0381         #tr = np.array( [0, xyz[2]] )
0382         tr = np.array( [xyz[0], xyz[2]] )
0383         shape = self.shape_kls(self.name, [left, right, tr], **kwa)
0384         return shape
0385
0386
0387 class Intersection(Boolean):
0388     operation = "intersection"
0389     shape_kls = SIntersectionSolid
0390
0391 class Subtraction(Boolean):
0392     operation = "difference"
0393     shape_kls = SSubtractionSolid
0394
0395 class Union(Boolean):
0396     operation = "union"
0397     shape_kls = SUnionSolid
0398
0399 class Primitive(Geometry):
0400     lunit = property(lambda self:self.att('lunit', 'mm', typ=str))
0401     aunit = property(lambda self:self.att('aunit', 'deg', typ=str))
0402     startphi = property(lambda self:self.att('startphi', 0, typ=float))
0403     deltaphi = property(lambda self:self.att('deltaphi',  360, typ=float))
0404     starttheta = property(lambda self:self.att('starttheta', 0, typ=float))
0405     deltatheta = property(lambda self:self.att('deltatheta', 180, typ=float))
0406     rmin = property(lambda self:self.att('rmin', 0, typ=float))
0407     rmax = property(lambda self:self.att('rmax', 0, typ=float))
0408
0409     x = property(lambda self:self.att('x', 0, typ=float))
0410     y = property(lambda self:self.att('y', 0, typ=float))
0411     z = property(lambda self:self.att('z', 0, typ=float))
0412
0413     @classmethod
0414     def fromstring(cls, st ):
0415         gg = GDML.fromstring("<gdml><solids>" + st + "</solids></gdml>")
0416         so = gg.solids(-1)   # pick last to allow composite booleans
0417         return so
0418
0419     @classmethod
0420     def deltaphi_segment_via_slab_DEPRECATED(cls, obj, phi0, phi1, dist):
0421         xyzw_ = lambda phi:(np.cos(phi*np.pi/180.), np.sin(phi*np.pi/180.),0,0)
0422
0423         slab_a = 0
0424         slab_b = dist
0425
0426         slab0 = CSG("slab", param=xyzw_(phi0+90),param1=[slab_a,slab_b,0,0] )
0427         slab1 = CSG("slab", param=xyzw_(phi1-90),param1=[slab_a,slab_b,0,0] )
0428
0429         # flipped signs to get segment from intended quadrant
0430
0431         obj_slab0 = CSG("intersection", left=obj, right=slab0 )
0432         obj_slab0_slab1 = CSG("intersection", left=obj_slab0, right=slab1 )
0433
0434         return obj_slab0_slab1
0435
0436     def __repr__0(self):
0437         depth = getattr(self, 'depth', 0)
0438         indent = " " * depth
0439         return "%d:%s %s %s %s %s rmin %s rmax %s  x %s y %s z %s  " % (depth, indent, self.gidx, self.typ, self.name, self.lunit, self.rmin, self.rmax, self.x, self.y, self.z)
0440
0441     def __repr__(self):
0442         grepr = Geometry.__repr__(self)
0443         return "%-45s : xyz %s,%s,%5.3f  " % (grepr, self.x, self.y, self.z)
0444
0445
0446 class Tube(Primitive):
0447     """
0448     G4Tubs is GDML serialized as tube
0449     """
0450     deltaphi_segment_enabled = True
0451
0452     def __repr__(self):
0453         hz = self.z/2.
0454         prepr = Primitive.__repr__(self)
0455         return "%-80s :  rmin %s rmax %6.3f hz %6.3f " % (prepr, self.rmin, self.rmax, hz )
0456
0457
0458     @classmethod
0459     def make_cylinder(cls, radius, z1, z2, name):
0460         cn = CSG("cylinder", name=name)
0461         cn.param[0] = 0
0462         cn.param[1] = 0
0463         cn.param[2] = 0
0464         cn.param[3] = radius
0465         cn.param1[0] = z1
0466         cn.param1[1] = z2
0467         return cn
0468
0469     @classmethod
0470     def make_disc(cls, x, y, inner, radius, z1, z2, name):
0471         cn = CSG("disc", name=name)
0472         cn.param[0] = x
0473         cn.param[1] = y
0474         cn.param[2] = inner
0475         cn.param[3] = radius
0476         cn.param1[0] = z1
0477         cn.param1[1] = z2
0478         return cn
0479
0480     def as_cylinder(self, nudge_inner=0.01):
0481         hz = self.z/2.
0482         has_inner = self.rmin > 0.
0483
0484         if has_inner:
0485             dz = hz*nudge_inner
0486             inner = self.make_cylinder(self.rmin, -(hz+dz), (hz+dz), self.name + "_inner")
0487         else:
0488             inner = None
0489         pass
0490         outer = self.make_cylinder(self.rmax, -hz, hz, self.name + "_outer" )
0491         tube = CSG("difference", left=outer, right=inner, name=self.name + "_difference" ) if has_inner else outer
0492
0493         has_deltaphi = self.deltaphi < 360
0494         if has_deltaphi and self.deltaphi_segment_enabled:
0495
0496              assert self.aunit == 'deg'
0497              phi0 = self.startphi
0498              phi1 = self.startphi + self.deltaphi
0499              sz  = self.z*1.01
0500              sr  = self.rmax*1.5
0501
0502              ## TODO: calulate how much the segmenting prism needs to poke beyind the radius
0503              ##       to avoid the outside plane from cutting the cylinder
0504
0505              segment = CSG.MakeSegment(phi0, phi1, sz, sr )
0506              log.info("as_cylinder doing phi0/phi1/sz/sr segmenting : name %s phi0 %s phi1 %s sz %s sr %s " % (self.name, phi0, phi1, sz, sr))
0507              tube_segment = CSG("intersection", left=tube, right=segment )
0508
0509              #tube_segment = self.deltaphi_slab_segment(tube, phi0, phi1, dist)
0510              #result.balance_disabled = True
0511
0512              result = tube_segment
0513         else:
0514              result = tube
0515         pass
0516         return result
0517
0518     def as_disc(self):
0519         hz = self.z/2.
0520         return self.make_disc( self.x, self.y, self.rmin, self.rmax, -hz, hz, self.name + "_disc" )
0521
0522     def as_ncsg(self, hz_disc_cylinder_cut=1.):
0523         hz = self.z/2.
0524         rmin = self.rmin
0525         rmax = self.rmax
0526         pick_disc = hz < hz_disc_cylinder_cut
0527         if pick_disc:
0528             log.debug("Tube.as_ncsg.CSG_DISC %s as hz < cut,  hz:%s cut:%s rmin:%s rmax:%s " % (self.name, hz, hz_disc_cylinder_cut, rmin, rmax))
0529         pass
0530         return self.as_disc() if pick_disc else self.as_cylinder()
0531
0532     def as_shape(self, **kwa):
0533         hz = self.z/2.
0534         shape = STubs(self.name, [self.rmax, hz], **kwa )
0535         return shape
0536
0537
0538 class Sphere(Primitive):
0539     deltaphi_slab_segment_enabled = False
0540
0541     def __repr__(self):
0542         prepr = Primitive.__repr__(self)
0543         return "%-80s :  rmin %s rmax %6.3f  " % (prepr, self.rmin, self.rmax )
0544
0545     def as_shape(self, **kwa):
0546         shape = SEllipsoid(self.name, [self.rmax, self.rmax], **kwa )
0547         return shape
0548
0549     def as_ncsg(self, only_inner=False):
0550         pass
0551         assert self.aunit == "deg" and self.lunit == "mm"
0552
0553         has_inner = not only_inner and self.rmin > 0.
0554         if has_inner:
0555             inner = self.as_ncsg(only_inner=True)  # recursive call to make inner
0556         pass
0557
0558         radius = self.rmin if only_inner else self.rmax
0559         assert radius is not None
0560
0561         startThetaAngle = self.starttheta
0562         deltaThetaAngle = self.deltatheta
0563
0564         x = 0
0565         y = 0
0566         z = 0
0567
0568         # z to the right, theta   0 -> z=r, theta 180 -> z=-r
0569         rTheta = startThetaAngle
0570         lTheta = startThetaAngle + deltaThetaAngle
0571
0572         assert rTheta >= 0. and rTheta <= 180.
0573         assert lTheta >= 0. and lTheta <= 180.
0574
0575         log.debug("Sphere.as_ncsg radius:%s only_inner:%s  has_inner:%s " % (radius, only_inner, has_inner))
0576
0577         zslice = startThetaAngle > 0. or deltaThetaAngle < 180.
0578
0579         if zslice:
0580             zmin = radius*math.cos(lTheta*math.pi/180.)
0581             zmax = radius*math.cos(rTheta*math.pi/180.)
0582             assert zmax > zmin, (startThetaAngle, deltaThetaAngle, rTheta, lTheta, zmin, zmax )
0583
0584             log.debug("Sphere.as_ncsg rTheta:%5.2f lTheta:%5.2f zmin:%5.2f zmax:%5.2f azmin:%5.2f azmax:%5.2f " % (rTheta, lTheta, zmin, zmax, z+zmin, z+zmax ))
0585
0586             cn = CSG("zsphere", name=self.name, param=[x,y,z,radius], param1=[zmin,zmax,0,0], param2=[0,0,0,0]  )
0587
0588             ZSPHERE_QCAP = 0x1 << 1   # zmax
0589             ZSPHERE_PCAP = 0x1 << 0   # zmin
0590             flags = ZSPHERE_QCAP | ZSPHERE_PCAP
0591
0592             cn.param2.view(np.uint32)[0] = flags
0593             pass
0594         else:
0595             cn = CSG("sphere", name=self.name)
0596             cn.param[0] = x
0597             cn.param[1] = y
0598             cn.param[2] = z
0599             cn.param[3] = radius
0600         pass
0601         if has_inner:
0602             ret = CSG("difference", left=cn, right=inner, name=self.name + "_difference"  )
0603         else:
0604             ret = cn
0605         pass
0606
0607         has_deltaphi = self.deltaphi < 360
0608         if has_deltaphi and not only_inner and self.deltaphi_slab_segment_enabled:
0609              assert self.aunit == 'deg'
0610              phi0 = self.startphi
0611              phi1 = self.startphi + self.deltaphi
0612              rmax = self.rmax + 1
0613              ret_segment = self.deltaphi_slab_segment(ret, phi0, phi1, rmax)
0614              result = ret_segment
0615         else:
0616              result = ret
0617         pass
0618
0619         return result
0620
0621
0622
0623
0624
0625 class Ellipsoid(Primitive):
0626     ax = property(lambda self:self.att('ax', 0, typ=float))
0627     by = property(lambda self:self.att('by', 0, typ=float))
0628     cz = property(lambda self:self.att('cz', 0, typ=float))
0629     zcut1 = property(lambda self:self.att('zcut1', 0, typ=float))
0630     zcut2 = property(lambda self:self.att('zcut2', 0, typ=float))
0631
0632     def _get_semi_axes(self):
0633         ax = self.ax
0634         by = self.by
0635         cz = self.cz
0636         a = np.array([ax,by,cz], dtype=np.float32)
0637         return a
0638     semi_axes = property(_get_semi_axes)
0639
0640     def as_ncsg(self):
0641         ax = self.semi_axes
0642         cn = CSG.MakeEllipsoid(axes=ax, name=self.name, zcut1=self.zcut1, zcut2=self.zcut2)
0643         return cn
0644
0645     def as_shape(self, **kwa):
0646         assert self.ax == self.by
0647         shape = SEllipsoid(self.name, [self.ax, self.cz], **kwa )
0648         return shape
0649
0650     def __repr__(self):
0651         prepr = Primitive.__repr__(self)
0652         return "%-80s :  ax/by/cz %6.3f/%6.3f/%6.3f  zcut1 %6.3f zcut2 %6.3f  " % (prepr, self.ax, self.by, self.cz, self.zcut1, self.zcut2 )
0653
0654
0655
0656
0657
0658 class Torus(Primitive):
0659     rtor = property(lambda self:self.att('rtor', 0, typ=float))
0660     def as_ncsg(self):
0661         cn = CSG.MakeTorus(R=self.rtor, r=self.rmax, name=self.name)
0662         return cn
0663
0664     def as_shape(self, **kwa):
0665         shape = STorus(self.name, [self.rmax, self.rtor], **kwa)
0666         return shape
0667
0668     def __repr__(self):
0669         prepr = Primitive.__repr__(self)
0670         return "%-80s :  rmin %6.3f rmax %6.3f rtor %6.3f  " % (prepr, self.rmin, self.rmax, self.rtor  )
0671
0672
0673
0674 class Box(Primitive):
0675     def as_ncsg(self):
0676         assert self.lunit == 'mm'
0677         cn = CSG("box3", name=self.name)
0678         cn.param[0] = self.x
0679         cn.param[1] = self.y
0680         cn.param[2] = self.z
0681         cn.param[3] = 0
0682         return cn
0683
0684 class Cone(Primitive):
0685     """
0686     GDML inner cone translated into CSG difference
0687
0688     rmin1: inner radius at base of cone
0689     rmax1: outer radius at base of cone
0690     rmin2: inner radius at top of cone
0691     rmax2: outer radius at top of cone
0692     z height of cone segment
0693
0694     """
0695     rmin1 = property(lambda self:self.att('rmin1', 0, typ=float))
0696     rmin2 = property(lambda self:self.att('rmin2', 0, typ=float))
0697
0698     rmax1 = property(lambda self:self.att('rmax1', 0, typ=float))
0699     rmax2 = property(lambda self:self.att('rmax2', 0, typ=float))
0700     z = property(lambda self:self.att('z', 0, typ=float))
0701     pass
0702
0703     @classmethod
0704     def make_cone(cls, r1,z1,r2,z2, name):
0705         cn = CSG("cone", name=name)
0706         cn.param[0] = r1
0707         cn.param[1] = z1
0708         cn.param[2] = r2
0709         cn.param[3] = z2
0710         return cn
0711
0712     def as_ncsg(self, only_inner=False):
0713         pass
0714         assert self.aunit == "deg" and self.lunit == "mm" and self.deltaphi == 360. and self.startphi == 0.
0715         has_inner = not only_inner and (self.rmin1 > 0. or self.rmin2 > 0. )
0716         if has_inner:
0717             inner = self.as_ncsg(only_inner=True)  # recursive call to make inner
0718         pass
0719
0720         r1 = self.rmin1 if only_inner else self.rmax1
0721         z1 = -self.z/2
0722
0723         r2 = self.rmin2 if only_inner else self.rmax2
0724         z2 = self.z/2
0725
0726         cn = self.make_cone( r1,z1,r2,z2, self.name )
0727
0728         return CSG("difference", left=cn, right=inner ) if has_inner else cn
0729
0730     def __repr__(self):
0731         return "%s %s z:%s rmin1:%s rmin2:%s rmax1:%s rmax2:%s  " % (self.gidx, self.typ, self.z, self.rmin1, self.rmin2, self.rmax1, self.rmax2 )
0732
0733     def plot(self, ax):
0734         zp0 = FakeZPlane(z=0., rmin=self.rmin1,rmax=self.rmax1)
0735         zp1 = FakeZPlane(z=self.z, rmin=self.rmin2,rmax=self.rmax2)
0736         zp = [zp0, zp1]
0737         Polycone.Plot(ax, zp)
0738
0739
0740
0741 class FakeZPlane(object):
0742     def __init__(self, z, rmin, rmax):
0743         self.z = z
0744         self.rmin = rmin
0745         self.rmax = rmax
0746
0747     def __repr__(self):
0748         return "%s %s %s %s " % (self.__class__.__name__, self.z, self.rmin, self.rmax)
0749
0750
0751 class ZPlane(Primitive):
0752     pass
0753
0754
0755 class Trapezoid(Primitive):
0756     """
0757     The GDML Trapezoid is formed using 5 dimensions:
0758
0759     x1: x length at -z
0760     x2: x length at +z
0761     y1: y length at -z
0762     y2: y length at +z
0763     z:  z length
0764
0765     The general form is ConvexPolyhedron, modelled by a list of planes
0766     defining half spaces which come together to define the polyhedron.
0767
0768     * 4:tetrahedron
0769     * 5:triangular prism
0770     * 6:trapezoid, cube, box
0771
0772     Q&A:
0773
0774     * where to transition from specific Trapezoid to generic ConvexPolyhedron ? Did this immediately with python input.
0775     * where to extract the bbox, thats easier before conversion to generic ? Again immediately in python.
0776
0777     * :google:`convex polyhedron bounding box`
0778
0779     Enumerating the vertices of a convex polyhedron defined by its planes
0780     looks complicated, hence need to extract and store the bbox prior
0781     to generalization of the shape.
0782
0783     * http://cgm.cs.mcgill.ca/~avis/doc/avis/AF92b.pdf
0784
0785     ::
0786
0787         In [1]: run gdml.py
0788
0789         In [2]: trs = g.findall_("solids//trd")
0790
0791         In [3]: len(trs)
0792         Out[3]: 2
0793
0794         In [13]: trs = g.findall_("solids//trd")
0795
0796         In [14]: trs
0797         Out[14]:
0798         [Trapezoid name:SstTopRadiusRibBase0xc271078 z:2228.5 x1:160.0 y1:20.0 x2:691.02 y2:20.0  ,
0799          Trapezoid name:SstInnVerRibCut0xbf31118 z:50.02 x1:100.0 y1:27.0 x2:237.2 y2:27.0  ]
0800
0801         In [4]: print trs[0]
0802         <trd xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" lunit="mm" name="SstTopRadiusRibBase0xc271078" x1="160" x2="691.02" y1="20" y2="20" z="2228.5"/>
0803
0804
0805         In [5]: print trs[1]
0806         <trd xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" lunit="mm" name="SstInnVerRibCut0xbf31118" x1="100" x2="237.2" y1="27" y2="27" z="50.02"/>
0807
0808
0809     Tips for finding LV with a solid... use filternodes_so to get some nodes
0810     then can use parent links to get idea of local tree around the solid.
0811
0812     ::
0813
0814         In [17]: sivr = t.filternodes_so("SstInnVerRib")
0815
0816         In [18]: sivr[0]
0817         Out[18]:
0818         Node 4473 : dig a60e pig 8be9 depth 11 nchild 0
0819         pv:PhysVol /dd/Geometry/AD/lvOIL#pvSstInnVerRibs#SstInnVerRibs#SstInnVerRibRot0xbf1abc0
0820          Position mm 2428.0 0.0 -40.0  None
0821         lv:Volume /dd/Geometry/AdDetails/lvSstInnVerRibBase0xbf31748 /dd/Materials/StainlessSteel0xc2adc00 SstInnVerRibBase0xbf30b50
0822            Subtraction SstInnVerRibBase0xbf30b50
0823              l:Box SstInnVerRibBox0xbf310d8 mm rmin 0.0 rmax 0.0  x 120.0 y 25.0 z 4875.0
0824              r:Trapezoid name:SstInnVerRibCut0xbf31118 z:50.02 x1:100.0 y1:27.0 x2:237.2 y2:27.0
0825            Material /dd/Materials/StainlessSteel0xc2adc00 solid : Position mm 2428.0 0.0 -40.0
0826
0827         In [26]: print list(map(lambda n:n.index, sivr))
0828         [4473, 4474, 4475, 4476, 4477, 4478, 4479, 4480, 6133, 6134, 6135, 6136, 6137, 6138, 6139, 6140]
0829
0830
0831         In [28]: sivr[0].parent
0832         Out[28]:
0833         Node 3155 : dig 8be9 pig a856 depth 10 nchild 520
0834         pv:PhysVol /dd/Geometry/AD/lvSST#pvOIL0xc241510
0835          Position mm 0.0 0.0 7.5  Rotation deg 0.0 0.0 -180.0
0836         lv:Volume /dd/Geometry/AD/lvOIL0xbf5e0b8 /dd/Materials/MineralOil0xbf5c830 oil0xbf5ed48
0837            Tube oil0xbf5ed48 mm rmin 0.0 rmax 2488.0  x 0.0 y 0.0 z 4955.0
0838            Material /dd/Materials/MineralOil0xbf5c830 solid
0839            PhysVol /dd/Geometry/AD/lvOIL#pvOAV0xbf8f638
0840
0841
0842
0843     """
0844     x1 = property(lambda self:self.att('x1', 0, typ=float))
0845     x2 = property(lambda self:self.att('x2', 0, typ=float))
0846     y1 = property(lambda self:self.att('y1', 0, typ=float))
0847     y2 = property(lambda self:self.att('y2', 0, typ=float))
0848     z = property(lambda self:self.att('z', 0, typ=float))
0849
0850     def __repr__(self):
0851         return "%s name:%s z:%s x1:%s y1:%s x2:%s y2:%s  " % (self.typ, self.name, self.z, self.x1, self.y1, self.x2, self.y2 )
0852
0853     def as_ncsg(self):
0854         assert self.lunit == 'mm'
0855         cn = CSG("trapezoid", name=self.name)
0856         planes, verts, bbox, src = make_trapezoid(z=self.z, x1=self.x1, y1=self.y1, x2=self.x2, y2=self.y2 )
0857         cn.planes = planes
0858         cn.param2[:3] = bbox[0]
0859         cn.param3[:3] = bbox[1]
0860
0861         cn.meta.update(src.srcmeta)
0862
0863         return cn
0864
0865
0866 class Polycone(Primitive):
0867     """
0868     ::
0869
0870         In [1]: run gdml.py
0871
0872         In [2]: pcs = gdml.findall_("solids//polycone")
0873
0874         In [3]: len(pcs)
0875         Out[3]: 13
0876
0877         In [6]: print pcs[0]
0878         <polycone xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" aunit="deg" deltaphi="360" lunit="mm" name="gds_polycone0xc404f40" startphi="0">
0879           <zplane rmax="1520" rmin="0" z="3070"/>              # really flat cone
0880           <zplane rmax="75" rmin="0" z="3145.72924106399"/>    # little cylinder
0881           <zplane rmax="75" rmin="0" z="3159.43963177189"/>
0882         </polycone>
0883
0884        In [7]: print pcs[1]
0885         <polycone xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" aunit="deg" deltaphi="360" lunit="mm" name="iav_polycone0xc346448" startphi="0">
0886           <zplane rmax="1565" rmin="0" z="3085"/>                 # flat cylinder (repeated rmax)
0887           <zplane rmax="1565" rmin="0" z="3100"/>                 # repeated z, changed r   <--- polycone abuse, zero height truncated cone
0888           <zplane rmax="1520.39278882354" rmin="0" z="3100"/>     #    big taper cone
0889           <zplane rmax="100" rmin="0" z="3174.43963177189"/>
0890         </polycone>
0891
0892
0893     Summarizing with zp, see that
0894
0895     * all using same rmin, thats polycone abuse : better to subtract a cylinder
0896     * common pattern repeated z with different rmax to model "lips"
0897
0898     ::
0899
0900         pcs[ 0].zp :          gds_polycone0xc404f40  3 z:      [3145.72924106399, 3159.43963177189, 3070.0] rmax:                     [1520.0, 75.0] rmin:               [0.0]
0901         pcs[ 7].zp :          SstTopHubBot0xc2635b8  2 z:                                  [-320.0, -340.0] rmax:                            [220.5] rmin:             [150.5]
0902         pcs[ 8].zp :         SstTopHubMain0xc263d80  2 z:                                     [-320.0, 0.0] rmax:                            [170.5] rmin:             [150.5]
0903
0904
0905         pcs[ 1].zp :          iav_polycone0xc346448  4 z:                [3100.0, 3085.0, 3174.43963177189] rmax:  [100.0, 1520.39278882354, 1565.0] rmin:               [0.0]
0906         pcs[ 2].zp :             IavTopHub0xc405968  4 z:         [0.0, 85.5603682281126, 110.560368228113] rmax:                     [100.0, 150.0] rmin:              [75.0]
0907         pcs[ 3].zp :       CtrGdsOflBotClp0xbf5dec0  4 z:                                 [0.0, 25.0, 30.0] rmax:                      [36.5, 150.0] rmin:              [31.5]
0908         pcs[ 4].zp :          OcrGdsPrtPln0xbfa1408  4 z:                               [0.0, 160.0, 185.0] rmax:                     [100.0, 150.0] rmin:              [75.0]
0909         pcs[ 5].zp :          lso_polycone0xc02a418  4 z:      [4076.62074383385, 4058.59604160589, 3964.0] rmax:              [1930.0, 50.0, 125.0] rmin:               [0.0]
0910         pcs[ 6].zp :          oav_polycone0xbf1c840  4 z:      [4094.62074383385, 3937.0, 4000.02470222796] rmax:            [2040.0, 1930.0, 125.0] rmin:               [0.0]
0911
0912         pcs[ 9].zp :             OavTopHub0xc2c9030  6 z:                          [0.0, 57.0, 90.0, 120.0] rmax:                [98.0, 68.0, 125.0] rmin:              [50.0]
0913         pcs[10].zp :       CtrLsoOflTopClp0xc178498  6 z:                         [0.0, 16.0, 184.0, 200.0] rmax:              [102.5, 112.5, 100.0] rmin:              [50.0]
0914         pcs[11].zp :       OcrGdsLsoPrtPln0xc104000  4 z:         [0.0, 214.596041605889, 184.596041605889] rmax:                       [98.0, 68.0] rmin:              [50.0]
0915         pcs[12].zp :       OcrCalLsoPrtPln0xc2fadd8  4 z:         [0.0, 214.596041605889, 184.596041605889] rmax:                       [98.0, 68.0] rmin:              [50.0]
0916
0917
0918     """
0919     pass
0920     zplane = property(lambda self:self.findall_("zplane"))
0921
0922     zp_num = property(lambda self:len(self.zplane))
0923     zp_rmin = property(lambda self:list(set(map(lambda _:_.rmin,self.zplane))))
0924     zp_rmax = property(lambda self:list(set(map(lambda _:_.rmax,self.zplane))))
0925     zp_z = property(lambda self:list(set(map(lambda _:_.z,self.zplane))))
0926     zp = property(lambda self:"%s %30s %2d z:%50r rmax:%35r rmin:%20r " % (self.gidx, self.name, self.zp_num, self.zp_z, self.zp_rmax, self.zp_rmin))
0927
0928     #def __repr__(self):
0929     #    return self.zp
0930
0931     def __repr__(self):
0932         prepr = Primitive.__repr__(self)
0933         return "%-80s :  zp_num %2d z:%r rmax:%r rmin:%r  " % (prepr, self.zp_num, self.zp_z, self.zp_rmax, self.zp_rmin )
0934
0935
0936     def zp_array(self):
0937         a = np.zeros( [len(self.zplane), 3] )
0938         for i in range(len(self.zplane)):
0939             zp = self.zplane[i]
0940             a[i] = [ zp.rmin, zp.rmax, zp.z ]
0941         pass
0942         return a
0943
0944     def prims(self):
0945         """
0946         Auto correct simple case of wrong z-order::
0947
0948                421     <polycone aunit="deg" deltaphi="360" lunit="mm" name="PMT_3inch_pmt_solid_cyl0x1c9da50" startphi="0">
0949                422       <zplane rmax="30.001" rmin="0" z="-15.8745078663875"/>
0950                423       <zplane rmax="30.001" rmin="0" z="-75.8755078663876"/>
0951                424     </polycone>
0952
0953         """
0954         zp = self.zplane
0955         zz = list(map(lambda _:_.z, zp))
0956
0957         if len(zp) == 2 and zz[0] > zz[1]:
0958            log.warning("Polycone swap misordered pair of zplanes for %s " % self.name)
0959            zp = list(reversed(zp))
0960         pass
0961
0962         prims = []
0963         for i in range(1,len(zp)):
0964             zp1 = zp[i-1]
0965             zp2 = zp[i]
0966
0967             r1 = zp1.rmax
0968             r2 = zp2.rmax
0969             z1 = zp1.z
0970             z2 = zp2.z
0971
0972             if z2 == z1:
0973                 log.debug("skipping z2 == z1 zp" )
0974             else:
0975                 #assert z2 > z1, (z2,z1)
0976                 if not z2 > z1:
0977                     raise ValueError("Polycone bad z-order expect z2>z1 : but z1 %s z2 %s " % (z1,z2 ))
0978                 pass
0979                 name = self.name + "_zp_%d" % i
0980                 pr = Tube.make_cylinder( r1, z1, z2, name ) if r1 == r2 else Cone.make_cone( r1, z1, r2, z2, name )
0981                 prims.append(pr)
0982             pass
0983         pass
0984         return prims
0985
0986     def inner(self):
0987         rmin = self.zp_rmin
0988         assert len(rmin) == 1
0989         rmin = rmin[0]
0990
0991         zp = self.zplane
0992         z = list(map(lambda _:_.z, zp))
0993         zmax = max(z)
0994         zmin = min(z)
0995
0996         has_inner = rmin > 0.
0997         if has_inner:
0998             inner = Tube.make_cylinder( rmin, zmin, zmax, self.name + "_inner_cylinder" )
0999         else:
1000             inner = None
1001         pass
1002         return inner
1003
1004
1005     def as_ncsg(self):
1006         assert self.aunit == "deg" and self.lunit == "mm" and self.deltaphi == 360. and self.startphi == 0.
1007         try:
1008             prims = self.prims()
1009         except ValueError as e:
1010             log.fatal("Polycone.as_ncsg failed ValueError : %r " % e )
1011             return None
1012         pass
1013         cn = TreeBuilder.uniontree(prims, name=self.name + "_uniontree")
1014         inner = self.inner()
1015         #return CSG("difference", left=cn, right=inner ) if inner is not None else cn
1016         return cn
1017
1018     def as_shape(self, **kwa):
1019         param = self.zp_array()
1020         shape = SPolycone(self.name, param, **kwa)
1021         return shape
1022
1023     def plot(self, ax):
1024         self.Plot(ax, self.zplane)
1025
1026     @classmethod
1027     def Plot(cls, ax, zp):
1028
1029         rmin = list(map(lambda _:_.rmin, zp))
1030         rmax = list(map(lambda _:_.rmax, zp))
1031
1032         z = list(map(lambda _:_.z, zp))
1033         zmax = max(z)
1034         zmin = min(z)
1035         zsiz = zmax - zmin
1036
1037         for i in range(1,len(zp)):
1038             zp0 = zp[i-1]
1039             zp1 = zp[i]
1040             ax.plot( [zp0.rmax,zp1.rmax], [zp0.z,zp1.z])
1041             ax.plot( [-zp0.rmax,-zp1.rmax], [zp0.z,zp1.z])
1042
1043             if i == 1:
1044                 ax.plot( [-zp0.rmax, zp0.rmax], [zp0.z, zp0.z] )
1045             pass
1046             if i == len(zp) - 1:
1047                 ax.plot( [-zp1.rmax, zp1.rmax], [zp1.z, zp1.z] )
1048             pass
1049         pass
1050
1051         for i in range(1,len(zp)):
1052             zp0 = zp[i-1]
1053             zp1 = zp[i]
1054             ax.plot( [zp0.rmin,zp1.rmin], [zp0.z,zp1.z])
1055             ax.plot( [-zp0.rmin,-zp1.rmin], [zp0.z,zp1.z])
1056
1057             if i == 1:
1058                 ax.plot( [-zp0.rmin, zp0.rmin], [zp0.z, zp0.z] )
1059             pass
1060             if i == len(zp) - 1:
1061                 ax.plot( [-zp1.rmin, zp1.rmin], [zp1.z, zp1.z] )
1062             pass
1063         pass
1064
1065         xmax = max(rmax)*1.2
1066         ax.set_xlim(-xmax, xmax )
1067         ax.set_ylim( zmin - 0.1*zsiz, zmax + 0.1*zsiz )
1068
1069
1070
1071
1072
1073
1074 class Volume(G):
1075     """
1076     Volume : GDML volume are Geant4 Logical Volumes (LV)
1077     ======================================================
1078
1079     As very many pv (placements of lv) share the same lv
1080     it does not make sense to simply jump from an lv to its pv.
1081     Conversely that is a question that must be asked of the whole
1082     geometry tree to select all the placements (pv) that have a particular lv.
1083
1084     ::
1085
1086         In [15]: for v in gdml.volumes.values():print v.material.shortname
1087         PPE
1088         MixGas
1089         Air
1090         Bakelite
1091         Air
1092         Bakelite
1093         Foam
1094         Aluminium
1095         Air
1096         ...
1097
1098     """
1099     materialref = property(lambda self:self.elem.find("materialref").attrib["ref"])
1100     solidref = property(lambda self:self.elem.find("solidref").attrib["ref"])
1101     solid = property(lambda self:self.g.solids[self.solidref])
1102     material = property(lambda self:self.g.materials[self.materialref])
1103
1104     physvol = property(lambda self:self.findall_("physvol"))
1105     children = property(lambda self:self.physvol)
1106
1107     def physvol_xyz(self, pfx="pLPMT_Hamamatsu", sub="position"):
1108         """
1109         :param pfx: pv name prefix
1110         :param sub: element name, must be either "position" or "rotation"
1111         :return xyz: array of shape (npv, 3)
1112         """
1113         assert sub in ["position", "rotation"]
1114         x = list(map(float,self.elem.xpath("./physvol[starts-with(@name,'%s')]/%s/@x" % (pfx,sub))))
1115         y = list(map(float,self.elem.xpath("./physvol[starts-with(@name,'%s')]/%s/@y" % (pfx,sub))))
1116         z = list(map(float,self.elem.xpath("./physvol[starts-with(@name,'%s')]/%s/@z" % (pfx,sub))))
1117         assert len(x) == len(y) == len(z)
1118         xyz = np.zeros( [len(x), 3])
1119         xyz[:,0] = x
1120         xyz[:,1] = y
1121         xyz[:,2] = z
1122         return xyz
1123
1124
1125
1126     def filterpv(self, pfx):
1127         return list(filter(lambda pv:pv.name.startswith(pfx), self.physvol))
1128
1129     def rdump(self, depth=0):
1130         print(self)
1131         for pv in self.physvol:
1132             lv = pv.volume
1133             lv.rdump(depth=depth+1)
1134         pass
1135
1136     def as_ncsg(self):
1137         """
1138         Hmm pv level transforms need to be applied to the
1139         csg top nodes of the solids, as the pv transforms
1140         are not being propagated GPU side, and even when they
1141         are will want to be able to handle a csg instance
1142         comprising a few solids only (eg the PMT)
1143         as is with its own transforms.
1144
1145         Essentially are collapsing a subtree for the
1146         handful of solids into the self contained instance
1147         of a list of csgnodes.
1148         """
1149         pass
1150
1151
1152     def __repr__(self):
1153         repr_ = lambda _:"   %r" % _
1154         pvs = list(map(repr_, self.physvol))
1155         line = "%s %s %s" % (self.gidx, self.typ, self.name )
1156         return "\n".join([line, "solid", repr(self.solid), "material", repr(self.material), "physvol %d" % len(self.physvol)] + pvs )
1157
1158 class Physvol(G):
1159     """
1160
1161     volume
1162          opticks.ana.pmt.gdml.Volume
1163
1164     transform
1165          4x4 np.ndarray homogeneous TRS matrix ie Translate-Rotate-Scale
1166          combining the position, rotation and scale attributes
1167
1168     position
1169     rotation
1170     scale
1171          dont use these, use transform that combines them
1172
1173
1174     """
1175     volumeref = property(lambda self:self.elem.find("volumeref").attrib["ref"])
1176     volume = property(lambda self:self.g.volumes[self.volumeref])
1177     children = property(lambda self:[self.volume])
1178
1179     position = property(lambda self:self.find1_("position"))
1180     rotation = property(lambda self:self.find1_("rotation"))
1181     scale = property(lambda self:self.find1_("scale"))
1182     transform = property(lambda self:construct_transform(self))
1183
1184     def __repr__(self):
1185         return "\n".join(["%s %s" % (self.typ, self.name)," %r %r " % ( self.position, self.rotation)])
1186
1187
1188
1189 class odict(collections.OrderedDict):
1190     """
1191     Used for GDML collections of materials, solids and volumes which are
1192     always keyed by name.
1193
1194     Call method gives integer index access to materials, solids and volumes (LV)
1195     The idx is set in GDML.init from the order the elements appear in the GDML.
1196
1197     ::
1198
1199         g.material(idx)
1200         g.solids(idx)
1201         g.volumes(idx)   # integer index
1202
1203     """
1204     def __call__(self, iarg):
1205         return list(self.items())[iarg][1]   # items returns (key,value) the key is the name
1206
1207
1208 class GDML(G):
1209     """
1210     Parsing and wrapping
1211     """
1212     kls = {
1213         "material":Material,
1214
1215         "tube":Tube,
1216         "torus":Torus,
1217         "sphere":Sphere,
1218         "ellipsoid":Ellipsoid,
1219         "box":Box,
1220         "cone":Cone,
1221         "polycone":Polycone,
1222         "zplane":ZPlane,
1223         "trd":Trapezoid,
1224
1225         "intersection":Intersection,
1226         "subtraction":Subtraction,
1227         "union":Union,
1228
1229         "position":Position,
1230         "rotation":Rotation,
1231         "scale":Scale,
1232
1233         "volume":Volume,
1234         "physvol":Physvol,
1235     }
1236
1237     @classmethod
1238     def parse(cls, path):
1239         """
1240         :param path: to GDML file
1241         """
1242         #if path is None:
1243         #    path = os.environ['OPTICKS_GDMLPATH']   # set within opticks_main
1244         #    log.info("gdmlpath defaulting to OPTICKS_GDMLPATH %s which is derived by opticks.ana.base from the IDPATH input envvar " % path )
1245         pass
1246         log.info("parsing gdmlpath %s " % (path) )
1247         gdml_elem = parse_(path)         # lxml parse and return root "gdml" element
1248         wgg = cls.wrap(gdml_elem, path=path)
1249         return wgg
1250
1251     @classmethod
1252     def fromstring(cls, st ):
1253         log.info("parsing string %s " % st )
1254         gg = ET.fromstring(st)
1255         wgg = cls.wrap(gg)
1256         return wgg
1257
1258     @classmethod
1259     def wrap(cls, gdml, path=None):
1260         """
1261         :param gdml: lxml root "gdml" element of parsed GDML document
1262         :param path: of the gdml as metadata
1263         """
1264         log.debug("wrapping gdml elements with python instances")
1265
1266         if not path is None:
1267             label = os.path.splitext(os.path.basename(path))[0]
1268         else:
1269             label = None
1270         pass
1271
1272         gg = cls(gdml)
1273         gg.g = gg
1274         gg.path = path
1275         gg.path_label = label
1276         gg.string = tostring_(gdml)
1277         gg.init()
1278         return gg
1279
1280     def get_traversed_volumes(self, lv_base, maxdepth=-1):
1281         """
1282         :param lv_base: base logical volume
1283         :param maxdepth: -1 for no limit
1284         :return tvol: odict of logical volumes within lv_base, obtained by recursion and keyed by index
1285
1286         Recursive traversal of the lv/pv GDML structure collecting
1287         and ordered dict keyed on the traversal order index
1288
1289         Note this was originally intended for assigning shortened local names
1290         within subtrees of associated volumes, eg those forming the PMT-mask and its PMT inside,
1291         that have associated names with a common prefix.
1292         """
1293         tvol = odict()
1294         def traverse_r(lv0, depth):
1295
1296             pvs = lv0.physvol
1297             indent = "   " * depth
1298             print("[%2d] %s %4d %s " % (depth, indent,len(pvs),lv0.name))
1299
1300             local_index = len(tvol)
1301             tvol[local_index] = lv0
1302             lv0.local_index = local_index
1303             lv0.local_depth = depth
1304
1305             if depth < maxdepth or maxdepth == -1:
1306                 for pv in pvs:
1307                     lv = pv.volume
1308                     traverse_r(lv, depth+1)
1309                 pass
1310             pass
1311         pass
1312         traverse_r(lv_base, 0)
1313         self.label_traversed_volumes(tvol)
1314         return tvol
1315
1316     def label_traversed_volumes(self, tvol ):
1317         """
1318         Assign shortened local_name to the lv
1319         by removing common prefix and pointer refs
1320
1321         NB this assumes that associated volumes have associated names
1322         with a common prefix
1323
1324         """
1325         names = list(map(lambda _:_.name, tvol.values()))
1326         unref_ = lambda _:_[:_.find("0x")] if _.find("0x") > -1 else _
1327
1328         prefix = os.path.commonprefix(names)
1329         for lv in tvol.values():
1330             lv.local_name = unref_( lv.name[len(prefix):] )
1331             lv.local_title = "%d:%s/%s" % ( lv.local_index, lv.local_name, unref_(lv.material.name) )
1332             lv.local_prefix = "%s/%s" % (lv.g.path_label, prefix)
1333         pass
1334
1335     def find_by_prefix(self, d, prefix):
1336         """
1337         :param d: dict
1338         :param prefix: str
1339         :return sel: list of objects from d with name starting with prefix
1340         """
1341         return list(filter(lambda v:v.name.startswith(prefix), d.values()))
1342
1343     def find_volumes(self, prefix):
1344         return self.find_by_prefix(self.volumes, prefix)
1345
1346     def find_solids(self, prefix):
1347         return self.find_by_prefix(self.solids, prefix)
1348
1349     def find_materials(self, prefix):
1350         return self.find_by_prefix(self.materials, prefix)
1351
1352     def smry(self):
1353         g = self
1354         ns = len(g.solids.keys())
1355         nv = len(g.volumes.keys())
1356         log.info(" g.path %s g.path_label %s ns:%d nv:%d (logical volumes) " % (g.path, g.path_label,ns,nv))
1357
1358     def find_one_volume(self, prefix):
1359         """
1360         :param prefix: str
1361         :return lv: first Volume instance with name matching starting with prefix
1362         """
1363         lvs = self.find_volumes(prefix)
1364         log.info("prefix argument %s matched %d volumes" % (prefix, len(lvs)))
1365         for i,lv in enumerate(lvs):
1366             log.info(" %2d : %s" % (i,lv.name) )
1367         pass
1368         return lvs[0] if len(lvs) > 0 else None
1369
1370
1371
1372
1373
1374
1375     def volume_summary(self):
1376         """
1377         1. get unique list of all volumeref/@ref lv names referenced from physvol
1378         2. for each lv name count the number of physvol that references it
1379         3. for each physvol get all parent lv names
1380
1381         25600 :                             PMT_3inch_log0x3a2d940 : lInnerWater0x30e6330
1382         12612 :                    NNVTMCPPMTlMaskVirtual0x32a49b0 : lInnerWater0x30e6330
1383          5000 :               HamamatsuR12860lMaskVirtual0x32906d0 : lInnerWater0x30e6330
1384          2400 :          mask_PMT_20inch_vetolMaskVirtual0x32a6960 : lOuterWaterPool0x30e5560
1385           590 :                                    lUpper0x31765d0 : lInnerWater0x30e6330
1386           590 :                                 lAddition0x31bd440 : lInnerWater0x30e6330
1387           590 :                                lFasteners0x312ddf0 : lInnerWater0x30e6330
1388           590 :                                    lSteel0x30e9b10 : lInnerWater0x30e6330
1389            64 :                                  lCoating0x45073a0 : lPanelTape0x4507210
1390            64 :                                lXJfixture0x320b950 : lTarget0x30e7010 lInnerWater0x30e6330
1391            63 :                                  lWallff_0x4506d50 : lAirTT0x4506ad0
1392            56 :                                 lXJanchor0x3205390 : lInnerWater0x30e6330
1393            36 :                                lSJFixture0x3210610 : lTarget0x30e7010
1394             8 :                               lSJReceiver0x3215340 : lTarget0x30e7010
1395             4 :                                    lPanel0x4507080 : lPlanef_0x4506f70
1396             2 :                                  lPlanef_0x4506f70 : lWallff_0x4506d50
1397             2 :                              lSJCLSanchor0x320f400 : lTarget0x30e7010
1398             1 :                             lLowerChimney0x4504230 : lInnerWater0x30e6330
1399             1 :                        lLowerChimneySteel0x4504880 : lLowerChimney0x4504230
1400             1 :                                   lTarget0x30e7010 : lAcrylic0x30e69a0
1401             1 :                             lUpperChimney0x4501f50 : lExpHall0x30dfe10
1402             1 :                           lLowerChimneyLS0x4504660 : lLowerChimney0x4504230
1403             1 :                                      lBar0x4507530 : lCoating0x45073a0
1404             1 :                      lLowerChimneyAcrylic0x4504450 : lLowerChimney0x4504230
1405             1 :                                  lBtmRock0x30e4820 : lWorld0x30d4f70
1406             1 :     HamamatsuR12860_PMT_20inch_inner2_log0x32b85a0 : HamamatsuR12860_PMT_20inch_body_log0x32b8250
1407             1 :                      PMT_3inch_inner2_log0x3a2db60 : PMT_3inch_body_log0x3a2d830
1408             1 :                               lInnerWater0x30e6330 : lReflectorInCD0x30e5cc0
1409
1410
1411         """
1412         lvns = list(set(self.elem.xpath("//physvol/volumeref/@ref")))
1413
1414         fmt = "%5d : %50s : %s"
1415         d = {}
1416         p = {}
1417         for lvn in lvns:
1418             d[lvn] = len(self.elem.xpath("//physvol[volumeref/@ref = '%s']" % lvn))
1419             p[lvn] = self.elem.xpath("//physvol[volumeref/@ref = '%s']/../@name" % lvn)
1420             log.debug(fmt  % (d[lvn], lvn, p[lvn] ))
1421         pass
1422         print("\n".join([fmt % (d[lvn], lvn, " ".join(p[lvn]) ) for lvn in sorted(d, reverse=True, key=lambda lvn:d[lvn])]))
1423
1424
1425
1426
1427     world = property(lambda self:self.volumes[self.worldvol])
1428
1429     def init(self):
1430         """
1431         Heart of the GDML parsing.
1432         """
1433         self.materials = odict()
1434         self.solids = odict()
1435         self.volumes = odict()
1436
1437         for i, e in enumerate(self.findall_("materials/material")):
1438             e.idx = i
1439             self.materials[e.name] = e
1440         pass
1441         # opticalsurface also in solids so switch from elem.findall to elem.xpath and exclude opticalsurface
1442         for i, e in enumerate(self.findall_("solids/*[local-name() != 'opticalsurface']")):
1443             e.idx = i
1444             self.solids[e.name] = e
1445         pass
1446         # skinsurface,bordersurface also in structure so : structure/* -> structure/volume
1447         for i, e in enumerate(self.findall_("structure/volume")):
1448             e.idx = i
1449             self.volumes[e.name] = e
1450         pass
1451         self.setup_world = self.elem.find("setup/world")
1452         self.worldvol = self.setup_world.attrib["ref"] if self.setup_world is not None else None
1453
1454         vv = self.volumes.values()
1455
1456         vvs = list(filter(lambda v:hasattr(v,'solid'), vv))
1457
1458         log.info("vv %s (number of logical volumes) vvs %s (number of lv with associated solid) " % (len(vv),len(vvs)))
1459         #for v in vvs:print repr(v)
1460         # getting the much larger number of physvol (ie placements of lv)
1461         # entails traversing the full geometry tree
1462
1463         self.lv2so = dict([(v.idx, v.solid.idx) for v in vvs])
1464
1465
1466
1467
1468
1469
1470
1471
1472 def oneplot(obj):
1473     fig = plt.figure()
1474     ax = fig.add_subplot(1,1,1, aspect='equal')
1475     obj.plot(ax)
1476     plt.show()
1477
1478 def multiplot(fig, objs, nx=4, ny=4):
1479     for i in range(len(objs)):
1480         ax = fig.add_subplot(nx,ny,i+1, aspect='equal')
1481         objs[i].plot(ax)
1482     pass
1483
1484
1485
1486
1487 def test_children(t):
1488
1489     lvn = "/dd/Geometry/PMT/lvPmtHemi0x"
1490
1491     l = t.filternodes_lv(lvn)  # all nodes
1492     assert len(l) == 672
1493
1494     n = l[0]
1495     sibs = list(n.siblings)
1496     assert len(sibs) == 192
1497
1498     cn = n.lv.solid.as_ncsg()
1499     cn.dump()
1500
1501     v = n.children[0]
1502     for c in v.children:
1503         print(c.pv.transform)
1504
1505
1506 def test_polycone(g):
1507     pcs = g.findall_("solids//polycone")
1508     for i,pc in enumerate(pcs):
1509         print("pcs[%2d].zp : %s " % (i, pc.zp))
1510         cn = pc.as_ncsg()
1511
1512 def test_gdml(g):
1513     print(g.world)
1514     g.world.rdump()
1515     g.volumes["/dd/Geometry/PMT/lvPmtHemi0xc133740"].rdump()
1516
1517 def test_trd(g):
1518     trs = g.findall_("solids//trd")
1519     cns = []
1520     for i,tr in enumerate(trs):
1521         cn = tr.as_ncsg()
1522         cns.append(cn)
1523
1524
1525 def test_gdml_fromstring():
1526     gg = GDML.fromstring(r"""<gdml><solids><tube aunit="deg" deltaphi="360" lunit="mm" name="AdPmtCollar0xc2c5260" rmax="106" rmin="105" startphi="0" z="12.7"/></solids></gdml>""")
1527     so = gg.solids(0)
1528     assert type(so) is Tube
1529     cn = so.as_ncsg()
1530     assert cn.typ == cn.DIFFERENCE
1531
1532 def test_primitive_fromstring():
1533     so = Primitive.fromstring(r"""<tube aunit="deg" deltaphi="360" lunit="mm" name="AdPmtCollar0xc2c5260" rmax="106" rmin="105" startphi="0" z="12.7"/>""")
1534     assert type(so) is Tube
1535     cn = so.as_ncsg()
1536     assert cn.typ == cn.DIFFERENCE
1537
1538
1539
1540 def parse_args(doc, **kwa):
1541     np.set_printoptions(suppress=True, precision=3, linewidth=200)
1542     parser = argparse.ArgumentParser(doc)
1543     parser.add_argument(     "idx",  type=int, nargs="*", help="Node index to dump.")
1544     parser.add_argument(     "--level", default="info", help="logging level" )
1545     parser.add_argument(  "-g","--gdmlpath", default=kwa.get("gdmlpath",None), help="Override the default key.gdmlpath:\n%(default)s " )
1546     args = parser.parse_args()
1547     fmt = '[%(asctime)s] p%(process)s {%(pathname)s:%(lineno)d} %(levelname)s - %(message)s'
1548     logging.basicConfig(level=getattr(logging,args.level.upper()), format=fmt)
1549     if len(args.idx) == 0:
1550         args.idx = [0]
1551     pass
1552     return args
1553
1554
1555 def define_matrix_values(elem, scale=1e6, startswith=None):
1556     log.info("define_matrix_values startswith:%s scale:%s " % (startswith, scale))
1557     mm = elem.findall("define/matrix")
1558     fmt = " {coldim:3s} {a.shape!r:8s} {dig:32s} {name:40s}   "
1559     d = odict()
1560     for m in mm:
1561         name = m.attrib["name"]
1562         if startswith is None or name.startswith(startswith):
1563             coldim = int(m.attrib["coldim"])
1564             values = m.attrib["values"]
1565             a = np_fromstring(values, coldim=coldim, scale=scale)
1566             d[name] = a
1567             dig = np_digest(a)
1568             print(fmt.format(**m.attrib,a=a, dig=dig))
1569
1570
1571             values2 = np_tostring(a, scale=scale)
1572             a2 = np_fromstring(values2, coldim=coldim, scale=scale)
1573             assert np.allclose( a, a2)
1574             ## number formatting changes, but should be no significant change in actual values
1575         pass
1576     pass
1577     if startswith is None:
1578         assert len(d) == len(mm)
1579     pass
1580     return d
1581
1582
1583
1584 if __name__ == '__main__':
1585
1586     key = key_()
1587     args = parse_args(__doc__, gdmlpath=key.gdmlpath)
1588
1589     g = GDML.parse(args.gdmlpath)
1590
1591     from opticks.analytic.treebase import Tree
1592     t = Tree(g.world)
1593
1594     #test_children(t)
1595     #test_polycone(g)
1596     #test_trd(g)
1597     #test_gdml_fromstring()
1598     #test_primitive_fromstring()
1599
1600     mv = define_matrix_values(g.elem, startswith="EFFICIENCY")
1601     #mv = define_matrix_values(g.elem)
1602
1603
1604