ana/pmt/geom.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 import os, logging, math
0022 import numpy as np
0023
0024 from opticks.sysrap.OpticksCSG import CSG_
0025
0026 #TYPECODE = {'Sphere':1, 'Tubs':2, 'Box':3 }  ## equivalent to pre-unified hardcoded and duplicitous approach
0027 TYPECODE = {'Sphere':CSG_.SPHERE, 'Tubs':CSG_.TUBS, 'Box':CSG_.BOX }
0028
0029
0030
0031 X,Y,Z = 0,1,2
0032
0033 log = logging.getLogger(__name__)
0034
0035 class Rev(object):
0036     def __init__(self, typ, name, xyz, radius, sizeZ=None, startTheta=None, deltaTheta=None, width=None):
0037         self.typ = typ
0038         self.name = name
0039         self.xyz = xyz
0040         self.radius = radius
0041         self.sizeZ = sizeZ
0042         self.startTheta = startTheta
0043         self.deltaTheta = deltaTheta
0044         self.width = width
0045
0046     def __repr__(self):
0047         return "Rev('%s','%s' xyz:%s, r:%s, sz:%s, st:%s, dt:%s wi:%s)" % \
0048             (self.typ, self.name, self.xyz, self.radius, self.sizeZ, self.startTheta, self.deltaTheta, self.width)
0049
0050
0051 class ZPlane(object):
0052     def __init__(self, name, z, y):
0053         """
0054         xy symmetry assumed
0055         """
0056         self.name = name
0057         self.z = z
0058         self.y = y
0059     def __repr__(self):
0060         return "ZPlane %s z:%s y:%s " % (self.name, self.z, self.y )
0061
0062 class Part(object):
0063
0064     @classmethod
0065     def ascending_bbox_zmin(cls, parts):
0066         return sorted(parts, key=lambda p:p.bbox.zmin)
0067
0068     @classmethod
0069     def intersect_tubs_sphere(cls, name, tubs, sphere, sign ):
0070         """
0071         :param name: identifier of ZPlane created
0072         :param tubs: tubs Part instance
0073         :param sphere: sphere Part instance
0074         :param sign: 1 or -1 sign of the sqrt
0075         :return zplane: ZPlane instance representing disc intersection with z-position and radius
0076
0077         Sphere at zp on Z axis
0078
0079             xx + yy + (z-zp)(z-zp) = RR
0080
0081         Cylinder along Z axis from -sizeZ/2 to sizeZ/2
0082
0083             xx + yy = rr
0084
0085         Intersection is a circle in Z plane
0086
0087             (z-zp) = sqrt(RR - rr)
0088
0089         """
0090         R = sphere.radius
0091         r = tubs.radius
0092
0093         RR_m_rr = R*R - r*r
0094         assert RR_m_rr > 0
0095
0096         iz = math.sqrt(RR_m_rr)
0097         assert iz < tubs.sizeZ
0098
0099         zp = sphere.xyz[Z]
0100
0101         return ZPlane(name, zp + sign*iz, r)
0102
0103
0104     def enable_endcap(self, tag):
0105         ENDCAP_P = 0x1 <<  0
0106         ENDCAP_Q = 0x1 <<  1
0107         pass
0108         if tag == "P":
0109             self.flags |= ENDCAP_P
0110         elif tag == "Q":
0111             self.flags |= ENDCAP_Q
0112         else:
0113             log.warning("tag is not P or Q, for the low Z endcap (P) and higher Z endcap (Q)")
0114
0115
0116     def __init__(self, typ, name, xyz, radius, sizeZ=0.):
0117         """
0118         :param typ: typ string eg Sphere, Tubs, Box
0119
0120         Canonical instanciations of Part arise in the as_part
0121         methods of dd.py:Primitive subclasses,
0122         currently only Sphere and Tubs.
0123
0124         Geometry defined here is used in::
0125
0126             oxrap/cu/intersect_analytic.cu
0127             oxrap/cu/intersect_ztubs.h etc..
0128
0129         """
0130         self.typ = typ
0131         self.name = name
0132         self.xyz = xyz
0133         self.radius = radius
0134         self.sizeZ = sizeZ   # used for Tubs
0135         self.bbox = None
0136         self.parent = None
0137         self.node = None
0138         self.material = None
0139         self.boundary = None
0140
0141         self.flags = 0
0142         # Tubs endcap control
0143
0144         typecode = TYPECODE.get(typ, -1)
0145         assert typecode > -1, ("unhandled typ ", typ)
0146         self.typecode = typecode
0147
0148
0149     @classmethod
0150     def make_container(cls, parts, factor=3. ):
0151         """
0152         create container box for all the parts
0153         optionally enlarged by a multiple of the bbox extent
0154         """
0155         bb = BBox([0,0,0],[0,0,0])
0156         for pt in parts:
0157             #print pt
0158             bb.include(pt.bbox)
0159         pass
0160         bb.enlarge(factor)
0161
0162         p = Part('Box', "make_container_box", bb.xyz, 0., 0. )
0163         p.bbox = bb
0164         log.info(p)
0165         return p
0166
0167
0168     def __repr__(self):
0169         return "Part %6s %12s %32s %15s r:%6s sz:%5s %40s %s" % (self.typ, self.material, self.name, repr(self.xyz), self.radius, self.sizeZ, repr(self.bbox), self.boundary)
0170
0171     def as_quads(self):
0172         quads = []
0173         quads.append( [self.xyz[0], self.xyz[1], self.xyz[2], self.radius] )
0174         quads.append( [self.sizeZ, 0, 0, 0] )
0175         for q in self.bbox.as_quads():
0176             quads.append(q)
0177         pass
0178         return quads
0179
0180
0181 class BBox(object):
0182     def __init__(self, min_, max_):
0183         self.min_ = np.array(min_)
0184         self.max_ = np.array(max_)
0185
0186     def include(self, other):
0187         """
0188         Expand this bounding box to encompass another
0189         """
0190         self.min_ = np.minimum(other.min_, self.min_)
0191         self.max_ = np.maximum(other.max_, self.max_)
0192
0193     def enlarge(self, factor):
0194         """
0195         Intended to duplicate ggeo-/GVector.hh/gbbox::enlarge
0196         """
0197         dim = self.max_ - self.min_
0198         ext = dim.max()/2.0
0199         vec = np.repeat(ext*factor, 3)
0200         self.min_ = self.min_ - vec
0201         self.max_ = self.max_ + vec
0202
0203     def _get_zmin(self):
0204         return self.min_[Z]
0205     def _set_zmin(self, val):
0206         self.min_[Z] = val
0207     zmin = property(_get_zmin, _set_zmin)
0208
0209     def _get_zmax(self):
0210         return self.max_[Z]
0211     def _set_zmax(self, val):
0212         self.max_[Z] = val
0213     zmax = property(_get_zmax, _set_zmax)
0214
0215     def _get_xymin(self):
0216         """
0217         xy symmetry assumed
0218         """
0219         assert self.min_[X] == self.min_[Y]
0220         return self.min_[X]
0221     def _set_xymin(self, val):
0222         self.min_[X] = val
0223         self.min_[Y] = val
0224     xymin = property(_get_xymin, _set_xymin)
0225
0226     def _get_xymax(self):
0227         assert self.max_[X] == self.max_[Y]
0228         return self.max_[X]
0229     def _set_xymax(self, val):
0230         self.max_[X] = val
0231         self.max_[Y] = val
0232     xymax = property(_get_xymax, _set_xymax)
0233
0234     x = property(lambda self:self.xyz[X])
0235     y = property(lambda self:self.xyz[Y])
0236     z = property(lambda self:self.xyz[Z])
0237     xyz = property(lambda self:(self.min_+self.max_)/2.)
0238
0239     def as_quads(self,scale=1):
0240         qmin = np.zeros(4)
0241         qmin[:3] = self.min_*scale
0242         qmax = np.zeros(4)
0243         qmax[:3] = self.max_*scale
0244         return qmin, qmax
0245
0246     def __repr__(self):
0247         xyz = self.xyz
0248
0249         assert xyz[X] == xyz[Y] == 0.
0250         return "BB %30s %30s z %6.2f" % (str(self.min_), str(self.max_), xyz[Z] )
0251
0252
0253
0254 if __name__ == '__main__':
0255
0256     bb = BBox([-100,-100,-100],[100,100,100])
0257     print bb
0258