qudarap/tests/rayleigh_scatter_align.py

0001 #!/usr/bin/env python
0002
0003 import os, numpy as np
0004 from opticks.ana.fold import Fold
0005 from opticks.ana.pvplt import pvplt_simple, pvplt_photon
0006 import pyvista as pv
0007 import matplotlib.pyplot as plt
0008
0009 GUI = not "NOGUI" in os.environ
0010 SIZE = np.array([1280, 720])
0011
0012
0013 FOLD = os.environ["FOLD"]
0014 TEST = os.environ["TEST"]
0015
0016 G4_FOLD = "/tmp/G4OpRayleighTest"
0017
0018
0019
0020 if __name__ == '__main__':
0021
0022     #t = Fold.Load(FOLD)
0023     t = Fold.Load(G4_FOLD)
0024
0025     print(t.p[:,:3])
0026     print(t.p[:,3].view(np.uint32))
0027
0028
0029     mom = t.p[:,1,:3]
0030     pol = t.p[:,2,:3]
0031
0032     mom0 = t.p0[0,1,:3]
0033     pol0 = t.p0[0,2,:3]
0034
0035     ct = np.sum( pol*pol0, axis=1 )
0036
0037
0038     transverse = np.sum( mom*pol, axis=1 )
0039     tr_min, tr_max, tr_abs = transverse.min(), transverse.max(), np.abs(transverse).max()
0040     print( " tr_min : %s  tr_max : %s tr_abs : %s " % ( tr_min, tr_max, tr_abs ))
0041     assert tr_abs < 1e-5
0042
0043
0044     mom_norm = np.abs(np.sum(mom*mom, axis=1 ) - 1.).max()
0045     pol_norm = np.abs(np.sum(pol*pol, axis=1 ) - 1.).max()
0046
0047     print(" mom_norm : %s pol_norm %s " % (mom_norm, pol_norm))
0048     assert mom_norm < 1e-5
0049     assert pol_norm < 1e-5
0050
0051     pl = pv.Plotter(window_size=SIZE*2 )  # retina 2x ?
0052     pl.show_grid()
0053
0054     #pvplt_simple(pl, mom, "%s.mom" % TEST )
0055
0056     pvplt_simple(pl, pol[:100000], "%s.pol" % TEST )
0057     pvplt_photon(pl,  t.p0 )
0058
0059     cp = pl.show() if GUI else None
0060
0061
0062     fig, ax = plt.subplots(1)
0063     ax.hist( ct*ct, bins=100 )
0064     fig.show()
0065
0066