qudarap/tests/QPMTTest.py

0001 #!/usr/bin/env python
0002
0003 import os, numpy as np, matplotlib as mp, textwrap
0004 from opticks.ana.fold import Fold
0005 import matplotlib.pyplot as plt
0006 SIZE = np.array([1280, 720])
0007 np.set_printoptions(edgeitems=16)
0008
0009 hc_eVnm = 1239.84198433200208455673
0010
0011 e2w_ = lambda e:hc_eVnm/e
0012 w2e_ = lambda w:hc_eVnm/w
0013
0014 PMTIDX = int(os.environ.get("PMTIDX","0")) # 0-based index into list of lpmtid
0015 SCRIPT = os.environ.get("SCRIPT", "unknown-SCRIPT")
0016
0017 class QPMTTest(object):
0018
0019     NAMES = "NNVT HAMA NNVTHiQE".split()
0020     S_NAMES = "SPMT".split()
0021
0022     def __init__(self, t):
0023         self.t = t
0024         self.init_energy_eV_domain()
0025         self.init_theta_radians_domain()
0026         self.init_mct_domain()
0027         self.init_costh_domain()
0028         self.title_prefix = "%s : %s " % ( SCRIPT, t.base )
0029
0030     def init_energy_eV_domain(self):
0031         e = self.t.qscan.energy_eV_domain
0032         #e0,e1 = 2.3, 3.3
0033         e0,e1 = 1.55, 4.3
0034         w0,w1 = e2w_(e0), e2w_(e1)
0035         ese = np.logical_and( e >= e0, e <= e1 )
0036
0037         self.e0 = e0
0038         self.e1 = e1
0039         self.w0 = w0
0040         self.w1 = w1
0041         self.ese = ese
0042
0043     def init_theta_radians_domain(self):
0044         h = self.t.qscan.theta_radians_domain
0045         h0 = h[0]
0046         h1 = h[-1]
0047         hse = np.logical_and( h >= h0, h <= h1 )
0048
0049         self.h0 = h0
0050         self.h1 = h1
0051         self.hse = hse
0052
0053     def init_mct_domain(self):
0054         mct = self.t.qscan.mct_domain
0055         self.mct = mct
0056
0057     def init_costh_domain(self):
0058         c = self.t.qscan.costh_domain
0059         c0 = c[-1]   ## flip the order as reversed
0060         c1 = c[0]
0061         cse = np.logical_and( c >= c0, c <= c1 )
0062
0063         self.c0 = c0
0064         self.c1 = c1
0065         self.cse = cse
0066
0067
0068
0069     def present_qeshape(self):
0070         t = self.t
0071         se = self.ese
0072         d = t.qscan.energy_eV_domain
0073         a = t.qscan.pmtcat_qeshape
0074         if a is None: return
0075
0076         prop_ni = t.qpmt.qeshape[:,-1,-1].view(np.int32)  ## last values from input prop arrays
0077
0078         v0,v1 = 0.0,0.38
0079
0080         assert len(a.shape) == 2, interp.shape
0081
0082         ni = a.shape[0]  # pmtcat
0083         nj = a.shape[1]  # energy
0084
0085         title = "%s : qeshape GPU interpolation lines and values " % self.title_prefix
0086
0087         fig, axs = mp.pyplot.subplots(1, ni, figsize=SIZE/100.)
0088         fig.suptitle(title)
0089
0090         for i in range(ni):
0091             ax = axs[i]
0092             ax.set_ylim( v0, v1 )
0093             v = a[i]
0094             name = self.NAMES[i]
0095             label = "%s qeshape" % name
0096             ax.set_xlabel("energy [eV]")
0097             ax.plot( d[se], v[se], label=label )
0098             ax.legend(loc=os.environ.get("LOC", "upper left")) # upper/center/lower right/left
0099
0100             p_e = t.qpmt.qeshape[i,:prop_ni[i],0]
0101             p_v = t.qpmt.qeshape[i,:prop_ni[i],1]
0102             p_s = np.logical_and( p_e >= self.e0, p_e <= self.e1 )
0103
0104             ax.scatter( p_e[p_s], p_v[p_s] )
0105         pass
0106         fig.show()
0107
0108
0109     def present_s_qeshape(self):
0110         t = self.t
0111         se = self.ese
0112         d = t.qscan.energy_eV_domain
0113         a = t.qscan.pmtcat_s_qeshape
0114         if a is None: return
0115
0116         prop_ni = t.qpmt.s_qeshape[:,-1,-1].view(np.int32)  ## last values from input prop arrays
0117
0118         v0,v1 = 0.0,0.30
0119
0120         assert len(a.shape) == 2, interp.shape
0121
0122         ni = a.shape[0]  # pmtcat
0123         nj = a.shape[1]  # energy
0124
0125         assert(ni == 1)
0126
0127         title = "%s : s_qeshape GPU interpolation lines and values " % self.title_prefix
0128
0129         fig, axs = mp.pyplot.subplots(1, ni, figsize=SIZE/100.)
0130         fig.suptitle(title)
0131
0132         for i in range(ni):
0133             ax = axs[i] if ni > 1 else axs
0134             ax.set_ylim( v0, v1 )
0135             v = a[i]
0136             name = self.S_NAMES[i]
0137             label = "%s s_qeshape" % name
0138             ax.set_xlabel("energy [eV]")
0139             ax.plot( d[se], v[se], label=label )
0140             ax.legend(loc=os.environ.get("LOC", "upper left")) # upper/center/lower right/left
0141
0142             p_e = t.qpmt.s_qeshape[i,:prop_ni[i],0]
0143             p_v = t.qpmt.s_qeshape[i,:prop_ni[i],1]
0144             p_s = np.logical_and( p_e >= self.e0, p_e <= self.e1 )
0145
0146             ax.scatter( p_e[p_s], p_v[p_s] )
0147         pass
0148         fig.show()
0149
0150
0151
0152
0153
0154
0155     def present_cetheta(self):
0156         t = self.t
0157         se = self.hse
0158         d = t.qscan.theta_radians_domain
0159         a = t.qscan.pmtcat_cetheta
0160         if a is None: return
0161
0162         prop_ni = t.qpmt.cetheta[:,-1,-1].view(np.int32)  ## last values from input prop arrays
0163
0164         v0,v1 = 0.0,1.1
0165
0166         assert len(a.shape) == 2, interp.shape
0167
0168         ni = a.shape[0]  # pmtcat
0169         nj = a.shape[1]  # theta
0170
0171         title = "%s : cetheta GPU interpolation lines and values " % self.title_prefix
0172
0173         fig, axs = mp.pyplot.subplots(1, ni, figsize=SIZE/100.)
0174         fig.suptitle(title)
0175
0176         for i in range(ni):
0177             ax = axs[i]
0178             ax.set_ylim( v0, v1 )
0179             v = a[i]
0180             name = self.NAMES[i]
0181             label = "%s cetheta" % name
0182             ax.set_xlabel("theta [radians]")
0183             ax.plot( d[se], v[se], label=label )
0184             ax.legend(loc=os.environ.get("LOC", "upper left")) # upper/center/lower right/left
0185
0186             ## input (domain,value) pairs used by the interpolation
0187             p_d = t.qpmt.cetheta[i,:prop_ni[i],0]
0188             p_v = t.qpmt.cetheta[i,:prop_ni[i],1]
0189             p_s = np.logical_and( p_d >= self.h0, p_d <= self.h1 )
0190
0191             ax.scatter( p_d[p_s], p_v[p_s] )
0192         pass
0193         fig.show()
0194
0195
0196     def present_cecosth(self):
0197         t = self.t
0198         se = self.cse
0199         d = t.qscan.costh_domain
0200         a = t.qscan.pmtcat_cecosth
0201         if a is None: return
0202
0203         prop_ni = t.qpmt.cecosth[:,-1,-1].view(np.int32)  ## last values from input prop arrays
0204
0205         v0,v1 = 0.0,1.1
0206
0207         assert len(a.shape) == 2, interp.shape
0208
0209         ni = a.shape[0]  # pmtcat
0210         nj = a.shape[1]  # theta
0211
0212         title = "%s : cecosth GPU interpolation lines and values " % self.title_prefix
0213
0214         fig, axs = mp.pyplot.subplots(1, ni, figsize=SIZE/100.)
0215         fig.suptitle(title)
0216
0217         for i in range(ni):
0218             ax = axs[i]
0219             ax.set_ylim( v0, v1 )
0220             v = a[i]
0221             name = self.NAMES[i]
0222             label = "%s cecosth" % name
0223             ax.set_xlabel("cosine_theta")
0224             ax.plot( d[se], v[se], label=label )
0225             ax.legend(loc=os.environ.get("LOC", "upper left")) # upper/center/lower right/left
0226
0227             ## input (domain,value) pairs used by the interpolation
0228             p_d = t.qpmt.cecosth[i,:prop_ni[i],0]
0229             p_v = t.qpmt.cecosth[i,:prop_ni[i],1]
0230             p_s = np.logical_and( p_d >= self.c0, p_d <= self.c1 )
0231
0232             ax.scatter( p_d[p_s], p_v[p_s] )
0233         pass
0234         fig.show()
0235
0236
0237     def present_rindex(self):
0238         t = self.t
0239
0240         a = t.qscan.pmtcat_rindex
0241         if a is None: return
0242         assert len(a.shape) == 4, a.shape
0243
0244         se = self.se
0245         e = t.qscan.energy_eV_domain
0246
0247         prop_ni = t.qpmt.rindex[:,-1,-1].view(np.int32)
0248         v0,v1 = -0.1,3.2
0249
0250         ni = a.shape[0]  # pmtcat
0251         nj = a.shape[1]  # layers
0252         nk = a.shape[2]  # props
0253
0254         title = "%s : PMT layer refractive index interpolations on GPU  " % self.title_prefix
0255
0256         fig, axs = mp.pyplot.subplots(1, ni, figsize=SIZE/100.)
0257         fig.suptitle(title)
0258
0259         for i in range(ni):
0260
0261             ax = axs[i]
0262             ax.set_ylim( v0, v1 )
0263
0264             name = self.NAMES[i]
0265             #ax.set_title(name)
0266             ax.set_xlabel("energy [eV]")
0267
0268             sax = ax.secondary_xaxis('top', functions=(e2w_, w2e_))
0269             sax.set_xlabel('%s   wavelength [nm]' % name)
0270             # secondary_xaxis w2e_ : RuntimeWarning: divide by zero encountered in true_divide
0271
0272             for j in range(nj):
0273                 if j in [0,3]: continue   # skip layers 0,3 Pyrex,Vacuum
0274                 for k in range(nk):
0275                     v = a[i,j,k]
0276                     iprop = i*nj*nk + j*nk + k
0277
0278                     label = "L%d %sINDEX" % ( j, "R" if k == 0 else "K" )
0279
0280                     ax.plot( e[se], v[se], label=label )
0281
0282                     p_ni = prop_ni[iprop]
0283                     p_e = t.qpmt.rindex[iprop,:p_ni,0]
0284                     p_v = t.qpmt.rindex[iprop,:p_ni,1]
0285
0286                     p_s = np.logical_and( p_e >= self.e0, p_e <= self.e1 )
0287                     ax.scatter( p_e[p_s], p_v[p_s] )
0288                 pass
0289             pass
0290             ax.legend(loc=os.environ.get("LOC", "lower right")) # upper/center/lower right/left
0291         pass
0292         fig.show()
0293
0294
0295     def present_atqc(self):
0296         """
0297         In [4]: t.qscan.atqc.shape
0298         Out[4]: (9, 900, 4)
0299
0300         In [5]: t.qscan.lpmtid
0301         Out[5]: array([    0,    10,    55,    98,   100,   137,  1000, 10000, 17611], dtype=int32)
0302
0303         t.qscan.mct_domain.shape
0304         Out[8]: (900,)
0305
0306         900 is scan over mct : minus_cos_theta of angle of incidence   (NOT landing position)
0307
0308         atqc[:,:,3]
0309               behaves as expected, fixed across the scan at high values ~0.91 - 0.97 depending on pmtid
0310
0311         """
0312         pass
0313         t = self.t
0314         mct = t.qscan.mct_domain
0315         atqc = t.qscan.atqc
0316
0317
0318
0319     def present_art(self):
0320         """
0321         In [2]: t.qscan.art.shape
0322         Out[2]: (9, 181, 4, 4)
0323         """
0324
0325         t = self.t
0326         lpmtid = t.qscan.lpmtid[PMTIDX]   # pick lpmtid by PMTIDX
0327
0328         all_art = t.qscan.art
0329         if all_art is None:
0330             print("present_art : ABORT t.qscan.art is None ")
0331             return
0332         pass
0333         art = all_art[PMTIDX]
0334         mct = t.qscan.mct_domain
0335
0336         consistent = len(art) == len(mct)
0337
0338         if not consistent:
0339             log.error("present_lpmtid_ART : INCONSISTENT : art.shape %s mct.shape %s " %
0340                      (str(art.shape), str(mct.shape)) )
0341             return
0342         pass
0343         assert consistent
0344
0345         As   = art[...,0,0]
0346         Ap   = art[...,0,1]
0347         Aa   = art[...,0,2]
0348         A_   = art[...,0,3]
0349
0350         Rs   = art[...,1,0]
0351         Rp   = art[...,1,1]
0352         Ra   = art[...,1,2]
0353         R_   = art[...,1,3]
0354
0355         Ts   = art[...,2,0]
0356         Tp   = art[...,2,1]
0357         Ta   = art[...,2,2]
0358         T_   = art[...,2,3]
0359
0360         SF    = art[...,3,0]
0361         wl    = art[...,3,1]
0362         ARTa  = art[...,3,2]
0363         mct   = art[...,3,3]
0364
0365
0366         opt = os.environ.get("OPT", "A_,R_,T_,As,Rs,Ts,Ap,Rp,Tp,Aa,Ra,Ta")
0367         title = "%s : PMTIDX %d lpmtid %d OPT %s " % (t.base, PMTIDX, lpmtid, opt)
0368         fig, ax = plt.subplots(1, figsize=SIZE/100.)
0369         fig.suptitle(title)
0370
0371         if "As" in opt:ax.plot(  mct, As, label="As" )
0372         if "Ap" in opt:ax.plot(  mct, Ap, label="Ap" )
0373         if "Aa" in opt:ax.plot(  mct, Aa, label="Aa" )
0374         if "A_" in opt:ax.plot(  mct, A_, label="A_" )
0375
0376         if "Rs" in opt:ax.plot(  mct, Rs, label="Rs" )
0377         if "Rp" in opt:ax.plot(  mct, Rp, label="Rp" )
0378         if "Ra" in opt:ax.plot(  mct, Ra, label="Ra" )
0379         if "R_" in opt:ax.plot(  mct, R_, label="R_" )
0380
0381         if "Ts" in opt:ax.plot(  mct, Ts, label="Ts" )
0382         if "Tp" in opt:ax.plot(  mct, Tp, label="Tp" )
0383         if "Ta" in opt:ax.plot(  mct, Ta, label="Ta" )
0384         if "T_" in opt:ax.plot(  mct, T_, label="T_" )
0385
0386         if "SF" in opt:ax.plot(  mct, SF, label="SF")
0387         if "wl" in opt:ax.plot(  mct, wl, label="wl" )
0388         if "ARTa" in opt:ax.plot(  mct, ARTa, label="ARTa" )
0389         if "mct" in opt:ax.plot(  mct, mct, label="mct" )
0390
0391
0392         ax.legend()
0393         fig.show()
0394
0395
0396
0397     def check_lpmtcat(self):
0398         t = self.t
0399         p = t.qpmt  # QPMT members
0400         s = t.qscan # QPMTTest scan results
0401
0402         expect_lpmtcat = p.src_lcqs[s.lpmtidx,0]
0403         lpmtcat = s.spec[:,:,0,3].astype(np.int32)
0404         assert( np.all( lpmtcat[:,0] == expect_lpmtcat ) )
0405
0406         lpmtid = s.lpmtid
0407
0408         ## across the mct scan the identity values all the same, so use np.max
0409         lpmtid_lpmtcat  = np.max(s.spec[:,:,0,3].astype(np.int32), axis=1)
0410         lpmtid_qe_scale = np.max(s.spec[:,:,1,3], axis=1)
0411         lpmtid_qe_shape = np.max(s.spec[:,:,2,3], axis=1)
0412         lpmtid_qe       = np.max(s.spec[:,:,3,3], axis=1)
0413         ## these were formerly lpmtid_stackspec now just spec
0414
0415         expr = "np.c_[lpmtid,lpmtid_lpmtcat,lpmtid_qe_scale,lpmtid_qe_shape,lpmtid_qe]"
0416         lpmtid_tab = eval(expr)
0417         print("lpmtid_tab:%s\n%s" % ( expr,  lpmtid_tab))
0418         print(" note the qe_shape factor depends only on lpmtcat, the others have lpmtid dependency ")
0419         print(" also note the qe_shape factor for lpmtcat 0:NNVT and 2:NNVT_HiQE are the same, diff from 1:HAMA  ")
0420
0421
0422 class SPMTID(object):
0423     """
0424     t.qscan.s_qescale.shape
0425     t.qscan.s_qescale
0426
0427     p.s_qescale[:,0].shape
0428     p.s_qescale[:,0]
0429
0430     t.qscan.spmtid.shape
0431     t.qscan.spmtid
0432
0433     np.all( t.qscan.s_qescale == p.s_qescale[:,0][t.qscan.spmtid - 20000] )
0434
0435     """
0436     @classmethod
0437     def EXPR(cls):
0438         return list(map(str.strip,textwrap.dedent(cls.__doc__).split("\n")))
0439
0440     def __init__(self, t, symbol="sp"):
0441         self.t = t
0442         self.symbol = symbol
0443
0444     def __repr__(self):
0445         t = self.t
0446         p = t.qpmt  # QPMT members
0447         s = t.qscan # QPMTTest scan results
0448
0449         lines = []
0450         lines.append("[%s" % self.symbol)
0451         for expr in self.EXPR():
0452             lines.append(expr)
0453             if expr == "" or expr[0] == "#": continue
0454             lines.append(repr(eval(expr)))
0455         pass
0456         lines.append("]%s" % self.symbol)
0457         return "\n".join(lines)
0458
0459
0460 if __name__ == '__main__':
0461     t = Fold.Load(symbol="t")
0462     print(repr(t))
0463     pt = QPMTTest(t)
0464
0465     pt.check_lpmtcat()
0466
0467     sp = SPMTID(t, symbol="sp")
0468     print(repr(sp))
0469
0470     p = t.qpmt
0471     s = t.qscan
0472
0473     #plot = "rindex"
0474     #plot = "qeshape"
0475     #plot = "cetheta"
0476     #plot = "cecosth"
0477     #plot = "art"
0478     plot = "s_qeshape"
0479
0480     PLOT = os.environ.get("PLOT", plot )
0481     if PLOT == "rindex":
0482         pt.present_rindex()
0483     elif PLOT == "qeshape":
0484         pt.present_qeshape()
0485     elif PLOT == "s_qeshape":
0486         pt.present_s_qeshape()
0487     elif PLOT == "cetheta":
0488         pt.present_cetheta()
0489     elif PLOT == "cecosth":
0490         pt.present_cecosth()
0491     elif PLOT == "art":
0492         pt.present_art()
0493     else:
0494         print("PLOT:%s not handled " % PLOT)
0495     pass
0496
0497
0498