sysrap/tests/stree_mat_test.py

0001 #!/usr/bin/env python
0002
0003 import numpy as np, textwrap
0004 from opticks.ana.fold import Fold
0005 import matplotlib.pyplot as plt
0006 SIZE = np.array([1280,720])
0007
0008
0009 def eprint(exprs):
0010     print("eprint(\"\"\"%s\"\"\")" % exprs  )
0011     for expr in list(filter(None,textwrap.dedent(exprs).split("\n"))):
0012         print(expr)
0013         if expr[0] in " #": continue
0014         print(eval(expr))
0015         print("\n")
0016     pass
0017
0018
0019
0020 if __name__ == '__main__':
0021     t = Fold.Load(symbol="t")
0022     print(repr(t))
0023
0024     a = t.gg    # old X4/GGeo NPFold
0025     b = t.st    # new U4/stree NPFold
0026
0027     # kludge to avoid obnoxious array presentation
0028     # when have very large values
0029     a.mat[np.where( a.mat == 1e9 )] = 1e6
0030     b.mat[np.where( b.mat == 1e9 )] = 1e6
0031     a.bnd[np.where( a.bnd == 1e9 )] = 1e6
0032     b.bnd[np.where( b.bnd == 1e9 )] = 1e6
0033
0034
0035     abn = np.array(a.bnd_names)
0036     aop = a.optical.reshape(-1,4,4).view(np.int32)
0037     assert len(abn) == len(aop)
0038
0039     bbn = np.array(b.bnd_names)
0040     bop = b.optical
0041     assert len(bbn) == len(bop)
0042
0043
0044     # tracking down where missing sur entry within bnd diff of 1. is coming from
0045     ab = a.bnd[:,1,0,0,0] - b.bnd[:,1,0,0,0]   # osur check
0046
0047
0048     eprint("""
0049     np.all( a.mat_names == b.mat_names )
0050     np.all( a.sur_names == b.sur_names )
0051     np.all( a.bnd_names == b.bnd_names )
0052
0053     a.mat.shape == b.mat.shape
0054     np.unique(np.where( np.abs(a.mat - b.mat) > 1e-3 )[0])
0055     a.mat_names[np.unique(np.where( np.abs(a.mat - b.mat) > 1e-3 )[0])]
0056
0057     #  RINDEX     ABSLENGTH  RAYLEIGH   REEMISSIONPROB   GROUPVEL
0058     np.max(np.abs(a.mat - b.mat), axis=2).reshape(-1,8)   # max deviation across wavelength domain
0059     np.c_[np.arange(a.mat_names.size),a.mat_names]
0060
0061     np.all( a.icdf == b.icdf )
0062     np.all( a.optical == b.optical )
0063     np.where( np.abs(a.sur - b.sur) > 1e-6 )
0064
0065     np.unique(np.where( np.abs( a.bnd - b.bnd ) > 1e-4 )[1])  ## only omat and imat
0066
0067     # payload group 0 has small discreps for ABSLENGTH and RAYLEIGH
0068     np.unique(np.where( np.abs( a.mat[:,0,:,:] - b.mat[:,0,:,:]) > 1e-4 )[2])
0069
0070     # payload group 1 has discreps in GROUPVEL
0071     np.unique(np.where( np.abs( a.mat[:,1,:,:] - b.mat[:,1,:,:]) > 1e-4 )[2])
0072
0073     # max absolute discrepancies in ABSLENGTH for each material
0074     np.c_[a.mat_names,  np.max( np.abs( a.mat[:,0,:,1]-b.mat[:,0,:,1]), axis=1 ) ]
0075
0076     # max relative difference from 1. in ABSLENGTH for each material
0077     np.c_[a.mat_names, np.max( np.abs( (a.mat[:,0,:,1]/b.mat[:,0,:,1]) - 1. ), axis=1 ) ]
0078
0079
0080     # max absolute discrepancies in RAYLEIGH for each material
0081     np.c_[a.mat_names,  np.max( np.abs( a.mat[:,0,:,2]-b.mat[:,0,:,2]), axis=1 ) ]
0082
0083     # max relative difference from 1. in RAYLEIGH for each material
0084     np.c_[a.mat_names, np.max( np.abs( (a.mat[:,0,:,2]/b.mat[:,0,:,2]) - 1. ), axis=1 ) ]
0085
0086
0087     # max absolute discrepancies in GROUPVEL for each material
0088     np.c_[a.mat_names,  np.max( np.abs( a.mat[:,1,:,0]-b.mat[:,1,:,0]), axis=1 ) ]
0089
0090     # max relative difference from 1. in GROUPVEL for each material
0091     np.c_[a.mat_names,  np.max( np.abs( a.mat[:,1,:,0]/b.mat[:,1,:,0] - 1.), axis=1 ) ]
0092
0093     # substantial diffs in GROUPVEL : comparing plots
0094     # seems to indicate a calculation difference wrt bin handling
0095     # new way is coming more directly from Geant4 so simply adopt the new
0096
0097
0098     """
0099     )
0100
0101
0102 if 0:
0103     BASE = "$HOME/.opticks/GEOM/$GEOM/CSGFoundry/SSim"
0104
0105     ## t.mat[np.where( t.mat == 300. )] = 299.792458  # GROUPVEL default kludge
0106
0107     rayleigh_Water_idx = np.where( np.array(t.rayleigh_names)  == "Water" )[0][0]
0108     rayleigh_Water_ = t.rayleigh[rayleigh_Water_idx]
0109     rayleigh_Water_ni = rayleigh_Water_[-1,-1].view(np.int64)
0110     rayleigh_Water = rayleigh_Water_[:rayleigh_Water_ni]
0111
0112     #wl = np.linspace(60.,820.,761)
0113     wl = t.wavelength
0114     en = t.energy
0115
0116     doms = {"wl":wl, "en":en }
0117     DOM = os.environ.get("DOM", "en")
0118     dom = doms[DOM[:2]]
0119
0120     if DOM == "ensel4":
0121         sel = np.where( dom < 4 )
0122     else:
0123         sel = slice(None)
0124     pass
0125
0126     plot = "RINDEX"
0127     #plot = "GROUPVEL"
0128     #plot = "RAYLEIGH"
0129     #plot = "ABSLENGTH"
0130
0131     PLOT = os.environ.get("PLOT", plot)
0132     qwns="mat oldmat".split()
0133
0134     MM = [4,11,14,17,18,19]
0135
0136     prop = {}
0137     prop["RINDEX"] = (0,0)
0138     prop["ABSLENGTH"] = (0,1)
0139     prop["RAYLEIGH"] = (0,2)
0140     prop["REEMISSIONPROB"] = (0,3)
0141     prop["GROUPVEL"] = (1,0)
0142
0143     for M in MM:
0144         MAT = t.mat_names[M]
0145
0146         title = "PLOT:%s DOM:%s M:%d MAT:%s " % (PLOT,DOM, M, MAT)
0147         print(title)
0148
0149         fig, ax = plt.subplots(1, figsize=SIZE/100.)
0150         fig.suptitle(title)
0151
0152         if PLOT.split("_")[0] == "DIFF":
0153             a = getattr(t, qwns[0], None)
0154             b = getattr(t, qwns[1], None)
0155             p = prop.get(PLOT.split("_")[1], None)
0156             aq = a[M,p[0],:,p[1]]
0157             bq = b[M,p[0],:,p[1]]
0158             abq = aq - bq
0159
0160             label = "%s %s" % (PLOT, "%s-%s" % (qwns[0],qwns[1]) )
0161             ax.plot( dom[sel], abq[sel] , label=label )
0162
0163             if PLOT == "DIFF_RAYLEIGH" and DOM[:2] == "en" and MAT == "Water":
0164
0165                 rayleigh_Water_en = rayleigh_Water[:,0]*1e6
0166                 if DOM == "ensel4":
0167                     rayleigh_Water_ensel = np.where(rayleigh_Water_en < 4 )
0168                 else:
0169                     rayleigh_Water_ensel = slice(None)
0170                 pass
0171                 for xc in rayleigh_Water_en[rayleigh_Water_ensel]:
0172                     ax.axvline(x=xc)
0173                 pass
0174             pass
0175         else:
0176             for qwn in qwns:
0177                 a = getattr(t, qwn, None)
0178                 assert not a is None
0179                 p = prop.get(PLOT, None)
0180                 assert not p is None
0181                 aq = a[M,p[0],:,p[1]]
0182                 label = "%s %s" % (PLOT, qwn)
0183
0184                 ax.plot( dom[sel], aq[sel], label=label )
0185
0186                 REL = "stree/material/%(MAT)s/%(PLOT)s.npy" % locals()
0187                 path = os.path.expandvars(os.path.join(BASE,REL))
0188
0189                 if DOM[:2] == "en" and os.path.exists(path):
0190                     meas = np.load(path)
0191                     meas_en = meas[:,0]*1e6
0192                     meas_va = meas[:,1]
0193                     if DOM == "ensel4":
0194                         meas_sel = np.where(meas_en < 4 )
0195                     else:
0196                         meas_sel = slice(None)
0197                     pass
0198                     ax.scatter( meas_en[meas_sel], meas_va[meas_sel], label="meas" )
0199                 pass
0200
0201                 RID = "stree/material/%(MAT)s/RINDEX.npy" % locals()
0202                 ri_path = os.path.expandvars(os.path.join(BASE,RID))
0203
0204                 if DOM[:2] == "en" and os.path.exists(ri_path):
0205                     ri = np.load(ri_path)
0206                     ri_en = ri[:,0]*1e6
0207                     ri_va = ri[:,1]
0208                     if DOM == "ensel4":
0209                         ri_ensel = np.where(ri_en < 4)
0210                     else:
0211                         ri_ensel = slice(None)
0212                     pass
0213                     for xc in ri_en[ri_ensel]:
0214                         ax.axvline(x=xc)
0215                     pass
0216                 pass
0217                 #if PLOT == "RAYLEIGH" and DOM[:2] == "en" and MAT == "Water":
0218                 #    ax.scatter( rayleigh_Water[:,0]*1e6, rayleigh_Water[:,1], label="rayleigh_Water")
0219                 #pass
0220             pass
0221         pass
0222         ax.legend()
0223         fig.show()
0224     pass
0225
0226 if 0:
0227     o = Fold.Load("/tmp/SBnd_test", symbol="o")
0228     print(repr(o))
0229
0230     s = Fold.Load("/tmp/blyth/opticks/U4TreeCreateTest/stree", symbol="s")
0231     print(repr(s))
0232
0233     Vacuum_kludge([t,o,s])
0234
0235
0236     exprs = """
0237     np.all( t.old_optical == o.optical )
0238     np.all( np.array( t.mat_names ) == np.array( o.mat_names ) )
0239     """
0240
0241     for expr in list(filter(None,textwrap.dedent(exprs).split("\n"))):
0242         print(expr)
0243         print(eval(expr))
0244     pass
0245
0246     assert len(t.mat) == len(o.mat)
0247
0248     print("ij")
0249     for i in range(len(o.mat)):
0250         for j in range(2):
0251             expr= "np.all( o.mat[%(i)d,%(j)d] == s.mat[%(i)d,%(j)d] )" % locals()
0252             print(" %s : %s " % (expr, eval(expr)))
0253         pass
0254     pass
0255
0256     print("ijk")
0257     for i in range(len(o.mat)):
0258         tname = t.mat_names[i]
0259         oname = o.mat_names[i]
0260         assert( tname == oname )
0261         print( "\n i : %(i)d  %(tname)s " % locals() )
0262         for j in range(2):
0263             print( " j : %(j)d " % locals() )
0264             for k in range(4):
0265                 expr= "len(np.where( np.abs( o.mat[%(i)d,%(j)d,:,%(k)d] - s.mat[%(i)d,%(j)d,:,%(k)d] ) > 1e-4)[0])" % locals()
0266                 print(" %s : %s " % (expr, eval(expr)))
0267             pass
0268         pass
0269     pass
0270 pass
0271