ana/debug/types.py

0001 #!/bin/env python
0002 #
0003 # Copyright (c) 2019 Opticks Team. All Rights Reserved.
0004 #
0005 # This file is part of Opticks
0006 # (see https://bitbucket.org/simoncblyth/opticks).
0007 #
0008 # Licensed under the Apache License, Version 2.0 (the "License");
0009 # you may not use this file except in compliance with the License.
0010 # You may obtain a copy of the License at
0011 #
0012 #   http://www.apache.org/licenses/LICENSE-2.0
0013 #
0014 # Unless required by applicable law or agreed to in writing, software
0015 # distributed under the License is distributed on an "AS IS" BASIS,
0016 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
0017 # See the License for the specific language governing permissions and
0018 # limitations under the License.
0019 #
0020
0021 """
0022 types.py : MOSTLY REPLACED BY OTHER MODULES
0023 ==============================================
0024
0025 ::
0026
0027     simon:ana blyth$ grep ana.types *.py
0028
0029     fresnel.py:from opticks.ana.types import *
0030     reflection.py:from opticks.ana.types import *
0031
0032
0033
0034 """
0035 import os, datetime, logging, json
0036 log = logging.getLogger(__name__)
0037 import numpy as np
0038
0039 from opticks.ana.base import json_, ffs_, ihex_
0040 from opticks.ana.nbase import count_unique
0041
0042 lsb2_ = lambda _:(_ & 0xFF)
0043 msb2_ = lambda _:(_ & 0xFF00) >> 8
0044 hex_ = lambda _:"0x%x" % _
0045
0046 idp_ = lambda _:os.path.expandvars("$IDPATH/%s" % _ )
0047
0048 pro_ = lambda _:load_("prop",_)
0049 ppp_ = lambda _:load_("photon",_)
0050 hhh_ = lambda _:load_("hit",_)
0051 ttt_ = lambda _:load_("test",_)
0052 stc_ = lambda _:load_("cerenkov",_)
0053 sts_ = lambda _:load_("scintillation",_)
0054 chc_ = lambda _:load_("opcerenkov",_)
0055 chs_ = lambda _:load_("opscintillation",_)
0056 oxc_ = lambda _:load_("oxcerenkov",_)
0057 auc_ = lambda _:load_("aucerenkov",_)
0058 oxs_ = lambda _:load_("oxscintillation",_)
0059 aus_ = lambda _:load_("auscintillation",_)
0060 rxc_ = lambda _:load_("rxcerenkov",_)
0061 rxs_ = lambda _:load_("rxscintillation",_)
0062 dom_ = lambda _:load_("domain",_)
0063 idom_ = lambda _:load_("idomain",_)
0064 seqc_ = lambda _:load_("seqcerenkov",_)
0065 seqs_ = lambda _:load_("seqscintillation",_)
0066 phc_ =  lambda _:load_("phcerenkov",_)
0067 phs_ =  lambda _:load_("phscintillation",_)
0068
0069 recsel_cerenkov_ = lambda _:load_("recsel_cerenkov", _)
0070 phosel_cerenkov_ = lambda _:load_("phosel_cerenkov", _)
0071 recsel_scintillation_ = lambda _:load_("recsel_scintillation", _)
0072 phosel_scintillation_ = lambda _:load_("phosel_scintillation", _)
0073
0074
0075
0076 g4c_ = lambda _:load_("gopcerenkov",_)
0077 g4s_ = lambda _:load_("gopscintillation",_)
0078 pmt_ = lambda _:load_("pmthit",_)
0079
0080
0081 def mat_(path="~/.opticks/GMaterialIndexLocal.json"):
0082     """
0083     Customized names to codes arranged to place
0084     more important materials at indices less than 0xF::
0085
0086         simon:~ blyth$ cat ~/.opticks/GMaterialIndexLocal.json
0087         {
0088             "ADTableStainlessSteel": "19",
0089             "Acrylic": "3",
0090             "Air": "15",
0091             "Aluminium": "18",
0092             "Bialkali": "5",
0093             "DeadWater": "8",
0094             "ESR": "10",
0095             "Foam": "20",
0096             "GdDopedLS": "1",
0097             "IwsWater": "6",
0098             "LiquidScintillator": "2",
0099             "MineralOil": "4",
0100
0101     """
0102     js = json_(path)
0103     return dict(zip(map(str,js.keys()),map(int,js.values())))
0104
0105
0106
0107 def imat_():
0108     """
0109     Inverse mat providing names for the custom integer codes::
0110
0111         In [90]: im[0xF]
0112         Out[90]: 'Air'
0113
0114         im = imat_()
0115         im_ = lambda _:im.get(_,'?')
0116
0117     """
0118     mat = mat_()
0119     return dict(zip(map(int,mat.values()),map(str,mat.keys())))
0120
0121 def cmm_(path="$IDPATH/../ChromaMaterialMap.json"):
0122     """
0123     Longname to chroma material code::
0124
0125         simon:~ blyth$ cat $IDPATH/../ChromaMaterialMap.json
0126         {"/dd/Materials/OpaqueVacuum": 18, "/dd/Materials/Pyrex": 21, "/dd/Materials/PVC": 20, "/dd/Materials/NitrogenGas": 16,
0127
0128     """
0129     js = json_(path)
0130     return dict(zip(map(lambda _:str(_)[len("/dd/Materials/"):],js.keys()),map(int,js.values())))
0131
0132 def icmm_():
0133     cmm = cmm_()
0134     return dict(zip(cmm.values(),cmm.keys()))
0135
0136
0137 def lmm_(path="$IDPATH/GBoundaryLibMetadataMaterialMap.json"):
0138     """
0139     Shortname to wavelength line number::
0140
0141         simon:~ blyth$ cat $IDPATH/GBoundaryLibMetadataMaterialMap.json
0142         {
0143             "ADTableStainlessSteel": "330",
0144             "Acrylic": "84",
0145             "Air": "12",
0146             "Aluminium": "24",
0147             "Bialkali": "126",
0148             "DeadWater": "42",
0149              ...
0150
0151     """
0152     js = json_(path)
0153     return dict(zip(map(str,js.keys()),map(int,js.values())))
0154
0155
0156
0157 def bnd_(path="$IDPATH/GBoundaryLibMetadata.json"):
0158     js = json_(path)
0159
0160     boundary = js['lib']['boundary']
0161
0162     bnd = {}
0163     bnd[0] = "MISS"
0164
0165     def shorten_sur_(sur):
0166         return sur if len(sur) < 2 else sur.split("__")[-1]
0167
0168     for i in sorted(map(int,boundary.keys())):
0169         b = boundary[str(i)]
0170         j = i + 1
0171         imat = b['imat']['shortname']
0172         omat = b['omat']['shortname']
0173         isur = shorten_sur_(b['isur']['name'])
0174         osur = shorten_sur_(b['osur']['name'])
0175         d = "%s/%s/%s/%s" % ( imat, omat, isur, osur )
0176
0177         bnd[j]   = "(+%0.2d) %s " % (j,d)
0178         bnd[-j] =  "(-%0.2d) %s " % (j,d)
0179     pass
0180     return bnd
0181
0182
0183 def c2g_():
0184     """
0185     From chroma material indice to ggeo customized
0186     """
0187     cmm = cmm_()
0188     gmm = mat_()
0189     return dict(zip(map(lambda _:cmm.get(_,-1),gmm.keys()),gmm.values()))
0190
0191 def c2l_():
0192     """
0193     Equivalent to G4StepNPY lookup
0194     providing mapping from chroma material index that is present in the gensteps
0195     to GGeo wavelength texture line number
0196     """
0197     cmm = cmm_()  # chroma material map
0198     lmm = lmm_()  # ggeo line numbers into wavelength texture
0199     return dict(zip(map(lambda _:cmm.get(_,-1),lmm.keys()),lmm.values()))
0200
0201 def check_c2l():
0202     """
0203     Rock is discrepant as not present in cmm
0204     """
0205     c2l = c2l_()    # int to int mapping
0206
0207     icmm = icmm_()  # chroma int to name
0208     mat = mat_()    # ggeo name to int   (customized(
0209     line2mat = line2mat_()   # wavelength buffer line
0210
0211     cnam = map(lambda _:icmm.get(_,None),c2l.keys())
0212     lnam = map(lambda _:line2mat.get(_,None),c2l.values())
0213     assert lnam[:-1] == cnam[:-1]
0214     print set(lnam) - set(cnam)
0215
0216
0217
0218 def check_gs():
0219     """
0220     In [71]: map(lambda _:icmm.get(_,None),np.unique(gsmat))
0221     Out[71]:
0222     ['Acrylic',
0223      'DeadWater',
0224      'GdDopedLS',
0225      'IwsWater',
0226      'LiquidScintillator',
0227      'MineralOil',
0228      'OwsWater']
0229     """
0230
0231     gs = stc_(1)
0232     icmm = icmm_()
0233     gsmat = gs.view(np.int32)[:,0,2]
0234     gnam = map(lambda _:icmm.get(_,None),np.unique(gsmat))
0235
0236
0237
0238 _line2g = {}
0239
0240 def line2g_():
0241     """
0242     ::
0243
0244         In [93]: l2g = line2g_()
0245
0246         In [96]: l2g
0247         Out[96]:
0248         {0: 13,
0249          1: 13,
0250          6: 16,
0251          7: 13,
0252          12: 15,
0253          13: 16,
0254          18: 17,
0255
0256
0257     """
0258     global _line2g
0259     if _line2g:
0260         return _line2g
0261     o = np.load(os.path.expandvars("$IDPATH/optical.npy")).reshape(-1,6,4)
0262     _line2g = {}
0263     for i,b in enumerate(o):
0264         _line2g[i*6+0] = b[0,0]
0265         _line2g[i*6+1] = b[1,0]
0266     return _line2g
0267
0268 def line2mat_():
0269     im = imat_()
0270     line2g = line2g_()
0271     return dict(zip(line2g.keys(),map(lambda _:im.get(_,None),line2g.values())))
0272
0273
0274 class Index(object):
0275     def table(self, cu, hex_=False):
0276         assert len(cu.shape) == 2 and cu.shape[1] == 2
0277         for msk,cnt in sorted(cu, key=lambda _:_[1], reverse=True):
0278             label = ihex_(msk) if hex_ else int(msk)
0279             log.debug("Index msk %d cnt %d label %s " % (msk, cnt, label))
0280             print "%20s %10d : %40s " % ( label, int(cnt), self(msk) )
0281         pass
0282
0283 class GFlags(Index):
0284     def __init__(self):
0285         self.mf = maskflags_()
0286         self.skip = ["TORCH"]
0287     def __call__(self, i):
0288         return "|".join(map(lambda kv:kv[1], filter(lambda kv:kv[1] not in self.skip,filter(lambda kv:i & kv[0], self.mf.items()))))
0289
0290 def gflags_table(cu):
0291     gf = GFlags()
0292     gf.table(cu)
0293
0294
0295 class SeqHis(Index):
0296     def __init__(self, abbrev=True):
0297         self.f  = abflags_() if abbrev else flags_()
0298         self.fi = iabflags_() if abbrev else iflags_()
0299     def __call__(self, i):
0300         x = ihex_(i)
0301         log.debug("seqhis %s " % x )
0302         return " ".join(map(lambda _:self.fi.get(int(_,16),'?%s?' % int(_,16) ), x[::-1] ))
0303     def seqhis_int(self, s):
0304         f = self.f
0305         return reduce(lambda a,b:a|b,map(lambda ib:ib[1] << 4*ib[0],enumerate(map(lambda n:f[n], s.split(" ")))))
0306
0307 def seqhis_int(s, abbrev=True):
0308     f  = abflags_() if abbrev else flags_()
0309     return reduce(lambda a,b:a|b,map(lambda ib:ib[1] << 4*ib[0],enumerate(map(lambda n:f[n], s.split(" ")))))
0310
0311 def seqhis_table(cu):
0312     sh = SeqHis()
0313     sh.table(cu, hex_=True)
0314
0315 def seqhis_table_():
0316     seq = load_("phtorch","1","dayabay")
0317     seqhis = seq[:,0,0]
0318     cu = count_unique(seqhis)
0319     sh = SeqHis()
0320     sh.table(cu, hex_=True)
0321
0322
0323
0324 class GFlags(Index):
0325     def __init__(self):
0326         self.mf = maskflags_()
0327         self.skip = ["TORCH"]
0328     def __call__(self, i):
0329         return "|".join(map(lambda kv:kv[1], filter(lambda kv:kv[1] not in self.skip,filter(lambda kv:i & kv[0], self.mf.items()))))
0330
0331 def gflags_table(cu):
0332     gf = GFlags()
0333     gf.table(cu)
0334
0335 def maskflags_():
0336     iaf = iabflags_()
0337     return dict(zip(map(lambda b:0x1 << (b-1),iaf.keys()),iaf.values()))
0338
0339
0340
0341
0342
0343
0344 def maskflags_string(i, skip=["TORCH"]):
0345     mf = maskflags_()
0346     return "|".join(map(lambda kv:kv[1], filter(lambda kv:kv[1] not in skip,filter(lambda kv:i & kv[0], mf.items()))))
0347
0348 def maskflags_int(s):
0349     mf = maskflags_()
0350     fm = dict(zip(mf.values(),mf.keys()))
0351     return reduce(lambda a,b:a|b, map(lambda n:fm[n], s.split("|")))
0352
0353
0354
0355
0356 _abbrev_mat = {
0357     "NitrogenGas":"NG",
0358     "ADTableStainlessSteel":"SS",
0359     "LiquidScintillator":"LS",
0360     "MineralOil":"MO" }
0361
0362 def abbrev_mat_(mat):
0363     return _abbrev_mat.get(mat,mat[0:2])
0364
0365
0366 def abmat_():
0367     mat = mat_()
0368     return dict(zip(map(abbrev_mat_,mat.keys()),mat.values()))
0369
0370 def iabmat_():
0371     mat =  abmat_()
0372     d = dict(zip(map(int,mat.values()),map(str,mat.keys())))
0373     d[0] = '??'
0374     return d
0375
0376 def seqhis_(i, abbrev=True):
0377     fi = iabflags_() if abbrev else iflags_()
0378     x = ihex_(i)
0379     return " ".join(map(lambda _:fi[int(_,16)], x[::-1] ))
0380
0381 def seqmat_(i, abbrev=True):
0382     mi = iabmat_() if abbrev else imat_()
0383     x = ihex_(i)
0384     return " ".join(map(lambda _:mi[int(_,16)], x[::-1] ))
0385
0386
0387
0388
0389
0390
0391 def load_(typ, tag, det="dayabay"):
0392     path = path_(typ,tag, det)
0393
0394     if os.path.exists(path):
0395         log.debug("loading %s " % path )
0396         os.system("ls -l %s " % path)
0397         return np.load(path)
0398     pass
0399     log.warning("load_ no such path %s " % path )
0400     return None
0401
0402
0403 gjspath_ = lambda _:os.path.expandvars("$IDPATH/%s" % _)
0404 geopath_ = lambda _:os.path.expandvars("$IDPATH/%s.npy" % _)
0405 geoload_ = lambda _:np.load(geopath_(_))
0406
0407
0408
0409
0410 global typs
0411 typs = "photon hit test cerenkov scintillation opcerenkov opscintillation gopcerenkov gopscintillation".split()
0412
0413 global typmap
0414 typmap = {}
0415
0416 class NPY(np.ndarray):
0417     shape2type = {
0418             (4,4):"photon",
0419             (6,4):"g4step",
0420                  }
0421
0422     @classmethod
0423     def from_array(cls, arr ):
0424         return arr.view(cls)
0425
0426     @classmethod
0427     def empty(cls):
0428         a = np.array((), dtype=np.float32)
0429         return a.view(cls)
0430
0431     @classmethod
0432     def detect_type(cls, arr ):
0433         """
0434         distinguish opcerenkov from opscintillation ?
0435         """
0436         assert len(arr.shape) == 3
0437         itemshape = arr.shape[1:]
0438         typ = cls.shape2type.get(itemshape, None)
0439         if typ == "g4step":
0440             zzz = arr[0,0,0].view(np.int32)
0441             if zzz < 0:
0442                 typ = "cerenkov"
0443             else:
0444                 typ = "scintillation"
0445             pass
0446         pass
0447         return typ
0448
0449     @classmethod
0450     def get(cls, tag):
0451         """
0452         # viewing an ndarray as a subclass allows adding customizations
0453           on top of the ndarray while using the same storage
0454         """
0455         a = load_(cls.typ, tag).view(cls)
0456         a.tag = tag
0457         return a
0458
0459     label = property(lambda self:"%s.get(%s)" % (self.__class__.__name__, self.tag))
0460
0461     @classmethod
0462     def summary(cls, tag):
0463         for typ in typs:
0464             path = path_(typ,tag)
0465             if os.path.exists(path):
0466                 mt = os.path.getmtime(path)
0467                 dt = datetime.datetime.fromtimestamp(mt)
0468                 msg = dt.strftime("%c")
0469             else:
0470                 msg = "-"
0471             pass
0472             print "%20s : %60s : %s " % (typ, path, msg)
0473
0474     @classmethod
0475     def mget(cls, tag, *typs):
0476         """
0477         Load multiple typed instances::
0478
0479             chc, g4c, tst = NPY.mget(1,"opcerenkov","gopcerenkov","test")
0480
0481         """
0482         if len(typs) == 1:
0483             typs = typs[0].split()
0484
0485         klss = map(lambda _:typmap[_], typs)
0486         arys = map(lambda kls:kls.get(tag), klss)
0487         return arys
0488
0489
0490 class Record(NPY):
0491     """
0492     OpenGL normalized shorts, a form of float -1.f:1.f compression
0493     """
0494     posx         = property(lambda self:self[:,0,0])
0495     posy         = property(lambda self:self[:,0,1])
0496     posz         = property(lambda self:self[:,0,2])
0497     time         = property(lambda self:self[:,0,3])
0498
0499
0500 class Photon(NPY):
0501     posx         = property(lambda self:self[:,0,0])
0502     posy         = property(lambda self:self[:,0,1])
0503     posz         = property(lambda self:self[:,0,2])
0504     time         = property(lambda self:self[:,0,3])
0505     position     = property(lambda self:self[:,0,:3])
0506
0507     dirx         = property(lambda self:self[:,1,0])
0508     diry         = property(lambda self:self[:,1,1])
0509     dirz         = property(lambda self:self[:,1,2])
0510     wavelength   = property(lambda self:self[:,1,3])
0511     direction    = property(lambda self:self[:,1,:3])
0512
0513     polx         = property(lambda self:self[:,2,0])
0514     poly         = property(lambda self:self[:,2,1])
0515     polz         = property(lambda self:self[:,2,2])
0516     weight       = property(lambda self:self[:,2,3])
0517     polarization = property(lambda self:self[:,2,:3])
0518
0519     aux0         = property(lambda self:self[:,3,0].view(np.int32))
0520     aux1         = property(lambda self:self[:,3,1].view(np.int32))
0521     aux2         = property(lambda self:self[:,3,2].view(np.int32))
0522     aux3         = property(lambda self:self[:,3,3].view(np.int32))
0523
0524     photonid     = property(lambda self:self[:,3,0].view(np.int32))
0525     spare        = property(lambda self:self[:,3,1].view(np.int32))
0526     flgs         = property(lambda self:self[:,3,2].view(np.uint32))
0527     pmt          = property(lambda self:self[:,3,3].view(np.int32))
0528
0529     history      = property(lambda self:self[:,3,2].view(np.uint32))   # cannot name "flags" as that shadows a necessary ndarray property
0530     pmtid        = property(lambda self:self[:,3,3].view(np.int32))
0531
0532     hits         = property(lambda self:self[self.pmtid > 0])
0533     aborts       = property(lambda self:self[np.where(self.history & 1<<31)])
0534
0535
0536     last_hit_triangle = property(lambda self:self[:,3,0].view(np.int32))
0537
0538     def history_zero(self):
0539         log.info("filling history with zeros %s " % repr(self.history.shape))
0540         self.history.fill(0)
0541
0542     #def _get_last_hit_triangles(self):
0543     #    if self._last_hit_triangles is None:
0544     #        self._last_hit_triangles = np.empty(len(self), dtype=np.int32)
0545     #        self._last_hit_triangles.fill(-1)
0546     #    return self._last_hit_triangles
0547     #last_hit_triangles = property(_get_last_hit_triangles)
0548     #
0549     #def _get_history(self):
0550     #    if self._last_hit_triangles is None:
0551     #        self._last_hit_triangles = np.empty(len(self), dtype=np.int32)
0552     #        self._last_hit_triangles.fill(-1)
0553     #    return self._last_hit_triangles
0554     #last_hit_triangles = property(_get_last_hit_triangles)
0555
0556
0557
0558
0559     def dump(self, index):
0560         log.info("dump index %d " % index)
0561         print self[index]
0562         print "photonid: ", self.photonid[index]
0563         print "history: ",  self.history[index]
0564         print "pmtid: ",    self.pmtid[index]    # is this still last_hit_triangle index when not a hit ?
0565
0566
0567 class G4CerenkovPhoton(Photon):
0568     """DsChromaG4Cerenkov.cc"""
0569     typ = "gopcerenkov"
0570     cmat = property(lambda self:self[:,3,0].view(np.int32)) # chroma material index
0571     sid = property(lambda self:self[:,3,1].view(np.int32))
0572     pass
0573 typmap[G4CerenkovPhoton.typ] = G4CerenkovPhoton
0574
0575 class G4ScintillationPhoton(Photon):
0576     """DsChromaG4Scintillation.cc"""
0577     typ = "gopscintillation"
0578     cmat = property(lambda self:self[:,3,0].view(np.int32)) # chroma material index
0579     sid = property(lambda self:self[:,3,1].view(np.int32))
0580     pdg = property(lambda self:self[:,3,2].view(np.int32))
0581     scnt = property(lambda self:self[:,3,3].view(np.int32))
0582     pass
0583 typmap[G4ScintillationPhoton.typ] = G4ScintillationPhoton
0584
0585 class ChCerenkovPhoton(Photon):
0586     typ = "opcerenkov"
0587     pass
0588 typmap[ChCerenkovPhoton.typ] = ChCerenkovPhoton
0589
0590 class OxCerenkovPhoton(Photon):
0591     typ = "oxcerenkov"
0592     pass
0593 typmap[OxCerenkovPhoton.typ] = OxCerenkovPhoton
0594
0595
0596
0597
0598
0599 class ChCerenkovPhotonGen(Photon):
0600     typ = "opcerenkovgen"
0601     pass
0602 typmap[ChCerenkovPhotonGen.typ] = ChCerenkovPhotonGen
0603
0604
0605 class ChScintillationPhoton(Photon):
0606     typ = "opscintillation"
0607     pass
0608 typmap[ChScintillationPhoton.typ] = ChScintillationPhoton
0609
0610 class ChScintillationPhotonGen(Photon):
0611     typ = "opscintillationgen"
0612     pass
0613 typmap[ChScintillationPhotonGen.typ] = ChScintillationPhotonGen
0614
0615 class OxScintillationPhoton(Photon):
0616     typ = "oxscintillation"
0617     pass
0618 typmap[OxScintillationPhoton.typ] = OxScintillationPhoton
0619
0620
0621
0622 class RxScintillationRecord(Record):
0623     typ = "rxscintillation"
0624     pass
0625 typmap[RxScintillationRecord.typ] = RxScintillationRecord
0626
0627 class RxCerenkovRecord(Record):
0628     typ = "rxcerenkov"
0629     pass
0630 typmap[RxCerenkovRecord.typ] = RxCerenkovRecord
0631
0632
0633
0634
0635 class TestPhoton(Photon):
0636     typ = "test"
0637     pass
0638 typmap[TestPhoton.typ] = TestPhoton
0639
0640
0641 class ChromaPhoton(Photon):
0642     typ = "chromaphoton"
0643
0644     GENERATE_SCINTILLATION = 0x1 << 16
0645     GENERATE_CERENKOV      = 0x1 << 17
0646
0647     cerenkov      = property(lambda self:self[np.where( self.history & self.GENERATE_CERENKOV )[0]])
0648     scintillation = property(lambda self:self[np.where( self.history & self.GENERATE_SCINTILLATION )[0]])
0649
0650     @classmethod
0651     def from_arrays(cls, pos, dir, pol, wavelengths, t, last_hit_triangles, flags, weights, hit=0):
0652         """
0653         #. NB a kludge setting of pmtid into lht using the map argument of propagate_hit.cu
0654         """
0655         nall = len(pos)
0656         a = np.zeros( (nall,4,4), dtype=np.float32 )
0657         pmtid = np.zeros( nall, dtype=np.int32 )
0658
0659         a[:,0,:3] = pos
0660         a[:,0, 3] = t
0661
0662         a[:,1,:3] = dir
0663         a[:,1, 3] = wavelengths
0664
0665         a[:,2,:3] = pol
0666         a[:,2, 3] = weights
0667
0668         assert len(last_hit_triangles) == len(flags)
0669
0670         SURFACE_DETECT = 0x1 << 2
0671         detected = np.where( flags & SURFACE_DETECT  )
0672         pmtid[detected] = last_hit_triangles[detected]             # sparsely populate, leaving zeros for undetected
0673
0674         a[:,3, 0] = np.arange(nall, dtype=np.int32).view(a.dtype)  # photon_id
0675         a[:,3, 1] = 0                                              # used in comparison againt vbo prop
0676         a[:,3, 2] = flags.view(a.dtype)                            # history flags
0677         a[:,3, 3] = pmtid.view(a.dtype)                            # channel_id ie PmtId
0678
0679         if hit:
0680             pp = a[pmtid > 0].view(cls)
0681         else:
0682             pp = a.view(cls)
0683         pass
0684         return pp
0685     pass
0686 typmap[ChromaPhoton.typ] = ChromaPhoton
0687
0688
0689
0690
0691
0692
0693
0694
0695
0696 class Prop(NPY):
0697     """
0698     See test_ScintillationIntegral from gdct-
0699     """
0700     typ = "prop"
0701     flat = property(lambda self:self[:,0])  # unform draw from 0 to max ScintillationIntegral
0702     wavelength = property(lambda self:1./self[:,1])# take reciprocal to give wavelength
0703     pass
0704 typmap[Prop.typ] = Prop
0705
0706 class G4Step(NPY):
0707     typ = "g4step"
0708     sid = property(lambda self:self[:,0,0].view(np.int32))    # 0
0709     parentId = property(lambda self:self[:,0,1].view(np.int32))
0710     materialIndex = property(lambda self:self[:,0,2].view(np.int32))
0711     numPhotons = property(lambda self:self[:,0,3].view(np.int32))
0712
0713     position      = property(lambda self:self[:,1,:3])
0714     time          = property(lambda self:self[:,1,3])
0715
0716     deltaPosition = property(lambda self:self[:,2,:3])
0717     stepLength    = property(lambda self:self[:,2,3])
0718
0719
0720     code = property(lambda self:self[:,3,0].view(np.int32))   # 3
0721
0722     totPhotons = property(lambda self:int(self.numPhotons.sum()))
0723     materialIndices = property(lambda self:np.unique(self.materialIndex))
0724
0725     def materials(self, _cg):
0726         """
0727         :param _cg: chroma geometry instance
0728         :return: list of chroma material instances relevant to this evt
0729         """
0730         return [_cg.unique_materials[materialIndex] for materialIndex in self.materialIndices]
0731     pass
0732 typmap[G4Step.typ] = G4Step
0733
0734
0735 class ScintillationStep(G4Step):
0736     """
0737     see DsChromaG4Scintillation.cc
0738     """
0739     typ = "scintillation"
0740     pass
0741 typmap[ScintillationStep.typ] = ScintillationStep
0742
0743
0744 class CerenkovStep(G4Step):
0745     """
0746     see DsChromaG4Cerenkov.cc
0747     """
0748     typ = "cerenkov"
0749     BetaInverse = property(lambda self:self[:,4,0])
0750     maxSin2 = property(lambda self:self[:,5,0])
0751     bialkaliIndex = property(lambda self:self[:,5,3].view(np.int32))
0752     pass
0753 typmap[CerenkovStep.typ] = CerenkovStep
0754
0755
0756 class PmtHit(Photon):
0757     typ = "pmthit"
0758 typmap[PmtHit.typ] = PmtHit
0759
0760 class G4PmtHit(Photon):
0761     typ = "g4pmthit"
0762 typmap[G4PmtHit.typ] = G4PmtHit
0763
0764
0765
0766 class VBOMixin(object):
0767     numquad = 6
0768     force_attribute_zero = "position_weight"
0769
0770     @classmethod
0771     def vbo_from_array(cls, arr, max_slots=None):
0772         if arr is None:return None
0773         assert max_slots
0774         a = arr.view(cls)
0775         a.initialize(max_slots)
0776         return a
0777
0778     def initialize(self, max_slots):
0779         self.max_slots = max_slots
0780         self._vbodata = None
0781         self._ccolor = None
0782         self._indices = None
0783
0784     def _get_ccolor(self):
0785         """
0786         #. initialize to red, reset by CUDA kernel
0787         """
0788         if self._ccolor is None:
0789             self._ccolor = np.tile( [1.,0.,0,1.], (len(self),1)).astype(np.float32)
0790         return self._ccolor
0791     ccolor = property(_get_ccolor, doc=_get_ccolor.__doc__)
0792
0793     def ccolor_from_code(self):
0794         ccolor = np.tile( [1.,1.,1.,1.], (len(self),1)).astype(np.float32)
0795         ccolor[np.where(self.code == 13),:] = [1,0,0,1]  #     mu:red
0796         ccolor[np.where(self.code == 11),:] = [0,1,0,1]  #      e:green
0797         ccolor[np.where(self.code == 22),:] = [0,0,1,1]  #  gamma:blue
0798         return ccolor
0799
0800     def _get_indices(self):
0801         """
0802         List of indices
0803         """
0804         if self._indices is None:
0805             self._indices = np.arange( len(self), dtype=np.uint32)
0806         return self._indices
0807     indices = property(_get_indices, doc=_get_indices.__doc__)
0808
0809     def _get_vbodata(self):
0810         if self._vbodata is None:
0811             self._vbodata = self.create_vbodata(self.max_slots)
0812         return self._vbodata
0813     vbodata = property(_get_vbodata)
0814
0815
0816
0817
0818 class VBOStep(G4Step, VBOMixin):
0819     def create_vbodata(self, max_slots):
0820         """
0821         """
0822         if len(self) == 0:return None
0823         dtype = np.dtype([
0824             ('position_time'   ,          np.float32, 4 ),
0825             ('direction_wavelength',      np.float32, 4 ),
0826             ('polarization_weight',       np.float32, 4 ),
0827             ('ccolor',                    np.float32, 4 ),
0828             ('flags',                     np.uint32,  4 ),
0829             ('last_hit_triangle',         np.int32,   4 ),
0830           ])
0831         assert len(dtype) == self.numquad
0832
0833         data = np.zeros(len(self)*max_slots, dtype )
0834         log.info( "create_data items %d max_slots %d nvert %d (with slot scaleups) " % (len(self),max_slots,len(data)) )
0835
0836
0837         def pack31_( name, a, b ):
0838             data[name][::max_slots,:3] = a
0839             data[name][::max_slots,3] = b
0840         def pack1_( name, a):
0841             data[name][::max_slots,0] = a
0842         def pack4_( name, a):
0843             data[name][::max_slots] = a
0844
0845         pack31_( 'position_time',        self.position ,    self.time )
0846         pack31_( 'direction_wavelength', self.deltaPosition, self.time )
0847
0848         # not setting, leaving at zero initially
0849         #pack31_( 'polarization_weight',  self.polarization, self.weight  )
0850         #pack1_(  'flags',                self.history )            # flags is used already by numpy
0851         #pack1_(  'last_hit_triangle',    self.last_hit_triangle )
0852
0853         # attempting to get gensteps sensitive to time selection
0854         data['flags'][::max_slots, 1] = self.time
0855         data['flags'][::max_slots, 2] = self.time
0856
0857         pack4_(  'ccolor',               self.ccolor_from_code())
0858
0859         return data
0860
0861
0862
0863 class VBOPhoton(Photon, VBOMixin):
0864     @classmethod
0865     def from_vbo_propagated(cls, vbo ):
0866         """
0867         Pulling the correct slot (-2) ?
0868         Hmm seems to be pulling all slots ?
0869         """
0870         r = np.zeros( (len(vbo),4,4), dtype=np.float32 )
0871         r[:,0,:4] = vbo['position_time']
0872         r[:,1,:4] = vbo['direction_wavelength']
0873         r[:,2,:4] = vbo['polarization_weight']
0874         r[:,3,:4] = vbo['last_hit_triangle'].view(r.dtype) # must view as target type to avoid coercion of int32 data into float32
0875         return r.view(cls)
0876
0877     def create_vbodata(self, max_slots):
0878         """
0879         :return: numpy named constituent array with numquad*quads structure
0880
0881         Trying to replace DAEPhotonsData
0882
0883         The start data is splayed out into the slots, leaving loadsa free slots
0884         this very sparse data structure with loadsa empty space limits
0885         the number of photons that can me managed, but its for visualization
0886         anyhow so do not need more than 100k or so.
0887
0888         Caution sensitivity to data structure naming:
0889
0890         #. using 'position' would use traditional glVertexPointer furnishing gl_Vertex to shader
0891         #. using smth else eg 'position_weight' uses generic attribute ,
0892            which requires force_attribute_zero for anythinh to appear
0893
0894         """
0895         if len(self) == 0:return None
0896         dtype = np.dtype([
0897             ('position_time'   ,          np.float32, 4 ),
0898             ('direction_wavelength',      np.float32, 4 ),
0899             ('polarization_weight',       np.float32, 4 ),
0900             ('ccolor',                    np.float32, 4 ),
0901             ('flags',                     np.uint32,  4 ),
0902             ('last_hit_triangle',         np.int32,   4 ),
0903           ])
0904         assert len(dtype) == self.numquad
0905
0906         data = np.zeros(len(self)*max_slots, dtype )
0907         log.info( "create_data items %d max_slots %d nvert %d (with slot scaleups) " % (len(self),max_slots,len(data)) )
0908
0909         def pack31_( name, a, b ):
0910             data[name][::max_slots,:3] = a
0911             data[name][::max_slots,3] = b
0912         def pack1_( name, a):
0913             data[name][::max_slots,0] = a
0914         def pack4_( name, a):
0915             data[name][::max_slots] = a
0916
0917
0918         pack31_( 'position_time',        self.position ,    self.time )
0919         pack31_( 'direction_wavelength', self.direction,    self.wavelength )
0920         pack31_( 'polarization_weight',  self.polarization, self.weight  )
0921
0922         # not setting leaving at zero
0923         #pack1_(  'flags',                self.history )            # flags is used already by numpy
0924         #pack1_(  'last_hit_triangle',    self.last_hit_triangle )
0925
0926         pack4_(  'ccolor',               self.ccolor)
0927
0928         return data
0929
0930
0931
0932
0933
0934 if __name__ == '__main__':
0935     pass
0936
0937
0938