eic-opticks/ana/gpmt.py

0001 #!/usr/bin/env python
0002 """
0003 ::
0004
0005     ipython -i gpmt.py
0006
0007
0008     mkdir -p /Users/blyth/simoncblyth.bitbucket.io/env/presentation/ana/gpmt/
0009     cp /tmp/fig/*.png /Users/blyth/simoncblyth.bitbucket.io/env/presentation/ana/gpmt/
0010
0011
0012 https://stackoverflow.com/questions/22959698/distance-from-given-point-to-given-ellipse
0013
0014 """
0015
0016 import os, sys, argparse, logging, textwrap
0017 import numpy as np, math
0018 import matplotlib.pyplot as plt
0019 import matplotlib.lines as mlines
0020
0021 specs_ = lambda s:filter(lambda s:s[0] != "#", filter(None,map(str.strip, textwrap.dedent(s).split("\n"))))
0022
0023 log = logging.getLogger(__name__)
0024 sys.path.insert(0, os.path.expanduser("~"))  # assumes $HOME/opticks
0025
0026 from opticks.analytic.GDML import GDML
0027 from opticks.ana.shape import ellipse_points, circle_points
0028 from opticks.ana.gplt import GPlot, add_line
0029 from opticks.ana.gargs import GArgs
0030
0031 from j.PMTEfficiencyCheck_ import PMTEfficiencyCheck_
0032
0033
0034 if __name__ == '__main__':
0035     args = GArgs.parse(__doc__)
0036     path = args.gdmlpath(4)    # 3:origin_CGDMLKludge_jun15  4:origin_CGDMLKludge_jun28
0037     path = os.path.expandvars(path)
0038
0039     if not os.path.exists(path):
0040         log.fatal("GDML path does not exist: %s " % path )
0041         sys.exit(1)
0042     else:
0043         log.info("parsing GDML path:%s " % path)
0044     pass
0045
0046     g = GDML.parse(path)
0047     g.smry()
0048
0049     NNVT = 1
0050     HAMA = 2
0051     pmt = NNVT
0052     #pmt = HAMA
0053
0054
0055     pec = None
0056     #pec = PMTEfficiencyCheck_()
0057     if not pec is None:
0058         ipec = { HAMA:0,NNVT:1 }
0059         pec_sli = 10000
0060     pass
0061
0062
0063     lvx = args.lvname(pmt)
0064     lv = g.find_one_volume(lvx)
0065
0066     if lv == None:
0067         log.fatal("failed to find lvx:[%s] " % (lvx))
0068     assert lv
0069
0070     #s = lv.solid
0071     #s.sub_traverse()
0072
0073     log.info( "lv %r" % lv )
0074
0075     lvs = g.get_traversed_volumes( lv, maxdepth=args.maxdepth )
0076
0077
0078     log.info( "lvs %r" % lvs )
0079     figpath_suffix  = lv.local_prefix.replace("/","_")
0080     log.info( "lv.local_prefix %s" % lv.local_prefix )
0081
0082
0083     plt.ion()
0084
0085     fig, ax = GPlot.MakeFig(plt, lv, args, recurse=True)  # all volumes together
0086     if not pec is None:
0087         pec.rz_plot(ax, ipec[pmt], pec_sli )
0088     pass
0089
0090     ax.set_aspect('equal')
0091     fig.show()
0092     path = args.figpath("CombinedFig_%s" % figpath_suffix )
0093     log.info("saving to %s " % path)
0094     fig.savefig(path)
0095
0096     #axs = GPlot.MultiFig(plt, lvs, args)
0097
0098     fig, axs = GPlot.SubPlotsFig(plt, [lvs], args)
0099     if not pec is None:
0100         pec.rz_plot(axs, ipec[pmt], pec_sli )
0101     pass
0102
0103     fig.show()
0104     path = args.figpath("SplitFig_%s" % figpath_suffix )
0105     log.info("saving to %s " % path)
0106     fig.savefig(path)
0107
0108     #scribble( axs[0,2] )
0109
0110
0111 def scribble(ax):
0112     mm = 1.
0113     deg = 2.*np.pi/360.
0114
0115     m4_torus_r = 80.
0116     m4_torus_angle = 45.*deg
0117     m4_r_2 = 254./2.
0118     m4_r_1 = (m4_r_2+m4_torus_r) - m4_torus_r*np.cos(m4_torus_angle)
0119
0120     m4_h = m4_torus_r*np.sin(m4_torus_angle) + 5.0       # full height of the tube
0121
0122     m4_h/2   #    tube to centerline torus offset : so torus centerline level with bottom of tube
0123
0124     neck_z = -210.*mm+m4_h/2.
0125     torus_z = neck_z - m4_h/2
0126
0127     torus_x = m4_r_2+m4_torus_r    # radial distance to center of torus circle
0128
0129     add_line(ax, [-300,torus_z], [300,torus_z] )
0130     add_line(ax, [torus_x, -300], [torus_x, 300] )
0131
0132     e = ellipse_points( xy=[0,-5.], ex=254., ez=190., n=1000000 )
0133
0134     #ax.scatter( e[:,0], e[:,1], marker="." )
0135
0136     tc = np.array([torus_x,torus_z])
0137     tr = m4_torus_r
0138     t = circle_points( xy=tc, tr=tr , n=100 )
0139     #ax.scatter( t[:,0], t[:,1], marker="." )
0140
0141     e_inside_t = np.sqrt(np.sum(np.square(e-tc),1)) - tr < 0.  # points on the ellipse that are inside the torus circle
0142
0143     ax.scatter( e[e_inside_t][:,0], e[e_inside_t][:,1], marker="." )
0144
0145
0146
0147