plots.py
Go to the documentation of this file.
1 import numpy
2 import matplotlib.pyplot as plt
3 import matplotlib as mpl
4 from matplotlib.collections import LineCollection
5 from matplotlib.ticker import AutoMinorLocator
6 from matplotlib.colors import LogNorm
7 
8 from collections import defaultdict
9 from wirecell import units
10 from wirecell.util.wires.schema import plane_face_apa;
11 
12 
13 class Hist2D(object):
14  def __init__(self, nx, xmin, xmax, ny, ymin, ymax):
15  self.arr = numpy.zeros((ny, nx))
16  self.nx = nx
17  self.rangex = (xmin, xmax)
18  self.ny = ny
19  self.rangey = (ymin, ymax)
20 
21  def xbin(self, x):
22  xmin,xmax=self.rangex
23  xrel = max(0, min(1.0, (x - xmin) / (xmax-xmin)))
24  return int(xrel * self.nx);
25 
26  def ybin(self, y):
27  ymin,ymax=self.rangey
28  yrel = max(0, min(1.0, (y - ymin) / (ymax-ymin)))
29  return int(yrel * self.ny);
30 
31  def fill(self, x, y, v):
32  xi = self.xbin(x)
33  yi = self.ybin(y)
34  self.arr[yi, xi] += v
35 
36  def extent(self):
37  return (self.rangex[0], self.rangex[1],
38  self.rangey[1], self.rangey[0])
39  def imshow(self, ax):
40  return ax.imshow(self.arr, extent=self.extent())
41 
42  def like(self):
43  return Hist2D(self.nx, self.rangex[0], self.rangex[1],
44  self.ny, self.rangey[0], self.rangey[1])
45 
46 def activity(cm):
47  '''
48  Given a ClusterMap, return a figure showing activity
49  '''
50  channels = set()
51  slices = set()
52  for snode in cm.nodes_oftype('s'):
53  sdat = cm.gr.nodes[snode]
54  for c in sdat['activity'].keys():
55  channels.add(int(c))
56  slices.add(sdat['ident'])
57 
58  cmin = min(channels);
59  cmax = max(channels);
60  smin = min(slices);
61  smax = max(slices);
62  print ("activity: c:[%d,%d], s:[%d,%d]" % (cmin, cmax, smin, smax))
63 
64  hist = Hist2D(smax-smin+1, smin, smax+1,
65  cmax-cmin+1, cmin, cmax+1)
66  for snode in cm.nodes_oftype('s'):
67  sdat = cm.gr.nodes[snode]
68  si = sdat['ident']
69  for c,v in sdat['activity'].items():
70  ci = int(c)
71  hist.fill(si+.1, ci+.1, v)
72 
73  return hist
74  # fig,ax = plt.subplots(nrows=1, ncols=1)
75  # im = hist.imshow(ax)
76  # plt.colorbar(im, ax=ax)
77  # return fig,ax,hist
78 
79 def blobs(cm, hist):
80  for bnode in cm.nodes_oftype('b'):
81  bdat = cm.gr.nodes[bnode]
82  for wnode in cm.gr[bnode]:
83  wdat = cm.gr.nodes[wnode]
84  if wdat['code'] != 'w':
85  continue;
86 
87  hist.fill(bdat['sliceid']+0.1, wdat['chid']+.1, 1)
88  return hist
89  # fig,ax = plt.subplots(nrows=1, ncols=1)
90  # im = hist.imshow(ax)
91  # plt.colorbar(im, ax=ax)
92  # return fig,ax,hist
93 
94 
95 def mask_blobs(a, b, sel=lambda a: a < 1, extent=None):
96  '''
97  Plot the activity with blobs masked out
98  '''
99  cmap = plt.get_cmap('gist_rainbow')
100  #cmap = plt.cm.gist_rainbow
101  cmap.set_bad(color='black')
102  arr = numpy.ma.masked_where(sel(b), a)
103  fig,ax = plt.subplots(nrows=1, ncols=1)
104  fig.set_size_inches(8.5,11.0)
105  im = ax.imshow(arr, cmap=cmap, interpolation='none', extent=extent)
106  minorLocator = AutoMinorLocator()
107  ax.yaxis.set_minor_locator(minorLocator)
108  ax.tick_params(which="both", width=1)
109  ax.tick_params(which="major", length=7)
110  ax.tick_params(which="minor", length=3)
111  plt.colorbar(im, ax=ax)
112  return fig, ax
113 
114 def wire_blob_slice(cm, sliceid):
115  '''
116  Plot slice information as criss crossing wires and blob regions.
117  '''
118  from . import converter
119 
120  snodes = cm.find('s', ident=sliceid)
121  if len(snodes) != 1:
122  print ("Unexpected number of slices with ID: %d, found %d" % (sliceid, len(snodes)))
123  return
124  snode = snodes[0]
125  by_face = defaultdict(list)
126  sdata = cm.gr.nodes[snode]
127  for cdat,cval in sdata["activity"].items():
128  chid = int(cdat)
129  cnode = cm.channel(chid)
130  cdata = cm.gr.nodes[cnode]
131  wnodes = cm.neighbors_oftype(cnode, 'w')
132  if not wnodes:
133  print("No wires for channel %d" % chid)
134  for wnode in wnodes:
135  wdat = cm.gr.nodes[wnode]
136  wpid = wdat['wpid']
137  p,f,a = plane_face_apa(wpid)
138  by_face[f].append((cval, wdat))
139 
140  blob_xs_byface = defaultdict(list)
141  blob_ys_byface = defaultdict(list)
142  blob_cs_byface = defaultdict(list)
143  for bnode in cm.neighbors_oftype(snode, 'b'):
144  bdata = cm.gr.nodes[bnode]
145  cx = list()
146  cy = list()
147  cpoints = converter.orderpoints(bdata['corners'])
148  for cp in cpoints + [cpoints[0]]:
149  cx.append(cp[2]/units.m)
150  cy.append(cp[1]/units.m)
151  faceid = bdata['faceid']
152  blob_xs_byface[faceid].append(cx)
153  blob_ys_byface[faceid].append(cy)
154  blob_cs_byface[faceid].append(bdata['value'])
155 
156 
157  cmap = plt.get_cmap('gist_rainbow')
158  linewidth=0.1
159  fig,axes = plt.subplots(nrows=1, ncols=len(by_face))
160 
161 
162 
163  for ind, (faceid, wdats) in enumerate(sorted(by_face.items())):
164  ax = axes[ind]
165  ax.set_title("face %d" % faceid)
166  ax.set_xlabel("Z [m]")
167  ax.set_ylabel("Y [m]")
168 
169  xs = list()
170  ys = list()
171  cs = list()
172  for cval, wdat in wdats:
173  cs.append(cval)
174  h = wdat['head']
175  t = wdat['tail']
176  c = [0.5*(h[i]+t[i]) for i in range(3)]
177 
178  p,f,a = plane_face_apa(wdat['wpid'])
179  wid = wdat['ident']
180  wip = wdat['index']
181 
182  chid = wdat['chid']
183  toffset = (chid%5) * 0.2
184 
185  if p == 4: # collection
186  if wdat['seg'] == 0:
187  ax.text(t[2]/units.m, t[1]/units.m + toffset, "C%d" %chid, fontsize=0.2, rotation=90, va='top')
188 
189  ax.text(c[2]/units.m, c[1]/units.m-toffset, "P%d WID%d WIP%d" %(p,wid,wip), fontsize=0.2, rotation=90, va='top')
190  ax.plot([c[2]/units.m,c[2]/units.m], [c[1]/units.m, c[1]/units.m-toffset],
191  color='black', linewidth = 0.1, alpha=0.5)
192 
193  else:
194  if wdat['seg'] == 0:
195  ax.text(h[2]/units.m, h[1]/units.m - toffset, "C%d" %chid, fontsize=0.2, rotation=90, va='bottom')
196  ax.plot([h[2]/units.m,h[2]/units.m], [h[1]/units.m - toffset, h[1]/units.m],
197  color='black', linewidth = 0.1, alpha=0.5)
198 
199  ax.text(c[2]/units.m+toffset, c[1]/units.m, "P%d WID%d WIP%d" %(p,wid,wip), fontsize=0.2, rotation=0, va='top')
200  ax.plot([c[2]/units.m+toffset,c[2]/units.m], [c[1]/units.m, c[1]/units.m],
201  color='black', linewidth = 0.1, alpha=0.5)
202 
203  xs.append([t[2]/units.m, h[2]/units.m])
204  ys.append([t[1]/units.m, h[1]/units.m])
205 
206 
207  segments = [numpy.column_stack([x,y]) for x,y in zip(xs, ys)]
208  lc = LineCollection(segments, cmap=cmap, linewidth=linewidth, alpha=0.5, norm=LogNorm())
209  lc.set_array(numpy.asarray(cs))
210  ax.add_collection(lc)
211  ax.autoscale()
212  fig.colorbar(lc, ax=ax)
213 
214  segments = [numpy.column_stack([x,y]) for x,y in zip(blob_xs_byface[faceid], blob_ys_byface[faceid])]
215  lc = LineCollection(segments, linewidth=2*linewidth, alpha=0.5)
216  lc.set_array(numpy.asarray(blob_cs_byface[faceid]))
217  ax.add_collection(lc)
218 
219 
220  # norm = mpl.colors.Normalize(vmin=cmin, vmax=cmax)
221  # cb = mpl.colorbar.ColorbarBase(axes[-1], cmap=cmap, norm=norm)
222  # cb.set_label("signal charge")
223  return fig, axes
224 
def activity(cm)
Definition: plots.py:46
def xbin(self, x)
Definition: plots.py:21
def blobs(cm, hist)
Definition: plots.py:79
auto enumerate(Iterables &&...iterables)
Range-for loop helper tracking the number of iteration.
Definition: enumerate.h:69
def mask_blobs(a, b, sel=lambda a:a< 1, extent=None)
Definition: plots.py:95
def wire_blob_slice(cm, sliceid)
Definition: plots.py:114
def ybin(self, y)
Definition: plots.py:26
static int max(int a, int b)
def imshow(self, ax)
Definition: plots.py:39
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55
def plane_face_apa(wpid)
Definition: schema.py:240
auto zip(Iterables &&...iterables)
Range-for loop helper iterating across many collections at the same time.
Definition: zip.h:295
def __init__(self, nx, xmin, xmax, ny, ymin, ymax)
Definition: plots.py:14
def fill(self, x, y, v)
Definition: plots.py:31