qudarap/tests/propagate_at_multifilm.py

0001 #!/usr/bin/env python
0002 """
0003 propagate_at_multifilm.py
0004 =================================
0005
0006 As everything is happening at the origin it would be
0007 impossible to visualize everything on top of each other.
0008 So represent the incoming hemisphere of photons with
0009 points and vectors on the unit hemi-sphere with mom
0010 directions all pointing at the origin and polarization
0011 vector either within the plane of incidence (P-polarized)
0012 or perpendicular to the plane of incidence (S-polarized).
0013
0014
0015
0016                    .
0017             .              .
0018         .                     .
0019
0020      .                          .
0021
0022     -------------0---------------
0023
0024
0025
0026 Notice with S-polarized that the polarization vectors Z-component is zero
0027
0028 ::
0029
0030    EYE=-1,-1,1 LOOK=0,0,0.5 PARA=1 ./QSimTest.sh ana
0031
0032
0033 """
0034 import os, numpy as np
0035 from opticks.ana.fold import Fold
0036 from opticks.ana.pvplt import *
0037 import pyvista as pv
0038 import matplotlib.pyplot as plt
0039
0040 #FOLD = os.environ["FOLD"]
0041 #TEST = os.environ["TEST"]
0042 GUI = not "NOGUI" in os.environ
0043
0044 '''
0045 In [13]: prd=np.load("prd.npy")
0046
0047 In [14]: prd.shape
0048 Out[14]: (1, 2, 4)
0049 '''
0050
0051 class MultiFilmPropagate():
0052     def __init__(self, src , mutate_src , pmttype = 0 , boundary = 0):
0053         self.mutate_src_path = mutate_src;
0054         self.t0 = Fold.Load(src)
0055         self.t1 = Fold.Load(mutate_src)
0056         self.src_photon = self.t0.p
0057         self.mutate_src_photon = self.t1.p
0058         self.prd = self.t1.prd
0059         self.pmttype = pmttype
0060         self.boundary = boundary
0061
0062         pmtStringName={}
0063         pmtStringName[0] = "kPMT_NNVT"
0064         pmtStringName[1] = "kPMT_Hama"
0065         pmtStringName[2] = "kPMT_NNVT_HiQE"
0066
0067         self.pmtStringName = pmtStringName
0068
0069         direction = {}
0070         direction[0] = "glass_to_vacuum"
0071         direction[1] = "vacuum_to_glass"
0072
0073         self.direction = direction
0074         flag = {}
0075         flag[1024] = "boundary_reflect"
0076         flag[2048] = "boundary_refract"
0077         self.flag = flag
0078
0079         print(" src_photon.shape = {}  {}".format(self.src_photon.shape,self.src_photon[:100,:,:]))
0080         print(" mutate_src_photon = {}".format(self.mutate_src_photon.shape))
0081         print(" prd.shape = {} ".format(self.prd.shape))
0082
0083     def draw_fix_wv(self):
0084         total_c1_ = self.src_photon[:,1,2] #mom.z
0085         total_c1 = np.absolute(total_c1_)
0086         reflect_c1_ = self.src_photon[self.mutate_src_photon[:,3,-1].view(np.uint32) == 1024, 1, 2]
0087         reflect_c1 = np.absolute(reflect_c1_)
0088         transmit_c1_ = self.src_photon[self.mutate_src_photon[:,3,-1].view(np.uint32) == 2048, 1, 2]
0089         transmit_c1 = np.absolute(transmit_c1_)
0090         print("total_c1 = {} , reflect_c1 = {}".format(total_c1.shape,reflect_c1.shape))
0091         print("total_c1 = {} , reflect_c1 = {}".format(total_c1.shape,reflect_c1.shape))
0092
0093         num_bins = 200
0094         bins = np.linspace(0.0, 1.0, num_bins+1)
0095
0096         count_total,   bin_array,  patch = plt.hist(total_c1 , bins)
0097         count_reflect , bin_array, patch = plt.hist( reflect_c1 , bins)
0098         count_transmit , bin_array, patch = plt.hist( transmit_c1 , bins)
0099         plt.cla()
0100
0101         reflect_ratio = count_reflect/count_total
0102         transmit_ratio = count_transmit/count_total
0103
0104         print("count_reflect.shape = {} bin ={}".format(count_reflect.shape,bins.shape))
0105         plt.scatter(bins[:-1],reflect_ratio, s=0.5,label="reflect")
0106         plt.scatter(bins[:-1],transmit_ratio,s=0.5,label="Transmit")
0107         plt.legend()
0108         plt.xlabel("cos_theta")
0109         plt.ylabel("probability")
0110         title = self.mutate_src_path +"\n" +self.pmtStringName[self.pmttype]+ "/" + self.direction[self.boundary]
0111         plt.title(title)
0112         plt.show()
0113
0114
0115         #save data to file which will be used by other programm
0116         data=np.zeros((3,num_bins),dtype = float )
0117         data[0,:] = bins[:-1]
0118         data[1,:] = reflect_ratio
0119         data[2,:] = transmit_ratio
0120
0121         tmp_file_name = self.mutate_src_path+"pmttype{}_boundary{}_R_T.npy".format(self.pmttype,self.boundary)
0122         np.save(tmp_file_name,data)
0123         print("the data save in {}".format(tmp_file_name))
0124
0125     def draw_with_pv(self, flag):
0126
0127         '''
0128         1024 BF reflect
0129         2048 BR
0130         '''
0131
0132         p = self.mutate_src_photon
0133         p = p[p[:,3,-1].view(np.uint32) == flag , :, :]
0134
0135         num_point = 1000
0136
0137         prd = self.prd
0138         lim = slice(0,num_point)
0139
0140         #print( " TEST : %s " % TEST)
0141         #print( " FOLD : %s " % FOLD)
0142         print( "p.shape %s " % str(p.shape) )
0143         print( "prd.shape %s " % str(prd.shape) )
0144         print(" using lim for plotting %s " % lim )
0145
0146         mom = p[:,1,:3]   # hemisphere of photons all directed at origin
0147         pol = p[:,2,:3]   # S or P polarized
0148         pos = mom          # illustrative choice of position on unit hemisphere
0149
0150         normal = prd[:,0,:3]  # saved from qprd
0151         #point =  prd[:,1,:3]  # not really position but its all zeros... so will do
0152         point = np.array( [0,0,0], dtype=np.float32 )
0153
0154         print("mom\n", mom)
0155         print("pol\n", pol)
0156         print("pos\n", pos)
0157
0158         label = "pvplt_polarized"+"_"+self.flag[flag]+"_"+self.mutate_src_path
0159         pl = pvplt_plotter(label=label)
0160
0161
0162         pvplt_viewpoint( pl )
0163         pvplt_polarized( pl, pos[lim], mom[lim], pol[lim] )
0164
0165
0166         pos = np.zeros((num_point,3),dtype=float)
0167         pvplt_lines(     pl, pos[lim], mom[lim] )
0168
0169
0170         pvplt_arrows( pl, point, normal )
0171         FOLD=self.mutate_src_path
0172
0173         label = "pvplt_polarized"+"_"+self.flag[flag]
0174         outpath = os.path.join(FOLD, "figs/%s.png" % label )
0175         outdir = os.path.dirname(outpath)
0176         if not os.path.isdir(outdir):
0177             os.makedirs(outdir)
0178         pass
0179
0180         print(" outpath: %s " % outpath )
0181         cp = pl.show(screenshot=outpath) if GUI else None
0182
0183 if __name__ == '__main__':
0184     src = os.path.expandvars("/tmp/${USER}/opticks/QSimTest/hemisphere_s_polarized/")
0185     #src = os.path.expandvars("/tmp/${USER}/opticks/QSimTest/hemisphere_p_polarized/")
0186     #src = os.path.expandvars("/tmp/${USER}/opticks/QSimTest/hemisphere_x_polarized/")
0187     mutate_src = os.path.expandvars("/tmp/${USER}/opticks/QSimTest/propagate_at_multifilm_s_polarized/")
0188     #mutate_src = os.path.expandvars("/tmp/${USER}/opticks/QSimTest/propagate_at_multifilm_p_polarized/")
0189     #mutate_src =os.path.expandvars("/tmp/${USER}/opticks/QSimTest/propagate_at_multifilm_x_polarized/")
0190     a=MultiFilmPropagate(src,mutate_src, pmttype = 0,boundary = 0)
0191     #a.draw_fix_wv()
0192     a.draw_with_pv(flag=1024)