ana/pmt/ddpart.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 ddpart.py
0023 ===========
0024
0025 Approach:
0026
0027 * kicking off from Dddb.logvol_ wraps results of
0028   lxml find into do-everything Elem subclass appropriate to the tags of the elements
0029
0030 * primary purpose is the Elem.parts method, trace back from there to understand
0031
0032
0033 Crucial Methods
0034 -----------------
0035
0036 ElemPartitioner.parts
0037      nexus of control of the partitioning, which gets manually mixed into Elem
0038
0039
0040
0041
0042 Need to play some mixin tricks
0043
0044 * http://stackoverflow.com/questions/8544983/dynamically-mixin-a-base-class-to-an-instance-in-python
0045
0046
0047 """
0048 import os, re, logging, math
0049 from opticks.ana.base import opticks_main, manual_mixin
0050
0051
0052 from ddbase import Dddb
0053 from ddbase import Elem, Sphere, Tubs, Primitive, Physvol, Logvol, PosXYZ
0054
0055 from geom import Part, BBox, ZPlane, Rev
0056 from gcsg import GCSG
0057
0058 X,Y,Z = 0,1,2
0059
0060
0061 log = logging.getLogger(__name__)
0062
0063
0064 class Parts(list):
0065    def __repr__(self):
0066        return "Parts " + list.__repr__(self)
0067
0068    def __init__(self, *args, **kwa):
0069        list.__init__(self, *args, **kwa)
0070
0071
0072 class Uncoincide(object):
0073     def __init__(self, top="CONTAINING_MATERIAL", sensor="SENSOR_SURFACE"):
0074          self.top = top
0075          self.sensor = sensor
0076          pass
0077     def face_bottom_tubs(self, name):
0078         """
0079         Only face and bottom have coincident surfaces in literal
0080         translation that need fixing.
0081
0082         ::
0083
0084                tubs                      bottom                    face
0085
0086                OM///Pyrex                OM///Pyrex               OM///Pyrex
0087
0088                Pyrex///Vacuum            Pyrex///OpaqueVacuum     Pyrex/SENSOR//Bialkali
0089
0090                Vacuum///OpaqueVacuum     OpaqueVacuum///Vacuum    Bialkali///Vacuum
0091
0092
0093         """
0094         log.warning("UNCOINCIDE boundary setting for %s " % name)
0095         face = bottom = tubs = None
0096         if name == "pmt-hemi":
0097             pass
0098             boundary = "%s///Pyrex" % self.top
0099             face = boundary
0100             bottom = boundary
0101             tubs = boundary
0102             pass
0103         elif name == "pmt-hemi-vac":
0104             pass
0105             face = "Pyrex/%s//Bialkali" % self.sensor
0106             bottom = "Pyrex///OpaqueVacuum"
0107             tubs = "Pyrex///Vacuum"
0108             pass
0109         elif name == "pmt-hemi-cathode":
0110             pass
0111             face = "Bialkali///Vacuum"
0112             bottom = None
0113             tubs = None
0114             pass
0115         elif name == "pmt-hemi-bot":
0116             pass
0117             face = None
0118             bottom = "OpaqueVacuum///Vacuum"
0119             tubs = None
0120             pass
0121         elif name == "pmt-hemi-dynode":
0122             pass
0123             face = None
0124             bottom = None
0125             tubs = "Vacuum///OpaqueVacuum"
0126             pass
0127         else:
0128             assert 0
0129         pass
0130         return face, bottom, tubs
0131
0132 UNCOINCIDE = Uncoincide()
0133
0134
0135 # intended to be a mixin to Elem
0136 class ElemPartitioner(object):
0137     def partition_intersection_3spheres(self, spheres, material=None):
0138         """
0139         :param spheres:  list of three *Sphere* in ascending center z order, which are assumed to intersect
0140         :return parts:  list of three sphere *Part* instances
0141
0142         Consider splitting the lens shape made from the boolean CSG intersection
0143         of two spheres along the plane of the intersection.
0144         The left part of the lens comes from the right Sphere
0145         and the right part comes left Sphere.
0146
0147         Extending from two sphere intersection to three spheres
0148         numbered s1,s2,s3 from left to right, with two ZPlane intersections z23, z12.
0149
0150         left
0151             from s3, bounded by z23 on right
0152
0153         middle
0154             from s2, bounded by z23 on left, z12 on right
0155
0156         right
0157             from s1, bounded by z12 on left
0158
0159         """
0160         s1, s2, s3 = spheres
0161
0162         assert s1.z < s2.z < s3.z
0163
0164         z12 = Sphere.intersect("z12",s1,s2)   # ZPlane of s1 s2 intersection
0165         z23 = Sphere.intersect("z23",s2,s3)   # ZPlane of s2 s3 intersection
0166
0167         assert z23.z < z12.z
0168
0169         p1 = s3.part_zleft(z23)
0170         p2 = s2.part_zmiddle(z23, z12)
0171         p3 = s1.part_zright(z12)
0172
0173         assert p1.bbox.z < p2.bbox.z < p3.bbox.z
0174
0175         pts = Parts([p3,p2,p1])
0176         return pts
0177
0178     def partition_intersection_2spheres(self, spheres, material=None, unoverlap=True):
0179         """
0180         :param spheres: list of two *Sphere* in ascending center z order
0181         :return parts: list of two sphere *Part* instances
0182
0183         Used for the very thin photocathode
0184
0185         For pmt-hemi-cathode-face_part pmt-hemi-cathode-belly_part the
0186         bbox of the Sphere parts from just the theta angle ranges are overlapping
0187         each other very slightly with the plane of sphere intersection in between.
0188         Because of the thin photocathode this creates an dangerous sharp protuding edge.
0189
0190         Suspect overlapping bbox may be problematic, so unoverlap by adjusting z
0191         edges of bbox, and adjust in xy.
0192         """
0193         s1, s2 = spheres
0194         assert s1.z < s2.z
0195
0196         # parts from just the theta ranges
0197         p1, p2 = Part.ascending_bbox_zmin([s1.as_part(),s2.as_part()])
0198         assert p1.bbox.zmin < p2.bbox.zmin
0199         p12 = Sphere.intersect("p12",s1,s2)
0200
0201         if unoverlap:
0202             p1.bbox.zmax = p12.z
0203             p2.bbox.zmin  = p12.z
0204             p2.bbox.xymax = p12.y
0205             p2.bbox.xymin = -p12.y
0206
0207         log.warning("material %s name %s " % (material, self.name))
0208         log.warning("p1 %s " % p1)
0209         log.warning("p2 %s " % p2)
0210
0211         # pluck inner radii, for the insides
0212         i1, i2 = Part.ascending_bbox_zmin([s1.as_part(inner=True),s2.as_part(inner=True)])
0213         assert i1.bbox.zmin < i2.bbox.zmin
0214         i12 = Sphere.intersect("p12",s1,s2, inner=True)
0215
0216         if unoverlap:
0217             i1.bbox.zmax = i12.z
0218             i2.bbox.zmin  = i12.z
0219             i2.bbox.xymax = i12.y
0220             i2.bbox.xymin = -i12.y
0221
0222         log.warning("i1 %s" % i1)
0223         log.warning("i2 %s" % i2)
0224         if UNCOINCIDE and self.name == "pmt-hemi-cathode":
0225             face, bottom, tubs = UNCOINCIDE.face_bottom_tubs(self.name)
0226             i2.boundary = face
0227             i1.boundary = face
0228             ret = [i2,i1]    # skipping the coincidents
0229         else:
0230             #i1.parent = p1
0231             #i2.parent = p2
0232             ret = [p2,i2,p1,i1]
0233         pass
0234         return Parts(ret)
0235
0236     def partition_intersection(self, material=None):
0237         #log.info(self)
0238
0239         spheres = []
0240         comps = self.findall_("./*")
0241         self.link_prior_posXYZ(comps)
0242
0243         other = []
0244         for c in comps:
0245             if type(c) is Sphere:
0246                 spheres.append(c)
0247             elif type(c) is PosXYZ:
0248                 pass
0249             else:
0250                 other.append(c)
0251             pass
0252
0253         assert len(other) == 0, ("only 2 or 3-sphere intersections handled", other )
0254
0255         for i,s in enumerate(spheres):
0256             log.debug("s%d: %s %s " % (i, s.desc, s.outerRadius.value))
0257
0258         if len(spheres) == 3:
0259             pts = self.partition_intersection_3spheres(spheres, material=material)
0260         elif len(spheres) == 2:
0261             pts = self.partition_intersection_2spheres(spheres, material=material)
0262         else:
0263             assert 0
0264
0265
0266         pts.gcsg = GCSG(self, spheres) # so pts keep hold of reference to the basis Elem and shapes list they came from
0267
0268         return pts
0269
0270
0271     def partition_union_intersection_tubs(self, comps, material=None, verbose=False):
0272         """
0273         """
0274         log.info("material %s name %s " % (material, self.name))
0275
0276         ipts = comps[0].partition_intersection()
0277         sparts = Part.ascending_bbox_zmin(ipts)
0278         assert len(sparts) == 3
0279
0280         tpart = comps[1].as_part()
0281         ts = Part.intersect_tubs_sphere("ts", tpart, sparts[0], -1)   # -ve root for leftmost
0282
0283
0284         if verbose:
0285             log.info(self)
0286             log.info("ts %s " % repr(ts))
0287             for i, s in enumerate(sparts):
0288                 log.info("sp(%s) %s " % (i, repr(s)))
0289             log.info("tp(0) %s " % repr(tpart))
0290         pass
0291
0292         # CRUCIAL setting of same zmax for ztubs and zmin for zsphere
0293         # done via the bbox...
0294         # TODO: store z-range in parameters (2*float), not the costly bbox (6*float)
0295
0296         tpart.bbox.zmax = ts.z
0297         sparts[0].bbox.zmin = ts.z
0298
0299         tpart.enable_endcap("P")  # smaller Z endcap
0300
0301
0302         if UNCOINCIDE:
0303             if self.name == "pmt-hemi" or self.name == "pmt-hemi-vac":
0304                 face, bottom, tubs = UNCOINCIDE.face_bottom_tubs(self.name)
0305             else:
0306                 assert 0
0307             pass
0308             tpart.boundary = tubs
0309             sparts[0].boundary = bottom
0310             sparts[1].boundary = face
0311             sparts[2].boundary = face
0312         pass
0313
0314         rparts = Parts()
0315         rparts.extend(sparts)
0316         rparts.extend([tpart])
0317
0318         # had CSG list/tree confusion here, from thinking ahead to serialization
0319         # dont do that, use the best representation for the immediate problem at hand
0320         # of grabbing the CSG tree into an in-memory representation
0321
0322         rparts.gcsg = GCSG(self, [ipts.gcsg, comps[1]] )
0323
0324         return rparts
0325
0326
0327     def partition_union(self, material=None, verbose=False):
0328         """
0329         union of a 3-sphere lens shape and a tubs requires:
0330
0331         * adjust bbox of the abutting part Sphere to the intersection z of tubs and Sphere
0332         * avoid a surface at the interface of tubs endcap and part Sphere
0333
0334         """
0335         comps = self.findall_("./*")
0336         self.link_prior_posXYZ(comps)
0337
0338         rparts = Parts()
0339         xret = None
0340
0341         if len(comps) == 3 and comps[0].is_intersection and comps[1].is_tubs and comps[2].is_posXYZ:
0342             # pmt-hemi-glass-bulb, pmt-hemi-bulb-vac
0343             xret = self.partition_union_intersection_tubs(comps[0:3], material=material)
0344         elif len(comps) == 3 and comps[0].is_sphere and comps[1].is_sphere and comps[2].is_posXYZ:
0345             xret = self.partition_intersection_2spheres(comps[0:2], material=material)
0346             if not hasattr(xret, 'gcsg'):
0347                 xret.gcsg = GCSG(self, comps[0:2])
0348         else:
0349             xret = self.parts()
0350
0351         if xret is not None:
0352             rparts.extend(xret)
0353             if hasattr(xret, 'gcsg'):
0354                 rparts.gcsg = xret.gcsg
0355         pass
0356         return rparts  ;
0357
0358     def parts_sphere_with_inner(self, c):
0359         log.info("name %s " % c.name)
0360
0361         p = c.as_part()
0362         i = c.as_part(inner=True)
0363         log.info("    part   %s " % p )
0364         log.info("    inner  %s " % i )
0365
0366         if UNCOINCIDE and c.name == "pmt-hemi-bot":
0367             face, bottom, tubs = UNCOINCIDE.face_bottom_tubs(c.name)
0368             i.boundary = bottom
0369             ret = [i]
0370         else:
0371             #i.parent = p
0372             ret = [p,i]
0373         pass
0374         pts = Parts(ret)
0375         pts.gcsg = GCSG(c)
0376         return pts
0377
0378
0379     def parts_primitive(self, c):
0380         log.info("name %s " % c.name)
0381         p = c.as_part()
0382         if UNCOINCIDE and c.name == "pmt-hemi-dynode":
0383             face, bottom, tubs = UNCOINCIDE.face_bottom_tubs(c.name)
0384             p.boundary = tubs
0385         pass
0386         pts = Parts([p])
0387         pts.gcsg = GCSG(c)
0388         return pts
0389
0390
0391     def parts(self, analytic_version=0):
0392         """
0393         :return: list of Part instances
0394
0395         Provides parts from a single LV only, ie not
0396         following pv refs. Recursion is needed
0397         in order to do link posXYZ transforms with geometry
0398         and skip them from the parts returned.
0399         """
0400         if type(self) is Physvol:
0401             return []
0402
0403         if type(self) is Logvol and not self.is_logvol:
0404             log.warning("inconsistent LV %s " % repr(self))
0405
0406         if self.is_logvol:
0407             base = self.posXYZ
0408             material = self.material
0409         else:
0410             base = None
0411             material = None
0412         pass
0413
0414         comps = self.findall_("./*")  # one lev only
0415         self.link_prior_posXYZ(comps, base)
0416
0417         if not base is None:
0418             log.info("&&&&&&&&&&&&&&&& name:%s  ncomp:%d c0:%r ElemPartitioner.parts base translate %r " % (self.name, len(comps), comps[0],  base) )
0419         pass
0420
0421         rparts = Parts()
0422
0423         gcsg_ = []
0424
0425         for c in comps:
0426             xret = []
0427             if c.is_sphere and c.has_inner():
0428                 log.info("-> sphere with non-zero inner radius  %s " % c.name)
0429                 xret = self.parts_sphere_with_inner(c)
0430             elif c.is_primitive:
0431                 log.info("-> primitive %s " % c.name)
0432                 xret = self.parts_primitive(c)
0433             elif c.is_intersection:
0434                 log.info("-> intersection %s " % c.name)
0435                 xret = c.partition_intersection(material=material)
0436             elif c.is_union:
0437                 log.info("-> union %s " % c.name)
0438                 xret = c.partition_union(material=material)
0439             elif c.is_composite:
0440                 log.info("-> composite %s " % c.name )
0441                 xret = c.parts()
0442             elif c.is_posXYZ:
0443                 pass
0444             else:
0445                 log.warning("skipped component %s " % repr(c))
0446             pass
0447
0448             if len(xret) > 0:
0449                 rparts.extend(xret)
0450                 if hasattr(xret, 'gcsg'):
0451                     gcsg_.append(xret.gcsg)
0452                 pass
0453             pass
0454         pass
0455
0456         log.info("%s : %s comps yield %s rparts %s gcsg " % (repr(self), len(comps), len(rparts), len(gcsg_)))
0457
0458         if self.is_logvol:
0459             for p in rparts:
0460                 if p.material is None:
0461                     p.material = self.material
0462                 pass
0463             pass
0464             for cn in gcsg_:
0465                 # associate lv to the first gcsg node
0466                 if cn.lv is None:
0467                     cn.lv = self
0468                     #log.info("cn lv : %s " % repr(cn.lv))
0469
0470
0471         rparts.gcsg = gcsg_
0472
0473         return rparts
0474
0475     def geometry(self):
0476         return filter(lambda c:c.is_geometry, self.components())
0477
0478
0479
0480
0481
0482
0483 class SpherePartitioner(object):
0484
0485     @classmethod
0486     def intersect(cls, name, a_, b_, inner=False):
0487         """
0488         Find Z intersect of two Z offset spheres
0489
0490         * http://mathworld.wolfram.com/Circle-CircleIntersection.html
0491
0492         :param name: identifer passed to ZPlane
0493         :param a_: *Sphere* instance
0494         :param b_: *Sphere* instance
0495
0496         :return zpl: *ZPlane* instance with z and y attributes where:
0497                      z is intersection coordinate
0498                      y is radius of the intersection circle
0499         """
0500
0501         if inner:
0502             R = a_.innerRadius.value
0503             r = b_.innerRadius.value
0504         else:
0505             R = a_.outerRadius.value
0506             r = b_.outerRadius.value
0507         pass
0508         assert R is not None and r is not None
0509
0510         a = a_.xyz
0511         b = b_.xyz
0512
0513         log.debug(" R %s a %s " % ( R, repr(a)) )
0514         log.debug(" r %s b %s " % ( r, repr(b)) )
0515
0516         dx = b[X] - a[X]
0517         dy = b[Y] - a[Y]
0518         dz = b[Z] - a[Z]
0519
0520         assert dx == 0
0521         assert dy == 0
0522         assert dz != 0
0523
0524         d = dz             # use Sphere a_ frame
0525         dd_m_rr_p_RR = d*d - r*r + R*R
0526         z = dd_m_rr_p_RR/(2.*d)
0527         yy = (4.*d*d*R*R - dd_m_rr_p_RR*dd_m_rr_p_RR)/(4.*d*d)
0528         y = math.sqrt(yy)
0529
0530         # add a[Z] to return to original frame
0531         return ZPlane(name, z+a[Z], y )
0532
0533
0534     def as_gcsg(self):
0535         R = self.outerRadius.value
0536         if self.has_inner():
0537             r = self.innerRadius.value
0538         else:
0539             r = -1.0
0540         pass
0541
0542         sta = self.startThetaAngle.value
0543         dta = self.deltaThetaAngle.value
0544         if sta is None:
0545             sta = 0.
0546         if dta is None:
0547             dta = 180.
0548
0549         gcsg = []
0550         gcsg.append( [self.xyz[0], self.xyz[1], self.xyz[2], R] )
0551         gcsg.append( [sta,dta,0,r] )
0552         gcsg.append( [0,0,0,0] )
0553         gcsg.append( [0,0,0,0] )
0554         return gcsg
0555
0556     def as_part(self, inner=False):
0557         """
0558         Spheres sliced via startThetaAngle and deltaThetaAngle
0559         are handled by setting of the bbox
0560         """
0561         if inner:
0562             radius = self.innerRadius.value
0563         else:
0564             radius = self.outerRadius.value
0565         pass
0566         if radius is None:
0567             return None
0568
0569         sta = self.startThetaAngle.value
0570         dta = self.deltaThetaAngle.value
0571
0572         assert self.xyz[2] == self.z
0573         z = self.z
0574
0575         p = Part('Sphere', self.name + "_part", self.xyz, radius )
0576
0577         if sta is None and dta is None:
0578             thetacut = False
0579             bb = self.bbox(z-radius, z+radius, -radius, radius)
0580         else:
0581             # above case is equivalent to sta=0 dta=180
0582             if sta is None:
0583                 sta = 0.
0584             rta = sta
0585             lta = sta + dta
0586             log.debug("Sphere.as_part %s  leftThetaAngle %s rightThetaAngle %s " % (self.name, lta, rta))
0587             assert rta >= 0. and rta <= 180.
0588             assert lta >= 0. and lta <= 180.
0589             zl = radius*math.cos(lta*math.pi/180.)
0590             yl = radius*math.sin(lta*math.pi/180.)
0591             zr = radius*math.cos(rta*math.pi/180.)
0592             yr = radius*math.sin(rta*math.pi/180.)
0593             ym = max(abs(yl),abs(yr))
0594             bb = self.bbox(z+zl, z+zr, -ym, ym)
0595         pass
0596         p.bbox = bb
0597
0598         return p
0599
0600     def part_zleft(self, zpl):
0601         radius = self.outerRadius.value
0602         z = self.xyz[2]
0603         ymax = zpl.y
0604         p = Part('Sphere', self.name + "_part_zleft", self.xyz, radius )
0605         p.bbox = self.bbox(z-radius, zpl.z, -ymax, ymax)
0606         return p
0607
0608     def part_zright(self, zpr):
0609         radius = self.outerRadius.value
0610         z = self.xyz[2]
0611         ymax = zpr.y
0612         p = Part('Sphere', self.name + "_part_zright", self.xyz, radius )
0613         p.bbox = self.bbox(zpr.z,z+radius, -ymax, ymax)
0614         return p
0615
0616     def part_zmiddle(self, zpl, zpr):
0617         p = Part('Sphere', self.name + "_part_zmiddle", self.xyz, self.outerRadius.value )
0618         ymax = max(zpl.y,zpr.y)
0619         p.bbox = self.bbox(zpl.z,zpr.z,-ymax,ymax )
0620         return p
0621
0622     def asrev(self):
0623         xyz = self.xyz
0624         ro = self.outerRadius.value
0625         ri = self.innerRadius.value
0626         st = self.startThetaAngle.value
0627         dt = self.deltaThetaAngle.value
0628         sz = None
0629         wi = None
0630         if ri is not None and ri > 0:
0631             wi = ro - ri
0632
0633         return [Rev('Sphere', self.name,xyz, ro, sz, st, dt, wi)]
0634
0635
0636
0637 class TubsPartitioner(object):
0638     def as_part(self):
0639         sizeZ = self.sizeZ.value
0640         radius = self.outerRadius.value
0641         z = self.xyz[2]
0642         p = Part('Tubs', self.name + "_part", self.xyz, radius, sizeZ )
0643         p.bbox = self.bbox(z-sizeZ/2, z+sizeZ/2, -radius, radius)
0644         return p
0645
0646     def as_gcsg(self):
0647         sizeZ = self.sizeZ.value
0648         outer = self.outerRadius.value
0649         inner = self.innerRadius.value
0650         gcsg = []
0651         gcsg.append( [self.xyz[0], self.xyz[1], self.xyz[2], outer] )
0652         gcsg.append( [0,0,sizeZ, inner] )
0653         gcsg.append( [0,0,0,0] )
0654         gcsg.append( [0,0,0,0] )
0655         return gcsg
0656
0657     def asrev(self):
0658         sz = self.sizeZ.value
0659         r = self.outerRadius.value
0660         xyz = self.xyz
0661         return [Rev('Tubs', self.name,xyz, r, sz )]
0662
0663
0664 class PrimitivePartitioner(object):
0665     def bbox(self, zl, zr, yn, yp ):
0666         assert yn < 0 and yp > 0 and zr > zl
0667         return BBox([yn,yn,zl], [yp,yp,zr])
0668
0669
0670
0671 def ddpart_manual_mixin():
0672     """
0673     Using manual mixin approach to avoid changing
0674     the class hierarchy whilst still splitting base
0675     functionality from partitioner methods.
0676     """
0677     manual_mixin(Tubs, TubsPartitioner)
0678     manual_mixin(Sphere, SpherePartitioner)
0679     manual_mixin(Elem, ElemPartitioner)
0680     manual_mixin(Primitive, PrimitivePartitioner)
0681
0682
0683
0684 if __name__ == '__main__':
0685     args = opticks_main(apmtidx=2)
0686
0687     ddpart_manual_mixin()  # add methods to Tubs, Sphere, Elem and Primitive
0688
0689     g = Dddb.parse(args.apmtddpath)
0690
0691     g.dump_context('PmtHemi')
0692
0693     lv = g.logvol_("lvPmtHemi")
0694
0695     pts = lv.parts()
0696
0697     for pt in pts:
0698         log.info(pt)
0699
0700
0701