eic-opticks/ana/simtrace_positions.py

0001 #!/usr/bin/env python
0002
0003 import os, numpy as np
0004 import logging
0005 from opticks.ana.axes import X, Y, Z
0006
0007 log = logging.getLogger(__name__)
0008
0009
0010 class SimtracePositions(object):
0011     """
0012     Transforms global intersect positions into local frame
0013
0014     HMM: the local frame positions are model frame in extent units
0015     when using tangential ... not so convenient would be better
0016     with real mm dimensions in local : kludge this with local_extent_scale=True
0017
0018     The isect, gpos used here come from qevent::add_simtrace
0019     """
0020
0021     @classmethod
0022     def Check(cls, p):
0023         num_photon = len(p)
0024         distance = p[:,0,3]
0025         landing_count = np.count_nonzero( distance )
0026         landing_msg = "ERROR NO PHOTON LANDED" if landing_count == 0 else ""
0027         print(" num_photon: %d : landing_count : %d   %s " % (num_photon, landing_count, landing_msg) )
0028
0029
0030     def __init__(self, simtrace, gs, frame, local=True, mask="pos", symbol="t_pos" ):
0031         """
0032         :param simtrace: "photons" array
0033         :param gs: FrameGensteps instance, used for gs.lim position masking
0034         :param frame: sframe instance
0035         :param local:
0036         :param mask:
0037
0038
0039         The simtrace array is populated by:
0040
0041         1. cx/CSGOptiX7.cu:simtrace
0042         2. CPU version of this ?
0043
0044         271 static __forceinline__ __device__ void simtrace( const uint3& launch_idx, const uint3& dim, quad2* prd )
0045         272 {
0046         ...
0047         274     sevent* evt  = params.evt ;
0048         280     const quad6& gs     = evt->genstep[genstep_id] ;
0049         281
0050         282     qsim* sim = params.sim ;
0051         283     RNG rng = sim->rngstate[idx] ;
0052         284
0053         285     quad4 p ;
0054         286     sim->generate_photon_simtrace(p, rng, gs, idx, genstep_id );
0055         287
0056         288     const float3& pos = (const float3&)p.q0.f  ;
0057         289     const float3& mom = (const float3&)p.q1.f ;
0058         290
0059         291     trace(
0060         292         params.handle,
0061         293         pos,
0062         294         mom,
0063         295         params.tmin,
0064         296         params.tmax,
0065         297         prd
0066         298     );
0067         299
0068         300     evt->add_simtrace( idx, p, prd, params.tmin );
0069         301
0070         302 }
0071
0072         410 SEVENT_METHOD void sevent::add_simtrace( unsigned idx, const quad4& p, const quad2* prd, float tmin )
0073         411 {
0074         412     float t = prd->distance() ;  // q0.f.w
0075         413     quad4 a ;
0076         414
0077         415     a.q0.f  = prd->q0.f ;
0078         416
0079         417     a.q1.f.x = p.q0.f.x + t*p.q1.f.x ;
0080         418     a.q1.f.y = p.q0.f.y + t*p.q1.f.y ;
0081         419     a.q1.f.z = p.q0.f.z + t*p.q1.f.z ;
0082         420     a.q1.i.w = 0.f ;
0083         421
0084         422     a.q2.f.x = p.q0.f.x ;
0085         423     a.q2.f.y = p.q0.f.y ;
0086         424     a.q2.f.z = p.q0.f.z ;
0087         425     a.q2.u.w = prd->boundary() ; // was tmin, but expecting bnd from CSGOptiXSimtraceTest.py:Photons
0088         426
0089         427     a.q3.f.x = p.q1.f.x ;
0090         428     a.q3.f.y = p.q1.f.y ;
0091         429     a.q3.f.z = p.q1.f.z ;
0092         430     a.q3.u.w = prd->identity() ;  // identity from __closesthit__ch (( prim_idx & 0xffff ) << 16 ) | ( instance_id & 0xffff )
0093         431
0094         432     simtrace[idx] = a ;
0095         433 }
0096         """
0097         local_extent_scale = frame.coords == "RTP"  ## KINDA KLUDGE DUE TO EXTENT HANDLING BEING DONE BY THE RTP TRANSFORM
0098         isect = simtrace[:,0]
0099
0100         gpos = simtrace[:,1].copy()              # global frame intersect positions
0101         gpos[:,3] = 1.
0102
0103         lpos = np.dot( gpos, frame.w2m )   # local frame intersect positions
0104
0105         if local and local_extent_scale:
0106             extent = frame.ce[3]
0107             lpos[:,:3] *= extent
0108         pass
0109
0110         upos = lpos if local else gpos   # local usually True
0111
0112         poslim = {}
0113         poslim[X] = np.array([upos[:,X].min(), upos[:,X].max()])
0114         poslim[Y] = np.array([upos[:,Y].min(), upos[:,Y].max()])
0115         poslim[Z] = np.array([upos[:,Z].min(), upos[:,Z].max()])
0116
0117         self.poslim = poslim
0118         self.gs = gs
0119         self.frame = frame
0120
0121         self.isect = isect
0122         self.gpos = gpos
0123         self.lpos = lpos
0124         self.local = local
0125
0126         ## NB when applying the mask the below are changed
0127         ## note that need to exclude the nan for comparisons to work howver
0128         ##      np.all( t_pos.simtrace[:,:3] == t.simtrace[:,:3] )
0129         self.simtrace = simtrace
0130         self.upos = upos
0131
0132         #self.make_histogram()
0133
0134         if mask == "pos":
0135             self.apply_pos_mask()
0136         elif mask == "t":
0137             self.apply_t_mask()
0138         else:
0139             pass
0140         pass
0141         self.symbol = symbol
0142
0143
0144     def __repr__(self):
0145         symbol = self.symbol
0146         return "\n".join([
0147                    "SimtracePositions",
0148                    "%s.simtrace %s " % (symbol, str(self.simtrace.shape)),
0149                    "%s.isect %s " % (symbol, str(self.isect.shape)),
0150                    "%s.gpos %s " % (symbol, str(self.gpos.shape)),
0151                    "%s.lpos %s " % (symbol, str(self.lpos.shape)),
0152                    ])
0153
0154
0155
0156     def apply_mask(self, mask):
0157         """
0158         applying the mask changes makes a selection on the self.p and self.upos arrays
0159         which will typically decrease their sizes
0160         """
0161         self.mask = mask
0162         self.simtrace = self.simtrace[mask]
0163         self.upos = self.upos[mask]
0164         self.upos2simtrace = np.where(mask)[0]   # map upos indices back to simtrace indices before mask applied
0165
0166     def apply_pos_mask(self):
0167         """
0168         pos_mask restricts upos positions to be within the limits defined by
0169         the source genstep positions
0170         """
0171         lim = self.gs.lim
0172
0173         xmin, xmax = lim[0]
0174         ymin, ymax = lim[1]
0175         zmin, zmax = lim[2]
0176
0177         upos = self.upos
0178
0179         xmask = np.logical_and( upos[:,0] >= xmin, upos[:,0] <= xmax )
0180         ymask = np.logical_and( upos[:,1] >= ymin, upos[:,1] <= ymax )
0181         zmask = np.logical_and( upos[:,2] >= zmin, upos[:,2] <= zmax )
0182
0183         xy_mask = np.logical_and( xmask, ymask )
0184         xyz_mask = np.logical_and( xy_mask, zmask )
0185         mask = xyz_mask
0186
0187         log.info("apply_pos_mask")
0188         self.apply_mask(mask)
0189
0190     def apply_t_mask(self):
0191         """
0192         t_mask restricts the intersect distance t to be greater than zero
0193         this excludes misses
0194         """
0195         log.info("apply_t_mask")
0196         #t = self.simtrace[:,2,2]
0197         t = self.simtrace[:,0,3]
0198         mask = t > 0.
0199         self.apply_mask( mask)
0200
0201     def make_histogram(self):
0202         """
0203         TODO: use  3d histo like this to sparse-ify gensteps positions,
0204         to avoiding shooting rays from big voids
0205         """
0206         lim = self.gs.lim
0207         nx = self.frame.nx
0208         ny = self.frame.ny
0209         nz = self.frame.nz
0210         upos = self.upos
0211
0212         # bizarrely some python3 versions think the below are SyntaxError without the num=
0213         #      SyntaxError: only named arguments may follow *expression
0214
0215         binx = np.linspace(*lim[X], num=2*nx+1)
0216         biny = np.linspace(*lim[Y], num=max(2*ny+1,2) )
0217         binz = np.linspace(*lim[Z], num=2*nz+2)
0218
0219         bins = ( binx, biny, binz )
0220
0221         h3d, bins2 = np.histogramdd(upos[:,:3], bins=bins )
0222
0223         self.h3d = h3d
0224         self.bins = bins
0225         self.bins2 = bins2
0226
0227
0228 if __name__ == '__main__':
0229     pass
0230