eic-opticks/ana/x018.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 Aiming for this to be generated, so keep it simple
0023 """
0024 import numpy as np
0025 from opticks.ana.shape import X, SEllipsoid, STubs, STorus, SCons, SSubtractionSolid, SUnionSolid, SIntersectionSolid
0026
0027
0028 class x018(X):
0029     """
0030     G4VSolid* make_solid()
0031     {
0032         G4ThreeVector A(0.000000,0.000000,-23.772510);
0033         G4ThreeVector B(0.000000,0.000000,-195.227490);
0034         G4ThreeVector C(0.000000,0.000000,-276.500000);
0035         G4ThreeVector D(0.000000,0.000000,92.000000);
0036
0037
0038         G4VSolid* d = new G4Ellipsoid("PMT_20inch_inner_solid_1_Ellipsoid0x4c91130", 249.000000, 249.000000, 179.000000, -179.000000, 179.000000) ; // 3
0039         G4VSolid* g = new G4Tubs("PMT_20inch_inner_solid_2_Tube0x4c91210", 0.000000, 75.951247, 23.782510, 0.000000, CLHEP::twopi) ; // 4
0040         G4VSolid* i = new G4Torus("PMT_20inch_inner_solid_2_Torus0x4c91340", 0.000000, 52.010000, 97.000000, -0.000175, CLHEP::twopi) ; // 4
0041
0042         G4VSolid* f = new G4SubtractionSolid("PMT_20inch_inner_solid_part20x4cb2d80", g, i, NULL, A) ; // 3
0043         G4VSolid* c = new G4UnionSolid("PMT_20inch_inner_solid_1_20x4cb30f0", d, f, NULL, B) ; // 2
0044
0045         G4VSolid* k = new G4Tubs("PMT_20inch_inner_solid_3_EndTube0x4cb2fc0", 0.000000, 45.010000, 57.510000, 0.000000, CLHEP::twopi) ; // 2
0046                                                                          //    rmin,    rmax,      dz,        sphi,     dphi
0047
0048         G4VSolid* b = new G4UnionSolid("PMT_20inch_inner_solid0x4cb32e0", c, k, NULL, C) ; // 1
0049         G4VSolid* m = new G4Tubs("Inner_Separator0x4cb3530", 0.000000, 254.000000, 92.000000, 0.000000, CLHEP::twopi) ; // 1
0050
0051         G4VSolid* a = new G4IntersectionSolid("PMT_20inch_inner1_solid0x4cb3610", b, m, NULL, D) ; // 0
0052         return a ;
0053     }
0054
0055
0056                            a                                            Intersection
0057                                                                        /            \
0058                       b             m(D)                          Union             m:Tubs
0059                                                                   /    \
0060                 c          k(C)                             Union       Tubs
0061                                                            /     \
0062             d     f(B)                            Ellipsoid   Subtraction
0063                                                               /          \
0064                 g(B)  i(B+A)                                Tubs         Torus
0065
0066
0067
0068     * intersection with m:Tubs blows away the rest of the tree leaving
0069       just the top half of the ellipsoid
0070
0071     * x019 is almost identical to this just with "Intersection" -> "Subtraction"
0072       so for that root.left becomes root
0073
0074
0075     """
0076     def __init__(self, mode=0):
0077
0078         d = SEllipsoid("d", [249.000, 179.000 ] )
0079         g = STubs("g", [75.951247,23.782510] )
0080         i = STorus("i", [ 52.010000, 97.000000] )
0081
0082         A = np.array( [0, -23.772510] )
0083         f = SSubtractionSolid("f", [g,i,A] )
0084         B = np.array( [0, -195.227490] )
0085         c = SUnionSolid("c", [d,f,B] )
0086
0087         k = STubs("k", [45.010000, 57.510000] )
0088
0089         C = np.array( [0, -276.500000] )
0090         b = SUnionSolid("b", [c,k,C] )
0091         m = STubs("m", [254.000000, 92.000000] )
0092
0093         D = np.array( [0, 92.000000] )
0094         a = SIntersectionSolid("a", [b,m,D] )
0095
0096         X.__init__(self, a )   # sets root
0097
0098
0099 if __name__ == '__main__':
0100     x = x018()
0101     print(repr(x))
0102
0103
0104     i = x.find_one("STorus")
0105     r = i.param[0]
0106     R = i.param[1]
0107
0108     d = x.find_one("SEllipsoid")
0109     ex = d.param[0]
0110     ez = d.param[1]
0111
0112     print("r %s R %s ex %s ez %s " % (r,R,ex,ez))
0113
0114
0115     print(" SEllipsoid d.xy %s " % repr(d.xy) )
0116     print(" STorus     i.xy %s " % repr(i.xy) )
0117
0118     z0 = i.xy[1]
0119
0120     torus_rhs = [R,z0]
0121
0122
0123     print(" z0  %s " % z0 )
0124     print(" torus_rhs  %s " % repr(torus_rhs) )
0125
0126     from opticks.ana.shape import ellipse_closest_approach_to_point
0127     p = ellipse_closest_approach_to_point( ex, ez, torus_rhs  )
0128
0129     pr, pz = p    # at torus/ellipse closest point : no guarantee of intersection
0130     print(" p %s " % repr(p) )
0131
0132
0133     r2 = pr
0134     r1 = R - r
0135     mz = (z0 + pz)/2.   # mid-z cone coordinate (ellipsoid frame)
0136     hz = (pz - z0)/2.   # cons galf height
0137
0138
0139     f = SCons( "f", [r1,r2,hz] )
0140     B = np.array( [0, mz] )
0141
0142     print(" replacment SCons %s offset %s " % (repr(f),repr(B)))
0143
0144     # hmm how to do tree surgery
0145
0146