g4cx/tests/G4CXTest_raindrop.py

0001 #!/usr/bin/env python
0002 """
0003 G4CXTest_raindrop.py
0004 ======================
0005
0006
0007 """
0008 import os, logging, textwrap, numpy as np
0009 from opticks.ana.fold import Fold
0010 from opticks.sysrap.sevt import SEvt, SAB
0011 from opticks.ana.nbase import chi2_pvalue
0012
0013 MODE = int(os.environ.get("MODE","3"))
0014 PICK = os.environ.get("PICK","B")
0015
0016 if MODE in [2,3]:
0017     try:
0018         import opticks.ana.pvplt as pvp
0019     except ImportError:
0020         pvp = None
0021     pass
0022 else:
0023     pvp = None
0024 pass
0025
0026
0027 def eprint(expr, l, g ):
0028     print(expr)
0029     try:
0030        val = eval(expr,l,g)
0031     except AttributeError:
0032        val = "eprint:AttributeError"
0033     pass
0034     print(val)
0035 pass
0036
0037
0038 if __name__ == '__main__':
0039     logging.basicConfig(level=logging.INFO)
0040
0041     a = SEvt.Load("$AFOLD", symbol="a")
0042     print(repr(a))
0043     b = SEvt.Load("$BFOLD", symbol="b")
0044     print(repr(b))
0045     ee = {'A':a, 'B':b}
0046
0047     print("[--- ab = SAB(a,b) ----")
0048
0049     if not a is None and not b is None:
0050         ab = SAB(a,b)
0051     else:
0052         ab = None
0053     pass
0054     print("]--- ab = SAB(a,b) ----")
0055
0056     print("[----- repr(ab) ")
0057     print(repr(ab))
0058     print("]----- repr(ab) ")
0059
0060
0061     EXPR_ = r"""
0062     np.c_[np.unique(a.q, return_counts=True)]
0063     np.c_[np.unique(b.q, return_counts=True)]
0064     """
0065     EXPR = list(filter(None,textwrap.dedent(EXPR_).split("\n")))
0066     for expr in EXPR:eprint(expr, locals(), globals() )
0067
0068
0069     #select = "WL"            ## select on wavelength to avoid unfilled zeros
0070     #select = "TO BT BT SA"   ## select on photon history
0071     select = "TO BT SA,TO BR BT SA"
0072     MSELECT = os.environ.get("SELECT", select )
0073
0074     sli = slice(0,100000)   # restrict to 1M to stay interactive
0075
0076     speed_ = lambda r,i:np.sqrt(np.sum( (r[:,i+1,0,:3]-r[:,i,0,:3])*(r[:,i+1,0,:3]-r[:,i,0,:3]), axis=1))/(r[:,i+1,0,3]-r[:,i,0,3])
0077
0078     for Q in PICK:
0079         e = ee.get(Q,None)
0080         if e is None:continue
0081
0082         for SELECT in MSELECT.split(","):
0083             context = "PICK=%s MODE=%d SELECT=\"%s\" ~/opticks/g4cx/tests/G4CXTest_raindrop.sh " % (Q, MODE, SELECT )
0084             print(context)
0085
0086             if SELECT == "WL":
0087                 sel = np.where(e.f.record[:,:,2,3] > 0)
0088                 sel_p, sel_r = sel  # 2-tuple selecting photon and record points
0089                 sel_p = sel_p[sli]
0090                 sel_r = sel_r[sli]
0091             else:
0092                 sel_p = e.q_startswith(SELECT)
0093                 sel_p = sel_p[sli]
0094
0095                 SELECT_ELEM = SELECT.split()
0096                 SELECT_POINT = len(SELECT_ELEM)
0097
0098                 sel_r = slice(0,SELECT_POINT)
0099                 r = e.f.record[sel_p,sel_r]
0100
0101                 if "SAVE_SEL" in os.environ:
0102                     sel_dir = os.path.join(e.f.base, SELECT.replace(" ", "_") )
0103                     if not os.path.isdir(sel_dir):
0104                         os.makedirs(sel_dir)
0105                     pass
0106                     sel_path = os.path.join( sel_dir, "record.npy" )
0107                     print("REC=%s ~/o/examples/UseGeometryShader/run.sh" % sel_dir)
0108                     np.save( sel_path, e.f.record[sel_p] )
0109                 pass
0110
0111                 for i in range(SELECT_POINT-1):
0112                     speed = speed_(r,i)
0113                     speed_min = speed.min() if len(speed) > 0 else -1
0114                     speed_max = speed.max() if len(speed) > 0 else -1
0115                     fmt = "speed len/min/max for : %d -> %d : %s -> %s : %8d %7.3f %7.3f "
0116                     print(fmt % (i,i+1,SELECT_ELEM[i],SELECT_ELEM[i+1],len(speed),speed_min, speed_max ))
0117                 pass
0118
0119                 if Q == "B":
0120                     kludge = e.f.NPFold_meta.U4Recorder__PreUserTrackingAction_Optical_UseGivenVelocity_KLUDGE
0121                     print("e.f.NPFold_meta.U4Recorder__PreUserTrackingAction_Optical_UseGivenVelocity_KLUDGE:%s " % kludge )
0122                     print("e.f.NPFold_meta.G4VERSION_NUMBER:%s " % e.f.NPFold_meta.G4VERSION_NUMBER )
0123                 pass
0124             pass
0125
0126             _pos = e.f.photon[sel_p,0,:3]        ## end position
0127
0128             _beg = e.f.record[sel_p,0,0,:3]      ## begin positions
0129             _poi = e.f.record[sel_p,sel_r,0,:3]  ## all positions
0130
0131             print("_pos.shape %s " % str(_pos.shape))
0132             print("_beg.shape %s " % str(_beg.shape))
0133             print("_poi.shape %s " % str(_poi.shape))
0134
0135             pos = _pos.reshape(-1,3)
0136             beg = _beg.reshape(-1,3)
0137             poi = _poi.reshape(-1,3)
0138
0139             elabel = "%s : %s " % ( e.symbol.upper(), e.f.base )
0140             label = context + " ## " + elabel
0141
0142             if MODE == 3 and not pvp is None:
0143                 if len(poi) == 0:
0144                     print("FAILED TO SELECT ANY POI : SKIP PLOTTING" )
0145                 else:
0146                     pl = pvp.pvplt_plotter(label)
0147                     pvp.pvplt_viewpoint(pl) # sensitive EYE, LOOK, UP, ZOOM envvars eg EYE=0,-3,0
0148                     pl.add_points( poi, color="green", point_size=3.0 )
0149                     pl.add_points( pos, color="red", point_size=3.0 )
0150                     pl.add_points( beg, color="blue", point_size=3.0 )
0151                     if not "NOGRID" in os.environ: pl.show_grid()
0152                     cp = pl.show()
0153                 pass
0154             pass
0155         pass
0156     pass
0157
0158