eic-opticks/CSGOptiX/cxs_min.py

0001 #!/usr/bin/env python
0002
0003 import os, logging, numpy as np
0004 log = logging.getLogger(__name__)
0005
0006 print("[from opticks.sysrap.sevt import SEvt, SAB")
0007 from opticks.sysrap.sevt import SEvt, SAB
0008 print("]from opticks.sysrap.sevt import SEvt, SAB")
0009
0010
0011
0012
0013 TEST = os.environ.get("TEST","")
0014 PLOT = os.environ.get("PLOT","scatter")
0015 GLOBAL = int(os.environ.get("GLOBAL","0")) == 1
0016 MODE = int(os.environ.get("MODE","3"))
0017 SEL = int(os.environ.get("SEL","0"))
0018
0019 print("before pvplt import MODE:%d " % MODE)
0020
0021 if MODE in [2,3]:
0022     from opticks.ana.pvplt import *
0023     # HMM this import overrides MODE, so need to keep defaults the same
0024 pass
0025
0026 if __name__ == '__main__':
0027     logging.basicConfig(level=logging.INFO)
0028     print("GLOBAL:%d MODE:%d SEL:%d" % (GLOBAL,MODE, SEL))
0029
0030     a = SEvt.Load("$AFOLD", symbol="a")
0031     print(repr(a))
0032
0033     if "BFOLD" in os.environ:
0034         b = SEvt.Load("$BFOLD", symbol="b")
0035         print(repr(b))
0036         ab = SAB(a,b)
0037         print(repr(ab))
0038     pass
0039
0040
0041     e = a
0042
0043     qtab = e.minimal_qtab()  # requires seq array saving to be useful
0044     print("qtab")
0045     print(qtab)
0046
0047
0048     if hasattr(e.f, 'record'):
0049         pp = e.f.record[:,:,0,:3].reshape(-1,3)
0050         found = "record"
0051     elif hasattr(e.f, 'hit'):
0052         pp = e.f.hit[:,0,:3]
0053         found = "hit"
0054     elif hasattr(e.f, 'photon'):
0055         pp = e.f.photon[:,0,:3]
0056         found = "photon"
0057     elif hasattr(e.f, 'inphoton'):
0058         pp = e.f.inphoton[:,0,:3]
0059         found = "inphoton"
0060     elif hasattr(e.f, 'genstep'):
0061         pp = e.f.genstep[:,1,:3]
0062         found = "genstep"
0063     else:
0064         pp = None
0065         found = "NONE"
0066         pass
0067     pass
0068
0069     assert not pp is None
0070
0071     gpos = np.ones( [len(pp), 4 ] )
0072     gpos[:,:3] = pp
0073     lpos = np.dot( gpos, e.f.sframe.w2m )
0074     upos = gpos if GLOBAL else lpos
0075
0076     if MODE in [0,1]:
0077         print("not plotting as MODE %d in environ" % MODE )
0078     elif MODE == 2:
0079         label = "%s : " % ( e.f.base.replace("/data/blyth/opticks/","") )
0080         label +=  " : %s " % PLOT
0081
0082         expl = ""
0083         if PLOT.startswith("seqnib") and hasattr(e.f, "seqnib"):
0084             expl = "Photon History Step Counts Occurrence in single 1M photon event"
0085         elif PLOT.startswith("thit") and hasattr(e.f, "hit"):
0086             expl = "Histogram hit times[ns] of all(and step tail) : from  1M photon TORCH event"
0087         pass
0088         title = "\n".join([label, expl])
0089         pl = mpplt_plotter(label=title)
0090         fig, axs = pl
0091         assert len(axs) == 1
0092         ax = axs[0]
0093     elif MODE == 3:
0094
0095         label = "%s : found %s " % ( e.f.base, found )
0096         print(label)
0097         pl = pvplt_plotter(label)
0098         pvplt_viewpoint(pl)   # sensitive to EYE, LOOK, UP envvars
0099         pvplt_frame(pl, e.f.sframe, local=not GLOBAL )
0100     pass
0101
0102
0103     H,V = 0,2  # X, Z
0104
0105     if SEL == 1:
0106         sel = np.logical_and( np.abs(upos[:,H]) < 500, np.abs(upos[:,V]) < 500 )
0107         spos = upos[sel]
0108     else:
0109         spos = upos
0110     pass
0111
0112     if MODE == 2:
0113
0114         if PLOT.startswith("scatter"):
0115             ax.scatter( spos[:,H], spos[:,V], s=0.1 )
0116         elif PLOT.startswith("seqnib") and hasattr(e.f, "seqnib"):
0117             seqnib = e.f.seqnib
0118
0119             bounce = np.arange(len(seqnib))
0120             ax.set_aspect('auto')
0121             ax.plot( bounce, seqnib )
0122             ax.scatter( bounce, seqnib )
0123             ax.set_ylabel("Photon counts for step points (0->31)", fontsize=20 )
0124
0125             cs_seqnib = np.cumsum(seqnib)
0126             ax2 = ax.twinx()
0127             ax2.plot( bounce, cs_seqnib, linestyle="dashed" )
0128             ax2.set_ylabel( "Cumulative counts rising to total of 1M photons", fontsize=20 )
0129         elif PLOT.startswith("thit") and hasattr(e.f, "hit"):
0130             hit_time = e.f.hit[:,0,3]
0131             mn_hit_time = hit_time.min()
0132             mx_hit_time = hit_time.max()
0133             bins = np.linspace( mn_hit_time, mx_hit_time, 100 )
0134
0135             hit_time_n,_  = np.histogram(hit_time, bins=bins )
0136             ax.set_aspect('auto')
0137             ax.plot( bins[:-1], hit_time_n, drawstyle="steps-mid", label="Simulated Hit time [ns]" )
0138
0139             hit_idx = e.f.hit.view(np.uint32)[:,3,2] & 0x7fffffff
0140             hit_nib = a.f.seqnib[hit_idx]   # nibbles of the hits
0141
0142             cut = 23
0143             hit_time_sel = hit_time[hit_nib>cut]
0144             sel_time_n, _ = np.histogram(hit_time_sel, bins=bins)
0145
0146             ax.plot( bins[:-1], sel_time_n, drawstyle="steps-mid", label="Simulated Hit time [ns] for step points > %d " % cut  )
0147
0148             ax.set_ylabel("Photon counts in simulated time[ns] bins", fontsize=20)
0149             ax.set_xlabel("Simulated time[ns]", fontsize=20)
0150
0151             ax.set_yscale('log')
0152             ax.legend()
0153             # need the seqnib of each photon so can check times of the big seqnib
0154             # PLOT=thit MODE=2 ~/o/cxs_min.sh
0155         else:
0156             print("PLOT:%s unhandled" % PLOT)
0157         pass
0158     elif MODE == 3:
0159         if PLOT.startswith("scatter"):
0160             pl.add_points(spos[:,:3])
0161         else:
0162             pass
0163         pass
0164     else:
0165         pass
0166     pass
0167
0168     if MODE == 2:
0169         fig.show()
0170     elif MODE == 3:
0171         pl.show()
0172     pass
0173 pass
0174
0175
0176
0177