ana/pmt/ncsgconverter.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 ana/pmt/ncsgconverter.py
0023 ===========================
0024
0025 Translation of detdesc XML defined anaytic PMT into NCSG input serialization.
0026 This probably predated GDML parsing, and it not much more than a curiosity.
0027 Although this may allow direct comparison of the same PMT geometry
0028 using different generations of analytic approaches.
0029
0030 ::
0031
0032     2645 tboolean-dd(){          TESTCONFIG=$(tboolean-dd- 2>/dev/null)     tboolean-- $* ; }
0033     2646 tboolean-dd-()
0034     2647 {
0035     2648     python $(tboolean-dir)/tboolean_dd.py \
0036     2649           --csgpath $TMP/$FUNCNAME \
0037     2650           --container $(tboolean-container)  \
0038     2651           --testobject $(tboolean-testobject)
0039     2652
0040     2653     # got too long for here-string  so broke out into script
0041     2654 }
0042     2655 tboolean-dd-check(){ tboolean-dd- 2> /dev/null ; }
0043     2656 tboolean-dd-edit(){ vi $(tboolean-dir)/tboolean_dd.py  ; }
0044     2657 tboolean-dd-scan(){ SCAN="0,0,127.9,0,0,1,0,0.1,0.01" NCSGScanTest $TMP/tboolean-dd-/1 ; }
0045     2658
0046
0047 ::
0048
0049     delta:pmt blyth$ opticks-find NCSGConverter
0050     ./ana/pmt/ncsgconverter.py:class NCSGConverter(object):
0051     ./ana/pmt/ncsgconverter.py:    cn = NCSGConverter.ConvertLV( tr.root.lv )
0052     ./tests/tboolean_dd.py:from opticks.ana.pmt.ncsgconverter import NCSGConverter
0053     ./tests/tboolean_dd.py:        obj = NCSGConverter.ConvertLV( lv )
0054     delta:opticks blyth$
0055
0056
0057 """
0058 import os, logging, sys, math, numpy as np
0059 log = logging.getLogger(__name__)
0060
0061 from opticks.ana.base import opticks_main
0062 from opticks.analytic.csg import CSG
0063
0064 from ddbase import Dddb, Sphere, Tubs, Intersection, Union, Subtraction
0065
0066 from opticks.analytic.treebase import Tree
0067
0068
0069
0070 class NCSGConverter(object):
0071     """
0072     Translate single volume detdesc primitives and CSG operations
0073     into an NCSG style node tree
0074     """
0075     @classmethod
0076     def ConvertLV(cls, lv ):
0077         """
0078         :param lv: Elem
0079         :return cn: CSG node instance
0080         """
0081         lvgeom = lv.geometry()
0082         assert len(lvgeom) == 1, "expecting single CSG operator or primitive Elem within LV"
0083
0084         cn = cls.convert(lvgeom[0])
0085
0086         if lv.posXYZ is not None:
0087             assert cn.transform is None, cn.transform
0088             translate  = "%s,%s,%s" % (lv.xyz[0], lv.xyz[1], lv.xyz[2])
0089             cn.translate = translate
0090             log.info("TranslateLV posXYZ:%r -> translate %s  " % (lv.posXYZ, translate) )
0091         pass
0092         return cn
0093
0094
0095     @classmethod
0096     def convert(cls, node):
0097         """
0098         :param node: instance of ddbase.Elem subclass
0099         :return cn: CSG node
0100         """
0101         assert node.is_operator ^ node.is_primitive, "node must either be operator or primitive "
0102         cn = cls.convert_primitive(node) if node.is_primitive else cls.convert_operator(node)
0103         return cn
0104
0105     @classmethod
0106     def convert_Sphere(cls, en, only_inner=False):
0107         """
0108         :param en: source element node
0109         :param use_slab: alternative approach using intersect with a slab rather than the nascent zsphere
0110         :param only_inner: used to control/distinguish internal recursive call handling the inner sphere
0111
0112         Prior to implementing zsphere with caps, tried using infinite slab
0113         in boolean intersection with the sphere to effect the zslice.
0114         But this approach unavoidably yields a cap, as switching off the
0115         slab caps causes the intersection with the slab to yield nothing.
0116         Hence proceeded to implement zsphere with cap handling.
0117
0118         * z-slice sphere primitive OR intersect with a slab ?
0119         * r-range sphere primitive OR difference two spheres ?
0120
0121         * doing z and r both at once is problematic for param layout
0122         """
0123         outerRadius = en.outerRadius.value
0124         innerRadius = en.innerRadius.value
0125         x = en.xyz[0]
0126         y = en.xyz[1]
0127         z = en.xyz[2]
0128
0129         has_inner = not only_inner and innerRadius is not None
0130         if has_inner:
0131             inner = cls.convert_Sphere(en, only_inner=True)  # recursive call to make inner sphere
0132         pass
0133
0134         radius = innerRadius if only_inner else outerRadius
0135         assert radius, (radius, innerRadius, outerRadius, only_inner)
0136
0137
0138         startThetaAngle = en.startThetaAngle.value
0139         deltaThetaAngle = en.deltaThetaAngle.value
0140
0141         log.info("convert_Sphere outerRadius:%s innerRadius:%s radius:%s only_inner:%s  has_inner:%s " % (outerRadius,innerRadius,radius, only_inner, has_inner))
0142
0143         zslice = startThetaAngle is not None or deltaThetaAngle is not None
0144
0145         if zslice:
0146             if startThetaAngle is None:
0147                 startThetaAngle = 0.
0148
0149             if deltaThetaAngle is None:
0150                 deltaThetaAngle = 180.
0151
0152
0153
0154             # z to the right, theta   0 -> z=r, theta 180 -> z=-r
0155             rTheta = startThetaAngle
0156             lTheta = startThetaAngle + deltaThetaAngle
0157
0158             assert rTheta >= 0. and rTheta <= 180.
0159             assert lTheta >= 0. and lTheta <= 180.
0160             zmin = radius*math.cos(lTheta*math.pi/180.)
0161             zmax = radius*math.cos(rTheta*math.pi/180.)
0162             assert zmax > zmin, (startThetaAngle, deltaThetaAngle, rTheta, lTheta, zmin, zmax )
0163
0164             log.info("convert_Sphere 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 ))
0165
0166             cn = CSG("zsphere", name=en.name, param=[x,y,z,radius], param1=[zmin,zmax,0,0], param2=[0,0,0,0]  )
0167
0168             ZSPHERE_QCAP = 0x1 << 1   # zmax
0169             ZSPHERE_PCAP = 0x1 << 0   # zmin
0170             flags = ZSPHERE_QCAP | ZSPHERE_PCAP
0171
0172             cn.param2.view(np.uint32)[0] = flags
0173             pass
0174         else:
0175             cn = CSG("sphere", name=en.name)
0176             cn.param[0] = x
0177             cn.param[1] = y
0178             cn.param[2] = z
0179             cn.param[3] = radius
0180         pass
0181         if has_inner:
0182             ret = CSG("difference", left=cn, right=inner )
0183             #ret = inner
0184             #ret = cn
0185         else:
0186             ret = cn
0187         pass
0188         return ret
0189
0190
0191
0192     @classmethod
0193     def convert_Tubs(cls, en):
0194         cn = CSG("cylinder", name=en.name)
0195
0196         cn.param[0] = en.xyz[0]
0197         cn.param[1] = en.xyz[1]
0198         cn.param[2] = 0    ## hmm this z is ignored in NCylinder
0199         cn.param[3] = en.outerRadius.value
0200
0201         zoff = en.xyz[2]
0202
0203         hz = en.sizeZ.value/2.
0204         cn.param1[0] = -hz + zoff
0205         cn.param1[1] =  hz + zoff   # direct offsetting or use transform ?
0206
0207
0208         PCAP = 0x1 << 0
0209         QCAP = 0x1 << 1
0210         flags = PCAP | QCAP
0211         cn.param1.view(np.uint32)[1] = flags
0212
0213         return cn
0214
0215
0216     @classmethod
0217     def convert_primitive(cls, en):
0218         convert_method_name = "convert_%s" % en.__class__.__name__
0219         convert_method = getattr(cls, convert_method_name, None )
0220         assert convert_method, "missing convert method: %s " % convert_method_name
0221         #log.info("convert_primitive with %s " % convert_method_name )
0222         cn = convert_method(en)
0223         cn.elem = en   # <-- temporary during dev, not used downstream
0224         return cn
0225
0226     @classmethod
0227     def convert_operator(cls, en):
0228         """
0229         Source Elem xml tree CSG operator nodes with three children
0230         have to be divided up to fit into binary CSG tree::
0231
0232
0233                    1
0234                   / \
0235                  10  11
0236                 /  \
0237                100 101
0238
0239         """
0240         op = en.__class__.__name__.lower()
0241         assert op in ["intersection", "union", "difference"]
0242
0243         children = en.geometry()
0244         nchild = len(children)
0245
0246         if nchild == 2:
0247
0248             cn = CSG(op, name=en.name)
0249             cn.left = cls.convert(children[0])
0250             cn.right = cls.convert(children[1])
0251
0252         elif nchild == 3:
0253
0254             cn = CSG(op, name=en.name)
0255
0256             ln = CSG(op, name=en.name + "_split3")
0257             ln.left = cls.convert(children[0])
0258             ln.right = cls.convert(children[1])
0259
0260             cn.left = ln
0261             cn.right = cls.convert(children[2])
0262
0263         else:
0264             assert 0, "CSG operator nodes must have 2 or 3 children"
0265         pass
0266         return cn
0267
0268
0269
0270
0271
0272 if __name__ == '__main__':
0273
0274     args = opticks_main(apmtidx=2)
0275
0276     g = Dddb.parse(args.apmtddpath)
0277
0278
0279     lvn = "lvPmtHemi"
0280     #lvn = "lvPmtHemiwPmtHolder"
0281
0282     lv = g.logvol_(lvn)
0283
0284     tr = Tree(lv)
0285
0286     cn = NCSGConverter.ConvertLV( tr.root.lv )
0287
0288     cn.dump()
0289
0290
0291
0292
0293
0294
0295
0296
0297
0298
0299
0300