eic-opticks/analytic/gdmlchk.py

0001 #!/usr/bin/env python
0002 """
0003 gdmlchk.py
0004 ===========
0005
0006 ::
0007
0008     ipython --pdb -i gdmlchk.py
0009
0010 """
0011 import os, sys, logging
0012 from collections import OrderedDict as odict
0013 from fnmatch import fnmatch
0014
0015 log = logging.getLogger(__name__)
0016
0017 import numpy as np
0018 import lxml.etree as ET
0019 import lxml.html as HT
0020
0021 #np.seterr(all='raise')
0022
0023 try:
0024     import matplotlib.pyplot as plt
0025 except ImportError:
0026     plt = None
0027 pass
0028
0029
0030 tostring_ = lambda _:ET.tostring(_)
0031 exists_ = lambda _:os.path.exists(os.path.expandvars(_))
0032 parse_ = lambda _:ET.parse(os.path.expandvars(_)).getroot()
0033 fparse_ = lambda _:HT.fragments_fromstring(file(os.path.expandvars(_)).read())
0034 pp_ = lambda d:"\n".join([" %30s : %f " % (k,d[k]) for k in sorted(d.keys())])
0035 values_ = lambda _:np.fromstring(_, sep=" ")
0036
0037
0038
0039
0040 class P(object):
0041     def __init__(self, property_, m):
0042         name = property_.attrib["name"]
0043         ref = property_.attrib["ref"]
0044         self.e = property_
0045         self.name = name
0046         self.ref = ref
0047         self.m = m
0048         mname = m.name
0049         mname = mname[:mname.find("0x")]
0050         self.mname = mname
0051         self.title = "%s.%s" % (self.mname,self.name)
0052
0053     def __str__(self):
0054         return "%30s :  %s" % (self.title, repr(self))
0055
0056     def __repr__(self):
0057         return "    P %30s : %30s : %s  " % ( self.name, self.ref, self.m.name  )
0058
0059 class M(object):
0060     def __init__(self, material):
0061         name = material.attrib["name"]
0062         self.name = name
0063         pp = odict()
0064         for property_ in material.xpath(".//property"):
0065             p = P(property_, self)
0066             pp[p.name] = p
0067         pass
0068         self.e = material
0069         self.pp = pp
0070
0071     def find_properties_with_ref(self, ref):
0072         return list(filter(lambda p:p.ref == ref, self.pp.values()))
0073
0074     def num_properties_with__ref(self, ref):
0075         return len(self.find_properties_with_ref(ref))
0076
0077     def __repr__(self):
0078         hdr = "M %30s : %d " % (self.name, len(self.pp))
0079         return "\n".join([hdr] + list(map(str,self.pp.values())) + [""])
0080
0081
0082 class MM(object):
0083     def __init__(self, material_elements):
0084         """
0085         :param material_elements: from lxml parsed GDML tree
0086         """
0087         d = odict()
0088         for material in material_elements:
0089             m = M(material)
0090             d[m.name] = m
0091         pass
0092         self.d = d
0093
0094     def find_materials_with_matrix_ref(self, matrix_ref):
0095         found = list(filter( lambda m:m.num_matrix_ref(matrix_ref) > 0, self.d.values()))
0096         return found
0097
0098     def find_material_properties_with_matrix_ref(self, matrix_ref):
0099         mp = []
0100         for m in self.d.values():
0101             pp = m.find_properties_with_ref(matrix_ref)
0102             mp.extend(pp)
0103         pass
0104         return mp
0105
0106     def __repr__(self):
0107         return "\n".join(map(str,self.d.values()))
0108
0109     def __call__(self, arg):
0110         """
0111         :param arg: integer or name
0112         :return m: M object
0113         """
0114         try:
0115             iarg = int(arg)
0116         except ValueError:
0117             iarg = None
0118         pass
0119         items = self.d.items()
0120         values = list(map(lambda kq:kq[1], items))
0121         return self.d[arg] if iarg is None else values[iarg]
0122
0123
0124
0125
0126 class Q(object):
0127     """
0128     Using approx hc eV to nm conversion of 1240 as it seems that was done upstream,
0129     hence using this approx value will actually be better as it should
0130     give the measurement nm from LS group.
0131
0132     Rather than using the more precise 1239.8418754199977 which
0133     will give nm slightly off
0134     """
0135     headings = tuple(["name", "ptitle", "shape", "eV[0]", "eV[-1]", "nm[0]", "nm[-1]", "len(vals)", "msg" ])
0136     fmt = "%40s : %40s : %10s : %10s : %10s : %10s : %10s : %10s : %s "
0137
0138     #hc_eVnm=1239.8418754199977  # h_Planck*c_light*1e12
0139     hc_eVnm=1240.0
0140
0141     def __init__(self, name, v, sv, msg, mm):
0142
0143         self.name = name
0144         self.v = v
0145         self.sv = sv
0146         self.mm = mm
0147         self.len = len(v)
0148
0149         vals = v[:,1]
0150         mev = v[:,0]      # CLHEP SystemOfUnits megaelectronvolt = 1.  g4-cls SystemOfUnits
0151         ev = 1.e6*v[:,0]
0152         is_optical = ev.min() >= 1.0 and ev.max() <= 20.0    # reasonal energy range for optical property
0153         is_zero_in_domain = np.any(ev == 0.)
0154
0155         if is_optical: msg += " optical, "
0156         if is_zero_in_domain: msg += " zero_in_domain, "
0157
0158         if not is_zero_in_domain:
0159             nm = self.hc_eVnm/ev
0160         else:
0161             nm = np.zeros(len(v), dtype=v.dtype)
0162         pass
0163
0164         a = np.zeros( v.shape, dtype=v.dtype )
0165         a[:,0] = nm[::-1]
0166         a[:,1] = vals[::-1]
0167         self.a = a
0168
0169         self.msg = msg
0170         self.mev = mev
0171         self.ev = ev
0172         self.nm = nm
0173         self.is_optical = is_optical
0174         self.is_zero_in_domain = is_zero_in_domain
0175         self.property = self.get_property(self.name)
0176         self.ptitle = self.property.title if not self.property is None else "-"
0177
0178     def get_property(self, matrix_ref):
0179         props = self.mm.find_material_properties_with_matrix_ref(matrix_ref)
0180         assert len(props) < 2
0181         return props[0] if len(props) == 1 else None
0182
0183     def row(self):
0184         fmt_ = lambda f:"%10.4f" % f
0185         return tuple([self.name, self.ptitle, str(self.v.shape), fmt_(self.ev[0]), fmt_(self.ev[-1]), fmt_(self.nm[0]), fmt_(self.nm[-1]), len(self.sv),self.msg])
0186
0187     @classmethod
0188     def hdr(cls):
0189         return cls.fmt % cls.headings
0190
0191     def __str__(self):
0192         return self.fmt % self.row()
0193
0194     def __repr__(self):
0195         return self.fmt % self.row()
0196
0197     def plot(self):
0198         fig, ax = plt.subplots(1)
0199         a = self.a
0200         x = a[:,0]
0201         y = a[:,1]
0202         ax.plot(x, y)
0203         ax.text(x.max(),y.max(), self.ptitle, ha='right' )
0204         print(self)
0205         print(a.shape)
0206         print(a)
0207         fig.show()
0208         return ax
0209
0210
0211
0212 class QQ(object):
0213     def __init__(self, matrix_elements, mm):
0214         """
0215         :param g: lxml parsed GDML tree
0216         """
0217         self.mm = mm
0218         d = odict()
0219         for matrix in matrix_elements:
0220             q = self.parse_matrix(matrix)
0221             d[q.name] = q
0222         pass
0223         self.d = d
0224         self._select_keys = None
0225
0226     def parse_matrix(self, matrix):
0227         coldim = matrix.attrib["coldim"]
0228         assert coldim == "2"
0229         name = matrix.attrib["name"]
0230         #name = name[:name.find("0x")]
0231
0232         sv = matrix.attrib["values"]
0233         v = values_(sv)
0234
0235         if len(v) % 2 == 0:
0236             msg = ""
0237         else:
0238             v = v[:-1]
0239             msg = "truncated+trimmed?"
0240         pass
0241         v = v.reshape(-1,2)
0242         q = Q(name, v, sv, msg, self.mm)
0243         return q
0244
0245     def __call__(self, arg):
0246         """
0247         :param arg: integer or name
0248         :return q: Q object
0249         """
0250         try:
0251             iarg = int(arg)
0252         except ValueError:
0253             iarg = None
0254         pass
0255         items = self.sorted_items()
0256         values = list(map(lambda kq:kq[1], items))
0257         return self.d[arg] if iarg is None else values[iarg]
0258
0259     def unselect(self):
0260         self._select_keys = None
0261
0262     def select(self, keys_include=[], keys_exclude=[], optical=True):
0263         """
0264         :param keys_include: list of fnmatch patterns to include from key selection
0265         :param keys_exclude: list of fnmatch patterns to exclude from key selection
0266         :param optical: when true requires properties to have energy in optical range 1-20 eV
0267         :return select_keys_: list of selected keys
0268
0269         Has side effect of setting self._select_keys which inflences
0270         other presentation of contents.
0271         """
0272         filter_q_optical = lambda q: q.is_optical == optical
0273         filter_k_include = lambda k: any(fnmatch(k, include_pattern) for include_pattern in keys_include)
0274         filter_k_exclude = lambda k: any(fnmatch(k, exclude_pattern) for exclude_pattern in keys_exclude)
0275         filter_func = lambda kq:filter_q_optical(kq[1]) and filter_k_include(kq[0]) and not filter_k_exclude(kq[0])
0276
0277         items = self.sorted_items()
0278         _select_items = list(filter(filter_func, items))
0279         _select_keys  = list(map(lambda kq:kq[0], _select_items))
0280         self._select_keys = _select_keys
0281         return self._select_keys
0282
0283     def num_select(self):
0284         return len(self._select_keys)
0285
0286     def nm_min_(self):
0287         items = self.selected_items()
0288         return np.array(list(map(lambda kq:kq[1].nm.min(), items)))
0289     def nm_min(self):
0290         return self.nm_min_().min()
0291     def nm_max_(self):
0292         items = self.selected_items()
0293         return np.array(list(map(lambda kq:kq[1].nm.max(), items)))
0294     def nm_max(self):
0295         return self.nm_max_().max()
0296
0297     def find(self, arg):
0298         """
0299         :param arg: eg RINDEX
0300         :return qq: list of all Q objects starting with arg
0301         """
0302         items = self.sorted_items()
0303         keys   = list(map(lambda kq:kq[0], items))
0304         ksel = filter(lambda k:k.startswith(arg), keys)
0305         return list(map(self, ksel))
0306
0307     def sorted_items(self):
0308         items = list(self.d.items())
0309         items = list(sorted(items,key=lambda kq:kq[1].len, reverse=True))
0310         return items
0311
0312     def selected_items(self):
0313         items = self.sorted_items()
0314         if not self._select_keys is None:
0315             items = list(filter(lambda kq:kq[0] in self._select_keys, items))
0316         pass
0317         return items
0318
0319     def __repr__(self):
0320         hdr =  Q.hdr()
0321         items = self.selected_items()
0322         values = list(map(lambda kq:kq[1], items))
0323         nm_min = "nm_min:"+str(self.nm_min())
0324         nm_max = "nm_max:"+str(self.nm_max())
0325         return "\n".join(map(str,[hdr]+values+[hdr,nm_min,nm_max]))
0326
0327
0328
0329     def plot(self):
0330         num = self.num_select()
0331         print("num:%s" % num)
0332         print("\n".join(self._select_keys))
0333         fig, axs = plt.subplots(num, sharex=True)
0334
0335         nm_max = self.nm_max()
0336
0337         for i in range(num):
0338             q = self(i)
0339             a = q.a
0340             x = a[:,0]
0341             y = a[:,1]
0342             ymid = (y.min() + y.max())/2.
0343             ax = axs[i]
0344             ax.plot(x, y)
0345             ax.text(nm_max,ymid, q.ptitle, ha='right' )
0346         pass
0347         fig.show()
0348
0349
0350
0351 if __name__ == '__main__':
0352
0353    np.set_printoptions(suppress=True)
0354
0355    #path = sys.argv[1]
0356    path = "/usr/local/opticks/opticksaux/export/juno2102/tds_ngt_pcnk_sycg_202102_v0.gdml"
0357    g = parse_(path)
0358
0359    mm = MM(g.xpath("//material"))
0360    qq = QQ(g.xpath("//matrix"), mm)
0361
0362    print("all optical properties")
0363    qq.select(keys_include=["*"],optical=True)
0364    print(qq)
0365
0366    print("all properties without selection")
0367    qq.unselect()
0368    print(qq)
0369
0370    wild = "*ABSLENGTH*"
0371    print("all optical properties matching wildcard %s " % wild )
0372    qq.select(keys_include=[wild], optical=True)
0373    print(qq)
0374
0375    qq.plot()
0376
0377    #q = qq.find("PPOABSLENGTH")[0]
0378    #a = q.a
0379
0380
0381